# Generalized Method of Moments (GMM)
```{r}
#| include: false
library(MASS)
library(ivreg)
library(gmm)
library(sandwich)
```
## The Method of Moments (MOM)
A population moment $\gamma$ can be defined as the expectation of
some continuous function $g$ of a random variable $x$:
$$
\gamma={\mathrm{E}} [g(x)]
$$ {#eq-gmm-1}
On the other hand, a sample moment is the sample version of the
population moment in a particular sample:
$$
\hat \gamma=\frac{1}{n} \sum [g(x)]
$$ {#eq-gmm-2}
## OLS as a moment problem
Consider the simple linear regression
$$
{\symbf{y=X\beta+u, \quad u \sim IID(0, \sigma^2).
}}$$ {#eq-gmm-3}
If the model is correctly specified, then
$$
{\symbf{E (X'u)=0.
}}$$ {#eq-gmm-4}
The MOM principle suggests that we replace the left-hand side with
its sample analog $\frac{1}{n} {\symbf{X'(y-X\beta)}}$.
Since we know that the true ${\symbf{\beta}}$ sets the population moment
equal to zero in expectation, it seems reasonable to assume that a
good choice of ${\symbf{\hat \beta}}$ would be one that sets the sample
moment to zero. The MOM procedure suggests an estimate of ${\symbf{\beta}}$ that solves
$$
\frac{1}{n} {\symbf{X'(y-X \hat \beta)=0.
}}$$ {#eq-gmm-5}
The MOM estimator is
$$
{\symbf{\hat \beta=(X'X)^{-1}X'y,
}}$$ {#eq-gmm-6}
which is the same as the OLS estimator.
## IV as a moment problem
Consider the simple linear regression
$$
{\symbf{y=X\beta+u, \quad u \sim IID(0, \sigma^2).
}}$$ {#eq-gmm-7}
If one or more of the regressors is endogenous, then
$$
{\symbf{E (X'u)\neq 0.
}}$$ {#eq-gmm-8}
We have to find an instrumental variable ${\symbf{Z}}$ which is
$$
{\symbf{E (Z'u)= 0.
}}$$ {#eq-gmm-9}
Or,
$$
{\symbf{E (Z'(y-X\beta))= 0.
}}$$ {#eq-gmm-10}
The sample analogy of this is
$$
\frac{1}{n} {\symbf{Z'(y-X \hat \beta)=0.
}}$$ {#eq-gmm-11}
When the model is *exactly identified* -- the number of instruments
equals the number of regressors, so ${\symbf{Z'X}}$ is square and
nonsingular -- the sample moment condition has a unique solution, the
simple IV estimator
$$
{\symbf{\hat \beta=(Z'X)^{-1}Z'y.
}}$$ {#eq-gmm-12}
When the model is *overidentified* -- more instruments than regressors
-- ${\symbf{Z'X}}$ is no longer square and we cannot set all sample moments
to zero simultaneously. We instead minimize a quadratic form in the
moments, which gives the two-stage least squares (2SLS) estimator
$$
{\symbf{\hat \beta=(X'P_Z X)^{-1}X'P_Z y, \qquad P_Z = Z(Z'Z)^{-1}Z',
}}$$ {#eq-gmm-13}
the special case of the GMM estimator below with weighting matrix
$({\symbf{Z'Z)^{-1}}}$.
## The Generalized Method of Moments
The expectation ${{\mathrm{E}}}(Y^r)$ for any $r=1,2, \dots$ is called the
$r^{th}$ (raw) moment of $Y$. The expectation ${{\mathrm{E}}} [(Y-{{\mathrm{ E}}}(Y))^r]$ is called the $r^{th}$ centered moment of $Y$.
The mean is the first raw moment.
The variance is the second centered moment.
The third centered moment measures the skewness of the
distribution.
The fourth centered moment measures the kurtosis of the
distribution. Interpreted as a measure of "fatness of tails".
The standardized kurtosis is
$$k=\frac{E[(Y-E(Y))^4]}{E[(Y-E(Y))^2]^2}.$$ {#eq-gmm-14}
For a normal distribution, $k=3.$
For a $t$ distribution with $v > 4$ degrees of freedom,
$k=3+6/(v-4) > 3$, i.e., the $t$ distribution has fatter tails
than a normal distribution (the excess kurtosis $6/(v-4)$ shrinks toward 0 as $v$ grows).
The distribution function of a random variable captures all
information about the random variable. If the moment generating
function exists in a neighborhood of zero, the full set of moments
also determines the distribution (without such a condition this can
fail: the lognormal distribution is not determined by its moments).
This matters directly for GMM: GMM only ever matches a finite set of
moment conditions, and if the underlying distribution is not uniquely
pinned down by its (even infinite) sequence of moments, then no amount
of moment-matching -- however many moment conditions we add -- can
recover the full data-generating distribution. This is precisely the
sense in which GMM asks for less than MLE: MLE assumes we know the
full distributional family and estimates its parameters (recovering
the whole distribution when correctly specified), while GMM only
requires correctly specified moment conditions and only ever identifies
the parameters entering those moments, not the full distribution. When
the distributional assumption behind MLE is correct, MLE is efficient
and gives you more (the whole distribution); when it is not, GMM's
weaker requirements make it more robust.
This distinction underlies the relative strengths and weaknesses
of ML and GMM.
## GMM
The statistical model takes the general form
$$
E[m(Y_i; \theta_0)]=0
$$ {#eq-gmm-15}
where
- $Y_1, \cdots, Y_n$ are random variables from which the sample $y_1, \cdots, y_n$ is drawn,
- $m(Y, \theta)$ is a function specifying the model,
- $\theta_0$ is the "true value" of the parameter.
$E[m(Y_i; \theta_0)]=0$ are called the population moment
conditions.
Two ideas behind GMM:
1. Replace the population mean $E[.]$ with the sample mean
calculated from the observed sample $y_1, \cdots, y_n$.
2. Since $E[m(Y_i; \theta_0)]=0$, choose $\hat \theta_{GMM}$ to make
$\frac{1}{n}\sum_{i=1}^{n}m(y_i; \hat \theta_{GMM})$ as close to
zero as possible.
Define the notation
$$
\bar m(\theta)=\frac{1}{n} \sum_{i=1}^n m(y_i; \theta).
$$ {#eq-gmm-16}
$\hat \theta_{GMM}$ is chosen to make $\bar m(\theta)'\bar
m(\theta)$ as close to
zero as possible.
More generally, $\hat \theta_{GMM}$ is chosen to minimize $\bar
m(\theta)'W \bar m(\theta)$ for some weighting matrix $W$.
The choice of $W$ matters only under overidentification. With as many moment
conditions as parameters, all sample moments can be driven to zero exactly and
$W$ drops out -- which is why the just-identified IV estimator above did not
need one. When there are more moments than parameters, the efficient choice is
$W = S^{-1}$, where $S$ is the asymptotic variance of the sample moments,
$S = \lim_{n \to \infty} \mbox{Var} \bigl( \sqrt n \, \bar m (\theta_0) \bigr)$.
With that weighting,
$$
\sqrt n (\hat \theta_{GMM} - \theta_0) \xrightarrow{d} N \bigl( 0, (G'S^{-1}G)^{-1} \bigr),
\qquad G = {{\mathrm{E}}} \left[ \frac{\partial m (Y_i ; \theta_0)}{\partial \theta'} \right],
$$ {#eq-gmm-17}
which is the GMM counterpart of the information-matrix bound of the previous
chapter. Since $S$ depends on $\theta_0$, in practice one starts from any
consistent weighting (often the identity, or $({\symbf{Z'Z)^{-1}}}$ -- scalar
multiples of ${\symbf{W}}$ leave the minimizer unchanged, so with $\bar m$ defined as a
sample average this can equally be written $[n^{-1}{\symbf{Z'Z}}]^{-1}$), estimates $\hat
S$ from the resulting residuals, and re-minimizes: two-step GMM. Stopping at
$W = ({\symbf{Z'Z)^{-1}}}$ is exactly 2SLS, which is therefore efficient GMM only
under homoskedasticity.
### An example
Let's see an example with GMM, using the same simulated data as
before. We have the same situation as before, $X$ is endogenous. We
are doing GMM version of 2sls.
Here I use R's "gmm" library, which makes things easy. Its first two arguments
are `g` and `x`: with the formula interface, `g` is the model formula whose
residual forms the moment function, and `x` is a formula giving the instruments.
The moment condition in this example is $E({\symbf{z' u) = 0}}$, where ${\symbf{z}}$ is the
instrument -- the variable `w` in the simulated data -- and $u$ is the residual
from the structural equation. Note that ${\symbf{W}}$ above denotes the *weighting
matrix* in $\bar m' W \bar m$, which is a different object entirely, so the
instrument is written ${\symbf{z}}$ here.
```{r}
## DGP: data$y <- data$x + data$z + data$u
set.seed(66)
nobs=10000
nDim = 3
sdxx = 1
sdww=1
sdzz=1
## here we have three variables x,z,w.
## z is the omitted variable,x and z are correlated, w is the instrument, which is correlated with x, but not z. u is independent of everything else.
crxz=.6
crzw=0
crxw=.5
covarMat = matrix( c(sdxx^2, crxz, crxw, crxz, sdzz^2, crzw, crxw, crzw, sdww^2 ) , nrow=nDim , ncol=nDim )
covarMat
data = data.frame(mvrnorm(n=nobs, mu=rep(0,nDim), Sigma=covarMat ))
names(data) <- c('x','z','w')
data$u <- rnorm(nobs,0,1)
# dgp
data$y <- data$x + data$z + data$u
gmm.fit = gmm(y ~ x, x = ~w, data = data)
summary(gmm.fit)
# It returns the same estimates as the 2sls results.
tsls.model <- ivreg(y ~ x | w, data=data)
summary(tsls.model)
```
Two things about that agreement are worth being explicit about, because it is easy
to over-generalise from it. First, GMM and 2SLS coincide *numerically* here because
the model is just-identified: one instrument (`w`) for one endogenous regressor
(`x`), so ${\symbf{Z'X}}$ is square and the sample moments can be set to zero exactly,
whatever the weighting matrix. The output says so -- the $J$-test reports
"degrees of freedom is 0", and its value ($\sim 10^{-26}$) is numerical noise
rather than evidence. Second, `gmm()` is using an optimal two-step weighting with a
heteroskedasticity-and-autocorrelation-consistent kernel (note `Method: twoStep` and
`Kernel: Quadratic Spectral` in the output). Under *over*-identification that
choice would generally **not** reproduce 2SLS, for exactly the reason given above:
$W = ({\symbf{Z'Z)^{-1}}}$ is the efficient weighting only under homoskedasticity.
In OLS case, it would be $E(X u) = 0$.
```{r}
gmm.ols = gmm(y ~ x, x = ~x, data = data)
summary(gmm.ols)
# It returns the same estimates as the OLS results.
ols <- lm(y ~ x, data=data)
summary(ols)
```
### Overidentification, the weighting matrix, and the $J$ test
The example above cannot show what GMM adds, because it is just-identified. One
instrument for one endogenous regressor leaves nothing for the weighting matrix to
do and nothing for a specification test to check. Two things have to change. We
need more instruments than endogenous regressors, and we need errors that are not
homoskedastic -- under homoskedasticity $({\symbf{Z'Z)^{-1}}}$ already is the efficient
weighting, so 2SLS and efficient GMM agree and the comparison is empty.
The design keeps $y = x + z + u$ with $z$ omitted from the regression, so the
composite error is $z+u$ and $x$ is endogenous through its correlation with $z$.
There are now two instruments. Both $w_1$ and $w_2$ correlate $0.5$ with $x$ and
$0.2$ with each other, and neither correlates with $z$, so both are valid. The
error is drawn with standard deviation $\exp(0.8 w_1)$, which leaves
$E[u \mid w]=0$ intact while making the variance depend on an instrument. That is
exactly the case where the efficient weighting differs from
$({\symbf{Z'Z)^{-1}}}$. The true coefficient on $x$ is $1$.
```{r}
set.seed(66)
nobs <- 10000
# cor(x,z)=.6; cor(x,w1)=cor(x,w2)=.5; cor(z,w1)=cor(z,w2)=0; cor(w1,w2)=.2
covarMat <- matrix(c(1, .6, .5, .5,
.6, 1, 0, 0,
.5, 0, 1, .2,
.5, 0, .2, 1), nrow = 4)
d <- data.frame(mvrnorm(n = nobs, mu = rep(0, 4), Sigma = covarMat))
names(d) <- c("x", "z", "w1", "w2")
# the error variance depends on w1, so E[u|w] = 0 still holds but 2SLS is
# no longer the efficient GMM estimator
d$u <- rnorm(nobs, 0, exp(0.8 * d$w1))
d$y <- d$x + d$z + d$u
fit_ident <- gmm(y ~ x, x = ~ w1 + w2, data = d, wmatrix = "ident")
fit_eff <- gmm(y ~ x, x = ~ w1 + w2, data = d)
tsls <- ivreg(y ~ x | w1 + w2, data = d)
data.frame(
estimator = c("GMM, identity weighting", "GMM, efficient two-step",
"2SLS, classical SE", "2SLS, robust SE"),
b_x = round(c(coef(fit_ident)["x"], coef(fit_eff)["x"],
coef(tsls)["x"], coef(tsls)["x"]), 4),
se = round(c(sqrt(diag(vcov(fit_ident)))["x"], sqrt(diag(vcov(fit_eff)))["x"],
sqrt(diag(vcov(tsls)))["x"],
sqrt(diag(vcovHC(tsls, type = "HC0")))["x"]), 4),
row.names = NULL)
specTest(fit_eff)
```
Three readings. The weighting matrix now matters: identity weighting gives
$1.015$ and the efficient two-step estimator gives $0.989$, where in the
just-identified case every weighting gave the same number. The efficiency gain is
real but modest, and it shows up only against the *robust* 2SLS standard error,
$0.042$ against $0.046$. And the classical 2SLS standard error, $0.033$, is far
too small -- the earlier remark about homoskedasticity, arriving as a number
instead of a caveat.
The $J$ statistic is now a test rather than a formality. It is a joint specification
test of the moment conditions, not a certificate of instrument validity. It has one degree of
freedom -- three moment conditions, from the constant and the two instruments,
less the two parameters -- and its value of $1.65$ with a $p$-value of $0.20$
gives no reason to doubt the moment conditions.
To see it earn its keep, make $w_2$ invalid by correlating it with the omitted
$z$. Nothing else changes.
```{r}
covarMat2 <- covarMat
covarMat2[2, 4] <- covarMat2[4, 2] <- 0.4 # w2 now correlates with omitted z
set.seed(66)
d2 <- data.frame(mvrnorm(n = nobs, mu = rep(0, 4), Sigma = covarMat2))
names(d2) <- c("x", "z", "w1", "w2")
d2$u <- rnorm(nobs, 0, exp(0.8 * d2$w1))
d2$y <- d2$x + d2$z + d2$u
fit_bad <- gmm(y ~ x, x = ~ w1 + w2, data = d2)
coef(fit_bad)["x"]
specTest(fit_bad)
```
The coefficient moves to $1.60$ and the $J$ statistic goes to $143$, with a
$p$-value of $5 \times 10^{-33}$. This is what overidentification buys. With one
instrument there is no way to notice that it is invalid, because a single moment
condition can always be satisfied exactly. With two, the instruments have to
agree about $\beta$, and when one of them is contaminated they do not. Two limits
are worth keeping in view. The test does not say *which* instrument is at fault,
and it has no power against violations the instruments share: if $w_1$ and $w_2$
were both correlated with $z$ in the same way they would agree on the same wrong
answer, and $J$ would pass.
## A few concepts of conditioning
### Independence
If $X$ and $Y$ are independent then
$$
f(x,y)=f(x)f(y)
$$ {#eq-gmm-18}
and hence
$$
f(y|x)=f(y).
$$ {#eq-gmm-19}
If $X$ and $Y$ are independent then
$$
E[g(X)h(Y)]=E[g(X)]\cdot E[h(Y)]
$$ {#eq-gmm-20}
and hence
$$
Cov[g(X),h(Y)]=0.
$$ {#eq-gmm-21}
i.e. all functions of $X$ and $Y$ are uncorrelated.
### Law of Iterated Expectations
$$
E[Y]=E[E(Y|X)].
$$ {#eq-gmm-22}
### Dependence Concepts
$X$, $Y$ independent:
$$
Cov[g(X),h(Y)]=0
$$ {#eq-gmm-23}
$X$, $Y$ uncorrelated:
$$
Cov[X,Y]=0
$$ {#eq-gmm-24}
$E[Y|X]=0$:
$$
Cov[g(X),Y]=0
$$ {#eq-gmm-25}
These are listed in decreasing strength, and the implications run one way only:
independence $\Rightarrow$ mean independence $\Rightarrow$ uncorrelatedness. The
middle condition is what a *model of the conditional mean* needs. Neither converse holds: $X$ and
$Y$ can be uncorrelated with $E[Y \mid X] \neq 0$ (take $Y = X^2 - 1$ with
$X \sim N(0,1)$, so $Cov[X,Y] = E[X^3] = 0$ while $E[Y \mid X] = X^2 - 1$), and
mean independence permits dependence through higher moments (any pure
heteroskedasticity, $Y = X\epsilon$ with $\epsilon$ independent of $X$ and mean
zero). Keeping the three apart matters because the moment conditions below need only the
*weakest* of them. $E[u_i]=0$ and $E[X_i u_i]=0$ are statements about means and
covariances, nothing more. Conditional mean zero is sufficient for them and is not
necessary, which is exactly why GMM asks for less than a conditional-mean model
does.
### Regression
A regression model is a model of $E[Y_i|X_i]$. For example,
$$
Y_i=\beta_0+\beta_1X_i+u_i
$$ {#eq-gmm-26}
where $E[u_i|X_i]=0$.
### GMM regression
The regression model
$$
Y_i=\beta_0+\beta_1X_i+u_i, \quad E[u_i|X_i]=0
$$ {#eq-gmm-27}
implies the moment condition
$$
E[u_i]=0 \quad \mbox{and} \quad E[X_i u_i]=0
$$ {#eq-gmm-28}
That is,
$$
E[Y_i-\beta_0-\beta_1X_i]=0
$$ {#eq-gmm-29}
$$
E[X_i(Y_i-\beta_0-\beta_1X_i)]=0
$$ {#eq-gmm-30}
The sample moment conditions are
$$
\frac{1}{n}\sum_{i=1}^{n}(y_i-\hat \beta_0-\hat \beta_1x_i)=0
$$ {#eq-gmm-31}
$$
\frac{1}{n}\sum_{i=1}^{n}x_i(y_i-\hat \beta_0-\hat \beta_1x_i)=0
$$ {#eq-gmm-32}
These are just normal equations for OLS.
A characteristic of GMM: the specification of the model generates
the estimator. i.e. only $E[Y_i|X_i]=\beta_0+\beta_1 X_i$ is
assumed.
Note there are no assumptions that $u_i$ is homoscedastic, not
autocorrelated or normally distributed. These properties affect
the statistical properties of the GMM estimator, not its
definition.