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

$A_i=\left(\frac{S_{i+1}+S_i}{2}\right)-\left(\frac{S_{i-1}+S_i}{2}\right)$ .

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

$S_{i+1/2}=\frac{S_i+S_{i+1}}{2}$

$S_{i-1/2}=\frac{S_i+S_{i-1}}{2}$

After the integration of equation (1) over the $i$ th cell, we obtain

$\int_{A_i}V_{\tau}ds=\int_{A_i}\frac{\sigma^2}{2}S^2V_{SS}ds+\int_{A_i}rSV_Sds-\int_{A_i}rVds$ .

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

$V(S_i,\tau_n)=V_i^n$ .

Using approximations we get

$\int_{A_i}\frac{\sigma^2}{2}S^2V_{SS}ds$

$\approx\frac{\sigma^2}{2}S^2_i\int_{A_i}V_{SS}ds$

$=\frac{\sigma^2}{2}S^2\left((V_S)_{i+1/2}-(V_S)_{i-1/2}\right)$

$\approx\frac{\sigma^2}{2}S^2\left(\frac{V_{i+1}-V_i}{S_{i+1}-S_i}+\frac{V_{i-1}-V_i}{S_i-S_{i-1}}\right)$

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]$

$=rS_i\left[\frac{V_{i+1}-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

$A_i\left[\frac{V_i^{n+1}-V_i^n}{\Delta\tau}\right]=\frac{\sigma^2}{2}S_i^2\left(\frac{V_{i+1}-V_i}{S_{i+1}-S_i}+\frac{V_{i-1}-V_i}{S_i-S_{i-1}}\right)+rS_i\left[\frac{V_{i+1}-V_{i-1}}{2}\right]-rV_iA_i$

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_0^{n+1}=V_0^n(1-r\Delta\tau)$
$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=100$ , $T=2$ , $r=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
clear;
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;
end

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);
end

Analyze results

Finally, we use Matlab to present the computations:

%% plot the option values
figure(1);
mesh(S,0:dtau:T,V');
title('V(S,t)','FontWeight','bold');
xlabel('asset price S','FontWeight','bold');
ylabel('time','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

figure(2);
index= find((S>60) & (S<120));
mesh(S(index),0:dtau:T,V(index,:)');
title('V(S,t) zoomed on S in [60,120]','FontWeight','bold');
xlabel('asset price S','FontWeight','bold');
ylabel('time','FontWeight','bold');
zlabel('option value V','FontWeight','bold');

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$ .

V_S=(V(3:end,time_steps+1)-V(1:end-2,time_steps+1))./(Sp-Sm)';
figure3 = figure(3);
axes4 = axes('Parent',figure3,'YGrid','on','XGrid','on','FontWeight','bold','FontSize',12);
box(axes4,'on');
hold(axes4,'all');
plot(S(2:length(S)-1),V_S,'LineWidth',2);
title('Delta: V_S(S,t_0)','FontWeight','bold');
xlabel('asset price S','FontWeight','bold');
ylabel('option value V','FontWeight','bold');

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

V_SS=2*(V_S(3:end)-V_S(1:end-2))./(Sp(3:end)-Sm(1:end-2))';
figure4 = figure(4);
axes4 = axes('Parent',figure4,'YGrid','on','XGrid','on','FontWeight','bold','FontSize',12);
box(axes4,'on');
hold(axes4,'all');
plot(S(3:length(S)-2),V_SS,'LineWidth',2);
title('Gamma V_S_S(S,t_0)','FontWeight','bold');
xlabel('asset price S','FontWeight','bold');
ylabel('option gamma V_S_S','FontWeight','bold');

Conclusion

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

Computer Aided Finance – Excel, Matlab, Theta Suite etc.

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

All posts in this series:

Discretizing the Black-Scholes PDE using a Finite Volume Method

Implicit Euler Discretization

Translation into M-Code

Generation of the mesh

Start parameters

Assemble discretization matrix

Solve for all time-steps

Analyze results

Conclusion

Update 15.4.2012:

Leave a Reply Cancel reply

All posts in this series:

Discretizing the Black-Scholes PDE using a Finite Volume Method

Implicit Euler Discretization

Translation into M-Code

Generation of the mesh

Start parameters

Assemble discretization matrix

Solve for all time-steps

Analyze results

Conclusion

Update 15.4.2012:

Related

Leave a Reply Cancel reply