IV estimator: efficient implementation? I would like to implement (in R) an instrumental variable (IV) estimator, that takes the most general form (here not 2SLS or GMM!):
$$ \beta_{IV} = (Z'X)^{-1}Z'Y $$
I could code this in the naive way, inverting $Z'X$, and computing the rest, but I remember that this is very inefficient in term of computing, and that in the similar case of OLS $(X'X)^{-1}X'Y$, one should use rather a QR decomposition rather than inverting $X'X$. 
However, in the case of my IV estimator, I guess the standard QR won't apply, so I am not sure whether there is any similar technique/decomposition that could be applied in this case? 
Any idea of an efficient technique to compute $(Z'X)^{-1}Z'Y$ without inverting $(Z'X)^{-1}$? 
Thanks!
 A: OLS and IV estimators are both solutions to systems of equations, which can be computed using QR decompositions. 
Linear regression model
The linear regression model is
$$
\begin{align}
\boldsymbol{Y} &= \mathbf{X}\boldsymbol{\beta} + \boldsymbol{\varepsilon}\\
\mathbb{E}(\boldsymbol{\varepsilon} \mid \mathbf{X}) &= \boldsymbol{0}
\end{align}
$$
The OLS estimator ($\hat{\boldsymbol{\beta}}_{OLS}$) is the solution to the system of equations:
$$
\left(\mathbf{X}'\mathbf{X}\right)\hat{\boldsymbol{\beta}}_{OLS}  = \mathbf{X}'\boldsymbol{Y}
$$
This computation can be easily performed using a call to qr.solve. I have added formula handling because that makes the interface neat.
#==========================================================
# download and pre-process the data
#==========================================================
library(AER)
library(sandwich)
library(Formula)

# download the Mroz dataset
download.file(url = 'http://people.stern.nyu.edu/wgreene/Text/tables/TableF4-1.txt', 
              destfile = 'mroz.txt')

# read in the file
dfMroz = read.table('mroz.txt', header = TRUE, skip = 36)
names(dfMroz) = tolower(names(dfMroz))
summary(dfMroz)

#==========================================================
# linear regression of wife's wages
#==========================================================
# using lm()
lmMroz = lm(log(ww) ~ we + ax. + I(ax.^2), 
             data = dfMroz[dfMroz$ww > 0, ])
summary(lmMroz)

# computation by hand
## define the formula on which the OLS estimate is based
formulaMrozLM = log(ww) ~ we + ax. + I(ax.^2)

## get the model matrices
mfMrozLM = model.frame(formulaMrozLM, data = dfMroz, subset = ww > 0)
vYMrozLM = model.response(mfMrozLM)
mXMrozLM = model.matrix(formulaMrozLM, data = mfMrozLM)

# solve using solve() and qr.solve()
# solve(crossprod(mXMrozLM), crossprod(mXMrozLM, vYMrozLM))
vBetaLMHat = qr.solve(mXMrozLM, vYMrozLM)

Instrumental variables model
The instrumental variables model is
$$
\begin{align}
\boldsymbol{Y} &= \mathbf{X}\boldsymbol{\beta} + \boldsymbol{\varepsilon}\\
\mathbf{X} &= \mathbf{Z}\mathbf{\Pi} +  \mathbf{V}\\
\mathbb{E}(\boldsymbol{\varepsilon} \mid \mathbf{Z}) &= \boldsymbol{0} 
\end{align}
$$
where I have included the instrument exogeneity condition, but excluded the relevance condition.
The IV estimator for the case that there are as many instruments as there are endogenous variables, is defined as the solution to the system of equations
$$
\left(\mathbf{Z}'\mathbf{X}\right)\hat{\boldsymbol{\beta}}_{IV} = \mathbf{Z}'\boldsymbol{Y}
$$
In the general case, where there are more instrumental variables than endogenous variables in the system, we have the 2SLS estimator which is an IV estimator, for $\mathbf{Z} = \hat{\mathbf{X}}$, where $\hat{\mathbf{X}}$ is the linear projection of the matrix $\mathbf{X}$ onto the columns of the matrix $\mathbf{Z}$, that is, $\hat{\mathbf{X}} = \mathbf{Z}\hat{\mathbf{\Pi}}$. This projection is computed by way of a solution of the auxiliary system of equations:
$$
\left(\mathbf{Z}'\mathbf{Z}\right)\hat{\mathbf{\Pi}} = \mathbf{Z}'\mathbf{X}
$$ 
#==========================================================
# intrumental variables regression of wife's wages
#==========================================================
# computing using ivreg() in AER
ivMroz = ivreg(log(ww) ~ ax. + I(ax.^2) + kl6 + k618 + we |
                 ax. + I(ax.^2) + kl6 + k618 + wmed + wfed + he, 
               data = dfMroz[dfMroz$ww > 0, ])
summary(ivMroz)

# NOTE this is provided by the 'Formula' package of Achim Zeileis
#   and this is used for brevity only, and these matrices can 
#   easily created by hand
formulaMrozIV = as.Formula(log(ww) ~ ax. + I(ax.^2) + kl6 + k618 + we |
                        ax. + I(ax.^2) + kl6 + k618 + wmed + wfed + he)
mfMrozIV = model.frame(formulaMrozIV, data = dfMroz, subset = ww > 0)
vYMrozIV = model.response(mfMrozIV)
mXMrozIV = model.matrix(formulaMrozIV, data = mfMrozIV, rhs = 1)
mZMrozIV = model.matrix(formulaMrozIV, data = mfMrozIV, rhs = 2)

# compute the mXHat
# NOTE in the case of the pure IV estimator, where the number 
#   of endogenous variables is the same as the number of exogenous
#   variables, this computation is unnecessary
mXHat = mZMrozIV %*% qr.solve(mZMrozIV, mXMrozIV)

# compute the 2SLS (IV) estimator
vBetaIVHat = qr.solve(mXHat, vYMrozIV)

A: The IV estimator in your equation seems hard to understand. The length of the instrument vector $Z$ may not equal that of the regressor vector $X$. If they are not equal, the matrix $Z'X$ is not a square matrix; its inversion would be ill defined. 
One possible alternative is to introduce a positive-definite weighting matrix $S$ such that $\beta_{IV} = ((X'Z)S(Z'X))^{-1} (X'Z)S(Z'y)$. This returns to the GMM framework. The weighting matrix $S$ could be an identity matrix or a more efficient alternative. In this case, the inversion of the matrix could be done through the usual QR decomposition. I would suggest the textbook "Econometrics" by Fumio Hayashi, which has a more detailed description of the IV estimators. Good luck!
