new; cls;
/* Deterministic Binomial European Call & Put Options (checking by Binomial Tree Analyses) */
/* Notice: This version will not go into overflow in bsampler2 itself. */
/*         However, its probability pb2 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 eurocall(S0,K,r,sig,T,nn);
call europut(S0,K,r,sig,T,nn);

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


/*
**  beuro.txt - Deterministic Binomial European Call & Put Options.
** (C) Copyright 2005 Yosuke Amijima. All Rights Reserved.
**
**  Purpose:    Calculates european option prices by deterministic binomial simulation.
**
**  Format:     C=Bcall(S0,K,r,sig,T,nn);
**              P=Bput(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 binomial tree analysis.
*/

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

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

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

/*
**  bsampler2.txt - Probability of the Deterministic Binomial Path.
** (C) Copyright 2005 Yosuke Amijima. All Rights Reserved.
**
**  Purpose:    Gets the probability of the path for given binary index vector x
**              without overflow problem.
**
**  Format:     p=pb2(S0,r,sig,T,x);
**
**  Input:      S0       scalar,    initial value
**
**              r        scalar,    risk-free interest rate
**
**              sig      scalar,    volatility
**
**              T        scalar,    maturity
**
**              x        vector,    nn x 1 of binary index vector
**
**
**  Output:     p        scalar,    probability of the path
**
*/
proc pb2(S0,r,sig,T,x);
   local nn,delt,a,u,d,p,S;
   nn=rows(x);
   delt=T/nn;
   a=exp(r*delt);
   u=exp(sig*sqrt(delt));
   d=1/u;
   p=(a-d)/(u-d);
   p=p*(x.==1)+(1-p)*(x.==0);
   p=prodc(p);
   retp(p);
endp;


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






/* by comparison */

proc eurocall(S0,K,r,sig,T,nn);
   local delt,u,d,a,p,m1,m2,i,j;
   delt=T/nn;
   u=exp(sig*sqrt(delt));
   d=exp(-sig*sqrt(delt));
   a=exp(r*delt);
   if u==1 and d==1;
      p=1;                             
   else;
      p=(a-d)/(u-d);
   endif;
   print "== European Call Option ==";
   print "delt=" delt delt*365 "days";
   print "   u=" u;
   print "   d=" d;
   print "   p=" p;
   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-1-(i-1)/2)*d^((i-1)/2);
         i=i+2;
      endo;
      j=j+1;
   endo;
   print "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 */
   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+2;
   endo;
   /* backward */
   j=nn;
   do while j>=1;
      i=1;
      do while i<=2*j-1;
         m2[nn+1-j+i,j]=(p*m2[nn+1-j+i-1,j+1]+(1-p)*m2[nn+1-j+i+1,j+1])*exp(-r*delt); 
         i=i+2;
      endo;
      j=j-1;
   endo;   
   print "Option Price Tree:";
   print/rz seqa(0,1,nn+1)';
   print/rz m2;
   print;
   print/lz "C=" m2[nn+1,1];
   print;
   retp(m2[nn+1,1]);
endp;

proc europut(S0,K,r,sig,T,nn);
   local delt,u,d,a,p,m1,m2,i,j;
   delt=T/nn;
   u=exp(sig*sqrt(delt));
   d=exp(-sig*sqrt(delt));
   a=exp(r*delt);
   if u==1 and d==1;
      p=1;                             
   else;
      p=(a-d)/(u-d);
   endif;
   print "== European Put Option ==";
   print "delt=" delt delt*365 "days";
   print "   u=" u;
   print "   d=" d;
   print "   p=" p;
   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-1-(i-1)/2)*d^((i-1)/2);
         i=i+2;
      endo;
      j=j+1;
   endo;
   print "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 */
   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+2;
   endo;
   /* backward */
   j=nn;
   do while j>=1;
      i=1;
      do while i<=2*j-1;
         m2[nn+1-j+i,j]=(p*m2[nn+1-j+i-1,j+1]+(1-p)*m2[nn+1-j+i+1,j+1])*exp(-r*delt); 
         i=i+2;
      endo;
      j=j-1;
   endo;   
   print "Option Price Tree:";
   print/rz seqa(0,1,nn+1)';
   print/rz m2;
   print;
   print/lz "P=" m2[nn+1,1];
   print;
   retp(m2[nn+1,1]);
endp;