• C++ Programming for Financial Engineering
    Highly recommended by thousands of MFE students. Covers essential C++ topics with applications to financial engineering. Learn more Join!
    Python for Finance with Intro to Data Science
    Gain practical understanding of Python to read, understand, and write professional Python code for your first day on the job. Learn more Join!
    An Intuition-Based Options Primer for FE
    Ideal for entry level positions interviews and graduate studies, specializing in options trading arbitrage and options valuation models. Learn more Join!

Assistance with C++ and Projected SOR Algorithm

Joined
12/10/10
Messages
6
Points
11
Hello,

I am getting incorrect results while trying to price American option using Projected SOR. This is for a class project, and I have been instructed to implement it by transforming Black-Scholes into Heat equation for the algorithm.

For this assignment:
K = 25;
r = 0.02;
volatility = 0.05;
Dividend Rate = 0.03;
T=1;

(Delta)X = 0.05
(Delta)t = 0.00125
so lambda = (delta)t/(delta)x^2 = 0.5

For the algorithm, the range for x is -2.5 <= x <= 2.5 after converting to heat equation.

I am to find prices for values of S=5, 25, and 125.

Variables: jMax = iMax = 101 = size of vectors/arrays since C++ uses zero-based vector indexing.

qDiv = (2 * (r - dividend yield)) / volatility^2
q = (2 * r) / volatility^2

The implementation is per class notes that I can attach if needed. The code segment for the American Call/Implicit FDM scheme is below.

Code:
//Runs projected SOR for American Call with implicit scheme
void ProjSORCallImp(int jMax, int iMax, double lambda, double qDiv, double q)
{

   //Heading for output of values
   cout << endl;
   cout << "Implicit Scheme Calculation for American Call";
   cout << endl;
   cout << "Using Projected SOR Iterative Algorithm." << endl;
   cout << endl;

   int i, j;
   double tau;
   double b[jMax];
   double vPrev[jMax];
   double vNext[jMax];
   double w[jMax];
   double g[jMax][iMax];
   double terminateCondition;
   double compareValue;


   //populate g matrix   
   for (j = 0; j < jMax; j++)
   {
      for (i = 0; i < iMax; i++)
      {
         compareValue=exp(((qDiv+1)/2)*(-2.5+(j*0.05)))-exp(((qDiv-1)/2)*(-2.5+(j*0.05)));
         g[j][i]=exp(((pow(qDiv-1,2)/4)+q)*(i*.00125))*max(compareValue, (double)0);

      }
   }

   //set up iteration vector
   for (j = 0; j < jMax; j++)
   {
      w[j] = g[j][0];
   }
  

   //Begin i loop (time increments) for prototype core algorithm 
   for(i = 0; i <= iMax - 1; i++)
   {

      tau = (double)i * 0.00125;

      //set up b vector
      b[1]=w[1]+lambda*g[0][i+1];
      
      for (j = 2; j <= jMax - 3; j++)
      {
         b[j]=w[j];
      }
      
      b[jMax-2]=w[jMax-2]+lambda*g[jMax-1][i+1];
     
      //set v = max(w, g+1) 
      for(j = 0; j <= jMax - 1; j++)
      {
         vPrev[j] = max(w[j],g[j][i+1]);
      }

      //Populate with initial values
      for (j = 0; i < jMax; i++)
      {
         vNext[i] = 0;
      }

      //do loop while not within desired tolerance
      do
      {

         for (j = 1; j <= jMax - 2; j++)
         {
            compareValue = vPrev[j]+(1.3/(1+(2*lambda)))*(b[j]+lambda*vNext[j-1]-(1+(2*lambda))*vPrev[j]+lambda*vPrev[j+1]);
           
            vNext[j] = max(g[j][i+1], compareValue);
         }

         terminateCondition = vPrev[82]; //for comparison for tolerance

         for (i = 0; i < jMax; i++)
         {
            vPrev[i] = vNext[i]; //sets vPrev vector for next iteration
         }

      } 
      while (sqrt(pow(vNext[82] - terminateCondition, 2)) > 0.00000001);

      //set new iteration vector
      for(j = 0; j < jMax; j++)
      {
         w[j] = vNext[j];
      }
   }

   double Stock[jMax];
   double Value[jMax];
   for(j = 0; j < jMax; j++)
   {
      Stock[j] = 25.0*exp((j*0.05)-2.5);
      Value[j] = 25.0*exp((((qDiv-1)*((j*0.05)-2.5))/-2)-((pow(qDiv-1,2)/4)+q)*.125)*w[j];
   }


   cout << "S = 5: " << Value[18] << endl;

   cout << "S = 25: " << Value[50] << endl;

   cout << "S = 125: " << Value[82] << endl;


}
Results from this are below:

