| Time | Block |
|---|---|
| 0:00 – 0:20 | 74.1 The partially linear regression (PLR) model |
| 0:20 – 0:43 | 74.2 Frisch–Waugh–Lovell: partialling out, the linear way |
| 0:43 – 0:50 | Break |
| 0:50 – 1:13 | 74.3 Partialling out with machine learning |
| 1:13 – 1:40 | 74.4 Two biases of naive ML: regularisation and overfitting |
| 1:40 – 1:48 | Break |
| 1:48 – 2:08 | 74.5 Neyman orthogonality: why residualising both sides works |
| 2:08 – 2:30 | 74.6 Mini-project: a Monte Carlo of three estimators |
74.1 The Partially Linear Regression (PLR) Model (20 min)
The PLR model (Robinson, 1988; Chernozhukov et al., 2018) has two equations:
Y = θ·D + g₀(X) + ε, E[ε | D, X] = 0 D = m₀(X) + V, E[V | X] = 0
| Symbol | Our question | Role |
|---|---|---|
| Y | housing-starts growth over h months | outcome |
| D | change in the overnight rate (pp) | treatment |
| X | pre-decision economic situation (inflation, labour market, activity, global rates, past policy, past housing) | controls (confounders) |
| θ | effect per 1pp rate increase | target parameter: the causal effect under unconfoundedness |
| g₀(X) | how the economic situation affects housing, in any form | nuisance function |
| m₀(X) | the Bank's reaction function: expected rate change given the situation | nuisance function |
| V | rate change not explained by the situation: a policy surprise (Day 73.4) | the identifying variation |
A nuisance parameter is something we must estimate in order to get θ, but don't care about for its own sake. Here the nuisances are whole functions, potentially complex and non-linear. That's where machine learning comes in.
The model is partially linear: linear in the treatment (one number, θ, the effect per pp) but unrestricted in the controls. That's a sensible compromise for our question. We want one interpretable policy elasticity per horizon, and we want to make no functional-form assumptions about how a dozen correlated macro variables affect housing and policy.
What θ means: the average effect on housing growth of a 1pp larger rate change, comparing months with the same pre-decision economic situation. Under the Day 73 assumptions (unconfoundedness given X, overlap, no anticipation), that's a causal effect.
74.2 Frisch–Waugh–Lovell: Partialling Out, the Linear Way (23 min)
You know this theorem from econometrics. In the linear regression Y = θD + Xβ + ε, the OLS estimate of θ equals the coefficient from regressing:
- the residuals of Y on X, on
- the residuals of D on X.
"Partialling out" X from both sides isolates the variation in D that's unrelated to X, and asks how Y's unexplained variation co-moves with it.
import numpy as np
import statsmodels.api as sm
from causal.simulate import confounded_world
df = confounded_world() # Day 71: true effect -5
X = sm.add_constant(df[["cycle"]])
full = sm.OLS(df["housing_growth"], sm.add_constant(df[["rate_change", "cycle"]])).fit()
y_res = df["housing_growth"] - sm.OLS(df["housing_growth"], X).fit().fittedvalues
d_res = df["rate_change"] - sm.OLS(df["rate_change"], X).fit().fittedvalues
fwl = sm.OLS(y_res, d_res).fit()
print(f"full regression θ = {full.params['rate_change']:.4f}")
print(f"FWL partialling θ = {fwl.params.iloc[0]:.4f}") # identical
The two numbers are identical. FWL is the template for everything this week. DML is FWL where "regress on X" means any flexible ML model.
74.3 Partialling Out with Machine Learning (23 min)
Taking expectations of the PLR equations given X:
ℓ₀(X) ≡ E[Y | X] = θ·m₀(X) + g₀(X)
Subtracting gives the residual-on-residual equation:
Y − ℓ₀(X) = θ · (D − m₀(X)) + ε
So, just like FWL:
- estimate ℓ̂(X) = predict Y from X (any ML model);
- estimate m̂(X) = predict D from X (any ML model): the reaction function;
- regress the Y-residuals on the D-residuals:
θ̂ = Σ (Yᵢ − ℓ̂(Xᵢ))(Dᵢ − m̂(Xᵢ)) / Σ (Dᵢ − m̂(Xᵢ))²
This is Robinson's (1988) estimator. With linear models for ℓ and m it's exactly FWL; with random forests or XGBoost it handles non-linear confounding, which defeated linear controls in Day 71's exercise 3.
A non-linear world
Add to causal/simulate.py:
def nonlinear_world(T: int = 500, true_effect: float = -5.0, seed: int = 0) -> pd.DataFrame:
"""Confounding through non-linear functions of several controls (true effect known)."""
rng = np.random.default_rng(seed)
X = rng.normal(size=(T, 5))
m = 0.8 * np.tanh(2 * X[:, 0]) + 0.4 * X[:, 1] ** 2 - 0.3 * X[:, 2] * X[:, 3] # reaction fn
g = 6 * np.sin(1.5 * X[:, 0]) + 3 * X[:, 1] ** 2 + 2 * np.abs(X[:, 4]) # direct effects
d = m + rng.normal(0, 0.6, T)
y = true_effect * d + g + rng.normal(0, 2, T)
out = pd.DataFrame(X, columns=[f"x{i}" for i in range(5)])
out["rate_change"], out["housing_growth"] = d, y
return out
Compare linear partialling out with ML partialling out (for now without cross-fitting, which comes on Day 75):
from sklearn.ensemble import RandomForestRegressor
from sklearn.linear_model import LinearRegression
from causal.simulate import nonlinear_world
df = nonlinear_world(T=2000)
Xc = df[[f"x{i}" for i in range(5)]].to_numpy()
y, d = df["housing_growth"].to_numpy(), df["rate_change"].to_numpy()
def partial_out(learner_y, learner_d):
y_res = y - learner_y.fit(Xc, y).predict(Xc)
d_res = d - learner_d.fit(Xc, d).predict(Xc)
return (d_res @ y_res) / (d_res @ d_res)
print("linear controls:", round(partial_out(LinearRegression(), LinearRegression()), 2))
print("random forests (in-sample!):",
round(partial_out(RandomForestRegressor(200, min_samples_leaf=5, random_state=0),
RandomForestRegressor(200, min_samples_leaf=5, random_state=0)), 2))
The linear version is biased because the confounding is non-linear. The in-sample random-forest version is closer, but still biased. Why? Next block.
74.4 Two Biases of Naive ML: Regularisation and Overfitting (28 min)
Chernozhukov et al. (2018) identify two distinct problems with plugging ML into causal estimation.
Bias 1 — Regularisation bias (from a non-orthogonal score)
A tempting "naive" approach estimates g with ML and then regresses the remainder on D without residualising D:
θ̂_naive = Σ Dᵢ (Yᵢ − ĝ(Xᵢ)) / Σ Dᵢ²
ML estimators of g are regularised: they shrink, smooth and trade bias for variance, so ĝ has errors. Because D is correlated with X (that's confounding), those errors in ĝ leak directly into θ̂. The estimate converges too slowly and is centred in the wrong place.
Bias 2 — Overfitting bias (from using the same data twice)
Even with residualisation on both sides, if ℓ̂ and m̂ are fitted on the same observations whose residuals we then use, flexible learners partly fit the noise. In-sample residuals are too small and correlated with the noise in ways that bias θ̂. In 74.3, the random forests predicted the very points they were trained on.
| Estimator | Residualise D? | Out-of-sample nuisance predictions? | Bias |
|---|---|---|---|
| naive plug-in | ✗ | ✗ | regularisation and overfitting |
| orthogonal, in-sample | ✓ | ✗ | overfitting |
| DML (orthogonal + cross-fitting) | ✓ | ✓ | neither, asymptotically |
Here's the naive plug-in estimator. Fitting g requires an initial step; a common heuristic is to fit ML for Y on X and D jointly, then remove D's part. For teaching, we use an even simpler naive variant that captures the issue: fit ĝ on Y − θ₀D iteratively starting from θ₀ = 0:
def naive_plugin(learner, iters: int = 5):
theta = 0.0
for _ in range(iters):
g_hat = learner.fit(Xc, y - theta * d).predict(Xc)
theta = (d @ (y - g_hat)) / (d @ d) # no residualisation of D
return theta
print("naive plug-in:", round(naive_plugin(RandomForestRegressor(200, min_samples_leaf=5,
random_state=0)), 2))
Compare the three numbers (naive, orthogonal in-sample, and DML from Day 75) in the Monte Carlo below.
74.5 Neyman Orthogonality: Why Residualising Both Sides Works (20 min)
The DML estimator solves an equation of the form "average score = 0". For the partialling-out PLR, the score is:
ψ(W; θ, ℓ, m) = [Y − ℓ(X) − θ(D − m(X))] · (D − m(X))
Neyman orthogonality means: small mistakes in the nuisance functions ℓ and m have no first-order effect on the estimating equation. Mathematically, the derivative of the expected score with respect to the nuisances, at the truth, is zero.
Intuition
In the naive score, D · (Y − g(X)), an error δ in ĝ shifts the score by E[D · δ(X)]. Since D depends on X, that's not zero: errors leak through. In the orthogonal score, the error is multiplied by (D − m(X)) = V, the policy surprise, which by construction is uncorrelated with any function of X. So errors in ℓ have no first-order effect. Errors in m enter only multiplied by errors in ℓ, a product of two small numbers.
Why that matters for ML
ML nuisance estimates converge slower than parametric ones (they're flexible). Orthogonality means θ̂ only needs the product of the two nuisance errors to be small: each converging faster than n^(−1/4) is enough, rather than the usual n^(−1/2). That's a mild requirement that random forests, boosting and lasso can meet in many settings. Combined with cross-fitting (Day 75), θ̂ is approximately normal, centred on the truth, with valid standard errors.
One sentence to remember: DML predicts what the economic situation implies for both the outcome and the policy, subtracts it from both, and estimates the effect from what's left, using predictions made on data the models didn't train on.
74.6 Mini-Project: A Monte Carlo of Three Estimators (23 min)
Create causal/montecarlo.py:
"""Monte Carlo comparisons of causal estimators (bias, spread, coverage)."""
import numpy as np
import pandas as pd
from sklearn.base import clone
from sklearn.ensemble import RandomForestRegressor
from sklearn.model_selection import KFold
from causal.simulate import nonlinear_world
RF = RandomForestRegressor(n_estimators=100, min_samples_leaf=5, random_state=0, n_jobs=2)
def naive(y, d, X, learner, iters=5):
theta = 0.0
for _ in range(iters):
g = clone(learner).fit(X, y - theta * d).predict(X)
theta = (d @ (y - g)) / (d @ d)
return theta
def orthogonal_insample(y, d, X, learner):
yr = y - clone(learner).fit(X, y).predict(X)
dr = d - clone(learner).fit(X, d).predict(X)
return (dr @ yr) / (dr @ dr)
def dml(y, d, X, learner, n_folds=5, seed=0):
yr, dr = np.empty_like(y), np.empty_like(d)
for tr, te in KFold(n_folds, shuffle=True, random_state=seed).split(X):
yr[te] = y[te] - clone(learner).fit(X[tr], y[tr]).predict(X[te])
dr[te] = d[te] - clone(learner).fit(X[tr], d[tr]).predict(X[te])
return (dr @ yr) / (dr @ dr)
def run(n_sims=100, T=500, true_effect=-5.0):
rows = []
for s in range(n_sims):
df = nonlinear_world(T=T, true_effect=true_effect, seed=s)
X = df[[f"x{i}" for i in range(5)]].to_numpy()
y, d = df["housing_growth"].to_numpy(), df["rate_change"].to_numpy()
rows += [{"sim": s, "estimator": "naive plug-in", "theta": naive(y, d, X, RF)},
{"sim": s, "estimator": "orthogonal, in-sample", "theta": orthogonal_insample(y, d, X, RF)},
{"sim": s, "estimator": "DML (cross-fitted)", "theta": dml(y, d, X, RF, seed=s)}]
res = pd.DataFrame(rows)
summary = res.groupby("estimator")["theta"].agg(mean="mean", sd="std")
summary["bias"] = summary["mean"] - true_effect
return res, summary
if __name__ == "__main__":
res, summary = run()
res.to_parquet("data/causal/mc_three_estimators.parquet", index=False)
print(summary.round(3))
python -m causal.montecarlo # a few minutes: 100 simulations × 3 estimators
Plot the three distributions of θ̂ (histograms with a vertical line at −5) in notebooks/day74.ipynb. You should see the naive plug-in shifted, the in-sample orthogonal estimator less shifted, and DML centred on −5. This reproduces, on your own simulation, the key figure in Chernozhukov et al. (2018).
Day 74 checklist
- I can write the PLR model, identify θ, the nuisance functions ℓ₀, m₀, g₀, and the identifying variation V
- I can show that FWL partialling-out equals the full linear regression
- I can implement Robinson's residual-on-residual estimator with any ML learner
- I can explain regularisation bias and overfitting bias, and which fix addresses which
- I can explain Neyman orthogonality intuitively, and why it lets slower ML nuisance estimates work
- My Monte Carlo shows naive and in-sample estimators biased, and cross-fitted DML centred on the truth
git add causal notebooks && git commit -m "Day 74: PLR, nuisance functions, ML biases, orthogonality, Monte Carlo" && git push