Limited Information Maximum Likelihood (LIML) estimation in R? Curious whether anyone knows a package, or has written an implementation themselves, for conducting instrumental variables regressions using LIML in R. All of the R packages I have seen for IV regressions seem to resort to 2SLS/GMM as opposed to LIML, which may have more desirable finite sample qualities compared to 2SLS (Hayashi 2000)
For additional context, stata's ivregress command includes options to use LIML estimation, and hoping someone has already implemented something similar in R so I don't have to write it myself.
 A: There are a number of ways to motivate the LIML estimator, the primary one being that it is a member of the family of IV estimators known as the $k$-class estimators.
Also, to my knowledge, a direct implementation of the LIML estimator is not available in R (see the sem package though).
$k$-class IV estimators
Consider the linear instrumental variables regression model:
$$
\begin{align}
\boldsymbol{Y}_1 &= \mathbf{Y}_2\boldsymbol{\beta} + \mathbf{Z}_1\boldsymbol{\delta} + \boldsymbol{\varepsilon}\\
&= \mathbf{X}\boldsymbol{\gamma} + \boldsymbol{\varepsilon}\\
\mathbf{Y}_2 &= \mathbf{Z}_1\mathbf{\Pi}_1 + \mathbf{Z}_2\mathbf{\Pi}_2+  \mathbf{V}\\
&= \mathbf{Z}\mathbf{\Pi} +  \mathbf{V}\\
\mathbb{E}(\boldsymbol{\varepsilon} \mid \mathbf{Z}) &= \boldsymbol{0} 
\end{align}
$$
where $\mathbf{Z}_1$ is the matrix of included exogenous covariates, and $\mathbf{Z}_2$ is the matrix of excluded exogenous covariates. $\mathbf{Y}_2$ is the matrix of endogenous covariates.
Then the $k$-class estimators, defined by Theil, can be written as solutions to the set of equations
$$
\mathbf{X}'(\boldsymbol{i}_N - k\mathbf{M}_{\mathbf{Z}})\mathbf{X}\hat{\boldsymbol{\beta}}_{KC} = \mathbf{X}'(\boldsymbol{i}_N - k\mathbf{M}_{\mathbf{Z}})\boldsymbol{Y}
$$
where $\mathbf{M}_{\mathbf{Z}} = \boldsymbol{i}_N - \mathbf{Z}\left(\mathbf{Z}'\mathbf{Z}\right)^{-1}\mathbf{Z}$ is the orthogonal projection matrix. Various choices of $k$ yield special cases, for example, the 2SLS corresponds to $k=1$.
A very crude implementation of a $k$-class estimators in R is as follows:
library(foreign)
library(zoo)
library(AER)

#==========================================================
# load and pre-process the data
#==========================================================
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)

#==========================================================
# k-class estimator (example: k = 0.9)
#==========================================================
# 2SLS estimator; to check that the formula is okay
ivreg(log(ww) ~ ax. + I(ax.^2) +  we |
        ax. + I(ax.^2) + kl6 + k618 + wa, data = dfMroz, subset = ww > 0)
# OK

# IV regression formula
formulaMrozKC = as.Formula(log(ww) ~ ax. + I(ax.^2) +  we |
                             ax. + I(ax.^2) + kl6 + k618 + wa)

## get the model matrices
mfMrozKC = model.frame(formulaMrozKC, data = dfMroz, subset = ww > 0)
vY = model.response(mfMrozKC)
mX = model.matrix(formulaMrozKC, data = mfMrozKC, rhs = 1)
mZ = model.matrix(formulaMrozKC, data = mfMrozKC, rhs = 2)

mMZ = diag(428) - mZ %*% solve(t(mZ) %*% mZ) %*% t(mZ)

# k-class estimator (k = 0.9)
dK = 0.9
solve(a = t(mX) %*% (diag(428) - dK*mMZ) %*% mX,
            b =  t(mX) %*% (diag(428) - dK*mMZ) %*% vY, tol = 1e-10)

