# Applying Frisch-Waugh-Lovell theorem to IV regression in R

I am estimating an instrumental variables linear regression that has a large number of indicator (factor) variables. I don't particularly care about the coefficient estimates on those indicator variables. In Stata's ivreg2 package there is a "partial" option that applies the Frisch-Waugh-Lovell theorem to orthogonalize the dependent and exogenous variables to the indicator variables. After this transformation the indicator variables are not estimated because they do not affect the coefficients on the variables I am interested in.

My question is, is there something like this in R? It doesn't have to be part of an IV regression package but I am looking to orthogonalize one set of variables to another set of variables. This seems like something that would have already been implemented. Thanks.

-

It is possible to implement this yourself, as long as you are able to keep the notation straight. I will work through this systematically by first laying out the notation, and then translating that notation transparently into R code. Due to the length of this note, bugs might exist. Additionally, the reference to what you have described is Baum, Schaffer and Stillman (2007), pg 484--.

## Instrumental variables model

The instrumental variables models is \begin{align} \boldsymbol{Y} &= \mathbf{X}_1\boldsymbol{\beta}_1 + \mathbf{X}_2\boldsymbol{\beta}_2 + \boldsymbol{\epsilon} \\ \mathbf{X}_1 &= \mathbf{Z}_1\mathbf{\Pi}_1 + \mathbf{Z}_2\mathbf{\Pi}_2 + \mathbf{V} \end{align} where $\mathbf{X}_1$ is the matrix of endogenous variables, $\mathbf{X}_2$ is the matrix of included exogenous variables, $\mathbf{Z}_1$ is the matrix of excluded exogenous variables ($\mathbf{Z}_2$ and $\mathbf{X}_2$ are identical).

## Partial projections

