Barraquand moment matching

A finite sample of normals never has mean exactly 0 and variance exactly 1, and the payoff inherits those sampling errors. Jérôme Barraquand (1995), pricing multi-asset options, simply rescaled the sample,

\[\tilde z_i = \frac{z_i - \bar z}{s_z},\]

so the first two moments are exact (“quadratic resampling”). The draws become dependent, so error bars come from independent replications.

import matplotlib.pyplot as plt
import numpy as np

from financialkit.monte_carlo import moment_match, moment_matching_estimate, sample_mean, terminal_prices, vanilla_payoff
from financialkit.monte_carlo.visualizers import plot_estimator_comparison
from financialkit.options import black_scholes_price

S, K, T, r, sigma = 100.0, 100.0, 1.0, 0.05, 0.2
exact = black_scholes_price(S, K, T, r, sigma)


def call(z):
    return np.exp(-r * T) * vanilla_payoff(terminal_prices(S, T, r, sigma, z), K)

Matched samples have exact moments

z = np.random.default_rng(0).standard_normal(500)
w = moment_match(z)
print(f"raw:     mean {z.mean():+.4f}, std {z.std():.4f}")
print(f"matched: mean {w.mean():+.4f}, std {w.std():.4f}")

raw:     mean -0.0269, std 1.0136
matched: mean -0.0000, std 1.0000

Twenty replications of 2000 draws each

n, reps = 2000, 20
gen = np.random.default_rng(1)
crude = sample_mean([call(gen.standard_normal(n)).mean() for _ in range(reps)], method="crude", n_evaluations=n * reps)
matched = moment_matching_estimate(call, n, n_replications=reps, rng=2)
print(f"crude {crude.estimate:.4f} +/- {crude.std_error:.4f}; matched {matched.estimate:.4f} +/- {matched.std_error:.4f}; exact {exact:.4f}")

fig, axes = plt.subplots(1, 2, figsize=(11, 4))
plot_estimator_comparison([crude, matched], exact=exact, ax=axes[0])

crude 10.2772 +/- 0.0496; matched 10.4502 +/- 0.0072; exact 10.4506

<Axes: ylabel='estimate (95% CI); variance reduction'>

The gain is largest for small samples

ns = np.array([100, 300, 1000, 3000, 10000])
ratio = []
for k in ns:
    g = np.random.default_rng(int(k))
    c = sample_mean([call(g.standard_normal(k)).mean() for _ in range(reps)], n_evaluations=k * reps)
    ratio.append(moment_matching_estimate(call, int(k), n_replications=reps, rng=int(k) + 1).variance_ratio(c))
axes[1].semilogx(ns, ratio, "o-")
axes[1].set_xlabel("draws per replication")
axes[1].set_ylabel("variance reduction factor")
fig.tight_layout()
plt.show()