TP2

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])

simulation d’un vecteur ayant pour copule la copule de student et marge uniforme

nu=4
X=A%*%rnorm(2)
S = rchisq(1,nu)
Y= X/sqrt(S/nu)
Y

##            [,1]
## [1,] -0.2875782
## [2,]  0.9603518

U = pt(Y,nu)
U

##           [,1]
## [1,] 0.3939768
## [2,] 0.8043690

simulation de M vecteurs ayant pour copule la copule de student et marge uniforme

A=t(chol(K))
nu = 1
M=1000
rho=0.8
K=matrix(c(1,rho,rho,1),2,2)


Y=matrix(rep(0,M*2),M,2)
for(m in 1:M)
{

X=A%*%rnorm(2)
S = rchisq(1, nu)
Y[m,] = X/sqrt(S/nu)
U = pt(Y,nu)
}
plot(U[,1],U[,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