How the correction streams¶
Every quantity in the KSS estimator can be computed with the operations the regression already uses: passes over the rows, sorted by the streamed dimension, and solves with the reduced system for the non-streamed effects. That is what lets the correction run inside the same bounded memory as the fit. Each step is described below, in the notation of the problem and the correction.
Pruning. The articulation points are found on the bipartite graph of distinct worker–firm matches, which is built by one streaming aggregation and whose size is the number of matches, not of rows. The search is a linear-time depth-first traversal (Tarjan, 1972), repeated until no articulation point remains, after which the largest connected component is kept. On every panel tested the pruned sample matched xhdfe’s row for row.
Random projections. Computing \(P_{ii}\) exactly requires a linear solve with \(S_{zz}\) for every observation. KSS approximate it by random projection, following Johnson and Lindenstrauss (their “JLA”). Since \(P_{ii}\) is the squared norm of row \(i\) of the projection matrix \(P = Z S_{zz}^- Z'\),
is unbiased for \(P_{ii}\) (Achlioptas, 2003), and each draw costs one solve with \(S_{zz}\) for the whole sample rather than one per observation. Applying \(P\) to a vector \(q\) means computing \(Z\,S_{zz}^-\,Z' q\), which is exactly what the regression machinery does for an outcome: \(Z' q\) is a sum over rows, the solve with \(S_{zz}\) reduces to a solve with the reduced matrix for the non-streamed effects, and \((Z u)_i\) is computed row by row. One block of draws therefore costs two passes over the sorted rows and one reduced solve with as many right-hand sides as draws in the block.
The Rademacher entries \(q_{si}\) are generated by a counter-based generator keyed on the row’s position in the sorted working copy and the draw’s index, so they are never stored and do not depend on how the rows are divided into chunks or threads (draws and reproducibility says what they do depend on). The squared projections are accumulated into five running sums per row, held in memory-mapped files rather than in memory, which is the source of the additional scratch disk of about 90 bytes per row while the correction runs.
Two refinements, from a later note by the same authors distributed with the reference code, are applied: the estimates of \(P_{ii}\) and \(1-P_{ii}\) are both formed from the same draws and normalized to sum to one, which keeps the estimated leverage in \([0,1]\); and the bias that arises from dividing by an estimated \(1-P_{ii}\) in \(\hat\sigma_i^2\) is removed to order \(1/p\), using the higher moments that the same running sums provide.
The bias term. Rather than estimating each \(B_{ii}\) and forming \(\sum_i B_{ii}\hat\sigma_i^2\), the package estimates the trace \(\operatorname{tr}\big(Q S_{zz}^{-} Z'\hat\Omega Z S_{zz}^{-}\big)\), with \(\hat\Omega = \operatorname{diag}(\hat\sigma_i^2)\), by Hutchinson’s stochastic trace estimator, with random vectors in the space of the parameters rather than of the observations. Each draw requires two applications of \(S_{zz}^-\), each a reduced solve with passes over the rows, and no row-sized array of \(B_{ii}\) is ever formed. Both routes are unbiased for the same quantity, and the trace route has lower Monte Carlo error at equal cost (methodological differences §4.3).
Leaving out a match. Following the reference code, the match-level estimator partials the covariates out using the full model, collapses the data to one row per match (the mean outcome, weighted by the spell’s length), and applies the correction to that weighted regression, whose leverages are those of the matches. The collapse reads only the fit’s residual file, and the collapsed regression is a second, smaller streaming fit.
Standard errors. KSS show that their estimator can be written as a quadratic form \(\tilde y' C\tilde y\) whose kernel \(C\) has a zero diagonal, and derive its variance. The package computes it the same way: the first term of the variance is a sum over rows once \(C\tilde y\) is available, which takes solves and row passes, and the second, a trace, is again estimated by random draws. Like the reference implementation, the package smooths the leave-out variance estimates instead of using the split-sample estimator of the original paper, which would require a per-observation graph search; KSS’s Lemma 5 implies that the resulting standard errors are conservative. The weak-identification diagnostic streams too: its leading eigenvalues come from Lanczos iterations, each one solve with no extra passes over the rows. Standard errors and weak identification has the details.
Memory. The leave-out passes hold about a dozen float64 columns per row of
the chunk they are working on. scratch_mb (default 32) bounds that, shrinking
the chunk below batch_rows as needed; it changes memory only, never the
answer. At 5 million workers (42.5 million rows) the whole decomposition peaked
at 4.7 GB, against xhdfe’s 20 GB; see the KSS benchmarks.