nb1.1 — LIE / LTV by Simulation

Companion to Ch 1.1, “Joint, Conditional, and the Two Laws That Run Everything.”

In the chapter we proved two identities with pencil and paper, on Sam’s calm/stormy regime table:

A proof is airtight, but it can still feel like a magic trick. This notebook does the thing that turns a trick into a fact you trust: we draw hundreds of thousands of (X, Y) pairs from Sam’s exact distribution, estimate everything from the raw sample by simple grouping, and watch the simulated numbers march toward the analytic answers as the sample grows. No appeal to the theorem — just data and arithmetic landing on the same place the algebra promised.

Setup

We fix a seed so every run of this notebook produces the identical sample — the reproducibility rule from CONVENTIONS.md. We also force matplotlib’s non-interactive Agg backend so the notebook runs headless (on a server, in CI, anywhere) without a display.

import matplotlib
matplotlib.use("Agg")  # headless, non-interactive backend

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt

rng = np.random.default_rng(42)  # one generator, seeded, for the whole notebook

print("numpy", np.__version__, "| pandas", pd.__version__)

Sam’s universe, as a table

Recall Sam’s joint distribution of (X, Y), where Y is the market regime (calm or stormy) and X is the day’s momentum-portfolio return in percent. Each cell is the joint probability p(x, y) = \Pr(X = x,\, Y = y).

X = -2\% X = 0\% X = +2\% row total
calm (Y=1) 0.05 0.35 0.30 0.70
stormy (Y=2) 0.15 0.10 0.05 0.30
col total 0.20 0.45 0.35 1.00

We encode it exactly as in the chapter’s code snippet: a 2\times 3 array of joint probabilities, with the three possible return values alongside.

# Return values (in %) and the joint probability table.
x_vals = np.array([-2.0, 0.0, 2.0])

# rows = regime (0: calm, 1: stormy); columns = return value
joint = np.array([[0.05, 0.35, 0.30],    # calm
                  [0.15, 0.10, 0.05]])   # stormy

assert np.isclose(joint.sum(), 1.0), "joint must sum to 1"

p_y = joint.sum(axis=1)   # marginal of Y (regime): [0.70, 0.30]
p_x = joint.sum(axis=0)   # marginal of X (return): [0.20, 0.45, 0.35]

print("p_y (calm, stormy) =", p_y)
print("p_x (-2, 0, +2)    =", p_x)

The analytic answers (the targets to hit)

Before we simulate anything, let’s recompute the exact numbers from the table the same way the chapter’s snippet did. These are the bullseyes our Monte Carlo darts should cluster around. Working with the true probabilities there is no sampling noise — the law holds to machine precision.

# Conditional distribution of X within each regime: divide each row by its total.
p_x_given_y = joint / p_y[:, None]

# Conditional mean and variance of X in each regime.
cond_mean = (p_x_given_y * x_vals).sum(axis=1)                        # E[X | Y]
cond_var  = (p_x_given_y * (x_vals - cond_mean[:, None])**2).sum(1)   # Var(X | Y)

# The two LTV pieces, averaged / varied over the regime distribution p_y.
within  = (p_y * cond_var).sum()                  # E[Var(X|Y)]
mean_X  = (p_y * cond_mean).sum()                 # E[X] via LIE
across  = (p_y * (cond_mean - mean_X)**2).sum()   # Var(E[X|Y])

# Direct total variance from the marginal of X.
var_X = (p_x * (x_vals - mean_X)**2).sum()

print(f"E[X | calm]         = {cond_mean[0]:.3f}   (chapter: 0.714)")
print(f"E[X | stormy]       = {cond_mean[1]:.3f}   (chapter: -0.667)")
print()
print(f"E[X]                = {mean_X:.3f}   (chapter: 0.300)")
print(f"within  E[Var(X|Y)] = {within:.3f}   (chapter: 1.710)")
print(f"across  Var(E[X|Y]) = {across:.3f}   (chapter: 0.400)")
print(f"within + across     = {within + across:.3f}   (chapter: 2.110)")
print(f"Var(X) direct       = {var_X:.3f}   (chapter: 2.110)")
assert np.isclose(within + across, var_X)   # the law, checked exactly

Now draw a sample

