Simulation d'EDP paraboliques par une méthode probabiliste¶
import numpy as np
import numpy.random as npr
import matplotlib.pyplot as plt
I) Equation de la chaleur¶
Nous allons simuler l'équation de la chaleur sur $\mathbb{R}$ suivante
$$ \begin{cases} \partial_t u(t,x) & = \frac{\sigma^2}{2} \Delta u(t,x),\\ u(t=0,x) & = f(x), \end{cases} \qquad (t,x) \in ]0,+\infty[\times \mathbb{R}, $$ où $f$ a de bonnes propriétés pour que l'équation soit bien posée. Dans la suite, on prend $f(x) = e^{-x^2}$ pour simplifier.
Cette equation admet la représentation probabiliste
$$ u(t,x) = \mathbb{E}[f(\sigma W_t + x)] \qquad (t,x)\in ]0,+\infty[\times \mathbb{R}, $$
qu'on peut donc approcher à l'aide d'une méthode de Monte-Carlo
$$ u(t,x) \simeq u_N(t,x) = \frac{1}{N} \sum_{j=1}^N f(W^j_t + x) $$
où $(W^j)_j$ est une famille de mouvement Brownien indépendants. On pourra comparer cette approximation avec la formule explicite
$$ u(t,x) = \int_{-\infty}^{+\infty} g_t(x-y) f(y) dy \qquad (t,x)\in ]0,+\infty[\times \mathbb{R}, $$
où
$$ g_t(x) = \frac{1}{\sqrt{2\pi \sigma^2 t}}e^{-x^2/(2\sigma^2 t)} $$
s'appelle le noyau de la chaleur, qu'on calculera par une méthode de quadrature. Il s'agit d'une convolution de $g$ par $f$.
# définition de la condition initiale
f = lambda x: np.exp(-x**2)
# approximation par Monte-Carlo
N = 4000 # nombre de trajectoire brownienne
sigma = 1. # coefficient de diffusion
h = 0.01 # pas de temps
n = 1000 # nombre de pas de temps, temps total T = nh
x = np.linspace(-10, 10, 100) # grille d'espace
u_N = np.zeros((n+1,len(x))) # vecteur représentant la solution par Monte-Carlo
for _ in range(N): # on répète N fois la construction d'une trajectoire
# simulation de la trajectoire Brownienne
V = npr.normal(0, 1, n)
W = np.sqrt(h)*np.cumsum(V)
W = np.append(0, W) # on ajoute le point de départ 0 au début de chaque trajectoire
WW, xx = np.meshgrid(W, x, indexing='ij')
# nécessaire pour évaluer la solution dans tout l'espace et en tout temps
u_N += f(sigma * WW + xx) / N
# calcul de la formule explicite par quadrature
from scipy.integrate import quad # import de la fonction pour le calcul des intégrales.
# Noyau de la chaleur G(t,x)
def heat_kernel(t, x, sigma):
return np.exp(-x**2 / (2 * sigma**2 * t)) / np.sqrt(2 * np.pi * sigma**2 * t)
# Solution u(t,x) calculée via une convolution de f avec g le noyau de la chaleur
def u(t, x, sigma):
integrand = lambda y: heat_kernel(t, x - y, sigma) * f(y)
result, _ = quad(integrand, -20, 20) # Intégration sur un intervalle suffisamment grand
return result
# évaluation de u(t,x) sur les grilles de temps et d'espace
t = np.array([1, 100, 1000]) * h
x_grid = np.linspace(-10, 10, 25)
u_values = np.array([np.array([u(t, x, sigma) for x in x_grid]) for t in t])
# comparaison de la solution explicite sous forme convolution
# et l'approximation par méthode de Monte-Carlo
cmap = plt.get_cmap("tab10") # définie des couleurs
plt.plot(x, u_N[1,:], color=cmap(0), label="t=0.01, Monte-carlo")
plt.scatter(x_grid, u_values[0,:], color=cmap(0), marker='o', label="t=0.01, convolution")
plt.plot(x, u_N[100,:], color=cmap(1), label="t=1, Monte-carlo")
plt.scatter(x_grid, u_values[1,:], color=cmap(1), marker='o', label="t=1, convolution")
plt.plot(x, u_N[1000,:], color=cmap(2), label="t=10, Monte-carlo")
plt.scatter(x_grid, u_values[2,:], color=cmap(2), marker='o', label="t=10, convolution")
plt.legend()
plt.xlabel(r'$x$', fontsize =15)
plt.show()
On a une bonne adéquation entre la solution calculée par la méthode de Monte-Carlo et la formule explicite où la convolution est évaluée par quadrature.
II) Réduction de variance et théorème de Girsanov¶
On va utiliser le théorème de Girsanov comme une méthode d'échantillonnage préférentiel pour simuler numériquement l'équation de la chaleur. On va s'intéresser à la solution au point $x=5$. On rappelle que pour un échantillon $(X_n)_n$ de $X$, on a
$$ Var(\bar X_n) = \frac{1}{n}Var(X),\qquad \text{avec}\qquad Var(X) = \mathbb{E}[X^2]- \mathbb{E}[X]^2. $$
# définition de la condition initiale centrée en -5
f = lambda x: np.exp(-0.5*(x+5)**2/(2)**2)
# illustration de la condition initiale
x = np.linspace(-20, 10, 400)
plt.plot(x, f(x))
plt.xlabel(r'$x$', fontsize =15)
plt.show()
Commençons par simuler notre EDP avec cette condition initiale en utilisant les méthodes présentées ci-dessus.
# paramètres du problème
sigma = 4. # coefficient de diffusion
h = 0.005 # pas de temps
n = 4000 # nombre de pas de temps, temps total T = nh
x = 5 # position d'observation
t = np.linspace(0, n*h, n+1) # grille des temps d'observation
# approximation par Monte-Carlo
N = 6000 # nombre de trajectoire brownienne
uN = np.zeros(n+1) # vecteur pour l'approximation de u par Monte-Carlo
u2N = np.zeros(n+1) # vecteur pour l'évaluation de la variance
for _ in range(N): # on répète N fois la construction d'une trajectoire
# simulation de la trajectoire Brownienne
V = npr.normal(0, 1, n)
W = np.sqrt(h)*np.cumsum(V)
W = np.append(0, W) # on ajoute le point de départ 0 au début de chaque trajectoire
uN += f(sigma * W + x) / N # calcul de la moyenne pour évaluer E[X]
u2N += f(sigma * W + x)**2 / N # calcul pour évaluer E[X^2]
var = u2N - uN**2 # evaluation Var(X)
from scipy.integrate import quad # import de la fonction pour le calcul des intégrales.
# Noyau de la chaleur g_t(x)
def heat_kernel(t, x, sigma):
return np.exp(-x**2 / (2 * sigma**2 * t)) / np.sqrt(2 * np.pi * sigma**2 * t)
# Solution u(t,x) calculée via une convolution de f avec g le noyau de la chaleur
def u(t, x, sigma):
integrand = lambda y: heat_kernel(t, x - y, sigma) * f(y)
result, _ = quad(integrand, -20, 20) # Intégration sur tout R
return result
# évaluation de u(t,x) sur les grilles de temps et d'espace
u_values = np.array([u(t, x, sigma) for t in t[1:]])
plt.plot(t, uN, label="Monte-Carlo")
plt.plot(t[1:], u_values, lw=3, label="convolution")
plt.xlabel(r'$t$', fontsize =15)
plt.legend(loc="lower right")
plt.show()
Illustrons maintenant le rôle que peut jouer le théorème de Girsanov dans une méthode d'échantillonnage préférentiel.
On considère $X$ la solution de $dX_t = \sigma dW_t$ avec $X_0 = x$, autrement dit $X_t = x + \sigma W_t$ pour tout $t\geq 0$. On peut réécrire cette équation
$$ dX_t = dW_t = \sigma dW^\ast_t + \theta(m-X_t)dt $$
où
$$ W^\ast_t = W_t - \int_0^t \frac{\theta}{\sigma}(m-X_s)ds = W_t- \int_0^t b(s) ds $$
est un mouvement Brownien sous $\mathbb{P}^\ast$, définie par
$$ \frac{d\mathbb{P}^\ast_{\Big| \mathcal{F}_t}}{d\mathbb{P}} = e^{\int_0^t b(s) ds - \frac{1}{2}\int_0^t b^2(s)ds} \qquad \text{avec} \qquad b(s) = \frac{\theta}{\sigma}(m-X_s). $$
On a donc \begin{align*} \mathbb{E}[f(\sigma W_t + x)] & = \mathbb{E}[f(X_t)] \\ & = \mathbb{E}\Big[f(X_t)e^{-\int_0^t b(s)dW_s +\frac{1}{2}\int_0^t b^2(s)ds} e^{\int_0^t b(s)dW_s -\frac{1}{2}\int_0^t b^2(s)ds} \Big] \\ & = \mathbb{E}^\ast\Big[f(X_t)e^{-\int_0^t b(s)dW_s +\frac{1}{2}\int_0^t b^2(s)ds}\Big] \\ & = \mathbb{E}^\ast\Big[f(X_t)e^{-\int_0^t b(s)dW^\ast_s-\frac{1}{2}\int_0^t b^2(s)ds}\Big]. \end{align*}
En posant $$ I_t = e^{-\int_0^t b(s)dW^\ast_s-\frac{1}{2}\int_0^t b^2(s)ds}, $$
on montre par la formule d'Itô que
$$ dI_t = -\frac{\theta}{\sigma}(m-X_t)I_t dW^\ast_t\qquad I_0 = 1. $$
Au final, on peut écrire
$$ \mathbb{E}[f(\sigma W_t + x)] = \mathbb{E}^\ast[f(X_t)I_t], $$
où les processus $X$ et $I$ sont tous deux solutions d'EDS qu'on discrétisera pour les simuler numériquement. On pourra remarquer que sous $\mathbb{P}^\ast$ le processus $X$ est un processus d'Orsntein-Ulhenbeck. On prendra $m$ (le paramètre de moyenne en temps long) comme étant le centre de la condition initiale $f$ (c'est à dire $m=-5$). Cela forcera $X$ à être autour du support de $f$. On pourra jouer avec le paramètre $\theta$ pour voir l'effet apporté.
Regardons maintenant si cette nouvelle représentation de la lsolution $u$ permet de réduire la variance.
def simulation_OU_int(T, n, m, theta, sigma, x0):
"""
T : durée de la simulation.
n : nombre de pas de temps.
theta : vitesse de retour vers m.
m : paramètre vers lequel le processus revient.
sigma : coefficient de diffusion.
x0 : condition initiale du processus.
"""
h = T / n # pas de temps
X = np.zeros(n) # vecteur qui enregistrera les differentes valeurs de X
I = np.zeros(n) # vecteur qui enregistrera les differentes valeurs de I
X[0] = x0 # on enregistre les données initiales
I[0] = 1.
for i in range(0, n-1): # discrétisation des EDS
V = np.sqrt(h) * npr.randn()
X[i+1] = X[i] + theta * (m - X[i]) * h + sigma * V
I[i+1] = (1 - theta * (m - X[i]) * V / sigma) * I[i]
return X, I
# approximation par Monte-Carlo avec un echantillonnage préférentiel
N = 6000 # nombre de trajectoire brownienne
m = -5. # paramètre pour l'EDS de X
theta = 0.1
vN = np.zeros(n) # vecteur pour l'évaluation de la variance
v2N = np.zeros(n) # vecteur pour l'approximation de u par Monte-Carlo
for _ in range(N): # on répète N fois la construction d'une trajectoire
# simulation de la trajectoire de X et I
X, I = simulation_OU_int(n*h, n, m, theta, sigma, x)
vN += f(X) * I / N # calcul pour évaluer E[X]
v2N += (f(X) * I)**2 / N # calcul pour évaluer E[X^2]
# illustration de la solution u(t,x=5)
plt.plot(t, uN, label="Monte-Carlo standard")
plt.plot(t[:-1], vN, label="Monte-Carlo Girsanov" )
plt.plot(t[1:], u_values, lw=3, label="convolution")
plt.xlabel(r'$t$', fontsize =15)
plt.legend(loc="lower right")
plt.show()
# Illustration des variances en fonction du temps
plt.plot(t, u2N, label="Monte-Carlo standard")
plt.plot(t[:-1], v2N, label="Monte-Carlo Girsanov" )
plt.xlabel(r'$t$', fontsize =15)
plt.legend(loc="lower right")
plt.show()
On voit donc que dans le cas d'un échantillonnage préférentiel la variance est environ deux fois plus petite qu'avec une méthode de Monte-Carlo standard.