from _api_doc_utils import *TwoSLS
Closed-form linear IV / two-stage least squares
1 Where it fits
Group: Causal inference
TwoSLS estimates linear instrumental-variables models. With endogenous regressors \(X_e\), exogenous controls \(X_c\), instruments \(Z\), and outcome \(y\), it runs the projection of the second-stage design on the instrument span and then estimates
\[ y = \alpha + X_e\beta_e + X_c\beta_c + u. \]
It supports multiple endogenous regressors and excluded instruments.
2 Criterion and first-stage projection
Let
\[ \tilde X=[\mathbf 1,\ X_{\mathrm{endog}},\ X_{\mathrm{exog}}], \qquad \tilde Z=[\mathbf 1,\ X_{\mathrm{exog}},\ Z_{\mathrm{excluded}}], \]
with the intercept omitted from both matrices only when intercept fitting is disabled internally. The implementation first solves \(\tilde Z\Pi\approx\tilde X\), forms \(\hat X=\tilde Z\hat\Pi=P_Z\tilde X\), and then solves \(\hat X\theta\approx y\). Thus the point estimate is
\[ \hat\theta = (\tilde X'P_Z\tilde X)^{-1}\tilde X'P_Zy. \]
Weighted fitting applies \(\sqrt{w_i}\) to \(y_i\), every regressor, and every instrument before constructing the projection. The sketched method applies one joint CountSketch to \(y\), \(\tilde X\), and \(\tilde Z\) and performs both stages in the sketched sample; it is an approximate point estimate.
3 Implementation walkthrough
The estimator implements 2SLS as two rectangular least-squares problems, avoiding an explicit \(n\times n\) projection matrix.
- The package concatenates endogenous regressors before exogenous controls for the reported second-stage design. The instrument design contains exogenous controls before excluded instruments. A common intercept column is prepended to both. It rejects fewer excluded instruments than endogenous regressors and any total instrument width smaller than the regressor width.
- A weighted fit scales \(\tilde X\), \(\tilde Z\), and \(y\) rowwise by the same \(\sqrt{w_i}\). It then solves the multi-response least-squares problem \(\tilde Z\Pi\approx\tilde X\) and forms \(\hat X=\tilde Z\hat\Pi\). No projection matrix \(P_Z\) is materialized.
- The second rectangular solve is \(\hat X\theta\approx y\). Algebraically it is the usual 2SLS estimator, but numerically both stages use the shared least-squares primitive rather than inverses of \(Z'Z\) and \(X'P_ZX\).
- The CountSketch path first constructs the full aligned designs, then hashes every observation once to a bucket and sign. That same draw is accumulated jointly into \(S\tilde X\), \(S\tilde Z\), and \(Sy\), preserving their row alignment. Both least-squares stages occur entirely in the sketched sample. The sketch width must be at least the larger design width.
- After either fit, the wrapper splits the intercept and slopes but stores the original endogenous, exogenous, instrument, outcome, and optional weight arrays. Full-sample structural residuals \(y-\tilde X\hat\theta\) drive every covariance and diagnostic, including after a sketched point estimate.
- Robust covariance is implemented through GMM-style parameter scores. It forms \(W=(Z'Z/n)^{-1}\), \(G=-Z'X/n\), transforms observation moments \(z_i\hat u_i\) by \(WG(G'WG)^{-1}/n\), and passes those parameter-score rows to the shared HC1, HAC, or cluster aggregator.
- The Anderson-Rubin method for one endogenous regressor constructs \(y_0=y-\beta_0x_e\), residualizes the null regression against included exogenous variables, and tests the excluded-instrument block jointly. This test is recomputed for the supplied null rather than inferred from the 2SLS Wald covariance.
The explicit fitted-regressor route is easy to inspect and stable relative to normal equations, but it materializes another dense \(n\times p\) matrix. The implementation checks order conditions, not statistical rank or instrument strength; near-collinearity is left to the numerical solve and diagnostics.
4 Inference and weak-instrument test
Define \(g_i=z_i\hat u_i\), \(G=-\tilde Z'\tilde X/n\), \(W=(\tilde Z'\tilde Z/n)^{-1}\), and \(A=G'WG\). The homoskedastic covariance is
\[ \widehat V_{\mathrm{vanilla}} = \frac{\hat\sigma^2}{n}A^{-1}, \qquad \hat\sigma^2=\frac{\hat u'\hat u}{n-p}. \]
The robust estimator replaces \(\hat\sigma^2W^{-1}\) with the empirical covariance \(\hat\Omega\) of \(g_i\):
\[ \widehat V = \frac1n A^{-1}G'W\hat\Omega WG A^{-1}. \]
The implemented \(\hat\Omega\) choices are HC1 iid scores, Bartlett Newey-West scores, and cluster sums with the same finite-sample corrections used by OLS. Wald tests use this covariance. The pairs bootstrap resamples all outcome, regressor, instrument, and weight rows together.
For one endogenous regressor, the Anderson-Rubin method tests \(H_0:\beta=\beta_0\) by regressing \(y-\beta_0x_{\mathrm{endog}}\) on the included exogenous variables and excluded instruments and jointly testing the excluded-instrument coefficients. It is unavailable with multiple endogenous regressors.
5 Performance and numerical behavior
The two least-squares stages scale with the instrument width as well as the regressor width, and inference forms dense instrument and parameter cross-products. Weak instruments do not necessarily cause a numerical error; they instead make \(\tilde X'P_Z\tilde X\) poorly conditioned and conventional Wald inference unreliable, which is why the scalar Anderson-Rubin test is exposed. CountSketch reduces the rows used in both point-estimation stages, but the fitted object retains the full sample and evaluates residuals and covariance there. Large instrument sets still incur quadratic storage and cubic dense-factorization costs in their width.
6 Python API
Constructor: cm.TwoSLS
Call fit(x_endog, x_exog, z, y). The predict(x) method expects a second-stage design matrix matching the fitted endogenous+exogenous columns. The robust linear covariance options match OLS: vanilla, hc1, newey_west, and cluster.
print(inspect.signature(cm.TwoSLS))()cls = cm.TwoSLS
display(HTML(html_table(["Public method"], public_methods(cls))))| Public method |
|---|
anderson_rubin_test(self, /, beta=0.0, vcov='hc1', lags=None, clusters=None) |
bootstrap(self, /, n_bootstrap, seed=None) |
fit(self, /, x_endog, x_exog, z, y) |
fit_sketch(self, /, x_endog, x_exog, z, y, sketch_size, seed=None) |
fit_weighted(self, /, x_endog, x_exog, z, y, sample_weight) |
predict(self, /, x) |
summary(self, /, vcov='hc1', lags=None, clusters=None) |
wald_test(self, /, r, q=None, vcov=None, lags=None, clusters=None) |
7 Minimal example
rng = np.random.default_rng(9)
n = 400
z = rng.normal(size=(n, 2))
x_exog = rng.normal(size=(n, 1))
v = rng.normal(size=(n, 1))
u = 0.6 * v[:, 0] + rng.normal(scale=0.4, size=n)
x_endog = z @ np.array([[0.8], [-0.4]]) + 0.2 * x_exog + v
y = 0.5 + 1.2 * x_endog[:, 0] - 0.7 * x_exog[:, 0] + u
model = cm.TwoSLS()
model.fit(x_endog, x_exog, z, y)
print(model.summary()['coef'])
print(model.summary(vcov='newey_west', lags=3)['coef_se'])[ 1.22456127 -0.6589091 ]
[0.04054546 0.03621148]8 summary() contract
The table below is generated by fitting the live class in this repository and then inspecting summary(). Shapes are shown because most values are plain NumPy arrays or scalars.
rng = np.random.default_rng(109)
n = 120
z = rng.normal(size=(n, 2))
x_exog = rng.normal(size=(n, 1))
v = rng.normal(size=(n, 1))
x_endog = z @ np.array([[0.8], [-0.4]]) + 0.2 * x_exog + v
y = 0.5 + 1.2 * x_endog[:, 0] - 0.7 * x_exog[:, 0] + 0.6 * v[:, 0] + rng.normal(size=n) * 0.4
model = cm.TwoSLS()
model.fit(x_endog, x_exog, z, y)
summary = model.summary()
display(HTML(html_table(["summary() key", "shape"], summary_shape_rows(summary))))| summary() key | shape |
|---|---|
intercept |
() |
coef |
(2,) |
intercept_se |
() |
coef_se |
(2,) |
vcov |
(3, 3) |
vcov_type |
() |