LIML estimator
The LIML estimator is had by setting $k$ in the $k$-class estimators to be the minimum eigenvalue of the matrix:
$$
\left(\mathbf{Y} \mathbf{M}_{\mathbf{Z}}\mathbf{Y}\right)^{-1/2}\mathbf{Y} \mathbf{M}_{\mathbf{Z}_1}\mathbf{Y} \left(\mathbf{Y} \mathbf{M}_{\mathbf{Z}}\mathbf{Y}\right)^{-1/2}
$$
where $\mathbf{Y} = [\boldsymbol{Y}_1, \mathbf{Y}_2]$. For a detailed discussion, see Davidson and MacKinnon (2001), pg. 539--.
R implementation
A very crude implementation of this in R is as follows (no guarantees about numerical efficiency are made):
#==========================================================
# LIML estimator (example: k = 0.9)
#==========================================================
# function to compute the inverse square root of a matrix
fnMatSqrtInverse = function(mA) {
  ei = eigen(mA)
  d = ei$values
      d = (d+abs(d))/2
      d2 = 1/sqrt(d)
      d2[d == 0] = 0
      return(ei$vectors %*% diag(d2) %*% t(ei$vectors))
}

mY2 = mX[, setdiff(colnames(mX), colnames(mZ))]
mZ1 = mZ[, intersect(colnames(mX), colnames(mZ))]
mMZ1 = diag(428) - mZ1 %*% solve(t(mZ1) %*% mZ1) %*% t(mZ1)
mYStar = cbind(vY, mY2)

mYZY = t(mYStar) %*% mMZ %*% mYStar
mLeftRight = fnMatSqrtInverse(mYZY)

dK.LIML = sort(eigen(b %*% (t(mYStar) %*% mMZ1 %*% mYStar) %*% b, 
                only.values = TRUE)$values)[1]

# LIML estimator
solve(a = t(mX) %*% (diag(428) - dK.LIML*mMZ) %*% mX,
      b =  t(mX) %*% (diag(428) - dK.LIML*mMZ) %*% vY, tol = 1e-10)

The computation of the matrix square root inverse is borrowed from here.
A: I made an implementation of the k-class in a package, currently hosted on github, RCompAngrist (dedicated to replicate examples in Angrist and Pishke's book). The function kclass() uses the same interface than ivreg, and let to pre-specify k (obtaining OLS, 2SLS, any estim) or estimating it. It seeks to use slightly more efficient computation, mainly through a QR decomposition, but is still under development.
With the actual development version, you can obtain the same results than above:
library(devtools)
install_github(repo = "RCompAngrist", username = "MatthieuStigler", subdir = "RcompAngrist")
library(RCompAngrist)

# run example of fg nu
download.file(url = 'http://people.stern.nyu.edu/wgreene/Text/tables/TableF4-1.txt',  destfile = 'mroz.txt')
dfMroz = read.table('mroz.txt', header = TRUE, skip = 36)
names(dfMroz) = tolower(names(dfMroz))

# do subset outside (currently arg subset within kclass does not work well)
dfMroz2 <- subset(dfMroz, ww>0)

# estimate 2SLS (k=1): same results than ivreg (different order though)
kclass(log(ww) ~ ax. + I(ax.^2) +  we |
    ax. + I(ax.^2) + kl6 + k618 + wa, data = dfMroz2,  k=1)

# estimate k-class with k= 0.9
k_0_9_kclass <- kclass(log(ww) ~ ax. + I(ax.^2) +  we |
    ax. + I(ax.^2) + kl6 + k618 + wa, data = dfMroz2,  k=0.9)

## Compare results with code of fg nu
k_0_9_manual <- solve(a = t(mX) %*% (diag(428) - dK*mMZ) %*% mX,
  b =  t(mX) %*% (diag(428) - dK*mMZ) %*% vY, tol = 1e-10)
all.equal(coef(k_0_9_kclass), k_0_9_manual[c(4,1,2,3),1], check.attributes = FALSE)

# obtain k value: leave arg k NULL:
k_0_9_kclass <- kclass(log(ww) ~ ax. + I(ax.^2) +  we |
                     ax. + I(ax.^2) + kl6 + k618 + wa, data = dfMroz2)

# same k than code of fg nu? Could not compare as missing b, to be done soon

A: The R package ivmodel implements the LIML estimator (link here). 
