Using ListConvolve for the Explicite Finite Difference Method

(c) Weber & Partner, financial technology, www.weber-und-partner.com

Developing an explicite finite difference method for the Black-Scholes equation

The partial differential equation from Black-Scholes for a derivate V on the security S is given as

pde = -(∂_t) V[S, t] == 
          (1/2) (σ^2) (S^2) (∂_(S, S))V[S, t] + 
            r S (∂_S) V[S, t] - r V[S, t];

Numerical approximation for the partial derivates follow for example from the Taylor series. The parital derivates in the equation above are replaced through the appropriate Taylor series

\(\(rls\  = \ {∂\_t V[S, 
                t] -> \ \(V[S, t - δt] - V[S, t]\)\/δt, 
	∂\_S V[S, 
                t] -> \ \(V[S + δS, t] - V[S - δS, t]\)\/\(2  δS\), 
\ \ \ \ \ \ \ \ \ ∂\_\(S, S\)V[S, 
                t] -> \ \(V[S + δS, t] - 2\ V[S, t] + V[S - δS, t]\)\/δS\^2};\)\)
prox  =  pde  /.  rls

In the next step the notation is changed to make it more consistent with a grid scheme.

prox  =  
      prox  /. {V[S, t]  ->  (Subscript[V, j, i]), 
          V[S, t - δt]  ->  (Subscript[V, j, i - 1]), 
          V[S + δS, t]  ->  (Subscript[V, j + 1, i]), 
          V[S - δS, t]  ->  (Subscript[V, j - 1, i])}

To better illustrate the structure of the equation some more notational adjustments are made. The new structur will later help to simplify the calculations

prox  =  
      prox  /.   {(((Subscript[V, j + 1, i]) - 2 (Subscript[V, j, i]) + (Subscript[V, j - 1, i]))/δ)(S^2) ->  (Subscript[Γ, i, j]), (((Subscript[V, j + 1, i]) - (Subscript[V, j - 1, i]))/(2  δS)) ->  (Subscript[Δ, i, j])}

Solving the last expression for (Subscript[V, i, j]) and simplifying leads to

diff  =  Solve[prox, (Subscript[V, j, i - 1])]  //  Simplify  //  
          First;
diff  //  TraditionalForm

The last expression can be rewritten as

\(TraditionalForm\`V\_\(j, i - 1\) = 
      V\_\(j, i\) + δt\ \((\(-\(1\/2\)\)\ S\^2\ Γ\_\(i, j\)\ σ\^2 + r\ V\_\(i, j\) - 
              r\ S\ Δ\_\(i, j\))\)\)

with \(TraditionalForm\`Γ\_j\%i = \(V\_\(j + 1\)\%i - \ 2  V\_j\%i + \ V\_\(j - 1\)\%i\)\/δS\^2\)and \(TraditionalForm\`Δ\_j\%i\  = \(V\_\(j + 1\)\%i - \ V\_\(j - 1\)\%i\)\/\(2  δS\)\).
The presented calculation scheme is a version of the explicit finite difference method. While applying this calculation schemebei the new values for the derivate (Subscript[V, j, i - 1])are stepwise calculated from (Subscript[V, j + 1, i]),  (Subscript[V, j, i]) and (Subscript[V, j - 1, i]). Picture  illustrates this approach.

Implementing the numerical method with ListConvolve

An efficient Mathematica function for the calculation of the differences needed in Δ and Γ is available through ListConvolve. To demonstrate this  ListConvolve is applied to a list of symbols.

Clear[V, δS];
v =  Table[(Subscript[V, j]), {j, 6}]

ListConvolve used for  Δ  results in the following expression

Δ  = ListConvolve[{1, 0, -1}, 
              v]/(2 δS) // TraditionalForm

The first list in ListConvolve, the kernel {1,0,1}, is applied piecewise to the second list and multiplicates the elements the elements of the second list and adds them up according to the values given in the kernel. This operation runs internally in Mathematica and is much faster as any loop written in Mathematica code.

The approach used for Δ can also be applied for the calculation of  Γ. ListConvolve can therefore replace loops which are common to the explicit approximation of partial differential equations.

Γ  = 
      ListConvolve[{1, -2, 1}, v]/δ(S^2)  //  
        TraditionalForm

To complete the specification of the finite difference method for the numerical solution of partial differential equations boundary and initial conditions are needed. These result from the given valuation problem.

Example: Numerical valuation of a call option

To illustrate the application of the method the following call option will be valued

S0  =  100. ;     
X  =  100. ;       
σ  =  0.2;         
τ  =  1.0;         
r    =  0.05;

The number of grid points for the discrete prices of the stock n can freely be chosen within a specific range. Increasing the number of time steps improves the accuracy but also increases the overall calculation time. For a first demonstration the number of discrete stock prices is set to 20.

n  =  20;

The grid points for the stock price should be placed in a range not to tight around the current stock price.  In this example the reange is choosen from zero  up to twice the strike price. From the chosen region results the step size ΔS for the discretization of the stock prices range. One way to generate the list of grid points is to use NestList .  #+ΔS& within NestList is a generic function defined for local use.
On the list of discrete stock prices the exercise function of the option can be applied. The resulting list provides the starting or initial values for the numerical method.

δS  = (( 2 X)/n);
S  =  NestList[# + δS &, 0, n];
V  =  Map[Max[# - X, 0] &, S];

The necessary number of time steps for the explicite FDM to converge depends on the step size for the discretization of the stock price, the volatility and the strike price (see Willmot [1998], p. 629). Using the formula given in  Willmot [1998] the number of time steps can be calculated as

nt  =  
        Floor[τ/(((δ(S/(2 X σ))))^2)] + 
          1;

Then the size of the time steps are

δt  =  (τ/nt);

For the calculation of Γ and Δ the Mathematica function ListConvolve is used as explained above. (Attention: Before executing the following Mathematica code make sure that V is initiated with numbers.)

Do[
   Δ  = 
        ListConvolve[{1, 0, -1}, V]/(2 δS); 
   Γ  = 
        ListConvolve[{1, -2, 1}, V]/δ(S^2); 
   s  =  Take[S, {2, -2}];  
   v  =  
        Take[V, {2, -2}]; 
   V  = 
        Join[{0},  
          v  + δt ((1/2) (σ^2) (s^2) Γ +  
                  r s Δ - r v), {Last[S] - 
              X (E^(-r i δt))}], {i, 
        nt}]

After the loop is finished V contains a list of option values. Each option value corresponds to a discrete stock price on the grid. Interplation on these numbers can be used to calculate the option price for the current price of the underlying (Subscript[S, 0]).

V0  =  Interpolation[Transpose[{S, V}]];
V0[S0]

Literature

Clewlow, L./Strickland, C., 1998, Implementing Derivatives Models, John Wiley & Sons Ltd.

Willmot, P., 1998,  Derivatives: The theory and practice of financial engineering, John Wiley & Sons Ltd.