Antithetische Variablen

Genauigkeit ist teuer

Der Standardfehler einer Monte-Carlo-Schätzung sinkt wie 1 / sqrt(N): Zehnmal mehr Genauigkeit braucht hundertmal mehr Pfade. Ist eine Simulation teuer, ist ein besserer Hebel als mehr Pfade weniger Rauschen pro Pfad. Das sind Varianzreduktionsverfahren, das einfachste sind antithetische Variablen. Sie nutzen dieselben Zufallszahlen zweimal: Zu einer Ziehung z wird auch der Pfad für -z gerechnet.

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

Warum die Paarung hilft

Ist z groß und positiv, endet der Pfad hoch, und ein Call zahlt viel; mit -z endet er tief, und der Call zahlt wenig. Ihre Auszahlungen sind negativ korreliert, der Paardurchschnitt ist also viel stetiger als jede einzelne. Bei einer monotonen Auszahlung wie der eines Calls hat der Paardurchschnitt (f(z) + f(-z)) / 2 die Varianz (1 + corr) / 2 mal die einer Auszahlung, mit negativem corr. Der Preis bleibt der Mittelwert einer unverzerrten Größe, es entsteht kein Bias.

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

Der Schätzer und sein Fehlerbalken

Behandeln Sie jedes Paar als eine Beobachtung: pair = (payoff(z) + payoff(-z)) / 2. Die Schätzung ist der abgezinste Mittelwert der Paar-Auszahlungen, der Standardfehler die abgezinste Standardabweichung (ddof=1) der Paare durch sqrt(Zahl der Paare). Gezählt werden die *Paare*, da die beiden Glieder eines Paares nicht unabhängig sind. Um den Gewinn zu beurteilen, vergleicht man mit einer einfachen Schätzung, die dieselbe Zahl von Auszahlungen aus unabhängigen Ziehungen nutzt.

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)

• Nutzen Sie jede Ziehung z zusammen mit -z; mitteln Sie beide Auszahlungen zu einer Paar-Auszahlung, die bei monotoner Auszahlung viel weniger rauscht.
• Der Standardfehler wird über die Zahl der Paare mit ddof=1 berechnet: discount * std(pairs) / sqrt(pairs).
• Vergleichen Sie mit einer einfachen Schätzung gleicher Kosten; bei in z symmetrischen Auszahlungen nutzlos, und kombinierbar mit Kontrollvariablen.

Als Nächstes: Kontrollvariablen, eine Größe mit bekanntem exakten Wert nutzen, um Rauschen aufzuheben.

Lab: Was bewertet wird

• Mit terminal_price Ziehungen in Preise verwandeln25 Pkt.
• Mit mc_estimate Schätzung und Fehler angeben25 Pkt.
• Mit antithetic_call_price jede Ziehung mit ihrem Gegenstück paaren25 Pkt.
• Mit variance_ratio den Gewinn messen25 Pkt.