Auditing a Monte Carlo: Method Notes from a Leveraged-ETF Simulation

Four Monte Carlo scripts, one leveraged ETF pair (LUNR and its 2× daily fund LUNL), 879 + 149 daily bars of data — and a headline probability that moved from 46.8% to 30.9% once the scripts were audited. These are the method notes: every technique used to take the simulations apart, find what they were actually measuring, and rebuild them, with the real numbers each step produced. The subject is one ticker pair; the techniques apply to any simulation you’re about to trust.

1. Establish what your data actually is

Every downstream number inherits whatever your data source really means. Three things get assumed constantly and are worth five minutes of checking each.

1.1 Verify field semantics against a finer dataset

I needed to know whether the daily bar’s open was the 9:30 regular-session open or the first pre-market print — the entire overnight-gap analysis hinges on it. Rather than read documentation, I pulled minute bars for the same day and compared:

2026-08-13: daily-agg open 14.2900 | first pre-mkt 17.1500 | 09:30 RTH open 14.2900
            daily-agg close 17.5600 | 15:59 RTH close 17.5600 | last tick   17.6799

Daily open equals the 09:30 bar exactly, and the close is the 16:00 close, not the 19:57 extended-hours tick. Technique: cross-validate a coarse dataset against a finer one you can reconcile. One request settles a question that guesswork would leave open all project. (This is the same discipline our market-data guide demands of timestamps.)

1.2 Know precisely what “adjusted” adjusts

The API’s adjusted=true handles splits only — not dividends. For a distributing fund, price-only returns then understate total return, and that gap would masquerade as a cost. Since I was about to claim a large unexplained drag, I had to rule this out, so I queried the dividends endpoint directly: none for either ticker. Only then was the drag attributable to costs rather than measurement.

1.3 Your sample period is a modelling assumption

The same ticker gives wildly different parameters depending on the window:

Full history spans the 2023 SPAC period with a +251% single day. Neither window is “right” — but choosing one silently chooses your tail behaviour, and that choice deserves to be an explicit, named parameter rather than a hardcoded start date.

2. Measure the relationship, don’t assume it

The scripts assumed the ETF returned exactly 2× the stock minus a stated expense ratio. That is a hypothesis, and there were 148 days of data sitting right there to test it against.

2.1 What an OLS fit actually buys you

yₜ = β·xₜ + α + εₜ

Four separate deliverables from one line of code. β is realised exposure. α is systematic drift the exposure doesn’t explain — costs, in this case. sd(ε) is tracking error. R² tells you whether the linear story is the whole story.

beta, alpha = np.polyfit(x, y, 1)
resid       = y - (beta * x + alpha)
r2          = np.corrcoef(x, y)[0, 1] ** 2

# beta 1.9846   alpha -0.001735/day   resid_sd 0.007105/day   R2 0.9978

β is essentially 2 and R² is 0.998, so the fund tracks superbly. But α annualises to −43.7%/yr against the 1.31%/yr the script assumed — a factor of 33.

2.2 Annualisation: ×252 for means, ×√252 for volatility

Means add over time, so a daily mean scales linearly. Variance adds over independent periods, so standard deviation scales with the square root. Mixing these up is the most common unit error in finance code.

mean_annual = r.mean() * 252
vol_annual  = r.std(ddof=1) * np.sqrt(252)   # sqrt, because VARIANCE is what adds

Use ddof=1 for a sample standard deviation — the divisor is n−1, since estimating the mean from the same data costs one degree of freedom.

2.3 Attack a surprising number before believing it

−43.7%/yr was implausible enough that reporting it unchallenged would have been negligent. Four independent attacks, each capable of killing the finding:

Technique: for each explanation that would make your finding boring, run the test that would reveal it. Surviving four such tests is what converts a number into a claim.

2.4 Variance decomposition

Splitting close-to-close returns into overnight (4pm→9:30) and intraday (9:30→4pm) legs localises where risk actually lives: for LUNR, 14.2% overnight, 72.6% intraday. Note the shares don’t sum to 100% — because var(a+b) = var(a) + var(b) + 2cov(a,b), and simple returns compound rather than add. Always work in logs when you want additivity, and expect a cross-term when you don’t.

2.5 The encompassing regression — how to settle an argument

Two competing hypotheses about when the fund holds exposure. The weak approach fits each separately and compares R². The strong approach nests both in one model and lets them compete for the same variance:

LUNLcc = b₁·LUNR_overnight + b₂·LUNR_intraday + a

A fund holding nothing overnight must produce b₁ = 0. Measured:

X = np.column_stack([np.ones(n), overnight, intraday])
coef, *_ = np.linalg.lstsq(X, y, rcond=None)
se = np.sqrt(np.diag(np.linalg.inv(X.T @ X) * resid.var(ddof=3)))

# b1 (overnight) = 2.0175 +/- 0.0450     b2 (intraday) = 1.9832 +/- 0.0236

The t-statistic against the null is t = (estimate − null) / SE = 2.0175/0.0230 = 87.8 standard errors from zero. That is not a close call, and it is a far more decisive instrument than eyeballing two R² values.