Now suppose we want to partial out the effect of some of the included exogenous variables, $\mathbf{X}_{21}$, such that $\mathbf{X}_2 = [\mathbf{X}_{21} \, \mathbf{X}_{22}]$, and are only interested in the effects $\boldsymbol{\beta}_{22}$ and $\boldsymbol{\beta}_1$. Then, define the orthogonal projection matrix, $$\mathbf{M}_{\mathbf{X}_{21}} = \boldsymbol{\iota}_N - \mathbf{X}_{21}\left(\mathbf{X}_{21}'\mathbf{X}_{21}\right)^{-1}\mathbf{X}_{21}'$$.

Now define the following set of partialled out matrices: \begin{align} \boldsymbol{Y}^P &= \mathbf{M}_{\mathbf{X}_{21}}\boldsymbol{Y} \\ \mathbf{X}_1^P &= \mathbf{M}_{\mathbf{X}_{21}}\mathbf{X}_1 \\ \mathbf{X}_{22}^P &= \mathbf{M}_{\mathbf{X}_{21}}\mathbf{X}_{22} \\ \mathbf{Z}_1^P &= \mathbf{M}_{\mathbf{X}_{21}}\mathbf{Z}_1 \end{align}

where $\cdot^P$ indicates the partialled (with respect to $\mathbf{X}_{21}$) versions of the respective matrices.

The new instrumental variables regression model is \begin{align} \boldsymbol{Y}^P &= \mathbf{X}_1^P\boldsymbol{\beta}_1 + \mathbf{X}_{22}^P\boldsymbol{\beta}_{22} + \boldsymbol{\epsilon}^P \\ \mathbf{X}_1^P &= \mathbf{Z}_1^P\mathbf{\Pi}_1 + \mathbf{Z}_{22}^P\mathbf{\Pi}_{22} + \mathbf{V}^P \end{align} And you can see the equivalence from the original model just by pre-multiplying through by $\mathbf{M}_{\mathbf{X}_{21}}$.

## R implementation

### Simulating an IV model

I begin by simulating the IV model that I have specified.

rm(list = ls())
set.seed(101)

#==========================================================
# simulate an IV model
#==========================================================
iK1 = 3 # number of endogenous variables
iK21 = 4 # number of included exogenous variables (partialled)
iK22 = 2 # number of included exogenous variables (not partialled)
## added one due to constant
iL = 4 # number of excluded exogenous variables

iN = 1000 # sample size

# first stage coefficients
mP1 = matrix(rnorm(iL*iK1), nrow = iL, ncol = iK1)
mP21 = matrix(rnorm(iK21*iK1), nrow = iK21, ncol = iK1)
mP22 = matrix(rnorm((iK22+1)*iK1), nrow = iK22 + 1, ncol = iK1)
# one extra dimension added due to the coefficient

# first stage matrices
mZ1 = matrix(rnorm(iL*iN), nrow = iN, ncol = iL)
mZ21 = matrix(rnorm(iK21*iN), nrow = iN, ncol = iK21)
mZ22 = cbind(matrix(rnorm(iK22*iN), nrow = iN, ncol = iK22))

# make assignments for consistent notation
mX21 = mZ21
mX22 = mZ22

# reduced form coefficient vectors
vBeta1 = rnorm(iK1)
vBeta21 = rnorm(iK21)
vBeta22 = rnorm(iK22)
dAlpha = rnorm(1)

# generate correlated draws from the MVN
library(MASS)
mCovar = matrix(sample(
c(1, 2, 3, 4,
1.5, 2.5, 3.5, 4.5,
.75, 1.75, 2.75, 3.75,
0.33, 0.66, 1.33, 2.33)),
nrow = 4, ncol = 4)
mCovar = crossprod(mCovar)

# correlated errors between the first stage and the reduced form equations
mVE = mvrnorm(iN, mu = rep(0, iK1 + 1), Sigma = mCovar)

# generate the endogenous variables
mX1 = mZ1 %*% mP1 + mZ21 %*% mP21 + cbind(mZ22, 1) %*% mP22 +
mVE[, -1]  # first stage error matrix

# generate the outcome variable
vY = dAlpha + mX1 %*% vBeta1 + mX21 %*% vBeta21 + cbind(mX22) %*% vBeta22 +
mVE[, 1]

# put it all together in a data frame
dfIVSim = data.frame(vY, mX21, mX22, mZ1, mX1)
names(dfIVSim) = c('y',
paste0('incExogP', 1:iK21),
paste0('incExogNP', 1:iK22),
paste0('excExog', 1:iL),
paste0('endog', 1:iK1))


### IV estimation without partialling

The next step is to estimate an instrumental varaibles models for the whole system without partialling out any of the exogenous regressors. This is easily done using the ivreg() function in the AER package.

#==========================================================
# estimate the IV model without partialling
#==========================================================
library(AER)
ivNonPartial = ivreg(as.formula(paste0('y ~',
paste0(
c(grep('incExogNP', names(dfIVSim), value = TRUE),
grep('incExogP', names(dfIVSim), value = TRUE),
grep('endog', names(dfIVSim), value = TRUE)),
collapse = '+'
), '|',
paste0(
c(grep('incExogNP', names(dfIVSim), value = TRUE),
grep('incExogP', names(dfIVSim), value = TRUE),
grep('excExog', names(dfIVSim), value = TRUE)),
collapse = '+'
),
sep = '')),
data = dfIVSim)


### IV estimation with partialling

The last step is to perform the same computations but after partialling out the matrix $\mathbf{X}_{21}$.

I will do these computations by hand, and if you are not familiar with the matrix formulae for constructing the 2SLS estimator, please see my answer here.

#==========================================================
# estimate the IV model after partialling out mX21
#==========================================================
# orthogonal projection matrix
mMX21 = diag(iN) - mX21 %*% solve(crossprod(mX21)) %*% t(mX21)

# matrix projections
vYP = mMX21 %*% vY
mX22P = mMX21 %*% cbind(mX22, 1)  # project intercept as well
mZ1P = mMX21 %*% mZ1
mX1P = mMX21 %*% mX1

# the instruments matrix
mZP = cbind(mZ1P, mX22P)
colnames(mZP) = c(paste0('excExog', 1:iL),
paste0('incExogNP', 1:iK22), '(Inercept)')

# the covariates matrix
mXP = cbind(mX1P, mX22P)
colnames(mXP) = c(paste0('endog', 1:iK1),
paste0('incExogNP', 1:iK22), '(Intercept)')

# projection of the covariates on the instruments
mXPHat = mZP %*% qr.solve(mZP, mXP)

# compute the 2SLS (IV) estimator
vBetaHat = qr.solve(mXPHat, vYP)


### Combined results

Once that is done, we can finally put everything together and see that the Frisch-Waugh-Lovell theorem does indeed hold in the case of the 2SLS estimator of the IV regression model.

# print out the combined results
mBeta = cbind(vBetaHat, coef(ivNonPartial)[
match(rownames(vBetaHat), names(coef(ivNonPartial)))], c(vBeta1, vBeta22, dAlpha))
colnames(mBeta) = c('Non-partial 2SLS', 'Partial 2SLS', 'True values')
print(mBeta)


This yields:

> print(mBeta)
Non-partial 2SLS Partial 2SLS True values
endog1            -0.4084523   -0.4084523  -0.2287651
endog2            -0.5572848   -0.5572848  -0.8119480
endog3            -0.7565501   -0.7565501  -0.5590752
incExogNP1         1.7839968    1.7839968   1.3362787
incExogNP2        -1.3708189   -1.3708189  -0.7752127
(Intercept)        0.7152798    0.7152798   0.9934898
`
-