Here is the bridge from “true distribution” to “data.” We have six joint cells, each with a known probability. Sampling a day means picking one of those six cells with the right odds, then reading off its regime and return. We flatten the 2\times 3 table into six outcomes, draw cell indices with rng.choice, and unflatten back into a (regime, return) pair for each draw.

This is exactly what a real dataset is: not the probabilities, but realized draws from them. Everything downstream pretends we never saw the table — we only get the sample.

def draw_sample(n, generator):
    '''Draw n (regime, return) pairs from Sam's joint distribution.

    Returns a tidy DataFrame with columns 'regime' (0=calm, 1=stormy) and 'ret' (%).
    '''
    flat_p = joint.ravel()                         # 6 cell probabilities
    # Index grid matching the flattened table: which regime row, which return column.
    regime_grid = np.repeat(np.arange(joint.shape[0]), joint.shape[1])  # [0,0,0,1,1,1]
    ret_grid    = np.tile(x_vals, joint.shape[0])                       # [-2,0,2,-2,0,2]

    cells = generator.choice(flat_p.size, size=n, p=flat_p)
    df = pd.DataFrame({"regime": regime_grid[cells],
                       "ret":    ret_grid[cells]})
    return df

sample = draw_sample(200_000, rng)
print(sample.head())
print("\nsample size:", len(sample))
print("\nrealized regime frequencies:")
print((sample["regime"].value_counts(normalize=True).sort_index()
       .rename({0: "calm", 1: "stormy"})))

Estimate everything from the raw sample by grouping

Now we throw away the table and act like data analysts. We groupby the regime and compute, from the sample, the within-regime means and variances — the empirical stand-ins for \mathbb{E}[X\mid Y] and \operatorname{Var}(X\mid Y). Then we recombine them exactly as the two laws prescribe.

Two small but important choices:

  • We use the population variance (ddof=0, dividing by n) so the estimates target the probability-weighted variances in the laws, not the bias-corrected sample variance.
  • The regime weights are the realized frequencies in the sample, not the true p_y — because in real life you never know p_y, you only count how often each regime showed up.
g = sample.groupby("regime")["ret"]

sim_cond_mean = g.mean()                  # E[X | Y], estimated
sim_cond_var  = g.var(ddof=0)             # Var(X | Y), estimated (population variance)
sim_weight    = g.size() / len(sample)    # realized regime frequencies

# LIE:  E[X] = E[ E[X|Y] ]  -> weight the conditional means by regime frequency
sim_mean_X = (sim_weight * sim_cond_mean).sum()

# LTV pieces
sim_within = (sim_weight * sim_cond_var).sum()                       # E[Var(X|Y)]
sim_across = (sim_weight * (sim_cond_mean - sim_mean_X)**2).sum()    # Var(E[X|Y])

# Direct total variance, straight off the pooled sample (ignoring regime entirely).
sim_var_X = sample["ret"].var(ddof=0)

print("--- estimated from a 200,000-draw sample ---")
print(f"E[X | calm]    = {sim_cond_mean[0]:+.4f}   (target  0.714)")
print(f"E[X | stormy]  = {sim_cond_mean[1]:+.4f}   (target -0.667)")
print()
print(f"E[X]                = {sim_mean_X:.4f}   (target 0.300)")
print(f"within E[Var(X|Y)]  = {sim_within:.4f}   (target 1.710)")
print(f"across Var(E[X|Y])  = {sim_across:.4f}   (target 0.400)")
print(f"within + across     = {sim_within + sim_across:.4f}   (target 2.110)")
print(f"Var(X) pooled       = {sim_var_X:.4f}   (target 2.110)")

Reveal the trick: LIE and LTV are not assumptions, they are bookkeeping

Look at the two pairs of lines above. The simulated E[X] from averaging the two regime means matches the simulated mean you’d get by ignoring regimes entirely. And the simulated within + across matches the simulated pooled Var(X). We never imposed either law — we computed the left side and the right side by completely different routes through the same sample, and they agree.

That agreement is the content of the laws. The “magic” is just that splitting the data into groups and recombining loses nothing, as long as you weight the groups by how often they occur. Let’s make the within/across split visible.

fig, ax = plt.subplots(figsize=(6, 4))

labels = ["within\nE[Var(X|Y)]", "across\nVar(E[X|Y])"]
sim_pieces = [sim_within, sim_across]
true_pieces = [within, across]

