Monte Carlo Methods

Monte Carlo methods

A standard error measures how much the estimator wobbles. It is blind to whatever the estimator is converging to. Variance reduction makes a wrong number precise faster; this skill is about the order of operations — get the target right, then reduce the variance.

../../../fin-models/skills/option-pricing-models/SKILL.md §3 covers the Euler-vs-exact discretisation of a GBM. This skill is the estimator machinery around it: reduction factors you can measure, American exercise, QMC, and a bias that survives 600,000 paths.

Every number marked ✅ Measured is printed by scripts/monte_carlo.py (numpy 2.2.6 + scipy 1.13.0, seed 20260909, about 45 s). ✅ source-verified means read in the cited paper or in the installed scipy source.

1. ✅ The rate, and three reductions with their factors

The test payoff is a discretely monitored arithmetic-average Asian call, S=K=100, T=1, r=5%, σ=25%, 12 monitoring dates.

✅ Measured — 4× the paths halves the error bar, three times over: se 0.1071 → 0.0538 → 0.0267 → 0.0133, ratios 1.989, 2.014, 2.008 against the 2.000 that 1/sqrt(N) predicts. To cut the error by 10 you pay 100×, and none of it touches a bias.

✅ Measured — at a fixed 100,000-path budget:

• 🔑 The control variate is not in the same league as the other two. The geometric-average Asian is exactly lognormal, so it has a closed form (⚠️ Kemna & Vorst 1990) — 6.99073 here — and it correlates ~0.999 with the arithmetic payoff. One extra exponential per path buys 859× the variance, i.e. the same accuracy as 86 million plain paths.
• 🚨 Antithetic sampling halves the number of independent units. The error bar must be computed over pair averages. Treating 100,000 antithetic draws as 100,000 independent samples understates the error and flatters the method — the same trap ../../../fin-models/skills/option-pricing-models/SKILL.md §6 finds in QuantLib's requiredSamples.
• 🚨 A stratified estimator with one draw per stratum has no honest standard error. The marginal law is unchanged, so the naive sample sd reports the plain-MC number and the 4.0× disappears. asian_mc uses two draws per stratum and pools the within-stratum variances, which is the only way to see the reduction you actually got.

2. ✅ Longstaff-Schwartz, reproduced from the paper

Longstaff & Schwartz (2001), Valuing American Options by Simulation: A Simple Least-Squares Approach, RFS 14(1), 113–147.

2a. The eight-path example of §1, exactly

✅ Measured against the paper's own printed output:

🔑 The paper never prints the price matrix — only the in-the-money X and Y columns. The matrix in LS_PATHS is the unique one consistent with all of them, and getting both regressions and both prices back out is the proof that it is right. ⚠️ The regressions are in the raw state, not in moneyness; lsm_price(..., scale=1.0) is what reproduces them.

2b. Rows of Table 1 — and a contract mismatch that dwarfs the Monte Carlo error

✅ source-verified — K=40, r=6%, "exercisable 50 times per year", "100,000 (50,000 plus 50,000 antithetic) paths", basis "a constant and the first three Laguerre polynomials".

Three unrelated methods in the same place: the paper's finite difference, an independent CRR tree, and LSM.

✅ And the tree itself is not new code taken on trust — option_models_cross_check() runs it against the separate implementation in ../../../fin-models/skills/option-pricing-models/SKILL.md, on that skill's own case (S=K=100, T=1, r=5%, q=2%, σ=20%, American put): 6.659528281882490 here against 6.659528281882506 there at 800 steps, a difference of −1.6e-14, and that skill documents 6.659528. It returns None on a bare install, so the module still runs standalone.

🚨 The last column is the trap. It is the same tree with exercise allowed at all 2,000 steps — a continuously exercisable American, worth 0.003 to 0.009 more. That is a different contract, not an error, and it is several times the paper's own standard error. A tree that allows exercise at every node will always look like it beats an LSM Bermudan. Compare a Bermudan against a Bermudan (bermudan_tree_put restricts exercise to the same 50 dates a year), or the comparison measures the exercise schedule instead of the method.

2c. 🚨 Footnote 9: max(intrinsic, fitted continuation) is upward biased

✅ source-verified — LS footnote 9: "if the American option were valued by taking the maximum of the immediate exercise value and the estimated continuation value, and discounting this value back, the resulting American option value could be severely upward biased. This bias arises since the maximum operator is convex; measurement error in the estimated continuation value results in the maximum operator being upward biased."

✅ Measured — same paths, same seed, same regressions; only the value carried backwards changes, so the difference is the rule and nothing else:

🚨 Positive at every strike, and no path count fixes it — it is a property of the estimator, not of the sample. max(a, b̂) with a noisy b̂ is biased up even when b̂ is unbiased. Carry the realised discounted cash flow backwards, which is what the paper does.

2d. 🚨 Two more one-line ways to lose a cent — or six

✅ Measured, S=36, σ=.20, T=1, against the Bermudan tree's 4.4780:

