Suppose we have a panel-data model with an individual effect \(c_i\):
\[y_{it} = X_{it}\beta + c_i + u_{it}\]
The fixed-effects estimator removes \(c_i\) by demeaning, so \(c_i\) can be correlated with \(X_{it}\). But the same transformation also removes time-invariant variables such as gender, race, or baseline schooling.
Mundlak (1978) instead writes \(c_i\) as a linear projection on the unit means of the time-varying covariates:
We can estimate this equation by pooled OLS or random effects. By the Frisch–Waugh–Lovell theorem, the coefficient on \(X_{it}\) is the same as the fixed-effects coefficient.
Testing \(H_0:\theta=0\) gives a regression-based version of the Hausman test and works with cluster-robust standard errors. We can also add time-invariant covariates directly. They do not change \(\hat\beta\) because within-unit deviations are orthogonal to every time-invariant variable.
8.2 A one-way panel example
I start with Stata’s nlswork data, using 19,229 person-year observations:
Code
clearwebuse nlsworkxtset idcode yeardropifunion==. | age==.bysort idcode: egen mean_age = mean(age)bysort idcode: egen mean_union = mean(union)* 1. Standard FExtreg ln_wage age union, fecluster(idcode)* 2. Pooled OLS with Mundlak meansreg ln_wage age union mean_age mean_union, cluster(idcode)* 3. Random Effects with Mundlak meansxtreg ln_wage age union mean_age mean_union, re cluster(idcode)
(National Longitudinal Survey of Young Women, 14-24 years old in 1968)
Panel variable: idcode (unbalanced)
Time variable: year, 68 to 88, but with gaps
Delta: 1 unit
(9,305 observations deleted)
Fixed-effects (within) regression Number of obs = 19,229
Group variable: idcode Number of groups = 4,150
R-squared: Obs per group:
Within = 0.0963 min = 1
Between = 0.0433 avg = 4.6
Overall = 0.0562 max = 12
F(2, 4149) = 331.52
corr(u_i, Xb) = 0.0127 Prob > F = 0.0000
(Std. err. adjusted for 4,150 clusters in idcode)
------------------------------------------------------------------------------
| Robust
ln_wage | Coefficient std. err. t P>|t| [95% conf. interval]
-------------+----------------------------------------------------------------
age | .0153507 .0006912 22.21 0.000 .0139956 .0167058
union | .1055274 .0098582 10.70 0.000 .0862001 .1248546
_cons | 1.248435 .0215661 57.89 0.000 1.206154 1.290716
-------------+----------------------------------------------------------------
sigma_u | .42353003
sigma_e | .26213464
rho | .72302816 (fraction of variance due to u_i)
------------------------------------------------------------------------------
Linear regression Number of obs = 19,229
F(4, 4149) = 234.77
Prob > F = 0.0000
R-squared = 0.0769
Root MSE = .44953
(Std. err. adjusted for 4,150 clusters in idcode)
------------------------------------------------------------------------------
| Robust
ln_wage | Coefficient std. err. t P>|t| [95% conf. interval]
-------------+----------------------------------------------------------------
age | .0153507 .0006912 22.21 0.000 .0139956 .0167059
union | .1055274 .0098587 10.70 0.000 .0861991 .1248556
mean_age | -.0064208 .001548 -4.15 0.000 -.0094556 -.0033859
mean_union | .1900785 .0223172 8.52 0.000 .1463249 .2338321
_cons | 1.405264 .042734 32.88 0.000 1.321482 1.489046
------------------------------------------------------------------------------
Random-effects GLS regression Number of obs = 19,229
Group variable: idcode Number of groups = 4,150
R-squared: Obs per group:
Within = 0.0963 min = 1
Between = 0.0750 avg = 4.6
Overall = 0.0768 max = 12
Wald chi2(4) = 1023.99
corr(u_i, X) = 0 (assumed) Prob > chi2 = 0.0000
(Std. err. adjusted for 4,150 clusters in idcode)
------------------------------------------------------------------------------
| Robust
ln_wage | Coefficient std. err. z P>|z| [95% conf. interval]
-------------+----------------------------------------------------------------
age | .0153507 .0006912 22.21 0.000 .0139959 .0167055
union | .1055274 .0098587 10.70 0.000 .0862047 .12485
mean_age | -.0061464 .0014254 -4.31 0.000 -.0089401 -.0033527
mean_union | .212183 .0210098 10.10 0.000 .1710046 .2533614
_cons | 1.362101 .0383784 35.49 0.000 1.286881 1.437322
-------------+----------------------------------------------------------------
sigma_u | .38505558
sigma_e | .26213464
rho | .68331727 (fraction of variance due to u_i)
------------------------------------------------------------------------------
All three models yield identical coefficients: \(\hat{\beta}_{\text{age}} = 0.0153507\) (\(\text{SE} = 0.00069\)) and \(\hat{\beta}_{\text{union}} = 0.1055274\) (\(\text{SE} = 0.00986\)).
8.2.1 Adding a time-invariant covariate
We add race (constant within individual) to the Mundlak Pooled OLS model:
Code
reg ln_wage age union mean_age mean_union i.race, cluster(idcode)
Linear regression Number of obs = 19,229
F(6, 4149) = 189.66
Prob > F = 0.0000
R-squared = 0.1069
Root MSE = .4422
(Std. err. adjusted for 4,150 clusters in idcode)
------------------------------------------------------------------------------
| Robust
ln_wage | Coefficient std. err. t P>|t| [95% conf. interval]
-------------+----------------------------------------------------------------
age | .0153507 .0006913 22.21 0.000 .0139955 .0167059
union | .1055274 .0098592 10.70 0.000 .0861981 .1248566
mean_age | -.0070345 .0015342 -4.59 0.000 -.0100424 -.0040266
mean_union | .2227511 .0220026 10.12 0.000 .1796142 .265888
|
race |
Black | -.1790995 .0143624 -12.47 0.000 -.2072575 -.1509415
Other | .0882461 .0715315 1.23 0.217 -.0519939 .2284861
|
_cons | 1.466452 .0426087 34.42 0.000 1.382916 1.549988
------------------------------------------------------------------------------
The coefficients on age and union remain unchanged to seven digits, while the coefficient on race is identified (Black women earn \(-17.9\%\) relative to white women, \(\text{SE} = 0.0144\)).
8.3 Two-way Mundlak and balanced panels
Wooldridge (2021) extends Mundlak’s idea to Two-Way Fixed Effects (TWFE):
\[y_{it} = X_{it}\beta + c_i + f_t + u_{it}\]
The Two-Way Mundlak (TWM) estimator adds both unit means (\(\bar{X}_{i\cdot}\)) and time means (\(\bar{X}_{\cdot t}\)):
The simple equivalence between TWM and TWFE requires the unit and time factors to be orthogonal. This holds automatically in a balanced panel, where every unit is observed in every period.
In an unbalanced panel, separate unit and time means do not span the two-way dummy-variable space. To recover the equivalence, we also include the unit means of the time dummies (\(\bar d_{it}\)). These means describe each unit’s observation pattern, \(s_{it}/T_i\).
Code
preservewebuse nlswork, clearkeepif !missing(ln_wage, age, union)xtset idcode yearbysort idcode: egen m_age = mean(age)bysort idcode: egen m_union = mean(union)egen t_age = mean(age), by(year)egen t_union = mean(union), by(year)* Unit meansof time dummieslevelsofyear, local(yrs)foreachyoflocal yrs {bysort idcode: egen my_`y' = mean(year==`y')}quietlyxtreg ln_wage age union i.year, fescalar b_twfe = _b[age]* Naive TWM with unit and time means onlyquietlyregress ln_wage age union m_age m_union t_age t_unionscalar b_twm = _b[age]* Correct TWM with unit meansof time dummiesquietlyregress ln_wage age union m_age m_union i.year my_*scalar b_fix = _b[age]di"TWFE age = " %9.7f b_twfedi"Naive TWM (unit + time means) age = " %9.7f b_twmdi"TWM + unit means of dummies age = " %9.7f b_fixrestore
(National Longitudinal Survey of Young Women, 14-24 years old in 1968)
(9,305 observations deleted)
Panel variable: idcode (unbalanced)
Time variable: year, 70 to 88, but with gaps
Delta: 1 unit
70 71 72 73 77 78 80 82 83 85 87 88
TWFE age = 0.0276645
Naive TWM (unit + time means) age = 0.0117968
TWM + unit means of dummies age = 0.0276646
TWFE gives \(\hat\beta_{\text{age}}=0.0276645\). The simple TWM regression gives only \(0.0117968\). Once we add the unit means of the time dummies, the estimate is \(0.0276646\), equal to TWFE to six decimal places.
8.4 More than two fixed effects
With three or more fixed-effect factors, such as individual \(i\), year \(t\), and group \(g\), one set of means for each dimension is not generally enough. It works only when the factors are mutually orthogonal, as in a Latin-square design. A safer approach is to eliminate the high-dimensional factor, usually the individual, keep dummies for the other factors, and add the within-individual means of every time-varying regressor, including those dummies.
Code
preservecaptureprogramdrop twm3programdefine twm3argslabelquietly {capturedrop x_i x_t x_g mt_* mg_* x y ft hgsort tby t: gendouble ft = rnormal() if_n==1by t: replace ft = ft[1]sort gby g: gendouble hg = rnormal() if_n==1by g: replace hg = hg[1]gendouble x = 0.6*ci + 0.5*ft + 0.4*hg + rnormal()gendoubley = 1*x + ci + ft + hg + rnormal()bysort i: egendouble x_i = mean(x)bysort t: egendouble x_t = mean(x)bysort g: egendouble x_g = mean(x)levelsof t, local(ts)foreachkoflocal ts {bysort i: egendouble mt_`k' = mean(t==`k') }levelsof g, local(gs)foreachkoflocalgs {bysort i: egendouble mg_`k' = mean(g==`k') }aregy x i.t i.g, absorb(i)scalar b_fe = _b[x]regressy x x_i x_t x_gscalar b_A = _b[x]regressy x x_i i.t i.g mt_* mg_*scalar b_B = _b[x] }di"`label'"di" Three-way FE = " %9.6f b_fedi" A: One set of means per factor = " %9.6f b_A " gap: " %9.6f (b_A-b_fe)di" B: Means over i, dummies included = " %9.6f b_B " gap: " %9.6f (b_B-b_fe)endclearsetseed 20260731setobs 300gen i = _ngendoubleci = rnormal()expand 10bysort i: gen t = _n* Non-orthogonal design (random assignment of g)gen g = 1 + int(runiform()*5)twm3 "Non-orthogonal factor assignment:"* Orthogonal design (Latin square)replace g = mod(i+t,5)+1twm3 "Orthogonal factor assignment (Latin square):"captureprogramdrop twm3restore
Number of observations (_N) was 0, now 300.
(2,700 observations created)
Non-orthogonal factor assignment:
Three-way FE = 0.994424
A: One set of means per factor = 1.007648 gap: 0.013224
B: Means over i, dummies included = 0.994424 gap: 0.000000
(2,435 real changes made)
Orthogonal factor assignment (Latin square):
Three-way FE = 0.996471
A: One set of means per factor = 0.996471 gap: 0.000000
B: Means over i, dummies included = 0.996471 gap: 0.000000
Construction B (means over \(i\) including factor dummies) is exact in both orthogonal and non-orthogonal designs.
8.5 Difference-in-differences with common timing
In standard DiD with common treatment timing: \[y_{it} = \beta w_{it} + c_i + f_t + u_{it}\] where \(w_{it} = d_i \cdot \text{post}_t\). Because \(d_i\) is time-invariant and \(\text{post}_t\) is unit-invariant: - Unit mean: \(\bar{w}_{i\cdot} = d_i \cdot \bar{p}\) - Time mean: \(\bar{w}_{\cdot t} = \bar{d} \cdot \text{post}_t\)
So we can span the Mundlak projection with \(d_i\) and the post-treatment indicator, or with the full set of time dummies. We can also interact treatment with time to allow the treatment effect to vary:
Code
clearsetobs 100setseed 1234567gen id = _ngen w1 = rnormal()gen w2 = rnormal()expand 4bysort id: genyear=2011+_ngenpost=(year>2013)gend=(id>50)genw=d*postgen x1 = rnormal() + w1 + 2*w2gen u = rnormal()gen f2014=(year==2014)gen f2015=(year==2015)* DGP with growing treatment effect: 2.5 in 2014, 3.5 in 2015 (ATT = 3.0)gen y1 = 1 + 2*x1 + 3*w1 + ugen logy = y1 + 2.5*d*f2014 + 3.5*d*f2015xtset id year
Number of observations (_N) was 0, now 100.
(300 observations created)
Panel variable: id (strongly balanced)
Time variable: year, 2012 to 2015
Delta: 1 unit
Code
* 1. TWFE with year-specific treatment effectsxtreg logy c.w#c.f2014 c.w#c.f2015 i.year, fevce(cluster id)* 2. Pooled OLS with Mundlak controls (d and year dummies)reg logy c.w#c.f2014 c.w#c.f2015 d f2014 f2015, vce(cluster id)lincom (c.w#c.f2014 + c.w#c.f2015)/2* 3. Adding time-invariant controls (x1 and interaction with d)reg logy c.w#c.f2014 c.w#c.f2015 d f2014 f2015 x1 c.d#c.x1, vce(cluster id)
The pooled Mundlak regression reproduces the TWFE estimates of 2.5 and 3.5 in the two post-treatment years. lincom gives their average, 3.0. We can add time-invariant controls and group interactions without changing these treatment-effect estimates.
8.6 What to remember
For a one-way model, adding unit means to pooled OLS or random effects recovers the within coefficient and lets us include time-invariant variables. For two-way models, the simple result holds in a balanced panel; an unbalanced panel also needs unit means of the time dummies. With more fixed-effect dimensions, I would eliminate the high-dimensional factor and retain explicit dummies for the others. The same formulation is useful in DiD because it accommodates time-varying effects and time-invariant moderators in an ordinary regression.
Source Code
---title: "The Mundlak device in panel data and difference-in-differences"date: "2021-10-23"---```{r}#| label: setup#| include: falselibrary(Statamarkdown)stataexe <-find_stata()knitr::opts_chunk$set(engine.path =list(stata = stataexe))```## The basic ideaSuppose we have a panel-data model with an individual effect $c_i$:$$y_{it} = X_{it}\beta + c_i + u_{it}$$The fixed-effects estimator removes $c_i$ by demeaning, so $c_i$ can be correlated with $X_{it}$. But the same transformation also removes time-invariant variables such as gender, race, or baseline schooling.Mundlak (1978) instead writes $c_i$ as a linear projection on the unit means of the time-varying covariates:$$c_i = \bar{X}_i \theta + \eta_i, \qquad E[\eta_i \mid X_i] = 0$$Substituting into the regression yields:$$y_{it} = X_{it}\beta + \bar{X}_i \theta + \eta_i + u_{it}$$We can estimate this equation by pooled OLS or random effects. By the Frisch--Waugh--Lovell theorem, the coefficient on $X_{it}$ is the same as the fixed-effects coefficient.Testing $H_0:\theta=0$ gives a regression-based version of the Hausman test and works with cluster-robust standard errors. We can also add time-invariant covariates directly. They do not change $\hat\beta$ because within-unit deviations are orthogonal to every time-invariant variable.## A one-way panel exampleI start with Stata's `nlswork` data, using 19,229 person-year observations:```{stata}*| label: stata1*| echo: true*| collectcode: trueclearwebuse nlsworkxtset idcode yeardrop if union==. | age==.bysort idcode: egen mean_age = mean(age)bysort idcode: egen mean_union = mean(union)* 1. Standard FExtreg ln_wage age union, fe cluster(idcode)* 2. Pooled OLS with Mundlak meansreg ln_wage age union mean_age mean_union, cluster(idcode)* 3. Random Effects with Mundlak meansxtreg ln_wage age union mean_age mean_union, re cluster(idcode)```All three models yield identical coefficients: $\hat{\beta}_{\text{age}} = 0.0153507$ ($\text{SE} = 0.00069$) and $\hat{\beta}_{\text{union}} = 0.1055274$ ($\text{SE} = 0.00986$).### Adding a time-invariant covariateWe add `race` (constant within individual) to the Mundlak Pooled OLS model:```{stata}*| label: stata2*| echo: true*| collectcode: truereg ln_wage age union mean_age mean_union i.race, cluster(idcode)```The coefficients on `age` and `union` remain unchanged to seven digits, while the coefficient on `race` is identified (Black women earn $-17.9\%$ relative to white women, $\text{SE} = 0.0144$).## Two-way Mundlak and balanced panelsWooldridge (2021) extends Mundlak's idea to Two-Way Fixed Effects (TWFE):$$y_{it} = X_{it}\beta + c_i + f_t + u_{it}$$The Two-Way Mundlak (TWM) estimator adds both unit means ($\bar{X}_{i\cdot}$) and time means ($\bar{X}_{\cdot t}$):$$y_{it} = X_{it}\beta + \bar{X}_{i\cdot}\lambda + \bar{X}_{\cdot t}\theta + u_{it}$$### What changes in an unbalanced panelThe simple equivalence between TWM and TWFE requires the unit and time factors to be orthogonal. This holds automatically in a balanced panel, where every unit is observed in every period.In an unbalanced panel, separate unit and time means do not span the two-way dummy-variable space. To recover the equivalence, we also include the unit means of the time dummies ($\bar d_{it}$). These means describe each unit's observation pattern, $s_{it}/T_i$.```{stata}*| label: twm_unbalanced*| echo: true*| collectcode: truepreservewebuse nlswork, clearkeep if !missing(ln_wage, age, union)xtset idcode yearbysort idcode: egen m_age = mean(age)bysort idcode: egen m_union = mean(union)egen t_age = mean(age), by(year)egen t_union = mean(union), by(year)* Unit means of time dummieslevelsof year, local(yrs)foreach y of local yrs { bysort idcode: egen my_`y' = mean(year==`y')}quietly xtreg ln_wage age union i.year, fescalar b_twfe = _b[age]* Naive TWM with unit and time means onlyquietly regress ln_wage age union m_age m_union t_age t_unionscalar b_twm = _b[age]* Correct TWM with unit means of time dummiesquietly regress ln_wage age union m_age m_union i.year my_*scalar b_fix = _b[age]di "TWFE age = " %9.7f b_twfedi "Naive TWM (unit + time means) age = " %9.7f b_twmdi "TWM + unit means of dummies age = " %9.7f b_fixrestore```TWFE gives $\hat\beta_{\text{age}}=0.0276645$. The simple TWM regression gives only $0.0117968$. Once we add the unit means of the time dummies, the estimate is $0.0276646$, equal to TWFE to six decimal places.## More than two fixed effectsWith three or more fixed-effect factors, such as individual $i$, year $t$, and group $g$, one set of means for each dimension is not generally enough. It works only when the factors are mutually orthogonal, as in a Latin-square design. A safer approach is to eliminate the high-dimensional factor, usually the individual, keep dummies for the other factors, and add the within-individual means of every time-varying regressor, including those dummies.```{stata}*| label: twm_three_way*| echo: true*| collectcode: truepreservecapture program drop twm3program define twm3 args label quietly { capture drop x_i x_t x_g mt_* mg_* x y ft hg sort t by t: gen double ft = rnormal() if _n==1 by t: replace ft = ft[1] sort g by g: gen double hg = rnormal() if _n==1 by g: replace hg = hg[1] gen double x = 0.6*ci + 0.5*ft + 0.4*hg + rnormal() gen double y = 1*x + ci + ft + hg + rnormal() bysort i: egen double x_i = mean(x) bysort t: egen double x_t = mean(x) bysort g: egen double x_g = mean(x) levelsof t, local(ts) foreach k of local ts { bysort i: egen double mt_`k' = mean(t==`k') } levelsof g, local(gs) foreach k of local gs { bysort i: egen double mg_`k' = mean(g==`k') } areg y x i.t i.g, absorb(i) scalar b_fe = _b[x] regress y x x_i x_t x_g scalar b_A = _b[x] regress y x x_i i.t i.g mt_* mg_* scalar b_B = _b[x] } di "`label'" di " Three-way FE = " %9.6f b_fe di " A: One set of means per factor = " %9.6f b_A " gap: " %9.6f (b_A-b_fe) di " B: Means over i, dummies included = " %9.6f b_B " gap: " %9.6f (b_B-b_fe)endclearset seed 20260731set obs 300gen i = _ngen double ci = rnormal()expand 10bysort i: gen t = _n* Non-orthogonal design (random assignment of g)gen g = 1 + int(runiform()*5)twm3 "Non-orthogonal factor assignment:"* Orthogonal design (Latin square)replace g = mod(i+t,5)+1twm3 "Orthogonal factor assignment (Latin square):"capture program drop twm3restore```Construction B (means over $i$ including factor dummies) is exact in both orthogonal and non-orthogonal designs.## Difference-in-differences with common timingIn standard DiD with common treatment timing:$$y_{it} = \beta w_{it} + c_i + f_t + u_{it}$$where $w_{it} = d_i \cdot \text{post}_t$. Because $d_i$ is time-invariant and $\text{post}_t$ is unit-invariant:- Unit mean: $\bar{w}_{i\cdot} = d_i \cdot \bar{p}$- Time mean: $\bar{w}_{\cdot t} = \bar{d} \cdot \text{post}_t$So we can span the Mundlak projection with $d_i$ and the post-treatment indicator, or with the full set of time dummies. We can also interact treatment with time to allow the treatment effect to vary:```{stata}*| label: stata3*| echo: true*| collectcode: trueclearset obs 100set seed 1234567gen id = _ngen w1 = rnormal()gen w2 = rnormal()expand 4bysort id: gen year=2011+_ngen post=(year>2013)gen d=(id>50)gen w=d*postgen x1 = rnormal() + w1 + 2*w2gen u = rnormal()gen f2014=(year==2014)gen f2015=(year==2015)* DGP with growing treatment effect: 2.5 in 2014, 3.5 in 2015 (ATT = 3.0)gen y1 = 1 + 2*x1 + 3*w1 + ugen logy = y1 + 2.5*d*f2014 + 3.5*d*f2015xtset id year``````{stata}*| label: stata4*| echo: true*| collectcode: true* 1. TWFE with year-specific treatment effectsxtreg logy c.w#c.f2014 c.w#c.f2015 i.year, fe vce(cluster id)* 2. Pooled OLS with Mundlak controls (d and year dummies)reg logy c.w#c.f2014 c.w#c.f2015 d f2014 f2015, vce(cluster id)lincom (c.w#c.f2014 + c.w#c.f2015)/2* 3. Adding time-invariant controls (x1 and interaction with d)reg logy c.w#c.f2014 c.w#c.f2015 d f2014 f2015 x1 c.d#c.x1, vce(cluster id)```The pooled Mundlak regression reproduces the TWFE estimates of 2.5 and 3.5 in the two post-treatment years. `lincom` gives their average, 3.0. We can add time-invariant controls and group interactions without changing these treatment-effect estimates.## What to rememberFor a one-way model, adding unit means to pooled OLS or random effects recovers the within coefficient and lets us include time-invariant variables. For two-way models, the simple result holds in a balanced panel; an unbalanced panel also needs unit means of the time dummies. With more fixed-effect dimensions, I would eliminate the high-dimensional factor and retain explicit dummies for the others. The same formulation is useful in DiD because it accommodates time-varying effects and time-invariant moderators in an ordinary regression.