xpos = np.arange(len(labels))
w = 0.38
ax.bar(xpos - w/2, true_pieces, w, label="analytic", color="#4C72B0")
ax.bar(xpos + w/2, sim_pieces, w, label="simulated", color="#DD8452")

for i, (t, s) in enumerate(zip(true_pieces, sim_pieces)):
    ax.text(i - w/2, t + 0.02, f"{t:.3f}", ha="center", va="bottom", fontsize=9)
    ax.text(i + w/2, s + 0.02, f"{s:.3f}", ha="center", va="bottom", fontsize=9)

ax.axhline(var_X, ls="--", color="gray", lw=1)
ax.text(1.45, var_X + 0.02, f"Var(X) = {var_X:.3f}", ha="right", va="bottom",
        fontsize=9, color="gray")
ax.set_xticks(xpos)
ax.set_xticklabels(labels)
ax.set_ylabel("variance contribution")
ax.set_title("LTV split: within + across = total")
ax.legend()
fig.tight_layout()
fig.savefig("ltv_split.png", dpi=110)
print("saved ltv_split.png")
print(f"within is {100*within/var_X:.0f}% of total risk; across is {100*across/var_X:.0f}%")

Convergence: watch the estimates settle onto the targets

A single 200,000-draw estimate being close could be luck. The real claim is that as the sample grows, the estimates converge to the analytic values — the law of large numbers doing its job underneath the two identities. We sweep the sample size over several orders of magnitude and, for each size, recompute the full LTV decomposition from a fresh draw.

def ltv_estimates(n, generator):
    '''Estimate (E[X], within, across, within+across, Var(X)) from an n-draw sample.'''
    df = draw_sample(n, generator)
    gg = df.groupby("regime")["ret"]
    cm = gg.mean()
    cv = gg.var(ddof=0)
    wt = gg.size() / len(df)
    mX = (wt * cm).sum()
    win = (wt * cv).sum()
    acr = (wt * (cm - mX)**2).sum()
    return mX, win, acr, win + acr, df["ret"].var(ddof=0)

sizes = [100, 300, 1_000, 3_000, 10_000, 30_000, 100_000, 300_000, 1_000_000]

rows = []
conv_rng = np.random.default_rng(7)   # separate stream so the sweep is self-contained
for n in sizes:
    mX, win, acr, tot, vX = ltv_estimates(n, conv_rng)
    rows.append({"n": n, "E[X]": mX, "within": win, "across": acr,
                 "within+across": tot, "Var(X)": vX})

conv = pd.DataFrame(rows).set_index("n")
pd.set_option("display.float_format", lambda v: f"{v:7.4f}")
print(conv)
print("\nanalytic targets: E[X]=0.300  within=1.710  across=0.400  total/Var(X)=2.110")

The wobble in the small-n rows is real sampling noise — with only 100 draws, a couple of unlucky stormy days can throw the across-regime term off badly. But by the time we reach hundreds of thousands of draws, every column has locked onto its target to three decimals. The plot below makes the funnel-shaped convergence obvious.

fig, axes = plt.subplots(1, 2, figsize=(11, 4))

# Left: E[X] converging to 0.300
axes[0].axhline(mean_X, ls="--", color="gray", label=f"analytic E[X] = {mean_X:.3f}")
axes[0].plot(conv.index, conv["E[X]"], "o-", color="#4C72B0")
axes[0].set_xscale("log")
axes[0].set_xlabel("number of draws (log scale)")
axes[0].set_ylabel("E[X] estimate")
axes[0].set_title("LIE: E[E[X|Y]] -> E[X]")
axes[0].legend()

# Right: the variance pieces converging to 1.710, 0.400, 2.110
axes[1].axhline(within, ls="--", color="#4C72B0")
axes[1].axhline(across, ls="--", color="#DD8452")
axes[1].axhline(var_X,  ls="--", color="#55A868")
axes[1].plot(conv.index, conv["within"], "o-", color="#4C72B0", label="within (-> 1.710)")
axes[1].plot(conv.index, conv["across"], "o-", color="#DD8452", label="across (-> 0.400)")
axes[1].plot(conv.index, conv["within+across"], "s-", color="#55A868",
             label="within+across (-> 2.110)")
