/*------------------------------------------------------------------*/
/*                       NPDE                                       */
/*===================================================================*/
/* Solves Differential Equations for the Upper-Lower Bound of a call
   option backwards from CT = max(ST-X,0) to find option price bounds
   at other dates in a model with:
       -Stochastic Volatility
       -Stochastic Interest Rates.
======================================================================

NOTE: It uses an EXPLICIT METHOD.
      Therefore, the time grid must be much finer than price grid or
      second deriv. becomes unstable. Condition: (dt)/(dS)^2 < 1/2
      It allows to implement  changes of variable for the S, V and r
*/

format /RD 6,2;
library pgraph;
graphset;
output file=difq.out reset;


/* -------------Characterization of Changes of Variables ------------*/

chvar=zeros(3,1);
c=zeros(3,1);

/*---------------- Definition of Grids for Variables ---------*/

u1g = seqb(10,500,5);       /* grid for u1 --> S */
Sg = U1toS(u1g,1);          /* grid for stock prices */

u2g = seqb(0.1,0.3,0.01);       /* grid for u2 --> V */
Vg = U1toS(u2g,2);          /* grid for volatility process */

u3g = seqb(0.01,0.12,0.02);       /* grid for u3 --> r */
rg = U1toS(u3g,3);          /* grid for interest rate */

dt = 0.00005;                 /* time grid for iterations */

X = 100;                  /* strike price */
_shr = 0.24;              /* parameter of mu S. see fn below: It implies
                             a constant Sharpe Ratio for the Stock */
A = 0.48;                  /* global sharpe ratio */
lambda= 0.;                /* market price of volatility risk */
maxt = 1;                 /* maximum time: Time to expiration */

                          /* type of smoothing for 2nd derivatives */
wind = 2;                 /* 1st deriv is 1/2 fwd + 1/2 backwd
                             wind=0: 2nd deriv is 1/2 fwd 1/2 backwd on 1st d
                             wind=1: average this point and 2 neighbors
                             wind=2: do it again ... */

nS = rows(Sg);
nV = rows(Vg);
nr = rows(rg);

/* -------------  specify stochastic processes ---------- */

      /* Parameter values */

_alphav=1.4;          /* volatility mean reversion */
_alphar=1.214;          /* interest rate mean reversion */
_V0=0.149;               /* long term mean volatility */
iVi = indexcat(Vg,(_V0-0.001)|(_V0+0.001));    /* change to index */
iV=vg[ivi];
_v0=iv;
ivi;
_r0=0.07;              /* long term mean interest rate */
r = indexcat(rg,(_r0-0.01)|(_r0+0.01));    /* change to index */
_r0=rg[r];
_r0;
sigmaV=0.2;            /* total volatility of volatility */
corrSV=0.1;             /* correlation coeff. between S and V */
_betaV=corrSV*sigmaV;
_betas=sqrt(sigmaV^2-_betaV^2);
sigmar=0.1;            /* total volatility of interest rates*/
corrSr=0.1;             /* correlation coeff. between S and r */
corrVr=0.1;             /* correlation coeff. between V and r */
_sigrS=corrSr*sigmar;
_sigrV=corrVr*sigmar;
_sigrr=sqrt(sigmar^2-_sigrV^2-_sigrS^2);
_betas=0.001;
_betaV=0.5;
_sigrs=0.001;
_sigrV=0.001;
_sigrr=0.18;


/*---------------- dynamics --------------------------------------------*/

proc(1) = mus(v,r);      /* dS = S mus(.) dt + S (V)^(1/2) dz */
    retp(_shr*V^(1/2)+r);
endp;

proc(1) = muv(v);    /* dV = muv(.) dt + sigvz(.)dz + sigvV(.) dwv */
    retp(_alphav*(_V0-V));
endp;

proc(1) = sigvz(v);
    retp(_betas*V^(1/2));
endp;

proc(1) = sigvV(v);
    retp(_betav*V^(1/2));
endp;

proc(1) = mur(r);    /* dr = mur(.)dt + sigrz(.)dz + sigrV(.)dwv + sigrr(.)dwr*/
    retp(_alphar*(_r0-r));
endp;

proc(1) = sigrz(r);
    retp(_sigrs*r^(1/2));
endp;

proc(1) = sigrV(r);
    retp(_sigrv*r^(1/2));
endp;

proc(1) = sigrr(r);
    retp(_sigrr*r^(1/2));
endp;

/*-----Initializations--------------------------------------------*/

uC=zeros(nr*nS,nV);
dC=zeros(nr*nS,nV);
SVC=zeros(nS,nV);
BSC=zeros(nS,1);
Cu1=zeros(nr*nS,nV);
Cu1u1=zeros(nr*nS,nV);
Cu2=zeros(nr*nS,nV);
Cu2u2=zeros(nr*nS,nV);
Cu3=zeros(nr*nS,nV);
Cu3u3=zeros(nr*nS,nV);
Cu1u2=zeros(nr*nS,nV);
Cu1u3=zeros(nr*nS,nV);
Cu2u3=zeros(nr*nS,nV);


/*------ Terminal Condition: option value at expiration ----------*/

