Here is the canonical DD equation with two groups and two periods:
$$ y_{ist} = \alpha + \gamma T_{s} + \lambda d_{t} + \beta(T_{s} \cdot d_{t}) + \epsilon_{ist}, $$
where, for example, we may observe individual/entity $i$, in state $s$, at time period $t$. In this setting, $T_{s}$ indexes only those states exposed to treatment, 0 otherwise. The variable $d_{t}$ indexes periods after treatment in both treatment and control groups. Because $d_{t}$ is the same across all $s$, this model is used when treated states enter into the treatment condition at precisely the same time.
The generalization of this equation would include dummies for each state and each time period but is otherwise unchanged. For example,
$$ y_{ist} = \gamma_{s} + \lambda_{t} + \beta D_{st} + \epsilon_{ist}, $$
where $D_{st}$ is equal to unity for treated states during periods when treatment is in effect. $\gamma_{s}$ denotes state (unit) fixed effects; $\lambda$ denotes year (time) fixed effects. Note, these fixed effects replace $T_{s}$ and $d_{t}$, respectively, in the former equation. $D_{st}$ is the same as before $(T_{s} \cdot d_{t})$. Instead of doing this interaction manually, we code this dummy explicitly to reflect early/late adopter states, or possibly ones experiencing intermittent treatment exposure. It is for these reasons that researchers estimate the equation you referenced. Please review this post which details the coding of the treatment dummy. You can find further insights here.
Why is it that so many papers use separate group and time fixed effects?
This is a requirement. DD performs a double-difference across units and across time. At a basic level, it is an interaction model. Put more simply, it assesses the before-and-after change in units exposed to treatment versus the before-and-after change in units unexposed to treatment. The more general case you referenced in your question is a 'two-way' fixed effects estimator, and it accommodates treatment exposure in multiple groups and multiple times periods. Once again, the variable $D_{st}$ is your interaction term. In practice, treatment exposure is often staggered and doesn’t always follow a precise pattern for some treated entities. Because of this, we regress the outcome on unit-specific effects, time-specific effects, and a treatment dummy. The main causal parameter of interest is akin to a weighted combination of all possible two-group/two-period DD estimators that can be constructed from your panel.
Why not use group times time fixed effects?
Multiplying fixed effects will often fail in practice. In most applications, there is not enough degrees of freedom to multiply the unit and time fixed effects. This equation attempts estimation of main effects for units (i.e., dummies for states), main effects for time (e.g., dummies for all years), and each pairwise interaction between unit and time. Thus, you would ‘chew up’ all your degrees of freedom. In other words, the model would be perfectly fit and you wouldn't be able to estimate your standard errors.
Suppose you observe 10 states over 10 years. Your total sample size ($N \times T$ = 100). Interacting the state effects with a discretized version of year results in the estimation of 99 dummies (i.e., 9 state dummies, 9 year dummies, and 81 state-year dummies). In addition to a constant, a treatment dummy, and possibly some time-varying covariates, you would have more parameters to estimate than observations.
In some DD applications, however, researchers interact state-specific effects (i.e., state dummies) with a linear time index. To be clear, this is not a discretized version of time. It is a continuous linear time trend variable (e.g., $t = 1, 2, 3, 4,…,T$). This is not equivalent to interacting the state-specific effects with individual year dummies.
For instance, in the case of firms, why not include industry-by-year and state-by-year fixed effects instead? Since these nest both year and state fixed effects, the coefficient of interest 𝛽3 still has the same interpretation as far as I can tell. However, the higher dimensionality of fixed effects would tighten the identification.
Time effects adjust for those "common shocks" affecting all states. Put differently, you're adjusting for potential effects that are constant across all states within a year. To address your question, estimating dummies for a concatenated version of 'state-year' in a state-year panel would estimate dummies for all state-year observations, which is more parameters than degrees of freedom.
If you're working with micro-data, then you have multiple $i$ (e.g., individuals/entities) nested within states. In this setting, I don't see how a single 'state-year' fixed effect would return the same estimate of $\beta$ (i.e., DD coefficient) without the accompanying state and time effects.
Please let me know if this addresses your concerns. I will adjust this answer if necessary. Or, maybe some else could shed further light on how this could work using data at the individual level.
This post may also be of interest to you.
Based upon our discussion, I simulated a three-level panel dataset with individual firms. Two firms embedded in 2 states observed over 3 years. A ‘state-year’ effect would absorb a treatment dummy when that treatment is implemented at a higher level of aggregation, most particularly when it affects level $s$ and only varies at the ‘state-year’ level. And it will not return the same estimate of $\beta$ if estimated with separate state and time effects.
If I misperceived your question, please let me know and I will adjust my response accordingly.
state year firm state_yr state_fe time_19 time_20 ny_19 ny_20 ca_18 ca_19 ca_20
NY 2018 1 NY-2018 0 0 0 0 0 0 0 0
NY 2019 1 NY-2019 0 1 0 1 0 0 0 0
NY 2020 1 NY-2020 0 0 1 0 1 0 0 0
NY 2018 2 NY-2018 0 0 0 0 0 0 0 0
NY 2019 2 NY-2019 0 1 0 1 0 0 0 0
NY 2020 2 NY-2020 0 0 1 0 1 0 0 0
CA 2018 1 CA-2018 1 0 0 0 0 1 0 0
CA 2019 1 CA-2019 1 1 0 0 0 0 1 0
CA 2020 1 CA-2020 1 0 1 0 0 0 0 1
CA 2018 2 CA-2018 1 0 0 0 0 1 0 0
CA 2019 2 CA-2019 1 1 0 0 0 0 1 0
CA 2020 2 CA-2020 1 0 1 0 0 0 0 1