The example question for these solutions can be found on my website (click here).

#### 6.1 Explicit Finite Difference For Option Pricing

In this example we are going to price a European call option with explicit finite difference.

##### Initial Setup

For the explicit method we shall need:

- All parameters for the option, such as X and S
_{0}etc. - The number of divisions in stock, jMax, and divisions in time iMax
- The size of the divisions ΔS = S
_{max}∕jMax and Δt = T∕iMax - A vector to store the stock prices, and two to store the option values at the current and previous time level.

It will make things easier if you just declare these at the top of the program before you start doing anything with them.

The program structure will look something like this:

Use the code below to as a template to get started.

#include<iostream>#include<fstream>#include<cmath>#include<vector>usingnamespacestd;/* Template code for the Explicit Finite Difference*/intmain() {// declare and initialise Black Scholes parameters// declare and initialise grid paramaters// declare and initialise local variables (ds,dt)// create storage for the stock price and option price (old and new)// setup and initialise the stock price// setup and initialise the final conditions on the option price// loop through time levels, setting the option price at each grid point, and also on the boundaries// output the estimated option price}

Inside the template the comments list the rough layout of your program. First declare the
parameters for the problem, local grid variables and the vectors required. Set the values of the
parameters to the following: σ = 0.4, r = 0.05, X = 2, T = 1, and iMax = jMax = 4. Next we
need to note that in option pricing problems we must for S over the semi inifinite domain,
therefore numerically we need to choose an appropriate S_{max}. In general, you should set this
according to mathematical reasoning if you want your code to be generic and work in all
cases. Often we just choose S_{max} = K × X where K is some constant (say 5) but it
should really be linked to the probability distribution of the underlying asset. For
simplicity here we choose S_{max} = 2X so that small grids still work, for any decent
level of accuracy S_{max} will need to set much higher. Next assuming we have set the
correct level of storage to your vector, calculate the values of ΔS and Δt for your grid
and using them to assign initial values to the stock price vector, and the option value
vectors.

It is a good idea at this stage to output the stock values and option values to check that the results are as you would expect them to be. The lines added to our template might look something like this

// declare and initialise Black Scholes parametersdouble S0=100.,X=100.,T=1.,r=0.06,sigma=0.2;// declare and initialise grid paramatersint iMax=4,jMax=4;// declare and initialise local variables (ds,dt)double S_max=2*X; double dS=S_max/jMax; double dt=T/iMax;// create storage for the stock price and option price (old and new)vector<double>S(jMax+1),vOld(jMax+1),vNew(jMax+1);// setup and initialise the stock pricefor(int j=0;j<=jMax;j++) { S[j] = j*dS; }// setup and initialise the final conditions on the option pricefor(int j=0;j<=jMax;j++) { vOld[j] =max(S[j]-X,0.); vNew[j] =max(S[j]-X,0.); cout << iMax << " " << j << " " << S[j] << " " << vNew[j] << " " << vOld[j] << endl; }/* OUTPUT4 0 0 0 04 1 1 0 04 2 2 0 04 3 3 1 14 4 4 2 2*/

##### Timestep Calculation

Next we must setup a loop to count backwards through time, and at each timestep another loop
to go through all stock price values (except the boundaries). We wish to use two time levels of
storage for V , which we shall call vNew and vOld. Here let V ^{i} be represented by the vector vNew,
and V ^{i+1} be represented by the vector vOld. Then at the end of each timestep we must overwrite
the old values with the new ones, please refer to the solution at the end to see this in action. This
allows us to move recursively through time with only two levels of storage. Since we can
often have timestep constraints, it can be extremely inefficient to store all levels of
time.

Firstly add in your code the timestep loop over i = iMax − 1 down through to i = 0 to step backwards through time. Now inside the loop, you must first implement the boundary condition at j = 0. For example, for a call we may write

and remember that vNew= V ^{i}. Next add in the condition at j = jMax, for a call we may
write

Check at this stage that you are implementing the boundary conditions correctly, so output them to screen. Your code for the loop may look like

// loop through time levels, setting the option price at each grid point, and also on the boundariesfor(int i=iMax-1;i>=0;i--) {// apply boundary condition at S=0vNew[0] = 0.; cout << i << " " << 0 << " " << S[0] << " " << vNew[0] << " " << vOld[0] << endl;// apply boundary condition at S=S_maxvNew[jMax] = S[jMax] - X*exp(-r*(T-i*dt)); cout << i << " " << jMax << " " << S[jMax] << " " << vNew[jMax] << " " << vOld[jMax] << endl; }

Now because we are only using five nodes in space and five nodes in time, this will allow us to check our calculations by hand. Make sure you are checking the values outputed from your code at every step.

Next we can put a loop in between the boundary conditions to calculate for the values 1 ≤ j ≤ jMax − 1. Inside the loop, to make things easier follow the notation from the lectures (click here) and declare variables to store the coefficients A, B and C. Then simply calculate the values of A, B and C for the current value of j, and write

Once all code is checked and working you should run grid checks (change the grid size and check the effect on the solution) to satisfy yourself that the code is correct. You will notice that if you increase jMax by a factor of 10 you will need to increase iMax by a factor of 100 to maintain stability. Please refer to the course notes and lectures for more details on this (click here).

The following code generates the results found in figure 1 and at the end you have estimates
for the option price at time t = 0 for a discrete set of values in S. The value estimates might not
though correspond to the value of S_{0} you need the option price at, so we will need to look at
interpolation in the next section.