i=1;
do while i <= nr;
 uC[nS*(i-1)+1:nS*i,1:nV] = maxc( (Sg-x)'|zeros(1,nS) )*ones(1,nV);
 dC[nS*(i-1)+1:nS*i,1:nV] = uC[nS*(i-1)+1:nS*i,1:nV];
                   /* s on rows, v on cols */
 i = i+1;
endo;
SVC= uC[1:nS,1:nV];
BSC= uC[1:nS,1];



         /* Setting Time Counter */

bigp = 0;                 /* intermediate print-outs along the way */

time0 = time;
timesec = 3600|60|1|.01;   /* converts all times to seconds */
time0 = timesec'time0;     /* seconds at start */

/*-----------  Option bounds at Grid Values -----------------*/

t = dt;
do while t <= maxt;        /* time before expiration */

                /* Evaluating Computer Time for Computations */
  if abs(trunc((t/dt)/10) - (t/dt)/10) < 0.0000001;
    output off;
    format /RD 4,0;
    timev = time;
    timev = timesec'timev; /* current time in seconds */
    "" ;
    "Running simulation" t/dt "of" maxt/dt;
     timprsim = (timev-time0)/(t/dt);
     sectogo = (maxt-t)*timprsim/dt;
    "  estimated  time per 100 sims h,m,s" sectohms(100*timprsim)';
    "          estimated time to go h,m,s" sectohms(sectogo)' ;
    output on;
  endif;

/* ====> Computing the Upper Bound */

    /* find derivatives Cu1 Cu2 Cu3 Cu1u1 Cu2u2 Cu3u3 Cu1u2 Cu1u3 Cu2u3
       by finite differences */
    /* take 1/2 of forward and backward differences */

   i=1;
   do while i <= nr;

    /* Cu1 and Cu1u1 */

    Cu1[1+nS*(i-1):nS*i-1,1:nV] =
           (uC[2+nS*(i-1):nS*i,1:nV]-
           uC[1+nS*(i-1):nS*i-1,1:nV])
           ./(u1g[2:nS]-u1g[1:nS-1]);
    Cu1[1+nS*(i-1):nS*i,1:nV] =
          ((Cu1[1+nS*(i-1),1:nV]|
            Cu1[1+nS*(i-1):nS*i-1,1:nV]) +
           (Cu1[1+nS*(i-1):nS*i-1,1:nV]|
            Cu1[nS*i-1,1:nV]))/2;

    Cu1u1[1+nS*(i-1):nS*i-1,1:nV] =
            (Cu1[2+nS*(i-1):nS*i,1:nV]-
            Cu1[1+nS*(i-1):nS*i-1,1:nV])
            ./(u1g[2:nS]-u1g[1:nS-1]);
    Cu1u1[1+nS*(i-1):nS*i,1:nV] =
            ((Cu1u1[1+nS*(i-1),1:nV]|
             Cu1u1[1+nS*(i-1):nS*i-1,1:nV]) +
            (Cu1u1[1+nS*(i-1):nS*i-1,1:nV]|
             Cu1u1[nS*i-1,1:nV]))/2;

    j = 1;
    do while j <= wind;
        Cu1u1[1+nS*(i-1):nS*i,1:nV] =
              ((Cu1u1[1+nS*(i-1),1:nV]|
                Cu1u1[1+nS*(i-1):nS*i-1,1:nV])+
               Cu1u1[1+nS*(i-1):nS*i,1:nV]
               +(Cu1u1[2+nS*(i-1):nS*i,1:nV]|
               Cu1u1[nS*i,1:nV]))/3;
        j = j+1;
    endo;

    /* Cu2 and Cu2u2 */

    Cu2[1+nS*(i-1):nS*i,1:nV-1] =
         (uC[1+nS*(i-1):nS*i,2:nV]
          -uC[1+nS*(i-1):nS*i,1:nV-1])
          ./((u2g[2:nV]-u2g[1:nV-1])');
    Cu2[1+nS*(i-1):nS*i,1:nV] =
          ((Cu2[1+nS*(i-1):nS*i,1]~
            Cu2[1+nS*(i-1):nS*i,1:nV-1]) +
           (Cu2[1+nS*(i-1):nS*i,1:nV-1]~
            Cu2[1+nS*(i-1):nS*i,nV-1]))/2;

    Cu2u2[1+nS*(i-1):nS*i,1:nV-1] =
         (Cu2[1+nS*(i-1):nS*i,2:nV]-
          Cu2[1+nS*(i-1):nS*i,1:nV-1])
          ./((u2g[2:nV]-u2g[1:nV-1])');
    Cu2u2[1+nS*(i-1):nS*i,1:nV] =
         ((Cu2u2[1+nS*(i-1):nS*i,1]~
           Cu2u2[1+nS*(i-1):nS*i,1:nV-1]) +
          (Cu2u2[1+nS*(i-1):nS*i,1:nV-1]~
           Cu2u2[1+nS*(i-1):nS*i,nV-1]))/2;

    /* Cu1u2 */

    Cu1u2[1+nS*(i-1):nS*i,1:nV-1] =
         (Cu1[1+nS*(i-1):nS*i,2:nV]-
          Cu1[1+nS*(i-1):nS*i,1:nV-1])
          ./((u2g[2:nV]-u2g[1:nV-1])');
    Cu1u2[1+nS*(i-1):nS*i,1:nV] =
         ((Cu1u2[1+nS*(i-1):nS*i,1]~
           Cu1u2[1+nS*(i-1):nS*i,1:nV-1]) +
          (Cu1u2[1+nS*(i-1):nS*i,1:nV-1]~
           Cu1u2[1+nS*(i-1):nS*i,nV-1]))/2;

    /* Cu3 */

    if i==1;
      Cu3[1+nS*(i-1):nS*i,1:nV] =
            (uC[1+nS*i:nS*(i+1),1:nV]-
             uC[1+nS*(i-1):nS*i,1:nV])
            ./(u3g[i+1]-u3g[i]);
    elseif i==nr;
      Cu3[1+nS*(i-1):nS*i,1:nV] =
            (uC[1+nS*(i-1):nS*(i),1:nV]-
             uC[1+nS*(i-2):nS*(i-1),1:nV])
            ./(u3g[i]-u3g[i-1]);
    else;
      Cu3[1+nS*(i-1):nS*i,1:nV] =
            (uC[1+nS*i:nS*(i+1),1:nV]-
             uC[1+nS*(i-1):nS*i,1:nV])
            ./(u3g[i+1]-u3g[i]);
      Cu3[1+nS*(i-1):nS*i,1:nV] =
           (Cu3[1+nS*(i-1):nS*i,1:nV]+
            (uC[1+nS*(i-1):nS*(i),1:nV]-
             uC[1+nS*(i-2):nS*(i-1),1:nV])
            ./(u3g[i]-u3g[i-1]))/2;
    endif;

    i = i+1;
  endo;



   i=1;
   do while i <= nr;

    /* Cu3u3  */

     if i==1;
      Cu3u3[1+nS*(i-1):nS*i,1:nV] =
            (Cu3[1+nS*i:nS*(i+1),1:nV]-
             Cu3[1+nS*(i-1):nS*i,1:nV])
            ./(u3g[i+1]-u3g[i]);
    elseif i==nr;
      Cu3u3[1+nS*(i-1):nS*i,1:nV] =
            (Cu3[1+nS*(i-1):nS*(i),1:nV]-
             Cu3[1+nS*(i-2):nS*(i-1),1:nV])
            ./(u3g[i]-u3g[i-1]);
    else;
      Cu3u3[1+nS*(i-1):nS*i,1:nV] =
            (Cu3[1+nS*i:nS*(i+1),1:nV]-
             Cu3[1+nS*(i-1):nS*i,1:nV])
            ./(u3g[i+1]-u3g[i]);
      Cu3u3[1+nS*(i-1):nS*i,1:nV] =
           (Cu3u3[1+nS*(i-1):nS*i,1:nV]+
            (Cu3[1+nS*(i-1):nS*(i),1:nV]-
             Cu3[1+nS*(i-2):nS*(i-1),1:nV])
            ./(u3g[i]-u3g[i-1]))/2;
    endif;

   /* Cu1u3 and Cu2u3 */

    Cu1u3[1+nS*(i-1):nS*i-1,1:nV] =
            (Cu3[2+nS*(i-1):nS*i,1:nV]-
            Cu3[1+nS*(i-1):nS*i-1,1:nV])
            ./(u1g[2:nS]-u1g[1:nS-1]);
    Cu1u3[1+nS*(i-1):nS*i,1:nV] =
            ((Cu1u3[1+nS*(i-1),1:nV]|
             Cu1u3[1+nS*(i-1):nS*i-1,1:nV]) +
            (Cu1u3[1+nS*(i-1):nS*i-1,1:nV]|
             Cu1u3[nS*i-1,1:nV]))/2;

    Cu2u3[1+nS*(i-1):nS*i,1:nV-1] =
         (Cu3[1+nS*(i-1):nS*i,2:nV]-
          Cu3[1+nS*(i-1):nS*i,1:nV-1])
          ./((u2g[2:nV]-u2g[1:nV-1])');
    Cu2u3[1+nS*(i-1):nS*i,1:nV] =
         ((Cu2u3[1+nS*(i-1):nS*i,1]~
           Cu2u3[1+nS*(i-1):nS*i,1:nV-1]) +
          (Cu2u3[1+nS*(i-1):nS*i,1:nV-1]~
           Cu2u3[1+nS*(i-1):nS*i,nV-1]))/2;

     i = i+1;
   endo;

   /* PDE Definition */
   i=1;
   do while i <= nr;
    uC[1+nS*(i-1):nS*i,1:nV] =
      uC[1+nS*(i-1):nS*i,1:nV]
      + dt*(
         - rg[i]*uC[1+nS*(i-1):nS*i,1:nV]
         + rg[i]*Sg.*
           (Cu1[1+nS*(i-1):nS*i,1:nV].*deriv1(u1g,Sg,1))
         + 1/2*(Sg^2).*Vg'.*
           (Cu1[1+nS*(i-1):nS*i,1:nV].*deriv2(u1g,Sg,1)+
            Cu1u1[1+nS*(i-1):nS*i,1:nV].*(deriv1(u1g,Sg,1))^2)
         + (muv(Vg')-(mus(Vg',rg[i])-rg[i])*_betas).*
           (Cu2[1+nS*(i-1):nS*i,1:nV].*deriv1(u2g',Vg',2))
         + 1/2*(sigvz(Vg')^2+sigvv(Vg')^2).*
           (Cu2[1+nS*(i-1):nS*i,1:nV].*deriv2(u2g',Vg',2)+
            Cu2u2[1+nS*(i-1):nS*i,1:nV].*(deriv1(u2g',Vg',2))^2)
         + (mur(rg[i])-(mus(Vg',rg[i])-rg[i]).*sigrz(rg[i])./sqrt(Vg')).*
           (Cu3[1+nS*(i-1):nS*i,1:nV].*deriv1(u3g[i],rg[i],3))
         + 1/2*(sigrz(rg[i])^2+sigrv(rg[i])^2+sigrr(rg[i])^2).*
           (Cu3[1+nS*(i-1):nS*i,1:nV].*deriv2(u3g[i],rg[i],3)+
            Cu3u3[1+nS*(i-1):nS*i,1:nV].*(deriv1(u3g[i],rg[i],3))^2)
         + _betas*Sg.*Vg'.*
           (Cu1u2[1+nS*(i-1):nS*i,1:nV].*deriv1(u1g,Sg,1).*
             deriv1(u2g',Vg',2))
         + Sg.*sqrt(Vg').*sigrz(rg[i]).*
           (Cu1u3[1+nS*(i-1):nS*i,1:nV].*deriv1(u1g,Sg,1).*
             deriv1(u3g[i],rg[i],3))
         + (sigvz(Vg').*sigrz(rg[i])+sigvv(Vg').*sigrv(rg[i])).*
           (Cu2u3[1+nS*(i-1):nS*i,1:nV].*deriv1(u2g',Vg',2).*
             deriv1(u3g[i],rg[i],3))
         + sqrt(A^2-((mus(Vg',rg[i])-rg[i])^2)./Vg').*sqrt(
           (sigvv(Vg').*
            (Cu2[1+nS*(i-1):nS*i,1:nV].*deriv1(u2g',Vg',2))
            + sigrv(rg[i]).*
           (Cu3[1+nS*(i-1):nS*i,1:nV].*deriv1(u3g[i],rg[i],3)) )^2
            + sigrr(rg[i])^2.*
           (Cu3[1+nS*(i-1):nS*i,1:nV].*deriv1(u3g[i],rg[i],3))^2 )
             );


    /* Boundary Conditions */

    i = i+1;
   endo;


    /* Boundary Conditions */


/*====> Computing the Lower Bound */

    /* find derivatives Cu1 Cu2 Cu3 Cu1u1 Cu2u2 Cu3u3 Cu1u2 Cu1u3 Cu2u3
       by finite differences */
    /* take 1/2 of forward and backward differences */

   i=1;
   do while i <= nr;

    /* Cu1 and Cu1u1 */

    Cu1[1+nS*(i-1):nS*i-1,1:nV] =
           (dC[2+nS*(i-1):nS*i,1:nV]-
           dC[1+nS*(i-1):nS*i-1,1:nV])
           ./(u1g[2:nS]-u1g[1:nS-1]);
    Cu1[1+nS*(i-1):nS*i,1:nV] =
          ((Cu1[1+nS*(i-1),1:nV]|
            Cu1[1+nS*(i-1):nS*i-1,1:nV]) +
           (Cu1[1+nS*(i-1):nS*i-1,1:nV]|
            Cu1[nS*i-1,1:nV]))/2;

    Cu1u1[1+nS*(i-1):nS*i-1,1:nV] =
            (Cu1[2+nS*(i-1):nS*i,1:nV]-
            Cu1[1+nS*(i-1):nS*i-1,1:nV])
            ./(u1g[2:nS]-u1g[1:nS-1]);
    Cu1u1[1+nS*(i-1):nS*i,1:nV] =
            ((Cu1u1[1+nS*(i-1),1:nV]|
             Cu1u1[1+nS*(i-1):nS*i-1,1:nV]) +
            (Cu1u1[1+nS*(i-1):nS*i-1,1:nV]|
             Cu1u1[nS*i-1,1:nV]))/2;

    j = 1;
    do while j <= wind;
        Cu1u1[1+nS*(i-1):nS*i,1:nV] =
              ((Cu1u1[1+nS*(i-1),1:nV]|
                Cu1u1[1+nS*(i-1):nS*i-1,1:nV])+
               Cu1u1[1+nS*(i-1):nS*i,1:nV]
               +(Cu1u1[2+nS*(i-1):nS*i,1:nV]|
               Cu1u1[nS*i,1:nV]))/3;
        j = j+1;
    endo;

    /* Cu2 and Cu2u2 */

    Cu2[1+nS*(i-1):nS*i,1:nV-1] =
         (dC[1+nS*(i-1):nS*i,2:nV]
          -dC[1+nS*(i-1):nS*i,1:nV-1])
          ./((u2g[2:nV]-u2g[1:nV-1])');
    Cu2[1+nS*(i-1):nS*i,1:nV] =
          ((Cu2[1+nS*(i-1):nS*i,1]~
            Cu2[1+nS*(i-1):nS*i,1:nV-1]) +
           (Cu2[1+nS*(i-1):nS*i,1:nV-1]~
            Cu2[1+nS*(i-1):nS*i,nV-1]))/2;

    Cu2u2[1+nS*(i-1):nS*i,1:nV-1] =
         (Cu2[1+nS*(i-1):nS*i,2:nV]-
          Cu2[1+nS*(i-1):nS*i,1:nV-1])
          ./((u2g[2:nV]-u2g[1:nV-1])');
    Cu2u2[1+nS*(i-1):nS*i,1:nV] =
         ((Cu2u2[1+nS*(i-1):nS*i,1]~
           Cu2u2[1+nS*(i-1):nS*i,1:nV-1]) +
          (Cu2u2[1+nS*(i-1):nS*i,1:nV-1]~
           Cu2u2[1+nS*(i-1):nS*i,nV-1]))/2;

    /* Cu1u2 */

    Cu1u2[1+nS*(i-1):nS*i,1:nV-1] =
         (Cu1[1+nS*(i-1):nS*i,2:nV]-
          Cu1[1+nS*(i-1):nS*i,1:nV-1])
          ./((u2g[2:nV]-u2g[1:nV-1])');
    Cu1u2[1+nS*(i-1):nS*i,1:nV] =
         ((Cu1u2[1+nS*(i-1):nS*i,1]~
           Cu1u2[1+nS*(i-1):nS*i,1:nV-1]) +
          (Cu1u2[1+nS*(i-1):nS*i,1:nV-1]~
           Cu1u2[1+nS*(i-1):nS*i,nV-1]))/2;

    /* Cu3 */

    if i==1;
      Cu3[1+nS*(i-1):nS*i,1:nV] =
            (dC[1+nS*i:nS*(i+1),1:nV]-
             dC[1+nS*(i-1):nS*i,1:nV])
            ./(u3g[i+1]-u3g[i]);
    elseif i==nr;
      Cu3[1+nS*(i-1):nS*i,1:nV] =
            (dC[1+nS*(i-1):nS*(i),1:nV]-
             dC[1+nS*(i-2):nS*(i-1),1:nV])
            ./(u3g[i]-u3g[i-1]);
    else;
      Cu3[1+nS*(i-1):nS*i,1:nV] =
            (dC[1+nS*i:nS*(i+1),1:nV]-
             dC[1+nS*(i-1):nS*i,1:nV])
            ./(u3g[i+1]-u3g[i]);
      Cu3[1+nS*(i-1):nS*i,1:nV] =
           (Cu3[1+nS*(i-1):nS*i,1:nV]+
            (dC[1+nS*(i-1):nS*(i),1:nV]-
             dC[1+nS*(i-2):nS*(i-1),1:nV])
            ./(u3g[i]-u3g[i-1]))/2;
    endif;

    i = i+1;
  endo;



   i=1;
   do while i <= nr;

    /* Cu3u3  */

     if i==1;
      Cu3u3[1+nS*(i-1):nS*i,1:nV] =
            (Cu3[1+nS*i:nS*(i+1),1:nV]-
             Cu3[1+nS*(i-1):nS*i,1:nV])
            ./(u3g[i+1]-u3g[i]);
    elseif i==nr;
      Cu3u3[1+nS*(i-1):nS*i,1:nV] =
            (Cu3[1+nS*(i-1):nS*(i),1:nV]-
             Cu3[1+nS*(i-2):nS*(i-1),1:nV])
            ./(u3g[i]-u3g[i-1]);
    else;
      Cu3u3[1+nS*(i-1):nS*i,1:nV] =
            (Cu3[1+nS*i:nS*(i+1),1:nV]-
             Cu3[1+nS*(i-1):nS*i,1:nV])
            ./(u3g[i+1]-u3g[i]);
      Cu3u3[1+nS*(i-1):nS*i,1:nV] =
           (Cu3u3[1+nS*(i-1):nS*i,1:nV]+
            (Cu3[1+nS*(i-1):nS*(i),1:nV]-
             Cu3[1+nS*(i-2):nS*(i-1),1:nV])
            ./(u3g[i]-u3g[i-1]))/2;
    endif;

   /* Cu1u3 and Cu2u3 */

    Cu1u3[1+nS*(i-1):nS*i-1,1:nV] =
            (Cu3[2+nS*(i-1):nS*i,1:nV]-
            Cu3[1+nS*(i-1):nS*i-1,1:nV])
            ./(u1g[2:nS]-u1g[1:nS-1]);
    Cu1u3[1+nS*(i-1):nS*i,1:nV] =
            ((Cu1u3[1+nS*(i-1),1:nV]|
             Cu1u3[1+nS*(i-1):nS*i-1,1:nV]) +
            (Cu1u3[1+nS*(i-1):nS*i-1,1:nV]|
             Cu1u3[nS*i-1,1:nV]))/2;

    Cu2u3[1+nS*(i-1):nS*i,1:nV-1] =
         (Cu3[1+nS*(i-1):nS*i,2:nV]-
          Cu3[1+nS*(i-1):nS*i,1:nV-1])
          ./((u2g[2:nV]-u2g[1:nV-1])');
    Cu2u3[1+nS*(i-1):nS*i,1:nV] =
         ((Cu2u3[1+nS*(i-1):nS*i,1]~
           Cu2u3[1+nS*(i-1):nS*i,1:nV-1]) +
          (Cu2u3[1+nS*(i-1):nS*i,1:nV-1]~
           Cu2u3[1+nS*(i-1):nS*i,nV-1]))/2;

     i = i+1;
   endo;

   /* PDE Definition */
   i=1;
   do while i <= nr;
    dC[1+nS*(i-1):nS*i,1:nV] =
      dC[1+nS*(i-1):nS*i,1:nV]
      + dt*(
         - rg[i]*dC[1+nS*(i-1):nS*i,1:nV]
         + rg[i]*Sg.*
           (Cu1[1+nS*(i-1):nS*i,1:nV].*deriv1(u1g,Sg,1))
         + 1/2*(Sg^2).*Vg'.*
           (Cu1[1+nS*(i-1):nS*i,1:nV].*deriv2(u1g,Sg,1)+
            Cu1u1[1+nS*(i-1):nS*i,1:nV].*(deriv1(u1g,Sg,1))^2)
         + (muv(Vg')-(mus(Vg',rg[i])-rg[i])*_betas).*
           (Cu2[1+nS*(i-1):nS*i,1:nV].*deriv1(u2g',Vg',2))
         + 1/2*(sigvz(Vg')^2+sigvv(Vg')^2).*
           (Cu2[1+nS*(i-1):nS*i,1:nV].*deriv2(u2g',Vg',2)+
            Cu2u2[1+nS*(i-1):nS*i,1:nV].*(deriv1(u2g',Vg',2))^2)
         + (mur(rg[i])-(mus(Vg',rg[i])-rg[i]).*sigrz(rg[i])./sqrt(Vg')).*
           (Cu3[1+nS*(i-1):nS*i,1:nV].*deriv1(u3g[i],rg[i],3))
         + 1/2*(sigrz(rg[i])^2+sigrv(rg[i])^2+sigrr(rg[i])^2).*
           (Cu3[1+nS*(i-1):nS*i,1:nV].*deriv2(u3g[i],rg[i],3)+
            Cu3u3[1+nS*(i-1):nS*i,1:nV].*(deriv1(u3g[i],rg[i],3))^2)
         + _betas*Sg.*Vg'.*
           (Cu1u2[1+nS*(i-1):nS*i,1:nV].*deriv1(u1g,Sg,1).*
             deriv1(u2g',Vg',2))
         + Sg.*sqrt(Vg').*sigrz(rg[i]).*
           (Cu1u3[1+nS*(i-1):nS*i,1:nV].*deriv1(u1g,Sg,1).*
             deriv1(u3g[i],rg[i],3))
         + (sigvz(Vg').*sigrz(rg[i])+sigvv(Vg').*sigrv(rg[i])).*
           (Cu2u3[1+nS*(i-1):nS*i,1:nV].*deriv1(u2g',Vg',2).*
             deriv1(u3g[i],rg[i],3))
         - sqrt(A^2-((mus(Vg',rg[i])-rg[i])^2)./Vg').*sqrt(
           (sigvv(Vg').*
            (Cu2[1+nS*(i-1):nS*i,1:nV].*deriv1(u2g',Vg',2))
            + sigrv(rg[i]).*
           (Cu3[1+nS*(i-1):nS*i,1:nV].*deriv1(u3g[i],rg[i],3)) )^2
            + sigrr(rg[i])^2.*
           (Cu3[1+nS*(i-1):nS*i,1:nV].*deriv1(u3g[i],rg[i],3))^2 )
             );


    /* Boundary Conditions */

    i = i+1;
   endo;


    /* Boundary Conditions */

/* ====> Computing the Stochastic Volatility-Cte Interest rate-Cte Risk Price */

    /* find derivatives Cu1 Cu2 Cu1u1 Cu2u2 Cu1u2
       by finite differences */
    /* take 1/2 of forward and backward differences */

   i=1;

    /* Cu1 and Cu1u1 */

    Cu1[1+nS*(i-1):nS*i-1,1:nV] =
           (SVC[2:nS,1:nV]-
           SVC[1:nS-1,1:nV])
           ./(u1g[2:nS]-u1g[1:nS-1]);
    Cu1[1+nS*(i-1):nS*i,1:nV] =
          ((Cu1[1+nS*(i-1),1:nV]|
            Cu1[1+nS*(i-1):nS*i-1,1:nV]) +
           (Cu1[1+nS*(i-1):nS*i-1,1:nV]|
            Cu1[nS*i-1,1:nV]))/2;

    Cu1u1[1+nS*(i-1):nS*i-1,1:nV] =
            (Cu1[2+nS*(i-1):nS*i,1:nV]-
            Cu1[1+nS*(i-1):nS*i-1,1:nV])
            ./(u1g[2:nS]-u1g[1:nS-1]);
    Cu1u1[1+nS*(i-1):nS*i,1:nV] =
            ((Cu1u1[1+nS*(i-1),1:nV]|
             Cu1u1[1+nS*(i-1):nS*i-1,1:nV]) +
            (Cu1u1[1+nS*(i-1):nS*i-1,1:nV]|
             Cu1u1[nS*i-1,1:nV]))/2;

    j = 1;
    do while j <= wind;
        Cu1u1[1+nS*(i-1):nS*i,1:nV] =
              ((Cu1u1[1+nS*(i-1),1:nV]|
                Cu1u1[1+nS*(i-1):nS*i-1,1:nV])+
               Cu1u1[1+nS*(i-1):nS*i,1:nV]
               +(Cu1u1[2+nS*(i-1):nS*i,1:nV]|
               Cu1u1[nS*i,1:nV]))/3;
        j = j+1;
    endo;

    /* Cu2 and Cu2u2 */

    Cu2[1+nS*(i-1):nS*i,1:nV-1] =
         (SVC[1:nS,2:nV]
          -SVC[1:nS,1:nV-1])
          ./((u2g[2:nV]-u2g[1:nV-1])');
    Cu2[1+nS*(i-1):nS*i,1:nV] =
          ((Cu2[1+nS*(i-1):nS*i,1]~
            Cu2[1+nS*(i-1):nS*i,1:nV-1]) +
           (Cu2[1+nS*(i-1):nS*i,1:nV-1]~
            Cu2[1+nS*(i-1):nS*i,nV-1]))/2;

    Cu2u2[1+nS*(i-1):nS*i,1:nV-1] =
         (Cu2[1+nS*(i-1):nS*i,2:nV]-
          Cu2[1+nS*(i-1):nS*i,1:nV-1])
          ./((u2g[2:nV]-u2g[1:nV-1])');
    Cu2u2[1+nS*(i-1):nS*i,1:nV] =
         ((Cu2u2[1+nS*(i-1):nS*i,1]~
           Cu2u2[1+nS*(i-1):nS*i,1:nV-1]) +
          (Cu2u2[1+nS*(i-1):nS*i,1:nV-1]~
           Cu2u2[1+nS*(i-1):nS*i,nV-1]))/2;

    /* Cu1u2 */

    Cu1u2[1+nS*(i-1):nS*i,1:nV-1] =
         (Cu1[1+nS*(i-1):nS*i,2:nV]-
          Cu1[1+nS*(i-1):nS*i,1:nV-1])
          ./((u2g[2:nV]-u2g[1:nV-1])');
    Cu1u2[1+nS*(i-1):nS*i,1:nV] =
         ((Cu1u2[1+nS*(i-1):nS*i,1]~
           Cu1u2[1+nS*(i-1):nS*i,1:nV-1]) +
          (Cu1u2[1+nS*(i-1):nS*i,1:nV-1]~
           Cu1u2[1+nS*(i-1):nS*i,nV-1]))/2;


   /* PDE Definition */


    SVC = SVC + dt*(- _r0*SVC
         + _r0*Sg.*
           (Cu1[1+nS*(i-1):nS*i,1:nV].*deriv1(u1g,Sg,1))
         + 1/2*(Sg^2).*Vg'.*
           (Cu1[1+nS*(i-1):nS*i,1:nV].*deriv2(u1g,Sg,1)+
            Cu1u1[1+nS*(i-1):nS*i,1:nV].*(deriv1(u1g,Sg,1))^2)
         + (muv(Vg')-sqrt(Vg')*lambda).*
           (Cu2[1+nS*(i-1):nS*i,1:nV].*deriv1(u2g',Vg',2))
         + 1/2*(sigvz(Vg')^2+sigvv(Vg')^2).*
           (Cu2[1+nS*(i-1):nS*i,1:nV].*deriv2(u2g',Vg',2)+
            Cu2u2[1+nS*(i-1):nS*i,1:nV].*(deriv1(u2g',Vg',2))^2)
         + _betas*Sg.*Vg'.*
           (Cu1u2[1+nS*(i-1):nS*i,1:nV].*deriv1(u1g,Sg,1).*
             deriv1(u2g',Vg',2))
             );

/* ====> Computing the Black-Scholes Price */

    /* find derivatives Cu1 Cu1u1
       by finite differences */
    /* take 1/2 of forward and backward differences */

   i=1;

    /* Cu1 and Cu1u1 */

    Cu1[1+nS*(i-1):nS*i-1,1] =
           (BSC[2:nS,1]-
           BSC[1:nS-1,1])
           ./(u1g[2:nS]-u1g[1:nS-1]);
    Cu1[1+nS*(i-1):nS*i,1] =
          ((Cu1[1+nS*(i-1),1]|
            Cu1[1+nS*(i-1):nS*i-1,1]) +
           (Cu1[1+nS*(i-1):nS*i-1,1]|
            Cu1[nS*i-1,1]))/2;

    Cu1u1[1+nS*(i-1):nS*i-1,1] =
            (Cu1[2+nS*(i-1):nS*i,1]-
            Cu1[1+nS*(i-1):nS*i-1,1])
            ./(u1g[2:nS]-u1g[1:nS-1]);
    Cu1u1[1+nS*(i-1):nS*i,1] =
            ((Cu1u1[1+nS*(i-1),1]|
             Cu1u1[1+nS*(i-1):nS*i-1,1]) +
            (Cu1u1[1+nS*(i-1):nS*i-1,1]|
             Cu1u1[nS*i-1,1]))/2;

    j = 1;
    do while j <= wind;
        Cu1u1[1+nS*(i-1):nS*i,1] =
              ((Cu1u1[1+nS*(i-1),1]|
                Cu1u1[1+nS*(i-1):nS*i-1,1])+
               Cu1u1[1+nS*(i-1):nS*i,1]
               +(Cu1u1[2+nS*(i-1):nS*i,1]|
               Cu1u1[nS*i,1]))/3;
        j = j+1;
    endo;


   /* PDE Definition */


    BSC = BSC + dt*(- _r0*BSC
         + _r0*Sg.*
           (Cu1[1+nS*(i-1):nS*i,1].*deriv1(u1g,Sg,1))
         + 1/2*iV*(Sg^2).*
           (Cu1[1+nS*(i-1):nS*i,1].*deriv2(u1g,Sg,1)+
            Cu1u1[1+nS*(i-1):nS*i,1].*(deriv1(u1g,Sg,1))^2)
         );


    if (bigp==1) and abs(trunc(t/0.1) - (t/0.1)) < 0.0000001;
/*        d1 = (ln(sg)-ln(x) + (_r0+(iV)/2)*t )
              /((iV*t)^0.5);
        d2 = d1 - (iV*t)^0.5;
        bsc = sg.*cdfn(d1) - X*exp(-(_r0)*t)*cdfn(d2); */

    _pstype = 1|9|12;
    _pltype = 3|6|6;
    _plctrl = 0|2|0;

    _plegstr ="Upper bound\000Lower bound\000volatility ";
    _plegctl = 1;
        xy(sg,uC[1+nS*(r-1):nS*r,iVi]-bsc[1:nS,1]~
           dC[1+nS*(r-1):nS*r,iVi]-bsc[1:nS,1]~
           SVC[1:nS,iVi]-bsc[1:nS,1]);
        ""; "TIME " t*1000;
     endif;



   t = t+dt;
endo;

output off;

title("Price differences with stochastic V and r");
xlabel("stock price");
ylabel("price difference");
termprc = maxc( (sg-x)'|zeros(1,nS) );  /* final price */;
lowbnd = maxc( (sg-x*exp(-_r0*maxt))'|zeros(1,nS) );   /* lower arb bound
d1 = (ln(sg)-ln(x) + (_r0+(iV)/2)*maxt )
        /((iV*maxt)^0.5);
d2 = d1 - (iV*maxt)^0.5;
bsc = sg.*cdfn(d1) - X*exp(-(_r0)*maxt)*cdfn(d2); */
    _pstype = 1|9|12|6|3|1|1;
    _pltype = 3|6|6|3|3|6|6;
    _plctrl = 0|2|2|2|2|0|0;

    _plegstr ="Upper bound\000Lower bound\000volatility ";
    _plegctl = 1;
    xy(sg,uC[1+nS*(r-1):nS*r,iVi]-bsc[1:nS,1]~
          dC[1+nS*(r-1):nS*r,iVi]-bsc[1:nS,1]~
          SVC[1:nS,iVi]-bsc[1:nS,1]);

uS = indexcat(Sg,(400-0.001)|(400+0.001));    /* change to index */
title("Price bounds with stochastic V and r");
xlabel("stock prices");
ylabel("price diferences");
    _pstype = 1|9|12|6|3|1|1;
    _pltype = 3|6|6|3|3|6|6;
    _plctrl = 0|2|2|2|2|0|0;

    _plegstr ="Upper bound\000Lower bound\0000 prices of risk ";
    _plegctl = 1;
    _pstype = 1|8|9@|1|1@|1;
    _pltype = 3|6|6@|6|6@|6;
    _plctrl = 0|0|10@|2|2@|0;

    _plegctl = 1;
    _pnum = 1|2;
    _pdate = 0;
    _pframe = 0|1;

    _protate = 1;
    begwind;
    makewind(9,4,0,3.2,1);
    _pnumht=0.2;

    xy(sg[1:uS],uC[1+nS*(r-1):uS+nS*(r-1),iVi]-bsc[1:uS,1]~
         dC[1+nS*(r-1):uS+nS*(r-1),iVi]-bsc[1:uS,1]~
        SVC[1:uS,iVi]-bsc[1:uS,1]);
    endwind;

@----------------------- Changes of Variable ------------------------------@

proc StoU1(Svalue,ii);  /* change of variable S->U1 */

    local uvalue;
    uvalue=zeros(rows(Svalue),1);
    if chvar[ii]==0; uvalue=Svalue; endif;
    if chvar[ii]==1; uvalue=c[ii]*Svalue./(1+c[ii]*Svalue); endif;
    if chvar[ii]==2; uvalue=(Svalue./(1+Svalue))^(1/c[ii]); endif;

    retp(uvalue);

endp;

proc U1toS(uvalue,ii);  /* inverted change of variable U1->S */

    local Svalue;
    Svalue=zeros(rows(uvalue),1);
    if chvar[ii]==0; Svalue=uvalue; endif;
    if chvar[ii]==1; Svalue=(1/c[ii])*(uvalue./(1-uvalue)); endif;
    if chvar[ii]==2;
       Svalue=uvalue^(c[ii]);
       Svalue=((Svalue)./(1-Svalue));
    endif;

    retp(Svalue);

endp;


proc deriv1(U,S,ii);

    local d;
    d=zeros(rows(S),1);

    if chvar[ii]==0; d=ones(rows(S),1); endif;
    if chvar[ii]==1;
      d=c[ii]*ones(rows(S),1)./((1+c[ii]*S).*(1+c[ii]*S));
    endif;
    if chvar[ii]==2;
      d=(ones(rows(S),1)./(1+S))^2;
      d=(1/c[ii])*d.*(U^(1-c[ii]));
    endif;

    retp(d);

endp;



proc deriv2(U,S,ii);

    local d2, d1;
    d2=zeros(rows(S),1);

    if chvar[ii]==1;
     d2=-2*c[ii]*c[ii]*ones(rows(S),1)./((1+c[ii]*S).*(1+c[ii]*S).*(1+c[ii]*S));
    endif;

    if chvar==2;
     d2=ones(rows(S),1)./(1+S);
     d2=(1/c[ii])*(d2^3).*(-2*(U^(1-c[ii]))+
          ((1-c[ii])/c[ii])*d2.*(U^(1-2*c[ii])));
    endif;

    retp(d2);

endp;

/* ------------------------------------------------------------------------- */
/*                             SEQB.G                                        */
/*
    An alternative to seqa.  specify instead min, max, and a step size.  If
    max is not min + integer number * step size, the next smaller number is
    used.

    Usage:

    x = seqb(min,max,step);

*/
/* ------------------------------------------------------------------------- */
proc seqb(min,max,step);
    local nsteps, x;
     if max < min ; "seqb: asking for max < min"; endif;
    nsteps = trunc((max-min)/step+1.000000000001);
    x = seqa(min,step,nsteps);
    retp(x);
endp;

proc(1) = sectohms(sec);  /* converts seconds to hrs, mins,secs */
    local hrs,mins,secs;
    hrs = trunc(sec/3600);
    mins = trunc(sec/60)-60*hrs;
    secs = sec-3600*hrs-60*mins;
    retp(hrs|mins|secs);
endp;