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.
