Finite Difference Methods
We learned last week that the price of a call option in the Black-Scholes models satisfies a partial differential equation (PDE) called the Black-Scholes PDE.
Once you know this you can compute the price of a call option by solving this PDE. Solving PDEs has been a big topic in mathematics since the 19th Century. For example Fourier and Poincaré are just two of the many famous mathematicians to have studied PDEs.
Once Black and Scholes reduced their problem to solving a PDE they could just appeal to known results on PDEs to compute the price.
Copyright By PowCoder代写 加微信 powcoder
In this topic we will solve the Black and Scholes PDE numerically by a “brute force” approach called a finite difference method.
Plan of campaign
• We will start by looking at how to approximate derivatives numerically
• We will then see how to solve ordinary differential equations (ODE) numerically • We will then see how to solve the Black-Scholes PDE numerically
Solving ODEs numerically isn’t quite as straightforward as you might expect. There is more than one method available corresponding to three ways of estimating derivatives.
• The Euler scheme based on a forward difference is an intuitive approach but can be unstable. This in effect means it gives terrible answers.
• A backward difference approach to can be more stable and so give better results.
• Using central differences is an even better approach as it is stable and converges more rapidly.
• Problem: Convergence of finite difference methods
• Key Paper: “Survey of the Stability of Linear Finite Difference Equations” • Keywords: Convergence, Stability
1 Approximating Derivatives
One of the simplest numerical methods is the forward difference method for approximating derivatives. If f : R → R is a differentiable function, one can estimate the derivative f′(x) using
f′(x)≈ f(x+h)−f(x) h
for a small value of h. By definition of the derivative
f′(x)= lim f(x+h)−f(x)
So this method converges.
This isn’t the only method to estimate derivatives. Two alternatives are the backward difference estimate and the central difference estimate.
Forward Backward Central
f′(x) ≈ f(x+h)−f(x) O(h) h
f′(x) ≈ f(x)−f(x−h) O(h)
f′(x) ≈ f(x+h)−f(x−h) O(h2)
Name Estimate Accuracy
If we choose f (x) with a known derivative at some point x, we can then plot a log-log plot of the error in the estimate against the true value.
I’ve taken f (x) = sin(x) and x = 1.2
Theorem: Suppose that f (x) is three times continuously differentiable on some interval [x − ε, x + ε] then f′(x) = f(x+h)−f(x−h) +O(h2)
for h ≤ ε.
Proof: Without loss of generality, we may take x = 0. By Taylor’s Theorem with remainder, we have that
f(h) = f(0) + f′(0)h + 1f′′(0)h2 + f′′′(ξ)h3 26
− f′(0) = 2h f(0) + f′(0)h + 2f′′(0)h2 + 6f′′′(ξh)h3
11 −f(0) − f′(0)(−h) − 2f′′(0)(−h)2 − 6f′′′(ξ−h))(−h)3
Since f ′′′ is continuous on the closed interval [x − ε, x + ε] it is bounded. So there is some constant C such 4
for some ξh with |ξh| < |h|. f(0+h)−f(0−h)
− f′(0) 11
= 12f′′′(ξh) + 12f′′′(ξ−h) h2
f (0 + h) − f (0 − h) 2h
− f′(0) ≤ Ch2
for all h < ε.
1.1 Higher Order Derivatives
One can approximate f′′(x) using the formula
f′′(x) ≈ f(x+h)−2f(x)+f(x−h).
h2 A heuristic proof is that by the forward estimate for f′:
f′(x) ≈ f(x+h)−f(x), f′(x−h) ≈ f(x)−f(x−h) hh
The backward estimate for f′′ is
2 Exercises 2.1 Exercise
f′′(x) ≈ f′(x) − f′(x − h)
f(x+h)−f(x) − f(x)−f(x−h) ≈hh
= f(x+h)−2f(x)+f(x−h)
Suppose f is smooth (i.e. infinitely differentiable) on some interval [x−ε, x+ε]. By considering a numerical example, find out the largest value of n for which this identity holds.
f′(x) = −f(x+2h)+8f(x+h)−8f(x−h)+f(x−2h) + O(hn) 12h
for h < ε. Store the value in a variable called n. 2.2 Exercise
Prove using Taylor’s theorem that
f′′(x) = f(x+h)−2f(x)+f(x−h) +O(h) h2
for 3 times continuously differentiable functions f defined on some interval [x − ε, x + ε].
3 Solving ODEs and Stability
Consider the Ordinary Differential Equation
dy = −20y dt
With initial condition y0 = 1. We can approximate this using the Euler scheme yt+δt = yt − 20ytδt
The true solution is easy to calculate it y(t) = e−20t.
When we have a large number of steps (100), this gives a good approximation. The grey arrows are a visual
representation of the ODE, showing the tangent direction at each point.
We can reduce the step size to 30, and although it is getting rather innacurate, at least it has the correct qualitiative features.
IF we reduce the step size further, we start to see oscillatory behaviour because our estimate heads to the wrong side of the axis. Still, we correctly capture the decay.
When the step size is 10, we see nothing but oscillation, the decay is gone.
When the step size is above 10, the solution starts to blow up exponentially fast. This will happen even if our initial condition is small.
For a general linear ODE of the form
dy =−λy dt
The true solution of the ODE decays exponentially.
The Euler scheme with step size δt gives us yt+δt = yt − λytδt = (1 − λδt)yt
ynδt = y0(1 − λδt)n.
This will grow exponentially if |1 − λδt| > 1, i.e. if δt > λ2 .
This means that if we choose a step size that is too large, small initial conditions will blow up. This is called instability of the method. The Euler scheme is unstable for large values of δt.
Given an n-dimensional ODE of the form dy =Ay
we can attempt to diagonalise A by finding its eigenvectors and eigenvalues.
If it has n independent eigenvectors v1 , . . . vn with eigenvalues λ1 , . . . λn then writing y out with respect
to the basis of eigenvalues, y = y1v1 + . . . + ynvn the ODE becomes dyi = λiyi
If all the eigenvalues are negative, then all the terms will decay. We will then need to choose a time step δt such that |1 − λiδt| < 1 for all the eigenvalues λi so that all the terms in the approximation decay too.
For a high dimensional system there may be a lot of large, negative eigenvalues, so this can force us to choose very low values of δt to solve a ODE using the Euler scheme.
This will be a problem when we solve PDEs because we will be approximating PDEs with high dimensional ODEs. This means using the Euler scheme can be a problem.
3.1 A stable scheme
An alternative approach is to use a different scheme. The Euler scheme can be justified using the forward estimate for the derivative.
TheODEdy =λyisapproximatedby dt
which we rearrange to get the scheme
However, it is possible to start with the backwards estimate
yt−yt−δt =λyt δt
yt+δt−yt =λyt δt
yt+δt = yt + λyt δt.
which we can rearrange to get
yt= 1−λδt yt−δt.
We will call this the backward difference scheme
n ynδt = 1 δt
For an ODE which does not grow exponentially, we have λ ≤ 0 so 1
≤ 1 and so this new scheme will
not grow exponentially either.
3.2 A Central Difference Scheme
Can we use central differences to get a more accurate approach? The most obvious choice of difference equation doesn’t quite work:
yt+δt−yt=λy 1 δt t+2δt
The problem is we have introduced a new unknown yt+ 12 δt. But we can estimate this in terms of yt+δt and yt to get a scheme that does work.
Rearranging
yt+δt − yt = λ(yt + yt+δt) δt 2
yt+δt = 1+ 12λδtyt 1 − 12 λδt
This central difference scheme will be stable so long as λ ≤ 0 and it gives much better results for large n - it has an error of O(δt2) wherease the Euler scheme has an error of order O(δt). With 100 steps, the central scheme and the true solution are visually indistinguishable.
The forward Euler scheme is more intuitive than the backward and central difference schemes, but can be unstable. The backward difference scheme gives the best qualitative results for large step sizes. The central
difference scheme has the best rate of convergence.
For high dimensional problems, one must be particularly careful when using the forward Euler scheme. The more rapid the decay of an eigenvector for the ODE, the smaller the step size required.
4 Exercises 4.1 Exercise
Write a function central_difference_error which computes the size of error of the central difference scheme for solving the ODE
dy = λyt dt
with initial conidition y0 on the interval [0, 1] with n steps. The error is defined to be the distance between the true value at time 1 and the approximation at time 1. Generate a plot that shows the error is of order O(n−2) when λ = −1.
4.2 Exercise
Prove that if λ ≤ 0, the central difference scheme is stable for all values of δt. 5 The Finite Difference Method Part 1
The PDE is
∂V +1σ2S2∂2V +rS∂V −rV =0 ∂t 2 ∂S2 ∂S
We will solve this PDE using finite difference methods. We will approximate the PDE by a large multi- dimensional family of ODE’s. Then we will solve the ODE.
Just as an ODE requires an initial condition, a PDE requires boundary conditions. In this section we will focus on pricing a European call option by PDE methods, so our boundary condition will be that V (T, S) = max{ST −K,0}.
If you consider more general processes SDEs for asset prices than just geometric Brownian motion, you will obtain similar PDEs. So solving the PDE numerically is a good example that generalizes quite easily. We know the analytic formula for the solution of the Black-Scholes PDE to price a call option, so this gives a good test our method works.
5.1 Choice of grid
We plan to solve the Black-Scholes PDE by approximating V (t, S ) on a finite grid of points (ti , Sj ) in a rectangular region [0, T ] × [Smin , Smax ].
Our grid will have N steps in the t direction and M steps in the S direction, giving a total of (N + 1)(M + 1) points.
The result will be a value for V at every point in our grid. We will then be able to generate a 3-d graph of V.
/tmp/ipykernel_43761/3424515936.py:6: MatplotlibDeprecationWarning: Calling
gca() with keyword arguments was deprecated in Matplotlib 3.4. Starting two
minor releases later, gca() will take no keyword arguments. The gca() function
should only be used to get the current axes, or if no axes exist, create new
axes with default keyword arguments. To create a new axes with non-default
arguments, use plt.axes() or plt.subplot().
ax = plt.gca(projection='3d')
5.2 Boundary conditions
As well as the PDE itself, we have one boundary condition. We know V (S, T ) = max{S − K, 0}.
When solving an ODE, the only boundary condition you normally need is the initial condition. For the
you would normally need to be given y0. For the Black-Scholes equation we have a final condition, so we
will be stepping backwards in time when we solve the equation.
If we want to have a well-defined PDE on our finite rectangular grid we will also need to know V (Smax , t) and V (Smin , t). Whatever values we take for V (Smax , t) and V (Smin , t) (within reason) it will be possible to solve the PDE on our grid, so the partial differential equation doesn’t make sense without these boundary conditions.
This seems like a chicken-and-egg problem as we don’t know the V (t, Smax) or V (t, Smin). Fortunately we can estimate them.
dy = f(y) dt
This is possible because if St is much larger than the strike K, the chance that ST will drop below K by maturity will be negligible. A derivivative that pays out ST −K (without any maximum) at time T has value St−e−r(T−t)K byasimplereplicationargument. Youcanreplicatethepayoutbypurchasingoneunitofthe stock and putting e−rT −tK into a risk-free account. So for large S we will have V (t, S) ≈ S − e−r(T −t)K.
To be more precise, we know that the probability of a draw from a normal distribution being more than 8 standard deviations from the mean is negligible.
We know that the change in the log of the stock price from time t to time T is normally distributed with
standard deviation σ we are given that
T − t. We will neglect the drift, which is reasonable over short time horizons. So if
thenthechancethatST
So if we take Smax = K exp(8σT ) we can estimate V (t, Smax) by the replication argument.
From this estimation argument, we will take as our top boundary condition
V (t, Smax) = Smax − e−r(T −t)K
The last argument is heuristic, and rather dodgy.
We haven’t yet got a rigorous theory connecting the probability distribution of the final price and the price of the option.
In fact, the probability distribution of the price depends upon μ, but the PDE doesn’t mention μ so our argument can’t be correct.
When we learn the Feynman-Kac theorem we will see that the argument can be fixed, but you need to choose the pricing measure in the argument, not the physical measure.
For a call option, the bottom boundary condition is simpler. If St = 0 then the SDE for geometric Brownian motionimpliesST =0andsothecalloptionwillbeworthlessattimeT.
We take as bottom boundary condition
V (t, Smin) = 0
To summarize where we have got to, we started with a PDE defined on the [0, T ] × [−∞, ∞]. We have imposed boundary conditions at the top and bottom to give a PDE on a finite region.
WewillapproximateV onagridof(N+1)×(M+1)pointslyingintheregion[0,T]×[Smin,Smax]where Smin = 0
Smax = K exp(8σ
Smax = K exp(8σ
We will take as boundary conditions
V (t, Smin) = 0
V (t, Smax) = Smax − e−r(T −t)K V(T,ST)=max{ST −K,0}
We use this to calculate V at all the other grid points by using a discrete approximation to the PDE ∂V +1σ2S2∂2V +rS∂V −rV =0.
∂t 2 ∂S2 ∂S 6 Finite Difference Methods Part II
Having defined the PDE we wish to solve on a rectangular region, we can now focus on the approximation on a grid.
6.1 Space Discretization
Let us now by discretize in the space (stock price) direction. Let us number the points in our stock price grid {Smin = S0,S1,…,SM = Smax}.
Si := Smin + iδS
δS := Smax − Smin M
Let us write Vti := V (t, Si ) for the corresponding values of V . Using the central difference estimate
∂V V(t,S +δS)−V(t,S −δS) V(t,S (t,Si)≈ i i =
)−V(t,S ) Vi+1 −Vi−1 i−1 = t t
∂S 2δS 2δS 2δS
(t,Si)≈ The Black-Scholes PDE is
) Vi+1 −2Vi +Vi−1 i−1 = t t t
)−2V(t,S )+V(t,S i
∂V +1σ2S2∂2V +rS ∂V −rVi =0
∂t 2 i ∂S2 i ∂S
So we may approximate V (t, S) at the grid points Si by solving the system of (M − 1) ODEs for V i
Vi+1 −2Vi +Vi−1 dt 2 (δS)2
Vi+1 −Vi−1 2δS
for i = 1,2,…,M − 1.
Rearranging the system of ODEs gives
So we have the ODE
i+1 + ciVt
σ 2 S i2 bi := r + (δS)2
ci :=−2(δS)2 −2δS
i−1 i i+1 ≈ aiVt + biVt + ciVt
which holds for i = 1,2,…,M − 1. Together with the boundary conditions
We can write our system of ODEs in matrix form
dV i i−1 i dt ≈ aiVt + biVt
σ 2 S i2 r S i ai := −2(δS)2 + 2δS,
VtM =Smax−e−r(T−t)K Vt0 = 0.
b1 c1 0 0 … 0 0
a 1 V t0 0 0
. 0
a2 b2 c2 0 … 0 0 a3 b3 c3 … 0
. , Wt:=
Λ := . . . 000
So that our ODE can now be written
0 … aM−1 bM−1 dV ≈ΛVt +Wt
0…bc M−2 M−2
… 0 0 V1 a1Vt0
… 0 0V2 0
… 0 0V3 0
V1 b1 c1 0 0
V2 a2 b2 c2 0 dV dV3 0 a3 b3 c3
dt=dt…≈. . . . M−2…. V 0 0 0 0
Let Vt denote the vector containing the components Vti for i = 1, . . . , M − 1. Define a matrix Λ and a
. . . . . + .
. . M − 2 . … bM−2 cM−2 V 0
… aM−1 bM−1 V M−1 cM−1VtM
Following are discussion of solving ODEs, we can find three alternative schemes for solving this ODE. Using backward differences we get the scheme
Vt −Vt−δt ≈ΛVt +Wt. δt
Equivalently, writing 1 for the identity matrix:
Vt−δt ≈ (1 − Λ δt)Vt − δtWt
This is called the explicit method for approximating the solution to the PDE. This is because this is a straight- forward equation describing Vt+δt which you can apply iteratively. Since we are working backwards in time, using backward differences is essentially the same as solving the ODE using the Euler method. So we can expect this method to be unstable unless δt is small.
Each vector Vt corresponds to a column in our grid of points. We start with the vector VT which is known from the final boundary condition.
In the diagram below, the blue dots denote the boundary conditions. We can compute the value at the red dot from three values at time T as shown by the red lines. The set of red lines is called the stencil. It is just a visualisation of which values are used in your finite difference scheme.
Once we’ve computed all the time points T − δt we can move back one time step
Using forward differences we get the scheme
Vt+δt −Vt ≈ΛVt +Wt.
Rearranging
This is called the implicit method for approximating the solution to the PDE. This is because to calculate
Vt+δt ≈(1+Λδt)Vt +δtWt. Vt we will have to solve this system of M − 1 linear equations.
We may write
Vt ≈ (1 + Λ δt)−1(Vt+δt − δtWt)
but this is an inefficient way of calculating Vt. When you solve a system of simultaneous equations you usually use the method of row-reduction, also known as Gaussian elimination. This is much more efficient than computing the inverse and also introduces less rounding errors.
To solve a general linear equation for x given by Ax = b in Python you should use numpy.linalg.solve. If the matrix is mostly zeros, you should use scipy.sparse.linalg.spsolve. This is a partic- ularly simple form of equation since Λ is tridiagonal, so for optimal performance you could use scipy.linalg.solve_banded which is designed to solve such equations efficiently.
Using central differences we get the scheme
Vt+δt −Vt ≈ 12Λ(Vt+δt +Vt)+ 21(Wt+δt +Wt)
Rearranging we get the Crank-Nicolson scheme:
(1+12Λδt)Vt ≈(1−12Λδt)Vt+δt−12δt(Wt+δt+Wt) It too is implicit as we must solve this linear equation at each stage of our iteration.
Again we may write:
Vt ≈(1+12Λ)−1(1−12Λ)Vt+δt−21δt(Wt+δt+Wt) but once again it is better to solve the equations using numpy.linalg.solve
6.2 Summary
To price a call option using finite difference methods we may use either the explicit, implicit or Crank- Nicolson schemes. These schemes may be written as
Vt−δt ≈ (1 − Λ δt)Vt − δtWt
Vt ≈ (1 + Λ δt)−1(Vt+δt − δtWt)
Vt ≈(1+12Λ)−1(1−12Λ)Vt+δt−21δt(Wt+δt+Wt)
When implementing the schemes use an efficient method such as scipy.sparse.linalg.spsolve rather than actually performing the matrix inverse.
The vector W encodes the top and bottom boundary conditions. The matrix Λ is a tridiagonal matrix. Explicit formulae for W and Λ are given above.
To start the iteration we use the final condition
VTi =max{Si−K,0}
which comes from the option payoff, we then step backwards in time to compute V at time 0. 7 Exercises
7.1 Exercise
Suppose that we want to price a European put option with strike K. Choose Smax and Smin = 0 as for a call option. What values should we choose for the top and bottom boundary? Write functions corresponding to your answer.
7.2 Exercise
A digital call option with strike K gives a payoff at maturity T of
1 ST>K 0 otherwise
Using the same values for Smin and Smax, compute appropriate boundary conditions for a digital call option. Write functions corresponding to your answer.
7.3 Exercise
The backwards heat equation for W (t, x) is
Suppose that you wish to solve this equation in a rectangular region and that you have been given top, bottomandfinalboundaryconditions.Discretizethex-coordinateatevenlyspacedpointsxi =xmin+iδx (1 ≤ i ≤ M − 1). Define Wti = W(t,xi). Write W for the vector with components Wti. Following the same method used for the Black-Scholes PDE, compute an ODE approximation to this equation of the form
dW =LW+zt dt
where L is a matrix and zt is a vector which will depend upon the top and bottom boundary conditions. Write a function heat_equation_matrix which returns L.
8 Implementation
We have all the pieces in place, so let us now implement the finite difference method to price a call option using the PDE.
We now wan
程序代写 CS代考 加微信: powcoder QQ: 1823890830 Email: powcoder@163.com