An American option is an option which the owner can exercise at any time during its lifetime. That means the option’s value cannot drop below the exercise value, i.e. the option value  of an American put option satisfies

. (1)

We use the above condition in the PDE solver (How Can I Price an Option with a PDE Method in Matlab?) to price an American put option with , , ,  and .

All posts in this series:

• Basics of a PDE solver in Matlab
• Pricing American options with a PDE solver
• Efficient pricing of barrier options

Implementation of the PDE Solver

Introducing the constraint (1) can be done by adding a single line of code within the solver loop:

% solve for V(t_0)
for tau=1:time_steps
  V(:,tau+1) = M\V(:,tau);
  V = max(K-S',V);
end

Analyzing the convergence of the solver, we will increase “refine” from 1 to 9 and then take a close look at the results. The complete solver with iteration of the refine is:

% init S with a structured mesh
clear;
refine = 9;

sigma = 0.4; r = 0.05; K = 100; T = 1;

S = [0 60:20:160 200 1000];

Sn = [];
for iterRefine=1:refine
 % refine mesh
 if iterRefine>1
   Sn(1:2:length(S)*2) = S;
   Sn(2:2:length(S)*2-1) = (S(2:length(S)) -S(1:length(S)-1))/2 + S(1:length(S)-1);
   S = Sn;
 end

 time_steps = round(256 *2^(iterRefine-1) * T);
 dtau = T/time_steps;

 V = max(K-S',0); % V(t=T) = payoff at maturity

 % Abbreviation
 Sm = S(1:end-2);
 Si = S(2:end-1);
 Sp = S(3:end);

 m_iminus = -(sigma^2*Si.^2*dtau)./((Sp-Sm).*(Si-Sm)) + Si.*r*dtau./(Sp-Sm);

 m_i = 1 + (sigma^2*Si.^2*dtau).*(1./( (Sp-Sm).*(Sp-Si) ) +1./( (Sp-Sm).*(Si-Sm) ) ) + r*dtau;
 m_i = [1+r*dtau m_i 1]; % boundary

 m_iplus =-(sigma^2*Si.^2*dtau)./((Sp-Sm).*(Sp-Si)) - Si.*r*dtau./(Sp-Sm);

 % Same but more efficient than: M = diag(m_iminus,-1) + diag(m_i) + diag(m_iplus,1);
 M = spdiags([ [m_iminus 0 0]; m_i; [0 0 m_iplus]]', -1:1, length(m_i), length(m_i));

 for tau=1:time_steps
   V = M\V;
   V = max(K-S',V);
 end

 % save value at S=100 for convergence study
 V_res(iterRefine,1) = V(S==100);
end

On my laptop, this code runs about 30 sec. Then, the option prices for for the different refine values are in V_res. An increase of the refinement doubles the nodes in the mesh for S and doubles the number of time steps. So, what happens to the error?

 
 refine        V_res  Difference       Factor
      1      13.1160           -            -
      2      13.5188      0.4028            -
      3      13.6265      0.1077       3.7399
      4      13.6555      0.0290       3.7185
      5      13.6637      0.0083       3.5074
      6      13.6663      0.0026       3.2000
      7      13.6671      0.0008       3.0796
      8      13.6674      0.0003       2.9366
      9      13.6675      0.0001       2.6696

In the above table, we see that increasing the “refine” by one (doubling the nodes in S mesh and doubling the time steps) reduces the error to about 1/3.

Possible further improvements

Instead of the implicit Euler time-stepping, one can increase the convergence combining implicit and explicit Euler to Crank-Nickolson. Then, using Crank-Nickolson and apply the constraint (1) in an implicit way, the increasing “refine” by one reduces the error to 1/4. That means, highly accurate results are obtained in less time (see e.g. Peter Forsyth: Quadratic convergence of a penalty method for valuing American options).

Conclusion

It is easy to impose the early exercise constraint in a PDE solver. The presented implementation only requires a single extra line of code compared to the European Option PDE solver.