new; cls;
/* Stochastic Trinomial Barrier Call & Put Options    */
S0=50;
K=50;
r=0.10;
sig=0.40;
T=5/12;
nn=500;
vH=60;
vL=15; 
times=5000;

print/lz "     nn=" nn;
print/lz "  times=" times;
print;
print "double barriers(vL<=St<=vH):";
print/lz "     vH=" vH;
print/lz "     vL=" vL;
print "C=" PTbarriersC(S0,K,r,sig,T,nn,vH,vL,times);
print "P=" PTbarriersP(S0,K,r,sig,T,nn,vH,vL,times);
print;

vB=60;
print "up and in:";
print/lz "     vB=" vB;
print "C=" PTupandinC(S0,K,r,sig,T,nn,vB,times);
print "P=" PTupandinP(S0,K,r,sig,T,nn,vB,times);
print;
print "up and out:";
print/lz "     vB=" vB;
print "C=" PTupandoutC(S0,K,r,sig,T,nn,vB,times);
print "P=" PTupandoutP(S0,K,r,sig,T,nn,vB,times);
print;
print "down and in:";
print/lz "     vB=" vB;
print "C=" PTdownandinC(S0,K,r,sig,T,nn,vB,times);
print "P=" PTdownandinP(S0,K,r,sig,T,nn,vB,times);
print;

vB=15;
print "down and out:";
print/lz "     vB=" vB;
print "C=" PTdownandoutC(S0,K,r,sig,T,nn,vB,times);
print "P=" PTdownandoutP(S0,K,r,sig,T,nn,vB,times);


/*
**  ptbarrier.txt - Stochastic Trinomial Double Barriers Call & Put Options.
** (C) Copyright 2005 Yosuke Amijima. All Rights Reserved.
**
**  Purpose:    Calculates Barrier option prices by stochastic trinomial simulation.
**
**  Format:     C=PTbarriersC(S0,K,r,sig,T,nn,vH,vL,times);
**              P=PTbarriersP(S0,K,r,sig,T,nn,vH,vL,times);
**
**  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
**
**              vH, vL   scalar,  barriers ( effective when vH <= St <= vL )
**                                if it is oneside barrier, set 1e8 or 1e-8 to either.
**
**              times    scalar,  number of simulations
**
**
**  Output:     C        scalar,  call option price
**              P        scalar,  put  option price
**
*/

