The Two-Way Mundlak Estimator

1. Executive Summary & Overview

The Empirical Challenge in Panel Data

When estimating panel data regressions with two-way unobserved heterogeneity:

y_{it} = X_{it}\beta + Z_i\gamma + W_t\delta + \mu_i + \lambda_t + \epsilon_{it}

• Two-Way Fixed Effects (TWFE):
  • Advantage: Controls for arbitrary correlation between regressors X_{it} and individual effects \mu_i or time effects \lambda_t.
  • Limitation: Drops time-invariant variables Z_i (e.g., distance, gender, education) and individual-invariant variables W_t (e.g., inflation, macro policy, oil prices).
• Two-Way Random Effects (RE-GLS):
  • Advantage: Retains Z_i and W_t and is efficient under orthogonality.
  • Limitation: Biased and inconsistent if \text{Cov}(X_{it}, \mu_i) \neq 0 or \text{Cov}(X_{it}, \lambda_t) \neq 0.

2. Econometric Methodology: Baltagi (2023)

The Two-Way Mundlak Transformation

Following Mundlak (1978), Baltagi (2023) and Wooldridge (2021) model the unobserved individual and time effects as linear functions of the group averages of all time-varying regressors:

\mu_i = \bar{X}_{i.}\theta_1 + u_i \quad \text{where } \bar{X}_{i.} = \frac{1}{T}\sum_{t=1}^T X_{it}

\lambda_t = \bar{X}_{.t}\theta_2 + v_t \quad \text{where } \bar{X}_{.t} = \frac{1}{N}\sum_{i=1}^N X_{it}

Substituting \mu_i and \lambda_t into the outcome equation yields the Two-Way Mundlak (TWM) Augmented Regression:

y_{it} = X_{it}\beta + Z_i\gamma + W_t\delta + \bar{X}_{i.}\theta_1 + \bar{X}_{.t}\theta_2 + (u_i + v_t + \epsilon_{it})

4. Step 1: Data Simulation Setup

suppressPackageStartupMessages({
  library(dplyr)
  library(fixest)
  library(plm)
  library(car)
})

set.seed(42)
N <- 100
T_per <- 10
n_obs <- N * T_per

# Create panel structure
id <- rep(1:N, each = T_per)
time <- rep(1:T_per, times = N)

# Draw unobserved individual (mu_i) and time (lambda_t) effects
mu_i <- rep(rnorm(N, mean = 2, sd = 1), each = T_per)
lambda_t <- rep(rnorm(T_per, mean = 0.5, sd = 0.5), times = N)

# Generate time-invariant Z_i and individual-invariant W_t
Z_i <- rep(rnorm(N, mean = 5, sd = 1.5), each = T_per)
W_t <- rep(rnorm(T_per, mean = 10, sd = 2), times = N)

# Generate endogenous regressor X_it correlated with mu_i and lambda_t
X_it <- 0.7 * mu_i + 0.4 * lambda_t + rnorm(n_obs, mean = 0, sd = 1)

# True parameters: beta = 1.5, gamma = 0.8, delta = -0.4
eps_it <- rnorm(n_obs, mean = 0, sd = 1)
y_it <- 1.5 * X_it + 0.8 * Z_i - 0.4 * W_t + mu_i + lambda_t + eps_it

# Construct Dataset and calculate individual and time means
df <- data.frame(id, time, y_it, X_it, Z_i, W_t) %>%
  group_by(id) %>%
  mutate(X_bar_i = mean(X_it)) %>%
  ungroup() %>%
  group_by(time) %>%
  mutate(X_bar_t = mean(X_it)) %>%
  ungroup()

head(df, 4)

# A tibble: 4 × 8
     id  time  y_it  X_it   Z_i   W_t X_bar_i X_bar_t
  <int> <int> <dbl> <dbl> <dbl> <dbl>   <dbl>   <dbl>
1     1     1 10.6   3.19  4.96  7.60    2.56    1.94
2     1     2  6.92  2.42  4.96 14.1     2.56    1.86
3     1     3  5.73  1.84  4.96 10.2     2.56    1.42
4     1     4  8.79  1.86  4.96  9.83    2.56    1.90

* Clear workspace and set random seed
clear all
set seed 42

* Define observation count (100 individuals x 10 time periods)
set obs 1000
gen id = ceil(_n / 10)
bysort id: gen time = _n