3. Learn the arithmetic of the instrument

3.1 Work in log returns

With r_log = ln(1 + r_simple), returns become additive over time, which makes drift, variance and horizon scaling all linear. It has a second benefit that matters here: generating in log space makes a return below −100% impossible by construction, so you never need a floor to patch up a broken sampler.

log_r = mu + sd * z            # generate here
simple = np.expm1(log_r)       # convert out; always > -1

3.2 Volatility drag, derived

For geometric Brownian motion, Itô’s lemma gives the log dynamics of the stock, and the fund holds L× that exposure:

d ln S = (μ − ½σ²) dt + σ dW d ln Λ = (Lμ − ½L²σ²) dt + Lσ dW

Subtract L× the first from the second and the μ and dW terms cancel, leaving pure drag:

drag = −½σ²(L² − L) = −σ² for L = 2

No μ term: drag is unrelated to direction, and quadratic in leverage. A 3× fund suffers three times the drag of a 2× fund, not 1.5×. At LUNR’s σ = 1.107 that is −1.225/yr at 2× and −3.676/yr at 3×.

3.3 The subtlety: it’s mean-square, not variance

The discrete expansion is −½·E[r²]·(L²−L), and E[r²] = σ² + μ². Those coincide only when μ ≈ 0, which is why the textbook writes σ². Ten consecutive −5% days have zero variance, yet:

ln(2x fund) = -1.05361     2 * ln(underlying) = -1.02587     gap -0.02774
predicted:  -0.5 * mean(r**2) * (L**2 - L) * 10 = -0.02500

The fund lags 2× the log return even with constant returns. Drag is driven by the mean square of the daily move, not by variance in the statistical sense.

3.4 Two benchmarks — the source of nearly all confusion

“Does a leveraged ETF beat 2×?” has no answer until you say which 2×:

Against L× the log return the fund always lags — drag is unavoidable. Against L× the simple return the fund wins in trends and loses in chop, because sustained moves de-lever you into losses. Both statements are true; they are answers to different questions.

3.5 Terminal value is a product, so ordering cannot matter

∏(1 + L·rᵢ) is invariant under permutation. Multiplication commutes. Testing it on the fund’s real 148 returns:

terms = [np.prod(1 + rng.permutation(y)) for _ in range(20000)]

# true order 0.4592145015x    min 0.4592145015    max 0.4592145015    spread 2.05e-15

The −½σ²(L²−L) formula contains only a variance term — no autocorrelation, no sequencing. That is the same fact stated analytically.

3.6 …but the path is entirely about ordering

So ordering drives drawdown timing, barrier touches and stop-outs — everything you would actually trade against — and nothing about where you land if you hold. Report both kinds of statistic, and know which kind each one is. (Our equity-curve simulator makes the same point interactively: the fan is the strategy; one path is one draw.)

4. Generate futures honestly

4.1 The generator menu

The stationary bootstrap (Politis–Romano) draws geometric block lengths rather than fixed ones, so no artificial boundary recurs every b days. Cheap to implement and strictly better behaved than fixed blocks.

4.2 Verify the sampler does what you claim

Don’t assert that block sampling preserves volatility clustering — measure the autocorrelation of |r|:

|r| autocorrelation:   iid -0.001      block 0.394

4.3 Clustering acts through dispersion, not sequence

Given §3.5, how can clustering matter at all if order doesn’t? Because it changes the multiset each path draws. Block sampling lets one path pull an entire calm stretch and another an entire violent one, widening the spread of per-path variance — and drag keys off variance:

                 per-path realised vol       LUNL median   LUNL p05
iid bootstrap    mean 1.094  sd 0.149          7.744        1.399
block bootstrap  mean 1.085  sd 0.201          8.290        1.040

Same marginal distribution, materially fatter loss tail. That is the real argument for block over iid — and it is a dispersion argument, not an ordering one.

4.4 A bootstrap cannot invent a tail it never saw

worst day in the 252d sample          -15.93%
worst day the block bootstrap drew    -15.93%   # identical, by construction
worst day Student-t (df=3) drew       -99.99%

P(any day <= -50%, a 2x wipeout):   bootstrap 0.0000    student-t 0.0112

Every resampling method reports exactly zero wipeout probability — a property of your sample, not of the world. For tail questions you need a parametric generator that can extrapolate. Judgement still required in the other direction: Student-t with df = 3 produced a −99.99% day, which is its own kind of nonsense. Neither is truth; know which failure mode you’ve chosen.

4.5 Drift calibration: mean and median are different targets

For GBM, E[S_T] = S₀·e^(μT) but median = S₀·e^((μ−½σ²)T). At σ = 1.107 over 60 days the median sits 13.6% below the mean. Setting μ = ln(target/S₀)/T therefore aims the mean at your target while the median — the number people actually read — lands well short.

The robust fix avoids the algebra entirely: simulate first, then shift every log return by a constant so the finished sample lands where you want. Exact for either statistic, and immune to whatever else is in your generator.

log_r = np.log1p(sim_rets)
shift = (np.log(target / S0) - np.median(log_r.sum(axis=1))) / n_days
tilted = np.expm1(log_r + shift)          # median now lands exactly on target

