[Edited for clarity and detail]
Summary: A random effects model should produced biased estimates for a dynamic panel, but the lmer function of the R lme4 package produces unbiased estimates and I don't know why.
I am comparing the results of the R packages plm and lme4 estimating a dynamic panel of the form
$$y_{it} = \alpha_i + \theta y_{it-1} + x_{it}\beta + \epsilon_{it}$$
with the $\alpha_i$ being unobserved effects. Under the usual assumptions (see Wooldridge 2002 p. 677 for instance), it is well known that this model cannot be estimated by traditional random and fixed effects.
A random effect estimation would consider the unobserved individual parameter $\alpha_i$ as part of the unobserved error term, $\nu_{it} = \alpha_i + \epsilon_{it}$. Because $\alpha_i$ influences $y_{it}$ for every $t$, we have $\text{Cor}(y_{t-1},\nu) \neq 0$ which leads to biased estimates.
A fixed effect estimation can avoid this problem by eliminating the individual effect by, say, differencing. However, it is also biased because the lagged first-difference will end being correlated with the lagged error. (Nickel 1981).
There are several approaches to this problem, but the traditional solutions involve using lags as instruments (Arellano Bond estimator) and there is a rich literature proposing variations of this idea.
But lmer does the job, apparently.
If I have a simulated dataset with two predictors x1 and x2 and the units are indexed by the variable id, this model would be estimated in lme4 as
lmer(y~ lag_y + x1 + x2 + (1|id), data = my_data)
and in plm, an excellent panel data package, as
plm(y ~ lag_y + x1 + x2, index = "id", model = "random", data = my_data)
As one comment pointed out, there is not one way of estimating random effects, but several ("amemiya", "nerlove", etc.). None of these methods recovers the parameters of the dynamic model without bias.
Let's simulate. I will create 100 datasets according to
$$y_{it} = a_i + 0.7 y_{it-1} + 0.5 x_{1,it} - x_{2, it} + e_{it}$$
with $x_1$ and $x_2$ drawn from uniform (-1/2, 1/2). $e$ is a standard normal disturbance. The $a_i$ can drawn be from a normal, but in this example I sampled from {-2, -1.5, ..., 1.5, 2} and I have found that other distributions make no difference. I set $T = 10$ and the number of observed units equal to 50, making a dataset of size $n=500$
The results for 100 simulations with true parameter 0.7 are in the table below. You can find the average across the 100 simulations and the 2.5% and 97.5% percentiles. Notice that one row says "IV" in which I use a random effect model with lagged dependent variables and lagged x as IV.
| Method | average est. (true = 0.7) | 0.025 q | 0.975 q |
|---|---|---|---|
FD fixest |
0.612 | 0.549 | 0.671 |
Amemiya plm |
0.676 | 0.614 | 0.735 |
Walhus plm |
1.001 | 0 968 | 1.027 |
Nerlove plm |
0.667 | 0.606 | 0.724 |
Swar plm |
1.016 | 0.969 | 1.040 |
lags as "IV" plm |
0.662 | 0.594 | 0.731 |
Mixed lme4 |
0.699 | 0.632 | 0.763 |
While some methods "cover" the true value, the mixed model is the only one that is always centered at the true value.
So the question is, why is the mixed model so good at estimating this model? Some of the other random effects models are clearly producing biased estimates.
Can someone explain to me or guide me to a literature that explains why the mixed model works well in dynamic panels?
Here is reproducible code if you would like to try it:
library(lme4)
library(plm)
gen_one_unit = function(TT, theta = 0.7, beta = c(0.5, -1), s = 1){
# This function generates observations for one "unit" of the panel
x1 = runif(TT) - .5
x2 = runif(TT) - .5
alpha_i = sample(seq(-2, 2, by = .5), size = 1) # arbitrary
y = rep(1, TT) # initialize
for(t in 2:TT){
y[t] = alpha_i + theta*y[t-1] +
beta[1]*x1[t] + beta[2]*x2[t] + rnorm(1, 0, s)
}
ind_df = data.frame(alpha_i = rep(alpha_i, TT), x1, x2,
y, lagy = c(NA,y[1:(TT-1)]))
return(ind_df)
}
H = 50 # number of units
TT = 20 # number o time periods
dyn_linear_list = list()
for(i in 1:H){
temp_df = gen_one_unit(TT = TT, theta= .5)
temp_df$id = rep(i, TT)
dyn_linear_list[[i]] = temp_df
}
dyn_df = do.call(rbind, dyn_linear_list)
dyn_df = dyn_df[!is.na(dyn_df$lagy),] # get rid of the first obs
reg_plm = plm(y ~ x1 + x2 + lagy,
model = "random", index = "id",
# random.method = "amemiya"
# random.method = "nerlove", etc.
data = dyn_df)
reg_lme = lmer(y ~ x1 + x2 + lagy + (1|id),
data = dyn_df)
summary(reg_plm) # look at the coefficient of lagy
summary(reg_lme) # look at the coefficient of lagy
References
Nickell, S. (1981). Biases in dynamic models with fixed effects. Econometrica: Journal of the econometric society, 1417-1426.
