new; cls;
/* Stochastic Binomial Asian Call & Put Options  */
S0=50;
K=50;
r=0.10;
sig=0.40;
T=5/12;
nn=500;
times=5000;

print/lz "     nn=" nn;
print/lz "  times=" times;
print "C=" PBasianC(S0,K,r,sig,T,nn,times);
print "P=" PBasianP(S0,K,r,sig,T,nn,times);
print;
print "Another version:";
print "C=" PBasianC2(S0,r,sig,T,nn,times);
print "P=" PBasianP2(S0,r,sig,T,nn,times);
print;
print "Geometric Mean Asian:";
print "C=" PBgasianC(S0,K,r,sig,T,nn,times);
print "P=" PBgasianP(S0,K,r,sig,T,nn,times);
print;
print "Another version of Geometric Mean:";
print "C=" PBgasianC2(S0,r,sig,T,nn,times);
print "P=" PBgasianP2(S0,r,sig,T,nn,times);



/*
**  pbasian.txt - Stochastic Binomial Asian call & Put Options.
** (C) Copyright 2005 Yosuke Amijima. All Rights Reserved.
**
**  Purpose:    Calculates Asian option prices by stochastic binomial simulation.
**
**  Format:     C=PBasianC(S0,K,r,sig,T,nn,times);
**              P=PBasianP(S0,K,r,sig,T,nn,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
**
**              times    scalar,  number of simulations
**
**
**  Output:     C        scalar,  call option price
**              P        scalar,  put  option price
**
*/

