Constant difference between `DuckMundlakEventStudy` and `feols`

[patrick-dd]

When calculating event study coefficients, the coefficients using DuckMundlakEventStudy differ from those in feols by a constant number.

I expected them to be the same. Perhaps this is my misunderstanding of the method.

Replication

• Use event study notebook
• reduce sample size to 100_000 and 50_000 treated (working on a small machine)
• add this to cell prior to staggered adoption header

mundlak_estimates['evstudy_coefs'] = evstudy_coefs
mundlak_estimates['diff'] = mundlak_estimates['est'] - mundlak_estimates['evstudy_coefs']
mundlak_estimates

notice constant difference. it's small for seed 42, which could be numerical. but big enough for 421, and big enough in my use case that prompted the investigation.

[s3alfisc]

Hi @patrick-dd , thanks for reporting! I can reproduce the issue & can confirm that it's independent of the seeding and sample sizes:

I vaguely recall that @apoorvalal and I discussed a similar phenomenon a while back but cannot recall the conclusion ...

[s3alfisc]

Add.: the observed difference is equal to the coefficient on the reference period set in the feols call. Here is from another example, where the reference period is set to "5":

[s3alfisc]

Probably something with reference levels?

Below the Mundlak specification from Wooldridge (eq 5.13) and TWFE, which should be equivalent. If I drop the same reference level to avoid multicollinearity I get results that are identical up to a numerical difference:

import numpy as np
import pyfixest as pf 
from pyfixest.utils.dgps import get_sharkfin

data = get_sharkfin(seed = 100)

# create year variables: 
for year in data["year"].unique():
    data[f"time{year}"] = (data["year"] == year).astype(int)

ref_year = 14
time_dummies = [f"time{year}" for year in data["year"].unique() if year != ref_year]

fit_mundlak = pf.feols(
    f"""
        Y ~ ever_treated*({"+".join(time_dummies)})
    """, 
    data = data
)

fit_twfe = pf.feols(f"Y ~ i(year, ever_treated, ref = {ref_year}) | unit + year", data = data)

fit_mundlak.coef()[31:].values - fit_twfe.coef().values
# array([2.52020627e-14, 2.46469511e-14, 2.52020627e-14, 2.49800181e-14,

If I don't drop the time14 variable by hand instead, pyfixest will do it for me, and I get something close to what you observe above @patrick-dd :

time_dummies = [f"time{year}" for year in data["year"].unique()]

fit_mundlak = pf.feols(
    f"""
        Y ~ ever_treated*({"+".join(time_dummies)})
    """, 
    data = data
)

fit_twfe = pf.feols(f"Y ~ i(year, ever_treated, ref = {ref_year}) | unit + year", data = data)
fit_mundlak.coef()[31:].values - fit_twfe.coef().values

# UserWarning: 
#            The following variables are collinear: ['time29', 'ever_treated:time29'].
#            The variables are dropped from the model.
#            
#  warnings.warn(
# array([0.52581268, 0.52581268, 0.52581268, 0.52581268, 0.52581268,
#        0.52581268, 0.52581268, 0.52581268, 0.52581268, 0.52581268,
#        0.52581268, 0.52581268, 0.52581268, 0.52581268, 0.52581268,
#        0.52581268, 0.52581268, 0.52581268, 0.52581268, 0.52581268,
#        0.44974662, 0.49758658, 1.08311083, 0.60277033, 0.53944695,
#        0.52658065, 0.54057793, 0.49257393, 0.52573303])

[s3alfisc]

Yes, after some "offline" discussion agreeing with @apoorvalal that this should be about reference levels. One of the time dummies will have to be dropped, feols() drops ref=14, while duckreg might drop any of the time dummies.

import numpy as np
import pandas as pd
import pyfixest as pf 
from pyfixest.utils.dgps import get_sharkfin

data = get_sharkfin()

# create year variables: 
for year in data["year"].unique():
    data[f"time{year}"] = (data["year"] == year).astype(int)

ref_year = 14
time_dummies = [f"time{year}" for year in data["year"].unique() if year != ref_year]

fit_mundlak = pf.feols(
    f"""
        Y ~ ever_treated*({"+".join(time_dummies)})
    """, 
    data = data
)

fit_twfe = pf.feols(f"Y ~ i(year, ever_treated, ref = {ref_year}) | year + unit", data = data)

mundlak_ref = 0
mundlak_beta = fit_mundlak.coef()[31:].values
twfe_beta = fit_twfe.coef().values

# raw estimates
np.c_[mundlak_beta, twfe_beta, mundlak_beta - twfe_beta]

# with offset 
np.c_[mundlak_beta , twfe_beta, (mundlak_beta) - twfe_beta]

Maybe it would make sense to add a function argument to set a reference level in the DuckMundlakEventStudy class @apoorvalal?