/************   GAUSS Applications in Finance *********************
Author: Thierry Roncalli  Thierry.Roncalli@montesquieu.u-bordeaux.f
Date: 3 Oct 96
Provided without guarantee for public non-commercial use.

Author's statement:
   Voici le code source des programmes de la premiere partie
   de l'article "Quelques applications de GAUSS en finance"
   qui sera utiliser pour la presentation de GAUSS lors du
   colloque international "Dynamique des Marches Financiers
   - Taux de change, Taux d'interet, Bourse" organise par
   EPEE-Universite d'Evry. Le code source de la deuxieme
   partie de l'article "Modelisation des series temporelles
   en finance" est en partie livre avec TSM.

   Cela evitera des erreurs de frappe, et surtout du temps perdu.

   Th. Roncalli

Rough translation (I hope):
	Here is the source code for the programs in the first part of the
	article ``Some Applications of GAUSS in Finance'' . . . (snip)
	The time series code from the second part of the article
	is included with TSM.


   Contenu de la partie "Methodes numeriques en finance" :
   ------------------------------------------------------

   poisson.src (page 2, simulation de processus de Poisson)
   poisson.prg (page 3-4, exemple)
   saut.src    (page 6, simulation de processus de diffusion avec saut)
   saut.prg    (page 7-8, exemple)
   lbo.src     (page 10, valorisation d'une option look-back par MC)
   lbo.prg     (page 12, exemple)
   asian.src   (page 13, valorisation des options asiatiques et europeennes
                         par le modele binomial de Cox, Ross et Rubinstein)
   asian.prg   (page 16, exemple)
   mcarlo.prg  (page 17, illustration de la convergence de l'estimateur
                         de Monte Carlo avec accelerateur)
   markov.src  (page 19, approximation d'une EDS par une chaine de Markov)
   markov.prg  (page 20, exemple)
   tridiag.src (page 27, algorithme tridiagonal pour la resolution d'un
                         systeme et l'inversion d'une matrice)
   edp2.src    (page 28, Resolution de l'equation differentielle partielle
                         parabolique du modele de structure par terme a
                         un seul facteur)
   coupon1.prg (page 32, premier exemple (modele de Vasicek : influence
                         du theta-schema)
   coupon2.prg (page 33, deuxieme exemple (modele de Vasicek : influence
                         du pas de discretisation k)
   hjm.src     (page 35, approximation numerique de l'integrale HJM)
   hjm.prg     (page 37, exemple page 91 d'Heath, Jarrow et Morton)
   volimp.src  (page 39, calcul de la volatilite implicite d'une option d'achat
                         sur future base sur le modele de Black)
   volimp.prg  (page 41, exemple)

****************************************************************************/


/*************************************************************
---> poisson.src
*************************************************************/

proc (2) = Poisson(lambda,t0,T,L,K);
  local h,ti,u,N;

  h = (T-t0)/L;
  ti = seqa(t0,h,L+1);

  u = rndu(L,K);
  u = u .< (lambda*h);
  u = zeros(1,K)Ýu;

  N = cumsumc(u);
  retp(ti,N);
endp;


/*************************************************************
---> poisson.prg
*************************************************************/

new;
library pgraph;
#include poisson.src;

rndseed 123;

L = 1000;
t0 = 0;
TT = 10;
lambda = 5;
{t,N} = Poisson(lambda,t0,TT,L,1);
{t,N2} = Poisson(2,t0,TT,L,1);

graphset;
  fonts("simplex simgrma");
  _pdate = "";
  _plwidth = 5; _pnum = 2;
  title("Simulation de processus de Poisson"\
        "\L__ \202l\201 = 5    _ _ \202l\201 = 2");
  xtics(t0,TT,1,0);
  ytics(0,50,10,0);
  xs = t[1:L];
  xe = t[2:L+1];
  ys = N[1:L];
  ye = N[1:L];
  e = ones(L,1);

  p1 = e~6*e~xs~ys~xe~ye~0*e~7*e~0*e;

  ys = N2[1:L];
  ye = N2[1:L];

  p2 = e~3*e~xs~ys~xe~ye~0*e~7*e~0*e;

  _pline = p1Ýp2;

  xlabel("t");
  ylabel("N]t[(\202w\201)");
  graphprt("-c=3 -cf=poiss1a.hgl -w=5");
  draw;


M = N - lambda*t;

graphset;
  fonts("simplex simgrma");
  _pdate = "";
  _plwidth = 5; _pnum = 2;
  title("Processus de Poisson precedent d'intensite \202l \201= 5"\
        "\LM]t[ = N]t[ - \202l\201t est une martingale");
  xtics(t0,TT,1,0);
  ytics(-11,1,1,0);
  xs = t[1:L];
  xe = t[2:L+1];
  ys = M[1:L];
  ye = M[1:L];
  e = ones(L,1);

  _pline = e~6*e~xs~ys~xe~ye~0*e~7*e~0*e;

  xlabel("t");
  ylabel("M]t[(\202w\201)");
  graphprt("-c=3 -cf=poiss1b.hgl -w=5");
  draw;


V = M^2 - lambda*t;

