McKean SDE and Particle method

Keywords: McKean Stochastic Differential Equations, simulation

1. Definitions

A McKean equation for an n-dimensional process X is an SDE in which the drift and volatility depend not only on the current value $X_t$ of the process, but also on the probability distribution $\mathbb{P}_t$ of $X_t$, i.e.

$dX_t = b(t, X_t, \mathbb{P}_t) dt + \sigma(t, X_t, \mathbb{P}_t) dW_t$, $X_0 \in \mathbb{R}^n$

where $W_t$ is a d-dimensional Brownian motion.

One of the most celebrated form of McKean SDE is the McKean-Vlasov SDE, where for $1 \le i \le n$ and $1 \le j \le d$,

$b^{i}(t, x, \mathbb{P}_t) = \int_{\Omega} b^i(t, x, \omega) d\mathbb{P}_t(\omega)$ (scalar)

$\sigma^{i}_j(t, x, \mathbb{P}_t) = \int_{\Omega} \sigma^i_j(t, x, \omega) d\mathbb{P}_t(\omega)$ (scalar)

1.1 The uniqueness and existence of the SDE

If the drift and volatility coefficients are Lipschitz-continuous functions of x and $\mathbb{P}_t$ wrt Wasserstein distance, the solution of the SDE has a strong unique solution.

Theorem Let $b : \mathbb{R}_{+} \times \mathbb{R}^n \times \mathcal{P}_2(\mathbb{R}^n) \xrightarrow{} \mathbb{R}^n$ and $\sigma: \mathbb{R}_{+} \times \mathbb{R}^n \times \mathcal{P}_2(\mathbb{R}^n) \xrightarrow{} \mathbb{R}^{n \times d}$ be Lipschitz-continuous functions and satisfy a linear growth condition

$|b(t, X, \mathbb{P}) - b(t, Y, \mathbb{Q})| + |\sigma(t, X, \mathbb{P}) - \sigma(t, Y, \mathbb{Q})| \le C(|X - Y| + d(\mathbb{P}, \mathbb{Q}))$

and

$|b(t, X, \mathbb{P})| + |\sigma(t, X, \mathbb{P})| \le C(1 + |X|)$ (C is a positive constant)

$\mathcal{P}_2(\mathbb{R}^n)$ denotes the probability measures with finite second order moment (notice the difference between the $L^2$-space). The Wasserstein distance is defined as:

$d(\mu, v) = \inf_{\tau \in \mathcal{P}(\mathbb{R}^n \times \mathbb{R}^n)} \Big(\int_{\mathbb{R}^n \times \mathbb{R}^n} |x - y| \tau(dx, dy)\Big)^{\frac{1}{2}}$

Note that the marginal distribution of $\mathcal{P}(\mathbb{R}^n \times \mathbb{R}^n)$ should be $\mu$ and $v$. The optimizer here is the possible coupling method between two random variables. Then the nonlinear SDE

$dX_t = b(t, X_t, \mathbb{P}_t) dt + \sigma(t, X_t, \mathbb{P}_t)dW_t$,

where $\mathbb{P}_{t}$ denotes the probability distribution of $X_t$, admits a unique solution such that $\mathbb{E}(\sup_{0 \le t \le T}|X_t|^2) < \infty$

Proof CRC-Nonlinear Option Pricing 251

The probability density function $p(t, y) dy = \mathbb{P}_t(dy)$ of $X_t$ is a solution to the Fokker-Planck PDE:

$-\partial_t p(t, x) - \sum_{i = 1}^n \partial_i (b^{i}(t, x, \mathbb{P}_t)p(t, x)) + \frac{1}{2} \sum_{i, j = 1}^n \partial_{i, j} \Big(\sum_{k = 1}^d \sigma^{i}_{k}(t, x, \mathbb{P}_t) \sigma_k^j(t, x, \mathbb{P}_t)p(t, x)\Big) = 0$

with initial condition

$\lim\limits_{t \xrightarrow{} 0} p(t, x) = \delta(x - X_0)$ (Dirac function)

(Note that it is a non-linear PDE since the drift and diffusion depends on the $p(x, t)$)

2. Particle method (McKean-Vlasov)

We can approximate the $\mathbb{P}_t^N$ with the empirical distribution $\frac{1}{N} \sum_{i = 1}^N \delta_{X_t^{i, N}}$, where the $(X^{i, N})_{1 \le i \le N}$ are solutions to the $(\mathbb{R}^n)^N$-dimensional classical (linear) SDE, i.e.

$dX_t^{i, N} = \Big(\int b(t, X_t^{i, N}, y) d\mathbb{P}_t^N(y)\Big)dt + \Big(\int \sigma(t, X_t^{i, N}, y) d\mathbb{P}_t^N(y)\Big)dW^i_t$, which is equivalent to

$dX_t^{i, N} = \sum_{j = 1}^N b(t, X_t^{i, N}, X_t^{j, N}) dt + \sum_{j = 1}^N \sigma^i(t, X_t^{i, N}, X_t^{j, N}) dW^i_t$ (*)

(Law $\Big(X_0^{i, N}\Big)$)

The particle method differs from classical Monte Carlo methods as it involves a system of N interacting particles.

2.1 Propagation of chaos and convergence of the particle method

Let’s consider scalar case first. Let $\mu_t^{N}$ the density of $(X_t^{1, N}, \dots, X_t^{N, N})$. We have the marginal laws that:

$\mu_t^{k}(x_1, \dots, x_k) = \int \mu_t^{(N)}(x_1, \dots, x_N)d x_{k + 1} \dots d x_{N}$

The Fokker-Planck PDE of (*) can be written as:

$\partial_t \mu_t^{N}(x_1, \dots, x_N) = -\frac{1}{N} \sum_{i, j = 1}^N \partial_{x_i}{b(x_i, x_j) \mu_t^{(N)}} + \frac{1}{2N^2}\sum_{i, p, q = 1}^N \partial^2_{x_i} {\sigma(x_i, x_p)\sigma(x_i, x_q)\mu_t^{(N)}}$ ($L^2$ adjoint see )