new; cls;
let m[8,3]=
1.09  1.08  1.34
1.16  1.26  1.54
1.22  1.07  1.03
0.93  0.97  0.92
1.11  1.56  1.52
0.76  0.77  0.90
0.92  0.84  1.01
0.88  1.22  1.34
;
K=1.10;
r=0.06;
delt=1;


print "P=" PameLSM(m,K,r,delt);
print "C=" CameLSM(m,K,r,delt);



/*
**  amelsm.txt - American Call & Put Options by LSM.
** (C) Copyright 2004 Yosuke Amijima. All Rights Reserved.
**
**  Reference:  Longstaff,F.A. and E.S.Schwartz 2001,
**              "Valuing American Options by Simulation: A Simple Least-Squares Approach,"
**              The Review of Financial Studies, Vol. 14, No.1, pp.113-147
**
**
**  Purpose:    Calculates american option prices by LSM approach of 
**              Longstaff & Schwartz(2001) from path data.
**
**  Format:     C=CameLSM(m,K,r,delt);
**              P=PameLSM(m,K,r,delt);
**
**  Input:      m        matrix,  row(# of paths) x col(# of period -1) of path data
**                                cutting the initial values at t=0.
**
**              K        scalar,  strike price
**             
**              r        scalar,  risk-free interest rate
**
**              delt     scalar,  period interval
**
**
**  Output:     C        scalar,  call option price
**              P        scalar,  put  option price
**
**  Notice:     This procedure requires lots of memory to run. If you use a large number
**              of paths or data periods, you need normal GAUSS + large memory setting.
**
*/

proc PameLSM(m,K,r,delt);
   local v,p,index,i,j,ii,X,Y,one,XX,YY,X1,b,v1,stoprule,v0,invXX;
   v=maxc(((K-m[.,cols(m)])~zeros(rows(m),1))'); 
   p=zeros(rows(m),cols(m));
   p[.,cols(m)]=v;
   index=zeros(rows(m),cols(m));
   index[.,cols(m)]=(v.>0);
   j=cols(m);
   do while j>=2;
      Y=miss(zeros(rows(m),1),0);
      X=miss(zeros(rows(m),1),0);
      one=miss(zeros(rows(m),1),0);
      i=1;
      do while i<=rows(m);
         if K-m[i,j-1]>0;
            Y[i]=v[i]*exp(-r*delt); 
            X[i]=m[i,j-1];
            one[i]=1;
         endif;
         i=i+1;
      endo;
      XX=delif(X,X[.,1].==miss(0,0));
      YY=delif(Y,Y[.,1].==miss(0,0));
      one=delif(one,one[.,1].==miss(0,0));
      if one/=miss(0,0);   
         X1=one~XX~XX^2;
         trap 1;
         invXX=inv(X1'X1);
         trap 0;
         if scalerr(invXX);
            v=maxc(((K-m[.,j-1])~zeros(rows(m),1))');
         else;
            b=inv(X1'X1)*X1'YY;
            v1=X1*b; 
            v=maxc(((K-m[.,j-1])~zeros(rows(m),1))');
            i=1;
            ii=1;
            do while i<=rows(m);
               if K-m[i,j-1]>0;
                  if v[i]<v1[ii];
                     v[i]=0;
                  endif;
                  ii=ii+1;
               endif;
               i=i+1;
            endo;
         endif;
         index[.,j-1]=(v.>0);
      else;
         index[.,j-1]=zeros(rows(m),1);
         v=maxc(((K-m[.,j-1])~zeros(rows(m),1))');
      endif;
      p[.,j-1]=v;
      j=j-1;
   endo;
   stoprule=zeros(rows(m),cols(m));
   i=1;
   do while i<=rows(m);
      j=1;
      do while j<=cols(m);
         if index[i,j]==1;
            stoprule[i,j]=1;
            break;
         endif;
         j=j+1;
      endo;
      i=i+1;
   endo;
   v0=sumc(sumc(stoprule.*p.*exp(-r*delt*seqa(1,1,cols(m))')))/rows(m);
   retp(v0);
endp;

proc CameLSM(m,K,r,delt);
   local v,p,index,i,j,ii,X,Y,one,XX,YY,X1,b,v1,stoprule,v0,invXX;
   v=maxc(((m[.,cols(m)]-K)~zeros(rows(m),1))'); 
   p=zeros(rows(m),cols(m));
   p[.,cols(m)]=v;
   index=zeros(rows(m),cols(m));
   index[.,cols(m)]=(v.>0);
   j=cols(m);
   do while j>=2;
      Y=miss(zeros(rows(m),1),0);
      X=miss(zeros(rows(m),1),0);
      one=miss(zeros(rows(m),1),0);
      i=1;
      do while i<=rows(m);
         if m[i,j-1]-K>0;
            Y[i]=v[i]*exp(-r*delt); 
            X[i]=m[i,j-1];
            one[i]=1;
         endif;
         i=i+1;
      endo;
      XX=delif(X,X[.,1].==miss(0,0));
      YY=delif(Y,Y[.,1].==miss(0,0));
      one=delif(one,one[.,1].==miss(0,0));
      if one/=miss(0,0);   
         X1=one~XX~XX^2;
         trap 1;
         invXX=inv(X1'X1);
         trap 0;
         if scalerr(invXX);
            v=maxc(((m[.,j-1]-K)~zeros(rows(m),1))');
         else;
            b=inv(X1'X1)*X1'YY;
            v1=X1*b; 
            v=maxc(((m[.,j-1]-K)~zeros(rows(m),1))');
            i=1;
            ii=1;
            do while i<=rows(m);
               if m[i,j-1]-K>0;
                  if v[i]<v1[ii];
                     v[i]=0;
                  endif;
                  ii=ii+1;
               endif;
               i=i+1;
            endo;
         endif;
         index[.,j-1]=(v.>0);
      else;
         index[.,j-1]=zeros(rows(m),1);
         v=maxc(((m[.,j-1]-K)~zeros(rows(m),1))');
      endif;
      p[.,j-1]=v;
      j=j-1;
   endo;
   stoprule=zeros(rows(m),cols(m));
   i=1;
   do while i<=rows(m);
      j=1;
      do while j<=cols(m);
         if index[i,j]==1;
            stoprule[i,j]=1;
            break;
         endif;
         j=j+1;
      endo;
      i=i+1;
   endo;
   v0=sumc(sumc(stoprule.*p.*exp(-r*delt*seqa(1,1,cols(m))')))/rows(m);
   retp(v0);
endp;