Pricing European Calls by Evolving the Forward Density

# European option priced with forward Kolmogorov and Finite Difference method


# European option priced with forward Kolmogorov and Finite Difference method












I explored the Kolmogorov equations, which describe the evolution of probability distributions in stochastic processes, and was thinking how they could be used to price European options. However, turning this idea into a practical solution was more complicated than expected.

Let me describe the approach I took to price an European Call option. I have added the code at the end.

Step 1 - Be $X$ the underlying asset of European option described by the SDE $dX_t=\mu X_tdt+\sigma X_tdW_t$. The Forward Kolmogorov equation for the probability density function $p(x, t)$ of the random variable $X_t$ is $$\text{(FKE)}\hspace{0.5cm}\frac{\partial p(x,t)}{\partial t} = -\frac{\mu(x,t)p(x, t)}{\partial x} + \frac{\partial^2 \sigma^2(x,t)p(x,t)}{\partial x^2}.$$

Step 2 - To solve this parabolic PDE numerically, define a grid and the following finite difference approximations for the derivatives:

- Set $p_i^j = p(x_i, t_j)$ with $i \in \{1, ..., N\}$, $j \in \{1, ..., M\}$.

- $\frac{\partial p(x,t)}{\partial t} = \frac{p_i^{j+1}-p_i^j}{\Delta t}$

- $\frac{\partial p(x,t)}{\partial x} = \frac{p_{i+1}^{j}-p_{i-1}^j}{2\Delta x}$

- $\frac{\partial^2 p(x,t)}{\partial x^2} = \frac{p{_i+1}^{j} - 2 p_i^j - p_{i-1}^j}{(\Delta x)^2}$

Step 3 - Next, define the time-stepping scheme, which needs to be solved for each time step $t=\{0,...T\}$, by inserting our finite difference approximations.

(FKE) $\Rightarrow \frac{p_{i}^{j+1}- p_{i}^{j}}{\Delta t} = -\mu(x_i) \frac{p^{i+1}_{j}-p^{i-1}_{j}}{2\Delta x} + \frac{1}{2}\sigma^2(x_i)\frac{p^{i+1}_{j}-2p^{i}_{j}+p^{i-1}_{j}}{(\Delta x)^2}$ $ \Leftrightarrow p^{i}_{j+1} = p^{i}_{j} +\Delta t (-\mu(x_i) \frac{p^{i+1}_{j}-p^{i-1}_{j}}{2\Delta x} + \frac{1}{2}\sigma^2(x_i)\frac{p^{i+1}_{j}-2p^{i}_{j}+p^{i-1}_{j}}{(\Delta x)^2})$

With initial condition starting values to be zero everywhere except near $x_0$, and whose integral over the entire line is equal to one. $$p^{0}_{i} = 0 \hspace{0.5cm}\forall i\in\{ 1, ..., N+1\} \hspace{0.5cm} \text{and} \hspace{0.5cm} p_{\lfloor x_0 /\text{dx} \rfloor}^{0} = 1 / \text{dx}, \hspace{0.3cm}\text{dx}=\frac{x_{\text{max}}-x_{\text{min}}}{N}$$ The boundary condition is $$p^{j}_{0} = p^{j}_{N+1} = 0, \hspace{0.5cm} \forall j\in \{1, ..., M\}. $$

Step 4 - Set parameters and run the program to solve for the probability density function $p(x, T)$. $$x=[0, 300], \hspace{0.4cm} N=100, \hspace{0.4cm}T=1, \hspace{0.4cm}M=1000, \hspace{0.4cm}r=0.05, \hspace{0.4cm}\mu=r, \hspace{0.4cm}\sigma=0.2, \hspace{0.4cm}x_0=100, \hspace{0.4cm}K=100. $$ For each timestep $j$ calculate the vector p$^j = [p_1^j, ..., p_{N+1}^j]^\intercal$ until one arrives at time $T$ with the probability density vector $p(x, T)$.

The graphs below show the evolution of the probability density function starting at time $t=0$ with the initial underlying value $X_0=x_0$ and the target probability function $p(x, T)$, which is used in the next step for pricing.

Step 5 - Combine the probability density function $p(x,T)$ and payout function $\max (x-K, 0)$ to receive the distribution of payoff and ultimately the expected value of the payoff under the risk-neutral measure $\mathbb{E}^{Q}[\max(X_T - K, 0)]$.

Doing a simpel calculation of the expected value $\mathbb{E}^\mathbb{Q} = \sum_{i=1}^{N+1} p(x_i, T) * \max(x_i - K,0) * dx$ returns an option value of $16.9786$

The option value for the BS model is $10.4506$.

My approach seems logical, but I’m unsure if I am missing something critical here. What should I focus on to refine this numerical approach and align it with the Black-Scholes pricing model?

fke_call_price.py

```
import numpy as np
import matplotlib.pyplot as plt
from scipy.stats import norm

#%% Functions
# Numeric Solver for Fokker-Planck-Gleichung of Ito process
def fokker_planck_solver(mu, sigma, S0, x_min, x_max, T, N, M):
    dx = (x_max - x_min) / N
    dt = T / M
    x = np.linspace(x_min, x_max, N+1)
    p = np.zeros((N+1, M+1))
    
    # initial condition
    p[:, 0] = np.zeros(N+1)
    idx = int(S0 / dx)
    p[idx, 0] = 1 / dx
    
    for j in range(M):
        for i in range(1, N):
            drift_term = -mu(x[i])*(p[i+1, j] - p[i-1, j]) / (2*dx)
            diffusion_term = sigma(x[i])**2*(p[i+1, j] - 2*p[i, j] + p[i-1, j]) / (dx**2)
            p[i, j+1] = p[i, j] + dt * (drift_term + .5 * diffusion_term)
        
        # boundary conditions (absorbing)
        p[0, j+1] = 0
        p[-1, j+1] = 0
    
    return x, p
#%% Main
# Parameter
x_min, x_max = 0, 300
N = 100
T = 1
M = 1000

r=0.05
mu = lambda x: r * x
sigma = lambda x: 0.2 * x

S0=100
K=100

# Calculation of FKE & BS price
dx = (x_max - x_min) / N
x, p = fokker_planck_solver(mu, sigma, S0, x_min, x_max, T, N, M)
FKE_price = np.sum(np.maximum(x - K, 0) * p[:, -1] * dx) * np.exp(-r * T)
print(f"FKE Preis: {FKE_price:.4f}")

d1 = (np.log(S0/K) + (r + 0.5*sigma(1)**2)*T) / (sigma(1)*np.sqrt(T))
d2 = d1 - sigma(1)*np.sqrt(T)
BS_price = S0*norm.cdf(d1) - K*np.exp(-r*T)*norm.cdf(d2)
print(f"Black-Scholes Preis: {BS_price:.4f}")
```