new; cls;
/* Deterministic Trinomial European Call & Put Options (checking by Trinomial Tree Analyses) */
/* Notice: This version will not go into overflow in tsampler2 itself. */
/*         However, its probability pt2 will become too small to zero. */
/*         This way of thinking can not be applied with probability.   */
/*         Use stochastic vesion of it instead if the number of time   */
/*         steps is large and then each probablity becomes zero.       */
S0=50;
K=50;
r=0.10;
sig=0.40;
T=5/12;

nn=10;
call eurocall3(S0,K,r,sig,T,nn);
call europut3(S0,K,r,sig,T,nn);

nn=10;
print "by this method(both results should be the same):";
print "C=" Tcall(S0,K,r,sig,T,nn);
print "P=" Tput(S0,K,r,sig,T,nn);


/*
**  teuro.txt - Deterministic Trinomial European Call & Put Options.
** (C) Copyright 2005 Yosuke Amijima. All Rights Reserved.
**
**  Purpose:    Calculates european option prices by deterministic trinomial simulation.
**
**  Format:     C=Tcall(S0,K,r,sig,T,nn);
**              P=Tput(S0,K,r,sig,T,nn);
**
**  Input:      S0       scalar,  initial value
**
**              K        scalar,  strike price
**             
**              r        scalar,  risk-free interest rate
**
**              sig      scalar,  volatility
**
**              T        scalar,  maturity
**
**              nn       scalar,  number of time steps
**
**
**  Output:     C        scalar,  call option price
**              P        scalar,  put  option price
**
**  Notice:     The result should be exactly the same as that of trinomial tree analysis.
*/


proc Tcall(S0,K,r,sig,T,nn);
   local x,S,ST,pp,i,C;
   x=zeros(nn,1);
   S=tsampler2(S0,sig,T,x);
   ST=S[nn+1];
   pp=pt2(S0,r,sig,T,x);
   C=exp(-r*T)*maxc((ST-K)|0)*pp;
   do while not x==2;
      x=rec3index(x);
      S=tsampler2(S0,sig,T,x);
      ST=S[nn+1];
      pp=pt2(S0,r,sig,T,x);
      C=C+exp(-r*T)*maxc((ST-K)|0)*pp;
   endo;
   retp(C);
endp;

proc Tput(S0,K,r,sig,T,nn);
   local x,S,ST,pp,i,P;
   x=zeros(nn,1);
   S=tsampler2(S0,sig,T,x);
   ST=S[nn+1];
   pp=pt2(S0,r,sig,T,x);
   P=exp(-r*T)*maxc((K-ST)|0)*pp;
   do while not x==2;
      x=rec3index(x);
      S=tsampler2(S0,sig,T,x);
      ST=S[nn+1];
      pp=pt2(S0,r,sig,T,x);
      P=P+exp(-r*T)*maxc((K-ST)|0)*pp;
   endo;
   retp(P);
endp;


/*
**  tsampler2.txt - Deterministic Trinomial Path Sampler.
** (C) Copyright 2005 Yosuke Amijima. All Rights Reserved.
**
**  Purpose:    Gets a deterministic path on trinomial tree for given ternary index vector x
**              without overflow problem.
**
**  Format:     S=tsampler2(S0,sig,T,x)
**
**  Input:      S0       scalar,    initial value
**
**              sig      scalar,    volatility
**
**              T        scalar,    maturity
**
**              x        vector,    nn x 1 of ternary index vector
**
**
**  Output:     S        vector,    (nn+1) x 1 of resulting values including S0
**
*/
proc tsampler2(S0,sig,T,x);
   local nn,delt,u,S;
   nn=rows(x);
   x=x-1*(x.==0)-1*(x.==1)-1*(x.==2);
   delt=T/nn;
   u=exp(sig*sqrt(3*delt));
   S=S0*u^cumsumc(x);
   S=S0|S;
   retp(S);
endp;

/*
**  tsampler2.txt - Probability of the Deterministic Trinomial Path.
** (C) Copyright 2005 Yosuke Amijima. All Rights Reserved.
**
**  Purpose:    Gets the probability of the path for given ternary index vector x
**              without overflow problem.
**
**  Format:     p=pt2(S0,r,sig,T,nn,i);
**
**  Input:      S0       scalar,    initial value
**
**              r        scalar,    risk-free interest rate
**
**              sig      scalar,    volatility
**
**              T        scalar,    maturity
**
**              x        vector,    nn x 1 of ternary index vector
**
**
**  Output:     p        scalar,    probability of the path
**
*/
proc pt2(S0,r,sig,T,x);
   local nn,delt,pd,pm,pu,p;
   nn=rows(x);
   delt=T/nn;
   pd=-sqrt(delt/(12*sig^2))*(r-1/2*sig^2)+1/6;
   pm=2/3;
   pu=sqrt(delt/(12*sig^2))*(r-1/2*sig^2)+1/6;
   p=pd*(x.==0)+pm*(x.==1)+pu*(x.==2);
   p=prodc(p);
   retp(p);
endp;

