Solving Heston 2 factor PDE in Python (with code)
In this article I will present a method of solving the Heston 2D PDE on the example of European Call and later up-and-out call. I will implement the equation with a finite difference, using Alternating Direction Implicit (ADI). I will also share full code and try to explain it as clearly as possible, as I had problems with understanding it myself. I will blend ADI explanation with code implementation so it is apparent what is going on. If you don’t know how to solve 1D PDE, go and check out my previous posts where I show it with the explicit and Crank-Nicholson methods.
About Heston Model
The Heston model is a widely used stochastic volatility model for option pricing. It extends the Black-Scholes model by allowing volatility to be stochastic, so it varies over time. The Heston model assumes that the asset price price and its variance follow stochastic differential equations (SDEs):
Where:
S_t — underlying asset price at time t
v_t — instantaneous variance (volatility squared)
r — risk‐free interest rate (annualized)
q — continuous dividend yield (if applicable)
κ — mean reversion speed of the variance process
θ — long‐run mean level of variance
σ — volatility of volatility (vol of vtv_tvt)
W_tS, W_tv — two Wiener processes (Brownian motions) with correlation ρ
This model can capture the volatility smile and skew seen in real market data and preserve the forward volatility term structure; however, it has problems with steep skews near maturity. Nevertheless, this is a very useful model, especially for pricing exotic options.
The most efficient method to price European options under the Heston model is using the characteristic function and Fourier inversion methods. The characteristic function of the log asset price is derived using Fourier transform techniques. I showed it in the past blog post.
For path-dependent options like barriers, binaries, and American options, a 2D Partial Differential Equation (PDE) approach is preferred. These options can be priced with 2D PDE with specific boundary conditions, which are fast and easy to implement. I will show it at the end on the example of Up and Out call.
The Heston PDE of option price u(t, S, v) can be written as follows:
For the implementation convenience it can be rearranged and writen in the following form:
Where
That is just a notational shift of all terms to the right‐hand side. In practice, for numerical methods, it is preferred to write it in the form. As you will see later, this representation will come in handy later, when we will implement ADI.
About ADI
ADI stands for alternating direction implicit. As the names suggest, it changes the directions when solving the equation; for example, you solve for one variable explicitly, then for the next implicitly, then again explicitly, and implicitly, and so on. You are effectively cycling between explicit and implicit direction as ADI splits the multi-dimensional operator into separate parts — each corresponding to a different spatial direction — and solves them alternately in an implicit manner. It is an efficient method, and you don’t need many steps for your solution to converge.
With ADI we are going to split our matrix form into 3 parts based on derivatives in respect to S, v and cross-derivatives. It is more efficient to have smaller, easier to solve pieces than inverting a large 2D matrix. Each sub-step involves matrices that are simpler (often tridiagonal or sparse) and much faster to invert.
The ADI method essentially “splits” the multi-dimensional problem into one-dimensional problems. Although this splitting introduces some splitting error, for many financial PDEs (like the Heston), the error is acceptable compared to the huge gains in efficiency.
Implementation
Let’s take care of a few crucial first steps to start working. This comes down to setting up grids and the rest. If you don’t understand this part, go to my previous post about solving Black-Scholes explicitly and with Crank-Nicholson. So let’s set it up:
# Heston model parameters
kappa = 5
theta = 0.04
sigma = 0.15
rho = -0.9
r = 0.02
q = 0.05
K = 100
T = 1.0
NAS = 100 # number of asset steps
NVS = 100 # number of volatility steps
NTS = 200 # number of time steps
Smin = 0.0
Smax = 2 * K
Vmin = 0.0
Vmax = 1
tau_min = 0.0
tau_max = T
ds = (Smax - Smin) / (NAS - 1)
dv = (Vmax - Vmin) / (NVS - 1)
dt = (tau_max - tau_min) / NTS
S = np.linspace(Smin, Smax, NAS)
V = np.linspace(Vmin, Vmax, NVS)
# Here I am creating 2d grid from 1d spot and 1d volatility gird
S2d, V2d = np.meshgrid(S, V, indexing='ij')
N = NAS * NVS
S_flat = S2d.ravel(order='F')
V_flat = V2d.ravel(order='F')Now something, which is different. I am using the ravel function with Fortran order, as the matrix calculations and ADI setup will be possible only with the vector set up this way. I am flattening it in a specific order so that the implementation is correct with matrix notations.
DS, D2S, DV, D2V, DSV = ADI.build_derivatives_2d(NAS, NVS, ds, dv)Here, the operator A1 (containing S-derivative terms) is multiplied by dt and scaled by theta (here, 0.5 as it corresponds to Crank-Nicholson). By subtracting this from the identity matrix, we are effectively setting up a system where the future (unknown) solution is tied to the S derivatives. In the time-stepping loop, this system is solved using an LU factorization (it is shown later). This means that the S-direction terms are handled implicitly — the new solution must be computed by solving a linear system that includes these terms. The same with V-direction.
This matrix represents the contribution of all the terms — including the cross derivative A0 and the parts of A1 and A2 not handled implicitly — from the previous time step. These terms are used directly (explicitly) in updating the solution without solving a new linear system for them.
theta_ADI = 0.5 # 0.5 for Crank Nicholson
A_S_mat = I - theta_ADI * dt * A1
A_V_mat = I - theta_ADI * dt * A2
B_mat = I + (1 - theta_ADI) * dt * (A0 + A1 + A2)
A_S_mat = A_S_mat.tocsc()
A_V_mat = A_V_mat.tocsc()
B_mat = B_mat.tocsc()
A_S_factor = spla.splu(A_S_mat)
A_V_factor = spla.splu(A_V_mat)Now we need to devide the PDE into three parts. In A0 we have cross-derivatives, A1 derivatives in respoect to S and in A2 derivatives in respect to V.
Sdiag = sp.diags(S_flat, 0)
S2diag = sp.diags(S_flat**2, 0)
Vdiag = sp.diags(V_flat, 0)
thetaV = sp.diags(theta - V_flat, 0)
I = sp.eye(N)
A0 = rho * sigma * Sdiag * Vdiag * DSV
A1 = (r - q) * Sdiag * DS + 0.5 * Vdiag * S2diag * D2S - 0.5 * r * I
A2 = kappa * thetaV * DV + 0.5 * sigma**2 * Vdiag * D2V- 0.5 * r * I Now let’s set up the payoff. It is as simple as that. The rest of the boundary conditions will be tackled later when solving, but I will post them here quickly as I discuss the payoff.
- for S=0, the value for the call is 0.
- for S=N, the value becomes Smax — K * exp(-r * T)
- for v = Vmin,
- for v = Vmax,
U_terminal_2d = np.maximum(S2d - K, 0)
U = U_terminal_2d.ravel(order='F')Now let’s look at the ADI implementation with Douglas Scheme. The steps are following:
U_n = U.copy()
for n in range(NTS):
Y0 = A_S_factor.solve(U_n)
Y1 = spla.spsolve(I - 0.5 * dt * A1, Y0 - 0.5 * dt * A1 @ U_n)
Y2 = spla.spsolve(I - 0.5 * dt * A2, Y1 - 0.5 * dt * A2 @ U_n)
U_n = Y2
t_current = T - (n + 1) * dt
for j in range(NVS):
idx_low = 0 + j * NAS
idx_high = (NAS - 1) + j * NAS
U_n[idx_low] = 0.0
U_n[idx_high] = Smax - K * np.exp(-r * (T - t_current))
# At v = 0 - zero slope: U_n[i,0] = U_n[i,1]
for i in range(NAS):
idx_v0 = i + 0 * NAS
idx_v1 = i + 1 * NAS
U_n[idx_v0] = U_n[idx_v1]
# At v = v_max - zero slope: U_n[i,NVS-1] = U_n[i,NVS-2]
for i in range(NAS):
idx_vm = i + (NVS - 1) * NAS
idx_vm_1 = i + (NVS - 2) * NAS
U_n[idx_vm] = U_n[idx_vm_1]
U_final_2d = U_n.reshape((NAS, NVS), order='F')As you can see, all the steps are there, and the boundary conditions are also handled. Now let’s see the results.
Volatility x spot x option values at t=0 (T till maturity)
Volatility x spot x option values at t=T (0 till maturity)
Now I know that it is hard to see, but here is an example of 1D Black-Scholes PDE solved with Crank-Nicholson and 2D Heston PDE fixed at a volatility equal to 20%. I just wanted to show that indeed the stochastic volatility will produce higher prices as there is a volvol parameter, which fattens the tails. Here are the plots. The price discrepancy is most likely too wide; however, they are not calibrated. I just did it as a sanity check.
Now how to price Up and Out calls? You just need to change Smax to the barrier level and change the payoff function. Here is the code and a plot. As you can see, it looks interesting. What’s even more important, it makes sense as the model has stochastic volatility — an extremely important feature in terms of exotic options. As you can see, the higher the volatility, the lower the price of up-and-out, because higher vol makes the option more likely to get kicked out.
B = 120
S_max = B
U_terminal_2d = np.where(S2d < B, np.maximum(S2d - K, 0), 0)
This shows how easy it is to further modify PDE to price more and more things.
Conclusion
As you can see, with ADI we can solve 2D PDEs with superior stability. Play with the code; try to do something better or faster. What’s cool about ADI and solving 2D PDEs is also that you can price many different things with it, e.g., callable bonds or some interest rate models
This post is mostly for educational purposes and the implementation might be proably better as it is my first time solving 2D PDE in Python. I treated this implementation as an exercise and a challenge when I was reading Paul Wilmott’s book, which ultimately directed me to “The Heston Model and Its Extensions”. I found there a practical implementation methods and ultimately based mine on author examples.
If you have any questions, write a comment or reach out to me on LinkedIn; I will try to help.
References:
- “The Heston Model and Its Extensions” — FABRICE DOUGLAS ROUAH
_ “A closed-form solution for options with stochastic volatility with applications to bond and currency options”, Heston
- “ADI Finite Difference Schemes for Option Pricing in the Heston Model with correlation” — KJ IN’T HOUT
Last function:
def build_derivatives_2d(NAS, NVS, ds, dv):
N = NAS * NVS
DS = lil_matrix((N, N))
D2S = lil_matrix((N, N))
DV = lil_matrix((N, N))
D2V = lil_matrix((N, N))
DSV = lil_matrix((N, N))
def idx(i, j):
# Flattening
return i + j * NAS
for j in range(NVS):
for i in range(NAS):
k = idx(i, j)
# First derivative S
if 1 <= i < NAS - 1:
DS[k, idx(i - 1, j)] = -0.5 / ds
DS[k, idx(i + 1, j)] = 0.5 / ds
elif i == 0:
DS[k, idx(i, j)] = -1.0 / ds
DS[k, idx(i + 1, j)] = 1.0 / ds
elif i == NAS - 1:
DS[k, idx(i - 1, j)] = -1.0 / ds
DS[k, idx(i, j)] = 1.0 / ds
# Second derivative S
if 1 <= i < NAS - 1:
D2S[k, idx(i - 1, j)] = 1.0 / ds**2
D2S[k, idx(i, j)] = -2.0 / ds**2
D2S[k, idx(i + 1, j)] = 1.0 / ds**2
elif i == 0:
D2S[k, idx(i, j)] = 1.0 / ds**2
D2S[k, idx(i + 1, j)] = -2.0 / ds**2
D2S[k, idx(i + 2, j)] = 1.0 / ds**2
elif i == NAS - 1:
D2S[k, idx(i - 2, j)] = 1.0 / ds**2
D2S[k, idx(i - 1, j)] = -2.0 / ds**2
D2S[k, idx(i, j)] = 1.0 / ds**2
# First derivative v
if 1 <= j < NVS - 1:
DV[k, idx(i, j - 1)] = -0.5 / dv
DV[k, idx(i, j + 1)] = 0.5 / dv
elif j == 0:
DV[k, idx(i, j)] = -1.0 / dv
DV[k, idx(i, j + 1)] = 1.0 / dv
elif j == NVS - 1:
DV[k, idx(i, j - 1)] = -1.0 / dv
DV[k, idx(i, j)] = 1.0 / dv
# Second derivative
if 1 <= j < NVS - 1:
D2V[k, idx(i, j - 1)] = 1.0 / dv**2
D2V[k, idx(i, j)] = -2.0 / dv**2
D2V[k, idx(i, j + 1)] = 1.0 / dv**2
elif j == 0:
D2V[k, idx(i, j)] = 1.0 / dv**2
D2V[k, idx(i, j + 1)] = -2.0 / dv**2
D2V[k, idx(i, j + 2)] = 1.0 / dv**2
elif j == NVS - 1:
D2V[k, idx(i, j - 2)] = 1.0 / dv**2
D2V[k, idx(i, j - 1)] = -2.0 / dv**2
D2V[k, idx(i, j)] = 1.0 / dv**2
# Cross derivative (central differences)
if (1 <= i < NAS - 1) and (1 <= j < NVS - 1):
DSV[k, idx(i - 1, j - 1)] = 1.0 / (4 * ds * dv)
DSV[k, idx(i + 1, j + 1)] = 1.0 / (4 * ds * dv)
DSV[k, idx(i - 1, j + 1)] = -1.0 / (4 * ds * dv)
DSV[k, idx(i + 1, j - 1)] = -1.0 / (4 * ds * dv)
return DS.tocsr(), D2S.tocsr(), DV.tocsr(), D2V.tocsr(), DSV.tocsr()
