Erreur forte et erreur faible
import math
import numpy as np
import torch
dev_cpu = torch.device("cpu")
rng_cpu = torch.Generator(device=dev_cpu)
dev_mps = torch.device("mps")
rng_mps = torch.Generator(device=dev_mps)
from scipy import stats
import pandas as pd
import matplotlib.pyplot as plt
import seaborn as sns
sns.set_theme()On commence par une première section sur la simulation exacte du processus d’Ornstein-Uhlenbeck. Il s’agit d’illustrer que le choix des instructions peut avoir un impact sur la performance des simulations et d’illustrer la différentiation automatique du module pytorch.
On considère \((X_t)_{t \in [0,T]}\) processus à valeurs dans \(\mathbf{R}\) solution de l’EDS \[ \operatorname{d}\! X_t = \lambda (\mu - X_t) \operatorname{d}\! t + \sigma \operatorname{d}\! W_t, \quad X_0 = x \in \mathbf{R}, \] où \(\lambda > 0\), \(\mu, \sigma \in \mathbf{R}\).
On pose \(Y_t = e^{\lambda t} (X_t - \mu)\) qui vérifie \[ \begin{aligned} \operatorname{d}\! Y_t &= \lambda e^{\lambda t} (X_t - \mu) \operatorname{d}\! t + e^{\lambda t} \operatorname{d}\! X_t, \\ &= \sigma e^{\lambda t} \operatorname{d}\! W_t, \end{aligned} \] et \(Y_0 = x-\mu\).
Le processus \((Y_t)_{t \in [0,T]}\) est une intégrale de Wiener, donc un processus de Markov et un processus Gaussien de fonction moyenne \(m_Y(t) = x-\mu\) et de fonction de covariance \[ \begin{aligned} K_Y(s, t) &= \operatorname{cov}(Y_s, Y_t) = \sigma^2 \mathbf{E}\bigg[\int_0^s e^{\lambda u} \operatorname{d}\! W_u \int_0^t e^{\lambda u} \operatorname{d}\! W_u \bigg], \\ &= \frac{\sigma^2}{2 \lambda} \big( e^{2 \lambda \min(s,t)} - 1 \big). \end{aligned} \] En particulier, \[ \forall 0 \le s \le t, \quad \mathcal{L}\big(Y_t| \sigma(Y_u, u \le s)\big) \sim \mathcal{N}\Bigl( Y_s ; \frac{\sigma^2}{2 \lambda} \bigl( e^{2 \lambda t} - e^{2 \lambda s}\bigr) \Bigr) \]
On en déduit que \((X_t)_{t \in [0,T]}\) est aussi un processus Gaussien de fonction moyenne \(m_X\) et de fonction de covariance \(K_X\) définies par \[ m_X(t) = \mu + (x-\mu) e^{-\lambda t} \quad \text{et} \quad K_X(s, t) = \frac{\sigma^2}{2 \lambda} e^{-\lambda (s+t)} \big( e^{2 \lambda \min(s,t)} - 1 \big). \] En particulier, \[ \forall 0 \le s \le t, \quad \mathcal{L}\big(X_t| \sigma(X_u, u \le s)\big) \simeq \mathcal{N}\Bigl( \mu + (X_s-\mu) e^{-\lambda(t-s)} ; \frac{\sigma^2}{2 \lambda} \bigl( 1 - e^{-2\lambda(t-s)} \bigr) \Bigr) \]
On illustre dans ce notebook un peu de programmation orientée objet (POO). En Python on peut facilement écrire un code sans avoir recourt à de la POO (définition de classe, etc.) mais dans certains cas la syntaxe POO est agréable. Voici un exemple de classe ornstein_uhlenbeck dont on définira les méthodes direct_simulation et recursive_simulation dans la suite.
class ornstein_uhlenbeck:
def __init__(self, lambd, mu, sigma):
self.lambd = lambd
self.mu = mu
self.sigma = sigma
def __repr__(self):
return f"Ornstein-Uhlenbeck de paramètres: " + \
f"lambda = {self.lambd}, mu = {self.mu}, sigma = {self.sigma}"
# déclaration pas nécessaire mais habitude du C++
def direct_simulation(self, x0, final_time, gaussians):
pass
def recursive_simulation(self, x0, final_time, gaussians):
passDans la suite on se donne les paramètres suivants: \(\lambda = 5\), \(\mu=1\), \(\sigma =0.3\), \(X_0 = 3\) et \(T = 1\) de sorte que la loi de \(X_T\) est \(\mathcal{N}\big( \mu + (x-\mu) e^{-\lambda T}; \frac{\sigma^2}{2\lambda} (1-e^{-2 \lambda T}) \big)\).
ou = ornstein_uhlenbeck(lambd=5, mu=1, sigma=0.3)
print(ou)
x0, T = 3.0, 1
print(f"Point initial x0 = {x0} et temps terminal T = {T}")
mean_T = ou.mu + (x0-ou.mu) * math.exp(-ou.lambd*T)
var_T = ou.sigma**2 / (2*ou.lambd) * (1 - math.exp(-2*ou.lambd*T))
print(f"mean: {mean_T} \t var: {var_T}") Ornstein-Uhlenbeck de paramètres: lambda = 5, mu = 1, sigma = 0.3
Point initial x0 = 3.0 et temps terminal T = 1
mean: 1.013475893998171 var: 0.008999591400632136On compare ici deux codes possibles pour simuler des trajectoires du processus d’Ornstein-Uhlenbeck:
On définit la méthode direct_simulation de la classe ornstein_uhlenbeck qui permet la construction d’un ensemble de trajectoires à partir d’un tenseur gaussians qui contient \(M \times N\) gaussiennes centrées réduites (\(M\) le nombe de trajectoires et \(N\) qui règle le pas de discrétisation \(h = T/N\)).
L’implémentation est la suivante: on simule les accroissements de \((Y_t)_{t \in [0,T]}\) aux instants \(t_k = k h\), puis on reconstruit les \(Y_{t_k}\) (via cumsum) et on obtient \(X_{t_k}\) à partir de \(Y_{t_k}\).
def ou_direct_simulation(self, x0, final_time, gaussians):
n_steps, n_paths = gaussians.shape
dt = final_time / n_steps
timesteps = (torch.arange(n_steps+1, device=gaussians.device)*dt)[:,None]
dY = torch.empty((n_steps+1, n_paths), device = gaussians.device)
dY[0] = x0 - self.mu
dY[1:] = torch.sqrt(self.sigma**2 / (2*self.lambd) *\
(torch.exp(2*self.lambd*timesteps[1:]) -\
torch.exp(2*self.lambd*timesteps[:-1]))) * gaussians
sample = self.mu+torch.exp(-self.lambd*timesteps)*torch.cumsum(dY,axis=0)
return sample
ornstein_uhlenbeck.direct_simulation = ou_direct_simulationOn définit la méthode recursive_simulation de la classe ornstein_uhlenbeck qui implémente la simulation des trajectoires en utilisant une boucle for pour construire itérativement \(X_{t_{k+1}}\) à partir de \(X_{t_k}\) (les \(M\) trajectoires en même temps).
def ou_recursive_simulation(self, x0, final_time, gaussians):
n_steps, n_paths = gaussians.shape
dt = final_time / n_steps
kap_ = math.exp(-self.lambd * dt)
sig_ = math.sqrt(self.sigma**2/(2*self.lambd)*(1 - kap_**2))
sample = torch.empty((n_steps+1, n_paths), device = gaussians.device)
sample[0] = x0
for k in range(1, n_steps+1):
sample[k] = self.mu + (sample[k-1] - self.mu) * kap_ +\
sig_*gaussians[k-1]
return sample
ornstein_uhlenbeck.recursive_simulation = ou_recursive_simulationLa boucle for sur les itérations en temps est acceptable, par contre en aucun cas il faut faire une boucle sur les trajectoires (avec les modules numpy ou pytorch).
n_paths = 100000
n_steps = 100
devs = { "cpu": { "device": dev_cpu, "generator":rng_cpu },
"mps": { "device": dev_mps, "generator":rng_mps }}
algos = { "direct": ou.direct_simulation,
"recursive": ou.recursive_simulation }
benchmarks = { "cpu": dict() , "mps": dict() }
def sampling(algo_function: callable, **kwargs):
gaussians = torch.randn((n_steps, n_paths), **kwargs)
sample = algo_function(x0, T, gaussians)
return sample
from itertools import product
for d, a in product(devs, algos):
time = %timeit -o sampling(algos[a], **devs[d])
benchmarks[d][a] = f"{time.average * 1000:.5} ms"
pd.DataFrame(benchmarks)65 ms ± 1.42 ms per loop (mean ± std. dev. of 7 runs, 10 loops each)
71.3 ms ± 1.28 ms per loop (mean ± std. dev. of 7 runs, 10 loops each)
4.97 ms ± 21.8 μs per loop (mean ± std. dev. of 7 runs, 1,000 loops each)
5.46 ms ± 203 μs per loop (mean ± std. dev. of 7 runs, 100 loops each)| cpu | mps | |
|---|---|---|
| direct | 65.008 ms | 4.9743 ms |
| recursive | 71.25 ms | 5.4635 ms |
Exercice: reprendre ces fonctions en inversant les axes. Que constatez-vous ?
n_steps, n_paths = 100, 200
timesteps = np.arange(n_steps+1) * (T/n_steps)
gaussians = torch.randn((n_steps, n_paths),device=dev_cpu,generator=rng_cpu)
paths = ou.direct_simulation(x0, T, gaussians)
fig, ax = plt.subplots()
for path in paths.T:
ax.plot(timesteps, path, color='C0', alpha=0.3)
ax.plot(timesteps, paths[:,0], color='C1', lw=2, label='One path')
ax.legend()
ax.set_title(f"{n_paths} OU paths with {n_steps} discretisation steps")
plt.show()
autogradLe module autograd de pytorch est un outil puissant qui implémente la différentiation automatique dite backward ou adjointe (AAD). C’est une méthode efficace pour calculer le gradient d’une fonction de \(\mathbf{R}^m \to \mathbf{R}\) ou la matrice Jacobienne d’une fonction de \(\mathbf{R}^m \to \mathbf{R}^n\).
Contrairement aux méthodes numériques classiques, comme les différences finies ou les méthodes des chocs,la différentiation automatique repose sur une approche analytique basée sur les règles de dérivation des fonctions composées.
PyTorch construit un graphe computationnel dynamique (coût assez cher) lors de l’exécution du code : chaque opération est enregistrée avec ses dépendances, ce qui permet de calculer les dérivées exactes de manière rétroactive à partir des résultats finaux.
Pour en savoir plus, consulter
Voici un exemple tiré du blog précédent.
x = torch.tensor([1.0, 2.0], requires_grad=True)
y = torch.empty(3)
y[0] = x[0] ** 2
y[1] = x[0] ** 2 + 5 * x[1] ** 2
y[2] = 3 * x[1]
v = torch.tensor([1.0, 2.0, 3.0])
y.backward(v) # Vector-Jacobian Product
print("y:", y)
print("x.grad:", x.grad)
## Manual computation.
dydx = torch.tensor(
[
[2 * x[0], 0],
[2 * x[0], 10 * x[1]],
[0, 3],
]
)
explicit_grad = v @ dydx
print("equality:", torch.equal(x.grad, explicit_grad))y: tensor([ 1., 21., 6.], grad_fn=<CopySlices>)
x.grad: tensor([ 6., 49.])
equality: TrueSoit \(f: \mathbf{R} \to \mathbf{R}\) la fonction définie par \[ f(x) = \mathbf{E} \bigl[ X^2_T | X_0 = x \bigr]. \] Par le calcul précédent on sait que la fonction \(f\) est quadratique et on connait sa dérivée \(f'\). Plus précisément on a \[ f(x) = \big((x-\mu) e^{-\lambda T} + \mu \big)^2 + \frac{\sigma^2}{2\lambda}(1 - e^{-2 \lambda T}) \] et donc \(f'(x) = 2 x e^{-\lambda T}\).
f_x0 = var_T + mean_T**2
df_x0 = 2 * mean_T * math.exp(-ou.lambd*T)
print(f"f({x0}) = {f_x0:.5} \t f'({x0}) = {df_x0:.5}")f(3.0) = 1.0361 f'(3.0) = 0.013657Le but est de faire un code de différentiation automatique qui calcule automatiquement la sensibilité de \(f\) par rapport à la variable \(x\) (sans utiliser la formule explicite…)
n_steps, n_paths = 100, 100000
## on crée un tensor x0_ avec l'attribut `requires_grad=True`
x0_ = torch.tensor(float(x0), device=dev_cpu, requires_grad=True)
## on calcule une quantité qui dépend de x0: évaluation forward
## dans cette phase il y a création du graphe de calcul
gaussians = torch.randn((n_steps, n_paths),device=dev_cpu,generator=rng_cpu)
sample = ou.direct_simulation(x0_, T, gaussians)
f_x0 = (sample[-1]**2).mean()
## évaluation rétrograde en utilisant le graphe de calcul
## pour évaluer la dérivée qui sera stockée dans le champ .grad de x0_
f_x0.backward()
df_x0 = x0_.grad
print(f"Approximation AAD: f({x0}) = {f_x0.item():.5} \t" + \
f"f'({x0}) = {df_x0.item():.5}")Approximation AAD: f(3.0) = 1.0367 f'(3.0) = 0.013661Comment obtenir un intervalle de confiance pour f'(3) ? En fait en échangeant dérivation et espérance et en notant \(X_T^x\) la solution issue de \(x\) en zéro (\(X_0 = x\)) on a aussi la représentation \[
f(x) = \mathbf{E} \Big[ (X^x_T)^2 \Big] \quad \text{et} \quad
f'(x) = \mathbf{E} \Big[ \frac{\partial}{\partial x} (X^x_T)^2 \Big].
\] Ainsi on a envie de dériver (mettre en oeuvre l’AAD) sous le signe somme c’est à dire sur les réalisations de \(((X^{(j), x}_T)^2)_{1 \le j \le M}\) utilisées dans l’estimateur Monte Carlo \(\frac{1}{M} \sum_{j=1}^M (X^{(j), x}_T)^2\).
Etant donné un échantillon i.i.d. de variables aléatoires réelles \((Z^{(j),x})_{1 \le j \le M}\) qui dépendent d’un paramètre \(x \in \mathbf{R}\) on veut calculer \(\partial_x Z^{(j),x}\). Pour utiliser le mécanisme VJP (vector jacobian product) de pytorch on utilise la représentation \[ \big(\partial_x Z^{(1),x}, \dots, \partial_x Z^{(M),x} \big) = \big( 1, \dots, 1 \big) J_{x}(g) \] où \(J_{x}(g)\) est la matrice Jacobienne de \(g:\mathbf{R}^M \to \mathbf{R}^M\) qui au point \((x,\dots,x) \in \mathbf{R}^M\) donne la valeur \((Z^{(1),x}, \dots, Z^{(M),x})\) (les dérivées croisées sont nulles).
Le code pytorch correspondant s’écrit alors:
n_steps, n_paths = 100, 10000
x0_ = torch.full((n_paths,), float(x0), device=dev_cpu, requires_grad=True)
gaussians = torch.randn((n_steps, n_paths),device=dev_cpu,generator=rng_cpu)
sample = ou.direct_simulation(x0_, T, gaussians)
payoffs = sample[-1]**2
v = torch.ones_like(x0_)
payoffs.backward(v)
grad_payoffs = x0_.grad
print(grad_payoffs)tensor([0.0127, 0.0132, 0.0133, ..., 0.0134, 0.0127, 0.0134])monte_carlodef monte_carlo(sample: np.ndarray, proba: float = 0.95) -> dict:
mean = np.mean(sample)
var = np.var(sample, ddof=1)
alpha = 1 - proba
quantile = stats.norm.ppf(1 - alpha / 2)
ci_size = quantile * np.sqrt(var / sample.size)
return {
"mean": mean,
"var": var,
"lower": mean - ci_size,
"upper": mean + ci_size
}Ainsi on peut obtenir un intervalle de confiance sur la sensibilité \(f'(3)\).
pd.DataFrame({
"f(3)": monte_carlo(payoffs.cpu().detach().numpy()),
"f'(3)": monte_carlo(grad_payoffs.cpu().detach().numpy())
})| f(3) | f'(3) | |
|---|---|---|
| mean | 1.038342 | 0.013673 |
| var | 0.036951 | 0.000002 |
| lower | 1.034575 | 0.013648 |
| upper | 1.042110 | 0.013698 |
Soit \((x_t)_{t \in [0,T]}\) solution d’une EDS \[ \operatorname{d}\! x_t = b(x_t) \operatorname{d}\! t + \sigma(x_t) \operatorname{d}\! W_t, \quad x_0 \in \mathbf{R^d}. \] On se donne des instants de discrétisation \(0=t_0 < t_1 < \dots < t_N = T\) et dans la suite on spécifie la grille homogène \(t_n = n \frac{T}{N} = n h\) avec \(h = \frac{T}{N}\) appelé le pas de discrétisation.
Schéma d’Euler
On rappelle le schéma d’Euler pour une diffusion: \[
X_{t_{n+1}} = X_{t_n} + b(X_{t_n}) h + \sigma(X_{t_n}) \sqrt{h} G_{n+1}, \quad X_0 = x_0.
\] où dans toute la suite \((G_1,\dots,G_N)\) est un vecteur Gaussien centré réduit (composantes indépendantes). La loi de \(\sqrt{h} G_{n+1}\) est celle de l’accroissement \(W_{t_{n+1}} - W_{t_n}\).
Schéma de Milstein
Avec les mêmes notations que précédemment on définit le schéma de Milstein pour une diffusion (en dimension 1): \[
X_{t_{n+1}} = X_{t_n} + b(X_{t_n}) h + \sigma(X_{t_n}) \sqrt{h} G_{n+1} + \frac{1}{2} (\sigma \sigma')(X_{t_n})(G_{n+1}^2 - 1) h.
\]
On peut aussi utiliser la version équivalente de Newton, schéma avec des propriétés théoriques équivalentes mais sans la dépendance en la dérivée \(\sigma'\) \[ X_{t_{n+1}} = X_{t_n} + b(X_{t_n}) h - \sigma(X_{t_n}) \sqrt{h} + \tilde \sigma_{n+1} \sqrt{h} (G_{n+1} + 1), \] avec \(\tilde \sigma_{n+1} = \sigma\left(X_{t_n} + \frac{1}{2} \sigma(X_{t_n}) \sqrt{h} \big(G_{n+1} - 1\big) \right)\).
L’équation différentielle stochastique (EDS) est donnée par \[ \operatorname{d}\! x_t = r x_t \operatorname{d}\! t + \sigma x_t \operatorname{d}\! W_t, \quad x_0 > 0 \]
Pour Black-Scholes on connaît la vraie solution de l’EDS et il n’est donc pas nécessaire de faire un schéma de discrétisation. On utilise cet exemple simplement pour illustrer numériquement le comportement des schémas numériques. On rappelle la solution \[ \forall t \in [0,T], \quad x_t = x_0 e^{\left(r - \frac{\sigma^2}{2} \right) t + \sigma W_t}. \]
class BlackScholes:
def __init__(self, r, sigma):
self.r = r
self.sigma = sigma
def exact_simulation(self, x0, final_time, dW):
n_steps, n_paths = dW.shape
dt = final_time / n_steps
sample = torch.empty((n_steps+1,n_paths), device=dW.device)
sample[0] = x0
for n in range(1, n_steps+1):
sample[n] = sample[n-1] * torch.exp(
(self.r - 0.5*self.sigma**2)*dt + self.sigma*dW[n-1])
# on pourrait utiliser np.cumprod à la place de cette boucle...
return sample
def euler_simulation(self, x0, final_time, dW):
n_steps, n_paths = dW.shape
dt = final_time / n_steps
sample = torch.empty((n_steps+1,n_paths), device=dW.device)
sample[0] = x0
for n in range(1, n_steps+1):
sample[n] = sample[n-1] + self.r*sample[n-1]*dt +\
self.sigma*sample[n-1]*dW[n-1]
return sample
def milstein_simulation(self, x0, final_time, dW):
n_steps, n_paths = dW.shape
dt = final_time / n_steps
sample = torch.empty((n_steps+1,n_paths), device=dW.device)
sample[0] = x0
for n in range(1, n_steps+1):
sample[n] = sample[n-1] + self.r*sample[n-1]*dt +\
self.sigma*sample[n-1]*dW[n-1] +\
0.5*self.sigma**2*sample[n-1]*(dW[n-1]**2-dt)
return sampleT = 5
X = BlackScholes(r=0.04, sigma=0.1)
n_steps, n_paths = 10, 30
dW = math.sqrt(T/n_steps) * torch.randn((n_steps, n_paths))tt = np.linspace(0, T, n_steps+1)
fig, ax = plt.subplots(layout='tight')
ax.plot(tt, X.exact_simulation(x0, T, dW), color='C0')
ax.plot(tt, X.euler_simulation(x0, T, dW), color='C0',linestyle='dotted')
ax.plot(tt, X.milstein_simulation(x0, T, dW), color='C0',linestyle='dashed')
fig.suptitle(fr'M={n_paths} trajectoires sur [0,{T}]' +\
fr'discrétisées avec N={n_steps} points')
plt.show()
On dit que le schéma est d’erreur forte \(\mathbf{L}^p\), \(p > 1\) d’ordre \(\alpha\) si \[ \left\lVert \sup_{t \in [0,T]} |x_{t} - X^h_{t}|\right\rVert_p = \mathcal{O}(h^\alpha) \] ou encore \[ \left\lVert \sup_{k \in \{0,\dots,N\}} |x_{t_k} - X^h_{t_k}|\right\rVert_p = \mathcal{O}(h^\alpha) \]
Il s’agit d’une erreur trajectorielle, similaire à l’erreur étudiée en analyse numérique pour un schéma de discrétisation d’une ODE.
Il faut calculer \(\sup_k |x_{t_k} - X^h_{t_k}|\) sur le même événement \(\omega\) i.e. la même trajectoire Brownienne \((W_t)_{t \in [0,T]}\). C’est pour cela qu’on a choisit dans le design du code de séparer la simulation des accroissements Browniens dW et la construction des trajectoires pour chacun des schémas (exact, euler et milstein).
On fait un code pour calculer l’error forte dans \(\mathbf{L}^2\). On utilise un estimateur Monte Carlo de l’erreur \(\mathbf{L}^2\)
algo = { "euler": X.euler_simulation,
"milstein": X.milstein_simulation, }
strong_error = { key: [] for key in algo.keys() }
T = 5
n_paths = 5000
n_steps_values = T*2**np.arange(0, 10)
for n_steps in n_steps_values:
dW = math.sqrt(T/n_steps) * torch.randn((n_steps, n_paths))
sample_exact = X.exact_simulation(x0, T, dW)
for scheme in algo.keys():
sample_scheme = algo[scheme](x0, T, dW)
error = sample_exact - sample_scheme
max_error, _ = torch.max(error**2, axis=0)
mc = monte_carlo(max_error.cpu().numpy())
strong_error[scheme].append({
"mean": np.sqrt(mc["mean"]),
"lower": np.sqrt(mc["lower"]),
"upper": np.sqrt(mc["upper"]),
})
se_euler = pd.DataFrame(strong_error["euler"])
se_milstein = pd.DataFrame(strong_error["milstein"])dt = T / n_steps_values
fig, ax = plt.subplots(layout='tight')
ax.plot(dt, se_euler["mean"], color='C0', label='Euler')
ax.fill_between(dt, se_euler["lower"], se_euler["upper"],
color='C0', alpha=0.2, label='Euler Confident Interval')
ax.plot(dt, se_milstein["mean"], color='C1', label='Milstein')
ax.fill_between(dt, se_milstein["lower"], se_milstein["upper"],
color='C1', alpha=0.2, label='Milstein Confident Interval')
ax.legend()
ax.set_xlabel(fr'Timestep $h$')
ax.set_ylabel('Strong error')
ax.set_title(fr'Strong error with respect to timestep $h$')
plt.show()
On peut vérifier les coefficients \(\alpha = 1/2\) pour Euler et \(\alpha = 1\) pour Milstein en faisant une régression linéaire en échelle logarithmique. Par exemple en utilisant scikit-learn cela donne
from sklearn.linear_model import LinearRegression
x = np.log(dt)[:, np.newaxis]
y = np.log(se_euler["mean"])
model = LinearRegression().fit(x, y)
print(f"Strong convergence order, Euler scheme: {model.coef_.item():.4}")
y = np.log(se_milstein["mean"])
model = LinearRegression().fit(x, y)
print(f"Strong convergence order, Milstein scheme: {model.coef_.item():.4}")Strong convergence order, Euler scheme: 0.4973
Strong convergence order, Milstein scheme: 0.9694Comment estimer l’erreur forte d’un schéma sans connaître la vraie solution ?
Si \(\left\lVert \sup_k |x_{t_k} - X^h_{t_k}| \right\rVert_p \to 0\) alors \[ \left\lVert \sup_k |x_{t_k} - X^h_{t_k}| \right\rVert_p = \mathcal{O}(h^\alpha) \quad \Longleftrightarrow \quad \left\lVert \sup_k |X^h_{t_k} - X^{h/2}_{t_{2k}}| \right\rVert_p = \mathcal{O}(h^\alpha) \]
Attention encore une fois à bien construire les schémas \((X^h_{kh})_{k=0,\dots,N}\) et \((X^{h/2}_{\frac{kh}{2}})_{k=0,\dots,2N}\) en utilisant la même trajectoire brownienne, de donc le même objet dW.
On fera de même pour mettre en oeuvre l’extrapolation de Richardson-Romberg qui permet de gagner un ordre pour l’erreur faible.
algo = { "euler": X.euler_simulation,
"milstein": X.milstein_simulation, }
strong_error = { key: [] for key in algo.keys() }
T = 5
n_paths = 5000
n_steps_values = T*2**np.arange(1, 10)
for n_steps in n_steps_values:
dW_fine = math.sqrt(T/n_steps) * torch.randn((n_steps, n_paths))
dW_coarse = dW_fine[::2,:] + dW_fine[1::2,:]
for scheme in algo.keys():
sample_fine = algo[scheme](x0, T, dW_fine)
sample_coarse = algo[scheme](x0, T, dW_coarse)
error = sample_fine[::2,:] - sample_coarse
max_error, _ = torch.max(error**2, axis=0)
mc = monte_carlo(max_error.cpu().numpy())
strong_error[scheme].append({ "mean": np.sqrt(mc["mean"]),
"lower": np.sqrt(mc["lower"]),
"upper": np.sqrt(mc["upper"]),
})
se_euler_app = pd.DataFrame(strong_error["euler"])
se_milstein_app = pd.DataFrame(strong_error["milstein"])dt = T / n_steps_values
x = np.log(dt)[:, np.newaxis]
y = np.log(se_euler_app["mean"])
model = LinearRegression().fit(x, y)
coeff_euler = model.coef_.item()
y = np.log(se_milstein_app["mean"])
model = LinearRegression().fit(x, y)
coeff_milstein = model.coef_.item()
fig, ax = plt.subplots(layout='tight')
ax.plot(dt, se_euler_app["mean"], color='C0',
label=f'Euler, coeff={coeff_euler:.4}')
ax.fill_between(dt, se_euler_app["lower"], se_euler_app["upper"],
color='C0', alpha=0.2, label='Euler Confident Interval')
ax.plot(dt, se_milstein_app["mean"], color='C1',
label=f'Milstein, coef={coeff_milstein:.4}')
ax.fill_between(dt, se_milstein_app["lower"], se_milstein_app["upper"],
color='C1', alpha=0.2, label='Milstein Confident Interval')
ax.legend()
ax.set_xlabel(fr'Timestep $h$')
ax.set_ylabel('Approximated Strong error')
ax.set_title(fr'Approximated Strong error with respect to timestep $h$')
plt.show()
On rappelle que l’erreur faible entre la diffusion \(\bigl(x_t\bigr)_{t \in [0,T]}\) et un schéma numérique \(\bigl(X^N_{k}\bigr)_{k=0,\dots, N}\) le long d’une fonction \(\varphi\) vérifie le développement suivant \[ \exists c_1, \dots, c_R, \qquad \mathbb{E} \bigl[ \varphi(x_T) \bigr] - \mathbb{E} \bigl[ \varphi(X^N_N) \bigr] = \frac{c_1}{N} + \frac{c_2}{N^2} + \cdots + \frac{c_R}{N^R} + \mathcal{O}\bigl( \frac{1}{N^R} \bigr) \] Il faut de la régularité sur la fonction \(\varphi\) ou bien sur la densité de la loi de \(x_T\) pour obtenir ce résultat.
On illustre les résultats de convergence de l’erreur faible du schéma d’Euler et de Milstein dans le cas d’un pricing d’un call de strike \(K\) à maturité \(T\). Ainsi on a \[ \operatorname{d}\! x_t = r x_t \operatorname{d}\! t + \sigma x_t \operatorname{d}\! W_t, \quad x_0 > 0, \] et \[ P_{Call}(x_0, T, K) = e^{-r T} \mathbb{E}\bigl[ (x_T - K)_+ \bigr] = x_0 \phi(d_1) - K e^{-r T} \phi(d_2) \] avec \(d_1 = \log(S_0 / K) + \frac{(r+\sigma^2/2) T}{\sigma \sqrt{T}}\), \(d_2 = d_1 - \sigma \sqrt{T}\), et \(\phi\) la fonction de répartition d’une loi \(\mathcal{N}(0,1)\).
def bs_call(X, x0, T, K):
d1 = (np.log(x0/K) + (X.r+0.5*X.sigma**2)*T) / (X.sigma*math.sqrt(T))
d2 = d1 - X.sigma*math.sqrt(T)
price = x0*stats.norm.cdf(d1) - K*math.exp(-X.r*T)*stats.norm.cdf(d2)
return priceX = BlackScholes(r=0.2, sigma=0.3)
x0 = 100
T = 1
K = 100
exact_price = bs_call(X, x0, T, K)from torch.nn.functional import relu
algo = { "exact": X.exact_simulation,
"euler": X.euler_simulation,
"milstein": X.milstein_simulation, }
weak_error = { key: [] for key in algo.keys() }
n_paths = 100000
n_steps_values = T*2**np.arange(0, 10)
for n_steps in n_steps_values:
dW = math.sqrt(T/n_steps) * torch.randn((n_steps, n_paths))
for scheme in algo.keys():
sample_scheme = algo[scheme](x0, T, dW)
payoffs = math.exp(-X.r*T) * relu(sample_scheme[-1] - K)
mc = monte_carlo(payoffs.cpu().numpy() - exact_price)
weak_error[scheme].append(mc)
we_exact = pd.DataFrame(weak_error["exact"])
we_euler = pd.DataFrame(weak_error["euler"])
we_milstein = pd.DataFrame(weak_error["milstein"])dt = T / n_steps_values
fig, ax = plt.subplots(layout='tight')
ax.plot(dt, we_euler["mean"], color='C0', label='Euler')
ax.fill_between(dt, we_euler["lower"], we_euler["upper"],
color='C0', alpha=0.2, label='Euler Confident Interval')
ax.plot(dt, we_milstein["mean"], color='C1', label='Milstein')
ax.fill_between(dt, we_milstein["lower"], we_milstein["upper"],
color='C1', alpha=0.2, label='Milstein Confident Interval')
ax.plot(dt, we_exact["mean"], color='grey', label='exact')
ax.fill_between(dt, we_exact["lower"], we_exact["upper"],
color='grey', alpha=0.2, label='exact Confident Interval')
ax.legend()
ax.set_xlabel(fr'Timestep $h$')
ax.set_ylabel('Weak error')
ax.set_title(fr'Weak error with respect to timestep $h$')
plt.show()
L’extrapolation de Richardson est une technique simple, très utilisé en analyse numérique, permettant de gagner un ordre de convergence. En effet, en considérant une combinaison linéaire entre un schéma de pas \(\frac{T}{N}\) et un autre de pas \(\frac{T}{2N}\) on peut améliorer l’ordre de convergence en éliminant le premier terme \(\frac{c_1}{N}\) qui apparait dans le développement de l’erreur faible: \[ \mathbb{E} \bigl[ \varphi(x_T) \bigr] - \mathbb{E} \bigl[ 2 \varphi(X^{2N}_{2N}) - \varphi(X^{N}_{N}) \bigr] = -\frac{c_2}{2 N^2} + \mathcal{O}\bigl( \frac{1}{N^2} \bigr) \]
weak_error = { key: [] for key in algo.keys() }
n_paths = 100000
n_steps_values = T*2**np.arange(1, 10)
for n_steps in n_steps_values:
dW_fine = math.sqrt(T/n_steps) * torch.randn((n_steps, n_paths))
dW_coarse = dW_fine[::2,:] + dW_fine[1::2,:]
for scheme in algo.keys():
sample_fine = algo[scheme](x0, T, dW_fine)
sample_coarse = algo[scheme](x0, T, dW_coarse)
payoffs = math.exp(-X.r*T) * (
2*relu(sample_fine[-1] - K) - relu(sample_coarse[-1]-K)
)
mc = monte_carlo(payoffs.cpu().numpy() - exact_price)
weak_error[scheme].append(mc)
we_exact = pd.DataFrame(weak_error["exact"])
we_euler = pd.DataFrame(weak_error["euler"])
we_milstein = pd.DataFrame(weak_error["milstein"])dt = T / n_steps_values
fig, ax = plt.subplots(layout='tight')
ax.plot(dt, we_euler["mean"], color='C0', label='Euler')
ax.fill_between(dt, we_euler["lower"], we_euler["upper"],
color='C0', alpha=0.2, label='Euler Confident Interval')
ax.plot(dt, we_milstein["mean"], color='C1', label='Milstein')
ax.fill_between(dt, we_milstein["lower"], we_milstein["upper"],
color='C1', alpha=0.2, label='Milstein Confident Interval')
ax.legend()
ax.set_xlabel(fr'Timestep $h$')
ax.set_ylabel('Weak error')
ax.set_title(fr'Weak error with Richardson Extrapolation' +\
fr'with respect to timestep $h$')
plt.show()