axes[1].set_xscale("log")
axes[1].set_xlabel("number of draws (log scale)")
axes[1].set_ylabel("variance estimate")
axes[1].set_title("LTV: pieces -> analytic values")
axes[1].legend(fontsize=8)

fig.tight_layout()
fig.savefig("ltv_convergence.png", dpi=110)
print("saved ltv_convergence.png")

A final consistency check

Let’s assert, in code, what the eye already sees: at the largest sample size the simulated total variance and the simulated within-plus-across agree, and both sit near the analytic 2.110. We allow a small tolerance because this is a finite sample — the point of the whole exercise is that the gap shrinks, not that it is zero.

mX, win, acr, tot, vX = ltv_estimates(1_000_000, np.random.default_rng(123))

# Internal identity holds exactly on any single sample (it's algebra on that sample):
assert np.isclose(tot, vX, atol=1e-9), "within+across must equal pooled Var(X) on a sample"

# And both land near the analytic targets within Monte Carlo error at 1e6 draws:
assert abs(mX  - mean_X) < 0.01, f"E[X] off: {mX}"
assert abs(win - within) < 0.02, f"within off: {win}"
assert abs(acr - across) < 0.02, f"across off: {acr}"
assert abs(vX  - var_X)  < 0.02, f"Var(X) off: {vX}"

print("All checks pass.")
print(f"  E[X]            {mX:.4f}  (target {mean_X:.3f})")
print(f"  within          {win:.4f}  (target {within:.3f})")
print(f"  across          {acr:.4f}  (target {across:.3f})")
print(f"  within+across   {tot:.4f}  == Var(X) {vX:.4f}  (target {var_X:.3f})")

Your Turn

The §7 warning in the chapter was that the within/across split is only as honest as your choice of conditioning variable Y. Make it concrete by adding a third regime.

Sam now distinguishes calm, transition, and stormy days. Here is a 3\times 3 joint table on the same return grid (-2, 0, +2):

            X=-2    X=0    X=+2
calm        0.03    0.22    0.25     (row total 0.50)
transition  0.06    0.13    0.06     (row total 0.25)
stormy      0.11    0.10    0.04     (row total 0.25)

Before you run anything, predict: the new “transition” regime has its own conditional mean. Compared with the two-regime world, will the across-regime variance term go up or down? (Hint: revisit Check-Your-Understanding question 2 in the chapter — does a regime whose mean sits between the others add or remove spread among the conditional means?)

Then fill in the cell below to test your prediction. The scaffold draws from the new table and reuses the exact same grouping logic — only joint3 changes.

# --- Your Turn: three-regime version ---------------------------------------
your_rng = np.random.default_rng(2025)

joint3 = np.array([[0.03, 0.22, 0.25],   # calm
                   [0.06, 0.13, 0.06],   # transition
                   [0.11, 0.10, 0.04]])  # stormy
assert np.isclose(joint3.sum(), 1.0)

def simulate_from(joint_table, n, generator):
    flat_p = joint_table.ravel()
    regime_grid = np.repeat(np.arange(joint_table.shape[0]), joint_table.shape[1])
    ret_grid    = np.tile(x_vals, joint_table.shape[0])
    cells = generator.choice(flat_p.size, size=n, p=flat_p)
    df = pd.DataFrame({"regime": regime_grid[cells], "ret": ret_grid[cells]})
    gg = df.groupby("regime")["ret"]
    cm, cv = gg.mean(), gg.var(ddof=0)
    wt = gg.size() / len(df)
    mX = (wt * cm).sum()
    win = (wt * cv).sum()
    acr = (wt * (cm - mX)**2).sum()
    return cm, win, acr, win + acr, df["ret"].var(ddof=0)

cm3, win3, acr3, tot3, vX3 = simulate_from(joint3, 500_000, your_rng)

print("three-regime conditional means:")
print(cm3.rename({0: "calm", 1: "transition", 2: "stormy"}))
print()
print(f"within E[Var(X|Y)] = {win3:.3f}")
print(f"across Var(E[X|Y]) = {acr3:.3f}   (two-regime was 0.400 -- did it move as you predicted?)")
print(f"within + across    = {tot3:.3f}  ==  Var(X) = {vX3:.3f}")
print()
print("TODO: change the regime probabilities yourself and predict the split BEFORE running.")