Variables antithétiques

La précision coûte cher

L'erreur standard d'une estimation Monte Carlo décroît comme 1 / sqrt(N) : dix fois plus de précision exige cent fois plus de trajectoires. Quand une simulation coûte cher, un meilleur levier que plus de trajectoires est moins de bruit par trajectoire. Ce sont les méthodes de réduction de variance, dont la plus simple est celle des variables antithétiques. Elles utilisent deux fois les mêmes nombres aléatoires : avec un tirage z, on fait aussi tourner la trajectoire de -z.

import numpy as np

def terminal_price(s0, rate, sigma, t, z):
    return s0 * np.exp((rate - 0.5 * sigma ** 2) * t + sigma * np.sqrt(t) * z)

z = np.array([-1.0, 0.0, 1.0])
print(np.round(terminal_price(100.0, 0.05, 0.2, 1.0, z), 3))         # the price for z = -1, 0, 1
print(np.round(terminal_price(100.0, 0.05, 0.2, 1.0, -z), 3))        # and for the opposite draws: mirror images

Pourquoi l'appariement aide

Si z est grand et positif, la trajectoire finit haut et un call rapporte beaucoup ; avec -z, elle finit bas et le call rapporte peu. Leurs gains sont négativement corrélés, donc la moyenne de la paire est bien plus stable que chacun. Pour un gain monotone comme celui d'un call, la moyenne de la paire (f(z) + f(-z)) / 2 a une variance égale à (1 + corr) / 2 fois celle d'un gain, avec corr négatif. Le prix reste la moyenne d'une quantité sans biais, donc il n'y a pas de biais.

import numpy as np

def terminal_price(s0, rate, sigma, t, z):
    return s0 * np.exp((rate - 0.5 * sigma ** 2) * t + sigma * np.sqrt(t) * z)

z = np.random.default_rng(1).standard_normal(20000)
up = np.maximum(terminal_price(100.0, 0.05, 0.2, 1.0, z) - 100.0, 0.0)       # call payoffs for z
down = np.maximum(terminal_price(100.0, 0.05, 0.2, 1.0, -z) - 100.0, 0.0)    # ... and for -z
print(round(float(np.corrcoef(up, down)[0, 1]), 3))                           # negative: the payoffs offset each other

L'estimateur et sa barre d'erreur

Considérez chaque paire comme une observation : pair = (payoff(z) + payoff(-z)) / 2. L'estimation est la moyenne actualisée des gains des paires, et l'erreur standard est l'écart-type actualisé (ddof=1) des paires divisé par sqrt(nombre de paires). Le décompte porte sur les *paires*, car les deux membres d'une paire ne sont pas indépendants. Pour juger le gain, comparez à une estimation simple utilisant le même nombre de gains issus de tirages indépendants.

import numpy as np

def terminal_price(s0, rate, sigma, t, z):
    return s0 * np.exp((rate - 0.5 * sigma ** 2) * t + sigma * np.sqrt(t) * z)

def mc_estimate(payoffs, discount):
    return discount * payoffs.mean(), discount * payoffs.std(ddof=1) / np.sqrt(len(payoffs))

def call_payoff(z):
    return np.maximum(terminal_price(100.0, 0.05, 0.2, 1.0, z) - 100.0, 0.0)

rng = np.random.default_rng(3)
z, w = rng.standard_normal(5000), rng.standard_normal(5000)
discount = np.exp(-0.05)

plain = mc_estimate(call_payoff(np.concatenate([z, w])), discount)           # 10,000 independent payoffs
anti = mc_estimate((call_payoff(z) + call_payoff(-z)) / 2, discount)         # 5,000 pairs: the same 10,000 evaluations
print([round(float(v), 4) for v in plain], [round(float(v), 4) for v in anti])
print(round(float(plain[1] / anti[1]), 2), "times smaller error")

# A call profits from pairing; a straddle (symmetric in z) does not
import numpy as np

def terminal_price(z):
    return 100.0 * np.exp((0.05 - 0.5 * 0.2 ** 2) + 0.2 * z)

rng = np.random.default_rng(4)
z = rng.standard_normal(20000)
for name, payoff in (("call", lambda p: np.maximum(p - 100.0, 0.0)), ("straddle", lambda p: np.abs(p - 100.0))):
    pairs = (payoff(terminal_price(z)) + payoff(terminal_price(-z))) / 2
    plain = payoff(terminal_price(np.concatenate([z, rng.standard_normal(20000)])))
    ratio = (plain.std(ddof=1) / np.sqrt(len(plain))) / (pairs.std(ddof=1) / np.sqrt(len(pairs)))
    print(name, round(float(ratio), 2))                       # the call gains (1.41); the straddle gets worse (0.78)

• Utilisez chaque tirage z avec -z ; moyennez les deux gains en un gain de paire, bien moins bruité pour un gain monotone.
• L'erreur standard se calcule sur le nombre de paires, avec ddof=1 : discount * std(pairs) / sqrt(pairs).
• Comparez à une estimation simple de même coût ; aucun gain pour les gains symétriques en z, et elle se combine avec les variables de contrôle.

À suivre : les variables de contrôle, utiliser une quantité dont on connaît la valeur exacte pour annuler le bruit.

Lab : ce qui est noté

• Transformer les tirages en prix avec terminal_price25 pts
• Donner l'estimation et son erreur avec mc_estimate25 pts
• Apparier chaque tirage à son opposé avec antithetic_call_price25 pts
• Mesurer le gain avec variance_ratio25 pts