Des dérivées sur une grille
Une dérivée est une limite, mais un ordinateur travaille avec des pas finis. Une différence finie remplace f'(x) par la pente sur un petit pas h : avant (f(x + h) - f(x)) / h, ou centrée (f(x + h) - f(x - h)) / (2h), plus précise car les erreurs des deux côtés se compensent (son erreur décroît comme h**2, celle de la différence avant seulement comme h). La même idée donne la dérivée seconde, la courbure : (f(x + h) - 2 f(x) + f(x - h)) / h**2. Pour un prix d'option en fonction du spot, c'est le gamma.
import math
def forward(f, x, h):
return (f(x + h) - f(x)) / h
def central(f, x, h):
return (f(x + h) - f(x - h)) / (2 * h)
exact = math.cos(1.0) # the derivative of sin at 1
for h in (0.1, 0.01, 0.001):
print(h, abs(forward(math.sin, 1.0, h) - exact), abs(central(math.sin, 1.0, h) - exact))# Gamma of a call from the second difference of its Black-Scholes price
import math
def norm_cdf(x):
return 0.5 * (1 + math.erf(x / math.sqrt(2)))
def bs_call(spot, strike, rate, sigma, t):
d1 = (math.log(spot / strike) + (rate + 0.5 * sigma ** 2) * t) / (sigma * math.sqrt(t))
return spot * norm_cdf(d1) - strike * math.exp(-rate * t) * norm_cdf(d1 - sigma * math.sqrt(t))
def second_difference(f, x, h):
return (f(x + h) - 2 * f(x) + f(x - h)) / h ** 2
price = lambda s: bs_call(s, 100, 0.05, 0.2, 1)
print(round(second_difference(price, 100, 0.5), 6)) # 0.018761
print(round(second_difference(price, 100, 0.01), 6)) # 0.018762D'une dérivée à une équation
L'intérêt de la méthode est de résoudre des équations aux dérivées partielles. Black–Scholes en est une : le prix de l'option évolue dans le temps selon sa courbure en fonction du spot. Après un changement de variables, elle devient l'équation de la chaleur, du/dt = D * d2u/dx2 : la chaleur se répand, et plus le profil de température est courbe, plus il évolue vite. On met les valeurs sur une grille, on remplace les deux dérivées par des différences, et l'équation devient une règle pour mettre à jour chaque point à partir de ses voisins.
# The explicit scheme: a spike of heat spreads to its neighbours, step by step
import numpy as np
def explicit_step(values, alpha):
new = values.copy()
new[1:-1] = values[1:-1] + alpha * (values[2:] - 2 * values[1:-1] + values[:-2])
return new
grid = np.zeros(9)
grid[4] = 1.0
for step in range(4):
print(step, np.round(grid, 3))
grid = explicit_step(grid, 0.25)En code, alpha = D * dt / dx**2 relie le pas de temps dt et le pas de grille dx. Une mise à jour, new[1:-1] = values[1:-1] + alpha * (values[2:] - 2 * values[1:-1] + values[:-2]), déplace d'un coup tous les points intérieurs grâce aux tranches ; les deux extrémités sont les conditions aux limites, fixées ici. Le schéma est simple, et le piège aussi.
La stabilité : le piège
Le schéma explicite n'est stable que si alpha <= 0.5, soit D * dt / dx**2 <= 0.5. Au-delà, les erreurs d'arrondi sont amplifiées à chaque pas et les valeurs explosent en dents de scie. La condition est dure : diviser le pas de grille par deux impose un pas de temps quatre fois plus petit, donc le coût augmente vite. C'est pourquoi le code professionnel utilise souvent des schémas implicites, qui résolvent un petit système linéaire à chaque pas et sont stables quel que soit dt ; l'explicite est le plus simple à écrire et à comprendre.
# The same spike with alpha = 0.4 (stable) and alpha = 0.6 (unstable) after 60 steps
import numpy as np
def explicit_step(values, alpha):
new = values.copy()
new[1:-1] = values[1:-1] + alpha * (values[2:] - 2 * values[1:-1] + values[:-2])
return new
for alpha in (0.4, 0.6):
grid = np.zeros(21)
grid[10] = 1.0
for _ in range(60):
grid = explicit_step(grid, alpha)
print(alpha, "largest value", float(np.abs(grid).max())) # 0.4: 0.0556, 0.6: about 3.1e7h les entrées d'un modèle et diviser est la manière dont les systèmes de risque calculent les grecques d'un valorisateur sans formule, et le choix de h est le même compromis qu'ici : trop grand, l'erreur de troncature apparaît ; trop petit, le bruit d'arrondi domine.- Une différence finie remplace une dérivée par des valeurs sur une grille : la centrée
(f(x + h) - f(x - h)) / (2h)a une erreur enh**2; la différence seconde(f(x + h) - 2 f(x) + f(x - h)) / h**2donne la courbure (gamma). - Les EDP deviennent des règles de mise à jour sur une grille : le pas explicite de l'équation de la chaleur est
values[i] + alpha * (values[i + 1] - 2 * values[i] + values[i - 1]), écrit avec des tranches. - Le schéma explicite exige
D * dt / dx**2 <= 0.5pour rester stable ; les schémas implicites lèvent cette limite au prix d'un système linéaire par pas.
À suivre : les matrices de covariance et de corrélation, première étape de l'optimisation de portefeuille.
Lab : ce qui est noté
- Estimer une courbure avec
second_difference25 pts - Avancer la grille d'un pas avec
explicit_step25 pts - Vérifier le pas de temps avec
is_stable25 pts - Avancer dans le temps avec
solve_heat25 pts