proc PTbarriersC(S0,K,r,sig,T,nn,vH,vL,times);
   local OV,i,S,C;
   OV=zeros(times,1);
   i=1;
   do while i<=times;
      S=ptsampler(S0,r,sig,T,nn);
      if S>=vL and S<=vH;          /* vector evaluation */
         OV[i]=S[nn+1]-K;   /* effective  */
      else;
         OV[i]=0;           /* value=0    */
      endif;
      i=i+1;
   endo;
   C=exp(-r*T)*meanc(maxc((OV~zeros(times,1))'));
   retp(C);
endp;

proc PTbarriersP(S0,K,r,sig,T,nn,vH,vL,times);
   local OV,i,S,P;
   OV=zeros(times,1);
   i=1;
   do while i<=times;
      S=ptsampler(S0,r,sig,T,nn);
      if S>=vL and S<=vH;          /* vector evaluation */
         OV[i]=S[nn+1]-K;   /* effective  */
      else;
         OV[i]=0;           /* value=0    */
      endif;
      i=i+1;
   endo;
   P=exp(-r*T)*meanc(maxc((OV~zeros(times,1))'));
   retp(P);
endp;


/*
**  ptbarrier.txt - Stochastic Trinomial Knock-in(out) Call & Put Options.
** (C) Copyright 2005 Yosuke Amijima. All Rights Reserved.
**
**  Purpose:    Calculates Knock-in and Knock-out option prices by stochastic trinomial simulation.
**
**  Format:     C=PTupandinC(S0,K,r,sig,T,nn,vB,times);     UP and IN
**              P=PTupandinP(S0,K,r,sig,T,nn,vB,times);
**
**              C=PTupandoutC(S0,K,r,sig,T,nn,vB,times);    UP and OUT
**              P=PTupandoutP(S0,K,r,sig,T,nn,vB,times);
**
**              C=PTdownandinC(S0,K,r,sig,T,nn,vB,times);   DOWN and IN
**              P=PTdownandinP(S0,K,r,sig,T,nn,vB,times);
**
**              C=PTdownandoutC(S0,K,r,sig,T,nn,vB,times);  DOWN and OUT
**              P=PTdownandoutP(S0,K,r,sig,T,nn,vB,times);
**
**
**  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
**
**              vB       scalar,  barrier
**
**              times    scalar,  number of simulations
**
**
**  Output:     C        scalar,  call option price
**              P        scalar,  put  option price
**
*/

proc PTupandinC(S0,K,r,sig,T,nn,vB,times);
   local OV,i,S,C;
   OV=zeros(times,1);
   i=1;
   do while i<=times;
      S=ptsampler(S0,r,sig,T,nn);
      if S<=vB;                  /* vector evaluation */
         OV[i]=0;          /* value=0    */
      else;
         OV[i]=S[nn+1]-K;  /* effective  */
      endif;
      i=i+1;
   endo;
   C=exp(-r*T)*meanc(maxc((OV~zeros(times,1))'));
   retp(C);
endp;

proc PTupandinP(S0,K,r,sig,T,nn,vB,times);
   local OV,i,S,P;
   OV=zeros(times,1);
   i=1;
   do while i<=times;
      S=ptsampler(S0,r,sig,T,nn);
      if S<=vB;                   /* vector evaluation */
         OV[i]=0;           /* value=0    */
      else;
         OV[i]=K-S[nn+1];   /* effective  */
      endif;
      i=i+1;
   endo;
   P=exp(-r*T)*meanc(maxc((OV~zeros(times,1))'));
   retp(P);
endp;

proc PTupandoutC(S0,K,r,sig,T,nn,vB,times);
   local OV,i,S,C;
   OV=zeros(times,1);
   i=1;
   do while i<=times;
      S=ptsampler(S0,r,sig,T,nn);
      if S<=vB;                   /* vector evaluation */
         OV[i]=S[nn+1]-K;  /* effective  */
      else;
         OV[i]=0;          /* value=0    */
      endif;
      i=i+1;
   endo;
   C=exp(-r*T)*meanc(maxc((OV~zeros(times,1))'));
   retp(C);
endp;

proc PTupandoutP(S0,K,r,sig,T,nn,vB,times);
   local OV,i,S,P;
   OV=zeros(times,1);
   i=1;
   do while i<=times;
      S=ptsampler(S0,r,sig,T,nn);
      if S<=vB;                   /* vector evaluation */
         OV[i]=K-S[nn+1];  /* effective  */
      else;
         OV[i]=0;          /* value=0    */
      endif;
      i=i+1;
   endo;
   P=exp(-r*T)*meanc(maxc((OV~zeros(times,1))'));
   retp(P);
endp;

proc PTdownandinC(S0,K,r,sig,T,nn,vB,times);
   local OV,i,S,C;
   OV=zeros(times,1);
   i=1;
   do while i<=times;
      S=ptsampler(S0,r,sig,T,nn);
      if S>=vB;                  /* vector evaluation */
         OV[i]=0;          /* value=0    */
      else;
         OV[i]=S[nn+1]-K;  /* effective  */
      endif;
      i=i+1;
   endo;
   C=exp(-r*T)*meanc(maxc((OV~zeros(times,1))'));
   retp(C);
endp;

proc PTdownandinP(S0,K,r,sig,T,nn,vB,times);
   local OV,i,S,P;
   OV=zeros(times,1);
   i=1;
   do while i<=times;
      S=ptsampler(S0,r,sig,T,nn);
      if S>=vB;                  /* vector evaluation */
         OV[i]=0;          /* value=0    */
      else;
         OV[i]=K-S[nn+1];  /* effective  */
      endif;
      i=i+1;
   endo;
   P=exp(-r*T)*meanc(maxc((OV~zeros(times,1))'));
   retp(P);
endp;

proc PTdownandoutC(S0,K,r,sig,T,nn,vB,times);
   local OV,i,S,C;
   OV=zeros(times,1);
   i=1;
   do while i<=times;
      S=ptsampler(S0,r,sig,T,nn);
      if S>=vB;                  /* vector evaluation */
         OV[i]=S[nn+1]-K;  /* effective  */
      else;
         OV[i]=0;          /* value=0    */
      endif;
      i=i+1;
   endo;
   C=exp(-r*T)*meanc(maxc((OV~zeros(times,1))'));
   retp(C);
endp;

proc PTdownandoutP(S0,K,r,sig,T,nn,vB,times);
   local OV,i,S,P;
   OV=zeros(times,1);
   i=1;
   do while i<=times;
      S=ptsampler(S0,r,sig,T,nn);
      if S>=vB;                  /* vector evaluation */
         OV[i]=K-S[nn+1];  /* effective  */
      else;
         OV[i]=0;          /* value=0    */
      endif;
      i=i+1;
   endo;
   P=exp(-r*T)*meanc(maxc((OV~zeros(times,1))'));
   retp(P);
endp;



/*
**  ptsampler.txt - Stochastic Trinomial Path Sampler.
** (C) Copyright 2005 Yosuke Amijima. All Rights Reserved.
**
**  Purpose:    Gets a stochastic path on trinomial tree.
**
**  Format:     S=ptsampler(S0,r,sig,T,nn)
**
**  Input:      S0       scalar,    initial value
**
**              r        scalar,    risk-free interest rate
**
**              sig      scalar,    volatility
**
**              T        scalar,    maturity
**
**              nn       scalar,    number of time steps
**
**
**  Output:     S        vector,    (nn+1) x 1 of resulting values including S0
**
**  Notice:     This procedure uses 'ptrisampler' inside.
**
*/
proc ptsampler(S0,r,sig,T,nn);
   local delt,pd,pm,pu,p,S,u;
   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=pu|pm|pd;
   S=ptrisampler(p,nn);
   S=S-1;
   u=exp(sig*sqrt(3*delt));
   S=S0*u^cumsumc(S);
   S=S0|S;
   retp(S);
endp;

/*
**  ptrisampler.txt - Stochastic Trinomial 0-1-2 Sampler.
** (C) Copyright 2005 Yosuke Amijima. All Rights Reserved.
**
**  Purpose:    Gets stochastic Trinomial 0-1-2 index numbers in a very easy way.
**
**  Format:     x=ptrisampler(p,nr);
**
**  Input:      p        vector,    3 x 1 of probability vector(p[1] for 2, p[2] for 1, p[3] for 0)
**
**              nr       scalar,    number of rows
**
**
**  Output:     x        vector,    nr x 1 of resulting 0-1 index vector
**
*/
proc ptrisampler(p,nr);
   local x;
   p=cumsumc(p);
   x=rndu(nr,1);
   x=2*(x.<=p[1])+1*(x.>p[1] .and x.<=p[2]);
   retp(x);
endp;