# Joint significance test in 2SLS (instrumental variables) regression

I'm trying to implement a joint test of the two coefficients comprising a quadratic term in a 2-stage least squares regression. The quadratic term is endogenous. I'm using AER in R, and ivreg's anova method is not giving me the same result as the manual Wald test that I'm checking it with. I'd basically like to know whether my own manual method is correct or not. If not, why not, and if so, what AER is doing differently.

1> rm(list=ls())
1> set.seed(1)
1> N <- 100
1> z <- rnorm(N) #The instrument
1> u <- rnorm(N) #The error term
1> x <- 1 + z - .1*z^2 + u + rnorm(N) # x is correlated with the error term u (endogeneity) and the instrument z
1> ex <- 1 + rnorm(N) #an exogenous variable
1> y <- 1 + x-.1*x^2  + ex + u
1> x2 <- x^2
1> z2 <- z^2
1> summary(lm(y~x+x2+ex))

Call:
lm(formula = y ~ x + x2 + ex)

Residuals:
Min       1Q   Median       3Q      Max
-1.91064 -0.57302  0.04697  0.43678  1.62413

Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept)  0.51137    0.13751   3.719 0.000337 ***
x            1.22946    0.05578  22.042  < 2e-16 ***
x2          -0.03512    0.02145  -1.637 0.104893
ex           0.97401    0.07933  12.278  < 2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.7699 on 96 degrees of freedom
Multiple R-squared:  0.8952,    Adjusted R-squared:  0.892
F-statistic: 273.5 on 3 and 96 DF,  p-value: < 2.2e-16


BIAS

 1> library(AER)

1> miv = ivreg(y~x+x2+ex|z+z2+ex)
1> summary(miv)

Call:
ivreg(formula = y ~ x + x2 + ex | z + z2 + ex)

Residuals:
Min      1Q  Median      3Q     Max
-2.1212 -0.7533 -0.2623  0.6458  3.5927

Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept)   1.4767     0.4785   3.086  0.00265 **
x             1.0678     0.1351   7.902 4.58e-12 ***
x2           -0.2185     0.1014  -2.154  0.03373 *
ex            0.8690     0.1409   6.167 1.65e-08 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 1.245 on 96 degrees of freedom
Multiple R-Squared: 0.7259, Adjusted R-squared: 0.7173
Wald test: 43.79 on 3 and 96 DF,  p-value: < 2.2e-16


MUCH BETTER WITH THE VALID INSTRUMENTS

Wald test using built-in method, versus manual Wald test:

1> mnull = ivreg(y~ex)
1> anova(miv,mnull,vcov=vcov(miv))
Wald test

Model 1: y ~ x + x2 + ex | z + z2 + ex
Model 2: y ~ ex
Res.Df Df      F    Pr(>F)
1     96
2     98 -2 31.512 3.011e-11 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

1> #Wald test
1> B = coef(miv)
1> r = matrix(c(0,0),2)
1> R = t(matrix(c(   0,1,0,0
1+      ,0,0,1,0),4))
1> V = vcov(miv)
1> W = t(R%*%B -r) %*% solve(R%*%V%*%t(R)) %*% (R%*%B -r)
1> W/2
[,1]
[1,] 31.5123


Why is my wald statistic twice the one given by the method?

    1> 1-pf(W/2,nrow(R),N-nrow(R))
[,1]
[1,] 2.706191e-11


And why are the p-values similar, but not identical?

The $W$ you're manually calculating is a $\chi^2$, not an $F$. Thus 1-pchisq(W,nrow(R)) would be the appropriate command.

In addition, note that $W \sim\chi^2$ with $q$ degrees of freedom (# of joint hypothesis) and ${1 \over q} W \sim F$ with $(q,d)$ degrees of freedom ($d$ being model degrees of freedom).

That's why you noticed your statistic being twice the anova $F$ statistic.

And you have a mistake in your manual calculation 1-pf(W/2,nrow(R),N-nrow(R)). That should have 96 degrees of freedom instead of 98 (from the full model, not the reduced model), in which case you match anova:

> 1-pf(W/2,nrow(R),miv\$df.residual)
[,1]
[1,] 3.011147e-11


And the chi-square version of the test:

> anova(miv,mnull,vcov=vcov(miv),test='Chisq')
Wald test

Model 1: y ~ x + x2 + ex | z + z2 + ex
Model 2: y ~ ex
Res.Df Df  Chisq Pr(>Chisq)
1     96
2     98 -2 63.025  2.062e-14 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
> 1-pchisq(W,nrow(R))
[,1]
[1,] 2.065015e-14