I'm trying to replicate in R the results of a paper that uses the Stata
reghdfe command but keep running into some issues. I've attached the dataset I'm using here (5.49 MB).
Ultimately, I want to obtain the same results that these Stata commands output:
import delimited sample_data.csv, case(preserve) xtset ID year reghdfe y l(0/1).delta_t l(0/1).delta_p, absorb(i.fact i.year i.ID i.ID#c.year i.ID#c.year_sqr)
The coefficients produced are as follows:
- delta_t = -0.00243
- lag of delta_t = -0.00652
- delta_p = 0.00175
- lag of delta_p = 0.0195.
In R, I'd like to use the
fixest package. I've posted an issue on the package's GitHub, but it got closed before I managed to fully test the answer. My R code is as follows:
data_fe <- data.table::fread('sample_data.csv') data_fe_p <- panel(data_fe, ~ ID + year) feols(y ~ l(delta_t, 0:1) + l(delta_p, 0:1) | fact + ID + year + ID[year, year_sqr], data = data_fe_p)
but this results in the following - quite different - coefficients:
- delta_t = -0.000833
- lag of delta_t = -0.004863
- delta_p = 0.002202
- lag of delta_p = 0.021593.
I can match the Stata output exactly if I remove the
year_sqr interaction in both softwares, i.e., if I do
reghdfe y l(0/1).delta_t l(0/1).delta_p, absorb(i.fact i.year i.ID i.ID#c.year) in Stata and
feols(y ~ l(delta_t, 0:1) + l(delta_p, 0:1) | fact + ID + year + ID[year], data = data_fe_p).
Any idea how to include the interaction and match the outputs?