TP3

Exercice 1.1 Solution 1

En utilisant la library(mvtnorm)

library(mvtnorm)
K = matrix(c(14,8,8,6),2,2)
m=c(-1,2)
M=10000
X=rmvnorm(M,m,K)
plot(X[,1],X[,2])

C=cov(X)
moy=apply(X,2,mean)
C

##           [,1]     [,2]
## [1,] 13.951971 8.016277
## [2,]  8.016277 6.035382

moy

## [1] -0.9290129  2.0431630

###Y|X=x
Esp.cond<-function(x){mean(X[,2]) +cov(X[,1],X[,2])/var(X[,1])*(x - mean(X[,1]))}
Esp.th <- function(x){2+8/14*(x + 1)}
plot(Esp.cond,xlim=c(-5,5))
plot(Esp.th,xlim=c(-5,5),add=TRUE,col="red")

plot(X[,1],X[,2]-Esp.cond(X[,1]))

cov(X[,1],X[,2]-Esp.cond(X[,1])) ###independance

## [1] -3.853903e-16

Solution 2

X1 = X[,1]
Y1 = X[,2]
N=length(X1)
# Espérance conditionnelle de Y sachant X

Phi=cbind(rep(1,N),X1,X1^2,X1^3)
colnames(Phi)=rep("",4)
A=t(Phi)%*%Phi
B=t(Phi)%*%Y1
alpha=solve(A,B)
alpha

##           [,1]
##   2.573186e+00
##   5.773291e-01
##   2.903252e-04
##  -4.918726e-05

alphath=c(moy[2]-C[1,2]/C[1,1]*moy[1],C[1,2]/C[1,1])
alphath

## [1] 2.5769388 0.5745623

# Fonction de régression
reg=lm(Y1~X1+I(X1^2)+I(X1^3))
alpha2=reg$coeff
alpha2

##   (Intercept)            X1       I(X1^2)       I(X1^3) 
##  2.573186e+00  5.773291e-01  2.903252e-04 -4.918726e-05

Espcond=Phi%*%alpha; #N simulations indép de l'esp.cond.

# Orthogonalité 
Z=Y1-Espcond
plot(X1,Z)

cov(cbind(X1,Z))

##               X1              
## X1  1.395197e+01 -3.467590e-14
##    -3.467590e-14  1.429417e+00

Exercice 1.3

En utilisant le package Copula

library(copula)
rho=0.5
M=10000
Sigma=matrix(c(1,rho,rho,1),2,2)
cop<-normalCopula(rho)
X=rCopula(M,cop)

A la main

M=10000
rho=0.5
K=matrix(c(1,rho,rho,1),2,2)
X=rmvnorm(M,c(0,0),K)
X=pnorm(X)
V=qexp(X[,1],2)
W=qnorm(X[,2])

On change le chemin courrant pour se placer là où on a sauvegarder le fichier polynomesorthogonaux

setwd("F://TP Monte Carlo")
source("TP3_codepolynomesorthogonaux.R")
# Espérance conditionnelle de V sachant W
Phi=cbind(rep(1,N),H1(W),H2(W),H3(W))
colnames(Phi)=rep("",4)

A=t(Phi)%*%Phi
B=t(Phi)%*%V
alpha=solve(A,B)
alpha

##        [,1]
##  0.49799890
##  0.23305838
##  0.06627553
##  0.01608963

Espcond=Phi%*%alpha; #N simulations indép de l'esp.cond.
mean(Espcond) ##environ 0.5 = E(V)

## [1] 0.4969892

# Orthogonalité 
Z=V-Espcond
cov(cbind(W,Z)) # elle est diagonale

##               W              
## W  9.977547e-01 -9.455410e-17
##   -9.455410e-17  1.960614e-01

Exercice 1.3

N=10000 # nombre de simulations

X=rexp(N,1)
Y=numeric()
for (i in seq(1,N)){    
    Y[i]=(max(X[i],1)-1/2)*runif(1)
}

Phi=cbind(rep(1,N), R1(X), R2(X),R3(X), R4(X), R5(X))
colnames(Phi)=rep("",6)

A=t(Phi)%*%Phi
B=t(Phi)%*%Y
alpha=solve(A,B)
alpha

##         [,1]
##   0.43423761
##  -0.35438334
##   0.05888911
##   0.09240941
##   0.01818180
##   0.05860211

Espcond=Phi%*%alpha
hat.f<-function(x){alpha[1]+R1(x)*alpha[2]+R2(x)*alpha[3]+R3(x)*alpha[4]+R4(x)*alpha[5]+R4(x)*alpha[5]}

x=sort(X)
y=1/2*(pmax(1,x)-1/2)
plot(X,Espcond,col="blue",type="p",pch=".")
lines(x,y,col="red",type="l")

plot(X,Y)
plot(hat.f,xlim=c(0,8),add=TRUE,col="red")

# approximation d'une fonction non dérivable par des polynômes...

Exercice 1.4

N=10000 # nombre de simulations
X=runif(N);
U=rep(0,N);

f=function(t){(1+t)*exp(-t)-1+runif(1)}

for (i in seq(1,N)){
    
   U[i] = uniroot(f,c(0,20))$root 
    
}
Y=X+U

# Base classique

Phi=cbind(rep(1,N), X, X^2, X^3,X^4)
colnames(Phi)=rep("",5)

A=t(Phi)%*%Phi
B=t(Phi)%*%Y
alpha=solve(A,B)
alpha

##        [,1]
##   2.0547043
##   0.1960284
##   3.6603225
##  -5.6354527
##   2.7801557

Espcond=Phi%*%alpha;
mean(Espcond)

## [1] 2.523408

###On verifie que E(E(Y|X)) = 2+E(X) = 2.5

x=sort(X)
plot(X,Espcond,type="p")
lines(x,x+2,col="red")

Phi=cbind(rep(1,N), R1(X), R2(X), R3(X),R4(X))
colnames(Phi)=rep("",5)

A=t(Phi)%*%Phi
B=t(Phi)%*%Y
alpha2=solve(A,B)
alpha2

##        [,1]
##    42.48234
##  -180.29385
##   306.22446
##  -233.08188
##    66.72364

Espcond2=Phi%*%alpha2;
mean(Espcond2)

## [1] 2.523408

x=sort(X)
plot(X,Espcond2,type="p")
lines(x,x+2,col="red")