Next: Simulation Up: Finite Differences Previous: European Options.   Contents   Index

Subsections

American Options.

We now compare the American versions of the same algoritms, the only difference being the check for exercise at each point.

Computer algoritm, implicit finite differences.

```#include <cmath>
#include "newmat.h"  // definitions for newmat matrix library
#include <vector>
#include <algorithm>

double option_price_put_american_finite_diff_implicit(
double S, double X, double  r, double sigma, double time,
int no_S_steps, int no_t_steps)
{
double sigma_sqr = sigma*sigma;
// need no_S_steps to be even:
int M; if ((no_S_steps%2)==1) { M=no_S_steps+1; } else { M=no_S_steps; };
double delta_S = 2.0*S/M;
double S_values[M+1];
for (int m=0;m<=M;m++) { S_values[m] = m*delta_S; };
int N=no_t_steps;
double delta_t = time/N;

BandMatrix A(M+1,1,1); A=0.0;
A.element(0,0) = 1.0;
for (int j=1;j<M;++j) {
A.element(j,j-1) = 0.5*j*delta_t*(r-sigma_sqr*j);    // a[j]
A.element(j,j)   = 1.0 + delta_t*(r+sigma_sqr*j*j);  // b[j];
A.element(j,j+1) = 0.5*j*delta_t*(-r-sigma_sqr*j);   // c[j];
};
A.element(M,M)=1.0;
ColumnVector B(M+1);
for (int m=0;m<=M;++m){ B.element(m) = max(0.0,X-S_values[m]); };
ColumnVector F=A.i()*B;
for(int t=N-1;t>0;--t) {
B = F;
F = A.i()*B;
for (int m=1;m<M;++m) { 	// now check for exercise
F.element(m) = max(F.element(m), X-S_values[m]);
};
};
return F.element(M/2);
};
```

Computer algoritm, explicit finite differences.

```// file: findiff_exp_am_put.cc
// author: Bernt A Oedegaard

#include <cmath>
#include <algorithm>
#include <vector>

double option_price_put_american_finite_diff_explicit( double S,
double X,
double r,
double sigma,
double time,
int no_S_steps,
int no_t_steps) {
double sigma_sqr = sigma*sigma;

int M;           // need M = no_S_steps to be even:
if ((no_S_steps%2)==1) { M=no_S_steps+1; } else { M=no_S_steps; };
double delta_S = 2.0*S/M;
vector<double> S_values(M+1);
for (int m=0;m<=M;m++) { S_values[m] = m*delta_S; };
int N=no_t_steps;
double delta_t = time/N;

vector<double> a(M);
vector<double> b(M);
vector<double> c(M);
double r1=1.0/(1.0+r*delta_t);
double r2=delta_t/(1.0+r*delta_t);
for (int j=1;j<M;j++){
a[j] = r2*0.5*j*(-r+sigma_sqr*j);
b[j] = r1*(1.0-sigma_sqr*j*j*delta_t);
c[j] = r2*0.5*j*(r+sigma_sqr*j);
};
vector<double> f_next(M+1);
for (int m=0;m<=M;++m) { f_next[m]=max(0.0,X-S_values[m]); };
vector<double> f(M+1);
for (int t=N-1;t>=0;--t) {
f[0]=X;
for (int m=1;m<M;++m) {
f[m]=a[m]*f_next[m-1]+b[m]*f_next[m]+c[m]*f_next[m+1];
f[m] = max(f[m],X-S_values[m]);  // check for exercise
};
f[M] = 0;
for (int m=0;m<=M;++m) { f_next[m] = f[m]; };
};
return f[M/2];
}
```