How Do Efficient PDE Solvers for Barrier Options Look Like?

· Pricing and Hedging Tutoriual

Besides European and American Options, another challenge in option pricing is the valuation of Barrier Options. We will see that simply applying the algorithms from the previous posts does not converge well. Especially, pricing a long-term up-and-out barrier option is hard, due to the discontinuity of the payoff. 

All posts in this series:

Barrier Options

In this post, we will take a look at the following barrier option:


  • Stock price S – Strike K if stock price was never above barrier level H.
  • 0 otherwise.

That means, the value V for S>H is zero. First, we try to introduce this in an explicit constraint. We use the  condition in the PDE solver (How Can I Price an Option with a PDE Method in Matlab?) and change the solver loop:

for tau=1:time_steps
   V = M\ V;

Analyzing the convergence of the solver with explicit constraint

Now, we take a look at a numerical example. In the following, we price an up-and-out Barrier option with S=100K=100, H=120, T=2r=0.05 and \sigma=0.4. Using different levels of “refine” with the PDE solver and the explicit constraint, we obtain the following values:

 refine        V_res   Difference       Factor
      1   0.05571209            0            0
      2   0.07059614   0.01488404            0
      3   0.07385265   0.00325651   4.57054309
      4   0.07434690   0.00049424   6.58891052
      5   0.07402224  -0.00032465  -1.52237808
      6   0.07343880  -0.00058344   0.55643959
      7   0.07281298  -0.00062581   0.93229160

This table presents the option values for different levels of “refine”. The column “Difference” presents the difference between two refine levels and the column “Factor” presents the factor between the “Difference” of two levels of refine.

We can see that absolute values of factor changes between 0.5 and 6.5, i.e. the error does not decrease smoothly for higher levels.

Improving convergence using an implicit barrier constraint

The idea for improving the convergence is an implicit introduction of the constraint, similar to the penalty method for American options. I.e. we introduce an additional equation into

M\cdot V^{n+1}=V^n.

We want to force numerically V^{n+1}_i = 0 ~~\forall i:S_i>H. This can be done e.g. by the introduction of a penalty in the equation.

100000\cdot V^{n+1}_i =V^{n}_i ~~\forall i:S_i>H

Then, the solver loop is:

    for tau=1:time_steps
        V = (M + diag(P))\V

Analyzing the convergence of the solver with implicit constraint

Analyzing the convergence of the solver, we increase “refine” from 1 to 7 and then take a close look at the results.

 refine        V_res   Difference       Factor
      1   0.05269513            0            0
      2   0.06639892   0.01370379            0
      3   0.06940017   0.00300124   4.56602960
      4   0.07011419   0.00071402   4.20327399
      5   0.07028354   0.00016934   4.21631141
      6   0.07032179   0.00003824   4.42771568
      7   0.07032934   0.00000754   5.06592917

In this table, the factor between different refine levels is between 4 and 5. This is a good result: It shows that doubling the nodes in the mesh for S and doubling the number of time steps reduces the error to 1/4. This is called quadratic convergence. By the way, the real option value can be obtained from an analytic solution of the Barrier Option is V_{analytic}=0.07032913.


Here is the Matlab PDE solver for Barrier options:

% init S with a structured mesh
refine = 7;
explicit = false; % change to true for explicit handling of constraint

sigma = 0.4;
r = 0.05;
K = 100;
T = 2;
H = 120;
LargePenalty = 100000;

S = [0 60:20:160 200 1000];
for iterRefine=1:refine
    if 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;

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

    V = max(S'-K,0); % V(t=T) = payoff at maturity
    V(S>=H) = 0;     % no value above barrier

    Sm = S(1:length(S)-2);
    Si = S(2:length(S)-1);
    Sp = S(3:length(S));

    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);
    m_iminus = -(sigma^2*Si.^2*dtau)./((Sp-Sm).*(Si-Sm)) + Si.*r*dtau./(Sp-Sm);

    M = spdiags([ [m_iminus 0 0]; m_i; [0 0 m_iplus]]', -1:1, length(m_i), length(m_i));
    P = zeros(length(S),1);

    for tau=1:time_steps
        if explicit
           V = M\ V;
           % equivalent but more efficient than V = (M + diag(P))\V
           V = (M + spdiags(P,0,length(V),length(V))) \ V;
    % save value at S=100 for convergence study
    V_res(iterRefine,1) = V(S==100);

V_diff   = V_res(2:end)-V_res(1:end-1);
V_factor = V_diff(1:end-1)./V_diff(2:end);

V_diff   = [0; V_diff];
V_factor = [0; 0; V_factor];

indexRefine = [1:refine]';
disp([indexRefine V_res V_diff V_factor])


Using a penalty method, the PDE solver can efficiently handle barrier options. The numerical example shows that an up-and-out barrier option can be priced with quadratic convergence.

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.