* Draw individual (mu_i) and time (lambda_t) unobserved effects
bysort id: gen temp_mu = rnormal(2, 1) if _n == 1
bysort id: egen mu_i = max(temp_mu)

bysort time: gen temp_lam = rnormal(0.5, 0.5) if _n == 1
bysort time: egen lambda_t = max(temp_lam)

* Generate Z_i (time-invariant) and W_t (individual-invariant)
bysort id: gen temp_z = rnormal(5, 1.5) if _n == 1
bysort id: egen Z_i = max(temp_z)

bysort time: gen temp_w = rnormal(10, 2) if _n == 1
bysort time: egen W_t = max(temp_w)

* Generate endogenous X_it and outcome y_it
gen X_it = 0.7*mu_i + 0.4*lambda_t + rnormal(0, 1)
gen eps_it = rnormal(0, 1)
gen y_it = 1.5*X_it + 0.8*Z_i - 0.4*W_t + mu_i + lambda_t + eps_it

* Set panel structure and calculate group means
xtset id time
bysort id: egen X_bar_i = mean(X_it)
bysort time: egen X_bar_t = mean(X_it)

5. Model Estimation & Equivalence Demonstration

* 1. Standard Two-Way Fixed Effects (TWFE)
reghdfe y_it X_it, absorb(id time)
estimates store twfe

* 2. Two-Way Mundlak via Pooled OLS
reg y_it X_it Z_i W_t X_bar_i X_bar_t, vce(cluster id)
estimates store twm_ols

* 3. Two-Way Mundlak via Random Effects GLS
xtreg y_it X_it Z_i W_t X_bar_i X_bar_t, re vce(cluster id)
estimates store twm_re

* Print comparison table across models
estimates table twfe twm_ols twm_re, b(%9.6f) se(%9.6f) ///
    keep(X_it Z_i W_t X_bar_i X_bar_t)

6. Hausman Specification Test

# Joint F-test on the coefficients of group means theta_1 and theta_2
# H0: theta_1 = 0 AND theta_2 = 0 (Unobserved effects uncorrelated with X)

f_test_res <- linearHypothesis(fit_twm_ols, c("X_bar_i = 0", "X_bar_t = 0"))
print(f_test_res)

Linear hypothesis test:
X_bar_i = 0
X_bar_t = 0

Model 1: restricted model
Model 2: y_it ~ X_it + Z_i + W_t + X_bar_i + X_bar_t

  Res.Df    RSS Df Sum of Sq      F    Pr(>F)    
1    996 1782.8                                  
2    994 1178.0  2    604.85 255.19 < 2.2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

* Run joint test for significance of the individual and time averages
reg y_it X_it Z_i W_t X_bar_i X_bar_t, vce(cluster id)

test X_bar_i X_bar_t

* Note: Rejecting H0 (p < 0.05) indicates that RE orthogonality fails.
* The TWM model successfully corrects for this correlation while preserving Z_i and W_t.

7.1 On Gravity

The gravity model is the canonical application of two-way fixed effects in trade economics. Bilateral trade depends on variables that are essentially constant over time (distance, common language, colonial history, shared borders) together with variables that vary over time (free trade agreements, GDP).

The standard two-way fixed effects estimator consistently estimates the effect of time-varying regressors but absorbs all time-invariant bilateral characteristics into the pair fixed effects.

The two-way Mundlak estimator augments a pooled regression with the pair mean and year mean of every time-varying regressor. As established earlier, this reproduces the two-way fixed-effects coefficients on the time-varying regressors while simultaneously identifying the time-invariant regressors.

7.2 Data

The example uses the CEPII Gravity dataset (Conte, Cotterlaz & Mayer), a standard bilateral panel of trade flows and gravity covariates spanning 1948–2021 for nearly all country pairs.

Starting from the raw file, a working panel is built by:

• using BACI bilateral trade flows (tradeflow_baci) — the trade-flow variable in the file with the best, most consistent coverage,
• dropping self-pairs,
• restricting to 2010–2019 and to the country pairs observed in every year of that window, so the panel is balanced (see the callout above),
• keeping GDP of the origin and destination country as time-varying controls, and distance, common language, contiguity, and colonial history as time-invariant characteristics.

This yields 189,020 observations on 18,902 country pairs over 10 years. The filtering is done once by prepare_gravity_subset.R, which writes the resulting panel to gravity_twm.csv.

library(fixest)
library(dplyr)
library(modelsummary)

trade <- read.csv("gravity_twm.csv")