• ✅ source-verified — LS §1: "We use only in-the-money paths since it allows us to better estimate the conditional expectation function in the region where exercise is relevant and significantly improves the efficiency of the algorithm." Regressing on all paths costs 12× the error of doing it their way, because a three-term fit spends its degrees of freedom on a region where the exercise decision is never close.
• 🚨 The weighted Laguerre basis exp(−X/2)·L_k(X) underflows on a raw stock price. At S = 40, exp(−S/2) = 2e-9, so all three columns are numerically zero beside the constant and the fit degenerates to an intercept — −0.0634, as bad as regressing on everything, and it does not warn. Feed it the moneyness S/K. This is why so many published LSM implementations quietly use plain polynomials.
• ✅ In-sample 4.4785 against out-of-sample 4.4818, a gap of +0.0034 — the paper's own Table 2 diagnostic (fit the stopping rule on one set of paths, value it on another), and the paper's conclusion holds: they are "virtually identical".
• 🔑 LSM's stopping rule is estimated, therefore suboptimal, therefore the out-of-sample value is a lower bound on the true price — IN EXPECTATION. ⚠️ It is not a bound on any single run: at 100,000 paths the Monte Carlo noise is several times the bias, which is why LS ran the diagnostic five times per row. An LSM number above the tree is not by itself evidence of anything, and a test that asserts a strict inequality will fail on a seed.

3. ✅ Quasi-Monte Carlo: scrambled Sobol

✅ source-verified in the installed scipy 1.13.0 — the constructor is Sobol(d, *, scramble=True, bits=None, seed=None, optimization=None), "If True, use LMS+shift scrambling", "Max dimensionality is 21201", and the class carries the warning: "Sobol' sequences are a quadrature rule and they lose their balance properties if one uses a sample size that is not a power of 2, or skips the first point, or thins the sequence". ✅ Measured — Sobol(d=2).random(5) raises UserWarning: The balance properties of Sobol' points require n to be a power of 2.

🚨 ✅ Measured — Sobol(scramble=False).random(2)[0] is exactly [0. 0. 0. 0.], and norm.ppf(0) is -inf. Turning scrambling off hands you an infinite first path. The scrambled first point is [0.973 0.286 0.229 0.689]. scramble=True is the default; leave it on.

✅ Measured — RMS error over 8 repetitions against a control-variate reference (not a QMC reference — a QMC reference would share its points with the thing being measured):

✅ Fitted convergence rates: MC N^-0.55, QMC N^-0.66 (theory: −0.50, and up to −1 for a smooth integrand). ⚠️ The advantage is real and it is not the textbook N^-1: this is a 12-dimensional integrand with a kink at the strike, and kinks are what cost QMC its rate.

🚨 A single QMC run has no error estimate at all — the points are deterministic, so there is nothing to take a variance of. The error bar comes from independent scrambles: ✅ Measured, 2^14 points × 16 scrambles gives 7.31185 ± 0.00153. That is the entire point of randomised QMC. Use random_base2(m); slicing or thinning the sequence throws away the balance property that was the reason to use it.

4. 🚨 The trap: a bias that 600,000 paths cannot see

A down-and-out call, S=K=100, H=90, T=1, r=5%, σ=25%. ✅ Continuous monitoring has a closed form: 9.11122 (against the vanilla call's 12.33600). Simulate it by checking the barrier at m dates — exact lognormal steps, so there is no Euler bias anywhere — and every path that dips under H between two dates and comes back survives.

✅ Measured:

• 🚨 Read the two m = 52 rows against each other. 4× the paths cuts the error bar 0.0392 → 0.0197 and the price stays 0.8 too high, so the bias goes from 21 to 44 standard errors. More paths make the reported error bar a worse description of how wrong the number is.
• ✅ Monitoring dates are what fix it: m = 12 → 252 cuts the bias 1.58 → 0.39, a factor of 4.1 against the 4.6 that 1/sqrt(m) predicts.
• ✅ The last two columns are Broadie-Glasserman-Kou's continuity correction: evaluate the continuous formula at the shifted barrier H·exp(−β·σ·sqrt(T/m)), with β = −ζ(1/2)/sqrt(2π) = 0.582597. It turns a 40-sigma bias into under one sigma at every m, out of one closed-form evaluation and no extra paths.
• 🚨 ✅ scipy.special.zeta(0.5, 1) returns nan — the Hurwitz form carries no analytic continuation below 1. scipy.special.zetac(0.5) + 1 does, and gives −1.4603545088095866 (its docstring: "For x < 1 the analytic continuation is computed"). A nan there does not raise: it flows through exp() into a barrier level and surfaces much later as a comparison that is inexplicably False. That is how it was found here.

The rule this section exists for: report the discretisation and the standard error separately, and pick the discretisation before you spend anything on variance reduction. A control variate at m = 52 would have given you 9.97 ± 0.001 — three digits of precision around a number that is wrong in the first one.

5. What the script gives you

scripts/monte_carlo.py — numpy + scipy only, importable, no file writes, no network:

⚠️ What is not here: importance sampling (the right tool when the payoff is a rare event — §4's barrier at H=50 would need it), multilevel Monte Carlo, Brownian-bridge and PCA path constructions for QMC (which are what recover the N^-1 rate in high dimensions), Greeks by pathwise or likelihood-ratio methods, and the upper-bound (duality) estimators that bracket LSM's lower bound.

Where this sits

• ../../../fin-models/skills/option-pricing-models/SKILL.md — the models themselves: Heston's branch cut, CRR's parity oscillation, SABR, and the Euler-vs-exact discretisation bias on a GBM. §4 here is the same lesson on a payoff whose discretisation is in the contract.
• ../../../fin-models/skills/risk-measures-var-cvar/SKILL.md — when the simulation output is a risk number rather than a price, and the tail quantile has its own error bar.
• ../copulas-and-dependence/SKILL.md — what you simulate from. Its §6 VaR numbers are Monte Carlo and inherit everything above.
• ../hawkes-processes/SKILL.md §5 — the same warning from the other side: a standard error is only as good as the independence assumption underneath it, and clustered arrivals break it by a measured factor.
• ../../../fin-core/skills/backtest-validation/SKILL.md — a backtest is a Monte Carlo with one path, and §4's moral applies hardest there.