graphset;
  fonts("simplex simgrma");
  _pdate = "";
  _plwidth = 5; _pnum = 2;
  title("Processus de Poisson precedent d'intensite \202l \201= 5"\
        "\LV]t[ = M]t[ - \202l\201t est une martingale");
  xtics(t0,TT,1,0);
  ytics(-40,80,10,0);
  xs = t[1:L];
  xe = t[2:L+1];
  ys = V[1:L];
  ye = V[1:L];
  e = ones(L,1);

  _pline = e~6*e~xs~ys~xe~ye~0*e~7*e~0*e;

  xlabel("t");
  ylabel("V]t[(\202w\201)");
  graphprt("-c=3 -cf=poiss1c.hgl -w=5");
  draw;


/*************************************************************
---> saut.src
*************************************************************/

proc (3) = saut(x0,mu,sigma,kappa,lambda,t0,TT,N,Ns);
  local h,t,xt,u,dN,i,xi,ti,dN_,lambda_;
  local mu:proc,sigma:proc,kappa:proc,lambda:proc;

  h = (TT-t0)/N;
  t = seqa(t0,h,N+1);

  xt = zeros(N+1,Ns);
  xt[1,.] = x0.*ones(1,Ns);
  u = rndn(N,Ns)*sqrt(h);
  dN = zeros(N+1,Ns);

  i = 1;
  do until i > N;
    xi = xt[i,.]; ti = t[i];
    lambda_ = lambda(ti,xi);
    dN_ = rndu(1,Ns) .< (lambda_*h);
    xt[1+i,.] = xi + mu(ti,xi)*h
                   + sigma(ti,xi).*u[i,.] + kappa(ti,xi).*dN_;
    dN[1+i,.] = dN_;
    i = i + 1;
  endo;

  retp(t,xt,dN);
endp;


/*************************************************************
---> saut.prg
*************************************************************/

new;
library pgraph;
#include saut.src;

x0 = 2;

proc mu(t,x);
  retp(0.0002*x);
endp;

proc sigma(t,x);
  retp(1.09*x);
endp;

proc kappa(t,x);
  retp(1);
endp;

proc lambda1(t,x);
  retp(0);
endp;

rndseed 123;
{t,x1,dN} = saut(x0,&mu,&sigma,&kappa,&lambda1,0,10,100,1);

proc lambda2(t,x);
  if t <= 5;
    retp(1);
  else;
    retp(0.5);
  endif;
endp;

rndseed 123;
{t,x2,dN} = saut(x0,&mu,&sigma,&kappa,&lambda2,0,10,100,1);
xs2 = indexcat(dN[.,1],1);

graphset;
  fonts("simplex simgrma");
  _pdate = ""; _pnum = 2;
  title("Mouvement Brownien Geometrique : "\
        "dX]t[ = \202m\201X]t[dt + \202s\201X]t[dW]t["\
        "\LProcessus de saut : "\
        "dX]t[ = \202m\201X]t[dt + \202s\201X]t[dW]t[ +\202k\201]t[dN]t[");
  xs = t[xs2];
  e = ones(rows(xs),1);
  xe = xs;
  ys = e*0;
  ye = e*100;

  _pline = e~6*e~xs~ys~xe~ye~0*e~7*e~0*e;

  xs = t[xs2-1];
  e = ones(rows(xs),1);
  xe = xs;
  ys = e*0;
  ye = e*100;

  _pline = _plineÝ(e~6*e~xs~ys~xe~ye~0*e~7*e~0*e);
   graphprt("-c=3 -cf=saut.hgl -w=5");
  xy(t,x1~x2);


/*************************************************************
---> lbo.src
*************************************************************/

declare matrix _print                       = 1;
declare matrix _mc_simulations              = {1000,1};
declare matrix _mc_nbpoints                 = 0;
declare matrix _mc_sauvegarde               = 0;
declare string _type_option                 ?= "call";

