Chapter 4: Statistical Predictive Models
Chapter Introduction
Every lending decision, every insurance quote, every marketing budget and every portfolio weight is, at bottom, a prediction that somebody has agreed to act on. A bank that lends $10 000 to a stranger is betting on three numbers it cannot observe: whether the borrower will default, how much will still be owed if they do, and how much of that can be recovered afterwards. The bank’s entire margin lives in the gap between the interest it charges and the product of those three numbers. If it estimates them well, it earns a steady spread; if it estimates them badly, the defaults eat the spread and the good borrowers take their business to a competitor who priced them correctly. This chapter is about building those numbers from data — and about the discipline that separates a prediction you can defend from one that merely fits.
The vehicle is a single real case that runs through the first three sections: 8 000 loans from LendingClub, the largest peer-to-peer lender in the United States, every one of which has finished its life, so that we know whether it was repaid. From the raw ledger columns you will define default, exposure at default (EAD) and loss given default (LGD); convert strings, dates and missing values into 25 usable regressors; fit multiple regressions; discover how little of LGD any origination feature explains; and learn to judge a model on data it never saw. Section 4.2 turns to the question every applied statistician faces after the first fit: which variables, in which form? You will see why R² is a dishonest scorekeeper, why two variables can be individually useless and jointly decisive, how an exhaustive best-subset search behaves, and how residual plots, t-tests and variance inflation factors tell you whether a coefficient means anything.
Section 4.3 moves beyond the line. A regression tree finds thresholds and interactions on its own — and memorises noise just as readily. You will watch a tree choose its first split by hand, overfit to a training error of exactly zero, and then be tamed by two ideas that dominate modern tabular prediction: averaging many deep trees (random forests) and adding many shallow trees fitted to residuals (gradient boosting). The same boosting machine is then pointed at ten years of the S&P 500 under a walk-forward protocol that refuses to peek into the future, with a result that is instructive precisely because it is so modest. The section closes with the third credit-risk factor — the probability of default — by logistic regression, judged by the ROC curve, and turned into a decision by a cutoff that reflects the lender’s gain and loss rather than the arbitrary number 0.5.
Section 4.4 asks a question that none of the earlier machinery can answer: not “given \(x\), what is \(y\)?” but “if I change \(x\), what happens to \(y\)?”. You will build a case in which the single best predictor of repayment is a completely useless policy lever, understand why, and meet three designs — regression adjustment, difference-in-differences, instrumental variables — that recover a causal effect when a naive regression cannot. Section 4.5 returns to prediction with a modern twist: the stock cross-section, standardised month by month, and a single attention layer written in twelve lines of numpy that replaces the index provider’s sector label with a peer group learned from the data. The layer is run through the same walk-forward protocol, and the honest result — that on 213 surviving large caps over 71 months nothing beats zero — is the point, not an embarrassment.
By the end of the chapter you will be able to take a raw business table to a defensible predictive model, report its out-of-sample performance in numbers a manager can act on, know which coefficients deserve to be called real, choose between linear and tree-based learners on evidence, and recognise the moment when the question has quietly changed from prediction to causation. Read the chapter alongside the Chapter 4 slide deck: sections, datasets, seeds and printed numbers are identical, so the two never disagree.
Table of Contents
- Credit Risk Data: Default, EAD, LGD, Feature Engineering, Multiple Regression
- Model Selection: Adjusted R², Interactions, Best Subset, Inference, Multicollinearity
- Nonlinear Models: Trees, Random Forests, Gradient Boosting, Logistic Classification
- Causal Analysis: Prediction vs Causation, Confounders, Difference-in-Differences, Instrumental Variables
- Cross-Sectional Attention Features: Transformers for the Stock Cross-Section
Credit Risk Data: Default, EAD, LGD, Feature Engineering, Multiple Regression
How a lending platform makes money
LendingClub is a marketplace. Borrowers arrive wanting money for debt consolidation, a large purchase or a medical bill; investors arrive wanting a return above what a savings account pays. The platform sits between them and does the one thing neither side can do alone: it underwrites. Each application is scored by a credit-risk model, assigned a grade from A to G, and offered a rate and a term (36 or 60 months) that reflect that grade. The platform then either keeps the loan on its own book or sells it to investors as a note. Its revenue is an origination fee charged to the borrower, a servicing fee charged to the investor, and — on loans it retains — the net interest margin.
Every one of those revenue lines depends on the model being right. Charge too little for the risk and the defaults eat the margin; charge too much and the creditworthy borrowers, who have alternatives, go elsewhere, leaving the platform with exactly the applicants it should have declined. The credit-risk model is not a back-office nicety; it is the business, as the pricing engine is at an insurer.
The quantity the model is ultimately trying to price is expected loss, and the industry decomposes it into three factors:
\[\text{EL} = \underbrace{\text{PD}}_{\text{probability of default}} \times \underbrace{\text{EAD}}_{\text{exposure at default}} \times \underbrace{\text{LGD}}_{\text{loss given default}}.\]
PD is the probability that the loan will not be repaid. EAD is the fraction of the principal still outstanding at the moment of default — a loan that fails in month 3 of 60 exposes almost everything; one that fails in month 55 exposes very little. LGD is the fraction of that exposure that is not recovered afterwards through collections, settlements or the sale of the debt. The decomposition matters for modelling because the three factors are different kinds of statistical object. PD is the probability of a 0/1 event, which makes it a classification problem; we will fit it in Section 4.3 with logistic regression. EAD and LGD are continuous fractions that exist only for loans that actually defaulted, which makes them regression problems on that subset. This section builds the two regressions; the classification waits.
Loading 8 000 finished loans
The file lending_club_sample.csv is a random sample of 8 000 rows from the roughly 150 000-loan dataset used in the original lecture notebook. Every loan in it has reached the end of its term, so its final status is known: either Fully Paid or Charged Off, the latter being the lender’s term for a loan that has been written off as uncollectable. Three ledger columns carry the story we need: funded_amnt is what was lent, total_rec_prncp is the principal repaid over the life of the loan, and recoveries is the cash clawed back after the charge-off.
The table has 8 000 rows and 29 columns, and the split is exactly 80 % Fully Paid, 20 % Charged Off — the sample was drawn to preserve the population default rate. The map call turns the status string into the 0/1 column default; the three lines after it are pure ledger arithmetic. EADamount is the dollar exposure, principal lent minus principal repaid; EAD is the same thing as a fraction of the original loan; and LGD is one minus the recovered share of that exposure. The 1e-4 in the denominator is a guard: one charged-off loan in the file has EADamount exactly 0, and without the guard its LGD would be a division by zero.
Look at the four printed loans. The first was funded at $10 000 and repaid $1 437 of principal before failing, so 85.6 % was exposed; collections recovered $1 240, leaving an LGD of 0.855. The third is the interesting one: $5 200 lent, $2 080 repaid, $3 120 exposed — and then $3 170 recovered, more than was outstanding, so its LGD is −0.016. Collectors also collect late fees and accrued interest; six loans show this pattern, which is why the minimum LGD is −0.177 rather than 0. We keep the raw numbers, because real ledgers do this and a model should see the data it will meet in production.
Among the 1 604 charged-off loans the mean EAD is 0.698: on average, 70 % of the principal was still outstanding when the borrower stopped paying. The mean LGD is 0.892: recoveries claw back barely a tenth of what is owed. Those two numbers already tell you most of what a lender needs to know about a default — it happens early, and it is nearly a total loss. The regressions that follow ask whether anything known at origination moves either number.
Feature engineering: term, employment, credit age
Raw columns are not regressors. term reads " 36 months" as a string; emp_length reads "10+ years", "3 years", "< 1 year" or is missing; earliest_cr_line reads "Jan-1986"; the FICO score is a range. Feature engineering is the sequence of judgements that turn these into numbers whose magnitude means something. Each judgement is small, and the ones you get wrong are silent — the regression will happily fit a column that says “under one year” is the same as “one year”.
str.extract(r"(\d+)") pulls the first run of digits out of each string, so " 36 months" becomes 36 and "10+ years" becomes 10. The term split is 6 047 loans at 36 months and 1 953 at 60. The employment column needs two corrections that the regular expression cannot make. First, "< 1 year" extracts as 1, which is the wrong end of the ordering — under a year should sit below one year, not equal to it — so .where(...) overrides it to 0. Second, 484 loans have no employment record at all, and fillna(0) puts them at zero years; but “we do not know” is not the same information as “under a year”, so the flag employed records which is which. With both columns in the model, a missing record can carry its own coefficient instead of being silently merged with the least-experienced applicants. The value counts show the resulting shape: 1 127 loans at zero, a fairly flat run from one to nine years, and 2 643 at the 10+ ceiling.
The buggy column has only 484 zeros — the missing records — because every "< 1 year" borrower has been extracted as 1.0. The regular expression (\d+) finds the digit in "< 1 year" and returns it, and nothing in the string tells str.extract that the < reverses its meaning. Six hundred and forty-three borrowers who have been in their job for a few months are now indistinguishable from those with a full year, and the ordering of the feature is wrong at exactly the end where default risk is highest. The fix is the .where(debt["emp_length"] != "< 1 year", 0) clause from the cell above: check the string, not the extracted digit. The general lesson is that text-to-number conversions need a value count afterwards, every time.
Four more judgements are packed into this cell. Credit age. pd.to_datetime(..., format="%b-%Y") parses Jan-1986 directly — no guessing, no warnings — and subtracting from the dataset’s end date (29 December 2018) gives the credit history in days: 12 050 for Jan-1986, and across the file from 1 277 days (three and a half years) to 20 757 (57 years), median 6 724. FICO. The bureau reports a range; the midpoint is the natural point estimate. Instalment burden. installment / annual_inc is the monthly payment as a share of annual income — a ratio of a few thousandths, comparable across a $30 000 earner and a $300 000 earner. Three applicants report zero income; replacing that zero with NaN before dividing keeps their ratio missing rather than infinite, and dropna() later removes them. Instalment balance. total_bal_il / annual_inc is the same idea for other instalment debt; it is NaN when the borrower has none, and there the honest fill really is 0.
The two “time since trouble” clocks deserve their own paragraph, because the missing-value decision is the opposite of the usual one. mths_since_last_delinq is missing for roughly half the loans. The reflex from Chapter 2 is to fill with 0 or the mean — but ask what the missing value means: the bureau has no delinquency on record, so the borrower has never had one. Filling with 0 would say the delinquency happened this month, the worst possible reading of the cleanest possible history; the mean would put the never-delinquent borrower among the delinquent ones. The value that preserves the ordering is a large number: 360 months, comfortably beyond the observed maximum of 170. The printed medians (360 for both clocks) show that more than half the sample now sits at “never”.
Dummies, assembly and the split
Three categorical columns remain: grade (A–G), home_ownership and verification_status. A regression needs a 0/1 indicator for each level, with one level left out. The reason is arithmetic rather than convention: if you include an indicator for every level and an intercept, the indicators sum to the column of ones, \(X^\top X\) is singular, and the coefficients are not identified. Dropping one level makes it the baseline against which the others are measured.
get_dummies(...).iloc[:, :-1] drops the last level (grade G, verif_Verified). Home ownership is handled differently, and deliberately so: the column also contains ANY and OTHER with one loan each, and the original notebook keeps the three big categories OWN, RENT and MORTGAGE explicitly. That looks harmless — two loans out of 8 000 are excluded rather than one level — but it means the three retained dummies sum to one on 7 998 of 8 000 rows, which is almost the dummy trap. Section 4.2 will show you what “almost” costs when we compute variance inflation factors.
The assembled frame has 25 features and three targets. dropna() removes six rows with a missing dti or revol_util, the three zero-income applicants among them. Then a step that is easy to skip and costly to omit: the rows are permuted with a seeded generator before the first 70 % are taken as training data. Cross-sectional rows have no natural order, but files almost always do — by issue date, grade or loan ID — and the first 70 % of a sorted file would train on one kind of loan and test on another. The seed 5650 makes the permutation reproducible, so every number here and in the slides is the same number. The result is 5 595 training and 2 399 test loans; of those, 1 146 and 458 defaulted, and only these rows enter the EAD and LGD regressions.
Multiple regression for EAD
The model is the one from your earlier statistics course, with 25 regressors instead of one:
\[\text{EAD}_i = \beta_0 + \beta_1 x_{i1} + \cdots + \beta_{25} x_{i,25} + \varepsilon_i, \qquad \hat\beta = (X^\top X)^{-1} X^\top y.\]
statsmodels needs the intercept column added explicitly with sm.add_constant; it does not assume one. The fitted object carries .params, .tvalues, .rsquared and .rsquared_adj, and we print the seven coefficients that matter to the story rather than the full summary table.
R² is 0.192: the 25 origination features explain 19 % of the variation in how much principal is still outstanding at default. Adjusted R² is 0.174, the same number after charging for the 25 columns. Read the coefficients with their t-statistics. term has the largest, t = 7.5, with a coefficient of 0.0047 per month: a 60-month loan has, other things equal, 24 × 0.0047 ≈ 0.11 more of its principal outstanding at default than a 36-month loan, which makes sense — longer loans amortise more slowly, so any given month of failure leaves more unpaid. int_rate (t = 3.5) pushes the same way: a higher rate means a larger share of each instalment goes to interest rather than principal. install, the burden ratio, has a negative coefficient (t = −2.7): a borrower carrying a heavier instalment relative to income has paid more principal down before failing. balance (t = 4.9) is positive — borrowers with more other instalment debt fail earlier. And fico (t = 0.75) and dti (t = 0.23) add nothing once the rest are in: their coefficients are indistinguishable from zero. That last observation is not “FICO does not matter for credit”; it is “FICO does not move the timing of default, given the other 23 columns”. Section 4.2 will make the distinction precise.
The LGD model and what an R² of 0.04 means
Now the same 25 features against LGD. Before running the cell, form an expectation. LGD is determined by what collectors recover after the charge-off — the legal process in the borrower’s state, whether the borrower has assets, whether the debt is sold and at what price, luck. Almost none of that is visible in the origination file.
The EAD regression reached R² = 0.192 with these 25 features. Will LGD land near 0.2 as well, near 0.5, or under 0.05?
R² = 0.039, adjusted R² = 0.018, and exactly one of the 25 coefficients has |t| above 2. The full 150 000-loan notebook found 0.012. For practical purposes LGD is a constant near 0.89, and a model that knows it is a constant is more useful than one that pretends otherwise: the honest forecast for every defaulted loan is “you will recover about 11 cents on the dollar”, and the honest uncertainty is the standard deviation of LGD, not a regression’s standard error. This is the first of several places in the chapter where the right answer to “can we predict this?” is “no, and here is the number that proves it”.
Judging on data the model never saw
R² on the training sample is the fraction of variance the model explained on the rows it was fitted to. It is not a forecast of how the model will do on the next 458 loans, because the fit has already used those training rows to choose 26 numbers. The only score that answers the deployment question is computed on held-out data, and it should be computed with the same formulas so the two are comparable:
\[R^2 = 1 - \frac{\sum (y_i - \hat y_i)^2}{\sum (y_i - \bar y)^2}, \qquad R^2_{\text{adj}} = 1 - (1 - R^2)\,\frac{n-1}{n-k-1}, \qquad \text{RMSE} = \sqrt{\tfrac{1}{n}\sum (y_i - \hat y_i)^2}.\]
Two details in the helper matter. has_constant="add" forces add_constant to add the intercept even when it thinks one is present — a quiet source of shape errors otherwise. And the adjusted R² uses the test sample’s own \(n\) (458), so the penalty per column is heavier than on the 1 146 training rows, as it should be.
EAD: training R² 0.192 becomes 0.156 on test, and adjusted R² falls from 0.174 to 0.107. The RMSE is 0.204 in-sample and 0.190 out-of-sample — a typical exposure forecast is off by 19 percentage points of principal, which, given that EAD itself has a standard deviation of about 0.23, is a modest improvement over guessing the mean. LGD: test adjusted R² is −0.080. A negative out-of-sample R² has a precise meaning: the model’s predictions are further from the truth than the training mean would have been. When you see it, the right model is the constant.
The last two lines put the three factors together. The training default rate is 0.205; the mean EAD and LGD among defaulted loans are 0.698 and 0.892; and their product is an expected loss of 12.7 cents per dollar lent. That number is the reason the average interest rate in this file is about 13 %: a lender who expects to lose 12.7 cents needs to charge at least that much just to break even before funding costs, and the spread above it is the margin. Every decision in the rest of the chapter is, in one way or another, an attempt to make that 12.7 smaller by declining the right loans.
Report the test-set adjusted R² and the test-set RMSE, in the units of the target. Training R² is a diagnostic, not a result. If test adjusted R² is negative, say so and recommend the mean — a manager will trust an analyst who reports “unpredictable” far more than one who reports a model that loses to a constant.
Zero — all principal came back, so funded minus repaid is 0. EAD and LGD are only meaningful for loans that actually defaulted (1 604 of 8 000 here; 1 146 in the training split), so the regressions use that subset. PD, the third factor, is fitted on all loans as a classification.
A missing value means the bureau has no delinquency on record — the borrower never had one. Filling with 0 says the delinquency happened this month (the opposite of the truth); the mean puts a clean history in the middle of the dirty ones. A value beyond the observed maximum (170 months) keeps the ordering right: “never” sits past “long ago”.
Model Selection: Adjusted R², Interactions, Best Subset, Inference, Multicollinearity
Rebuilding the features
Every section of this chapter starts by re-importing and rebuilding what it needs, so that a reader who opens the book at Section 4.2 can run it without scrolling up. The cell below repeats the feature engineering of Section 4.1 in compressed form — same conversions, same fills, same seed, same 70/30 split — and refits the 25-feature EAD model.
Same 1 146 and 458 defaulted rows, same adjusted R² of 0.174. If either number differs on your machine, a line of the feature engineering has drifted, and everything downstream will drift with it.
Why R² cannot be the scorekeeper
Ordinary least squares chooses \(\hat\beta\) to minimise the sum of squared residuals. Add any column to \(X\) — a real predictor, a column of random noise, the borrower’s shoe size — and the minimisation now has one more degree of freedom, so the minimum can only fall or stay put. R² therefore never decreases when a variable is added, whatever the variable is. This is not a statistical subtlety; it is a property of the minimisation, and it means R² is useless for comparing models with different numbers of columns: the bigger model always wins.
Adjusted R² charges an entrance fee. The formula \(1 - (1 - R^2)\,(n-1)/(n-k-1)\) multiplies the unexplained fraction by a factor that grows with \(k\), so a new column must reduce the residual sum of squares by more than the fee to raise the score. The fee is small when \(n\) is large relative to \(k\) — adjusted R² is a weak penalty compared with AIC or BIC — but on 1 146 rows and 25 columns it is not negligible. The cleanest way to see it work is to add columns that cannot possibly matter.
We append ten columns of pure Gaussian noise (seed 1) to the 25 features and refit. Will R² rise or fall? Will adjusted R² rise or fall?
R² rises from 0.1920 to 0.1959; adjusted R² falls from 0.1740 to 0.1706. Ten columns of random numbers “explained” an extra 0.4 % of the variance of EAD, because with 1 146 rows there is always some accidental correlation to be found, and OLS finds it. The adjustment charged 10 × (roughly 0.07 %) for the privilege and correctly reported that the model got worse. Never select a model on raw R². The same logic applies, with larger stakes, to any procedure that searches over many candidate columns: the more you search, the more accidental fit you will find, and the more the score you report needs to be computed on data the search never touched.
Expanding the pool with interactions
We take eight base features and form all \(\binom{8}{2} = 28\) pairwise products, giving 36 candidates. Squares are omitted deliberately: term takes only the values 36 and 60, so term**2 is a linear function of term and would duplicate it. The full notebook does this for all 25 features plus a set of ratios and ends up with 650 columns — a number to keep in mind when we count models. The cell also defines a scoring helper that fits on training rows and returns the training and test adjusted R², and uses it to compare the 36-column model with the plain 8-column one.
The 36-column model scores 0.195 on training data and 0.079 on test; the 8-column model scores 0.178 and 0.141. The interactions bought 1.7 points of in-sample adjusted R² and cost 6 points out of sample. This is the pattern you should expect whenever the number of columns is large relative to the signal: 36 columns on 1 146 rows is plenty of room to fit noise, and the adjustment’s entrance fee, while enough to catch ten pure-noise columns, is not enough to catch 28 columns that each carry a sliver of accidental correlation. The test set catches all of them.
Best subset: 255 models, one envelope
Best-subset selection is the exhaustive answer to “which columns?”: fit every non-empty subset of the \(k\) candidates, and keep the best by your criterion. For \(k\) candidates there are \(2^k - 1\) subsets — 255 for \(k = 8\), about \(6.9 \times 10^{10}\) for \(k = 36\), and \(2^{650}\) for the notebook’s full pool, a number with 196 digits. So exhaustive search over the whole pool is impossible, and the practical procedure has two stages: screen the pool down to a handful of candidates, then search the handful exhaustively. The screen here keeps the eight columns with the highest absolute correlation with EAD (the notebook keeps the top 15 of 650, giving \(2^{15} = 32\,768\) fits).
The screen’s top eight are dominated by variants of two variables: int_rate*term (0.379), int_rate*fico, term, int_rate, term*fico, then three more products of int_rate or term with something else. install and balance — which Section 4.1 showed to have t-statistics of −2.7 and 4.9 in the full model — do not make the cut, because their marginal correlation with EAD is small. Keep that in mind.
The search then fits all 255 subsets and records, for each size \(k\), the best adjusted R². The envelope climbs 0.143 → 0.152 → 0.158 → 0.161 at \(k = 4\) and then falls — 0.161, 0.160, 0.159, 0.158 — as a fifth, sixth, seventh and eighth candidate are added. Beyond four columns every additional candidate is charged more than it earns. The winner is int_rate, term*fico, int_rate*creditdays and term*revol_util: the rate, and three interactions of term or rate with something else.
Every grey dot is one of the 255 fits; the red line is the best of each size. Two things are visible that a table hides. The spread of the dots at each \(k\) is wide — at \(k = 3\) the worst subset scores about 0.11 and the best 0.158 — so which three columns matters far more than how many. And the envelope is flat from \(k = 4\) onwards: the search is not telling you that four is a magic number, it is telling you that after four the candidates it has left are redundant with the ones already in.
The chosen model on test data — and a surprise
The best-four model scores 0.130 on test — far above the 36-column model’s 0.079, with a tenth of the columns. Best subset did what it was supposed to do: found a compact model that generalises better than the kitchen sink.
But the plain eight base features, with no interactions and no search, score 0.141. The searched model loses to the unsearched one. Why? Look back at the screen. Correlation with EAD kept five variants of int_rate and term — which are largely the same information five times over — and dropped install and balance, which are weak marginally but useful jointly. That is the x1*x2 lesson in reverse: a one-at-a-time screen cannot see joint value, so it discarded the two columns that would have helped and kept redundant copies of the two it liked. The best-subset search then did an excellent job of choosing among the wrong candidates.
Screening is itself a modelling choice with consequences. Screening by marginal correlation favours redundant copies of strong variables; screening by t-statistic in the full model would have kept install and balance, and so would a forward-stepwise search that adds the column which most improves adjusted R² given what is already in. The honest report is the one printed above: all three models, training and test, side by side.
When the pool is too large to search exhaustively, three greedy procedures replace the search: forward selection starts empty and adds the column that most improves the criterion; backward elimination starts full and removes the least useful; stepwise alternates the two. None is guaranteed to find the best subset, but forward selection evaluates \(k(k+1)/2\) models instead of \(2^k\) and judges each candidate given the columns already chosen. The isom5650.regression package used in Colab implements all three.
Inference: residuals, t-tests, VIF
Everything so far asked “how close?” — a prediction question, answered by test-set R² and RMSE, which need no assumption about the errors at all. Inference asks a different question: “is this coefficient real?” — is \(\beta_j\) distinguishable from zero, and what range of values is consistent with the data? That question is answered by t-statistics, p-values and confidence intervals, and every one of them rests on assumptions about \(\varepsilon\): independent errors, constant variance, and either normality or a large enough \(n\) for the central limit theorem to stand in. Violate the assumptions and the numbers still print — statsmodels does not know — but their coverage is wrong: a “95 %” interval may cover 80 % of the time, or 99 %.
So before trusting a t-test, look at the residuals.
The plot has two hard diagonal edges. EAD lives in \([0, 1]\), so the residual \(y - \hat y\) can never exceed \(1 - \hat y\) or fall below \(-\hat y\); the cloud is sliced off by those two lines. The spread is visibly not constant — wider in the middle of the fitted range, pinched at the ends — and each vertical slice is skewed, because so many defaults happen early (EAD near 1) and pile up against the upper edge. Homoskedasticity and normality both fail, mildly. The model’s predictions and test RMSE are exactly as good as they were; what changes is that the t-statistics below are approximate. With 1 146 rows the approximation is not bad: read large t-statistics as evidence of direction and borderline ones as borderline.
R², adjusted R², RMSE and the predictions themselves are computed from \(y\) and \(\hat y\) alone — no distributional assumption enters. t-tests, p-values and confidence intervals need independent, constant-variance errors (and normality, or a large \(n\)). Validate the assumptions when you are about to call a coefficient significant; skip the validation when you only need a forecast. A residual plot is thirty seconds well spent in the first case and irrelevant in the second.
The t-statistic of a coefficient is \(t = \hat\beta_j / \text{se}(\hat\beta_j)\): the estimate in units of its own uncertainty. If \(|t| < 2\) the 95 % interval covers zero, and the data cannot distinguish the coefficient from zero given the other columns in the model. That last clause is the one students forget. A small \(|t|\) does not say the variable is unrelated to \(y\); it says the variable is redundant in this model — its information is already carried by other columns, or it carries none. Remove the redundant columns and the survivors’ t-statistics generally sharpen, because the standard errors were inflated by the redundancy.
19 of 25 coefficients have \(|t| < 2\). At the bottom, grade_F (0.01), employed (0.06) and grade_E (0.07) are as close to zero as a coefficient can be. At the top, term (7.53), balance (4.86) and int_rate (3.49) carry the model — the same three variables the best-subset search kept reaching for. Refitting on the six survivors (int_rate, term, revol_util, install, balance, creditdays) gives adjusted R² 0.178, above the 25-feature model’s 0.174. Nineteen columns were pure cost: they added nothing to the fit and inflated every standard error in the table.
Multicollinearity and the variance inflation factor
Redundancy has a name when it is severe: multicollinearity, one predictor being nearly a linear combination of the others. It does not bias OLS and barely hurts prediction, but it wrecks inference: the data cannot tell which of two near-identical columns deserves the coefficient, so both get a huge standard error and often opposite signs. The diagnostic is the variance inflation factor,
\[\text{VIF}_j = \frac{1}{1 - R^2_j},\]
where \(R^2_j\) is the R² of regressing \(x_j\) on all the other predictors. A VIF of 1 means \(x_j\) is orthogonal to the rest; a VIF of 10 means 90 % of its variance is already explained by the others and its standard error is \(\sqrt{10} \approx 3.2\) times what it would be in an orthogonal design. The conventional alarm threshold is 10. We warned in Section 4.1 that keeping all three home-ownership dummies was a trap; the VIF makes the cost visible, and the cell then follows the notebook’s remedy — drop the three home dummies and grades A–D — and refits.
home_RENT has a VIF of 290.6, home_MORTGAGE 286.2, home_OWN 97.0. The three dummies sum to one on all but two loans, so regressing any one of them on the other two plus the intercept gives an R² of about 0.997, and \(1/(1 - 0.997)\) explodes. This is the dummy trap in its “almost” form: not an exact singularity that statsmodels would refuse, but a near-singularity that it fits silently with standard errors seventeen times too large. The grade dummies are the second story: B, C and D sit at 19–31 because LendingClub sets int_rate from the grade, so once the rate is in the model the grade dummies are close to a function of it.
After dropping the seven columns the largest VIF is 10.8, shared by dti and install, which both have annual income in the denominator — a mild, explicable overlap that one might leave alone or resolve by keeping just one. The model has 18 features instead of 25, and its test adjusted R² rises from 0.107 to 0.121. That is the empirical version of the claim above: collinearity hurt inference badly (VIFs in the hundreds) and prediction a little (1.4 points of test R² recovered by removing it).
Adjusted R² and the test score answer does this model predict? The t-test answers is this coefficient distinguishable from zero given the others? The VIF answers is this column nearly a copy of the others? A variable can fail the t-test because it is useless or because its twin is also in the model; the VIF tells you which. Run all three before you write a sentence about what “drives” EAD.
R² never falls when a column is added because OLS can always use the extra degree of freedom to fit a little more of the training noise. Adjusted R² multiplies the unexplained fraction by \((n-1)/(n-k-1)\), so each column must earn more than its fee; ten useless columns paid the fee and earned nothing, so the score fell.
VIF \(= 1/(1 - R^2_j)\), so 291 means the other predictors explain \(R^2_j \approx 0.9966\) of that dummy’s variance and its standard error is \(\sqrt{291} \approx 17\) times larger than in an orthogonal design. Cause: OWN + RENT + MORTGAGE = 1 on all but two loans, so the three are almost a copy of the intercept — the dummy trap. Drop one (or all three) and the VIFs collapse.
Nonlinear Models: Trees, Random Forests, Gradient Boosting, Logistic Classification
Why trees?
A linear model adds up effects. Each column contributes \(\beta_j x_j\) regardless of what the other columns are doing, and the only way to express “the effect of the rate depends on the term” is to build int_rate * term by hand, as Section 4.2 did. Threshold effects are worse: “IF the term is 60 months AND the other instalment balance exceeds 1 % of income THEN exposure is high” is a rule a linear model can only approximate with hand-made indicators, and real business relationships are full of such rules — a limit that kicks in above a score, a risk that appears only when two conditions hold at once.
A regression tree builds such rules from the data. It partitions the feature space into boxes by a sequence of binary splits — first on one feature at one threshold, then, within each side, on another feature at another threshold — and predicts the mean of \(y\) in each box. The result is a piecewise-constant approximation with thresholds and interactions built in: two successive splits on different features are an interaction, and no product term was needed. The price is that a tree has no notion of smoothness or of a linear trend, so it approximates a straight line with a staircase, and — as we will see — it can keep splitting until every box holds one observation.
We rebuild a compact version of the credit data with eight numeric features, the same seed and the same split as before, and start with the mechanics of one split.
The split rule is simple to state and expensive to run: for every feature and every cut-point between two consecutive observed values, divide the rows into left and right, compute the sum of squared errors around each child’s mean, and keep the cut that most reduces \(\text{SSE}_{\text{left}} + \text{SSE}_{\text{right}}\) below the parent’s SSE. The loop above does this for one feature. The parent SSE of the 1 146 training loans is 59.11 (the standard deviation of EAD is about 0.227, and \(1146 \times 0.227^2 \approx 59\)). The best cut on int_rate is at 17.93 %, which drops the children’s SSE to 54.44 — a reduction of 4.67, about 8 % of the total. A tree does this for all eight features, takes the best, and then repeats inside each child until a stopping rule fires.
Two levels, then no limit
export_text prints the fitted rules. The root does not split on int_rate — term at 48 months (the midpoint between 36 and 60) reduces SSE more — and then each side splits on balance. Four leaves, four means: 36-month loans with little other instalment debt have EAD 0.572; 60-month loans with a lot of it, 0.860. The tree found the term × balance interaction that Section 4.2 had to construct by hand and then failed to keep through its correlation screen, and it did so with four numbers. Its test RMSE is 0.191, already equal to the 25-feature OLS model’s 0.190.
Now the failure mode. DecisionTreeRegressor() with default settings has no depth limit and a minimum leaf size of 1.
With no depth limit the tree keeps splitting until it cannot. On 1 146 training rows, what will the training RMSE be, and how many leaves?
1 145 leaves for 1 146 loans (two loans with identical features share one), training RMSE 0.0, test RMSE 0.279 — 46 % worse than the two-level tree. Every training point is predicted exactly because every training point has its own box. This is memorisation in its purest form, and the depth sweep shows the whole curve: training RMSE falls monotonically (0.215, 0.204, 0.199, 0.189, 0.159) while test RMSE bottoms out at depth 2–3 (0.191, 0.190) and climbs from depth 5 (0.201) to depth 8 (0.229). A shallow tree has high bias and low variance — it cannot express much, but what it expresses is stable. A deep tree is the reverse. The sweet spot is found on validation data, never on training data, because the training curve never turns.
Four hyper-parameters control the trade-off, all of them variance dials: max_depth caps successive splits (typically 2–7); min_samples_leaf forbids leaves smaller than a given size (5–50) — the last line shows that a fully grown tree with min_samples_leaf=50 has 19 leaves and a test RMSE of 0.195, nearly all the damage undone by one constraint; min_samples_split forbids splitting small nodes; max_leaf_nodes caps complexity directly. Set one or two, tune on a held-out fold.
Bagging and random forests: averaging fixes variance
If a deep tree’s problem is variance — refit it on a slightly different sample and it changes a lot — the classical remedy is to average. Bagging (bootstrap aggregating) draws \(B\) bootstrap samples of the training rows, fits one deep tree to each, and averages the \(B\) predictions. Each tree is noisy in its own way; the noise averages out while the signal, which all trees share, survives. Averaging keeps the low bias of deep trees and divides the variance.
A random forest adds one more source of randomness: at each split only a random subset of the features is offered. This seems perverse — why deny the tree its best split? — but it decorrelates the trees. If one feature dominates, plain bagging grows \(B\) near-identical trees whose average is barely less variable than one of them; forcing some splits onto other features makes the trees genuinely different, and the average of less correlated estimates has lower variance. The forest below uses 100 trees with a leaf size of 20; the cell compares it with OLS on the same eight features and plots its feature importances.
The forest’s test RMSE is 0.191 — identical to OLS on these eight features — with a lower training RMSE (0.185 against 0.205). With 1 146 rows and a weak signal, the nonlinearity buys nothing out of sample; the forest fits the training set better and generalises the same. On the full 150 000-loan notebook the picture changes: there the forest lifted test adjusted R² from 0.163 to 0.216, because with a hundred times more rows there is enough data to estimate the interactions and thresholds reliably. Flexible models need data in proportion to their flexibility.
The bar chart ranks features by impurity importance: for every split on feature \(j\) in every tree, record the SSE reduction; sum, and average across trees. term (0.289), balance (0.249) and int_rate (0.178) dominate — the same three the t-tests and the best-subset search flagged. That agreement is reassuring, but read the number for what it is.
Impurity importance is a description of the forest, not of the world. It is biased toward features with many distinct values — a continuous variable offers hundreds of cut-points, a 0/1 dummy offers one, so the continuous one gets more chances to be chosen whatever the truth. Correlated features share the credit unstably: two near-copies split the importance between them at random across trees. And a feature that matters only inside an interaction can score low. For “does this variable matter?”, use permutation importance on test data (shuffle the column, measure the loss in test score) or a t-test in a model you can interpret.
Boosting: fitting the residuals
Averaging attacks variance. Boosting attacks bias. Start with a constant prediction, the training mean. Fit a small tree — a one-split stump, say — not to \(y\) but to the current residuals \(y - f_0(x)\), and add a shrunken copy to the model: \(f_1 = f_0 + \nu T_1\), with a learning rate \(\nu\) well below 1. Compute the new residuals, fit another small tree, add it, repeat. After \(M\) rounds,
\[f_M(x) = \bar y + \nu \sum_{m=1}^{M} T_m(x),\]
and each tree has corrected what the previous ones left behind. For squared-error loss the residual is exactly the negative gradient of the loss with respect to the prediction, which is why the method is called gradient boosting and why it generalises to any differentiable loss — logistic loss for classification, quantile loss for a median. Five rounds by hand make the mechanism concrete, and GradientBoostingRegressor then does 100 rounds with depth-3 trees.
Five one-split stumps with \(\nu = 0.5\) take the test RMSE from 0.200 to 0.189 — the forest’s level, reached with five splits in total. The full GradientBoostingRegressor lands at 0.192 on test, the forest at 0.191, OLS at 0.191: three models within a thousandth of each other on a weak signal. When the signal is weak, model choice matters less than the analyst expects, and honest reporting matters more.
The knobs: n_estimators (more rounds, more capacity, eventually overfitting); learning_rate (smaller steps need more trees and generalise better — the recipe is a small rate, 0.01–0.1, with the number of trees chosen on a validation fold); max_depth of 1–3, so that each tree stays a weak learner; and subsample < 1, bagging’s row sampling inside boosting. Random forests are the robust default; boosting, tuned, usually wins on tabular data by a margin that depends on how much signal there is to find.
Predicting the stock market: why you cannot shuffle
The lecture notebook’s last example turns boosting on prices. The target is whether the S&P 500 rises tomorrow; the features are five lagged daily returns, the trailing 10-day volatility and the 20-day momentum, all lagged by at least one day so that they were known at yesterday’s close. The question that decides whether the exercise means anything is how to split the data.
With cross-sectional loans, a seeded shuffle was right. With a time series it is exactly wrong: a shuffled split puts 2023 rows in the training set and 2016 rows in the test set, so the model is trained on the future of the rows it is tested on — it learns the regime and drift of 2023 and is then asked about days that preceded them. TimeSeriesSplit instead trains on an earlier block and tests on the next, with an expanding window: each fold’s training set is everything before its test block. That is the only honest test of a forecast; the toy call on ten indices at the end of the cell shows the pattern.
2 495 trading days from 2015 to 2024, of which 53.8 % closed up. Seven features, nine columns including the target and the raw return. The toy split shows the expanding window: fold 1 trains on indices 0–1 and tests on 2–3; fold 4 trains on 0–7 and tests on 8–9. Now the bug that every student writes at least once.
The test AUC is 1.000, and no honest model of daily direction reaches it. The leak is range(0, 5): the first feature is lag0 = r.shift(0), today’s return — the very quantity whose sign is the target. The last three printed rows show it directly: lag0 is positive exactly when up is 1. The model has been handed the answer as a feature, and TimeSeriesSplit cannot protect against that, because the leak is within each row, not across time. The fix is range(1, 6), as in the cell above. The general rule: every feature at date \(t\) must be computable from information available at the close of \(t-1\). An AUC that looks too good is a bug until proven otherwise.
Walk-forward boosting on ten years
Four folds, each testing on roughly two years. Training AUC runs from 0.833 on the first fold down to 0.705 on the fourth — highest on the smallest training block, where memorising is easiest — while test AUC is 0.516, 0.494, 0.555, 0.514, a mean of 0.52. The model memorises the training block and forecasts the next one at barely better than a coin flip. The notebook’s much larger exercise — LightGBM on 125 MB of 5-minute bars with 140 features and weekly folds — lands at 0.538 ± 0.031: the same picture at a different scale. Daily direction is nearly unforecastable from its own past, and a boosting machine with 100 trees does not change that.
Is an AUC of 0.52 worth anything? On the 821 out-of-sample days when the model said \(p > 0.55\), the S&P averaged 0.089 % a day against 0.054 % across all test days; on its 296 bearish days, −0.027 %. The notebook’s version found 0.077 % against 0.023 %. A tiny edge, before costs, on a sample where the standard error of a daily mean is about 0.04 % — and the importance ranking (lag4 first) is the unstable kind the callout above warned about. The lesson is the protocol: shuffled, the test AUC would have been meaningless; unlagged, it would have been 1.0.
Pyodide has no lightgbm; the lecture notebook runs the walk-forward above on 5-minute bars with it. It is the same gradient-boosting idea engineered for speed and regularisation:
from lightgbm import LGBMClassifier, early_stopping
model = LGBMClassifier(n_estimators=5000, learning_rate=0.02, num_leaves=63,
subsample=0.8, colsample_bytree=0.8,
reg_alpha=1.0, reg_lambda=2.0, random_state=42, verbose=-1)
model.fit(X_train, y_train, eval_set=[(X_valid, y_valid)], eval_metric="auc",
callbacks=[early_stopping(200, verbose=False)])num_leaves = 63 is roughly depth 6; subsample and colsample_bytree are row and feature sampling inside boosting; reg_alpha and reg_lambda are L1 and L2 penalties on leaf values; early_stopping chooses n_estimators on the validation fold. Split with TimeSeriesSplit, never a shuffle.
Logistic classification for PD
The third factor. Default is 0/1, and a model of it must produce a probability. Why not run OLS on the 0/1 column? Because the fitted values of that linear probability model can fall below 0 or above 1, and because the error variance of a 0/1 outcome is \(p(1-p)\), which changes with \(p\), so the constant-variance assumption behind every standard error is violated by construction. Logistic regression fixes both by modelling the log-odds as linear:
\[\log\frac{p}{1-p} = \beta_0 + \beta^\top x \quad\Longleftrightarrow\quad p = \frac{1}{1 + e^{-(\beta_0 + \beta^\top x)}}.\]
The logistic function maps the whole real line into \((0, 1)\), so \(p\) is always a probability; the coefficients are fitted by maximum likelihood rather than least squares; and each \(\beta_j\) is read as the change in log-odds per unit of \(x_j\). The model is fitted on all 5 595 training loans, not the defaulted subset, because the question is now whether a loan defaults. The second half of the cell turns probabilities into decisions with a confusion table at two cutoffs.
int_rate (z = 9.8) and term (z = 5.7) raise the log-odds of default; fico (z = −4.8) lowers them; install (3.7) and balance (3.1) raise them; dti and open_acc are indistinguishable from zero. The coefficient on int_rate, 0.085, says each extra percentage point of rate multiplies the odds of default by \(e^{0.085} \approx 1.09\) — which is what you would expect, since the rate was set by the platform’s own risk assessment. The mean predicted PD is 0.205, exactly the training default rate: logistic regression fitted by maximum likelihood with an intercept is calibrated on average by construction.
A probability is not a decision. Flag a loan as “will default” when \(p\) exceeds a cutoff \(c\), and the confusion table counts the four outcomes on the 2 399 test loans. Two rates summarise it: the true positive rate TPR = P(flag | default) and the false positive rate FPR = P(flag | repaid). At the conventional cutoff of 0.5 the model flags 62 loans, catches 33 of the 458 real defaulters — a TPR of 7.2 % — and wrongly flags 1.5 % of the good loans. Predicted PDs average 0.2 and rarely exceed 0.5, so a 0.5 cutoff barely flags anyone. At cutoff 0.2 the model catches 62.7 % of defaulters but flags 34.6 % of good loans. Every cutoff is a trade between the two rates, and there is no cutoff that improves one without worsening the other.
The ROC curve and the cutoff that maximises profit
The ROC curve plots TPR against FPR for every cutoff at once, from \(c = 1\) (flag nobody: both rates 0) to \(c = 0\) (flag everybody: both 1). A model that ranks defaulters above good loans bows the curve toward the top-left corner; a coin flip gives the diagonal. The area under the curve (AUC) is the probability that a randomly chosen defaulter is ranked riskier than a randomly chosen good loan — a cutoff-free summary of ranking quality.
Test AUC is 0.702, training AUC 0.692 — no overfitting in a 9-parameter logit, and the test set happens to be marginally easier. A random defaulter is ranked riskier than a random good loan 70 % of the time. The 150 000-loan notebook reaches a similar 0.70; origination data simply do not separate the two classes more than this, because a great deal of what causes default (job loss, illness, divorce) happens after origination.
Which point on the curve should the lender choose? Not 0.5 — that number has no business content. Suppose approving a good loan earns \(g\) and approving a bad one loses \(\ell\), each weighted by its base rate. A loan is approved when it is not flagged, so the expected profit per applicant at cutoff \(c\) is
\[\Pi(c) = g\,(1 - \text{FPR}(c)) - \ell\,(1 - \text{TPR}(c)).\]
Moving along the ROC curve changes both rates; at the optimum \(d\Pi = 0\), which gives \(g\,d\text{FPR} = \ell\,d\text{TPR}\), or
\[\frac{d\,\text{TPR}}{d\,\text{FPR}} = \frac{g}{\ell}.\]
The optimal cutoff is where the slope of the ROC curve equals gain over loss. If losses per bad loan rise relative to gains, \(g/\ell\) falls, the optimum moves to where the curve is flatter — higher TPR, higher FPR — which is a lower cutoff: flag more applicants and accept more false alarms. The cell puts in the chapter’s own numbers: a gain of about 10 % interest on a good loan and a loss of EAD × LGD ≈ 0.70 × 0.89 ≈ 0.62 of principal on a bad one, and then asks whether boosting ranks better than the logit.
The slope target \(g/\ell\) is 0.68, so the optimum sits where the ROC curve has flattened to that gradient: cutoff 0.166, flagging 77.5 % of defaulters at the price of 49.7 % of good applicants. Approving everyone loses 3.75 cents per applicant; at the optimal cutoff the book earns 1.41 cents. The number 0.5 was never a candidate — at 0.5 the model flags 7 % of defaulters and the book still loses money. Where exactly the cutoff lands depends on \(g\) and \(\ell\), which are business inputs, not statistical ones; the statistician’s job is the curve, the manager’s is the slope.
Boosting reaches a training AUC of 0.745 against the logit’s 0.692, and a test AUC of 0.703 against 0.702. The flexible model gained five points in-sample and one tenth of a point out of sample. When the signal is weak, a flexible model mostly finds more noise — and the logit has the advantage of coefficients that can be explained to a regulator.
Overfitting: the deep tree memorised every training row (one leaf each) and its variance made it 46 % worse out of sample. max_depth (cap the number of successive splits) and min_samples_leaf (forbid small leaves; 50 gave 19 leaves and test RMSE 0.195) are the variance dials, set on validation data — never on training data, whose error never turns.
Profit per applicant is \(g(1-\text{FPR}) - \ell(1-\text{TPR})\); setting its derivative along the curve to zero gives \(d\text{TPR}/d\text{FPR} = g/\ell\). A larger \(\ell\) lowers \(g/\ell\), so the optimum moves to a flatter part of the curve — higher TPR and FPR — which is a lower cutoff: flag more applicants. In the chapter, \(g/\ell = 0.68\) put the cutoff at 0.17, not 0.5.
Causal Analysis: Prediction vs Causation, Confounders, Difference-in-Differences, Instrumental Variables
Two questions that look alike
Everything in the first three sections answered one question: given \(x\), what is \(y\)? The PD model says that a borrower with a 20 % rate and a 60-month term is more likely to default than one with a 7 % rate and 36 months, and it is right. A manager reads that table and asks a question that looks identical and is not: if I change \(x\), what happens to \(y\)? Suppose the platform observes that “borrowers who use our budgeting tool repay 20 % more”. Four actions are on the table. Flag tool users as low-risk when pricing their next loan. Add tool usage as a feature in the PD model. Report the correlation in the annual review. Or push every borrower to install the tool, expecting 20 % better repayment.
The first three need only an association: any variable correlated with repayment improves a prediction, and it does not matter why it is correlated. The fourth needs something much stronger — that changing tool usage changes repayment — and the same regression coefficient that serves the first three will, in general, answer the fourth wrongly.
The potential outcomes framework makes the gap precise. Each borrower \(i\) has two potential repayment outcomes, \(Y_i(1)\) with the tool and \(Y_i(0)\) without, and we observe exactly one. The causal effect \(Y_i(1) - Y_i(0)\) is never observed for anyone — causal inference is a missing-data problem. What we can compute is the observed gap between users and non-users, which decomposes as
\[\underbrace{E[Y \mid T=1] - E[Y \mid T=0]}_{\text{observed gap}} = \underbrace{E[Y(1)-Y(0)\mid T=1]}_{\text{causal effect on users}} + \underbrace{E[Y(0)\mid T=1] - E[Y(0)\mid T=0]}_{\text{selection bias}}.\]
The second term asks: would the users have repaid better than the non-users even without the tool? If the people who install a budgeting tool were already disciplined about money, the selection term is large and positive and the observed 20 % gap is mostly selection. Randomised assignment kills the selection term by construction — that is why an A/B test is the gold standard — but most business data are observational, and the rest of this section is about what to do then.
A great predictor that is a useless lever
The cleanest way to see the problem is to build a world in which we know the answer. A hidden trait \(U\) — financial discipline — drives both tool adoption \(T\) and the repayment score \(Y\). By construction, the tool itself has zero effect on repayment.
36 % of borrowers use the tool, and the correlation between usage and repayment is 0.48. The naive regression of \(Y\) on \(T\) gives a coefficient of 2.17 with a t-statistic of 24 and an R² of 0.23. By every standard in Sections 3.1–3.3 this is an excellent predictor: a strongly significant coefficient, a respectable R², and it would improve a PD model if added as a feature. As a policy, pushing the tool would move repayment by exactly 0, because the world was built that way. The 2.17 is entirely selection: disciplined borrowers adopt the tool, and disciplined borrowers repay. The t-statistic of 24 is not a lie — the association is real and would replicate in any sample — it just answers a different question from the one the manager asked.
Before interpreting any coefficient as an effect, draw the arrows. Here they are \(T \leftarrow U \rightarrow Y\): \(U\) is a confounder, a common cause of treatment and outcome, and the path through it is a “backdoor” that carries association without causation. A predictive model does not care which arrows exist; a causal claim depends on nothing else driving both \(x\) and \(y\).
Regression adjustment: control for the confounder
If the confounder is measured, the fix is to put it in the regression. Conditional on \(U\), the comparison is between tool users and non-users with the same discipline, the backdoor path is closed, and the coefficient on \(T\) estimates the causal effect.
Add \(U\) as a second regressor. Will the coefficient on \(T\) stay near 2.17, collapse toward 0, or flip sign?
The adjusted coefficient on \(T\) is 0.066 — zero within noise — and the coefficient on \(U\) is 1.94, close to the true 2.0. In a second version of the world where the tool genuinely adds 0.5, the naive regression says 2.67 and the adjusted one says 0.45: the adjustment recovers the effect both when it is absent and when it is present.
Two cautions. First, this works only because \(U\) was measured; in the real platform “financial discipline” is not a column. Second, the rule is include confounders, not include everything. A mediator — a variable on the path from \(T\) to \(Y\), such as “budget alerts received” — must be left out, because conditioning on it removes the effect you want to measure. A collider — a variable caused by both \(T\) and \(Y\), such as “selected for a testimonial” — must also be left out, because conditioning on it opens a spurious path (among professional basketball players height and skill are negatively correlated: being short and skilled is the only way a short player makes the league). Which variables to include is decided by the graph, not by their t-statistics.
When the confounder is not in your data
U (unobserved)
/ \
v v
T -> YRegression adjustment cannot close a backdoor through a variable you do not have, and no amount of data on \(T\) and \(Y\) alone can recover the effect — more rows only give a more precise estimate of the wrong number. If you can randomise, do: an A/B test that assigns the tool by lottery makes \(T\) independent of \(U\), and a difference of means is then a causal estimate with a standard error from a two-sample t-test. When you cannot randomise, two observational designs replace the missing control with a structure. Difference-in-differences uses a control group that shares the treated group’s time trend. Instrumental variables use a source of variation in \(T\) that has nothing to do with \(U\). Each replaces the untestable “no unmeasured confounder” with a different, and sometimes more defensible, assumption.
Difference-in-differences: the gap between gaps
Suppose the tool is rolled out to one region (the treated group) between period 0 and period 1, and not to another (the control). Two naive comparisons are available and both are wrong: treated versus control after the rollout picks up any pre-existing level difference, and the treated region before and after picks up whatever happened to everyone — a recession, a rate change, seasonality. DiD subtracts one from the other:
\[\hat\tau_{\text{DiD}} = (\bar Y^{\text{tr}}_{1} - \bar Y^{\text{tr}}_{0}) - (\bar Y^{\text{ctl}}_{1} - \bar Y^{\text{ctl}}_{0}).\]
The control group’s change stands in for what the treated group would have done without the policy. That works if, absent treatment, both groups would have moved in parallel — the parallel trends assumption. Level differences are fine, because they cancel in the first differences; different trends are fatal, because the control’s change is then not the treated group’s counterfactual. The simulation gives the treated group a level 1.0 higher, gives everyone a common drift of 0.5, and adds a policy effect of 0.8 to the treated group in period 1.
The four cell means are 5.013 and 5.466 for the control group, 6.001 and 7.315 for the treated. The naive post-period gap is 1.849 — the true 0.8 plus the 1.0 level difference. The naive before-after change in the treated group is 1.314 — the true 0.8 plus the 0.5 common drift. The DiD is $1.314 - 0.453 = $ 0.861, within a standard error of 0.8.
The regression version, \(y = \alpha + \beta\,\text{group} + \gamma\,\text{period} + \delta\,(\text{group} \times \text{period})\), returns exactly the same four numbers as coefficients — \(\hat\alpha = 5.013\) is the control’s baseline, \(\hat\beta = 0.988\) the level difference, \(\hat\gamma = 0.453\) the common drift, and $= $ 0.861 the DiD — now with a standard error of 0.051, so that the effect is 17 standard errors from zero and 1.2 from the truth. The interaction term is the whole design: the effect that exists only for the treated group and only after treatment. Parallel trends cannot be tested in the post period (the treatment is in the way), but with several pre-periods one can and should check that the two groups moved in parallel before the rollout — the standard “event-study” plot that accompanies every serious DiD.
Instrumental variables: borrow a lottery
An instrument \(Z\) is a variable that moves \(T\) but touches \(Y\) only through \(T\). Three conditions must hold. Relevance: \(Z\) shifts \(T\) — testable from the first-stage regression. Independence: \(Z\) is unrelated to the unobserved \(U\) — usually argued from the design (a lottery, a weather shock, an arbitrary administrative rule). Exclusion: \(Z\) has no path to \(Y\) except through \(T\) — never verifiable from data, defended on substantive grounds, and where most instrumental-variable arguments are attacked. Classic instruments are distance to the nearest college for schooling (Card 1995), the draft-lottery number for military service (Angrist 1990), and rainfall for agricultural income. On the lending platform, a randomised in-app prompt would do: it moves adoption, it was assigned by lottery, and it has no reason to affect repayment except through the tool.
Two-stage least squares turns the instrument into an estimate. Stage 1 regresses \(T\) on \(Z\) and keeps the fitted values \(\hat T\) — the part of \(T\) that the instrument explains, which by independence is free of \(U\). Stage 2 regresses \(Y\) on \(\hat T\). With a single instrument the whole procedure collapses to a ratio, \(\hat\beta_{\text{IV}} = \text{cov}(Z, Y)/\text{cov}(Z, T)\): how much \(Y\) moves per unit of \(Z\), divided by how much \(T\) moves per unit of \(Z\). The cell does both stages by hand with np.linalg.lstsq, then repeats the exercise with an instrument that barely moves \(T\).
The true effect of \(T\) on \(Y\) is 1.0. OLS says 2.013 — twice the truth — because \(U\) pushes \(T\) and \(Y\) in the same direction and OLS attributes all of that co-movement to \(T\). Two-stage least squares lands on 1.027, and the covariance ratio gives the identical 1.027, as it must with one instrument. The instrument has done what regression adjustment could not: it recovered the effect without ever observing \(U\).
The second half shows the price. When the instrument’s coefficient in the first stage is cut from 0.5 to 0.05, the first-stage F-statistic falls from 637 to 7.0, and the 2SLS estimate drifts to 0.72 with a standard error many times larger. A weak instrument explains so little of \(T\) that \(\hat T\) is mostly noise, and the second stage is regressing \(Y\) on noise; the estimate is biased toward the OLS value in finite samples and its confidence interval is unreliable. The rule of thumb is a first-stage F above 10 — otherwise do not trust the second stage, and say so.
| Design | Replaces the missing control with | Fails when |
|---|---|---|
| Regression adjustment | the measured confounder | a confounder is unmeasured |
| Difference-in-differences | a control group’s time trend | trends are not parallel |
| Instrumental variables | an exogenous shock to \(T\) | \(Z\) is weak or has its own path to \(Y\) |
Each row of the table trades one untestable assumption for another, and none of them is free. The discipline that matters is not choosing the cleverest design but stating, in the same sentence as the estimate, which assumption it rests on: “the tool raises repayment by 0.86, assuming the two regions would have moved in parallel”. A causal number without its assumption is a prediction wearing a disguise.
Zillow Offers used the company’s Zestimate models to buy homes directly, renovate and resell them. The models were good predictors of sale prices as observed. Once Zillow acted on them — bidding at the model’s price — the data-generating process changed: sellers accepted the offers the model had overpriced and walked away from the underpriced ones (adverse selection), and a 2021 market shift moved prices faster than the model updated. In Q3 2021 Zillow wrote down its inventory by $304 million, disclosed roughly 7 000 unsold homes, and on 2 November 2021 shut the business and cut about 25 % of its workforce. A model tuned on test RMSE answers “how close is \(\hat y\) to \(y\) under the world that produced the training data?”. The moment predictions become actions, the question is causal, and the selection that follows your action is a confounder the test set never contained.
The hidden trait \(U\) (discipline) drives both adoption and repayment: \(T \leftarrow U \rightarrow Y\). The association is real and replicable — it is a fine predictor — but it is entirely selection bias, \(E[Y(0)\mid T=1] - E[Y(0)\mid T=0]\). Adding \(U\) to the regression closes the backdoor and the coefficient collapses to 0.066. A t-statistic measures the strength of an association, never its direction of causation.
DiD: absent treatment, the treated and control groups would have followed parallel trends — untestable after treatment, checkable with pre-treatment periods. IV: relevance (testable: first-stage F > 10), independence of \(Z\) from unobserved confounders (argued from design), and exclusion — \(Z\) affects \(Y\) only through \(T\) — which can never be verified from data alone.
Cross-Sectional Attention Features: Transformers for the Stock Cross-Section
One month, one matrix
Every month a systematic equity fund sees the same object. At month-end \(t\) there is a matrix \(X_t \in \mathbb{R}^{N \times K}\) — \(N\) stocks by \(K\) characteristics — and a vector \(y_t \in \mathbb{R}^N\) of next month’s returns, the target. The problem is cross-sectional: rank the \(N\) stocks within the month so that the top outperforms the bottom. The single-series machinery of Section 4.3 is replaced by a stack of monthly snapshots, and the modelling questions change with it: how to standardise a snapshot without using other snapshots, what “peer group” a stock should be judged against, and how to score a ranking rather than a point forecast.
25 335 rows: 119 months × 213 stocks in 11 sectors, from January 2015 to November 2024. Five price-based characteristics are the features: last month’s return ret_1_0, the 12-month return skipping the most recent month ret_12_1 (the standard momentum signal), 21-day realised volatility, log dollar volume over 126 days, and the price as a fraction of its 252-day high. The target ret_next at month \(t\) is the return over month \(t+1\) — the check on AAPL confirms it is exactly ret_1_0 at the next row, 0.100776 in both places. The descriptives show why raw characteristics cannot go straight into a model: log_dvol_126d averages 19.5 while prc_highprc_252d is a fraction near 0.9, and ret_12_1 runs from −0.93 to +9.34 (a stock that rose tenfold in a year).
The universe is today’s large caps, tracked back to 2015. Nothing that was large in 2015 and has since failed, shrunk or been delisted is in the file. That is survivorship bias, and it flatters any backtest: the average stock in the file did better than the average stock a 2015 investor could have bought, and the losers that would have populated the short leg of a long-short portfolio are simply absent. The production rule is to build the universe as of each month from a point-in-time index membership file. This file is fine for learning the mechanics — never quote its Sharpe ratio as evidence.
Standardise per snapshot, never across time
A z-score of rvol_21d computed over the whole 2015–2024 file uses the mean and standard deviation of nine years that, for the 2015 row, had not yet happened. A model run live in January 2015 would have seen a different number, so the backtest no longer reproduces what the model could have known: a look-ahead hidden inside a preprocessing step. The transform is monotone, which is why the bug is easy to miss — the ordering within a month is unchanged — but the scale now depends on the future, and any tree that splits on a threshold will split differently.
The fix is to standardise within each month, using only the cross-section observed at \(t\): subtract the month’s median and divide by the month’s standard deviation. The median rather than the mean, because a single stock that rose tenfold should not shift the centre; and clip at ±3, so that one extreme month cannot dominate a tree split. In production the same transform is applied within each month × sector, so that a stock is judged against the peers it is actually compared with on the day — the _s columns below.
groupby(...).transform returns a frame aligned to the original rows, so the subtraction and division are row-by-row against each row’s own month (or month × sector). The largest sector median after standardising is exactly 0, as it must be; the clipped columns have standard deviation 0.91 (market-relative) and 0.97 (sector-relative), slightly under 1 because of the clipping, and a mean of 0.07 because the distributions are right-skewed and the median sits below the mean. Ten features so far: five market-relative, five sector-relative.
Why a sector median is a crude peer group
“Sector” is a label assigned by an index provider. It is stable, cheap and often wrong for the purpose at hand. In October 2022 NVIDIA’s behavioural peers — high dollar volume, high volatility, far below the 52-week high after a year of falling — included AMD, which shares its sector, but also Netflix, Alibaba and Tesla, none of which is “Technology”. A sector-relative z-score compares NVIDIA with Apple, Microsoft and Oracle, which in that month looked nothing like it.
Attention builds the peer group from the data instead. The construction is the one at the core of the transformer, and it is worth understanding from first principles because it is far simpler than its reputation. Start with the matrix \(X\) of \(N\) stocks by \(K\) standardised characteristics. Three linear maps produce three views of each stock: a query \(q_i = W_q^\top x_i\), “what am I looking for in a peer?”; a key \(k_j = W_k^\top x_j\), “what do I offer as a peer?”; and a value \(v_j = W_v^\top x_j\), “what information do I contribute once chosen?”. The score of stock \(j\) as a peer of stock \(i\) is the inner product \(q_i \cdot k_j\), scaled by \(\sqrt{d}\) (the dimension of the keys, so that the scale of the scores does not grow with \(d\)) and by a temperature \(\tau\). The scores in row \(i\) are then passed through a softmax,
\[A_{ij} = \frac{\exp(q_i \cdot k_j / \sqrt{d}\,\tau)}{\sum_{l} \exp(q_i \cdot k_l / \sqrt{d}\,\tau)},\]
which turns them into positive weights summing to one — a probability distribution over peers. Finally the context of stock \(i\) is the weighted average of the values, \(c_i = \sum_j A_{ij} v_j\). In matrix form, with \(Q = XW_q\), \(K = XW_k\), \(V = XW_v\):
\[A = \text{softmax}_{\text{rows}}\!\left(\frac{QK^\top}{\sqrt{d}\,\tau}\right), \qquad C = AV.\]
Now set \(W_q = W_k = W_v = I\). The score becomes \(x_i \cdot x_j / \sqrt{K}\) and the weight becomes proportional to \(\exp(x_i \cdot x_j / \sqrt{K})\): a similarity kernel in characteristic space. Stocks that look like \(i\) — the same sign and size of momentum, volatility, volume, drawdown — get high weight; stocks that look unlike it get weight near zero; and the sector label never enters. The context \(c_i\) is then the characteristic vector of \(i\)’s data-driven peer group, and the deviation \(x_i - c_i\) says how \(i\) differs from that group. That is precisely what a sector-relative z-score was trying to be, with the peer group chosen by the data instead of by an index committee.
Twelve lines. The only numerical care is the subtraction of the row maximum before exponentiating, which prevents overflow and does not change the softmax (a constant in the exponent cancels between numerator and denominator). NVIDIA’s attention row in October 2022 gives AMD 0.353, NVIDIA itself 0.263, Netflix 0.154, Alibaba 0.077, Tesla 0.055, Amazon 0.028 — six stocks carry more than 1 % of the weight, and they come from three sectors. The context vector is the peer-weighted characteristic vector: NVIDIA’s own 12-1 momentum is −1.49 standard deviations against −1.38 for its peers, its dollar volume 3.0 (clipped) against 2.52, its distance from the 52-week high −2.67 against −2.59. In that month NVIDIA largely was its peer group — a cluster of large, liquid, volatile stocks in deep drawdown — and no sector median would have seen it.
Temperature: the dial between two crude peer groups
The temperature \(\tau\) divides the scores before the softmax, and it controls how concentrated the weights are. As \(\tau \to 0\) the scores are magnified, the softmax approaches a hard maximum, and all the weight goes to the single most similar stock. As \(\tau \to \infty\) every score goes to 0, every weight goes to \(1/N\), and the context collapses to the cross-sectional mean.
With identity weights and \(\tau = 10^6\), what does NVIDIA’s context vector become? And at \(\tau = 0.05\), does NVIDIA attend mostly to itself?
At \(\tau = 10^6\) the context is (0.02, 0.10, 0.19, 0.00, −0.14) — exactly the cross-sectional mean of \(X\), which per-month standardisation has already removed, so the feature carries nothing new. At \(\tau = 0.05\) the weight on AMD is 0.997 and NVIDIA’s self-weight is 0.003: NVIDIA is not even its own nearest point, because the score is an inner product rather than a distance, and AMD’s vector, pointing the same way with a larger norm, scores higher. At \(\tau = 5\) the weights are nearly uniform (self 0.025, largest other 0.027). \(\tau = 1\) sits in between — a handful of behavioural peers — and is the setting used below; it is a hyper-parameter, and the exercise at the end of the section is about what happens when it moves.
What Kelly, Kuznetsov, Malamud and Xu actually train
The idea of putting attention over the stock cross-section rather than over words is due to Kelly, Kuznetsov, Malamud and Xu, “Artificial Intelligence Asset Pricing Models” (2025). Their architecture differs from the toy above in three ways. Each layer uses linear attention, \(A = XW_qW_k^\top X^\top\) with no softmax, followed by a feed-forward block and a residual connection, \(X^{(\ell+1)} = X^{(\ell)} + \text{FFN}(AX^{(\ell)}W_v)\), stacked \(L\) times. The final layer maps each stock’s row to a portfolio weight, and every \(W\) is trained end-to-end to maximise the Sharpe ratio of that portfolio — not to minimise a squared error. And they train on thousands of stocks over six decades. Attention is what lets a stock’s weight depend on the whole cross-section: the peer group is learned, for the objective the fund actually cares about.
The weights stay fixed at \(W = I\) (a seeded random \(W\) works the same way and is the “random features” view of the same construction). The attention output is used only as feature engineering: for every stock-month, the context \(C\) — peer-weighted characteristics — and the deviation \(X - C\), “how I differ from my attention-weighted peers”. The learner on top is the boosting machine from Section 4.3, fitted to next month’s return by squared error. Nothing is trained on a Sharpe objective, and the comparison is between three feature sets under one protocol, not between architectures.
attention_features runs one attention pass per month — only month-\(t\) rows enter month \(t\)’s attention matrix, so there is no look-ahead — and returns ten columns per stock-month, aligned to the original frame by index rather than by position (the groupby does not preserve row order, and assigning by position would scramble stocks across months). Fifteen features now: five market-relative characteristics, five peer contexts, five deviations. NVIDIA’s October 2022 row reads: momentum −1.49 against a peer context of −1.38, a deviation of −0.10; price-to-high −2.67 against −2.59, a deviation of −0.08. With a self-weight of 0.26 and AMD at 0.35, NVIDIA’s deviations from its peer group are small, which is the correct description of that month.
The rest of the cell sets up the protocol. The test period starts in January 2019, giving 71 test months; the model is refitted every 12 months on all rows dated strictly before the block and predicts the next 12 — 6 refits, the first on 10 212 rows (a monthly refit would be 71 fits at ten times the runtime, with the same conclusion). The learner is HistGradientBoostingRegressor with depth-2 trees, learning rate 0.05, 200 rounds and a minimum leaf of 50 — Section 4.3’s recipe in its faster, histogram-based form.
One subtlety deserves a sentence. The row dated 2018-12-31 has ret_next equal to January 2019’s return, which is realised on 31 January 2019 — the same instant at which the 2019-01-31 features are observed. So when we predict the 2019-01-31 row, the 2018-12-31 target is already known, and training on all rows dated strictly before the test block is legal. A two-month-ahead target would need a one-month gap between the last training row and the first test row; a one-month-ahead target does not.
Two scores are computed per month. The rank IC is the Spearman correlation between the model’s prediction and the realised return across the 213 stocks — the natural score for a ranking problem, immune to the scale of the predictions. The top-20 minus bottom-20 spread is the equal-weight return of the 20 highest-ranked stocks minus the 20 lowest — the return of the simplest long-short portfolio the ranking could produce, before costs.
Three feature sets, one protocol
Form an expectation before running. The raw five features are the smallest set and the hardest to overfit; the sector-relative set is the production invariant; the attention set has data-driven peers. Which will win — and will the difference be detectable?
Mean rank IC: −0.015 (t = −1.04) for the raw features, −0.014 (t = −0.99) with sector features, −0.010 (t = −0.67) with attention. All three are zero within noise — the standard error of a mean IC over 71 months is about 0.015 — and all three are, if anything, slightly negative. The long-short spread tells a slightly different story: +7.7 %, +9.3 % and +2.8 % a year, with Sharpe ratios of 0.47, 0.62 and 0.17. The extremes of the ranking do a little better than the middle even though the ranking as a whole does not, which is common: the tails of a prediction are where the model is most confident and the signal, if any, is most concentrated.
Attention features did not help here. The attention set has the least negative IC and the lowest Sharpe; the ranking between the sets flips with the score you read; none of the differences approaches its standard error. That is the honest result, and the honest result is the deliverable. Five price-based characteristics on 213 large caps over 71 months is a small, noisy problem — KKMX have thousands of stocks and six decades — and a feature that needs that much data to show its value cannot be expected to show it here. What the protocol delivers is the ability to tell: had the split been shuffled, or the standardisation done over the full sample, the attention set could easily have “won”, and the win would have been an artefact.
Is a Sharpe ratio of 0.6 evidence?
The plot shows three cumulative spreads drifting upward with large swings, and the temptation is to read the 0.62 Sharpe as a modest real edge. The arithmetic says otherwise. The t-statistic of a monthly spread with annualised Sharpe \(S\) observed for \(T\) months is \(S\sqrt{T/12}\): with \(T = 71\), the sector set’s 0.62 gives t = 1.51, the raw set’s 0.47 gives 1.15, attention’s 0.17 gives 0.40. To reach t = 2 the sector set would need 124 months of the same performance, the raw set 18 years. A Sharpe ratio is a point estimate with standard error roughly \(\sqrt{12/T} \approx 0.41\) here — wider than the differences between the three lines — and the file is survivorship-biased in the spread’s favour on top of that.
Changing the peer group
The peer-group definition has two dials: the temperature, and which characteristics form the keys. The cell rebuilds the attention features using only momentum and volatility as keys — so that “similar” means “similar momentum and volatility”, ignoring volume and drawdown — and reruns the walk-forward on the attention set.
Momentum-and-volatility keys move the mean IC to −0.025 (t = −1.64), from −0.010. Try the other settings and the pattern is the same: \(\tau = 3\) leaves the IC at −0.010, \(\tau = 10\) gives −0.008, and \(\tau = 0.3\) gives −0.041 with t = −2.60 — the one “significant” number, and it appears on the fourth try. The sign of the change flips with the setting, and a t-statistic of 2.6 found after four searches is multiple testing, not a signal: with four independent tries, the chance that at least one exceeds |t| = 2 by luck is about 18 %. In KKMX the peer group is learned on a Sharpe objective with thousands of stocks and six decades; here it is a fixed kernel on 213 survivors, and the data cannot say which kernel is better. The chapter ends where it began: the protocol is the product, and the number it reports is only as good as the discipline behind it.
\(A_{ij} \propto \exp(x_i \cdot x_j / \sqrt{K}\,\tau)\) — a similarity kernel in characteristic space, independent of sector labels. As \(\tau \to 0\) all weight goes to the single most similar stock (at \(\tau = 0.05\) NVIDIA gave AMD 0.997); as \(\tau \to \infty\) every weight is \(1/N\) and the context is the cross-sectional mean, which per-month standardisation already removed.
The 2015 z-score would use the 2015–2024 mean and standard deviation, numbers a live model in 2015 could not have seen; the ordering within a month is unchanged but the scale — and therefore every tree threshold — depends on the future. Standardise within each month (median / std across the stocks observed at \(t\), clipped at ±3), and within month × sector when the peer group is the sector.
Chapter Wrap-up
You began with a ledger and ended with a protocol. Along the way the chapter built the three numbers a lender lives by — PD by logistic regression with AUC 0.70, EAD by regression with test adjusted R² of 0.11 to 0.14 depending on the columns, LGD as a constant near 0.89 because nothing at origination predicts recoveries — and combined them into an expected loss of 12.7 cents per dollar lent, which is why the loans in the file carry a 13 % rate. That is what a predictive model is for: numbers a business can price with.
The methodological lessons transfer to every other table you will meet. R² is not a scorekeeper; adjusted R² is a weak one; the test set is the only honest one, and a negative test R² means “use the mean”. Variables can be individually useless and jointly decisive, so interactions belong in the pool — searched by a procedure that judges columns given the others, which a marginal-correlation screen is not. t-tests and VIFs answer questions about coefficients, not predictions, and need assumptions the predictions never did. Trees find thresholds and interactions on their own and memorise noise just as readily; forests fix the variance by averaging, boosting fixes the bias by fitting residuals, and on a weak signal all three land within a thousandth of OLS. A probability is not a decision until a cutoff is chosen, and the right cutoff is where the ROC slope equals gain over loss — 0.17 here, not 0.5.
Then two lessons no predictive skill can substitute for. A coefficient with t = 24 can have a causal effect of exactly zero because a hidden trait drove both treatment and outcome; regression adjustment, difference-in-differences and instrumental variables each recover the effect under an assumption that must be named in the same sentence as the estimate. And a walk-forward test on the stock cross-section — standardised per month, never shuffled, refitted on the past only — reported that neither sector features nor an attention peer group beat zero on 213 surviving stocks over 71 months. The result is negative; the protocol that produced it is the skill.
Chapter 5 changes the frame. Everything here treated the coefficients as fixed unknowns to be estimated and tested. Bayesian methods treat them as uncertain quantities with distributions, combine the data with what was known before, and replace the p-value with a posterior — which, as you will see, makes the Sharpe ratio of Section 4.5 a rather different object. The slides for this chapter, with their multiple-choice checkpoints and in-class practice cells, are at the “Slides” link in the navigation bar; the sections, seeds and numbers there are the same as here.
← Chapter 2 · Contents · Chapter 5: Bayesian Methods →