Variables antitéticas

La precisión es cara

El error estándar de una estimación Monte Carlo cae como 1 / sqrt(N): diez veces más precisión exige cien veces más trayectorias. Cuando una simulación es costosa, una palanca mejor que más trayectorias es menos ruido por trayectoria. Son los métodos de reducción de varianza, y el más simple son las variables antitéticas. Usan dos veces los mismos números aleatorios: con un sorteo z también se ejecuta la trayectoria 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

Por qué ayuda emparejar

Si z es grande y positivo, la trayectoria termina alta y un call paga mucho; con -z termina baja y el call paga poco. Sus pagos están correlacionados negativamente, así que el promedio del par es mucho más estable que cada uno. Para un pago monótono como el de un call, el promedio del par (f(z) + f(-z)) / 2 tiene una varianza (1 + corr) / 2 veces la de un pago, con corr negativa. El precio sigue siendo el promedio de una cantidad insesgada, así que no hay sesgo.

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

El estimador y su barra de error

Trata cada par como una observación: pair = (payoff(z) + payoff(-z)) / 2. La estimación es la media descontada de los pagos de los pares, y el error estándar es la desviación típica descontada (ddof=1) de los pares entre sqrt(número de pares). El recuento es de *pares*, pues los dos miembros de un par no son independientes. Para juzgar la ganancia, compara con una estimación simple que use el mismo número de pagos de sorteos independientes.

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)

• Use cada sorteo z junto con -z; promedie los dos pagos en un pago de par, mucho menos ruidoso para un pago monótono.
• El error estándar se calcula sobre el número de pares, con ddof=1: discount * std(pairs) / sqrt(pairs).
• Compare con una estimación simple del mismo coste; no ayuda con pagos simétricos en z, y se combina con variables de control.

A continuación: las variables de control, usar una cantidad cuyo valor exacto se conoce para cancelar el ruido.

Lab: qué se califica

• Convertir sorteos en precios con terminal_price25 pts
• Informar la estimación y su error con mc_estimate25 pts
• Emparejar cada sorteo con su opuesto con antithetic_call_price25 pts
• Medir la ganancia con variance_ratio25 pts