How Can I Price an Option with a PDE Method in Matlab?


In this article, we build a very simple PDE solver for the Black-Scholes Equation. Using the Finite Volume Discretization Method, we derive the equations required for an efficient implementation in Matlab. The implicit Euler time-stepping of the solver guarantees a stable behavior and convergence.

All posts in this series:

Discretizing the Black-Scholes PDE using a Finite Volume Method

The Black-Scholes PDE is a Cauchy-Problem in backwards time where the initial values are given by the payoff at maturity. We denote the asset price by S, the value of the option by V, the volatility of the asset by \sigma and the risk free rate by r. Furthermore, a partial derivative is denoted by a subscript as usual for numerical PDE solutions. The solver is an implicit finite volume discretization of

V_{\tau}=\frac{\sigma^2}{2}S^2V_{SS}+rSV_S-rV (1)

working backwards in time from maturity to present time t_0. This equation is defined using the backwards time \tau=T-t. We integrate equation (1) over a finite volume A_i using the discrete values S_i, where


That means that the cell boundaries are placed half way between the nodes at



After the integration of equation (1) over the ith cell, we obtain


In the following, we denote the value of the option at time \tau_n and asset price S_i as


Using approximations we get





Furthermore, we get for the other terms

\int_{A_i}rVds\approx rV_iA_i

\int_{A_i}rSV_S ds

\approx r S_i\left[V_{i+1/2}-V_{i-1/2}\right]


On the left hand side of equation (1), we get

\int_{A_i}V_{\tau}ds\approx A_i\left[\frac{V_i^{n+1}-V_i^n}{\Delta\tau}\right].

All the above equations finally lead to


where we have to define the time level of the right hand side.

Implicit Euler Discretization

Choosing a time level n for the right hand side, we get a so-called explicit discretization; choosing a level of n+1, the discretization is called implicit. The implicit discretization is more stable so that we will use

A_i\left[\frac{V_i^{n+1}-V_i^n}{\Delta\tau}\right]=\frac{\sigma^2}{2}S_i^2\left(\frac{V^{n+1}_{i+1}-V^{n+1}_i}{S_{i+1}-S_i}+\frac{V^{n+1}_{i-1}-V^{n+1}_i}{S_i-S_{i-1}}\right)+rS_i\left[\frac{V^{n+1}_{i+1}-V^{n+1}_{i-1}}{2}\right]-rV^{n+1}_iA_i (2)

in the following. The final algorithm of the PDE solver is obtained by rearranging equation (2) in a matrix M such that linear equation

M\cdot V^{n+1}=V^n (3)

with the known vector V^n=(V^n_1\ldots V^n_{m}) and the unknown V^{n+1}=(V^n_1\ldots V^n_{m}) can be solved by Matlab.

The boundary, i.e. the first and the last element of V can be assumed as

V_m^i =0~~\forall i

Translation into M-Code

Now, we write a program which implements the PDE solver. We use equally distant time steps in \tau and a structured mesh in the S direction.

In this implementation, we compute the price of a put option with S=100, K=100T=2r=0.05 and \sigma=0.4. Furthermore, we use a start mesh with

S = [0:20:60 65:5:90 92:1:110 115:5:150 160:20:200 300 500 1000];

and set the initial values at maturity time to

V = max(K-S',0);" % value of V at maturity time

T, i.e. \tau=0.

Generation of the mesh

We use a resolution of 200 time steps. Then, we use an algorithm for refine the mesh and the time steps such that we can easily control the accuracy by changing a single parameter: “refine”.

% Pricing a European Option with a PDE method (finite volumes)

% init S with a structured mesh
refine = 4;
S = [0:20:60 65:5:90 92:1:110 115:5:150 160:20:200 300 500 1000];
for i=1:refine-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;

Start parameters

% define parameters for Option
sigma = 0.4;
r = 0.05;
K = 100;
T = 2;
time_steps = round(400 *2^(refine-1)* T);

% calculate size of time step
dtau = T/time_steps;

% V(t=T) = payoff at maturity
V = zeros(length(S),time_steps+1);
V(:,1) = max(K-S',0);

% Define values for S_i_minus, S_i and S_i_plus for convenient access
Sm = S(1:length(S)-2);
Si = S(2:length(S)-1);
Sp = S(3:length(S));

Assemble discretization matrix

Rearranging equation (2), we find that matrix M is a tri-diagonal matrix with the following properties:

% compute values for m_i_minus, m_i and m_i_plus according to the formula given
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 condition

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

% assemble matrix M from (sub-) diagonals
M = diag(m_i) + diag(m_iminus,-1) + diag(m_iplus,1);

Solve for all time-steps

Now, the solver according to equation (3) is simple. It solves for V^{n+1} at each time step n:

% solve for V(t_0)
for tau=1:time_steps
 V(:,tau+1) = M\V(:,tau);

Analyze results

Finally, we use Matlab to present the computations:

%% plot the option values
xlabel('asset price S','FontWeight','bold');
zlabel('option value V','FontWeight','bold');

And we can zoom in on S=100, the strike price of the option:

Values of European Put Option computed using a PDE solver

Values of European Put Option computed using a PDE solver

index= find((S>60) & (S<120));
title('V(S,t) zoomed on S in [60,120]','FontWeight','bold');
xlabel('asset price S','FontWeight','bold');
zlabel('option value V','FontWeight','bold');
Values of European Put Option computed using a PDE solver

Values of European Put Option computed using the PDE solver zoomed on the option's strike price.

We can analyze the hedge parameter delta, V_S, the first derivative of V w.r.t. S.

figure3 = figure(3);
axes4 = axes('Parent',figure3,'YGrid','on','XGrid','on','FontWeight','bold','FontSize',12);
title('Delta: V_S(S,t_0)','FontWeight','bold');
xlabel('asset price S','FontWeight','bold');
ylabel('option value V','FontWeight','bold');

European Option delta computed by PDE solver

We can analyze the hedge parameter gamma, V_{SS}, the second derivative of V w.r.t. S.

figure4 = figure(4);
axes4 = axes('Parent',figure4,'YGrid','on','XGrid','on','FontWeight','bold','FontSize',12);
title('Gamma V_S_S(S,t_0)','FontWeight','bold');
xlabel('asset price S','FontWeight','bold');
ylabel('option gamma V_S_S','FontWeight','bold');


This article presents a simple yet efficient PDE solver for the Black-Scholes Equation. This solver uses a finite volume discretization and an implicit Euler time-stepping.

Update 15.4.2012:

  • Improved accuracy for V_{SS}
  • Changed formula for time_steps to represent convergence
  • General code cleaning

Leave a Reply

Fill in your details below or click an icon to log in: Logo

You are commenting using your account. Log Out /  Change )

Facebook photo

You are commenting using your Facebook account. Log Out /  Change )

Connecting to %s

This site uses Akismet to reduce spam. Learn how your comment data is processed.