Estimation par Monte-Carlo

Répondre par des moyennes de tirages aléatoires

De nombreuses questions financières n'ont pas de formule simple : la probabilité qu'un portefeuille perde 20 % en un an, le prix d'une option dont le gain dépend de toute la trajectoire. Monte Carlo y répond par la force brute. Écrivez la quantité comme une moyenne, E[f(X)], simulez de nombreux tirages de X, appliquez f à chacun et prenez la moyenne. Par la loi des grands nombres, la moyenne d'échantillon converge vers l'espérance vraie. Une probabilité n'est que l'espérance d'un indicateur oui/non : « la fraction de simulations où l'événement a eu lieu » l'estime.

import math
import numpy as np

# Estimate pi: the quarter circle x^2 + y^2 <= 1 covers pi/4 of the unit square
rng = np.random.default_rng(42)
x, y = rng.random(100000), rng.random(100000)
inside = x ** 2 + y ** 2 <= 1                        # True/False for every point
print(round(4 * inside.mean(), 4))                    # 3.1378   close to 3.1416, with 100,000 points

Simuler des prix et des probabilités

Pour simuler une action, tirez des chocs normaux centrés réduits et faites-les passer par la formule log-normale du sujet précédent : les n prix terminaux sortent d'une seule expression sur tableau. Ensuite, toute question sur la distribution des prix devient une moyenne sur le tableau. La probabilité de finir sous 80 est np.mean(prices < 80), le prix espéré est prices.mean() et un quantile des pertes est np.quantile(...). Passez le générateur rng à la fonction au lieu de le créer à l'intérieur, pour que l'appelant garde le contrôle de la reproductibilité.

def simulate_terminal_prices(rng, s0, mu, sigma, t, n):
    z = rng.standard_normal(n)
    return s0 * np.exp((mu - 0.5 * sigma ** 2) * t + sigma * math.sqrt(t) * z)

prices = simulate_terminal_prices(np.random.default_rng(1), 100, 0.05, 0.2, 1, 200000)
print(round(float(np.mean(prices < 80)), 4))          # 0.1027   P(price below 80 after a year), about 10%
print(round(float(np.mean(prices > 130)), 4))         # 0.1216   P(price above 130)

Valoriser une option

Un call européen paie max(S - K, 0) à l'échéance. Dans le monde risque-neutre, tout actif croît au taux sans risque rate, donc le juste prix est le gain espéré actualisé : simulez les prix terminaux avec mu = rate, calculez le gain de chacun, faites la moyenne et multipliez par exp(-rate * t). Pour cette option simple, une formule fermée existe (Black-Scholes, 10,4506 pour les chiffres ci-dessous) et la simulation s'en approche. L'intérêt de Monte Carlo est que les mêmes quelques lignes fonctionnent encore pour des gains sans aucune formule.

def mc_call_price(rng, s0, strike, rate, sigma, t, n):
    prices = simulate_terminal_prices(rng, s0, rate, sigma, t, n)    # risk-neutral: the price grows at the risk-free rate
    payoff = np.maximum(prices - strike, 0)                           # max(S - K, 0) for every simulated path
    return math.exp(-rate * t) * payoff.mean()                        # discounted average payoff

price = mc_call_price(np.random.default_rng(1), 100, 100, 0.05, 0.2, 1, 200000)
print(round(price, 4))        # 10.3993   the exact Black-Scholes price is 10.4506: the difference is Monte Carlo error

# More points, less error (1 / sqrt(n)): the estimate of pi with the same seed
for n in (100, 10000, 1000000):
    g = np.random.default_rng(42)
    x, y = g.random(n), g.random(n)
    print(n, round(4 * float((x * x + y * y <= 1).mean()), 4))
# 100 3.36
# 10000 3.1488
# 1000000 3.1424

• Monte Carlo : écrivez la quantité comme une moyenne E[f(X)], simulez de nombreux tirages, appliquez f, prenez la moyenne.
• Une probabilité est la moyenne d'une comparaison : np.mean(prices < 80) ; passez le générateur en argument pour que l'exécution soit reproductible.
• Prix d'option = exp(-rate * t) * mean(max(S - K, 0)) avec les prix simulés à mu = rate ; donnez aussi l'erreur.

À suivre : l'erreur de Monte Carlo et la convergence, la précision d'une estimation et comment obtenir plus pour moins.

Lab : ce qui est noté

• Lancer des fléchettes sur un carré pour trouver pi avec estimate_pi25 pts
• Compter la fréquence d'un événement avec estimate_probability25 pts
• Simuler où une action pourrait finir avec simulate_terminal_prices25 pts
• Valoriser une option par simulation avec mc_call_price25 pts