/* A Deterministic Path on Binomial Tree   */
/* 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.       */
new; cls;
S0=50;
r=0.10;
sig=0.40;
T=5/12;
x={1,1,0,1,0,1,0,1,0,0,1,1,1,0,0,1,1,0,1,0,1,0,1,0,0,1,1,1,0,0,1,1,0,
1,0,1,0,1,0,0,1,0,1,0,0,1,1,0,0,0,1,0,1,0,0,1,1,1,0,0,1};   /* 0 or 1 you could change */
S=bsampler2(S0,sig,T,x);


print S;
print "p=" pb2(S0,r,sig,T,x);

library pgraph;
graphset;
title("A Path on Binomial Tree");
ylabel("S");
xlabel("t");
xy(seqa(0,T/(rows(S)-1),rows(S)),S);



/*
**  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;