13.5 Instrumental variables in software

The function tsls below computes 2SLS with the correct standard errors and returns the residuals \(Y - Xb\); it is used in all four panels.

Software 13.5.1 (Simultaneity in a market, and the second-stage trap). A market with demand \(q = 100 - 2p + u_d\) and supply \(q = 20 + 3p - 4w + u_s\), where \(w\) shifts supply (weather, or input costs) but not demand, so it is a valid instrument for price in the demand equation:

> tsls <- function(y, X, Z) {               # 2SLS with correct standard errors
+   Xh <- Z %*% solve(crossprod(Z), crossprod(Z, X))
+   b <- solve(crossprod(Xh), crossprod(Xh, y))
+   u <- y - X %*% b; s2 <- sum(u^2) / (length(y) - ncol(X))
+   list(b = drop(b), se = sqrt(diag(s2 * solve(crossprod(Xh)))), u = drop(u))
+ }
> set.seed(4); n <- 200
> w  <- rnorm(n, 5, 1)                      # cost shifter (weather)
> ud <- rnorm(n, 0, 4); us <- rnorm(n, 0, 4)
> p  <- (80 + 4 * w + ud - us) / 5          # equilibrium price
> q  <- 100 - 2 * p + ud                    # demand: slope -2
> X <- cbind(1, p); Z <- cbind(1, w)
> iv <- tsls(q, X, Z)
> naive <- lm(q ~ fitted(lm(p ~ w)))        # second stage by hand
> rbind(OLS = coef(summary(lm(q ~ p)))[2, 1:2], IV = c(iv$b[2], iv$se[2]),
+       naive2nd = coef(summary(naive))[2, 1:2])
           Estimate Std. Error
OLS      -0.3910552  0.1727943
IV       -1.9211272  0.3387597
naive2nd -1.9211272  0.2563027

Reading the output. Least squares estimates the demand slope as \(-0.39\), nowhere near the true \(-2\): price is correlated with the demand shock, because a positive \(u_d\) raises the equilibrium price, so the regression traces out a mixture of the demand and supply curves. Its probability limit, from the variances of the simulated shocks, is \(-2 + \Cov (p, u_d)/\Var (p) = -2 + 3.2/1.92 = -0.33\). IV recovers \(-1.92\). The second stage run by hand gives the same estimate but a standard error of \(0.256\) instead of the correct \(0.339\): it understates the uncertainty by a quarter, because it computes the residuals from \(\hat p\) instead of \(p\).

Stata: ivregress 2sls q (p = w), with estat firststage for the first stage. EViews: tsls q c p @ c w, the instruments listed after @.

Software 13.5.2 (Weak instruments). \(2000\) samples of size \(100\) from \(y = x + u\), \(x = \pi (z_1 + \dots + z_5) + v\) with \(\Corr (u, v) = 0.5\), for five values of the first-stage coefficient \(\pi \):

> set.seed(5)
> weak <- function(pi, n = 100, l = 5, R = 2000) {
+   r <- replicate(R, {
+     Z <- matrix(rnorm(n * l), n); v <- rnorm(n)
+     u <- 0.5 * v + sqrt(0.75) * rnorm(n)  # Corr(u, v) = 0.5
+     x <- Z %*% rep(pi, l) + v; y <- x + u # true beta = 1
+     fs <- lm(x ~ Z); F <- summary(fs)$fstatistic[1]
+     xh <- fitted(fs); b <- sum(xh * y) / sum(xh * x)
+     e <- y - b * x; se <- sqrt(sum(e^2) / (n - 1) / sum(xh^2))
+     c(ols = sum(x * y) / sum(x^2), tsls = b, F = F,
+       reject = abs((b - 1) / se) > 1.96)
+   })
+   rowMeans(r)
+ }
> tab <- sapply(c(0.02, 0.05, 0.1, 0.2, 0.5), weak)
> colnames(tab) <- paste0("pi=", c(0.02, 0.05, 0.1, 0.2, 0.5)); round(tab, 3)
        pi=0.02 pi=0.05 pi=0.1 pi=0.2 pi=0.5
ols       1.499   1.495  1.476  1.418  1.222
tsls      1.479   1.409  1.255  1.098  1.015
F.value   1.030   1.263  2.075  5.150 26.328
reject    0.228   0.223  0.166  0.116  0.066

Reading the output. The rows give the mean least squares and 2SLS estimates (the truth is \(1\)), the mean first-stage \(F\), and the rejection rate of a nominal \(5\%\) \(t\) test of the true value. With \(\pi = 0.02\), \(F \approx 1\) and 2SLS, at \(1.48\), is no better than least squares at \(1.50\); the \(t\) test rejects the truth \(23\%\) of the time. As the instruments strengthen, the bias falls roughly as \(1/F\): with \(F = 5.2\) the 2SLS bias is \(0.098\), about a quarter of the least squares bias \(0.418\), and with \(F = 26\) it is \(0.015\), and the test is close to its nominal size. The rule \(F > 10\) is the boundary between these regimes.

