//*********************************************
// LE COUSSIN FINANCIER
// Texte Scopos vol 16, p. 233.
//*********************************************

///// Simulation du coussin en temps continu avec sauts ////////

// Données (celles du texte)
mu=0.1;
r=0.05;
v=0.2;
S0=1;
V0=100;
P0=90;
lambda=2;
m=5;
// Pas pour le calcul des variations du coussin entre 0 et 1
d=0.05;      // pas assez grand pour que les sauts soient visibles...
t=[0:d:1]';
// On considère M trajectoires pour avoir une estimation de E(C1)
// valeur moyenne du coussin en fin de période.
M=100;

// Simulation des instants de sauts :
// -- Le nombre de sauts sur [0,1] suit une loi de Poisson de paramètre lambda,
//    on simule donc N v.a. de loi de Poisson de paramètre lambda :
N=grand(1,M,'poi',lambda);
Nmax=max(N);
Nmin=min(N);
// -- on crée une matrice à Nmax lignes et M colonnes contenant les instants
//    de sauts en utilisant le fait que ceux-ci sont uniformément répartis
//    sur [0,1] connaissant N (cf. cours sur les processus de Poisson).
T=2*ones(Nmax,M); // les instants de saut hors de [0,1] sont fixés à 2.
for j=1:M,
  T([1:N(j)],j)=rand(N(j),1,'uniform');
end
T=gsort(T,'r','i');

// Simulation de l'amplitude  des instants de sauts (loi normale N(a,sigma))
a=0;sigma=0.05;
// Pour les instants de saut hors de [0,1], on fait X=0.
X=(a+sigma*rand(T,'normal')) .* bool2s(T<=1);

// Calcul des valeurs du coussin. 
// L'application de la formule 11 est difficile car il faut simuler
// les trajectoires d'un mouvement brownien (les 'W_t' sont des
// normales réduites mais *pas* indépendantes : voir leur obtention
// comme sommes de X_{n,k}).
// Comme on devra au final discrétiser ces trajectoires, le plus simple
// est d'utiliser directement la version discrète (formule 6) pour la
// partie sans saut de la formule 11.
C=(V0-P0)*ones(t)*ones(1,M);      // valeur à t=0
[lastl,lastc]=size(C);
U=rand(C,'uniform');
Z=(-1+2*bool2s(U<0.5*ones(U))) / sqrt(lastl);
Y=(mu/lastl)*ones(U)+ v*Z;
C=C .* cumprod(ones(C) + m*Y + ((1-m)*r/lastl)*ones(C), 'c'); // formule 6
// Partie due aux sauts
// La ligne i de P contient les produits (1 + m Xn) pour les n: T_n<=t(i) 
P=ones(C);
for i=1:lastl,
  Xt= X .* bool2s(T<=t(i));
  P(i,:)=prod(ones(X)+m*Xt,'r');
end
D=C .* P;

// Représentation graphique de 2 trajectoires :
// -- une qui a le nombre min de sauts (en bleu),
// -- une qui a le nombre max de sauts (en rouge).
// N.B. : si Nmin=Nmax, on prend 2 trajectoires différentes.
// Choix des trajectoires à afficher :
J1=min(find(N==Nmin));
J2=max(find(N==Nmax));
xbasc();
xselect();
printf("\nNmin = %d", Nmin)
printf("\nNmax = %d\n", Nmax)
// Pour que plot2d montre les sauts, il faut dédoubler les abscisses des points
// où il y a un saut et mettre en ordonnées les valeurs à gauche et à droite.
// Le plus simple est de dédoubler tous les points d'abscisse et pour
// les ordonnées de donner en chaque point la valeur à gauche et la valeur
// à droite (souvent égales).
x=gsort([t;t],'g','i');
y1=zeros(x);
y2=zeros(x);
// Les valeurs à droite sont données par la matrice D
y1(2:2:length(x))=D(:,J1);
y2(2:2:length(x))=D(:,J2);
// Les valeurs à gauche sont obtenues en décalant P 
// (descendre les lignes d'un cran et mettre des '1' dans la 1ère ligne)
G=C .* [ones(1,M);P(1:lastl-1,:)];
y1(1:2:length(x)-1)=G(:,J1);
y2(1:2:length(x)-1)=G(:,J2);
plot2d(x,[y1,y2],style=[2,5])

// Estimation de la valeur moyenne du coussin en fin de période (t=1).
moyenne=sum(D(lastl,:))/M;
ecart=(sum(D(lastl,:).^2)-M*moyenne^2)/(M-1);
ecart=sqrt(ecart);
printf("\nValeur moyenne du coussin à t=1 : C1 = %f", moyenne)
printf("\nÉcart type = %f\n", ecart)