Fixed effects regression with state specific trends

1
votes

I'm working with panel data and I want to estimate a fixed effects regression with state specific trends.

In Stata, I could accomplish this by the following,

xi i.state i.year i.state*time
reg y x _I*

The above will create state dummies, year dummies, and 50 (state x time) dummies where time numerically identifies the trend (i.e. 1, 2, 3...)

In R, I can run a fixed effects model with plm or lm, for example,

plm(y ~ x, index = c("state", "year"), effect = "twoways", data = df)
 lm(y ~ x + factor(state) + factor(year), data = df)

How would I include the 50 (state x time) dummies the way xi does in Stata?

I know interaction() is not what I want because that creates a new factor variable with n levels, where n = (num of states) x (num of time periods). What I'm trying to do is create 50 (state x time) variables such that state1xtime is 1,2,3... when state == 1 and zero otherwise, repeat for state2xtime, where state == 2, etc.

3
rstata
I'm not sure why you use factor(year). I think just lm(y ~ x + factor(state) + year + factor(state):year, data = df) -- where I'm assuming you're centering your year variable. - Alex W

3 Answers

2
votes

You would do this simply interacting state with year. The correct operator for this is :, which only includes the interactions terms.

Note that there is a subtle difference between lm and plm:

  • for lm, just do state:year
  • for plm, year has been converted implicitly to a factor, so do state:as.integer(year) (doing state:year would give you all combinations of state and year).

check:

library(plm)

data("Produc", package = "plm")
produc_plm <- pdata.frame(Produc, index = c("state","year"))

### simple state-specific time trend ###
fe1_Ttrend_lm <- lm(log(gsp) ~ -1 + log(pcap) + log(pc) + log(emp) + unemp + state +state:year, data = Produc)
fe1_Ttrend_plm <- plm(log(gsp) ~ log(pcap) + log(pc) + log(emp) + unemp + state : as.integer(year), data = produc_plm)

summary(fe1_Ttrend_lm)$coef[1:4,]
summary(fe1_Ttrend_plm)$coef[1:4,]
2
votes

To complete Matifou's answer, you could also do it with the package fixest:

library(fixest)
data("Produc", package = "plm")
fe_fixest = feols(log(gsp) ~ log(pcap) + log(pc) + log(emp) + unemp + year::state | state, data = Produc)
# Notice the double colon (on the left the numeric variable, on the right the factor). The alternative also works,
# fe_fixest = feols(log(gsp) ~ log(pcap) + log(pc) + log(emp) + unemp + interact(year, state) | state, data = Produc)
# fe_fixest = feols(log(gsp) ~ log(pcap) + log(pc) + log(emp) + unemp + i(year, state) | state, data = Produc)

# Requesting ("standard" standard-errors, otherwise, clustered at state level by default)
coeftable(fe_fixest, se = "standard")[1:4, ]

Note that if you don't care about the standard-errors of the interacted coefficients, you can project them out using the following syntax:

fe_fixest_bis = feols(log(gsp) ~ log(pcap) + log(pc) + log(emp) + unemp | state[year], data = Produc)
# state[year] means 'state' fixed-effects and 'state x year' interactions
# The interacted terms are projected out, and the estimation is faster

coeftable(fe_fixest_bis, se = "s")
1
votes

Could this be what you are looking for:

dta <- data.frame(state = rep(LETTERS[1:3], 10), 
              time = rep(1:3, each = 10))
dta <- cbind(dta, model.matrix( ~ state - 1, data = dta) * dta$time)

head(dta, 1)
#     state time stateA stateB stateC
# 1     A    1      1      0      0

tail(dta, 1)
#      state time stateA stateB stateC
# 30     C    3      0      0      3