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.

Weights are analytic, not frequency weights. Zero-weight rows and zero-only clusters do not count toward summary, Wald, or Anderson-Rubin inference; HAC uses the active-row ordering. Prediction validates finite inputs and trained width, and failed refits clear old state. Sketched summaries record method, size, seed, original sample count, and inference_approximate=True; covariance does not account for sketch randomization.

3 Implementation walkthrough

The estimator implements 2SLS as two rectangular least-squares problems, avoiding an explicit n\times n projection matrix.

  1. 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.
  2. 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.
  3. 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.
  4. 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.
  5. 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.
  6. 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.
  7. 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.

Code 
print(inspect.signature(cm.TwoSLS))
()
Code 
cls = cm.TwoSLS
display(HTML(html_table(["Public method"], public_methods(cls))))
Public method
TwoSLS()
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.

Code 
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 ()
sketch_size ()
sketch_seed ()
inference_approximate ()
nobs ()
original_nobs ()
sketch_method ()