Extracted main text — title through conclusion, appendix excluded. This is what our citation measures are computed over, published so the extraction can be checked by eye.
64,862 characters · 10 sections · 102 citation commands
Continuous difference-in-differences with double/debiased machine learning
\allowdisplaybreaks
Difference-in-differences (DiD) is one of the most widely used research designs in empirical work. While conventional DiD settings typically focus on binary or discrete multi-valued treatments, there is growing interest in extending DiD to continuous treatments. The motivation for continuous DiD is clear: the treatment group rarely receives interventions at a constant level, and treatment effects can vary with the intensity or “dose” of the treatment. Thus, rather than comparing treated and control groups before and after an intervention at an aggregate level, one can further investigate how outcomes vary across different treatment intensities within the treated group.
Continuous treatments are prevalent in many empirical settings. For instance, individuals may experience varying levels of exposure to policy interventions, marketing campaigns, or environmental pollutants, all of which can be modeled as continuous treatments. Several recent studies have explored DiD with continuous treatments, including ZDS2022 on the impact of shutting down online advertising sites, CJLR23 on racial discrimination in public accommodations, and AGHP22 on the effects of the expanded child tax credit.
Despite its widespread use in empirical research, the theoretical foundation for continuous DiD remains relatively underdeveloped, particularly in comparison to the extensive body of literature on DiD with binary or discrete treatments (see RSBP23, DD2023, CALL2023 for recent overviews). A few recent studies have begun to bridge this gap, notably CDPV2022, DHS2021, and CGS2024. For instance, DHS2021 extend the change-in-changes model of AI2006 to accommodate continuous treatments, while CDPV2022 examine the average slope of stayers in the continuous DiD setting. Our paper is closely related to CGS2024, which studies continuous DiD in the commonly used two-way fixed effect (TWFE) regression framework. CGS2024 demonstrate that, under TWFE, the regression parameter of interest can be decomposed as weighted integrals of either the average treatment parameters across treatment intensities with potentially negative weights or average causal responses with selection bias but nonnegative weights. They also provide data-driven non-parametric estimators for these causal parameters that are rate optimal.
In this paper, we focus on the average treatment effect on the treated (ATT) for any given continuous treatment intensity. Although this parameter is one of several investigated in CGS2024, our primary contribution is to incorporate covariates non-parametrically into both the identification and estimation procedures. Specifically, we modify the parallel trends assumption in CGS2024 by conditioning on covariates in a manner analogous to the “conditional parallel trends" assumption used in DiD for binary or discrete treatments; see HIT97,HIST98, Abadie2005, Chang2020, and SZ2020, for example. As noted in Abadie2005, an unconditional parallel trends assumption can be restrictive if covariates that influence outcome dynamics have different distributions across treatment and control groups. By conditioning on such covariates, we obtain a more robust framework for identifying and estimating the ATT in continuous treatment settings.
We first establish identification results analogous to those in Abadie2005, adapted to the continuous treatment setting. Based on these identification results, a naive estimator for the ATT can be constructed in two steps. First, one estimates several nuisance parameters from the identification results, including the conditional density of the continuous treatment. In the second step, the nuisance estimates are substituted into a simple average to obtain the estimator of the causal parameter. However, for potentially high-dimensional controls, while one may employ machine learning methods to estimate the nuisance parameters, doing so can introduce substantial bias in the causal parameter estimation (see CCDDHNR and the references therein). Moreover, reusing the same sample for both nuisance and causal parameter estimation can result in additional overfitting bias. To address these concerns, we adopt the double/debiased machine learning (DML) framework studied in CCDDHNR, which uses orthogonalization and cross-fitting to reduce the influence of nuisance parameter estimation on causal estimates.
Previous studies have adopted similar strategies in related settings. For instance, Chang2020 considers the DML framework for DiD with binary or discrete treatments, and SZ2020 proposes efficient doubly robust estimators for DiD with binary treatment. We contribute to and extend this literature to the continuous treatment setting. In particular, in place of the usual propensity score for the treated group, our setting requires the conditional density of the continuous treatment, which poses additional difficulties for directly applying DML methods, often involving only conditional mean functions as the nuisance parameters. To circumvent this, we introduce an approximate causal parameter $ATT_h$ using a kernel function. As the kernel bandwidth shrinks, $ATT_h$ converges to the true $ATT$. Importantly, by focusing on $ATT_h$, we can replace the conditional density with a conditional mean, which allows us to apply the existing DML results. We then derive orthogonal scores for both panel and repeated cross-sectional cases and construct corresponding DML estimators. Building on CCK2014a, CCK2014b, CCDDHNR, and FHLZ22, we establish the asymptotic normality of these estimators and show that the asymptotic bias becomes negligible under an appropriate undersmoothing kernel bandwidth. Additionally, we provide consistent variance estimators via cross-fitting and develop uniform confidence bands for the treatment curve using a multiplier bootstrap procedure. The results from our carefully designed simulation studies suggest that our estimators perform well.
To illustrate the usefulness of our method, we revisit AF2008, which examines the impact of the 1983 Medicare payment system (PPS) reform on the healthcare industry. Since the PPS reform affected hospitals with varying proportions of Medicare inpatients differently, the share of Medicare inpatients can be interpreted as a continuous treatment variable. This makes AF2008 an exemplary case for applying our methods. Thus, we non-parametrically estimate the ATTs of the PPS reform in a continuous DiD context, providing a more detailed understanding of the effects of this policy reform. In particular, contrasting with the linear estimates from AF2008, our results suggest significant heterogeneity in the impact of the PPS reform across hospitals with different shares of Medicare inpatients.
We note that the kernel smoothing has been previously considered in the causal inference literature with continuous treatment. For example, KMMS17 studies average potential outcomes under a continuous treatment, proposing a doubly robust signal and a two-step estimation procedure involving a pseudo-outcome and local kernel linear regression. Along similar lines, SC2021 employs series methods to establish uniform asymptotic results. HLM25 recently adopted a similar framework as KMMS17 to establish identification and estimation results on the average dose effect on treated. This causal parameter differs from ours in that it relies on different sets of parallel trends assumptions and it is an average dose-response on the entire treated group, akin to the average potential outcome. It is important to emphasize that while KMMS17 and HLM25 also employ the kernel techniques, their approach differs from ours in non-trivial ways. We use kernels primarily to approximate the original causal parameters, facilitating the construction of orthogonal scores, after which the final estimation proceeds as a simple average; see BL2017 for a more general discussion on this method. This contrasts with their approach, which uses kernel regressions to estimate the conditional mean of a pseudo-outcome. In this respect, our work is also related to KZ2018, SUZ19, and CYYL25, all of which consider continuous treatments and employ kernel-based moment functions to study the average potential outcomes and partial effects.
The remainder of this paper is organized as follows. Section 2 introduces continuous DiD and demonstrates the identification of the causal parameter. Section 3 provides the orthogonal scores. In Section 4, we present our estimators and establish their asymptotic properties. Section 5 showcases the simulation results, followed by a detailed empirical example in Section 6. Section 7 concludes.
\setcounter{equation}{0} In this section, we formally set up the difference-in-differences with continuous treatment following Abadie2005 and CGS2024. First, using the potential outcome notation (e.g. Rubin74), let $Y_{i,t}(0)$ denote the potential outcome of individual $i$ in period $t$ when receiving no treatment, and similarly let $Y_{i,t}(d)$ denote the potential outcome of individual $i$ in period $t$ when receiving treatment with intensity $d$.
The treatment variable $D$ is modeled as a random variable with a mixture distribution: a probability mass at $0$ and a continuous distribution on an interval $[d_L,d_H]$ excluding $0$. Specifically, the control group consists of individuals who receive treatment $D=0$, and we need a relatively large number of individuals in the control group so that the comparison with the treated is meaningful. On the other hand, the treated individuals can receive varied treatments, each with a potentially different treatment dose/intensity $D=d\in[d_L,d_H]$. We restrict our attention to the two-period $(t-1, t)$ models and suppress the time notation in treatment $D_i$ in the panel setting. Let $X_i$ denote the set of individual-level covariates. We make the following assumptions:
Assumptions (ref) and (ref) are analogous to those in the DiD literature with a discrete treatment. While Assumption (ref) requires a balanced panel, Assumption (ref) allows for repeated cross-sections but imposes stationarity of $(D,X)$ and hence rules out compositional changes.\footnote{For DiD with compositional changes, see HONG2013, ZIMMERT2020, and SX2025 for detailed discussions in the discrete treatment setting, and HHZ2024 in the continuous treatment setting.} Assumption (ref) is the strong overlap assumption, ensuring sufficient support for both treated and untreated individuals, which is crucial for identification. Assumption (ref) formalizes the requirement that there is no anticipated treatment effect prior to the treatment. Assumption (ref), a generalization of the discrete case in HIT97,HIST98, is the key identifying condition for the causal parameter. This assumption essentially states that, conditional on covariates, the unobserved counterfactual trend of the treated at each given treatment intensity is the same as the observed trend of the control group.
Next, we describe our target parameter. The causal parameter we are interested in is the average treatment effect on the treated (ATT for short) at any given treatment intensity $d\in [d_L,d_H]$:
The interpretation of this parameter is analogous to the cases with discrete treatment variables: the expected effect of treatment with intensity $d$ for those who actually received treatment with intensity $d$. See also CGS2024 Section 3 for a comprehensive discussion on ((ref)) and an alternative parallel trends assumption under which the average treatment effect $ATE(d) := E[Y_t(d) - Y_t(0)]$ can be identified. The following theorem presents the main results of this section, in which we establish the identification of $ATT(d)$ for both panel and repeated cross-sectional settings.
With Theorem (ref), one can build estimators for $ATT(d)$ using the estimated sample analogs. For potentially high-dimensional covariates, machine learning methods can be employed to estimate the nuisance parameters, including the conditional density $f_{D|X}(d|X)$ and the conditional probability $P(D=0|X)$. However, the use of machine learning methods can often result in non-trivial first-order biases in the estimation of the causal parameter, see e.g. CCDDHNR and references therein for a detailed discussion. Therefore, we consider alternative estimating equations that reduce the influence of the nuisance parameters.
\setcounter{equation}{0} In this section, we focus on the panel case for illustration as the repeated cross-sectional case only requires minor modifications. We begin by introducing Neyman orthogonality. Let $\theta_0(d) \in \Theta\subset \mathbf{R}$ be the low-dimensional parameter of interest, e.g., $ATT(d)$, and let $\rho_0(d) \in\mathcal{H}(d)$ denote the true low-dimensional nuisance parameters, e.g., $\rho_0(d) = f_D(d)$. The true infinite-dimensional nuisance parameters $\eta_0(d)\in\mathcal{T}(d)$ include $f_{D|X}(d|X)$ and $P(D=0|X)$ with the estimated $\hat\eta(d)$ in the realization set $T_N(d)\subset \mathcal{T}(d)$ with high probability.\footnote{New infinite-dimensional nuisance parameters can arise when constructing the orthogonal scores. We also explicitly index the nuisance parameters and nuisance function spaces by treatment intensity $d$.} Let $Z$ be the observable random vector, e.g. $Z =(Y_{t-1},Y_t, D, X)$ in the panel setting, and let $\psi: (Z,\theta(d),\rho(d),\eta(d))\mapsto \mathbf{R}$ denote a score function.\footnote{We say $\psi$ is a score function if at the true nuisance parameters $(\rho_0(d),\eta_0(d))$ and the true $\theta_0(d)$, the moment condition $E[\psi(Z,\theta_0(d),\rho_0(d),\eta_0(d)] = 0$ holds.} With these notations, following CCDDHNR and Chang2020, we formally define the Neyman orthogonality with respect to the infinite-dimensional nuisance parameters.
In the above definition, (a) says that $\psi$ identifies the parameter of interests while (b) ensures the first-order bias from estimating the infinite-dimensional nuisance parameters is zero. Recall that in the panel case,
where $\Delta Y:= Y_t - Y_{t-1}$. First, given the continuous nature of the treatment intensity, $\theta_0(d)$ cannot be estimated non-parametrically at root-$N$ rate. This relates to a class of non-regular parameters involving continuous treatment variables; see GW2015, KMMS17, SUZ19, SC2021, FHLZ22, and CYYL25 for example. Moreover, a score based on the above expression does not satisfy Neyman orthogonality, and an adjustment term has to be added.
To this end, we approximate the non-regular $ATT(d)$ with a family of smoothed regular parameters that are tractable. We note that this approach has been discussed extensively in BL2017 and CYYL25, and specifically we rely on the following observation (e.g., FYT96):
where $K(\cdot)$ is a kernel function. Replacing $E[\Delta Y|D=d]$ and $f_{D|X}(d|x)$ by their kernel counterparts, we can define $ATT_h(d)$ as follows:
which is an expression that consists of only conditional expectations. Notably, it can be shown that
which suggests that we can work with $ATT_h(d)$ instead. In particular, define the bias $B_h(d):= ATT_d - ATT_h(d)$, one can show that $B_h(d) = O(h^2)$, and we defer the formal result to the next section. For notation simplicity, we now formally define $ATT_h(d)$ in both settings.
Our goal is to construct scores that satisfy Neyman orthogonality for each $h$, and then take the limit as $h\to 0$. The next lemma presents such scores. To simplify the expressions, denote: $g(X) := P(D=0|X)$; $f_h(d|X):= E[K_h(D-d)|X] $; $\mathcal{E}_{\Delta Y}(X) := E[\Delta Y|X,D=0]$; $\mathcal{E}_{\lambda Y}(X) := E\big[\frac{T-\lambda}{\lambda(1-\lambda)}Y\big|X, D=0\big]$ with $\lambda = P(T=1)$.
The proof is provided in the appendix, where we construct the adjustment term and verify the Neyman orthogonality conditions from Definition (ref). We also provide an alternative derivation showing $\psi_h^{(1)}$ and $\psi_h^{(2)}$ as the efficient influence functions for the smoothed parameter $ATT_h(d)$ using the method proposed in HDDV22. The assumption on the existence of integrable functions $M_h^{(1)}$ and $M_h^{(2)}$ is mild and it justifies interchanging expectation and differentiation. For simplicity, we omit superscripts on $\psi^{(1)}_h$ and $\psi^{(2)}_h$ whenever the context is clear. The infinite-dimensional nuisance parameters in these new scores include $f_h(d|X)$, $g(X)$, $\mathcal{E}_{\Delta Y}(X)$, and $\mathcal{E}_{\lambda Y}(X)$, with the latter two introduced by the adjustment terms. Notably, the estimating moments for $ATT_h(d)$ based on these orthogonal scores remain robust to the first-order biases introduced by the nuisance estimates. In the next section, we construct DML estimators of $ATT(d)$ using these scores and establish their asymptotic properties.
\setcounter{equation}{0} As mentioned in the introduction, constructing DML estimators involves two main steps. In the previous section, we established scores that satisfy Neyman orthogonality (Lemma (ref)). These scores are then used alongside a cross-fitting procedure, further reducing estimation bias. With these key components in place, we construct DML estimators following the procedure proposed by CCDDHNR.
First, we partition the sample $I_N$ into $K\geq 2$ disjoint subsets $\{I_k\}_{k=1}^K$ of equal size $n=N/K$. For each $k\in\{1,\cdots,K\}$, we use the auxiliary sample $I_k^c:= I_N\setminus I_k$ to estimate the nuisance parameters. We then compute sample averages according to ((ref)) and ((ref)) using these estimates, evaluated at $I_k$, to obtain $\widehat{ATT}_k(d)$. Finally, we average across the $K$ estimates to obtain the final estimator $\widehat{ATT}(d)$. We note that at each $k=1,\cdots, K$, the nuisance parameters and $\widehat{ATT}_k(d)$ are estimated using disjoint subsamples, which reduces the overfitting bias and significantly simplifies the asymptotic analysis. Moreover, since $K$ is fixed, it does not affect the asymptotic properties of the estimator. In practice, we recommend using $K=5$ as a rule of thumb and leave the optimal choice of $K$ to future research. The detailed algorithms are deferred to Appendix A.
Next, we outline the regularity conditions required to establish the asymptotic properties of our DML estimators. We focus on the panel case and present the analogous results for the repeated cross-sections in Appendix B. For notational simplicity, let $\mathcal{D}$ denote a closed sub-interval of $(d_L, d_H)$ whose boundary points can be chosen arbitrarily close to $d_L$ and $d_H$, and let $\mathcal{X}$ and $\Delta\mathcal{Y}$ denote the supports of $X$ and $\Delta Y$, respectively.
The kernel function is central to our analysis. In addition to its well-established theoretical properties for estimating the density $f_D(d)$, we also use it to approximate the point mass at $D=d$ and the conditional density $f_{D|X}(d|X)$. Assumption (ref) imposes the standard regularity conditions on the kernel function, which are essential for establishing the asymptotic normality of our estimator. Assumption (ref) requires smoothness and boundedness of the outcome variable and relevant distributions, while Assumption (ref) specifies conditions on the kernel bandwidth and the quality of the non-parametric nuisance estimators.
The following lemma characterizes the bias of using kernels to approximate $ATT(d)$.
The proof is given in the companion supplement. This lemma suggests that, for an undersmoothing bandwidth, the bias does not affect the asymptotic distribution of our estimators. The next theorem is the main result of this section that establishes the asymptotic normality of our estimator for $ATT(d)$.
The proof builds on the DML framework of CCDDHNR, modified to accommodate kernel smoothing. The asymptotic variance has two components, both depending inversely on the kernel bandwidth $h$: one arising from the kernels in the orthogonal score $\psi_h$, and the other from the linear expansion of the estimator with respect to the kernel density estimator $\hat{f}_D(d)$. Since $h$ is a function of the sample size $N$ under our assumptions, we index the asymptotic variance by $N$ to reflect this dependence. Therefore, our estimator $\widehat{ATT}(d)$ attains a convergence rate of $\sqrt{Nh}$, which, though slower than the parametric rate $\sqrt{N}$, is comparable to the optimal rate for one-dimensional non-parametric regression estimation.
Next, following CCDDHNR and Chang2020, we consider a cross-fitted variance estimator. For notation simplicity, denote $\hat{\theta}_h(d):= \widehat{ATT}(d)$ and $E_{n,k} f(Z_i):= n^{-1}\sum_{i\in I_k}f(Z_i)$ as the empirical average of a function $f$ evaluated at $Z_i$'s in the subsample $I_k$. For the panel case, define
Then, with this variance estimator, the $1-\alpha$ confidence interval can be constructed as $[\widehat{ATT}(d) - z_{1-\alpha/2}\hat{\sigma}_N(d)/\sqrt{N}, \widehat{ATT}(d) + z_{1-\alpha/2}\hat{\sigma}_N(d)/\sqrt{N}]$ where $z_{1-\alpha/2}$ denotes the $1-\alpha/2$-th quantile of the standard normal random variable. The following theorem establishes the consistency of the cross-fitted variance estimator.
Alternatively, we can consider a multiplier bootstrap procedure to construct confidence intervals. Such procedure has been discussed extensively in recent studies, see, e.g., CCK2014b, BCFH17, SUZ19, CJ21, FHLZ22, and CYYL25. First, we make the following assumption on the multiplier.
In practice, let $\{\xi_i\}_{i=1}^N$ be an i.i.d. sequence of random variables that satisfies Assumption (ref). Then for each $b=1,\cdots, B$, we independently draw such a sequence $\{\xi_i\}_{i=1}^{N}$ and construct estimates based on the following expression. For the panel case, define
Let $\hat{c}_{\alpha}$ denote the $\alpha$-th quantile of $\{\widehat{ATT}(d)_{b}^*- \widehat{ATT}(d)\}_{b=1}^B$, a $1-\alpha$ confidence interval can be constructed as $[\widehat{ATT}(d) - \hat{c}_{1-\alpha/2}, \widehat{ATT}(d) -\hat{c}_{\alpha/2}]$.
Moreover, we can establish valid uniform inference results based on the bootstrap estimator proposed here. The following assumption strengthens Assumption (ref).
This assumption differs from the pointwise case in two key ways. First, we require that, uniformly over $\mathcal{D}$, the nuisance estimator $\hat{f}_h(d|X)$ remains bounded and has rate $h^{-1/2}\varepsilon_N$. Second, we assume that the estimated density and conditional density to have bounded derivatives with probability tending to $1$, ensuring that the score functions are Lipschitz continuous on $\mathcal{D}$. These additional assumptions are mild and can be enforced during estimation procedures. With these modified assumptions, the linear expansion of the bootstrap estimators holds uniformly over $d\in\mathcal{D}$.
This theorem is the basis for establishing uniform inference theory using the multiplier bootstrap estimator. We consider the following procedure, see CCK2014b and FHLZ22 for example, to establish valid uniform confidence bands.
With Assumption (ref), we can easily adapt our proof of Theorem (ref) to establish the uniform consistency of our cross-fitted variance estimator ((ref)) over $\mathcal{D}$. Then, with Theorem (ref), we can show that the proposed uniform confidence band achieves asymptotic coverage of $1-\alpha$, using results from CCK2014a (Proposition 3.2 and Theorem 3.2) and CCK2014b (Corollary 3.1). Since this argument is well established in the literature, e.g., see the discussion of Theorem 4.2 in FHLZ22, we do not include the formal theoretical discussion here. Instead, we focus on presenting the new results in Theorem (ref) and defer its proof to the appendix.
\setcounter{equation}{0} Data-generating process (a) $p = 100$ dimensional covariates $X \sim N(0.2,\Sigma)$, where $\Sigma$ has variances $1$ on the diagonal and covariances $0.1$ off-diagonal; (b) the control group propensity score follows $P(D=0|X) = 1/(1+\exp(-X'\gamma))$, with $\gamma_j = 0.5j^{-2}$; (c) for $D>0$, the continuous treatment is generated as $D= (1+\exp(X'\alpha))^{-1} + V$, where $V \perp X$, $V \sim Beta(2,2)$, $\alpha_j = 0.3j^{-2}$; (d) the potential outcomes are given by $Y_{t-1}(0) = \epsilon_1$, $Y_t(0) = Y_{t-1}(0) + X'\beta + 1 + \epsilon_2$, $Y_t(D) = Y_t(0) - 0.5D^2 + \epsilon_3$, where $\beta_j = 0.5/j$ for $j = 1, \cdots, 6$ and $0$ otherwise, and $(\epsilon_1, \epsilon_2, \epsilon_3)\sim N(0, I_3)$. For the panel setting, the generated data are $(Y_{i,t-1}, Y_{i,t}, X_i, D_i)$, with $Y_{t-1} = Y_{t-1}(0)$ and $Y_t = 1\{D>0\}Y_t(D) + 1\{D=0\}Y_t(0)$. Additionally, for the repeated cross-sectional setting, the generated data are $(Y_i, T_i, X_i, D_i)$, with time indicator $T \sim \text{Bern}(0.5)$ and $Y = TY_t + (1-T)Y_{t-1}$, $Y_{t-1} = Y_{t-1}(0)$, $Y_t = 1\{D>0\}Y_t(D) + 1\{D=0\}Y_t(0)$.
In our simulations, the nuisance parameters $P(D=0|X)$, $f_h(d|X) = E[K_h(D-d)|X]$, $E[Y_t - Y_{t-1}|X,D=0]$, and $E\big[\frac{T-\lambda}{\lambda(1-\lambda)}Y\big|X,D=0\big]$ are estimated non-parametrically using random forests each with 200 trees of maximum depth 20. Throughout our simulations, we also use an undersmoothing kernel bandwidth $h = 1.06\hat{\sigma}_{\tilde{D}}N^{-1/4}$, where $\hat{\sigma}_{\tilde{D}}$ is the estimated standard deviation of positive treatment intensities. We consider sample sizes $N = 2000$ and $10000$ for both panel and repeated cross-sectional settings, and we conduct $B = 500$ simulations in each setting. The DGP implies the true $ATT(d) = -0.5d^2$, and we focus on a specific treatment intensity $d = 0.9$. Notably, the continuous treatment variable is dependent on the correlated high-dimensional covariates in a nonlinear way. Additionally, the DGPs suggest that the effective sample size should be small at the target intensity, which adds another layer of difficulty for estimation.
Despite these challenges, the simulation results suggest that our estimators perform well. The histograms of these simulation estimates are shown in Figure (ref), where the red lines indicate the true ATT. We see that as the sample size increases, both bias and variance decrease. The simulation estimates appear to follow a normal distribution in each case, which is consistent with our asymptotic theory. Moreover, in Table (ref), we report the bias, the standard deviation of estimated ATTs (Std), the root-mean-squared error (RMSE), the average standard deviations (AVSE), and the coverage probability of 95 percent confidence intervals. In both settings, bias, standard deviation, and RMSE decrease as the sample size increases. The standard deviations of the simulation estimates are very close to the average estimated standard errors, suggesting that our variance estimators perform well. The coverage of the estimated confidence intervals is close to 95 percent, although there is a slight under-coverage in the panel setting.
\setcounter{equation}{0}
The Medicare Prospective Payment System (PPS) reform, introduced in 1983, shifted Medicare hospital reimbursements from a full-cost model to a fixed payment per diagnosis. However, for the first three years, capital costs continued to be reimbursed based on actual expenses.\footnote{As noted in AF2008, Medicare’s capital cost reimbursements remained unchanged until 1991 due to delays.} This created a relative increase in labor costs for hospitals treating Medicare inpatients. AF2008 highlights this feature, showing that the PPS reform significantly increased hospitals’ capital-labor ratios and encouraged technology adoption.
Theoretically, AF2008 predicts that PPS reform would lead to a higher capital-labor ratio and, if capital-labor substitution is sufficiently elastic, an increased demand for capital and technology. Since only hospitals with Medicare inpatients were affected, these effects likely varied with Medicare inpatient share. To test these predictions, AF2008 uses data from the 1980–1986 Annual American Hospital Association (AHA) survey, which provides hospital information including expenditures, employment, and technology adoption. Their baseline specification is a linear regression:
where $Y_{i,t}$ is the capital-labor ratio or total number of medical facilities for hospital $i$ in year $t$, $D_i$ is the pre-reform Medicare inpatient share, and $\text{Post}_t$ is a treatment-timing indicator. $X_{i,t}$ represents covariates, and $\alpha_i$ and $\gamma_t$ are hospital and year fixed effects, respectively. AF2008 argues that $\beta$ captures the causal effect of PPS reform on capital-labor ratios and technology adoption, relying on a parallel trends assumption: in the absence of the PPS reform, hospitals with different shares $D_i$ should have experienced similar changes in outcomes over time.
Recent work by CGS2024 examines the same empirical setting in detail and finds suggestive evidence that the parallel trends assumption may be too strong. This underscores the importance of incorporating covariates to improve the plausibility of the identifying assumption. By conditioning on covariates, our approach refines the parallel trends assumption, ensuring that hospitals are compared based on more similar characteristics. In this way, our analysis complements CGS2024, offering an alternative perspective on the effects of the PPS reform.
Regression (ref) resembles a Two-Way Fixed Effects (TWFE) design but differs in that $D_i$ is continuous. As shown by CGS2024, with continuous treatment, the coefficient $\beta$ in (ref) can be viewed as a weighted average of $ATT(d)$ with possible negative weights, which complicates interpretation.\footnote{See Proposition 10 in CGS2024. They do not incorporate covariates, but the issue persists.} Our continuous DiD framework addresses this by reframing AF2008’s design as follows: \setcounter{bean}{0}
We identify the causal effect at intensity $d$ as: \[ ATT(d) = E[Y_t(d) - Y_t(0)|D=d]. \] Unlike the constant $\beta$ in ((ref)), the causal effect curve $ATT(d)$ can be used to study the policy impact at a much more granular level. For example, if the PPS reform raised the capital-labor ratio, $ATT(d)$ should be positive for all $d>0$. Moreover, $ATT(d)$ should increase in $d$ if the impact of PPS reform is larger for hospitals with higher shares of Medicare inpatients. We apply our panel estimator and, for comparability with AF2008, average pre-treatment outcomes ($Y_{t-1}$) over 1980–1983 and post-treatment outcomes ($Y_t$) over 1984–1986 (capital-labor ratio) or 1984–1985 (technology adoption). Our data source is the cleaned data file from AF2008.
First, we examine the results for capital-labor ratio. All estimated ATTs are positive, mirroring the findings in AF2008 and suggesting that the PPS reform led to an increase in the capital-labor ratio. For comparison, AF2008 reports an estimate of $1.13$, which exceeds most of our estimates. Moreover, our estimates vary across treatment intensities and do not exhibit a strictly increasing trend, contradicting the theoretical prediction that hospitals with higher Medicare shares would see greater increases in the capital-labor ratio. At low and high treatment intensities, we note that the small effective sample sizes lead to noisier estimates, as reflected in the wider confidence intervals. For completeness, an effect curve estimated without covariates using a kernel method shows a similar pattern to our DML estimates.
Next, we present evidence of increased technological adoption following the PPS reform. Specifically, we consider the total number of specialized medical facilities per hospital as a proxy for technological adoption. All of our estimated ATTs for this outcome are positive, aligning with AF2008’s prediction that the PPS reform would incentivize technological adoption. The estimated treatment curve initially rises with treatment intensity but then declines at higher intensities, again diverging from the theoretical prediction that hospitals with larger Medicare shares would invest more. For comparison, an effect curve estimated using a kernel method without covariates again shows a similar pattern to our DML estimates. As with the capital-labor ratio, estimates are especially noisy where data are sparse, an issue amplified by our undersmoothing bandwidth.
This paper studies difference-in-differences models with continuous treatments. Our identification results are based on a conditional parallel trends assumption, allowing researchers to account for covariates non-parametrically. Under the double/debiased machine learning framework, we develop non-parametric estimators for the average treatment effect on the treated at each continuous treatment intensity and establish their asymptotic properties. Monte Carlo simulations demonstrate that our estimators perform well despite the highly non-linear relationship between the continuous treatment and the high-dimensional covariates. To demonstrate the empirical relevance of our methodology, we re-examine the research questions posed in AF2008 by applying our estimators to their dataset and obtaining new empirical insights. The extension of difference-in-differences models to the continuous treatment setting has important implications for empirical research. Our methods provide researchers with new tools for examining the impacts of continuous treatment variables.