proc mc_lbo(S,r,d,sigma,tau);
  local Ns,Nr,n,h,hsqrt,mu,x,i,j,u;
  local somme,g,k,nb_simulations,prime;
  local fh,varnames,l,old_cursor;

  old_cursor = csrtype(0);
  Ns = _mc_simulations[1];  /* Nombre de simulations pour chaque replication */
  Nr = _mc_simulations[2];  /* Nombre de replications */
  somme = 0;
  nb_simulations = 0;       /* Nombre de simulations totales */

  if _mc_nbpoints == 0;
    N = trunc(tau*365*5);
  else;
    N = _mc_nbpoints;       /* Nombre de points de discretisation */
  endif;

  h = tau/n;  hsqrt = sqrt(h);  /* h : pas de discretisation */
  mu = r-d;

  if _print == 1;
    cls;
    locate 2,6;  print chrs(218~196*ones(1,72)~191);
    l = 1; do until l>3;
      locate 2+l,6; print chrs(179~zeros(1,72)~179);
      l = l+1;
    endo;
    locate 6,6;  print chrs(192~196*ones(1,72)~217);
    locate 4,8;
    if lower(_type_option) $== "put";
      print "Valorisation d'une option de vente LOOK-BACK -- PUT -- par"\
            " Monte Carlo";
    else;
      print "Valorisation d'une option d'achat LOOK-BACK -- CALL -- par"\
            " Monte Carlo";
    endif;
    locate 8,20;  print chrs(218~196*ones(1,40)~191);
    l = 1; do until l>6;
      locate 8+l,20; print chrs(179~zeros(1,40)~179);
      l = l+1;
    endo;
    locate 14,20;  print chrs(192~196*ones(1,40)~217);
    locate 9,22; print "Prix du sous-jacent :";
    locate 9,50; print ftos(s,"%lf",10,5);
    locate 10,22; print "Taux d'interet :     ";
    locate 10,50; print ftos(r,"%lf",10,5);
    locate 11,22; print "Taux de rendement :  ";
    locate 11,50; print ftos(d,"%lf",10,5);
    locate 12,22; print "Volatilite :         ";
    locate 12,50; print ftos(sigma,"%lf",10,5);
    locate 13,22; print "Maturite :           ";
    locate 13,50; print ftos(tau,"%lf",10,5);
  endif;

  if type(_mc_sauvegarde) == 13;
    varnames = { g } ;
    create fh = ^_mc_sauvegarde with varnames,0,8;
  endif;

  j = 1;
  do until j > nr;
    x = s.*ones(1,ns);
    k = x;
    i = 1;
    do until i > N;
      u = rndn(1,Ns);
      u = hsqrt*u;
      x = x+ (mu*x)*h + (sigma*x).*u;
      if lower(_type_option) $== "put";
        K = maxc(KÝx); K = K';
      else;
        K = minc(KÝx); K = K';
      endif;
      i=i+1;
    endo;
    nb_simulations = nb_simulations + Ns;
    if lower(_type_option) $== "put";
      g = K-x;                             /* Pay Off d'un put */
    else;
      g = x-K;                             /* Pay Off d'un call */
    endif;

    g = gÝzeros(1,ns);
    g = maxc(g);
    if type(_mc_sauvegarde) == 13;
      gosub sauvegarde;
    endif;
    somme = somme + sumc(g);
    prime = exp(-r*tau)*(somme/nb_simulations);  /* Formule de Feynman-Kac */
    if _print == 1;
      gosub affichage;
    endif;
    if key == 1079;
      retp(prime);
    endif;
    j = j+1;
  endo;
  locate 23,1;
  call csrtype(old_cursor);
  retp(prime);

  affichage:
    locate 16,20;  print chrs(218~196*ones(1,40)~191);
    l = 1; do until l>5;
      locate 16+l,20; print chrs(179~zeros(1,40)~179);
      l=l+1;
    endo;
    locate 21,20; print chrs(192~196*ones(1,40)~217);
    locate 17,22; print "Nombre de points :   ";
    locate 17,50; print ftos(n,"%lf",10,0);
    locate 18,22; print "Pas :                ";
    locate 18,50; print ftos(h,"%lf",10,8);
    locate 19,22; print "Nombre de simulations :";
    locate 19,50; print ftos(nb_simulations,"%lf",10,0);
    locate 20,22; print "Valeur de la prime : ";
    locate 20,50; print ftos(prime,"%lf",10,6);
  return;

  sauvegarde:
    call writer(fh,g);
  return;

endp;


/*************************************************************
---> lbo.prg
*************************************************************/

new;
#include lbo.src;

tt = hsec;

rndseed 123;
_mc_simulations = {100,10};  /* 1000 simulations */
_type_option = "call";
_print = 1;
prime1 = mc_lbo(100,ln(1+0.10),0,0.2,1);

rndseed 123;
_mc_nbpoints = 100;
prime2 = mc_lbo(100,ln(1+0.10),0,0.2,1);

output file = lbo.out reset;

print ftos(prime1,"prime1 : %lf",10,5);
print ftos(prime2,"prime2 : %lf",10,5);

print;
print ftos((hsec-tt)/100,"Temps de calcul : %lf secondes",5,3);

output off;


/*************************************************************
---> asian.src
*************************************************************/

declare matrix _print = 1;
declare matrix _accelerateur = 1;

proc (2) = simulation(_pi,u,d,n,m);
  local w,b1,b2,p1,p2;
  w = rndu(n,m);     /* variable aleatoire uniforme      */
  b1 = w.<=_pi;      /* variable aleatoire de Bernouilli */
  w = 1 - w;
  b2 = w.<=_pi;      /* variable aleatoire antithetique  */

  p1 = b1*u + (1-b1)*d;
  p1 = cumprodc(p1);

  p2 = b2*u + (1-b2)*d;
  p2 = cumprodc(p2);

  retp(p1,p2);
endp;


proc (3) = parametre_option(S,tau,sigma,r,n);
  local u,d,_pi;
  u = exp(sigma*sqrt(tau/n));
  d = 1/u;
  _pi = (exp(r*tau/n)-d)/(u-d);
  retp(_pi,u,d);
endp;


proc PayOFF(Stilde,K,r,tau);
  local N,ST,Sm,G,Prime;

  Prime = zeros(6,1);
  N = rows(Stilde);

  ST = Stilde[N,.]';          /* Valeurs finales du sous-jacent */
  Sm = meanc(Stilde);         /* valeurs moyennes prises par le sous-jacent */

  /* Floating Strike Options --- CALL */
  G = (ST-Sm);
  G = G.*(G.>0);
  Prime[1] = meanc(G);

  /* Floating Strike Options --- PUT  */
  G = (Sm-ST);
  G = G.*(G.>0);
  Prime[2] = meanc(G);

  /* Fixed Strike Options --- CALL    */
  G = (Sm-K);
  G = G.*(G.>0);
  Prime[3] = meanc(G);

  /* Fixed Strike Options --- PUT     */
  G = (K-Sm);
  G = G.*(G.>0);
  Prime[4] = meanc(G);

  /* European CALL                    */
  G = (ST-K);
  G = G.*(G.>0);
  Prime[5] = meanc(G);

  /* European PUT                     */
  G = (K-ST);
  G = G.*(G.>0);
  Prime[6] = meanc(G);

  Prime = exp(-r*tau)*Prime;

  retp(prime);