glimpse(trade)

Rows: 189,020
Columns: 12
$ pair            <chr> "AFG_AGO", "AFG_AGO", "AFG_AGO", "AFG_AGO", "AFG_AGO",…
$ iso3_o          <chr> "AFG", "AFG", "AFG", "AFG", "AFG", "AFG", "AFG", "AFG"…
$ iso3_d          <chr> "AGO", "AGO", "AGO", "AGO", "AGO", "AGO", "AGO", "AGO"…
$ year            <int> 2010, 2011, 2012, 2013, 2014, 2015, 2016, 2017, 2018, …
$ log_trade       <dbl> 6.4547283, 5.6148092, 6.9101554, 4.5845694, 4.3829387,…
$ fta             <int> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
$ log_distance    <dbl> 8.934982, 8.934982, 8.934982, 8.934982, 8.934982, 8.93…
$ common_language <int> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
$ contiguity      <int> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
$ colony          <int> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
$ log_gdp_o       <dbl> 16.58414, 16.69864, 16.83627, 16.82661, 16.81375, 16.7…
$ log_gdp_d       <dbl> 18.22796, 18.46101, 18.56341, 18.63723, 18.65794, 18.4…

* Built once from the raw CEPII Gravity.csv by prepare_gravity_subset.R
import delimited "gravity_twm.csv", clear
describe

7.3 Two-Way Fixed Effects Benchmark

The benchmark specification estimates

\ln(\text{Trade}_{ijt}) = \beta\, FTA_{ijt} + \psi_1 \ln(GDP)_{it} + \psi_2 \ln(GDP)_{jt} + \alpha_{ij} + \gamma_t + \varepsilon_{ijt}

Distance and the other time-invariant bilateral characteristics are absorbed into the pair effects \alpha_{ij} and dropped from the output.

twfe <- feols(log_trade ~ fta + log_gdp_o + log_gdp_d | pair + year,
              cluster = ~pair, data = trade)

etable(twfe)

                              twfe
Dependent Var.:          log_trade
                                  
fta               0.0411. (0.0239)
log_gdp_o       0.2649*** (0.0308)
log_gdp_d       0.4635*** (0.0284)
Fixed-Effects:  ------------------
pair                           Yes
year                           Yes
_______________ __________________
S.E.: Clustered           by: pair
Observations               189,020
R2                         0.92125
Within R2                  0.00451
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

reghdfe log_trade fta log_gdp_o log_gdp_d, absorb(pair year) vce(cluster pair)

7.4 Constructing the Two-Way Mundlak Variables

For every time-varying regressor (FTA, \ln GDP_o, \ln GDP_d), the pair mean \bar x_{i\cdot} and year mean \bar x_{\cdot t} are computed and added to the pooled regression. Omitting the means for any one time-varying regressor breaks exact equivalence for all coefficients, not only that regressor’s own.

trade_m <- trade |>
  group_by(pair) |>
  mutate(
    fta_pair       = mean(fta),
    log_gdp_o_pair = mean(log_gdp_o),
    log_gdp_d_pair = mean(log_gdp_d)
  ) |>
  ungroup() |>
  group_by(year) |>
  mutate(
    fta_year       = mean(fta),
    log_gdp_o_year = mean(log_gdp_o),
    log_gdp_d_year = mean(log_gdp_d)
  ) |>
  ungroup()

egen fta_pair       = mean(fta),       by(pair)
egen log_gdp_o_pair = mean(log_gdp_o), by(pair)
egen log_gdp_d_pair = mean(log_gdp_d), by(pair)

egen fta_year       = mean(fta),       by(year)
egen log_gdp_o_year = mean(log_gdp_o), by(year)
egen log_gdp_d_year = mean(log_gdp_d), by(year)

7.5 Two-Way Mundlak Estimator

The augmented regression is

\ln(\text{Trade}) = \beta\, FTA + \theta_1\overline{FTA}_{ij\cdot} + \theta_2\overline{FTA}_{\cdot t} + \psi_1 \ln GDP_i + \dots + \delta_1 \ln\text{Distance} + \delta_2\,\text{Language} + \delta_3\,\text{Contiguity} + \delta_4\,\text{Colony} + \varepsilon

No pair or year fixed effects are included here — the Mundlak means substitute for both, which is exactly what makes distance, language, contiguity, and colony identifiable.

