The exact contents of citations.db main_text.text for this paper — one flattened LaTeX string, title through conclusion, appendix excluded, unmodified except for removing email addresses. This is what our citation measures are computed over.
64,862 characters
Continuous difference-in-differences with double/debiased machine learning
\maketitle
\allowdisplaybreaks
\begin{abstract}
This paper extends difference-in-differences to settings with continuous treatments. Specifically, the average treatment effect on the treated (ATT) at any level of treatment intensity is identified under a conditional parallel trends assumption. Estimating the ATT in this framework requires first estimating infinite-dimensional nuisance parameters, particularly the conditional density of the continuous treatment, which can introduce substantial bias. To address this challenge, we propose estimators for the causal parameters under the double/debiased machine learning framework and establish their asymptotic normality. Additionally, we provide consistent variance estimators and construct uniform confidence bands based on a multiplier bootstrap procedure. To demonstrate the effectiveness of our approach, we revisit a previous study on the 1983 Medicare Prospective Payment System reform, reframing it as a DiD with continuous treatment and non-parametrically estimating its effects.
\textbf{Keywords:} Difference-in-differences, causal inference, continuous treatment, machine learning
\end{abstract}
\newpage
\section{Introduction}
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 \cite{ZDS2022} on the impact of shutting down online advertising sites, \cite{CJLR23} on racial discrimination in public accommodations, and \cite{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 \cite{RSBP23}, \cite{DD2023}, \cite{CALL2023} for recent overviews). A few recent studies have begun to bridge this gap, notably \cite{CDPV2022}, \cite{DHS2021}, and \cite{CGS2024}. For instance, \cite{DHS2021} extend the change-in-changes model of \cite{AI2006} to accommodate continuous treatments, while \cite{CDPV2022} examine the average slope of stayers in the continuous DiD setting. Our paper is closely related to \cite{CGS2024}, which studies continuous DiD in the commonly used two-way fixed effect (TWFE) regression framework. \cite{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 \cite{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 \cite{CGS2024} by conditioning on covariates in a manner analogous to the ``conditional parallel trends" assumption used in DiD for binary or discrete treatments; see \cite{HIT97,HIST98}, \cite{Abadie2005}, \cite{Chang2020}, and \cite{SZ2020}, for example. As noted in \cite{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 \cite{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 \cite{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 \cite{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, \cite{Chang2020} considers the DML framework for DiD with binary or discrete treatments, and \cite{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 \cite{CCK2014a, CCK2014b}, \cite{CCDDHNR}, and \cite{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 \cite{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 \cite{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 \cite{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, \cite{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, \cite{SC2021} employs series methods to establish uniform asymptotic results. \cite{HLM25} recently adopted a similar framework as \cite{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 \cite{KMMS17} and \cite{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 \cite{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 \cite{KZ2018}, \cite{SUZ19}, and \cite{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.
\section{Setup and Identification}
\setcounter{equation}{0}
In this section, we formally set up the difference-in-differences with continuous treatment following \cite{Abadie2005} and \cite{CGS2024}. First, using the potential outcome notation (e.g. \cite{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:
\begin{assumption}[Panel]\label{didas1}
The observed data $\{Y_{i,t-1}, Y_{i,t}, D_i, X_i\}_{i=1}^N$ are independently and identically distributed.
\end{assumption}
\begin{assumption}[Repeated Cross-Sections]\label{didas2} (a) For each individual $i$ in the pooled sample, $T_i$ is a time indicator $=1$ if observation $i$ belongs to the post-treatment sample and $=0$ otherwise, and $Y_i = (1-T_i)Y_{i,t-1} + T_iY_{i,t}$; (b) $(D,X)\perp T$ and the following holds: (i) conditional on $T=0$, data are i.i.d. from the distribution of $(Y_{t-1},D,X)$; (ii) conditional on $T=1$, data are i.i.d. from the distribution of $(Y_t,D,X)$.
\end{assumption}
\begin{assumption}[Support]\label{didas3} (a) The support of $D$ is $\{0\}\sqcup[d_L, d_H]$ with $0<d_L<d_H<\infty$; (b) there exists a constant $0<\kappa <\frac{1}{2}$ such that, almost surely, $\kappa <P(D=0|X)<1-\kappa$ and $f_{D|X}(d|X) > \kappa$ for all $d\in [d_L, d_H]$.
\end{assumption}
\begin{assumption}[No Anticipation]\label{didas5} $Y_t = Y_t(D)$, $Y_{t-1} = Y_{t-1}(0)$.
\end{assumption}
\begin{assumption}[Conditional Parallel Trends]\label{didas4}
For all $d\in [d_L, d_H]$, the following holds
\begin{align}
E[Y_t(0)-Y_{t-1}(0)|X, D=d] = E[Y_{t}(0)-Y_{t-1}(0)|X, D=0].
\end{align}
\end{assumption}
Assumptions \ref{didas1} and \ref{didas2} are analogous to those in the DiD literature with a discrete treatment. While Assumption \ref{didas1} requires a balanced panel, Assumption \ref{didas2} allows for repeated cross-sections but imposes stationarity of $(D,X)$ and hence rules out compositional changes.\footnote{For DiD with compositional changes, see \cite{HONG2013}, \cite{ZIMMERT2020}, and \cite{SX2025} for detailed discussions in the discrete treatment setting, and \cite{HHZ2024} in the continuous treatment setting.} Assumption \ref{didas3} is the strong overlap assumption, ensuring sufficient support for both treated and untreated individuals, which is crucial for identification. Assumption \ref{didas5} formalizes the requirement that there is no anticipated treatment effect prior to the treatment. Assumption \ref{didas4}, a generalization of the discrete case in \cite{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 \textit{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) \textit{at any given treatment intensity} $d\in [d_L,d_H]$:
\begin{equation}\label{att}
ATT(d) := E[Y_t(d) - Y_t(0)|D=d].
\end{equation}
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 \cite{CGS2024} Section 3 for a comprehensive discussion on (\ref{att}) 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.
\begin{theorem}[Identification of ATT]\label{thm:cdid1}
(a) (Panel) If Assumptions \ref{didas1}, \ref{didas3}, \ref{didas5}, and \ref{didas4} hold, then, for any $d\in[d_L, d_H]$,
\begin{align}
ATT(d) = E[Y_t-Y_{t-1}| D=d] - E\bigg[(Y_t-Y_{t-1})\mathbf{1}\{D=0\}\frac{f_{D|X}(d|X)}{f_{D}(d)P(D=0|X)}\bigg];
\end{align}
(b) (Repeated Cross-Sections) if Assumptions \ref{didas2}, \ref{didas3}, \ref{didas5}, and \ref{didas4} hold, then, for any $d\in[d_L, d_H]$,
\begin{align}
ATT(d) = E\bigg[\frac{T-\lambda}{\lambda(1-\lambda)}Y \bigg| D=d\bigg]- E\bigg[\frac{T-\lambda}{\lambda(1-\lambda)}Y\mathbf{1}\{D=0\}\frac{f_{D|X}(d|X)}{f_D(d)P(D=0|X)}\bigg]
\end{align}
where $\lambda := P(T=1)$.
\end{theorem}
With Theorem \ref{thm:cdid1}, 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. \cite{CCDDHNR} and references therein for a detailed discussion. Therefore, we consider alternative estimating equations that reduce the influence of the nuisance parameters.
\section{Orthogonal Scores}
\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 \textit{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 \cite{CCDDHNR} and \cite{Chang2020}, we formally define the Neyman orthogonality with respect to the infinite-dimensional nuisance parameters.
\begin{definition}[Neyman Orthogonality]\label{neyman} A score $\psi$ satisfies the Neyman orthogonality condition at $(\theta_0(d),\rho_0(d),\eta_0(d))$ with respect to a nuisance realization set $T_N(d)\subset \mathcal{T}(d)$ if (a) $\theta_0(d)$ satisfies the moment condition
\begin{align}
E_P[\psi(Z,\theta_0(d),\rho_0(d),\eta_0(d))] = 0;
\end{align}
(b) for $r\in[0,1)$ and $\eta(d)\in T_N(d)$, the Gateaux (directional) derivative satisfies
\begin{align}
\partial_r E_P[\psi(Z,\theta_0(d),\rho_0(d),\eta_0(d) + r(\eta(d)-\eta_0(d)))]|_{r=0} = 0.
\end{align}
\end{definition}
In the above definition, (a) says that $\psi$ identifies the parameter of interests while (b) ensures the first-order bias from estimating the \textit{infinite-dimensional} nuisance parameters is zero. Recall that in the panel case,
\begin{align}
\theta_0(d) = ATT(d) = E[\Delta Y|D=d] - E\bigg[\Delta Y\mathbf{1}\{D=0\}\frac{f_{D|X}(d|X)}{f_D(d)P(D=0|X)}\bigg].
\end{align}
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 \cite{GW2015}, \cite{KMMS17}, \cite{SUZ19}, \cite{SC2021}, \cite{FHLZ22}, and \cite{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 \cite{BL2017} and \cite{CYYL25}, and specifically we rely on the following observation (e.g., \cite{FYT96}):
\begin{equation}\label{cd:kernel}
f_{D|X}(d|x) = \lim_{h\to 0} E[K_h(D-d)|X=x],\quad K_h(u) := \frac{1}{h}K\Big(\frac{u}{h}\Big)
\end{equation}
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:
\begin{align}\label{atth}
ATT_h(d) :=& E\bigg[\Delta Y\frac{K_h(D-d)}{f_D(d)}\bigg] - E\bigg[\Delta Y\mathbf{1}\{D=0\}\frac{E[K_h(D-d)|X]}{f_D(d)P(D=0|X)}\bigg] \notag \\
=& E\bigg[\Delta Y\frac{K_h(D-d)P(D=0|X) - \mathbf{1}\{D=0\}E[K_h(D-d)|X]}{f_D(d)P(D=0|X)} \bigg],
\end{align}
which is an expression that consists of only conditional expectations. Notably, it can be shown that
\begin{align*}
ATT(d) = \lim_{h\to 0} ATT_h(d),
\end{align*}
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.
\begin{definition}[Panel]
\begin{align}\label{eq:atth1}
ATT_h(d) = E\bigg[\Delta Y \frac{K_h(D-d)P(D=0|X) - \mathbf{1}\{D=0\}E[K_h(D-d)|X]}{f_D(d)P(D=0|X)} \bigg]
\end{align}
where $\Delta Y = Y_t - Y_{t-1}$.
\end{definition}
\begin{definition}[Repeated Cross-Sections]
\begin{align}\label{eq:atth2}
ATT_h(d) =& E\bigg[Y^\lambda\frac{K_h(D-d)P(D=0|X) - \mathbf{1}\{D=0\}E[K_h(D-d)|X]}{f_D(d)P(D=0|X)} \bigg]
\end{align}
where $Y^\lambda := \frac{T-\lambda}{\lambda(1-\lambda)}Y$ with $\lambda = P(T=1)$.
\end{definition}
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)$.
\begin{lemma}\label{lm:scores}
Define (a) for the panel setting,
\begin{equation}\label{psi1}
\psi_h^{(1)} := \frac{K_h(D-d)g(X) - \mathbf{1}\{D=0\}f_h(d|X)}{f_D(d)g(X)}\bigg(\Delta Y - \mathcal{E}_{\Delta Y}(X)\bigg) -ATT_h(d),
\end{equation}
and (b) for the repeated cross-sectional setting,
\begin{equation}\label{psi2}
\psi_h^{(2)} := \frac{K_h(D-d)g(X) - \mathbf{1}\{D=0\}f_h(d|X)}{f_D(d)g(X)} \bigg(\frac{T-\lambda}{\lambda(1-\lambda)}Y - \mathcal{E}_{\lambda Y}(X) \bigg) -ATT_h(d).
\end{equation}
Suppose there exist $M_h^{(1)}\in L^1(P_{Y_{t-1},Y_t,D,X})$ and $M_h^{(2)}\in L^1(P_{Y,T,D,X})$ such that $|\psi_h^{(1)}|\leq M_h^{(1)}$ and $|\psi_h^{(2)}|\leq M_h^{(2)}$ almost surely. Then the scores $\psi_h^{(1)}$ and $\psi_h^{(2)}$ satisfy Neyman orthogonality defined in (\ref{neyman}).
\end{lemma}
The proof is provided in the appendix, where we construct the adjustment term and verify the Neyman orthogonality conditions from Definition \ref{neyman}. 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 \cite{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.
\section{Estimation and Inference}
\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{lm:scores}). 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 \cite{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{psi1}) and (\ref{psi2}) 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.
\begin{assumption}[Kernel]\label{cdidas1} The kernel function $K(\cdot)$ satisfies: (a) $K(\cdot)$ is bounded and differentiable; (b) $\int K(u) du = 1$, $\int uK(u)du = 0$, $0<\int u^2 K(u) du <\infty$. Moreover, for notation simplicity, define $K_h(u):= h^{-1}K(u/h)$.
\end{assumption}
\begin{assumption}[Bounds and Smoothness, Panel]\label{cdidas2} (a) There exist constants $c>0$ and $0<C<\infty$ such that $\sup_{d\in\mathcal{D}} f_D^0(d)>c$, $|Y_{t-1}| < C$, $|Y_{t}| < C$, $c <f_h^0(d|X)<C$ $\forall d\in\mathcal{D}$, and $|\mathcal{E}_{\Delta Y}^0(X)|<C$ almost surely; (b) $f_D^0(d) \in C^2(\mathcal{D})$ and $\sup_{d\in\mathcal{D}}|\partial_d^2 f_D^{0}(d)| < \infty$; (c) $f_{D|X}^0(d|x) \in C^2(\mathcal{D})$ $\forall x\in \mathcal{X}$ and $\sup_{d,x \in \mathcal{D},\mathcal{X}} |\partial^2_d f_{D|X}^0(d|x)| < \infty$; (d) $f_{\Delta Y, D}(t, d) \in C^2(\Delta\mathcal{Y})$ and $\sup_{t,d \in \Delta\mathcal{Y},\mathcal{D}} |\partial_t^2 f_{\Delta Y, D}(t, d)| < \infty$.
\end{assumption}
\begin{assumption}[Rates, Panel]\label{cdidas3} (a) The kernel bandwidth $h = h_N\to 0$ satisfies $Nh\to\infty$ and $\sqrt{Nh^5} = o(1)$; (b) there exists a sequence $\varepsilon_N\to 0$ such that $h^{-1}\varepsilon_N^2 = o(1)$; (c) with probability tending to $1$, $\|\hat{f}_h(d|X) - f_h^0(d|X)\|_{P,2}\leq h^{-1/2}\varepsilon_N$, $\|\hat{g}(X) - g_0(X)\|_{P,2}\leq \varepsilon_N$, $\|\hat{\mathcal{E}}_{\Delta Y}(X) - \mathcal{E}_{\Delta Y}^0(X)\|_{P,2}\leq \varepsilon_N$; (d) with probability tending to $1$, $\kappa<\hat{g}(X) < 1-\kappa$ and $c <\hat{f}_h(d|X)<C$ almost surely, and $\|\hat{\mathcal{E}}_{\Delta Y}(X)\|_{P,\infty}<C$.
\end{assumption}
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{cdidas1} imposes the standard regularity conditions on the kernel function, which are essential for establishing the asymptotic normality of our estimator. Assumption \ref{cdidas2} requires smoothness and boundedness of the outcome variable and relevant distributions, while Assumption \ref{cdidas3} specifies conditions on the kernel bandwidth and the quality of the non-parametric nuisance estimators.
\begin{remark}
\textnormal{Assumption \ref{cdidas3} (a) and (b) give $h = o(N^{-1/5})$ and $h = \omega(\varepsilon_N^2)$. However, consistency of our variance estimator additionally requires $h^{-2}\varepsilon_N^2 + h^{-3}N^{-1} = o(1)$, which imposes a more restrictive lower bound on $h$. Moreover, whereas the standard DML literature assumes the nuisance estimators to converge at rate $\varepsilon_N = o(N^{-1/4})$, we allow the conditional density $\hat{f}_h$ to converge at a slower rate $h^{-1/2}\varepsilon_N$. This relaxation does not contradict the existing DML results for regular parameters; our target parameter is non-regular and cannot be estimated non-parametrically at $\sqrt{N}$ rate because of the continuous treatment. Finally, although our results assume a deterministic kernel bandwidth, they should extend to data-driven choices (for example, the adaptive procedure in \cite{BL2017}), which we leave to future work.}
\end{remark}
The following lemma characterizes the bias of using kernels to approximate $ATT(d)$.
\begin{lemma}[Bias of $ATT_h(d)$, Panel]\label{lm:bias}
Suppose Assumptions \ref{cdidas1}, \ref{cdidas2}, \ref{cdidas3} hold. Then $B_h(d) := ATT(d) - ATT_h(d)$ satisfies $B_h(d) = O(h^2)$ for any $d\in \mathcal{D}$.
\end{lemma}
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)$.
\begin{theorem}[Asymptotic Normality, Panel]\label{thm:asymp}\hspace{-2.9pt} Suppose assumptions \ref{didas1}, \ref{didas3}, \ref{didas5}, \ref{didas4}, \ref{cdidas1}, \ref{cdidas2}, and \ref{cdidas3} hold. Then, for $d\in\mathcal{D}$, if $\varepsilon_N = o(N^{-1/4})$,
\[
\frac{\widehat{ATT}(d) - ATT(d)}{\sigma_{N}(d)/\sqrt{N}}\quad \to^d\quad N(0,1)
\]
where
\begin{align}\label{var1}
\sigma_{N}^2(d):= E\bigg[\bigg(\psi_h^{(1)}(Z,\theta_{0h}(d),f_D^0(d),\eta_0(d)) - \frac{\theta_{0h}(d)}{f_D^0(d)}\big(K_h(D-d)-E[K_h(D-d)]\big)\bigg)^2\bigg]
\end{align}
for $\theta_{0h}(d):= ATT_h(d)$ defined in (\ref{eq:atth1}) and $\psi_h^{(1)}$ defined in (\ref{psi1}).
\end{theorem}
The proof builds on the DML framework of \cite{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 \cite{CCDDHNR} and \cite{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
\begin{equation}\label{asvar}
\hat{\sigma}_{N}^2(d) := \frac{1}{K}\sum_{k=1}^K E_{n,k}\bigg[
\bigg(
\psi_h^{(1)}(Z,\hat{\theta}_h(d),\hat{f}_k(d), \hat{\eta}_k(d)) - \frac{\hat{\theta}_h(d)}{\hat{f}_k(d)}\big(K_h(D-d)-\hat{f}_k(d)\big)
\bigg)^2
\bigg].
\end{equation}
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.
\begin{theorem}[Consistency of Variance Estimator, Panel]\label{thm:asvar}
Suppose the conditions of Theorem \ref{thm:asymp} hold and assume that $h^{-2}\varepsilon_N^2 + h^{-3}N^{-1} = o(1)$. Then, for $d\in\mathcal{D}$,
\begin{align*}
\hat{\sigma}_{N}^2(d) = \sigma_{N}^2(d) + o_p(1)
\end{align*}
where $\hat{\sigma}_{N}^2(d)$ is defined in (\ref{asvar}) and $\sigma_{N}^2(d)$ is defined in (\ref{var1}).
\end{theorem}
Alternatively, we can consider a multiplier bootstrap procedure to construct confidence intervals. Such procedure has been discussed extensively in recent studies, see, e.g., \cite{CCK2014b}, \cite{BCFH17}, \cite{SUZ19}, \cite{CJ21}, \cite{FHLZ22}, and \cite{CYYL25}. First, we make the following assumption on the multiplier.
\begin{assumption}[Sub-exponential Multiplier]\label{as:subexp}
The random variable $\xi$ satisfies: (a) $\xi$ has a sub-exponential distribution; (b) $E[\xi]=Var(\xi) = 1$; (c) $\xi$ is independent of $(Y_{t-1}, Y_t, D, X)$ for the panel case and independent of $(Y, T, D, X)$ for the repeated cross-sectional case.
\end{assumption}
In practice, let $\{\xi_i\}_{i=1}^N$ be an i.i.d. sequence of random variables that satisfies Assumption \ref{as:subexp}. 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
\begin{align}\label{bootstrap}
\widehat{ATT}(d)_{b}^* := \frac{1}{N}\sum_{k=1}^K\sum_{i\in I_k} \xi_i& \frac{K_h(D_i-d)\hat{g}_k(X_i) - \mathbf{1}\{D_i=0\}\hat{f}_{h,k}(d|X_i)}{\hat{f}_k(d)\hat{g}_k(X_i)}\notag \\
\times& \big(\Delta Y_i - \hat{\mathcal{E}}_{\Delta Y,k}(X_i)\big).
\end{align}
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{cdidas3}.
\begin{assumption}[Uniform Inference Rates, Panel]\label{cdidas4} \sloppy
(a) The kernel bandwidth $h = h_N\to 0$ satisfies $Nh\to\infty$ and $\sqrt{Nh^5} = o(1)$; (b) there exists a sequence $\varepsilon_N\to 0$ such that $h^{-1}\varepsilon_N^2 = o(1)$; (c) with probability tending to $1$, $\sup_{d\in\mathcal{D}}\|\hat{f}_h(d|X) - f_h^0(d|X)\|_{P,2}\leq h^{-1/2}\varepsilon_N$, $\|\hat{g}(X) - g_0(X)\|_{P,2}\leq \varepsilon_N$, $\|\hat{\mathcal{E}}_{\Delta Y}(X) - \mathcal{E}_{\Delta Y}^0(X)\|_{P,2}\leq \varepsilon_N$; (d) with probability tending to $1$, $\kappa<\hat{g}(X)<1-\kappa$ and $c <\hat{f}_h(d|X)<C$ almost surely $\forall d\in\mathcal{D}$, $\sup_{d\in \mathcal{D}}|\hat{f}_D^{(1)}(d)|<C$, $\sup_{d\in\mathcal{D}}\|\partial_d \hat{f}_h(d|X)\|_{P,\infty}< C$, and $\|\hat{\mathcal{E}}_{\Delta Y}(X)\|_{P,\infty}<C$.
\end{assumption}
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}$.
\begin{theorem}[Uniform Linear Expansion, Panel]\label{thm:unifexp} Suppose assumptions \ref{didas1}, \ref{didas3}, \ref{didas5}, \ref{didas4}, \ref{cdidas1}, \ref{cdidas2}, \ref{as:subexp}, and \ref{cdidas4} hold. Then, for $d\in\mathcal{D}$, if $\varepsilon_N = o(N^{-1/4})$,
\begin{align}
&\widehat{ATT}(d) - \widehat{ATT}(d)^* \notag \\ &= \frac{1}{N}\sum_{i=1}^N \dot{\xi}_i\Bigg[\psi_h^{(1)}(Z_i,\theta_{0h}(d),f_D^0(d),\eta_0(d)) - \frac{\theta_{0h}(d)}{f_D^0(d)}\big(K_h(D_i-d)-E[K_h(D-d)]\big)\Bigg]\notag \\ &+ R^{(1)}(d)
\end{align}
where $\dot{\xi}_i := \xi_i - 1$ and $\sup_{d\in\mathcal{D}} |R^{(1)}(d)| = o_p( (Nh)^{-1/2})$.
\end{theorem}
This theorem is the basis for establishing uniform inference theory using the multiplier bootstrap estimator. We consider the following procedure, see \cite{CCK2014b} and \cite{FHLZ22} for example, to establish valid uniform confidence bands.
\begin{itemize}
\item[1.] Construct $\widehat{ATT}(d)$ and $\hat{\sigma}_N(d)$ on a finite grid of values $d\in \bar{\mathcal{D}} \subset \mathcal{D}$.
\item[2.] For each $b=1, \cdots, B$, draw an i.i.d. sequence of multipliers $\{\xi\}_{i=1}^N$ from a $N(1,1)$ distribution, and construct $\widehat{ATT}(d)_{b}^*$ for all $d\in \bar{\mathcal{D}}$.
\item[3.] Compute $\hat{c}(1-\alpha)$, which we denote as the $(1-\alpha)$-th quantile of
\[
\Bigg\{\max_{d\in \bar{\mathcal{D}}} \frac{\sqrt{N}|\widehat{ATT}(d) - \widehat{ATT}(d)_{b}^*|}{\hat{\sigma}_N(d)} \Bigg \}_{b=1}^B.
\]
\item[4.] For all $d\in\mathcal{D}$, construct the $1-\alpha$ uniform confidence band as
\[
[\widehat{ATT}(d) - \hat{c}(1-\alpha)\hat{\sigma}_N(d)/\sqrt{N}, \quad \widehat{ATT}(d) + \hat{c}(1-\alpha)\hat{\sigma}_N(d)/\sqrt{N}].
\]
\end{itemize}
With Assumption \ref{cdidas4}, we can easily adapt our proof of Theorem \ref{thm:asvar} to establish the uniform consistency of our cross-fitted variance estimator (\ref{asvar}) over $\mathcal{D}$. Then, with Theorem \ref{thm:unifexp}, we can show that the proposed uniform confidence band achieves asymptotic coverage of $1-\alpha$, using results from \cite{CCK2014a} (Proposition 3.2 and Theorem 3.2) and \cite{CCK2014b} (Corollary 3.1). Since this argument is well established in the literature, e.g., see the discussion of Theorem 4.2 in \cite{FHLZ22}, we do not include the formal theoretical discussion here. Instead, we focus on presenting the new results in Theorem \ref{thm:unifexp} and defer its proof to the appendix.
\begin{remark}
\textnormal{A natural extension of our framework is to develop a test for the conditional parallel trends assumption, akin to the approach in \cite{CS2018}, Section 4, which examines differences between the not-yet-treated and the never-treated in the pre-treatment period. Extending such a test to the continuous treatment setting requires a multi-period generalization of the methods considered in this paper, and we suspect that stronger parallel trends assumptions would be necessary for a valid test. While our companion study, \cite{HHZ2024}, proposes estimators that could aid in this analysis, a formal testing procedure remains an open question. Additionally, drawing on insights from \cite{SZ2020}, we recognize that more efficient estimators may exist in the repeated cross-sectional settings than those considered in this paper (Appendix B), and we defer a detailed investigation of such estimators to \cite{HHZ2024}.}
\end{remark}
\section{Simulation}
\setcounter{equation}{0}
\noindent \textbf{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{figs1}, 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{tab01}, 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.
\begin{figure}[htbp!]
\centering
\includegraphics[width=\textwidth]{simulation_histograms.png}
\caption{The simulation results, true $ATT(d) = -0.405$.}
\label{figs1}
\end{figure}
\begin{table}[htbp!]
\caption{Monte Carlo simulation results.}\label{tab01}
\vspace{-10pt}
\begin{center}
\begin{tabular*}{\textwidth}{@{\extracolsep{\fill}}lccccc@{}}
\hline\hline
Setting and Sample Size & Bias & Std & RMSE & AVSE & Coverage \\
\hline
panel, n=2000 & -0.0720 & 0.2725 & 0.2819 & 0.2524 & 0.9180 \\
panel, n=10000 & -0.0198 & 0.1261 & 0.1277 & 0.1262 & 0.9300 \\
cross-sections, n = 2000 & -0.0340 & 0.5745 & 0.5755 & 0.5578 & 0.9440 \\
cross-sections, n = 10000 & 0.0290 & 0.2754 & 0.2769 & 0.2710 & 0.9500 \\
\hline\hline
\end{tabular*}
\end{center}
\end{table}
\section{Empirical Example}
\setcounter{equation}{0}
\subsection{Background}
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 \cite{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. \cite{AF2008} highlights this feature, showing that the PPS reform significantly increased hospitals’ capital-labor ratios and encouraged technology adoption.
Theoretically, \cite{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, \cite{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:
\begin{equation}\label{eq:af}
Y_{i,t} = \alpha_i + \gamma_t + X_{i,t}'\eta + \beta\cdot(D_i\cdot \text{Post}_t) + \varepsilon_{i,t},
\end{equation}
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. \cite{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 \cite{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 \cite{CGS2024}, offering an alternative perspective on the effects of the PPS reform.
\subsection{Setup as a continuous DiD}
Regression \eqref{eq:af} resembles a Two-Way Fixed Effects (TWFE) design but differs in that $D_i$ is continuous. As shown by \cite{CGS2024}, with continuous treatment, the coefficient $\beta$ in \eqref{eq:af} can be viewed as a weighted average of $ATT(d)$ with possible negative weights, which complicates interpretation.\footnote{See Proposition 10 in \cite{CGS2024}. They do not incorporate covariates, but the issue persists.} Our continuous DiD framework addresses this by reframing \cite{AF2008}’s design as follows:
\setcounter{bean}{0}
\begin{list}
{(\alph{beana})}{\usecounter{beana}}
\item No Treatment Pre-PPS: Before the PPS reform, no hospital was treated.
\item Control Group: Hospitals with $D_i = 0$ (no Medicare patients) are the control.
\item Treatment Group: Hospitals with positive Medicare shares (treatment intensities) $D_i > 0$.
\item Outcomes: $Y$ includes the capital-labor ratio or measures of technological adoption.
\item Covariates: $X$ includes number of beds, metro status, private status, number of medical staff, and state dummies.\footnote{We exclude some additional characteristics in \cite{AF2008}—e.g., general, short-term, or federal status—to avoid conditioning on PPS exemption criteria.} For the capital-labor ratio, we also add binary indicators of specialized capital equipments (CT, MRI, etc.).
\item \textit{Conditional} parallel trends:
\[
E[Y_{t}(0) - Y_{t-1}(0)|X,D=d] = E[Y_{t}(0) - Y_{t-1}(0)|X,D=0],
\]
i.e., absent the PPS reform, hospitals with share $D=d$ would have experienced similar changes over time as hospitals with no Medicare inpatients (shares $D = 0$), conditional on hospital-specific covariates $X$ determined before the PPS reform.
\end{list}
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{eq:af}), 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 \cite{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 \cite{AF2008}.
\subsection{Results}
First, we examine the results for capital-labor ratio. All estimated ATTs are positive, mirroring the findings in \cite{AF2008} and suggesting that the PPS reform led to an increase in the capital-labor ratio. For comparison, \cite{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 \cite{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.
\begin{figure}[hbtp!]
\centering
\includegraphics[width=\textwidth, height=0.65\textheight]{dmlcdid_all.png}
\vspace{-10pt}
\caption{$\widehat{ATT}(d)$, panel data}
\label{fig:cdid_avg_panel}
\end{figure}
\begin{remark}
\textnormal{
We apply 5-fold cross-fitting, shuffling the data before sample splitting to prevent over-representation in subsamples. A second-order Gaussian kernel with an undersmoothing bandwidth $h = 1.06\times\hat{\sigma}_{\tilde{D}}N^{-1/4}$ is used to estimate both the density $f_D(d)$ and the conditional mean $E[K_h(D-d)|X]$ (see \cite{S2018}). The infinite-dimensional nuisance parameters are estimated using the Random Forest (RF) from the Python scikit-learn package, with 200 trees of maximum depth 20 and fixed minimum leaf size 5. The RF is chosen for its flexibility to accommodate both continuous and discrete covariates, though other ML methods, such as deep neural networks, can be similarly employed. The standard errors are obtained from the cross-fitted estimator defined in (\ref{asvar}) and used to construct the 95-percent pointwise confidence intervals and the bootstrap uniform confidence bands. For the bootstrap CIs, we use Gaussian multipliers $\{\xi_i\}_{i=1}^N$ drawn from a normal distribution with $E[\xi_i]= Var[\xi_i] = 1$, with $B=1000$ repetitions.
}
\end{remark}
\section{Conclusion}
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 \cite{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.
\bibliographystyle{chicago}
\begin{thebibliography}{}
\bibitem[\protect\citeauthoryear{Abadie}{Abadie}{2005}]{Abadie2005}
Abadie, A. (2005).
\newblock Semiparametric difference-in-differences estimators. \newblock {\em Review of Economic Studies\/} {\em 72}, 1--19.
\bibitem[\protect\citeauthoryear{Acemoglu and Finkelstein}{Acemoglu and Finkelstein}{2008}]{AF2008} Acemoglu, D. and Finkelstein, A. (2008). \newblock Input and technology choices in regulated industries: evidence from the health care sector.
\newblock {\em Journal of Political Economy\/} {\em 116}, 837--880.
\bibitem[\protect\citeauthoryear{Ananat et al.}{Ananat et al.}{2022}]{AGHP22} Ananat, E., Glasner, B., Hamilton, C., and Parolin, Z. (2022).
\newblock Effects of the expanded child tax credit on employment outcomes: evidence from real-world data from April to December 2021.
\newblock Technical Working Paper 29823, National Bureau of Economic Research.
\bibitem[\protect\citeauthoryear{Athey and Imbens}{Athey and Imbens}{2006}]{AI2006} Athey, S. and Imbens, G. W. (2006).
\newblock Identification and inference in nonlinear difference‐in‐differences models.
\newblock {\em Econometrica\/} {\em 74}, 431--497.
\bibitem[\protect\citeauthoryear{Belloni et al.}{Belloni et al.}{2017}]{BCFH17} Belloni, A., Chernozhukov, V., Fernández‐Val, I., and Hansen, C. (2017).
\newblock Program evaluation and causal inference with high‐dimensional data.
\newblock {\em Econometrica\/} {\em 85}, 233--298.
\bibitem[\protect\citeauthoryear{Bibaut and van der Laan}{Bibaut and van der Laan}{2017}]{BL2017} Bibaut, A. F. and van der Laan, M. J. (2017).
\newblock Data-adaptive smoothing for optimal-rate estimation of possibly non-regular parameters, arXiv:1706.07408.
\bibitem[\protect\citeauthoryear{Callaway and Sant'Anna}{Callaway and Sant'Anna}{2018}]{CS2018} Callaway, B. and Sant'Anna, P. H. (2018).
\newblock Difference-in-differences with multiple time periods and an application on the minimum wage and employment, arXiv:1803.09015v2.
\bibitem[\protect\citeauthoryear{Callaway and Sant'Anna}{Callaway and Sant'Anna}{2021}]{CS2021}
Callaway, B. and Sant'Anna, P. H. (2021).
\newblock Difference-in-differences with multiple time periods. \newblock {\em Journal of Econometrics\/} {\em 225}, 200--230.
\bibitem[\protect\citeauthoryear{Callaway}{Callaway}{2023}]{CALL2023}
Callaway, B. (2023).
\newblock Difference-in-differences for policy evaluation.
\newblock {\em Handbook of Labor, Human Resources and Population Economics\/}, 1--61.
\bibitem[\protect\citeauthoryear{Callaway, Goodman-Bacon, and Sant'Anna}{Callaway et al.}{2024}]{CGS2024}
Callaway, B., Goodman-Bacon, A., and Sant'Anna, P. H. (2024).
\newblock Difference-in-differences with a continuous treatment.
\newblock Technical Working Paper 32117, National Bureau of Economic Research.
\bibitem[\protect\citeauthoryear{Cattaneo and Jansson}{Cattaneo and Jansson}{2021}]{CJ21} Cattaneo, M. D. and Jansson, M. (2021).
\newblock Average density estimators: efficiency and bootstrap consistency.
\newblock {\em Econometric Theory\/} {\em 38}, 1140--1174.
\bibitem[\protect\citeauthoryear{Chang}{Chang}{2020}]{Chang2020}
Chang, N. C. (2020).
\newblock Double/debiased machine learning for difference-in-differences models.
\newblock {\em Econometrics Journal\/} {\em 23}, 177--191.
\bibitem[\protect\citeauthoryear{Chernozhukov et al.}{Chernozhukov et al.}{2014a}]{CCK2014a}
Chernozhukov, V., Chetverikov, D., and Kato, K. (2014).
\newblock Gaussian approximation of suprema of empirical processes.
\newblock {\em Annals of Statistics\/} {\em 42}, 1564--1597.
\bibitem[\protect\citeauthoryear{Chernozhukov et al.}{Chernozhukov et al.}{2014b}]{CCK2014b}
Chernozhukov, V., Chetverikov, D., and Kato, K. (2014).
\newblock Anti-concentration and honest, adaptive confidence bands.
\newblock {\em Annals of Statistics\/} {\em 42}, 1787--1818.
\bibitem[\protect\citeauthoryear{Chernozhukov et al.}{Chernozhukov et al.}{2018}]{CCDDHNR} Chernozhukov, V., Chetverikov, D., Demirer, M., Duflo, E., Hansen, C., Newey, W., and Robins, J. (2018).
\newblock Double/debiased machine learning for treatment and structural parameters.
\newblock {\em Econometrics Journal\/} {\em 21}, C1--C68.
\bibitem[\protect\citeauthoryear{Colangelo and Lee}{Colangelo and Lee}{2025}]{CYYL25} Colangelo, K. and Lee, Y. Y. (2025).
\newblock Double debiased machine learning non-parametric inference with continuous treatments.
\newblock {\em Journal of Business \& Economic Statistics\/}, 1--26.
\bibitem[\protect\citeauthoryear{Cook et al.}{Cook et al.}{2023}]{CJLR23} Cook, L. D., Jones, M. E., Logan, T. D., and Ros\'{e}, D. (2023).
\newblock The evolution of access to public accommodations in the United States.
\newblock {\em The Quarterly Journal of Economics\/} {\em 138}, 37--102.
\bibitem[\protect\citeauthoryear{de Chaisemartin et al.}{de Chaisemartin et al.}{2022}]{CDPV2022} de Chaisemartin, C., D'Haultfoeuille, X., Pasquier, F., and Vazquez-Bare, G. (2022). \newblock Difference-in-differences estimators for treatments continuously distributed at every period, arXiv:2201.06898.
\bibitem[\protect\citeauthoryear{de Chaisemartin and D'Haultfoeuille}{de Chaisemartin and D'Haultfoeuille}{2023}]{DD2023}
de Chaisemartin, C. and D'Haultfoeuille, X. (2023).
\newblock Two-way fixed effects and differences-in-differences with heterogeneous treatment effects: a survey.
\newblock {\em Econometrics Journal\/} {\em 26}, C1--C30.
\bibitem[\protect\citeauthoryear{D'Haultfoeuille et al.}{D'Haultfoeuille et al.}{2023}]{DHS2021} D'Haultfoeuille, X., Hoderlein, S., and Sasaki, Y. (2023).
\newblock non-parametric difference-in-differences in repeated cross-sections with continuous treatments.
\newblock {\em Journal of Econometrics\/} {\em 234}, 664--690.
\bibitem[\protect\citeauthoryear{Fan et al.}{Fan et al.}{1996}]{FYT96}
Fan, J., Yao, Q., and Tong, H. (1996).
\newblock Estimation of conditional densities and sensitivity measures in nonlinear dynamical systems.
\newblock {\em Biometrika\/} {\em 83}, 189--206.
\bibitem[\protect\citeauthoryear{Fan and Yao}{Fan and Yao}{2003}]{FY2003} Fan, J. and Yao, Q. (2003).
\newblock {\em Nonlinear Time Series: non-parametric and Parametric Methods\/} (Vol. 20).
\newblock New York: Springer.
\bibitem[\protect\citeauthoryear{Fan et al.}{Fan et al.}{2022}]{FHLZ22}
Fan, Q., Hsu, Y. C., Lieli, R. P., and Zhang, Y. (2022).
\newblock Estimation of conditional average treatment effects with high-dimensional data.
\newblock {\em Journal of Business and Economic Statistics\/} {\em 40}, 313--327.
\bibitem[\protect\citeauthoryear{Galvao and Wang}{Galvao and Wang}{2015}]{GW2015}
Galvao, A. F. and Wang, L. (2015).
\newblock Uniformly semiparametric efficient estimation of treatment effects with a continuous treatment.
\newblock {\em Journal of the American Statistical Association\/} {\em 110}, 1528--1542.
\bibitem[\protect\citeauthoryear{Haddad et al.}{Haddad et al.}{2024}]{HHZ2024} Haddad, M. F., Huber, M., and Zhang, L. Z. (2024).
\newblock Difference-in-differences with time-varying continuous treatments using double/debiased machine learning, arXiv:2410.21105.
\bibitem[\protect\citeauthoryear{H\"{a}rdle}{H\"{a}rdle}{1990}]{Hardle90} H\"{a}rdle, W. (1990).
\newblock {\em Applied non-parametric regression\/} (No.19).
\newblock United Kingdom: Cambridge University Press
\bibitem[\protect\citeauthoryear{Hettinger et al.}{Hettinger et al.}{2025}]{HLM25}
Hettinger, G., Lee, Y., and Mitra, N. (2025).
\newblock Multiply robust difference-in-differences estimation of causal effect curves for continuous exposures.
\newblock {\em Biometrics\/} {\em 81}, ujaf015.
\bibitem[\protect\citeauthoryear{Heckman et al.}{Heckman et al.}{1997}]{HIT97}
Heckman, J. J., Ichimura, H., and Todd, P. E. (1997).
\newblock Matching as an econometric evaluation estimator: Evidence from evaluating a job training programme.
\newblock {\em Review of Economic Studies\/} {\em 64}, 605--654.
\bibitem[\protect\citeauthoryear{Heckman et al.}{Heckman et al.}{1998}]{HIST98}
Heckman, J., Ichimura, H., Smith, J., and Todd, P. (1998).
\newblock Characterizing selection bias using experimental data.
\newblock {\em Econometrica\/} {\em 66}, 1017--1098.
\bibitem[\protect\citeauthoryear{Hines et al.}{Hines et al.}{2022}]{HDDV22}
Hines, O., Dukes, O., Diaz-Ordaz, K., and Vansteelandt, S. (2022).
\newblock Demystifying statistical learning based on efficient influence functions.
\newblock {\em The American Statistician\/} {\em 76}, 292--304.
\bibitem[\protect\citeauthoryear{Hirano and Imbens}{Hirano and Imbens}{2004}]{HI2004} Hirano, K. and Imbens, G. W. (2004).
\newblock The propensity score with continuous treatments.
\newblock {\em Applied Bayesian Modeling and Causal Inference from Incomplete-Data Perspectives\/} {\em 226164}, 73--84.
\bibitem[\protect\citeauthoryear{Hong}{Hong}{2013}]{HONG2013} Hong, S. H. (2013).
\newblock Measuring the effect of napster on recorded music sales: difference‐in‐differences estimates under compositional changes.
\newblock {\em Journal of Applied Econometrics\/} {\em 28}, 297--324.
\bibitem[\protect\citeauthoryear{Kallus and Zhou}{Kallus and Zhou}{2018}]{KZ2018}
Kallus, N. and Zhou, A. (2018).
\newblock Confounding-robust policy improvement.
\newblock {\em Advances in Neural Information Processing Systems\/} {\em 31}.
\bibitem[\protect\citeauthoryear{Kennedy et al.}{Kennedy et al.}{2017}]{KMMS17} Kennedy, E. H., Ma, Z., McHugh, M. D., and Small, D. S. (2017).
\newblock Non‐parametric methods for doubly robust estimation of continuous treatment effects.
\newblock {\em Journal of the Royal Statistical Society\/}, Series B (Statistical Methodology) {\em 79}, 1229--1245.
\bibitem[\protect\citeauthoryear{Kennedy}{Kennedy}{2024}]{KENNEDY24} Kennedy, E. H. (2024).
\newblock Semiparametric doubly robust targeted double machine learning: a review.
\newblock {\em Handbook of Statistical Methods for Precision Medicine\/}, 207--236.
\bibitem[\protect\citeauthoryear{Li and Racine}{Li and Racine}{2007}]{LR07}
Li, Q. and Racine, J.S. (2007).
\newblock {\em non-parametric Econometrics: Theory and Practice\/}.
\newblock Princeton, NJ: Princeton University Press.
\bibitem[\protect\citeauthoryear{Roth et al.}{Roth et al.}{2023}]{RSBP23}
Roth, J., Sant’Anna, P. H., Bilinski, A., and Poe, J. (2023).
\newblock What’s trending in difference-in-differences? A synthesis of the recent econometrics literature.
\newblock {\em Journal of Econometrics\/} {\em 235}, 2218--2244.
\bibitem[\protect\citeauthoryear{Rubin}{Rubin}{1974}]{Rubin74} Rubin, D. B. (1974).
\newblock Estimating causal effects of treatments in randomized and nonrandomized studies.
\newblock {\em Journal of Educational Psychology\/} {\em 66}, 688--701.
\bibitem[\protect\citeauthoryear{Sant'Anna and Xu}{Sant'Anna and Zhao}{2025}]{SX2025} Sant'Anna, P. H. and Xu, Q. (2025).
\newblock Difference-in-Differences with compositional changes, arXiv:2304.13925v2.
\bibitem[\protect\citeauthoryear{Sant'Anna and Zhao}{Sant'Anna and Zhao}{2020}]{SZ2020} Sant'Anna, P. H. and Zhao, J. (2020).
\newblock Doubly robust difference-in-differences estimators. \newblock {\em Journal of Econometrics\/} {\em 219}, 101--122.
\bibitem[\protect\citeauthoryear{Semenova and Chernozhukov}{Semenova and Chernozhukov}{2021}]{SC2021} Semenova, V. and Chernozhukov, V. (2021).
\newblock Debiased machine learning of conditional average treatment effects and other causal functions.
\newblock {\em Econometrics Journal\/} {\em 24}, 264--289.
\bibitem[\protect\citeauthoryear{Silverman}{Silverman}{2018}]{S2018} Silverman, B. W. (2018).
\newblock {\em Density Estimation for Statistics and Data Analysis\/}.
\newblock Routledge.
\bibitem[\protect\citeauthoryear{Su et al.}{Su et al.}{2019}]{SUZ19}
Su, L., Ura, T., and Zhang, Y. (2019).
\newblock Non-separable models with high-dimensional data. \newblock {\em Journal of Econometrics\/} {\em 212}, 646--677.
\bibitem[\protect\citeauthoryear{van der Vaart and Wellner}{van der Vaart and Wellner}{1996}]{VW96}
van der Vaart, A.W. and Wellner, J.A. (1996).
\newblock {\em Weak Convergence and Empirical Processes: With Applications to Statistics\/}.
\newblock New York: Springer.
\bibitem[\protect\citeauthoryear{Zeng et al.}{Zeng et al.}{2022}]{ZDS2022} Zeng, H. S., Danaher, B., and Smith, M. D. (2022).
\newblock Internet governance through site shutdowns: the impact of shutting down two major commercial sex advertising sites. \newblock {\em Management Science\/} {\em 68}, 8234--8248.
\bibitem[\protect\citeauthoryear{Zimmert}{Zimmert}{2020}]{ZIMMERT2020} Zimmert, M. (2020).
\newblock Efficient difference-in-differences estimation with high-dimensional common trend confounding, arXiv:1809.01643v5.
\end{thebibliography}
\newpage