endp;


proc CRR_asian(S,K,Tau,Sigma,r,N,M);
  local _pi,u,d,p1,p2,Stilde,prime1,prime2,prime,str;

  {_pi,u,d} = parametre_option(S,Tau,Sigma,r,N);

  {p1,p2} = simulation(_pi,u,d,n,m);

  Stilde = S*p1;            /* Premiere trajectoire simulee du sous-jacent */
  Prime1 = PayOff(Stilde,K,r,tau);

  if _accelerateur /= 1;
    retp(Prime1);
  endif;

  Stilde = S*p2;            /* Deuxieme trajectoire simulee du sous-jacent */
  Prime2 = PayOff(Stilde,K,r,tau);

  Prime = 0.5*(Prime1+Prime2);  /* Accelerateur de Monte carlo */

  if _print == 1;

    print "Floating Strike Option:";
    str = ftos(prime[1],"call = %lf",10,5);
    str = str$+ftos(prime[2],"     put = %lf",10,5);
    print str; print;

    print "Fixed Strike Option:";
    str = ftos(prime[3],"call = %lf",10,5);
    str = str$+ftos(prime[4],"     put = %lf",10,5);
    print str; print;

    print "European Option:";
    str = ftos(prime[5],"call = %lf",10,5);
    str = str$+ftos(prime[6],"     put = %lf",10,5);
    print str; print;

  endif;

  retp(Prime);
endp;


/*************************************************************
---> asian.prg
*************************************************************/

new;
#include asian.src;

output file = asian.out reset;

tt = hsec;

Prime = CRR_asian(100,98,90/365,0.2,0.08,3,10000);

print;
print ftos((hsec-tt)/100,"Temps de calcul : %lf secondes",5,3);

output off;


/*************************************************************
---> mcarlo.prg
*************************************************************/

new;
library ipg,pgraph;
ipgset;

output file = mcarlo.out reset;
tt = hsec;

#include asian.src;

MC1 = zeros(250,3);    /* avec accelerateur de Monte Carlo */
MC2 = zeros(250,3);    /* sans accelerateur de Monte Carlo */
_print = 0;

rndseed 123; _accelerateur = 1;

i = 1;
do until i > 250;
  Prime = CRR_asian(100,98,90/365,0.2,0.08,3,5);
  MC1[i,1] = Prime[5];
  Prime = CRR_asian(100,98,90/365,0.2,0.08,3,25);
  MC1[i,2] = Prime[5];
  Prime = CRR_asian(100,98,90/365,0.2,0.08,3,100);
  MC1[i,3] = Prime[5];
  i = i+1;
endo;

rndseed 123; _accelerateur = 0;

i = 1;
do until i > 250;
  Prime = CRR_asian(100,98,90/365,0.2,0.08,3,5);
  MC2[i,1] = Prime[5];
  Prime = CRR_asian(100,98,90/365,0.2,0.08,3,25);
  MC2[i,2] = Prime[5];
  Prime = CRR_asian(100,98,90/365,0.2,0.08,3,100);
  MC2[i,3] = Prime[5];
  i = i+1;
endo;


x1 = zeros(128,3); x2 = x1;
d1 = zeros(128,3); d2 = d1;

i = 1;
do until i > 3;
  {x1[.,i],d1[.,i],f} = noyau(MC1[.,i]);
  {x2[.,i],d2[.,i],f} = noyau(MC2[.,i]);
  i = i + 1;
endo;

print ftos((hsec-tt)/100,"Temps de calcul : %lf secondes",10,5);

output off;

graphset;
  TITLE("Fonction de densite (avec accelerateur de Monte Carlo)");
  _plwidth = 6;
  _pdate = "";
  _pline = 1~1~6.3957115~0~6.3957115~10~1~7~10;
  _plegstr = "Ns = 5\000Ns = 25\000Ns = 100";
  _plegctl = {2 7 6 3};
  xtics(0,18,2,0); ytics(0,1.2,0.4,0);
  graphprt("-c=3 -cf=mcarlo1.hgl -w=10");
  xy(x1,d1);

graphset;
  TITLE("Fonction de densite (sans accelerateur de Monte Carlo)");
  _plwidth = 6;
  _pdate = "";
  _pline = 1~1~6.3957115~0~6.3957115~10~1~7~10;
  _plegstr = "Ns = 5\000Ns = 25\000Ns = 100";
  _plegctl = {2 7 6 3};
  xtics(0,18,2,0); ytics(0,1.2,0.4,0);
  graphprt("-c=3 -cf=mcarlo2.hgl -w=10");
  xy(x2,d2);


/*************************************************************
---> markov.src
*************************************************************/