Implicit Scheme Calculation for American Call
Using Projected SOR Iterative Algorithm.

S = 5: 0
S = 25: 1.17982
S = 125: 96.8086

Any help would be appreciated. I am completely stuck on this and I received little help from the instructor when I asked him about it this morning (he didn't even look at my code!) and I only have until Wednesday night to correct my code.
 
Thanks, unfortunately, I cannot get the package to run (attempted FiniteDifference.cpp and received following error code):
Owner@your-09dedafe33 ~
$ g++ -Wall FiniteDifference.cpp -o FiniteDifference
/tmp/ccEgGENy.o:FiniteDifference.cpp:(.text+0x44c): undefined reference to `Full
Mtx::FullMtx(int, int, char)'
/tmp/ccEgGENy.o:FiniteDifference.cpp:(.text$_ZN7FullMtxD1Ev[FullMtx::~FullMtx()]
+0xd): undefined reference to `FullMtx::CoreDest()'
/usr/lib/gcc/i686-pc-cygwin/3.4.4/../../../libcygwin.a(libcmain.o):(.text+0xa9):
undefined reference to `_WinMain@16'
collect2: ld returned 1 exit status

I am running it in CygWin on WinXP.

In any event, I'm looking for where my code is incorrect so I can make changes to it. My professor specifically wants me to use the algorithm provided in class which is utilized in the code, but giving incorrect results.
 
"In any event, I'm looking for where my code is incorrect so I can make changes to it. My professor specifically wants me to use the algorithm provided in class which is utilized in the code, but giving incorrect results. "

You may have coded incorrectly. Normally, you debug to get correct results. It goes with the territory, as will become clear.

But what is more important is to tell how you set up the fd scheme. There are many risks:

1. The transformation to heat equation complicates the new domain.
2. Domain truncation and the imposition of numerical boundary coditions. (two issues here).
3. What time discretisation are you using?

It is possible to see what the algo is? Debugging the algos is less expensive than debugging code.


//
btw point 1 is a trick to get a symmetric matrix and is not needed. Non-symmetric version of PSOR exist (and well as faster ways to price this, but that is off-topic).<!-- / message -->
 
"In any event, I'm looking for where my code is incorrect so I can make changes to it. My professor specifically wants me to use the algorithm provided in class which is utilized in the code, but giving incorrect results. "

You may have coded incorrectly. Normally, you debug to get correct results. It goes with the territory, as will become clear.

But what is more important is to tell how you set up the fd scheme. There are many risks:

1. The transformation to heat equation complicates the new domain.
2. Domain truncation and the imposition of numerical boundary coditions. (two issues here).
3. What time discretisation are you using?

It is possible to see what the algo is? Debugging the algos is less expensive than debugging code.


//
btw point 1 is a trick to get a symmetric matrix and is not needed. Non-symmetric version of PSOR exist (and well as faster ways to price this, but that is off-topic).<!-- / message -->

I've been working for a week trying to debug this thing:/

I've attached a copy of the class notes that includes the algorithm... The initialization matrix g is covered on slide 9 and the algorithm itself is on slides 12-14.

I'm using the implicit scheme for this implementation.(Theta = 1 in the notes)
 

Attachments

I know the feeling; very annoying when it does not work :-)

I had a look through. The equations look correct. Personally I write PSOR in two steps which makes coding easier.

Attention area is boundary condition on sheets 6 and 8 (I cannot comment because I don't know your

1. q_delta = (r-delta) / (0.5 * sig^2)? [ q_delta < 1 ?? > 1?]
2. q = r / (0.5 * sig^2)?

Are you modelling a put? What value do you get? What value should you get?

What are the numeric values you use for r, sig, T etc. Do the values in S space and x space both look off?
 
For this assignment:
K = 25;
r = 0.02;
volatility (sigma) = 0.05;
Dividend Rate (delta) = 0.03;
T=1;

()X = 0.05
()t = 0.00125
so lambda = ()t/()x^2 = 0.5

For the algorithm, the range for x is -2.5 <= x <= 2.5 after converting to heat equation.

I am to find prices for values of S=5, 25, and 125.

Variables: jMax = iMax = 101 = size of vectors/arrays since C++ uses zero-based vector indexing.

1. q_delta = (2 * (r-delta)) / (sig^2)? [ q_delta < 1 ?? > 1?]
2. q =(2 * r) / (sig^2)?

I've posted the call, but I'm getting similar (incorrect) results for the put:
Implicit Scheme Calculation for American Put
Using Projected SOR Iterative Algorithm.

S = 5: 18.8648
S = 25: 0.342394
S = 125: 0

I don't have target results, only that they have to be equal or greater than the corresponding European values:

Implicit Scheme Calculation for European Call
Using SOR Iterative Algorithm.

S = 5: 0.00109004
S = 25: 4.68111
S = 125: 80.1774

Implicit Scheme Calculation for European Put
Using SOR Iterative Algorithm.

S = 5: 14.7837
S = 25: 4.9249
S = 125: 0.00663343
 
Please try a put

S = K = 100
T = 0.25
r = 0.1
d == q = 0
sig = 0.2

I have a P for this.

What is the answer you get?
 
Please try a put

S = K = 100
T = 0.25
r = 0.1
d == q = 0
sig = 0.2

I have a P for this.

What is the answer you get?

When I change the time parameter, the thing gets caught in an endless loop for the SOR do-while loop.

EDIT: Never mind, that was for a call. Though it is troubling the SOR loop never converges for the call only by changing the time parameter.

For put:
Implicit Scheme Calculation for American Put
Using Projected SOR Iterative Algorithm.

S = 100: 1.51764 by forcing q = 0.
S = 100: 0.0497831 if q = 5 by equation q =(2 * r) / (sig^2)
 
Okay, I'm an idiot... The problem was in this loop...

for (j = 0; i < jMax; i++)
{
vNext = 0;
}

So this loop would increment the i (time step) value during the first run, so I was essentially pricing it at the starting time. Once I changed all those i's to j's it gave me much more reasonable results... Not sure if they're correct, but at least they make sense when compared to my European option price results...

Thanks for all the help, and hopefully I won't have any more problems with the program!
 
When I change the time parameter, the thing gets caught in an endless loop for the SOR do-while loop.

EDIT: Never mind, that was for a call. Though it is troubling the SOR loop never converges for the call only by changing the time parameter.

For put:
Implicit Scheme Calculation for American Put
Using Projected SOR Iterative Algorithm.

S = 100: 1.51764 by forcing q = 0.


P = 3.0732

I see you also discuss LU and Brennan Schwartz. How does it compare to PSOR speed-wise.

BTW I have developed the ADE method for European and American options which is 40% faster than the most optimised LU. It is very easy to program and to parallelise. It could be a nice part of a comparison project.


http://papers.ssrn.com/sol3/papers.cfm?abstract_id=1552926


See other approaches here

http://www.wilmott.com/messageview.cfm?catid=34&threadid=60415&FTVAR_MSGDBTABLE=
 
Okay, I'm an idiot... The problem was in this loop...

for (j = 0; i < jMax; i++)
{
vNext = 0;
}

So this loop would increment the i (time step) value during the first run, so I was essentially pricing it at the starting time. Once I changed all those i's to j's it gave me much more reasonable results... Not sure if they're correct, but at least they make sense when compared to my European option price results...

Thanks for all the help, and hopefully I won't have any more problems with the program!



I alway use 'i' for the space direction and 'n' for time.

Most of the literature and books use 'i' and 'j' which are too close together for comfort.

I once was testing a FDM. I took me 2 days and nights to realise that I programmed the unary minus operator wrong:

v2 = -v1;

!
 
"For the algorithm, the range for x is -2.5 <= x <= 2.5 after converting to heat equation."


What is the rationale/motivation for this specific truncation?
 
Back
Top