Chapitre 5
Conclusion générale
Nous avons tout au long de ce mémoire intitulé
"Simulation de modèles de diffusion appliqués aux taux
d'intérêts" proposé des algorithmes de simulation de
processus de diffusion appliquées aux taux d'intérêts ainsi
que leurs programmes de simulation sous le logiciel R.
R est un langage de programmation et un environnement
mathématique utilisé pour le traitement de données et
l'analyse statistique. C'est un projet GNU fondé sur le langage S et sur
l'environnement développé dans les laboratoires Bell par John
Chambers et ses collègues.
R est considéré par ses créateurs comme
étant une exécution de S, avec la sémantique
dérivée du langage Scheme. R est un logiciel libre
distribué selon les termes de la licence GNU GPL et est disponible sous
GNU/Linux, FreeBSD, NetBSD, OpenBSD, Mac OS X et Windows. R représente
aujourd'hui l'un des objectifs techniques majeurs de la communauté
hacker GNU1.
Les processus stochastiques de taux que nous avons
étudié sont le modèle de Vasicek, le modéle de Cox,
Ingersoll & Ross ainsi que le modèle de Vasicek a 2 facteurs, qui
sont les modèles les plus utilisés pour la dynamique des taux
d'intérêts.
Des programmes en R pour la simulation du processus
d'Orstein-Uhlenbeck exponentiel ainsi que les modèles de Heston et
d'Eraker ont été ajoutés et proposés au lecteur.
Nous espérons avoir atteint l'objectif fixé, a
savoir initier le lecteur a la simulation des processus de diffusion
appliqués aux taux d'intérêts par l'utilisation de l'outil
informatique.
Chapitre 6
Annexes
6.1 Annexe I : Programmes en R
6.1.1 Code R pour simulation d'un mouvement brownien standard
SimulationBrownien<-function(n){
# n ... nombre de simulations
temps = seq(0,1,length=n+1) ## discretisation du temps
pas.temps = 1/n
B.acc = rnorm(n,sd=sqrt(pas.temps)) ## Simulation des
accroissements B.sim = c(0,cumsum(B.acc))} ## Simulation d'une trajectoire
6.1.2 Code R pour simulation de 30 trajectoires d'un mouvement
brownien
SimulTrajectoires<-function(n.sim,n.point){ # n.sim ... nombre
de simulation
# n.point ... points de discrétisation
temps = seq(0,1,length=n.point)
pas.temps = 1/(n.point-1)
B.acc =
matrix(rnorm((n.point-1)*n.sim,sd=sqrt(pas.temps)),nrow=n.sim) B.sim =
matrix(NA,ncol=n.point,nrow=n.sim)
for (i in 1 :n.sim)
{B.sim[i,] = c(0,cumsum(B.acc[i,]))}
B.mean = apply(B.sim,2,mean)} ## apply calcule la moyenne pour
chaque colonne
6.1.3 Code R pour covariance empirique du mouvement brownien
B.cov.emp = cov(B.sim) ## estimation de la matrice de
variance-covariance image(temps, temps,
B.cov.emp,col=terrain.colors(20),xlab="temps",ylab="temps", main="Covariance
empirique du mouvement Brownien")
contour(temps, temps,B.cov.emp, add=TRUE)
6.1.4 Code R pour simulation d'une marche aléatoire
SimulationRandomWalk<-function(n){
# n ... nombre de simulations
temps = seq(0,1,length=n-l-1) ## discretisation du temps
pas.temps = 1/n
M.tr = runif(30,-1,1)
M.sim = c(0,cumsum(M.tr))}
6.1.5 Code R pour simulation d'un mouvement brownien de dimension
2
SimulBrowDim2<-function(n,rho){
# n ... nombre de simulations
# rho ... coéffi cientde corrélation
R<-matrix(c(1 ,rho,rho,1),nrow=2) ## Matrice des
Corrélation
L<- t(chol(R)) ## decomposition de Cholesky
Z<- matrix(rnorm(n*2),nrow=2)
Z<- L%*%Z ## variables aléatoires
corrélées
Z2=t(Z) ## transposé de la matrice Z pas.temps = 1/n
B.acc1 = sqrt(pas.temps)*Z2[,1] B.sim1 = c(0,cumsum(B.acc1))
B.acc2 =sqrt(pas.temps)*Z2[,2] B.sim2 = c(0,cumsum(B.acc2))}
6.1.6 Code R pour simulation d'un mouvement brownien de dimension
3
SimulBrowDim3<-function(n,rho){
# n ... nombre de simulations
# rho ... coéffi cientde corrélation
R<-matrix(c(1 ,rho,rho,rho,1,rho,rho,rho,1),nrow=3)
L<- t(chol(R))
Z<- matrix(rnorm(n*3),nrow=3)
Z<- L%*%Z Z3=t(Z)
pas.temps = 1/n
B.acc1 = sqrt(pas.temps)*Z3[,1] B.sim1 = c(0,cumsum(B.acc1))
B.acc2 = sqrt(pas.temps)*Z3[,2] B.sim2 = c(0,cumsum(B.acc2)) B.acc3 =
sqrt(pas.temps)*Z3[,3] B.sim3 = c(0,cumsum(B.acc3))}
6.1.7 Code R pour simulation d'un mouvement brownien avec
derive
SimulationBrownienDrifté<-function(n,u){
# n ... nombre de simulations
# u ... Coéfficient de dérive
temps = seq(0,1,length=n+1) ## discretisation du temps
pas.temps = 1/n
B.acc =rnorm(n,sd=sqrt(pas.temps)) ## simulation des
accroissements y=u*pas.temps+(B.acc)
B.sim = c(0,cumsum(y))} ## simulation d'une trajectoire
6.1.8 Code R pour simulation de 30 trajectoires d'un brownien
avec derive
SimulTrajectoires<-function(n.sim,n.point,u){ # n.sim ...
nombre de simulation
# n.point ... points de discrétisation
# u ... Coéfficient de dérive
temps = seq(0,1,length=n.point)
pas.temps = 1/(n.point-1)
B.acc =
matrix(rnorm((n.point-1)*n.sim,sd=sqrt(pas.temps)),nrow=n.sim) B.sim =
matrix(NA,ncol=n.point,nrow=n.sim)
for (i in 1 :n.sim)
{B.sim[i,] = c(0,cumsum(u*pas.temps+(B.acc[i,])))}
B.mean = apply(B.sim,2,mean)} ## apply calcule la moyenne pour
chaque colonne
6.1.9 Code R pour simulation d'un mouvement brownien
géométrique
SimulBroGéo<-function(para,x0,dt,n){ # x0 ... valeur
intiale du processus
# n ... nombre de simulations
# u ... coéfficient de dérive
# sigma ... volatilité
u=para[1]
sigma=para[2] dt=1/n
x=c(0,n)
tvector=dt*0 :n dW=sqrt(dt)*rnorm(n)
x[1]=x0
for(i in 1 :n) {x[i+1]=x[i]+mu*x[i]*dt+sigma*x[i]*dW[i]}}
6.1.10 Code R pour estimation des paramètres du
modèle Vasicek
EstimationVasicek<-function(data,dt){
#data ... vecteur contenant les données des taux #dt ...
intervalle de temps entre les données
N=length(data)
rate=data[2 :N]
lagrate=data[1 :(N-1)]
alpha=(N*sum(rate*lagrate) - sum(rate)*sum(lagrate))/
(N*sum(lagrate^2)- (sum(lagrate)) 2)
k_hat = -log(alpha)/dt
theta_hat = sum(rate-alpha*lagrate ) / (N*(1-alpha))
v2hat<-sum((rate-lagrate*alpha-theta_hat*(1-alpha)) 2)/N
sigma_hat<-sqrt(2*k_hat*v2hat/(1-exp(-2*k_hat*dt)))
c(theta_hat,k_hat,sigma_hat)}
6.1.11 Code R pour simulation du modèle de Vasicek
VasicekSimulation<-function(para,r0,dt,n,m){ # r0 ... valeur
intiale du processus de taux
# n ... nombre de simulations
# dt ... pas du temps
# m ... nombre de trajectoires simulées k=para[1]
theta=para[2]
sigma=para[3]
r=matrix(0,m+1,n)
r[1,]=r0
for(j in 1 :n){
for(i in 2 :(m+1)){
dr=k*(theta-r[i-1,j])*dt + sigma*sqrt(dt)*rnorm(1,0,1)
r[i,j]=r[i-1,j] + dr }}}
6.1.12 Code R pour simulation de la solution du modèle de
Vasicek
VasicekSimulSolution<-function(para,r0,dt,n,m){
# r0 ... valeur intiale du processus de taux
# n ... nombre de simulations
# dt ... pas du temps
# m ... nombre de trajectoires simulées
k=para[1]
theta=para[2] sigma=para[3] r=matrix(0,m+1,n)
r[1,]=r0
for(j in 1 :n){
for(i in 2 :(m+1)){
dr=k*(1-exp(-sigma*dt))+(theta*sqrt((1-exp(-2*sigma*dt))/(2*sigma))*rnorm(1,0,1))
r[i,j]=r[i-1,j]*exp(-sigma*dt)+dr }}}
6.1.13 Code R pour le calcul du prix du bon dans le modèle
de Vasicek
VasicekPrix<-function(r,tau,para){
# r ... valeur actuelle du taux d'intérêt
# Para ... vecteur des paramètres du modèle de
Vasicek
# tau ... maturité
theta=Para[1]
k=Para[2]
sigma=Para[3]
B=(1-exp(-k*tau))/k
A=exp((B-tau)*(k 2*theta-0.5*sigma 2)/k 2 - sigma 2*B 2/(4*k))
P=A*(exp(-B*r))} # P : prix du bon du trésor
6.1.14 Code R pour estimation des paramètres du
modèle CIR
CIRloglike<-function(para,data){
#Calcul de la fonction log-vraissemblance du modèle CIR
theta=para[1]
k=para[2]
sigma= para[3]
N=length(data)
rate=data[1 :(N-1)]
lagrate=data[2 :N]
ncp= rate*((4*k*exp(-k*dt))/(sigma 2*(1-exp(-k*dt)))) d=
4*theta*k/sigma 2 #degr" de liberté
c= 4*k/(sigma 2*(1-exp(-k*dt)))
R=sum(dchisq(c*lagrate,df=d,ncp=ncp,log=TRUE)+log(c))}
#Code R pour l'optimisation Numérique
OptimCIR<-optim(par=c(0.1,0.1,0.1),fn=CIRloglike,method="L-BFGS-B",
lower=c(0.05,0.05,0.05),upper=c(0.7,0.7,0.7),data=read.table("dataUS.dat"))$par
6.1.15 Code R pour simulation du modèle de Cox, Ingersoll
& Ross
CIRsimulation<-function(para,r0,n,dt,m){
# r0 ... valeur intiale du processus de taux
# n ... nombre de simulations
# dt ... pas du temps
# m ... nombre de trajectoires simulées
k=para[1]
theta=para[2] sigma=para[3] r=matrix(0,m+1,n)
r[1,]=r0
for(j in 1 :n){
for(i in 2 :(m+1)){
dr=k*(theta-r[i-1,j])*dt +sigma*sqrt(
r[i-1,j]*dt)*rnorm(1,0,1)
r[i,j]=r[i-1,j] + dr }}}
6.1.16 Code R pour simulation de la solution exacte du
modèle CIR
CIRsimulSolution<-function(para,n,dt,m){ # r0 ... valeur
intiale du processus de taux # n ... nombre de simulations
# dt ... pas du temps
# m ... nombre de trajectoires simulées
# c... terme multiplicatif distribution khi-2 # df...
degré de liberté
# ncp... parametre de non centralité k=para[1]
theta=para[2]
sigma=para[3]
c= sigma 2*(1-exp(-k*dt))/(4*k)
df= 4*theta*k/sigma 2
r=c(0,m)
r[1]=r0
for(i in 1 :(m)){
ncp= r[i]*exp(-k*dt)/c
r[i+1]=c*rchisq(m,df=df,ncp=ncp)}}
6.1.17 Code R pour le calcul du prix du bon dans le
modèle CIR
CIRprix<-function(r,tau,para){
# r ... valeur actuelle du taux d'intérêt
# Para ... vecteur des paramètres du modèle CIR
# tau ... maturité
theta=Para[1]
k=Para[2]
sigma=Para[3]
h=sqrt(k 2+2*sigma 2)
B= 2*(exp(h*tau)-1) / (2*h+(k+h)*(exp(tau*h)-1))
A= ((2*h*exp((k+h)*(tau)/2)) / (2*h+(k+h)*(exp(tau*h)-1)))
(2*k*theta/sigma 2) P=A*(exp(-B*r))} # P : prix du bon du trésor
6.1.18 Code R pour estimation des paramètres du
modèle de Vasicek 2f
EstimationVasicek2f<-function(data,data2,dt){ #data ...
vecteur des données du taux court
#data2 ... vecteur des données du taux a long terme #dt
... intervalle de temps entre les données N=length(data)
rate=data[2 :N]
lagrate=data[1 :(N-1)]
alpha=(N*sum(rate*lagrate) - sum(rate)*sum(lagrate))/
(N*sum(lagrate 2)- (sum(lagrate)) 2) k_hat = -log(alpha)/dt
theta_hat = sum(rate-alpha*lagrate ) / (N*(1-alpha))
v2hat<-sum((rate-lagrate*alpha-theta_hat*(1-alpha)) 2)/N
sigma_hat<-sqrt(2*k_hat*v2hat/(1-exp(-2*k_hat*dt)))
c(theta_hat,k_hat,sigma_hat)
N=length(data2)
rate=data2[2 :N]
lagrate=data2[1 :(N-1)]
alpha2=(N*sum(rate*lagrate) - sum(rate)*sum(lagrate))/
(N*sum(lagrate^2)- (sum(lagrate)) 2)
k_hat2 = -log(alpha2)/dt
theta_hat2 = sum(rate-alpha2*lagrate ) / (N*(1-alpha2))
v2hat2<-sum((rate-lagrate*alpha2-theta_hat2*(1-alpha2)) 2)/N
sigma_hat2<-sqrt(2*k_hat2*v2hat2/(1-exp(-2*k_hat2*dt)))
c(theta_hat2,k_hat2,sigma_hat2)
corr=cor(T,T2)} #coéffi cient de corrélation
6.1.19 Code R pour simulation du modèle de Vasicek a 2
facteurs
Vasicek2fSimulation<-function(param,x0,y0,dt,n,d){ # x0 ...
valeur intiale du processus de taux x
# y0 ... valeur intiale du processus de taux y # n ... nombre de
simulations
# dt ... pas du temps
# m ... nombre de trajectoires simulées k_x =param[1]
theta_x =param[2]
sigma_x=param[3]
k_y =param[4]
theta_y =param[5] sigma_y=param[6] rho =param[7]
x=matrix(0,m+1,n) y=matrix(0,m+1,n) r=matrix(0,m+1,n)
x[1,]=x0
y[1,]=y0
r[1,]=x[1]+y[1]
for (j in 1 :(n)){
for(i in 2 :(m+1)){ norm1=rnorm(1,0,1) norm2=rnorm(1,0,1)
norm3=(rho*norm1)+(sqrt(1-rho 2)*norm2)
dx=k_x*(theta_x-x[i-1,j])*dt+sigma*sqrt(dt)*norm3
x[i,j]=x[i-1,j] + dx
dy=k_y*(theta_y-y[i-1,j])*dt + sigma_y*sqrt(dt)*norm1
y[i,j]=y[i-1,j] + dy r[i,j]=x[i,j]+y[i,j] }}}
6.1.20 Code R pour simulation de la solution exacte du Vasicek
2f
Vasicek2fSimulSolution<-function(param,x0,y0,dt,n,d){
# x0 ... valeur intiale du processus de taux x
# y0 ... valeur intiale du processus de taux y
# n ... nombre de simulations
# dt ... pas du temps
# m ... nombre de trajectoires simulées
k_x =param[1]
theta_x =param[2] sigma_x=param[3] k_y =param[4]
theta_y =param[5] sigma_y=param[6] rho =param[7]
x=matrix(0,m+1,n) y=matrix(0,m+1,n) r=matrix(0,m+1,n)
x[1,]=x0
y[1,]=y0
r[1,]=x[1,]+y[1,]
for (j in 1 :(n)){ for(i in 2 :(m+1)){ norm1=rnorm(1,0,1)
norm2=rnorm(1,0,1)
norm3=(rho*norm1)+(sqrt(1-rho 2)*norm2)
dx=k_x*(1-exp(-sigma_x*dt))+(theta_x*sqrt((1-exp(-2*sigma_x*dt))/(2*sigma_x))*norm3)
x[i,j]=x[i-1,j]*exp(-sigma_x*dt)+dx
dy=k*(1-exp(-sigma_y*dt))+(theta_y*sqrt((1-exp(-2*sigma_y*dt))/(2*sigma_y))*norm1)
y[i,j]=y[i-1,j]*exp(-sigma_y*dt)+dy
r[i,j]=x[i,j]+y[i,j] III
6.1.21 Code R pour le calcul du prix du bon dans Vasicek 2f
PrixVasicek2F<-function(t,para,x,y,tau){
# t ... valeur du taux d'intérêt (modele
vasicek2f)
# x ... valeur du processus y(t)
# y ... valeur du processus x(t)
# Para ... vecteur des parametres du modele Vasicek2f)
# tau ... maturité
k_x =para[1]
theta_x =para[2]
sigma_x=parm[3]
k_y =parm[4]
theta_y =para[5]
sigma_y=para[6]
rho =para[7]
R=theta_x+(1-exp(-k_x*tau))/(k_x*tau)*(x-theta_x)+theta_y+(1-exp(-k_y*tau))/
(k_y*tau)*(y-theta_y)-sigma_x-2/(2*k_x-2)*(1+(1-exp(-2*k_x*tau))/
(2*k_x*tau)-2*(1-exp(-k_x*tau))/(k_x*tau))-sigma_y-2/(2*k_y-2)*(1+(1-exp(-2*k_y*tau))/
(2*k_y*tau)-2*(1-exp(-k_y*tau))/(k_y*tau))-(rho*sigma_x*sigma_y)/(k_x*k_y)*(1-
(1-exp(-k_x*tau))/
(k_x*tau)-(1-exp(-k_y*tau))/(k_y*tau)+(1-exp(-(k_x+k_y)*tau))/((k_x+k_y)*tau))
P=exp(-R*t)} #P : prix du bon du trésor
6.1.22 Code R pour simulation d'un processus d'O-U exponentiel
SimulationOUexponentiel<-function(para,rO,dt,n,m){
# rO ... valeur intiale du processus OUexpo
# n ... nombre de simulations
# dt ... pas du temps
# m ... nombre de trajectoires simulées
k=para[1]
theta=para[2] sigma=para[3] r=matrix(O,m+1,n)
r[1,]=log(rO) for(j in 1 :n){
for(i in 2 :(m+1)){
dr=k*(theta-log(r))*dt + beta*sqrt(dt)*rnorm(1,O,1)
r[i,j]=r[i-1,j] + dr }}}
6.1.23 Code R pour simulation du modèle d'Eraker de
volatilité stochastique
SimulationEraker<-function(para,rO,vO,m,dt,n){ # vO ...
valeur intiale du processus de volatilité # rO ... valeur intiale du
modèle proposé par Eraker
# m ... nombre de simulations
# dt ... pas du temps
# n ... nombre de trajectoires simulées
k=para[1]
theta=para[2]
sigma=para[3]
rho=para[4]
v=matrix(0,m+1,n) r=matrix(0,m+1,n) v[1,]=v0
r[1,]=r0
for(j in 1 :n){ for(i in 2 :(m+1)){
norm1=rnorm(1,0,1) norm2=rnorm(1,0,1)
norm3=(rho*norm1)+(sqrt(1-rho 2)*norm2)
dv=k*v[i-1,j]*dt +sigma*sqrt(dt)*norm3
v[i,j]=v[i-1,j] + dv
dr=(theta+k*r[i-1,j])*dt +sigma*exp(0.5*v[i-1,j])*norm1
r[i,j]=r[i-1,j] + dr }}}
6.1.24 Code R pour simulation du modèle de Heston de
volatilité stochastique
Simulationlleston<-function(para,s0,v0,r,m,dt,n){ # v0 ...
valeur intiale du processus de volatilité
# s0 ... valeur intiale du modèle de Heston
# r ... taux d'intérêt sans risque
# m ... nombre de simulations
# dt ... pas du temps
# n ... nombre de trajectoires simulées
k=para[1]
theta=para[2]
sigma=para[3] #volatilité de la volatilité
rho=para[4]
v=matrix(0,m+1,n) s=matrix(0,m+1,n) v[1,]=v0
s[1,]=s0
for(j in 1 :n){
for(i in 2 :(m+1)){ norm1=rnorm(1,0,1) norm2=rnorm(1,0,1)
norm3=(rho*norm1)+(sqrt(1-rho 2)*norm2)
dv=k*(theta-v[i-1,j])*dt +sigma*sqrt(v[i-1,j]*dt)*norm3
v[i,j]=v[i-1,j] + dv
ds=r*s[i-1,j]*dt +s[i-1,j]*sqrt(v[i-1,j]*dt)*norm1
s[i,j]=s[i-1,j] + ds }}}
6.2 Annexe II : Les données
6.2.1 Daily Treasury Bill Rates Data
Nous présentons ici les données s'étallant
sur la période allant du 05/01/2012 au 21/05/2012.
| 