/*
**  rec3index.txt - Recursive Ternary Dimension Indexing.
** (C) Copyright 2005 Yosuke Amijima. All Rights Reserved.
**
**  Purpose:    Gets ternary index numbers for given dimension recursively 
**              without overflow problem.
**
**  Format:     x=rec3index(x);
**
**  Input:      x        vector,    dim x 1 of current ternary index vector
**
**
**  Output:     x        vector,    dim x 1 of next ternary index vector
**
**  Notice:     You could start from x={0,0,0,0,0}; for example.
*/
proc rec3index(x);
   local dim,j;
   dim=rows(x);
   x[dim]=x[dim]+1;
   j=dim;
   do until j==0;
      if x[j]==3;
         x[j]=0;
         x[j-1]=x[j-1]+1;
      endif;
      j=j-1;
   endo;
   retp(x);
endp;





/* by comparison */

proc eurocall3(S0,K,r,sig,T,nn);
   local delt,u,d,pd,pm,pu,m1,m2,i,j;
   delt=T/nn;
   u=exp(sig*sqrt(3*delt));
   d=exp(-sig*sqrt(3*delt));
   pd=1/6-sqrt(delt/(12*sig^2))*(r-sig^2/2);
   pm=2/3;
   pu=1/6+sqrt(delt/(12*sig^2))*(r-sig^2/2);
   print "== European Call Option ==";
   print "delt=" delt delt*365 "days";
   print "   u=" u;
   print "   d=" d;
   print "  pu=" pu;
   print "  pm=" pm;
   print "  pd=" pd;
   print;
 /* Stock Price Tree */
   m1=miss(zeros(2*nn+1,nn+1),0);
   j=1;
   do while j<=cols(m1);
      i=1;
      do while i<=2*j-1;
         m1[nn+1-j+i,j]=S0*u^(j-i);
         i=i+1;
      endo;
      j=j+1;
   endo;
   print "Trinomial Stock Price Tree:";
   print/rz seqa(0,1,nn+1)';
   print/rz m1;
   print "Node Time:";
   print delt*seqa(0,1,nn+1)';
   print;
/* Option Price Tree & Exercise Tree */
   m2=miss(zeros(2*nn+1,nn+1),0);
   /* the last column */
   i=1;
   do while i<=2*nn+1;
      m2[i,nn+1]=maxc(m1[i,nn+1]-K|0);    /* call setting */
      i=i+1;
   endo;
   /* backward */
   j=nn;
   do while j>=1;
      i=1;
      do while i<=2*j-1;
         m2[nn+1-j+i,j]=(pu*m2[nn+1-j+i-1,j+1]+pm*m2[nn+1-j+i,j+1]+pd*m2[nn+1-j+i+1,j+1])*exp(-r*delt); 
         i=i+1;
      endo;
      j=j-1;
   endo;  
   print "Option Price Tree:";
   print/rz seqa(0,1,nn+1)';
   print/rz m2;
   print/lz "C=" m2[nn+1,1];
   print;
   retp(m2[nn+1,1]);
endp;

proc europut3(S0,K,r,sig,T,nn);
   local delt,u,d,pd,pm,pu,m1,m2,i,j;
   delt=T/nn;
   u=exp(sig*sqrt(3*delt));
   d=exp(-sig*sqrt(3*delt));
   pd=1/6-sqrt(delt/(12*sig^2))*(r-sig^2/2);
   pm=2/3;
   pu=1/6+sqrt(delt/(12*sig^2))*(r-sig^2/2);
   print "== European Put Option ==";
   print "delt=" delt delt*365 "days";
   print "   u=" u;
   print "   d=" d;
   print "  pu=" pu;
   print "  pm=" pm;
   print "  pd=" pd;
   print;
 /* Stock Price Tree */
   m1=miss(zeros(2*nn+1,nn+1),0);
   j=1;
   do while j<=cols(m1);
      i=1;
      do while i<=2*j-1;
         m1[nn+1-j+i,j]=S0*u^(j-i);
         i=i+1;
      endo;
      j=j+1;
   endo;
   print "Trinomial Stock Price Tree:";
   print/rz seqa(0,1,nn+1)';
   print/rz m1;
   print "Node Time:";
   print delt*seqa(0,1,nn+1)';
   print;
/* Option Price Tree & Exercise Tree */
   m2=miss(zeros(2*nn+1,nn+1),0);
   /* the last column */
   i=1;
   do while i<=2*nn+1;
      m2[i,nn+1]=maxc(K-m1[i,nn+1]|0);    /* put setting */
      i=i+1;
   endo;
   /* backward */
   j=nn;
   do while j>=1;
      i=1;
      do while i<=2*j-1;
         m2[nn+1-j+i,j]=(pu*m2[nn+1-j+i-1,j+1]+pm*m2[nn+1-j+i,j+1]+pd*m2[nn+1-j+i+1,j+1])*exp(-r*delt); 
         i=i+1;
      endo;
      j=j-1;
   endo;  
   print "Option Price Tree:";
   print/rz seqa(0,1,nn+1)';
   print/rz m2;
   print/lz "P=" m2[nn+1,1];
   print;
   retp(m2[nn+1,1]);
endp;