Longstaff Schwartz method

Question

I try to implemente the LSM method with this algorithm but my price is always too low. By example for an American put option with the following parameters:

S0 = 36, Strike = 40, rate = 6%, T = 1 year, discrete path = 50, volatility = 20%

I got 4 dollars, but the Longstaff and Schwartz article lists 4.7 dollars. With a volatility of 40%, the error is bigger at 5 dollars for me vs. 7.3 dollars for L&S. But with my tree pricer I have the same result as the L&S article.

Could you help me to find the error please?

void LeastSquaresMC::calcLeastSquaresMC()
{

     mu_ = (rate_ - vol_*vol_*0.5)*dt; // drift
     voldt = vol_*sqrt(dt); // diffusion
          for (i = 0; i < M_; i++)
          {
               Paths(i,0) = 36;

               for (j = 1; j < b; j++)
                    {
                         // generate deviate
                         deviate = G();
                         Paths(i,j) =  Paths(i,j-1)*exp(mu_+voldt*deviate);                         
                    }
          }
     // initialize cash flow matrix by zero
     for (i = 0; i < z; i++)
          {
               for (j = 0; j < b; j++)
                    {
                         CashFlow(i,j,0);
                    }
          }

     for (i = 0; i < z; i++)
          {
               for (j = 0; j < b; j++)
                    {
                         Exercise(i,j) = MAX(strike_-Paths(i,j),0);
                    }
          }
     // compute cash flows at maturity
     for (i = 0; i < z; i++)
          {
               CashFlow(i,b-1,(Exercise(i,b-1)));

          }
     //cout <<CashFlow << endl;
     // recursion
     computeLSM(b-1, Paths, CashFlow, Exercise);

}

double LeastSquaresMC::computeLSM(int time, Matrix& Paths, Matrix& CashFlow, Matrix& Exercise)
{

     double disc = exp(-rate_*dt);     // discount factor
     vector<double> Y;               // vector of payoffs (dependent variables)
     vector<double> B;               // vector of regression coefficients
     vector<double> C;               // continuation
     vector<int> num;
     vector<double> stock;          
     vector<int>::iterator i = num.begin();


     /*long z = M_*2;*/

     for (j = 0; j < z; j++)
          {
               if(Exercise(j,time-1)>0)
                    {

                                   Y.push_back(MAX(CashFlow(j,time),0)*disc);
                                   num.push_back(j);
                                   stock.push_back(Paths(j,time-1));                    
                    }
          }

     if (time > 1)
          {
               if(num.empty()==false)
                    {
                         int size_l = Y.size();
                         Matrix X(size_l,3);    // 1 X X^2 (columns)

                         for (j = 0; j < size_l; j++)
                              {
                                   X(j,0,1);
                                   X(j,1,stock[j]);
                                   X(j,2,stock[j]*stock[j]);
                              }
                         B = ((X.transpose()*X).Inverse())*(X.transpose()*Y);
                         C = X*B;
                         j=0;
                         for(i = num.begin() ; i != num.end(); ++i)
                              {
                              if (Exercise(*i,time-1)>C[j])

                                   {

                                        CashFlow(*i,time-1,Exercise(*i,time-1));
                                        for (l = time; l < b; l++)
                                             {
                                                  CashFlow(*i,l,0);
                                             }
                              j++;
                              }
                         computeLSM(time-1, Paths, CashFlow, Exercise);
                    }
               else
                    {
                         computeLSM(time-1, Paths, CashFlow, Exercise);
                    }
               }
          else
               {
                    return computeValue(CashFlow);
               }

     return 0.0;     
}

double LeastSquaresMC::computeValue (Matrix& CashFlow)
{


     double discValue = 0.0; // discounted value
          for (i = 0; i < z; i++)
               {
                    for (j = 1; j < b; j++)
                         {
                              if (CashFlow(i, j) > 0)
                              {
                                   discValue = discValue + CashFlow(i, j)*exp(-0.06*j);
                              }
                         }
               }
          cout <<"prix:"<<discValue/z << endl;
 return discValue/z;
}

Answer 1

As noted by others, the code is very hard to read. What I spotted: is the discounting done right? I see you discount the continuation value only to calculate Y, but does the discounting enter the recursion?

(I have an implementation of the LS in Java here: http://www.finmath.net/java )

Answer 2

I think that you are simply discounting the cash flows incorrectly (j is an index):

Nearly at the end of your listing, instead of writing

discValue = discValue + CashFlow(i, j)*exp(-0.06*j);

you should write

discValue = discValue + CashFlow(i, j)*exp(-0.06*j*dt);

Answer 3

That bug sounds familiar to me, from when I implemented this myself in python. I can't figure out what your code is doing (it's too wordy, with too much whitespace, too much fragmentation into functions, has strange indents, etc.) But my guess is that either you're not actually using the regression results as quant_dev mentioned, or you're accidentally truncating a float to an int somewhere -- my recollection is that one of those two things is what bit me; using the same test data, I also got 4.0 the first time I ran it.