/* A Deterministic Path on Trinomial Tree */
new; cls;
S0=50;
r=0.10;
sig=0.40;
T=5/12;
nn=5;  /* Not so large in GAUSS Light. It could explode. */
i=floor(3^nn*rndu(1,1));
S=tsampler(S0,sig,T,nn,i);
print S;
print/ld "p=" pt(S0,r,sig,T,nn,i);

library pgraph;
graphset;
title("A Path on Trinomial Tree");
ylabel("S");
xlabel("t");
xy(seqa(0,T/nn,nn+1),S);


/*
**  tsampler.txt - Deterministic Trinomial Path Sampler.
** (C) Copyright 2005 Yosuke Amijima. All Rights Reserved.
**
**  Purpose:    Gets a deterministic path on trinomial tree for given i.
**
**  Format:     S=tsampler(S0,sig,T,nn,i)
**
**  Input:      S0       scalar,    initial value
**
**              sig      scalar,    volatility
**
**              T        scalar,    maturity
**
**              nn       scalar,    number of time steps
**
**              i        scalar,    deterministic order (0,1,2,...,3^nn-1)
**
**
**  Output:     S        vector,    (nn+1) x 1 of resulting values including S0
**
**  Notice:     This procedure uses 'trisampler' inside.
**
**              Parameter nn should not be large because a path could blow up if
**              you use Light version of GAUSS or memory is limited. It depends.
**              This comes from matrix transformation in procedure trisampler.
*/
proc tsampler(S0,sig,T,nn,i);
   local S,delt,u;
   S=trisampler(i,nn);
   S=S-1*(S.==0)-1*(S.==1)-1*(S.==2);
   delt=T/nn;
   u=exp(sig*sqrt(3*delt));
   S=S0*u^cumsumc(S);
   S=S0|S;
   retp(S);
endp;

/*
**  tsampler.txt - Probability of the Deterministic Trinomial Path.
** (C) Copyright 2005 Yosuke Amijima. All Rights Reserved.
**
**  Purpose:    Gets the probability of the path.
**
**  Format:     p=pt(S0,r,sig,T,nn,i);
**
**  Input:      S0       scalar,    initial value
**
**              r        scalar,    risk-free interest rate
**
**              sig      scalar,    volatility
**
**              T        scalar,    maturity
**
**              nn       scalar,    number of time steps
**
**              i        scalar,    deterministic order (0,1,2,...,2^nn-1)
**
**
**  Output:     p        scalar,    probability of the path
**
**  Notice:     This procedure uses 'trisampler' inside.
**
*/
proc pt(S0,r,sig,T,nn,i);
   local delt,pd,pm,pu,p,S;
   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;
   S=trisampler(i,nn);
   p=pd*(S.==0)+pm*(S.==1)+pu*(S.==2);
   p=prodc(p);
   retp(p);
endp;

/*
**  trisampler.txt - Deterministic Trinomial 0-1-2 Sampler.
** (C) Copyright 2005 Yosuke Amijima. All Rights Reserved.
**
**  Purpose:    Gets deterministic trinomial 0-1-2 index numbers for given i in a very easy way.
**
**  Format:     x=trisampler(i,nr);
**
**  Input:      i        scalar,    deterministic order (0,1,2,...,3^nr-1)
**
**              nr       scalar,    number of rows
**
**
**  Output:     x        vector,    nr x 1 of resulting 0-1-2 index vector
**
*/
proc trisampler(i,nr);
   local x,j;
   if i<0 or i>(3^nr-1) or i-floor(i)/=0;
      errorlog "ERROR: Parameter i must be an integer(0,1,2,...,3^nn-1).";
      retp(".");
   endif;
   x=zeros(nr,1);
   j=nr;
   do until j==1;
      x[j]=floor(i/(3^(j-1)));
      i=i-(x[j]*(3^(j-1)));
      j=j-1;
   endo;
   x[1]=i;
   retp(rev(x));
endp;