mundlak <- feols(
  log_trade ~ fta + fta_pair + fta_year +
    log_gdp_o + log_gdp_o_pair + log_gdp_o_year +
    log_gdp_d + log_gdp_d_pair + log_gdp_d_year +
    log_distance + common_language + contiguity + colony,
  cluster = ~pair, data = trade_m
)

etable(mundlak)

                           mundlak
Dependent Var.:          log_trade
                                  
Constant         -29.57*** (1.402)
fta               0.0411. (0.0239)
fta_pair        0.7498*** (0.0428)
fta_year        -1.136*** (0.2312)
log_gdp_o       0.2650*** (0.0309)
log_gdp_o_pair  0.9300*** (0.0316)
log_gdp_o_year       4.053 (3.184)
log_gdp_d       0.4635*** (0.0284)
log_gdp_d_pair  0.4376*** (0.0290)
log_gdp_d_year      -3.571 (3.161)
log_distance    -1.060*** (0.0186)
common_language 0.6393*** (0.0398)
contiguity      0.9199*** (0.0871)
colony          0.8934*** (0.0556)
_______________ __________________
S.E.: Clustered           by: pair
Observations               189,020
R2                         0.66630
Adj. R2                    0.66627
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

reg log_trade fta fta_pair fta_year ///
    log_gdp_o log_gdp_o_pair log_gdp_o_year ///
    log_gdp_d log_gdp_d_pair log_gdp_d_year ///
    log_distance common_language contiguity colony, vce(cluster pair)

7.6 Comparing Estimates

The key theoretical prediction is that

\hat\beta_{TWFE} = \hat\beta_{Mundlak}

modelsummary(
  list("TWFE" = twfe, "Two-Way Mundlak" = mundlak),
  coef_map = c(fta = "FTA", log_gdp_o = "log(GDP origin)",
               log_gdp_d = "log(GDP destination)",
               log_distance = "log(Distance)",
               common_language = "Common language",
               contiguity = "Contiguity", colony = "Colonial ties"),
  stars = TRUE
)

eststo clear

reghdfe log_trade fta log_gdp_o log_gdp_d, absorb(pair year) vce(cluster pair)
eststo twfe

reg log_trade fta fta_pair fta_year ///
    log_gdp_o log_gdp_o_pair log_gdp_o_year ///
    log_gdp_d log_gdp_d_pair log_gdp_d_year ///
    log_distance common_language contiguity colony, vce(cluster pair)
eststo mundlak

esttab twfe mundlak, keep(fta log_gdp_o log_gdp_d)

7.7 Interpretation

On the balanced 2010–2019 CEPII panel (189,020 observations, 18,902 pairs):

• \hat\beta_{\text{FTA}} is 0.0411 under TWFE and 0.0411 under the Two-Way Mundlak estimator — the same to four significant figures, and both only marginally significant (p<0.1) once GDP is controlled for.
• The GDP elasticities match exactly across estimators too (\hat\psi_{\text{origin}} \approx 0.265, \hat\psi_{\text{destination}} \approx 0.464 under both), confirming the equivalence holds for every time-varying regressor, not just the one of headline interest.
• The Mundlak estimator additionally recovers statistically significant coefficients on the time-invariant characteristics that TWFE cannot identify: log-distance (-1.06, close to the canonical gravity elasticity of -1), common language (+0.64), contiguity (+0.92), and colonial ties (+0.89) — all with the sign gravity theory predicts.
• The joint Wald test H_0: \theta_1=\theta_2=0 on the FTA pair and year means is rejected (p \approx 10^{-77}) — the Hausman-type diagnostic from Section 6 — confirming the unobserved pair and year effects are correlated with FTA status, exactly the situation where a random-effects model without the Mundlak correction would be biased.

The two-way Mundlak estimator preserves the identifying assumptions of the fixed-effects model while recovering coefficients fixed effects cannot identify — provided, as shown above, that every time-varying regressor’s own pair and year means are included and the panel is balanced.

References

• Baltagi, B. H. (2023). The Two-way Mundlak Estimator. Center for Policy Research Working Paper No. 256, Syracuse University.
• Mundlak, Y. (1978). On the pooling of time series and cross section data. Econometrica, 46(1), 69–85.
• Wooldridge, J. M. (2021). Two-way fixed effects, the two-way Mundlak regression, and Difference-in-Differences estimators. SSR Working Paper.
• Kang, S. (1985). A note on the equivalence of specification tests in the two-way error component model. Economics Letters, 17(4), 341–343.