#include<iostream>#include<fstream>#include<cmath>#include<vector>usingnamespacestd;/* Template code for the Explicit Finite Difference*/intmain() {// declare and initialise Black Scholes parametersdouble S0=1.,X=2.,T=1.,r=0.05,sigma=0.4;// declare and initialise grid paramatersint iMax=4,jMax=4;// declare and initialise local variables (ds,dt)double S_max=2*X; double dS=S_max/jMax; double dt=T/iMax;// create storage for the stock price and option price (old and new)vector<double>S(jMax+1),vOld(jMax+1),vNew(jMax+1);// setup and initialise the stock pricefor(int j=0;j<=jMax;j++) { S[j] = j*dS; }// setup and initialise the final conditions on the option pricefor(int j=0;j<=jMax;j++) { vOld[j] =max(S[j]-X,0.); vNew[j] =max(S[j]-X,0.); cout << iMax << " " << j << " " << S[j] << " " << vNew[j] << " " << vOld[j] << endl; }// loop through time levels, setting the option price at each grid point, and also on the boundariesfor(int i=iMax-1;i>=0;i--) {// apply boundary condition at S=0vNew[0] = 0.; cout << i << " " << 0 << " " << S[0] << " " << vNew[0] << " " << vOld[0] << endl;for(int j=1;j<=jMax-1;j++) { double A,B,C; A=0.5*sigma*sigma*j*j*dt+0.5*r*j*dt; B=1.-sigma*sigma*j*j*dt; C=0.5*sigma*sigma*j*j*dt-0.5*r*j*dt; vNew[j] = 1./(1.+r*dt)*(A*vOld[j+1]+B*vOld[j]+C*vOld[j-1]); cout << i << " " << j << " " << S[j] << " " << vNew[j] << " " << vOld[j] << endl; }// apply boundary condition at S=S_maxvNew[jMax] = S[jMax] - X*exp(-r*(T-i*dt)); cout << i << " " << jMax << " " << S[jMax] << " " << vNew[jMax] << " " << vOld[jMax] << endl;// set old values to newvOld=vNew; } }

#### Interpolation

##### Example - Linear Interpolation

We show here a simple example of how to take the code detailed above and output a value for the
option at a given stock price using linear interpolation. We use a style here which allows for a
more generic higher order Lagrange polynomial approximation to be derived. We assume that we
have the set of points S_{j} and vNew_{j} and we wish to find the value of the interpolated
function V (S). The linear approximation says that for points j∗ and j ∗ +1 we should
choose

There are three stages:

- Given the value S
_{0}, find j^{∗}such that S_{ 0}∈ [j ∗ ΔS, (j ∗ +1)ΔS]

`int jstar; jstar = S0/dS;` - Evaluate the Lagrange polynomial at S = S
_{0}

`double sum=0.; sum = sum + (S0 - S[jstar+1])/(S[jstar]-S[jstar+1])*vNew[jstar]; sum = sum + (S0 - S[jstar])/(S[jstar+1]-S[jstar])*vNew[jstar+1];` - Test and debug the code!!! You can check that the polynomial above is the equation of
a line passing through the points (S
_{j∗},vNew_{j∗}) and (S_{j∗+1},vNew_{j∗+1}). Try different values of S_{0}and check they look ok.

Once you are convinced the code is working move everything out into a function. You will now be
able to pass in S_{0} as an argument to the function, simply return the value from the
interpolation.

The final code might look like this

#include<iostream>#include<fstream>#include<cmath>#include<vector>usingnamespacestd;/* Template code for the Explicit Finite Difference*/doubleexplicitCallOption(double S0,double X,double T,double r,double sigma,int iMax,int jMax) {// declare and initialise local variables (ds,dt)double S_max=2*X; double dS=S_max/jMax; double dt=T/iMax;// create storage for the stock price and option price (old and new)vector<double>S(jMax+1),vOld(jMax+1),vNew(jMax+1);// setup and initialise the stock pricefor(int j=0;j<=jMax;j++) { S[j] = j*dS; }// setup and initialise the final conditions on the option pricefor(int j=0;j<=jMax;j++) { vOld[j] =max(S[j]-X,0.); vNew[j] =max(S[j]-X,0.); }// loop through time levels, setting the option price at each grid point, and also on the boundariesfor(int i=iMax-1;i>=0;i--) {// apply boundary condition at S=0vNew[0] = 0.;for(int j=1;j<=jMax-1;j++) { double A,B,C; A=0.5*sigma*sigma*j*j*dt+0.5*r*j*dt; B=1.-sigma*sigma*j*j*dt; C=0.5*sigma*sigma*j*j*dt-0.5*r*j*dt; vNew[j] = 1./(1.+r*dt)*(A*vOld[j+1]+B*vOld[j]+C*vOld[j-1]); }// apply boundary condition at S=S_maxvNew[jMax] = S[jMax] - X*exp(-r*(T-i*dt));// set old values to newvOld=vNew; }// get j* such that S_0 \in [ j*dS , (j*+1)dS ]int jstar; jstar = S0/dS; double sum=0.;// run 2 point Lagrange polynomial interpolationsum = sum + (S0 - S[jstar+1])/(S[jstar]-S[jstar+1])*vNew[jstar]; sum = sum + (S0 - S[jstar])/(S[jstar+1]-S[jstar])*vNew[jstar+1];returnsum; } intmain() {// declare and initialise Black Scholes parametersdouble S0=1.639,X=2.,T=1.,r=0.05,sigma=0.4;// declare and initialise grid paramatersint iMax=4,jMax=4; cout <<explicitCallOption(S0,X,T,r,sigma,iMax,jMax) << endl;/* OUTPUT0.194858*/}

## No comments:

## Post a Comment