TwoSLS

Closed-form linear IV / two-stage least squares

from _api_doc_utils import *

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 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.

4 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.

5 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)

6 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]

7 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	()