Diff-in-Diff time varying treatment I would like to estimate a Difference in Difference specification where the treatment (i.e., the policy change P) is off then on and then off again:
P=0 (0, …, T); P=1 (T+1, …, T+N); P=0 (T+N+1…) 

In this instance, the policy change is completely exogenous and isn't related to covariates or expectations or anticipations.
I'm wondering how to do this, since most DnD settings seem to have just two generic periods (off and on), whereas in my case there are three “regimes” and it’s also of interest to assess whether, even after the policy change is finally closed, there are still some lingering effects on the outcome.
Any papers that’ve looked at this, or methods suitable to such a configuration?
 A: Your setup is suitable for a "classical" difference-in-differences (DD) analysis. For example, you might wish to estimate the following:
$$
y_{ct} = \gamma T_{c} + \lambda P_{t} + \delta (T_{c} \times P_{t}) + \epsilon_{ct},
$$
where you observe cities $c$ across time periods $t$. $T_{c}$ equals 1 for your treated city, 0 for your non-treated city. $P_{t}$ indexes all post-treatment time periods in both your treated and non-treated city. You can extend this equation and investigate effects in all post-exposure periods. To do this, simply interact your treatment dummy with individual post-treatment time dummies. Suppose your policy is implemented in city $c$ in the year 2000. I assume for the purposes of this example your data is 'aggregated up' to the city-year level. However, due to budgetary constraints, the policy is repealed at the conclusion of 2001. Thus, your treatment is only in effect for two full years. Now suppose you continued to observe your outcome for several years after the policy was officially removed, in part because theory suggests effects might linger. You could then do the following:
$$
y_{ct} = \gamma T_{c} + \lambda_{1} (T_{c}*\mathbf{I}_{t = 2000}) + \lambda_{2} (T_{c}*\mathbf{I}_{t = 2001}) + \lambda_{3} (T_{c}*\mathbf{I}_{t = 2002}) + \lambda_{j} (T_{c}*\mathbf{I}_{t = j}) + \epsilon_{ct},
$$
where $\mathbf{I}$ is the indicator function. You may have seen this represented with a numeral (i.e., $\mathbf{1}_{\kappa}$). In words, set the value to 1 if condition $\kappa$ is satisfied, 0 otherwise. Put in context, if a city is in a particular post-treatment year then set it equal to 1, 0 otherwise. You could do this for all post-treatment years—even the periods after treatment is removed. Your interaction of $T_{c}$ with the time dummy for 2002 (i.e., $\mathbf{I}_{t = 2002}$) is the first year where the policy is officially 'off' (i.e., withdrawn). Effects may grow or fade over time in the "after" periods. To save space, I only included three post-treatment years. This could generalize to all $j$ periods (i.e., $\mathbf{I}_{t = j}$) after the initial treatment year (i.e., 2000). You could do these interactions manually, or you could leverage Stata's factor variable notation and let software do these interactions for you.
I should note a potential problem with this approach. You have two cities observed across 90 periods (e.g., weeks, months, years, etc.) for total of 180 city-year observations. You indicated you have 40 post-treatment periods; this includes 20 $t$ periods during which the policy is in effect and 20 $t$ periods after which the policy is terminated. The latter equation will result in 40 main effects for years and 40 city-year ($c\times t$) interactions. In addition to a main effect for treatment and other covariates, your model has over 80 parameters to be estimated. The ratio of observations-to-parameters is very small and so I am hesitant to advise you further using this approach.
I propose two recommendations. First, limit the number of interactions. You do not have to estimate separate interactions for all post-treatment periods. To be clear, post-treatment usually denotes your "during" and "after" treatment periods. What does theory tell you about the lingering effects of treatment? Is treatment likely to persist for many periods after it is removed? Or, is it likely to revert back to pre-policy levels in just a few post-policy periods? You could try limiting dynamic effects to the "after" period. For example, you could interact your treatment dummy with one post-treatment indicator indexing the entire "during" epoch. You would then add separate interactions of your treatment dummy with the individual "after" periods. Thus, you will only be assessing the dynamics of treatment in the 20 periods after treatment is removed.
Second, try to get more data! Right now you only have two cities. What about other untreated jurisdictions? Observing more cities over time will be helpful later on when you have to adjust your standard errors. In fact, it is well known that cluster-robust variance estimators perform poorly with a scanty number of clusters. If finding more units (cities) is difficult in your setting, it might also be worthwhile to look into synthetic control methods.
The "classical" DD approach can also be extended to exploit your different treatment regimes. For instance, you could 'discretize' your time index into before, during, and after treatment periods. Your equation would look something like the following:
$$
y_{ct} = \gamma T_{c} + \lambda_{1} \text{During}_{t} + \lambda_{2} \text{After}_{t} + \delta_{1} (T_{c} \times \text{During}_{t}) + \delta_{2} (T_{c} \times \text{After}_{t}) + \epsilon_{ct},
$$
where $\text{During}_{t}$ is a dummy for the all "during" periods and $\text{After}_{t}$ is a dummy for all "after" periods. Note, the reference period is the "before" period.  Feel free to switch the reference period to suit your needs. Remember, this works well when all units (cities) in your panel have well-defined before, during, and after treatment periods. The major advantage here is now instead of 80+ parameters to be estimated, you now have 6 (including the intercept).
To avoid any coding errors, I suggest creating one categorical "epoch" variable. It would look something like this:
epoch <- c("before", "during", "after)

This variable indexes the before, during, and after epoch for each city in your panel. Simply interact your treatment dummy with a factorized version of this variable and it will estimate each interaction and their constituent terms for free. It seems logical to force the "before" period to be the reference. But it doesn't have to work that way. You could also choose the "after" period to be your reference, where your estimate of $\delta_{1}$ is your  "before versus after" effect and your estimate of $\delta_{2}$ is your "during versus after" effect. It all depends on what comparisons are of substantive interest.