Stata: estat firststage reports the \(F\) and the Stock–Yogo critical values. EViews: View \(\to \) IV Diagnostics \(\to \) Weak Instrument Diagnostics.

Software 13.5.3 (Measurement error and the Hausman test). The survey data in R’s recommended package MASS record, for \(237\) statistics students, their height in centimetres and the spans of the writing and non-writing hands. Each span measures the same thing, hand size, with error. Regress height on the writing-hand span, using the other span as the instrument, in the spirit of Theorem 13.1.1:

> library(MASS)
> d <- na.omit(survey[, c("Height", "Wr.Hnd", "NW.Hnd", "Sex")]); n <- nrow(d)
> y <- d$Height; X <- cbind(1, d$Wr.Hnd)
> ols <- lm(Height ~ Wr.Hnd, data = d)
> iv1 <- tsls(y, X, cbind(1, d$NW.Hnd))
> c(n = n, F.first = summary(lm(Wr.Hnd ~ NW.Hnd, data = d))$fstatistic[[1]])
       n  F.first
 207.000 2637.647
> rbind(OLS = coef(summary(ols))[2, 1:2], IV = c(iv1$b[2], iv1$se[2]))
    Estimate Std. Error
OLS 3.121788  0.2895470
IV  3.147388  0.3005942
> H <- unname((iv1$b[2] - coef(ols)[2])^2 /
+             (iv1$se[2]^2 - vcov(ols)[2, 2]))  # Hausman
> v <- resid(lm(Wr.Hnd ~ NW.Hnd, data = d))  # control-function form
> c(H = H, p = pchisq(H, 1, lower.tail = FALSE),
+   t.v = coef(summary(lm(Height ~ Wr.Hnd + v, data = d)))["v", 3])
         H          p        t.v
 0.1005195  0.7512071 -0.3164351

Reading the output. With \(207\) complete records, the first-stage \(F\) is \(2638\): the two spans have correlation \(0.96\), so the instrument is very strong. The IV slope, \(3.15\) cm of height per cm of span, is barely above the least squares \(3.12\). The reliability of the writing-hand span, estimated by \(\Cov (X, Z)/\Var (X) = 1.00\), says that the hands were measured accurately, and Theorem 13.1.1 then predicts almost no attenuation. The Hausman statistic is \(0.10\) (\(p = 0.75\)), and the control-function \(t\) ratio is \(-0.316\), whose square is the same \(0.10\), as Theorem 13.4.2 says. Exogeneity is not rejected, and least squares, with its smaller standard error, is preferred. IV is a remedy for a problem that must first be shown to exist.

Stata: ivregress 2sls Height (Wr_Hnd = NW_Hnd), then estat endogenous. EViews: View \(\to \) IV Diagnostics \(\to \) Regressor Endogeneity Test.

Software 13.5.4 (An invalid instrument and the Sargan test). Add sex as a second instrument. It is certainly relevant, since men have larger hands, but it is not valid, because sex affects height directly, over and above hand size:

> Z2 <- cbind(1, d$NW.Hnd, d$Sex == "Male")   # add Sex as a second instrument
> iv2 <- tsls(y, X, Z2)
> S <- n * summary(lm(iv2$u ~ Z2[, -1]))$r.squared
> c(b = iv2$b[2], Sargan = S, p = pchisq(S, 1, lower.tail = FALSE))
           b       Sargan            p
3.172708e+00 4.716980e+01 6.509561e-12
> coef(summary(lm(Height ~ Wr.Hnd + Sex, data = d)))[, 1:2]
              Estimate Std. Error
(Intercept) 137.686951  5.7131129
Wr.Hnd        1.594446  0.3229395
SexMale       9.489814  1.2286579

Reading the output. The Sargan statistic is \(47.2\) on \(\ell - k = 1\) degree of freedom, \(p \approx 10^{-11}\): the two instruments disagree decisively. The last regression shows why: holding hand span fixed, men are \(9.5\) cm taller, so sex belongs in the equation, and the disturbance of the equation without it is correlated with sex. Here the test identifies the culprit because we know which instrument is the doubtful one; in general it does not. Including sex as a regressor and using the other hand’s span as the instrument would give a valid, exactly identified estimate.

Stata: estat overid after ivregress. EViews: the \(J\)-statistic is printed with every tsls output, and View \(\to \) IV Diagnostics \(\to \) Instrument Orthogonality Test tests subsets.

Questions on this section

Stuck on something here? Ask below and it stays attached to this topic.