Exercice 1.1
En utilisant la library(mvtnorm)
library(mvtnorm)
## Warning: package 'mvtnorm' was built under R version 3.4.1
M=3000
rho=0.8
K=matrix(c(1,rho,rho,1),2,2)
X=rmvnorm(M,c(0,0),K)
U=pnorm(X)
A la main
A=t(chol(K))
X=matrix(rep(0,M*2),M,2)
for(m in 1:M)
X[m,]=A%*%rnorm(2)
U=pnorm(X)
plot(U[,1],U[,2]) ##copule normale, marges uniformes

hist(U[,1],freq=FALSE)
x=seq(0,1,0.001)
lines(x,dunif(x),col="red")

hist(U[,2],freq=FALSE)
lines(x,dunif(x),col="red")

cor(U,method="spearman") ##Estimation du rho de Spearman
## [,1] [,2]
## [1,] 1.0000000 0.7966535
## [2,] 0.7966535 1.0000000
6/pi*asin(rho/2)
## [1] 0.7859393
Exercice 1.2
Copule normale, marges exponentielle et pareto
ppareto<-function(x){1-(1/x)^3}
qpareto<-function(x){(1-x)^(-1/3)}
U[,1] = qexp(U[,1],2) ###marge exponentielle
U[,2] = qpareto(U[,2]) ###marge pareto
cor(U,method="spearman")
## [,1] [,2]
## [1,] 1.0000000 0.7966535
## [2,] 0.7966535 1.0000000
plot(U[,1],U[,2]) ###Copule normale et marges exponentielle et pareto

On verifie que les marge sont bien respectivement exponentielle e pareto
plot(seq(0,4,0.01),dexp(seq(0,4,0.01),2),col="red",type='l')
hist(U[,1],freq=FALSE,add=TRUE)

x=seq(0,4,0.01)
plot(ecdf(U[,1]))
lines(x,pexp(x,2),col="red")

plot(ecdf(U[,2]))
lines(x,ppareto(x),col="red")

EXercice 1.3
Simulation du vecteur gaussien
M=1000
rho=-0.4
K=matrix(c(1,rho,rho,1),2,2)
X=rmvnorm(M,c(0,0),K)
head(X)
## [,1] [,2]
## [1,] 0.65125685 0.9101235
## [2,] -0.01854413 -1.5392342
## [3,] 0.06164945 -1.9503202
## [4,] -0.90066316 0.0532397
## [5,] -0.47882739 -1.4598861
## [6,] 1.58112531 -0.9203356
plot(X[,1],X[,2])

La mĆŖme mais vectoriellement (sans boucle for)
nu = 3
M=5000
rho=0.8
K=matrix(c(1,rho,rho,1),2,2)
X=rmvnorm(M,c(0,0),K)
S = rchisq(M,nu)
Y =X/sqrt(S/nu)
U=pt(Y,nu)
plot(U[,1],U[,2])

#ppareto<-function(x){1-(1/x)^3}
qpareto<-function(x){(1-x)^(-1/3)}
U[,1] = qexp(U[,1],2)
U[,2] = qpareto(U[,2])
plot(U[,1],U[,2])

PAckage Copula
##Exercice 1.4
library(copula)
## Warning: package 'copula' was built under R version 3.4.1
norm.cop <- normalCopula(rho)
gumb.cop <- gumbelCopula(3,2) ##gumbel 2-dimensionnelle de parametre theta=3
X=rCopula(M,norm.cop)
X[,1] = qexp(X[,1],2)
X[,2] = qpareto(X[,2])
plot(X[,1],X[,2])

t.cop <- tCopula(rho,3)
X=rCopula(M,t.cop)
X[,1] = qexp(X[,1],2)
X[,2] = qpareto(X[,2])
plot(X[,1],X[,2])

X=rCopula(M,gumb.cop)
X[,1] = qexp(X[,1],2)
X[,2] = qpareto(X[,2])
plot(X[,1],X[,2])

Copule multivariƩe
gumb.cop <- gumbelCopula(3,4) ##gumbel 4-dimensionnelle de parametre theta=3
X=rCopula(M,gumb.cop)
head(X)
## [,1] [,2] [,3] [,4]
## [1,] 0.25671821 0.10666880 0.378944259 0.15567088
## [2,] 0.62887948 0.69441422 0.726218558 0.77740402
## [3,] 0.07156137 0.03534631 0.058612380 0.16880067
## [4,] 0.18167689 0.05623169 0.461654598 0.21001083
## [5,] 0.23219228 0.11204877 0.009643894 0.02108046
## [6,] 0.13815131 0.31264631 0.498693733 0.22739622