Example - European Put

Now, assuming that our vectors are indexed from 0 to $N$ (and therefore have $N+1$ elements) I will go through an example for a European put with $\sigma=0.4$, $r=0.05$, $X=2$, $dS=1$, $dt=0.25$, and $N=4$. Note here that these calculations could have been done by hand, and trying to replicate an example of this scale with a code is a good way to start, and also to check for bugs/errors.

Matrix Setup

I tend to use a, b and c to represent the tridagonal matrix $A$, and d for the right hand side of the equation. Just remember that you must be consitent with your notation.

First setup the boundary condition

$$V(S=0) = X e^{-r(T-t)}$$

by writing the code

a[0] = 0.;   // this term multiplies V_new[-1] in the equation...
b[0] = 1.;   // this term multiplies V_new[0] in the equation...
c[0] = 0.;   // this term multiplies V_new[1] in the equation...
d[0] =   // fill this in...


Next write a for loop to assign each of the inner equations. Here we shall let $V_{new}[j]=V^i_j$ and $V_{old}[j]=V^{i+1}_j$. The inner equations may be written

$$a_jV^i_{j-1}+b_jV^i_{j}+c_j V^i_{j+1} = d_j$$

so inside the loop write something like:

a[j] =   // this term multiplies V_new[j-1] in the equation...
b[j] =   // this term multiplies V_new[j] in the equation...
c[j] =   // this term multiplies V_new[j+1] in the equation...
d[j] =   // this term will involve V_old...


After the loop setup the final boundary condition:

a[N] = 0.;   // this term multiplies V_new[N-1] in the equation...
b[N] = 1.;   // this term multiplies V_new[N] in the equation...
c[N] = 0.;   // this term multiplies V_new[N+1] in the equation...
d[N] =   // fill this in...


Now we should have finished setting up the matrix equations, they should look something like (at the first timestep)

j	a[j]	b[j]	c[j]	d[j]
0	0	1	0	1.97516
1	0.0275	-4.105	0.0525	-3.95
2	0.135	-4.345	0.185	-0.135
3	0.3225	-4.745	0.3975	-0
4	0	1	0	0

Matrix Solver

Still inside the time loop, we now need to solve the matrix equation. This can be done using either iterative methods (SOR) or direct methods (Thomas's solver). We shall go through Thomas's solver here but this cannot be used on American options so the SOR method is left as an exercise (it should be easier).

There are two stages in the direct solving of this matrix equation:

• Reduction to upper triangular matrix,
• Back-substitution to solve.

Reduction

In order to reduce the matrix to upper triangular, we must eliminate all $a[j]$ terms. So, use a for loop moving from the second equation ($j=1$) to the last updating each equation to eliminate $a[j]$. Inside the loop should look something like:

   // use a[j] and b[j-1] to eliminate a[j]...
b[j] = b[j] -   // update value of b[j]...
   // c[j] does not change...
d[j] = d[j] -   // update value of d[j]...
a[j] = 0.;   // can set this to zero at the end if you want, but not used again...


and carrying on from before the reduced matrix should be:

j	a[j]	b[j]	c[j]	d[j]
0	0	1	0	1.97516
1	0	-4.105	0.0525	-4.00432
2	0	-4.34327	0.185	-0.266689
3	0	-4.73126	0.3975	-0.0198024
4	0	1	0	0

Back-Substitute

We are now left we an upper triangulr matrix. In order to solve and find $V_{new}$ we must back-substitute. First calculate the value of $V_{new}[n]$ from the final equation:

   // find V_new[N] from equation " b[N] V_new[N] = d[N] "
V_new[N] =   // fill this in...


then use a for loop to iterate backwards from the $N-1$th equation back up to the 0th equation, finding the value of $V_{new}$ at each point. Inside the loop should look something like:

   // find V_new[j] from equation " b[j] V_new[j] + c[j] V_new[j+1] = d[j] "
V_new[j] = (d[j] -   // fill this in...


At the end of the loop the value of $V_{new}$ and $V_{old}$ should be:

i	j	V_old[j]	V_new[j]
3	0	2	1.97516
3	1	1	0.976261
3	2	0	0.061581
3	3	0	0.00418543
3	4	0	0

All that is left is to set the old values equal to the new. Done!!

Iterative methods

• The SOR method must be implemented at the solve stage, everything preceeding that may be kept the same
• Since both methods solve exactly the same equation, they should both post exactly the same results, provided the tolerance in the SOR method is taken small enough.
• Therefore the results below can be used to check your iterative method as well as your direct method.

Results

The full set of results for the example is below...

i	j	V_old[j]	V_new[j]
3	0	2	1.97516
3	1	1	0.976261
3	2	0	0.061581
3	3	0	0.00418543
3	4	0	0
2	0	1.97516	1.95062
2	1	0.976261	0.954845
2	2	0.061581	0.112606
2	3	0.00418543	0.01471
2	4	0	0
1	0	1.95062	1.92639
1	1	0.954845	0.935397
1	2	0.112606	0.155285
1	3	0.01471	0.0282984
1	4	0	0
0	0	1.92639	1.90246
0	1	0.935397	0.917626
0	2	0.155285	0.191233
0	3	0.0282984	0.0429639
0	4	0	0