proc PBasianC(S0,K,r,sig,T,nn,times);
   local AV,i,S,C;
   AV=zeros(times,1);
   i=1;
   do while i<=times;
      S=pbsampler(S0,r,sig,T,nn);
      AV[i]=meanc(S[2:nn+1]);
      i=i+1;
   endo;
   C=exp(-r*T)*meanc(maxc(((AV-K)~zeros(times,1))'));
   retp(C);
endp;

proc PBasianP(S0,K,r,sig,T,nn,times);
   local AV,i,S,P;
   AV=zeros(times,1);
   i=1;
   do while i<=times;
      S=pbsampler(S0,r,sig,T,nn);
      AV[i]=meanc(S[2:nn+1]);
      i=i+1;
   endo;
   P=exp(-r*T)*meanc(maxc(((K-AV)~zeros(times,1))'));
   retp(P);
endp;

/*
**  pbasian.txt - Another version of Stochastic Binomial Asian Call & Put Options.
** (C) Copyright 2005 Yosuke Amijima. All Rights Reserved.
**
**  Purpose:    Calculates Asian option prices by stochastic binomial simulation.
**
**  Format:     C=PBasianC2(S0,r,sig,T,nn,times);
**              P=PBasianP2(S0,r,sig,T,nn,times);
**
**  Input:      S0       scalar,  initial value
**             
**              r        scalar,  risk-free interest rate
**
**              sig      scalar,  volatility
**
**              T        scalar,  maturity
**
**              nn       scalar,  number of time steps
**
**              times    scalar,  number of simulations
**
**
**  Output:     C        scalar,  call option price
**              P        scalar,  put  option price
**
*/

proc PBasianC2(S0,r,sig,T,nn,times);
   local OV,i,S,C;
   OV=zeros(times,1);
   i=1;
   do while i<=times;
      S=pbsampler(S0,r,sig,T,nn);
      OV[i]=S[nn+1]-meanc(S[2:nn+1]);
      i=i+1;
   endo;
   C=exp(-r*T)*meanc(maxc((OV~zeros(times,1))'));
   retp(C);
endp;

proc PBasianP2(S0,r,sig,T,nn,times);
   local OV,i,S,P;
   OV=zeros(times,1);
   i=1;
   do while i<=times;
      S=pbsampler(S0,r,sig,T,nn);
      OV[i]=meanc(S[2:nn+1])-S[nn+1];
      i=i+1;
   endo;
   P=exp(-r*T)*meanc(maxc((OV~zeros(times,1))'));
   retp(P);
endp;

/*
**  pbasian.txt - Stochastic Binomial Geometric Asian call & Put Options.
** (C) Copyright 2005 Yosuke Amijima. All Rights Reserved.
**
**  Purpose:    Calculates Geometric Asian option prices by stochastic binomial simulation.
**
**  Format:     C=PBgasianC(S0,K,r,sig,T,nn,times);
**              P=PBgasianP(S0,K,r,sig,T,nn,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
**
**              times    scalar,  number of simulations
**
**
**  Output:     C        scalar,  call option price
**              P        scalar,  put  option price
**
*/

proc PBgasianC(S0,K,r,sig,T,nn,times);
   local AV,i,S,C;
   AV=zeros(times,1);
   i=1;
   do while i<=times;
      S=pbsampler(S0,r,sig,T,nn);
      AV[i]=geomean(S[2:nn+1]);
      i=i+1;
   endo;
   C=exp(-r*T)*meanc(maxc(((AV-K)~zeros(times,1))'));   /* We use arithmetic one here. */
   retp(C);
endp;

proc PBgasianP(S0,K,r,sig,T,nn,times);
   local AV,i,S,P;
   AV=zeros(times,1);
   i=1;
   do while i<=times;
      S=pbsampler(S0,r,sig,T,nn);
      AV[i]=geomean(S[2:nn+1]);
      i=i+1;
   endo;
   P=exp(-r*T)*meanc(maxc(((K-AV)~zeros(times,1))'));   /* we use arithmetic one here. */
   retp(P);
endp;

/*
**  pbasian.txt - Another version of Stochastic Binomial Geometric Asian call & Put Options.
** (C) Copyright 2005 Yosuke Amijima. All Rights Reserved.
**
**  Purpose:    Calculates Geometric Asian option prices by stochastic binomial simulation.
**
**  Format:     C=PBgasianC2(S0,r,sig,T,nn,times);
**              P=PBgasianP2(S0,r,sig,T,nn,times);
**
**  Input:      S0       scalar,  initial value
**             
**              r        scalar,  risk-free interest rate
**
**              sig      scalar,  volatility
**
**              T        scalar,  maturity
**
**              nn       scalar,  number of time steps
**
**              times    scalar,  number of simulations
**
**
**  Output:     C        scalar,  call option price
**              P        scalar,  put  option price
**
*/

proc PBgasianC2(S0,r,sig,T,nn,times);
   local OV,i,S,C;
   OV=zeros(times,1);
   i=1;
   do while i<=times;
      S=pbsampler(S0,r,sig,T,nn);
      OV[i]=S[nn+1]-geomean(S[2:nn+1]);
      i=i+1;
   endo;
   C=exp(-r*T)*meanc(maxc((OV~zeros(times,1))'));
   retp(C);
endp;

proc PBgasianP2(S0,r,sig,T,nn,times);
   local OV,i,S,P;
   OV=zeros(times,1);
   i=1;
   do while i<=times;
      S=pbsampler(S0,r,sig,T,nn);
      OV[i]=geomean(S[2:nn+1])-S[nn+1];
      i=i+1;
   endo;
   P=exp(-r*T)*meanc(maxc((OV~zeros(times,1))'));
   retp(P);
endp;

proc geomean(x);
   local n;
   n=rows(x);
   retp(exp(1/n*sumc(ln(x))));
endp;


/*
**  pbsampler.txt - Stochastic Binomial Path Sampler.
** (C) Copyright 2005 Yosuke Amijima. All Rights Reserved.
**
**  Purpose:    Gets a stochastic path on binomial tree.
**
**  Format:     S=pbsampler(S0,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 'pbisampler' inside.
**
*/
proc pbsampler(S0,r,sig,T,nn);
   local delt,a,u,d,p,S;
   delt=T/nn;
   a=exp(r*delt);
   u=exp(sig*sqrt(delt));
   d=1/u;
   p=(a-d)/(u-d);
   S=pbisampler(p,nn);
   S=S-1*(S.==0);
   S=S0*u^cumsumc(S);
   S=S0|S;
   retp(S);
endp;

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