proc (2) = EDS_to_MarkovChain(mu,sigma,cn);
  local mu:proc,sigma:proc;
  local X_min,X_max,Nx,k,h,X,i,j,Pij,mu_,sigma_;
  local div,num,num1,num2;
  X_min = cn[1];
  X_max = cn[2];
  Nx = cn[3];
  k = cn[4];
  h = (X_max-X_min)/Nx;
  X = seqa(X_min,h,Nx+1);
  i = seqa(0,1,Nx+1); j = i';
  Pij = zeros(Nx+1,Nx+1);
  mu_ = mu(X); sigma_ = sigma(X);
  div = sigma_*sqrt(k);
  num = (j[1,2:Nx]-i)*h-mu_*k;
  num1 = (num+0.5*h)./div;
  num2 = (num-0.5*h)./div;
  Pij[.,2:Nx] = cdfn(num1)-cdfn(num2);
  num = ((-i+0.5)*h-mu_*k)./div;
  Pij[.,1] = cdfn(num);
  num = ((Nx-i-0.5)*h-mu_*k)./div;
  Pij[.,Nx+1] = cdfnc(num);
  retp(X,Pij);
endp;


/*************************************************************
---> markov.prg
*************************************************************/

new;
#include markov.src;

hsec_ = hsec;

K = 98;
S0 = 100;
sig = 0.15;
t0 = 0;
TT = 95/365;
tau = TT - t0;
r = ln(1+0.08);

Prime = 5.3677;   /* Valeur de la prime de l'option d'achat dans
                     le modele de Black et Scholes */

proc mu(x);
  retp(r*x);
endp;

proc sigma(x);
  retp(sig*x);
endp;

S_min = 80;
S_max = 120;
Ns = 100;     /* Nombre de points pour la partition de S */
Nt = 50;      /* Nombre de points pour la partition du temps */

cn = S_minÝS_maxÝNsÝ(tau/Nt);

{S,P} = EDS_to_MarkovChain(&mu,&sigma,cn);

Prob = P;

i = 1;
do until i > Nt;
  Prob = Prob*P;
  i = i + 1;
endo;

indx = indexcat(S,S0);
Prob = Prob[indx,.]';

