Here we aim to solve the classic American put option problem using the PSOR method. We must solve the problem ∂P∂t+12σ2S2∂2P∂S2+rS∂P∂S−rP=0 with P(S,T)=max P(S,t) \geq X-S, and S_0 = 1.973 \text{ , } X=2 \text{ , } T=1 \text{ , } r=0.05 \text{ , } \sigma=0.4 ~. The method of PSOR is almost identical to the SOR method, but there are two important differences.
This means that we update the algorithm where v_j^q is the qth guess at the solution, and \epsilon_j is our current measure of the residual error, to
y = \left(d_j - a_jv^{q+1}_{j-1} - c_jv^q_{j+1}\right)/b_j
v^{q+1/2}_j=v^q_j+\omega(y-v^q_j) ,
v^{q+1}_j=\max\left(v^{q+1/2}_j, X - S_j\right) .
\epsilon_j = \left( v^{q+1}_j - v^q_j \right)^2
Then the exit condition on the algorithm becomes stop if
\sum_j \epsilon_j < \text{tol}^2
which should guarantee our solution is accurate to the level tol. The final solution is also on the web click here.
First include standard libraries:
#include <iostream>
#include <vector>
#include <chrono>
#include <string>
#include <fstream>
#include <cmath>
#include <algorithm>
using namespace std;
Next we have packaged the American put option solver up into a function. This code has been adapted from the European put option solver with SOR, so the inputs/outputs are the same. There is full flexibility on the choice of numerical parameters.
/**
@brief Return the value of an American Put option using the Crank-Nicolson method with PSOR solver
ON INPUT:
@param S0 -- initial stock price
@param X -- exercise (strike) price
@param T -- Time to expiry (years)
@param r -- interest rate (per annum)
@param sigma -- volatility (per annum^12)
@param iMax -- number of time steps
@param jMax -- number of space steps
@param SMax -- Maximum value of S
@param omega -- is the relaxation parameter
@param tol -- is the tolerance level
@param iterMax -- is maximum iterations
ON OUTPUT:
@return the value of a American put option at \f$ S=S0 \f$, \f$ t=0 \f$
**/
double americanPut_CN(double S0,double X,double T,double r,double sigma,int iMax,int jMax,double SMax,double omega,double tol,int iterMax)
{
double dS=SMax/jMax;
double dt=T/iMax;
// create storage for the stock price and option price (old and new)
std::vector<double> S(jMax+1),vOld(jMax+1),vNew(jMax+1);
// setup and initialise the stock price
for(int j=0;j<=jMax;j++)
{
S[j] = j*dS;
}
// reset initial condition
for(int j=0;j<=jMax;j++)
{
vOld[j] = std::max(X-S[j],0.);
vNew[j] = std::max(X-S[j],0.);
}
// run through timesteps
for(int i=iMax-1;i>=0;i--)
{
// declare vectors for matrix equations
std::vector<double> a(jMax+1),b(jMax+1),c(jMax+1),d(jMax+1);
// set up matrix equations a[j]=
a[0] = 0.;b[0] = 1.;c[0] = 0.;
d[0] = X*exp(-r*(iMax-i)*dt);
for(int j=1;j<jMax;j++)
{
a[j] = 0.25*(sigma*sigma*j*j-r*j);
b[j] =-0.5*sigma*sigma*j*j - 0.5*r - 1./dt;
c[j] = 0.25*(sigma*sigma*j*j+r*j);
d[j] =-a[j]*vOld[j-1] - (b[j]+2./dt)*vOld[j] - c[j]*vOld[j+1];
}
a[jMax] = 0.;b[jMax] = 1.;c[jMax] = 0.;
d[jMax] = 0.;
// solve with SOR method
int sor;
for(sor=0;sor<iterMax;sor++)
{
double error=0.;
// implement sor in here
{
double y = (d[0] - c[0]*vNew[1])/b[0];
y = vNew[0] + omega*(y-vNew[0]);
y = std::max(y,X-S[0]);
error+=(y-vNew[0])*(y-vNew[0]);
vNew[0] = y;
}
for(int j=1;j<jMax;j++)
{
double y = (d[j] - a[j]*vNew[j-1] - c[j]*vNew[j+1])/b[j];
y = vNew[j] + omega*(y-vNew[j]);
y = std::max(y,X-S[j]);
error+=(y-vNew[j])*(y-vNew[j]);
vNew[j] = y;
}
{
double y = (d[jMax] - a[jMax]*vNew[jMax-1])/b[jMax];
y = vNew[jMax] + omega*(y-vNew[jMax]);
y = std::max(y,X-S[jMax]);
error+=(y-vNew[jMax])*(y-vNew[jMax]);
vNew[jMax] = y;
}
// make an exit condition when solution found
if(error<tol*tol)
break;
}
if(sor>=iterMax)
{
std::cout << " Error NOT converging within required iterations\n";
std::cout.flush();
throw;
}
// set old=new
vOld=vNew;
}
int jStar=S0/dS;
double sum=0.;
sum+=(S0 - S[jStar])/dS * vNew[jStar+1];
sum+=(S[jStar+1] - S0)/dS * vNew[jStar];
return sum;
}
Now we declare our Black-Scholes parameters
// declare and initialise Black Scholes parameters
double S0,X,T,r,sigma;
// declare and initialise grid paramaters
int iMax,jMax;
// declare and initialise local variables (ds,dt)
double SMax;
and initialise them
// initialise Black Scholes parameters
S0=1.973;X=2.;T=1.;r=0.05;sigma=0.4;
// initialise grid paramaters
iMax=40;jMax=40;SMax=2*X;
before checking that our algorithm works and generates results.
cout << americanPut_CN(S0,X,T,r,sigma,iMax,jMax,SMax,1,1.e-6,1000);
Next we can look at getting some results for different values of n when iMax=n and jMax=n.
cout << " n | V_n \n";
for(int k=1;k<=4;k++)
{
int n=10*pow(2,k);
iMax = n;
jMax = n;
double value = americanPut_CN(S0,X,T,r,sigma,iMax,jMax,SMax,1.2,1.e-6,1000);
cout << n << " | " << value << "\n";
}
Obviously we do not have an analytical value for this option but there are many places that you can find an estimate for the value. Even on my website I have an online American put pricer written in javascript that estimates the value as V_\text{exact}(S=1.973,t=0) = 0.284193 so we can see that these estimates are not very far from the true value. If we are only interested in 4s.f. then just n=160 provides a good estimate. If we want something more accurate we will need to be more careful, including the interpolation and tolerance as before.
So for the next example (with extrapolation) we use S_0=X, so that interpolation errors are removed and also increase S_\text{max} and the tolerance. The new exact value we are aiming for is V_\text{exact}(S=X=2,t=0) = 0.273352 .
double valueOld=0;
S0 = X;
cout << " n | V_n | diff | V_extrap \n";
for(int k=1;k<=5;k++)
{
int n=10*pow(2,k);
iMax = n;
jMax = n;
double value = americanPut_CN(S0,X,T,r,sigma,iMax,jMax,5*X,1.2,1.e-8,1000);
double value_extrap = (4.*value - valueOld)/3.;
cout << n << " | " << value << " | " ;
cout << value - valueOld << " | ";
cout << value_extrap << "\n";
valueOld = value;
}
Again we are able to get quite smooth convergence although not quite as good as previously. This time the ratio of differences between the results is not perfect so the extrapolation doesn't perform quite as well. We are still able to improve our results though and the combination of n=160 and n=320 is able to give 5 almost 6s.f. of accuracy, whereas the standard method with n=320 gives only 3 or 4s.f. of accuracy.