4.6 Absorbing barriers

Flooring a daily return at −0.999 leaves a fund that lost 99.9% alive and able to rally back, inflating the right tail with paths that in reality would have been liquidated. Once dead, keep it dead:

dead = np.cumsum(r <= -1.0, axis=1) > 0
nav[:, 1:] = np.where(dead, 0.0, nav[:, 1:])

5. Quantify your own uncertainty

A Monte Carlo result is a sample statistic with its own sampling error. Printing “43.33%” without an error bar hides whether the third digit means anything.

5.1 Standard error of a simulated probability

SE = √( p(1−p) / n )

Verified empirically — 60 independent runs at n = 5,000:

observed sd of p_hat across runs   0.00674
formula sqrt(p(1-p)/n)             0.00706      # matches

The √n means halving your error costs 4× the paths. At 200,000 paths a probability near 50% carries ±0.22pp at two standard errors — which tells you a 3pp difference between scenarios is real and a 0.1pp difference is noise.

5.2 Quantiles have no closed form — bootstrap them

idx = rng.integers(0, len(x), size=(200, len(x)))
se  = np.std(np.quantile(x[idx], q, axis=1), ddof=1)

Resample your own output, recompute the statistic, take the spread. The same trick works for any statistic you can compute but not analyse.

5.3 Variance reduction: antithetic variates

For every draw z, also use −z. Symmetric error cancels in the pair rather than averaging down slowly, and it costs nothing:

half = (n + 1) // 2
z = rng.standard_normal((half, d))
return np.vstack([z, -z])[:n]

# sd of mean terminal, plain        0.01109
# sd of mean terminal, antithetic   0.00603    # 1.8x tighter, same cost

Its companion is common random numbers: when comparing scenarios, drive both with the same seed so the difference you measure is the assumption, not the sampling.

5.4 Report the sensitivity, not the point estimate

The most useful output of the whole exercise wasn’t a probability — it was finding which assumption owned it:

A 16-point swing from two assumptions the original scripts never surfaced. The bootstrap silently inherited the stock’s trailing +141%/yr drift — a strongly bullish input nobody had chosen on purpose. Any parameter that can move your answer that far belongs in a table, not buried in a default.

6. Prove the code does what you claim

6.1 Test invariants, not outputs

Asserting that a simulation returns 30.9% just freezes today’s bug. Assert properties that must hold for any correct implementation:

ok  view tilt lands the requested statistic exactly
ok  wipeout absorbing (dead path stays at 0 through later rallies)
ok  backtest aligned (exact 2x series reproduces to rmse=1.73e-16)
ok  all 5 generators match historical log-vol and stay above -100%
ok  |r| autocorr: iid=-0.001 vs block=0.394 (clustering preserved)
ok  antithetic normals cancel mean error exactly (sd 2.1e-19 vs 7.1e-04)

The third is the strongest pattern available: feed the model data it should reproduce perfectly — a synthetic exact-2× series — and demand machine precision. Any alignment or convention error shows up instantly.

6.2 Cross-check simulation against closed form

Wherever theory gives an answer, make the simulator reproduce it. Comparing measured drag against −½σ²(L²−L):

L=2.0: theory -1.2252/yr   simulated -1.1882/yr
L=3.0: theory -3.6757/yr   simulated -3.6774/yr

This validates the compounding, the sign conventions and the leverage handling in one shot.

6.3 Off-by-one errors hide in cumulative products

model_nav = np.r_[1.0, np.cumprod(1.0 + model_daily[1:])]   # both start at 1.0

Always verify that two series you compare start on the same date at the same value.

6.4 Vectorise — speed buys statistical quality

scalar Python double loop   ~28 s     (100,000 paths x 126 days)
vectorised numpy             ~0.2 s   # 140x

This isn’t cosmetic. One script ran 100 paths because more was too slow, and 100 paths cannot support a 5th percentile. Fast code is what lets you afford 200,000 paths, a sensitivity table and honest error bars.

6.5 Make the run reproducible

Seed the generator explicitly, cache fetched data to disk, keep configuration in one place, take secrets from the environment, and emit assumptions alongside results so the output can be re-derived months later.

The through-line

Nearly every finding came from the same reflex: asking what a number was actually measuring rather than what it was labelled. Four questions did most of the work.

1. Is this parameter assumed or measured? The 1.31%/yr cost was assumed; measuring it gave −43.7%.
2. What sample produced this, and what does that sample exclude? The window choice set the volatility, the tails, and whether wipeout risk existed at all.
3. Which statistic is this — mean, median, or something else? Drift calibration, benchmark comparison and drag all turn on this distinction.
4. How would I know if this were wrong? If no test could fail, it isn’t validation.

The scripts weren’t badly written. They were confidently answering a question slightly different from the one being asked — which is the normal failure mode of quantitative work, and why the audit is worth as much as the model. It’s the same standard we hold papers to, applied to our own code — and the reason the production checklist treats every unverified number as a lie until an independent mechanism confirms it.

From the corpus: risk measures beyond VaR → · trading simulation: replay vs agent-based vs generative →