PayOff = S - K;
PayOff = PayOff.*(PayOff.>0);
C = exp(-r*tau)*(Prob'PayOff);

output file = markov.out reset;

print ftos(Prime,"Valeur de la prime de l'option d'achat : %lf",5,4);
print ftos(C,"Valeur approchee par la chaine de Markov : %lf",5,4);

print;
print ftos((hsec-hsec_)/100,"Temps de calcul : %lf secondes",5,3);

output off;


/*************************************************************
---> tridiag.src
*************************************************************/

proc Tridiagonal(RR,r);
  local N,x,b_prime,r_prime,i,_aux;
  N = rows(RR);
  x = zeros(n,1); b_prime = zeros(N,1); r_prime = zeros(N,1);
  b_prime[n] = RR[n,n]; r_prime[n] = r[n];
  i = N-1;
  do while i>0;
    _aux = RR[i,i+1]/b_prime[i+1];
    b_prime[i]=RR[i,i]-_aux*RR[i+1,i];
    r_prime[i]=r[i]-_aux*r_prime[i+1];
    i=i-1;
  endo;
  x[1]=r_prime[1]/b_prime[1];
  i=2;
  do until i>n;
    x[i]=(r_prime[i]-RR[i,i-1]*x[i-1])/b_prime[i];
    i=i+1;
  endo;
  retp(x);
endp;

proc Tridiag_inv(A);
  local N,Id,A_inv,j;
  N = rows(A); Id = eye(N); A_inv = zeros(N,N);
  j=1;
  do until j>N;
    A_inv[.,j] = Tridiagonal(A,Id[.,j]);
    j=j+1;
  endo;
  retp(A_inv);
endp;


/*************************************************************
---> edp2.src
*************************************************************/

declare external _print;

proc (3) = EDP2(mu,sigma,lambda,alpha,theta,cn);
  local mu:proc,sigma:proc,lambda:proc,alpha:proc;
  local x_min,x_max,J,tau,N,Temps,h,k,old,t,x,r_j;
  local d0,d1,d2,_Upsilon,_Lambda;
  local sigma_j,mu_j,lambda_j,_s1,_s2,_s3,a_j,b_j,c_j,d_j,e_j,f_j;
  local M,u_n,u,nn,t1,t2,t_n;
  x_min = cn[1]; x_max = cn[2]; J = cn[3];
  tau = cn[4]; N = cn[5]; Temps = cn[6];
  k = tau/N; h = (x_max-x_min)/J;
  old = csrtype(0);
  if _print == 1;
    call EDP2_affichage(NÝJÝkÝhÝtheta);
  endif;
  if Temps /= 1;
    Temps = 0;
  endif;
  t = seqa(0,k,N+1); x = seqa(x_min,h,J+1);
  d0 = seqa(1,1,J+1); d1 = seqa(2,1,J); d2 = seqa(1,1,J);
  _Upsilon = zeros(J+1,J+1);
  _Lambda = zeros(J+1,J+1);

  if temps == 0;
    sigma_j = sigma(x); mu_j = mu(x); lambda_j = lambda(x); r_j = alpha(x);
    _s1 = (sigma_j)^2*k/(h^2); _s2 = (mu_j-lambda_j.*sigma_j)*k/h;_s3 = 1-theta;
    a_j = 0.5*(_s1+_s2);
    b_j = 1-_s1-r_j*k;
    c_j = 0.5*(_s1-_s2);
    d_j = a_j; e_j = b_j-2; f_j = c_j;
    clear sigma_j,mu_j,lambda_j;
    A_j = _s3*a_j; B_j = _s3*(b_j-1)+1; C_j = _s3*c_j;
    D_j = theta*d_j; E_j = theta*(e_j+1)-1; F_j = theta*f_j;
    E_j[1] = F_j[1]+E_j[1]; E_j[J+1] = E_j[J+1]+D_j[J+1];
    _Upsilon = Tridiag_create(F_j,E_j,D_j);
    _Lambda = Tridiag_create(C_j,B_j,A_j);
    clear A_j,B_j,C_j,D_j,E_j,F_j;
    M = -inv(_Upsilon)*_Lambda; /* M = -Tridiag_inv(_Upsilon)*_Lambda; */
    clear _Upsilon,_Lambda;
    u_n = ones(J+1,1);
    u = zeros(J+1,N+1); u[.,1] = u_n;
    nn = 1;
    t1 = hsec;
    do until nn>N;
      u_n = M*u_n;
      u[.,nn+1] = u_n;
      if _print == 1;
        t2 = hsec;
        call EDP2_print(t1,t2,nn,N);
      endif;
      nn = nn+1;
    endo;
  else;
    u_n = ones(J+1,1);
    u = zeros(J+1,N+1); u[.,N+1] = u_n;
    t1 = hsec;
    t_n = t[N+1];
    sigma_j = sigma(t_n,x); mu_j = mu(t_n,x); lambda_j = lambda(t_n,x);
    r_j = alpha(t_n,x);
    _s1 = (sigma_j)^2*k/(h^2); _s2 = (mu_j-lambda_j.*sigma_j)*k/h;
    _s3 = 1-theta;
    clear sigma_j,mu_j,lambda_j;
    a_j = 0.5*(_s1+_s2);
    b_j = 1-_s1-r_j*k;
    c_j = 0.5*(_s1-_s2);
    nn = N;
    do while nn>=1;
      A_j = _s3*a_j; B_j = _s3*(b_j-1)+1; C_j = _s3*c_j;
      _Lambda = Tridiag_create(C_j,B_j,A_j);
      t_n = t[nn];
      sigma_j = sigma(t_n,x); mu_j = mu(t_n,x); lambda_j = lambda(t_n,x);
      r_j = alpha(t_n,x);
      _s1 = (sigma_j)^2*k/(h^2); _s2 = (mu_j-lambda_j.*sigma_j)*k/h;
      a_j = 0.5*(_s1+_s2);
      b_j = 1-_s1-r_j*k;
      c_j = 0.5*(_s1-_s2);
      clear sigma_j,mu_j,lambda_j;
      D_j = theta*a_j; E_j = theta*(b_j-1)-1; F_j = theta*c_j;
      E_j[1] = F_j[1]+E_j[1]; E_j[J+1] = E_j[J+1]+D_j[J+1];
      _Upsilon = Tridiag_create(F_j,E_j,D_j);
      clear D_j,E_j,F_j;
      u_n = Tridiagonal(_Upsilon,-_Lambda*u_n);
      clear _Upsilon,_Lambda;
      u[.,nn] = u_n;
      if _print == 1;
        t2 = hsec;
        call EDP2_print(t1,t2,N-nn+1,N);
      endif;
      nn = nn-1;
    endo;
    u=rev(u')';
  endif;
 retp(t,x,u);
endp;


proc (1) = Tridiag_create(A0,A1,A2);
  local m,A,i;
  m = rows(A1); A = zeros(m,m);
  A[1,1] = A1[1]; A[1,2] = A2[1];
  i=2;
  do until i>m-1;
    A[i,i-1] = A0[i];
    A[i,i] = A1[i];
    A[i,i+1] = A2[i];
    i=i+1;
  endo;
  A[m,m-1] = A0[m]; A[m,m] = A1[m];
  retp(A);
endp;


proc (0) = EDP2_affichage(cn);
  local q0,q1,q2,w1,w2,w3,w4;
  cls;
  q0 = 196*ones(1,70); q1 = zeros(1,70); q2 = zeros(1,5);
  w1 = chrs(q2~218~q0~191); w2 = chrs(q2~179~q1~179); w3 = chrs(q2~192~q0~217);
  print w1; print w2; print w2; print w2; print w2; print w2; print w3;
  locate 3,17; print "Resolution de la valorisation des coupons zeros";
  locate 4,13;
  print "dans le modele financier d'arbitrage a une variable d'etat";
  locate 5,23; print "par la methode des theta-schemas.";
  locate 10,1;
  q0 = 196*ones(1,50); q1 = zeros(1,50); q2 = zeros(1,16);
  w1 = chrs(q2~218~q0~191); w2 = chrs(q2~179~q1~179); w3 = chrs(q2~192~q0~217);
  w4 = chrs(q2~195~q0~180);
  print w1; print w2; print w2; print w4; print w2; print w4;
  print w2; print w2; print w2; print w3;
  locate 11,22; print ftos(cn[1],"N : %lf    ",7,0);
  locate 11,45; print ftos(cn[2],"J : %lf    ",7,0);
  locate 12,22; print ftos(cn[3],"k : %lf    ",7,5);
  locate 12,45; print ftos(cn[4],"h : %lf    ",7,5);
  locate 14,22; print ftos(cn[5],"THETA : %lf    ",7,5);
  retp;
endp;


proc (0) = EDP2_print(t1,t2,nn,N);
  local t3;
  t3=(t2-t1)/6000;
  locate 16,22; print ftos(t3,"Temps ecoule : %lf minutes",5,2);
  t3=t3/nn*(N-nn);
  locate 17,22; print ftos(t3,"Temps restant (estimation) : %lf minutes",5,2);
  locate 18,22; print ftos(nn,"Iteration no %lf",5,0);
  retp;
endp;


/*************************************************************
---> coupon1.prg
*************************************************************/

new;
library pgraph;
#include tridiag.src;
#include edp2.src;

proc mu(x);
  retp(0.10-x);
endp;

proc sigma(x);
  retp(0.1);
endp;

proc lambda(x);
  retp(0.1);
endp;

proc alpha(x);
  retp(x);
endp;

_print =1;
cn = -1Ý1Ý100Ý5Ý500Ý0;

{t,x,u} = EDP2(&mu,&sigma,&lambda,&alpha,0,cn);
r = x[56]; sol=u[56,.]';
{t,x,u} = EDP2(&mu,&sigma,&lambda,&alpha,0.25,cn); sol=sol~u[56,.]';
{t,x,u} = EDP2(&mu,&sigma,&lambda,&alpha,0.5,cn); sol=sol~u[56,.]';
{t,x,u} = EDP2(&mu,&sigma,&lambda,&alpha,0.75,cn); sol=sol~u[56,.]';
{t,x,u} = EDP2(&mu,&sigma,&lambda,&alpha,1,cn); sol=sol~u[56,.]';

P = vasicek(t,1,0.10,0.1,r,0.1);

graphset;
  fonts("simplex simgrma");
  _pdate = ""; _plwidth = 5; _paxht = 0.2; _pnum = 2;
  title("Methode des \202Q\201-schemas");
  xlabel("\202t\201"); ylabel("P-u[n]]j[");
  _plegstr = "\202Q\201=0\000\202Q\201=0.25\000\202Q\201=0.5\000"\
             "\202Q\201=0.75\000\202Q\201=1";
  _plegctl = { 2 6 6 1.5};
  graphprt("-c=3 -cf=edp1.hgl -w=5");
  xy(t,P-sol);


proc vasicek(tau,a,b,sigma,r0,lamda);
  local b_prime,taux_long,x,c;
  b_prime=b-lamda.*sigma./a;
  taux_long=b_prime-(sigma^2)./(2*a^2);
  x=(1-exp(-a.*tau))./a;
  c=exp(-taux_long.*tau+(taux_long-r0).*x-(x^2).*(sigma^2)./(4*a));
  retp(c);
endp;


/*************************************************************
---> coupon2.prg
*************************************************************/

new;
library pgraph;
#include tridiag.src;
#include edp2.src;

proc mu(x);
  retp(0.10-x);
endp;

proc sigma(x);
  retp(0.1);
endp;

proc lambda(x);
  retp(0.1);
endp;

proc alpha(x);
  retp(x);
endp;

_print =1;

np = seqa(100,100,5);
axe1 = miss(zeros(501,5),0);
sol1 = miss(zeros(501,5),0);

i=1;
do until i>5;
  cn = -1Ý1Ý100Ý5Ýnp[i]Ý0;
  {t,x,u} = EDP2(&mu,&sigma,&lambda,&alpha,0.25,cn);
  r = x[56];
  P = vasicek(t,1,0.10,0.1,r,0.1);
  sol1[1:rows(P),i] = P-u[56,.]';
  axe1[1:rows(t),i] = t;
  i = i+1;
endo;

graphset;
  fonts("simplex simgrma");
  _pdate = ""; _plwidth = 5; _paxht = 0.2; _pnum = 2;
  title("Influence de k");
  xlabel("\202t\201"); ylabel("P-u[n]]j[");
  _plegstr = "k=0.05\000k=0.025\000k=0.0167\000"\
             "k=0.0125\000k=0.01";
  _plegctl = { 2 6 2 4};
  graphprt("-c=3 -cf=edp2.hgl -w=5");
  xy(axe1,sol1);


proc vasicek(tau,a,b,sigma,r0,lamda);
  local b_prime,taux_long,x,c;
  b_prime=b-lamda.*sigma./a;
  taux_long=b_prime-(sigma^2)./(2*a^2);
  x=(1-exp(-a.*tau))./a;
  c=exp(-taux_long.*tau+(taux_long-r0).*x-(x^2).*(sigma^2)./(4*a));
  retp(c);
endp;


/*************************************************************
---> hjm.src
*************************************************************/

declare matrix _n = 20; /* Nombre de points d'approximation
                           pour l'algorithme de SIMPSON */

/* Variables globales */

declare matrix _beta;  /* Fonction beta(t,T) */
declare matrix _A1;
declare matrix _A2;

proc _HJM_1(v);
  local f,res,u;
  u = _A1;
  f = _beta;
  local f:proc;
  res = f(u,v);
  retp(res);
endp;

proc _HJM_2(u,t);
  local h,v,s,I,j;
  _A1 = u;
  h = (t-u)/_n;
  v = seqa(u,h,_n+1);
  s = 1Ýreshape(4~2,1,_n-1)'Ý1;
  I = 0;
  j = 1;
  do until j > _n+1;
    I = I+h*(s[j]/3)*_HJM_1(v[j]);
    j = j+1;
  endo;
  retp(I);
endp;

proc _HJM_3(u);
  local res,t;
  t = _A2;
  res = _HJM_2(u,t);
  retp(res);
endp;

proc _HJM_4(u);
  local f,res,t;
  t = _A2;
  f = _beta;
  local f:proc;
  res = f(u,t)*_HJM_3(u);
  retp(res);
endp;

proc HJM(beta,t0,t);
  local beta:proc;
  _beta = &beta;
  local h,u,s,I,j;
  _A2 = t;
  h = (t-t0)/_n;
  u = seqa(t0,h,_n+1);
  s = 1Ýreshape(4~2,1,_n-1)'Ý1;
  I = 0;
  j = 1;
  do until j > _n+1;
    I = I+h*(s[j]/3)*_HJM_4(u[j]);
    j = j+1;
  endo;
  retp(I);
endp;


/*************************************************************
---> hjm.prg
*************************************************************/

new;
library pgraph;

#include heath.src;

sigma2 = 0.05;
lambda = 0.125;

proc beta(t,_T);
  local res;
  res = sigma2.*exp(-lambda/2.*(_T-t));
  retp(res);
endp;

proc solution(t);
  local y;
  y = -2*(sigma2/lambda)^2*( (1-exp(-lambda*t)) - 2*(1-exp(-lambda/2*t)) );
  retp(y);
endp;

N = 5Ý10Ý25Ý50;
ch = "Sol";

j = 1;
do until j > 4;

_N = N[j];
_T = 140;
t = seqa(0.1,0.5,_T);
r = zeros(_T,2);
i = 1;
do until i > _T;
  r[i,.] = HJM(&beta,0,t[i])~solution(t[i]);
  i = i+1;
endo;
call varput(r,ch$+ftos(j,"%lf",1,0));
j = j+1;
endo;

graphset;
begwind;
window(2,2,0);
setwind(1);
  _ptitlht = 0.25; _paxht = 0.25; _pnum = 2; _pnumht = 0.20; _plwidth = 5;
  title("N=5");
  _plegstr = "Solution Numerique\0Solution Exacte";
  _plegctl = { 2 6 3.5 2};
  xlabel("t");
  r = varget("Sol1");
  xy(t,r);
nextwind;
  title("N=10");
  r = varget("Sol2");
  xy(t,r);
nextwind;
  graphset;
  _ptitlht = 0.25; _paxht = 0.25; _pnum = 2; _pnumht = 0.20; _plwidth = 5;
  title("N=25");
  xlabel("t"); ylabel("ERREUR");
  r = varget("Sol3");
  xy(t,r[.,2]-r[.,1]);
nextwind;
  title("N=50");
  r = varget("Sol4");
  xy(t,r[.,2]-r[.,1]);
  graphprt("-c=3 -cf=hjm.hgl -w=5");
endwind;


/*************************************************************
---> volimp.src
*************************************************************/

declare matrix _volimp_K;
declare matrix _volimp_F0;
declare matrix _volimp_tau;
declare matrix _volimp_r;
declare matrix _volimp_Prime;
declare matrix _convergence = 0.001;


proc Volimp(Prime,K,F0,tau,r,a,b);

  _volimp_Prime = Prime;
  _volimp_K = K;
  _volimp_F0 = F0;
  _volimp_tau = tau;
  _volimp_r = r;

  retp( bisection(&_CALL_Future,a,b) );
endp;


proc CALL_future(K,F0,sigma,tau,r);
  local a,b,d1,d2,prime;

  a = sigma.*sqrt(tau);
  d1 = ln(F0./K)./a + 0.5*a;
  d2 = d1 - a;
  b = exp(-r.*tau);


  prime = F0.*b.*cdfn(d1) - K.*b.*cdfn(d2);

  retp(prime);
endp;


proc _CALL_Future(sigma);
   retp(Call_Future(_volimp_K,_volimp_F0,sigma,_volimp_tau,_volimp_r)
         - _volimp_Prime);
endp;


proc bisection(&f,a,b);
  local f:proc;
  local ya,yb,Nobs,c,yc,indx1,indx2;

  ya = f(a); yb = f(b);
  Nobs = rows(ya);

  if ( maxc(ya) > 0 ) or ( minc(yb) < 0 );
    ERRORLOG "erreur : mauvaise initialisation de a et b.";
    end;
  endif;

  do while maxc(abs(a-b)) > _convergence;
    c = (a+b)/2;
    yc = f(c);
    indx1 = yc.<0;
    indx2 = 1 - indx1;
    a = indx1.*c + indx2.*a;
    b = indx1.*b + indx2.*c;
  endo;
  c = (a+b)/2;
  retp(c);
endp;


/*************************************************************
---> volimp.prg
*************************************************************/

new;
#include volimp.src;

load data;

Prime = data[.,1];
K = data[.,2];
F0 = data[.,3];
r = data[.,4];
tau = data[.,5];

tt = hsec;

output file = volimp.out reset;

_convergence = 0.00001;

sigma = Volimp(Prime,K,F0,tau,r,0.001,0.5);

print ftos((hsec-tt)/100,"Temps de calcul : %lf secondes",6,2);

output off;