Nonlinear and Nonseparable Structural Functions in Fuzzy Regression Discontinuity Designs

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.

97,338 characters · 15 sections · 62 citation commands

Rendered from LaTeX for readability, not typeset faithfully. Citation keys are highlighted; maths is left as source; figures, tables and equation environments are summarised rather than reproduced; unrecognised commands are greyed out so nothing is silently dropped. Email addresses are removed.

abstractMany empirical examples of regression discontinuity (RD) designs concern a continuous treatment variable, but the theoretical aspects of such models are less studied. This study examines the identification and estimation of the structural function in fuzzy RD designs with a continuous treatment variable. The structural function fully describes the causal impact of the treatment on the outcome. We show that the nonlinear and nonseparable structural function can be nonparametrically identified at the RD cutoff under shape restrictions, including monotonicity and smoothness conditions. Based on the nonparametric identification equation, we propose a three-step semiparametric estimation procedure and establish the asymptotic normality of the estimator. The semiparametric estimator achieves the same convergence rate as in the case of a binary treatment variable. As an application of the method, we estimate the causal effect of sleep time on health status by using the discontinuity in natural light timing at time zone boundaries. {\bf Keywords:} Causal Inference, Continuous Treatment, Dual Monotonicity, Nonparametric Identification, Semiparametric Estimation, Asymptotic Normality.

Introduction

The regression discontinuity (RD) design is one of the most credible approaches to causal inference in non-experimental settings. In an RD design, the researcher is interested in the effect of a treatment $T$ on some outcome $Y$. The basic idea is that there is an observed running variable $R$ (also called score or index or forcing variable) such that the treatment varies discontinuously when the running variable crosses some cutoff (also called threshold) value $\bar{r}$. By utilizing this discontinuity, the researcher has the power to identify and estimate the causal impact of interest.

Most theoretical studies of the RD design assume that the treatment is a binary intervention. However, in empirical settings, researchers may be interested in a continuous treatment that takes value inside an interval. Such examples include sleep time, air pollution level, and medical spending. The goal of this study is to provide methods for examining the causal effect of a continuous treatment variable in an RD setting.

It takes a few steps to extend the idea of RD design from a binary treatment to a continuous one. With a binary treatment, the sharp design refers to the case where the running variable completely determines the treatment. In particular, the treatment changes from $0$ to $1$ when the running variable crosses the cutoff. The fuzzy design refers to the case where the treatment probability jumps at the cutoff. The jump can be smaller and need not be from $0$ to $1$. The sharp and fuzzy designs of a binary treatment are demonstrated in Figure (ref).

figure[figure omitted — 433 chars of source]

When the treatment variable is continuous, the representation of the RD becomes more complicated than the binary case. The reason is that the distribution of a binary variable can be completely summarized by the scalar treatment probability as in Figure (ref)(b), while a continuous variable contains much more information. Specifically, we can consider quantile regressions of the treatment on the running variable at different quantile levels. Each quantile level would deliver a different regression model with a different discontinuity. Eventually, we would obtain an infinite number of regression discontinuities based on all quantile levels of the treatment. This infinite set of regression discontinuities can be represented as the entire variation between the conditional quantile function of the treatment from just below and just above the cutoff. Figure (ref) provides a demonstration.

figure[figure omitted — 1,009 chars of source]

The exogenous variation of the treatment contained in the aforementioned set of regression discontinuities provides tremendous identification power on the causal effect of interest. To fully express the causal effect of the treatment $T$ on the outcome $Y$, we introduce the structural function

align*[align* omitted — 43 chars of source]

where $\varepsilon$ contains unobserved causal factors (for easy reference, $\varepsilon$ will be called the error term hereafter). The structural function $g^*$ specifies how the treatment $T$ determines the outcome $Y$ together with the running variable $R$ and error term $\varepsilon$.

Consider an empirical example for concreteness, where we are interested in the causal impact of sleep time on health. Figure (ref)(a) shows the histogram of sleep time based on the American Time Use Survey (ATUS) and demonstrates that sleep time is indeed a continuous treatment variable. The causal identification is based on exploiting the discontinuity in the timing of natural light at time zone boundaries. Individuals living on the late sunset side of the time zone boundary tend to go to bed at a later time, while in the morning, everyone gets up and goes to work at 8 am. This generates an exogenous variation in the sleep time across the time zone boundary.\footnote{This identification strategy is first proposed by Giuntella2019sunset within a linear model. See the empirical application in Section (ref) for more details.} Similar to the demonstration in Figure (ref)(b), we would expect the distribution of the sleep time for individuals living on the early sunset side to first-order stochastically dominate the distribution on the late sunset side. This relationship is supported by Figure (ref)(b) based on nonparametric estimates of the conditional quantiles of sleep time. In this example, the running variable is the distance to the time zone boundary, and the cutoff $\bar{r}$ is at the time zone boundary. The error term $\varepsilon$ may contain unobserved eating habits that correlate with both health and sleep time.

figure[figure omitted — 782 chars of source]

The goal of the RD design is to use the discontinuity to identify the structural function $g^*$ at the cutoff $\bar{r}$. When the treatment is binary, the information contained in the structural function can be reduced to a scalar treatment effect

align*[align* omitted — 74 chars of source]

which is the difference in outcome when the treatment is manipulated from $0$ to $1$. However, when the treatment is continuous, the structural function is an infinite-dimensional object and is much harder to identify.

In practice, empirical studies often use the two stage least squares (TSLS) method to estimate the following Wald ratio:

align*[align* omitted — 232 chars of source]

There are two motivations behind this procedure. First, in the binary treatment case, the Wald ratio would identify the treatment effect.\footnote{As shown in hahn2001identification, the Wald ratio identifies the average treatment effect in the sharp design and the local average treatment effect (for the compliers) in the fuzzy design.} Second, in the continuous treatment case, if the structural function is linear and separable in the treatment, that is, if the structural function can be decomposed as

align*[align* omitted — 76 chars of source]

then the Wald ratio would identify the slope coefficient $\beta$ of the treatment.\footnote{Such a linear specification of RD design with a continuous treatment can be found in Section 3.4.2 of lee2010regression.} However, the Wald ratio cannot identify the structural function in general because the structural function is infinite-dimensional while the Wald ratio is one-dimensional. Any attempt to condense the structural function into a scalar bears the risk of dampening the causal interpretation of the model.

The preceding discussion shows that the general identification of the structural function in RD designs remains an unsolved issue. It is desirable to know whether the structural function (at the cutoff) can be identified without the aforementioned linearity and separability conditions. This issue is a practical concern. For instance, in the time zone example, there are reasons for one to believe that the structural function is nonlinear and nonseparable in the treatment.\footnote{The nonlinearity can be due to the fact that both undersleeping and oversleeping are harmful to health. The nonseparability can be due to the effect heterogeneity caused by unobserved eating habits, which affect both sleep time and health.} The optimal sleep time can only be determined after the identification of the nonlinear structural function. From the theoretical perspective, it is wise to achieve identification in the nonparametric sense and avoid functional form restrictions such as linearity and separability that do not have economic theory foundations. As an advantage, the more general specification allows the treatment effect to be heterogeneous across different levels of the treatment and outcome.

The current study aims to precisely tackle the identification and estimation of the possibly nonlinear and nonseparable structural function. The nonparametric identification result is established based on shape restrictions, including monotonicity and smoothness conditions. The idea behind the identification result is that we are using the infinite set of regression discontinuities in Figure (ref)(b) to identify the infinite-dimensional structural function. The monotonicity condition restricts the structural function $g^*$ to be strictly increasing in the error term $\varepsilon$. This condition requires the error term to be one-dimensional, which is the potential restriction of the model. However, this condition is common in the nonparametric identification literature matzkin2003nonparametric and is satisfied by most, if not all, parametric models used in practice. The smoothness and other regularity conditions imposed in this paper are common in the RD literature.

A semiparametric estimation procedure is developed based on the nonparametric identification result. The structural function is parametrized while nonlinearity and nonseparability are maintained. One such parametrization could be

align*[align* omitted — 116 chars of source]

The relationship between the treatment and the running variable is left to be nonparametric. Under appropriate conditions, the semiparametric estimator of the structural parameter $\gamma=(\gamma_1,\gamma_2,\gamma_3)$ is shown to be consistent and asymptotically normal. As an interesting finding, the convergence rate of the semiparametric estimator, $n^{-2/5}$, is the same as in the binary treatment case. There is no loss in terms of convergence rate when extending the RD design from the binary treatment case to the continuous case. The faster convergence rate is due to the integral smoothing in the estimation of the criterion function constructed from the identification equation. To understand this phenomenon, one can consider the analogy in regular semiparametric estimation theory, where the first step is nonparametric while the second step recovers the parametric rate.

The rest of the paper is organized as follows. The remaining part of this section discusses the literature. Section (ref) introduces the RD model with a continuous treatment and presents the nonparametric identification result. Section (ref) proposes the semiparametric estimation procedure and derives its asymptotic properties. Section (ref) presents the empirical application and simulation studies. The technical proofs for the identification and estimation results are collected in Appendices (ref) and (ref), respectively.

commentThe regression discontinuity (RD) design is one of the most credible non-experimental methods for causal inference and program evaluation. In an RD design, there is an outcome of interest $Y$, a treatment variable $T$, and a running variable $R$ (also called score or index or forcing variable). Their causal relationship is modeled by a structural function \begin{align*} Y = g^*(T,R,\varepsilon), \end{align*} where $\varepsilon$ represents unobserved causal factors. The unique feature of an RD design is that there is a discontinuity in the treatment choice model. Specifically, there is a cutoff value $\bar{r}$ (also called threshold) of the running variable $R$ such that the treatment $T$ varies discontinuously when the running variable crosses that cutoff. In such a case, the RD design can be interpreted under the local average treatment effect (LATE) framework imbens1994identification. In practice, however, the treatment under consideration may well be continuous. Examples of continuous treatments in RD designs include air pollution measures, medical expenditure, and government spending. Table (ref) summarizes some of the empirical studies that involve a continuous treatment. The standard RD method, developed for binary treatments, cannot be readily applied when the treatment is continuous. The reason is that the main object of interest, the conditional average structural function (CASF): $t \mapsto \mathbb{E}[g^*(t,\bar{r},\varepsilon) | R=\bar{r}]$, is much more complicated when the treatment variable is continuous. In the binary case when $t \in \{0,1\}$, the CASF can be summarized by a scalar treatment effect that can be identified (under standard conditions) as a Wald ratio: \begin{align*} \mathbb{E}[g^*(1,\bar{r},\varepsilon) | R=\bar{r}] - \mathbb{E}[g^*(0,\bar{r},\varepsilon) | R=\bar{r}] = \frac{ \lim_{r \uparrow \bar{r}} \mathbb{E}[Y|R=r] - \lim_{r \downarrow \bar{r}} \mathbb{E}[Y|R=r] }{\lim_{r \uparrow \bar{r}} \mathbb{E}[T|R=r] - \lim_{r \downarrow \bar{r}} \mathbb{E}[T|R=r]}, \end{align*} which can be estimated by the two-stage least squares (TSLS) method in the neighborhood of the cutoff. However, in the case where the treatment is continuous, it is difficult to condense the CASF into a scalar, and any such attempt bears the risk of dampening the causal interpretation of the model. The usual solution, common in empirical studies, is to assume that the CASF is linear in the treatment and hence can be characterized by the slope coefficient, which is equal to the Wald ratio specified above.\footnote{Such a linear specification can be found, for example, in Section 3.4.2 of lee2010regression.} The drawback of this linear model is that it does not allow for heterogeneity in the treatment effect, which is the key of modern causal inference methods and the potential outcome framework. The current study aims to directly identify and estimate the possibly nonlinear and nonseparable structural function. This more general specification of the structural function allows the treatment effect to be heterogeneous across different levels of $T$ and $\varepsilon$. The main contribution of this paper includes the nonparametric identification of the CASF local to the RD cutoff and a semiparametric estimation procedure. The estimator converges at the $n^{-2/5}$-rate, which is the same rate for the estimator in the binary treatment case.\footnote{This rate is proved, for example, in the Appendix of hahn2001identification.} The asymptotic normal distribution of the estimator is also derived. There is a price to pay going from parametric identification to nonparametric identification. Besides smoothness and other regularity conditions, we need to assume a dual monotonicity condition for the model. The dual monotonicity condition (explained in detail later) restricts both the structural function and the treatment choice to be monotonic. Such a condition allows the structural function to be nonlinear with respect to the treatment $T$ and nonseparable between the treatment $T$ and the unobserved $\varepsilon$. Hence, it is satisfied by virtually all parametric models used in practice.

Relation to the literature

The RD method is first introduced by thistlethwaite1960regression into the literature. hahn2001identification establish the theoretical foundation of RD designs by using the potential outcome framework and show that the RD Wald ratio can be interpreted as the local average treatment effect (LATE) for compliers local to the cutoff. Early reviews of the RD design can be found in IMBENS2008regression and lee2010regression. For more recent reviews, see cattaneo2017regression and cattaneo2021regression.

There are many empirical papers that study the causal effect of a continuous treatment in an RD design, some of which are given in Table (ref). As explained earlier, these studies apply the TSLS method to estimate the Wald ratio. Hence, there is room for potential improvement in these settings by using the semiparametric estimator developed in the current study.

table[table omitted — 1,443 chars of source]

The theoretical literature on RD designs focuses on the case of a binary treatment variable. The one exception is the recent paper by dong2021regression, which studies RD designs specifically with a continuous treatment variable. Under simple conditions, they propose a way to identify and estimate the Quantile specific LATE. This parameter bears a causal interpretation as it is a weighted average of the derivative of the structural function dong2021regression. It can also be understood as the treatment effect given a particular quantile of the treatment. Their results are established under conditions weaker than the ones in our paper. In particular, they do not assume the monotonicity condition of the structural function. In certain situations, however, the policy design process may require information beyond the weighted average of the structural function. The current paper takes a different approach and aims to identify the structural function directly.

It has become common in the literature to identify a certain weighted average of the structural function as the causal estimand. One of the first examples is the 2SLS estimation with a multivalued treatment angrist-imbens1995. The reason for this trend is twofold: the direct identification of the structural function is difficult, and the researchers want the assumptions they make to be minimal. However, the weighted average only provides summary information on the structural function, which is not sufficient in optimal policy designs. In this paper, we make an effort to identify the structural function itself at the expense of making stronger assumptions. In the empirical application in Section (ref), we show that the estimated nonlinear structural function can help determine the optimal sleep time while the TSLS estimates cannot.

The identification in the RD design is related to that in the instrumental variables (IV) models of triangular systems. The control function approach described in imbens2009identification states that the variation in the treatment becomes exogenous after conditioning on the control function. As explained later, a similar phenomenon is also observed in the RD model with a continuous treatment. It explains the intuition behind the nonparametric identification equation. Another relevant literature is the one that studies instruments with small support torgovitsky2015identification,D2015id,TORGOVITSKY2017minimum. These papers examine a model with a discrete instrumental variable and a continuous treatment variable. Since RD can be interpreted as a local IV approach, our framework is related to the large body of this IV literature.

That said, this paper is not a straightforward extension of the results from the IV literature. The difference between the RD design and the IV approach includes the following. First, the identification in the IV model relies on the (conditional) independence of the IV with the error term, while the identification in RD designs is based on the discontinuity and does not depend on any independence assumption. This is one of the reasons that the RD method is considered to be more credible than IV for causal inference. Second, in an RD design, the running variable directly affects both the outcome and the treatment and hence does not satisfy the exclusion restriction typically required in IV models. This inclusion also gives rise to the unique issue of extrapolation away from the cutoff. Third, the estimation procedure in the RD design focuses on the local neighborhood of the cutoff. It is theoretically more challenging to derive the asymptotic properties of the estimator.

The problem studied by this paper is also related to the broad literature on the nonparametric identification of structural functions. Relevant papers include matzkin2003nonparametric and hoderlein2007identification,hoderlein2009identification. The identification there relies on the exogeneity of the treatment, which is not required in RD designs.

RD Design with a continuous treatment

This section describes the RD model with a continuous treatment, explains the assumptions of the model, and discusses the nonparametric identification of the structural function local to the cutoff.

The model

We study the following causal equation:

align[align omitted — 73 chars of source]

where $Y$ is the outcome of interest, $T$ is the treatment, and $R$ is the running variable. The scalar variable $\varepsilon$ represents unobserved causal factors in the outcome equation. We assume all the random variables are absolutely continuous. The function $g^*$ is the unknown true structural function.

The running variable $R$ partly determines the treatment $T$ by the following treatment choice function:

align[align omitted — 150 chars of source]

where $\bar{r}$ is the cutoff value, and $U_0$ and $U_1$ are scalar variables, representing other factors that are not observable to an econometrician. For easy reference, they will be referred to as the error terms hereafter. The important feature of the RD design is that the treatment varies discontinuously when the running variable crosses the cutoff $\bar{r}$. The functions $m_0$ and $m_1$ represent respectively the treatment choice mechanism when $R$ is below and above the cutoff.

It is important to point out that the variables $U_0,U_1,T$ and $R$ are allowed to be correlated with the error term $\varepsilon$. If we assume $\varepsilon$ to be independent of $(T,R)$, then we can follow matzkin2003nonparametric or hoderlein2007identification to identify the structural function. If we assume $\varepsilon \perp R$ and $R$ is excluded from $m_0,m_1$ and $g$, then we can follow torgovitsky2015identification to identify the structural function by treating the binary variable $\mathbf{1}\{R \geq \bar{r}\}$ as the instrument.

We make the following assumptions on the model imposed by ((ref)) - ((ref)). Let $\mathcal{G}$ be the set of candidate structural functions such that the true $g^*$ is contained in $\mathcal{G}$. That is, $\mathcal{G}$ is the infinite-dimensional parameter space where the structural function belongs to. Denote the conditional distribution function by $F_{\cdot|\cdot}(\cdot | \cdot)$, the conditional density function by $f_{\cdot|\cdot}(\cdot | \cdot)$, and the conditional quantile function by $F^{-1}_{\cdot|\cdot}(\cdot | \cdot)$.

assumption[Dual Monotonicity] \ \begin{enumerate}[label = (\roman*)] • Every $g \in \mathcal{G}$ satisfies that for each given $T=t$ and $R=r$, $g$ is strictly increasing in $\varepsilon$. • For each given $R=r$, $m_0$ is strictly increasing in $U_0$ and $m_1$ is strictly increasing in $U_1$. \end{enumerate}

assumption[Smoothness] \ \begin{enumerate} [label = (\roman*)] • The functions $m_0,m_1$ and every $g \in \mathcal{G}$ are continuous on their respective domains. • The conditional quantile functions $F^{-1}_{U_0 | R}(u | r)$ and $F^{-1}_{U_1 | R}(u | r)$ are strictly increasing in $u$ and continuous in $(u,r)$. • The conditional distribution functions $F_{\varepsilon | U_0,R}(e | u , r)$ and $F_{\varepsilon | U_1,R}(e | u , r)$ are strictly increasing in $e$ and are continuous in $r$ at $\bar{r}$. • The running variable $R$ is absolutely continuous, and its density is strictly positive around the cutoff $\bar{r}$. \end{enumerate}

assumption[Rank Similarity] $U_0 | (\varepsilon,R=\bar{r}^-) $ has the same distribution as $U_1 | (\varepsilon,R=\bar{r}^+) $. That is, \begin{align*} \lim_{r \uparrow \bar{r}} f_{U_0|\varepsilon,R}(u | e,r) = \lim_{r \downarrow \bar{r}} f_{U_1|\varepsilon,R}(u | e,r). \end{align*}

Assumption (ref) defines a one-to-one mapping between $(Y,T)$ and $(\varepsilon,U_0,U_1)$ for a given value of $R$. Assumption (ref) states that except for the discontinuity introduced in ((ref)), everything else is assumed to be reasonably smooth. Assumption (ref) is similar to Assumption 3 in dong2021regression. It imposes the rank similarity condition chernozhukov2005endogeneity on $(U_0,U_1)$.

The treatment choice functions $(m_0,m_1)$ are not identified. Rather than trying to identify them, it is more convenient to consider a normalization to a quantile representation. By using the monotonicity of $m_0$ and $m_1$ in Assumption (ref)(ii), we define

align[align omitted — 155 chars of source]

as the conditional rank of $T$ given $R$.\footnote{The second equality in Equation ((ref)) is proved in Lemma (ref)} Then the treatment choice model in ((ref)) can be written as

align*[align* omitted — 132 chars of source]

where

align*[align* omitted — 128 chars of source]

By using $[r_0,r_1]$ to denote the support of $R$, we can write the domains of $h_0$ and $h_1$ respectively as $[r_0,\bar{r}] \times [0,1]$ and $[\bar{r},r_1] \times [0,1]$.

The following lemma shows that the function $h$ defined above is the conditional quantile function of $T$ given $R$, and the quantile representation is a valid normalization in the sense that it preserves the monotonicity and smoothness conditions. Consequently, the function $h$ (including both $h_0$ and $h_1$) and the rank $U = h^{-1}(R,T)$ are identified from the data, where $h^{-1}$ denotes the inverse of $h$ with respect to the second argument $U$.

lemma[Quantile Representation] The following statements hold under Assumptions (ref) - (ref): \begin{enumerate} [label = (\roman*)] • $U \perp R$, $U | R \sim \text{Unif}[0,1]$, and $\mathbb{P}(T \leq h(R,u) | R) = u, u \in [0,1]$. • For each $R=r$, $h_0$ and $h_1$ are strictly increasing in $U$. • The functions $h_0$ and $h_1$ are continuous. • The conditional distribution function $F_{\varepsilon|U,R}(e | u,r)$ is strictly increasing in $e$ and is continuous in $r$ at $\bar{r}$, that is, \begin{align*} \lim_{r \uparrow \bar{r}}F_{\varepsilon|U,R}(e | u,r) = \lim_{r \downarrow \bar{r}}F_{\varepsilon|U,R}(e | u,r), for every (e,u). \end{align*} \end{enumerate}

After the normalization, $U$ is independent of $R$ but $U$ and $\epsilon$ are possibly correlated even after conditioning on $R$. Let $F_{Y | T, R}$ be the conditional distribution function of $Y$ given $T$ and $R$. We define

align*[align* omitted — 446 chars of source]

The above left and right limits exist in view of Assumptions (ref) and (ref). The following assumption states that the support of the unobserved $\varepsilon$ does not vary with $U$ or $R$. This invariance of the support is not strong since it still allows $\varepsilon$ to be correlated with $U$ or $R$ in any way.

assumption[Support Invariance] Supp$(\varepsilon | U=u, R=r)$ does not depend on $u$ or $r$ in the neighborhood of $\bar{r}$. This common support is denoted by $\mathcal{E}$.

Nonparametric identification

We derive an important implication of the model ((ref)) - ((ref)). This implication is the key to identification and estimation. A function $g \in \mathcal{G}$ is said to satisfy Condition ((ref)) if for every $e \in \mathcal{E}$ and $ u \in [0,1]$,

align[align omitted — 192 chars of source]

If the function $g$ in Condition ((ref)) is equal to the true $g^*$, then the left-hand side of ((ref)) is equal to the conditional distribution of $\varepsilon$ given $U$ evaluated from the left side of the cutoff $\bar{r}$. Symmetrically, the right-hand side of ((ref)) is equal to the conditional distribution of $\varepsilon $ given $U$ evaluated from the right side of the cutoff $\bar{r}$. Then the equality holds by the continuity of $F_{\varepsilon|U,R}(e|u,\cdot)$ stated in Lemma (ref)(iv). We summarize this result in the following lemma.

lemma[Local Control Function] Under Assumptions (ref) - (ref), $g^*$ satisfies Condition ((ref)) with both sides of the equation equal to $F_{\varepsilon|U,R}(e| u,\bar{r})$.

When $g = g^*$, Condition ((ref)) can be written as

align*[align* omitted — 155 chars of source]

This leads to another interpretation of Lemma (ref): $U$ can serve as a control function local to the cutoff. After fixing the value of $U$, the variation in the treatment $T$ becomes locally exogenous. This is because given $U$ and $R$, the treatment $T$ becomes deterministic. The only variation left in $T$ around the cutoff is due to the discontinuity in the treatment choice function. Lemma (ref) is essentially a version of Lemma 1(i) in dong2021regression. From the IV perspective, Lemma (ref) corresponds to Theorem 1 in imbens2009identification. It is also similar to Theorem 1 in torgovitsky2015identification in that it provides a (necessary) characterization of the identified set of the structural function.\footnote{The identified set can be defined as the subset of $\mathcal{G}$ that contains the functions $g$ that can generate the observed distribution of $(Y,T,R)$. However, it is rather a detour to formally define such a set because in Section (ref) we directly use Condition ((ref)) for estimation.}

For any $g \in \mathcal{G}$, Lemma (ref) can be used to verify whether $g = g^*$. In particular, if

align*[align* omitted — 170 chars of source]

then $g$ can not be the true structural function. We further introduce some regularity conditions below.

assumption\ \begin{enumerate} [label = (\roman*)] • (Fuzzy RD). The support of $T | R$ from just below and above the cutoff are intervals denoted respectively by $\text{Supp}(h_0(\bar{r},U)) = [t_0',t_0''] \text{ and } \text{Supp}(h_1(\bar{r},U)) = [t_1',t_1''].$ The two supports are overlapping: $[t_0',t_0''] \cap [t_1',t_1''] \ne \emptyset$.\footnote{Infinite intervals are also allowed. For example, $\text{Supp}(h_0(\bar{r},U))$ can be $(-\infty,t_0'']$, $[t_0',\infty)$, or $\mathbb{R}$. We use the notation $[t_0',t_0'']$ to represent all these cases. } • (Strong Discontinuity). Local to the cutoff, the functions $h_0$ and $h_1$ intersects and only intersects finitely many times. That is, the following set is nonempty and finite: \begin{align*} \{ h_0(\bar{r},u): h_0(\bar{r},u) = h_1(\bar{r},u) \in [t_0',t_0”] \cap [t_1',t_1”], u \in [0,1] \}. \end{align*} \end{enumerate}

Assumption (ref) imposes restrictions on the nature of the discontinuity. Assumption (ref)(i) requires that the RD design is fuzzy in that there are treatment levels that are taken both below and above the cutoff. Assumption (ref)(ii) imposes restrictions on the strength of the discontinuity. It requires that the conditional quantile functions $h_0(\bar{r},\cdot)$ and $h_1(\bar{r},\cdot)$ only intersects finitely many times. The two curves can intersect but not overlap. If the two functions $h_0(\bar{r},\cdot)$ and $h_1(\bar{r},\cdot)$ overlaps on some interval, then the structural function is not identified on that interval because there is no exogenous variation in the treatment inside that interval.\footnote{In that case, it is possible to partially identify the structural function. } In the extreme case where the two curves completely overlap, there is no discontinuity.

With the above assumptions, we present the main identification result of the paper. The following theorem shows that Condition ((ref)) identifies the true structural function up to a monotone transformation of the error term.

theorem[Nonparametric Identification] Let Assumptions (ref) - (ref) hold. If $g \in \mathcal{G}$ satisfies Condition ((ref)), then there exists a continuous and strictly increasing function $\lambda^g$ such that for every $t \in [t_0',t_0''] \cup [t_1',t_1''], e \in \mathcal{E},$ and $u \in [0,1]$, $g^*(t,\bar{r},e) = g(t,\bar{r},\lambda^g(e))$.

remarkTheorem (ref) is proved based on the sequencing approaching developed in the proof of Theorem 2 in torgovitsky2015identification.

Theorem (ref) is the best one can achieve in terms of identifying the nonseparable structural function because the error term is unobserved. Any $g \in \mathcal{G}$ that satisfies Condition ((ref)) is equally good as the true $g^*$. The only difference is that the error term is rescaled by the monotone transformation $\lambda^g$. Therefore, such a function $g$ can also be seen as a “version” of $g^*$.

Inspecting the conditions of Theorem (ref), we can see that no independence assumption is needed. This is why the RD design is often considered a more credible approach than instrumental variables for conducting causal inference. However, in many studies of RD designs, a local independence assumption is imposed, explicitly making the running variable exogenous around the cutoff. For example, Assumption A3(i) in hahn2001identification requires $(\varepsilon,U)$ to be jointly independent of $R$ conditioning on $R$ near $\bar{r}$. In the binary treatment case, dong2018alternative shows that this local independence condition is not needed to achieve identification.

Based on the observation made in Lemma (ref), we can recover the conditional distribution of $\varepsilon$ given $U$ and $R = \bar{r}$. For any $g \in \mathcal{G}$, if $g$ is the true structural function, then the corresponding conditional distribution of $\varepsilon$ is

align[align omitted — 229 chars of source]

In fact, the above conditional distribution $F_{\varepsilon | U,R}^g$ is a transformed version of the true conditional distribution $F_{\varepsilon | U,R}$, where the transformation is the $\lambda^g$ defined in Theorem (ref). This means that the conditional distribution of $\varepsilon | U,R=\bar{r}$ is identified up to the same monotone transformation as the structural function. To eliminate such inconvenience caused by the error term, we can integrate out $\varepsilon$ and obtain a unique conditional average structural function (CASF):

align*[align* omitted — 77 chars of source]

where the expectation is taken with respect to the true conditional distribution of $\varepsilon$ given $ R=\bar{r}$. The following corollary summarizes the above discussion.

corollary[CASF] Let Assumptions (ref) - (ref) hold. For any $g^{\circ} \in \mathcal{G}$ that satisfies Condition ((ref)), let $\lambda^g$ be the transformation defined in Theorem (ref). The following two statements hold true.\footnote{Since the conditional distribution of $\varepsilon$ given $U,R=\bar{r}$ is identified. We can also identify other structural parameters, including the conditional quantile structural function.} \begin{enumerate} [label = (\roman*)] • For every $e \in \mathcal{E}$ and $u \in [0,1]$, $F^{g^{\circ}}_{\varepsilon | U,R}(\lambda^{g^{\circ}}(e) | u,\bar{r}) = F_{\varepsilon | U,R}(e | u,\bar{r}).$ • For any $t \in [t_0',t_0''] \cup [t_1',t_1'']$, the CASF $\beta^*(t,\bar{r})$ is uniquely identified as \begin{align*} \beta^*(t) = \int {g^{\circ}}(t,\bar{r},e) dF^{g^{\circ}}_{\varepsilon | U,R}(e | u,\bar{r}) du. \end{align*} \end{enumerate}

The CASF $\beta^*(t)$ gives the average outcome the policy-maker can achieve when the treatment level for individuals with characteristic $R = \bar{r}$ is set to $t$. It is worth noting the difference between the CASF and the local average structural function (LASF) commonly seen in the LATE literature. The LASF represents the average outcome for the so-called compliers, an unobservable subpopulation. Therefore, the policy-maker cannot assign treatment to the compliers even when the LASF is identified. On the other hand, the identified CASF can directly guide the treatment assignment to the subpopulation with $R = \bar{r}$. The derivative of the CASF is not the causal effect specific to any subpopulation. Following the spirit of, for example, heckman2001policy, we may call CASF a policy-relevant parameter.

Semiparametric estimation

In this section, we consider parametrizations of the structural function that maintain nonlinearity and nonseparability. We propose a semiparametric estimation procedure and derive its large-sample properties. The estimator is semiparametric because the structural function is parametrically specified, while the treatment choice model is left nonparametrically specified.

We do not consider a fully nonparametric estimator since such a procedure can be too data-demanding for practical use, which is especially true for the RD design since the estimation is in the local neighborhood of the cutoff.\footnote{From the theoretical perspective, it can be challenging to construct a fully nonparametric estimator. If we follow the sieve approach, for example, we would need to consider a basis of functions that are strictly increasing in one of the arguments to accommodate the monotonicity of the structural function, which is a non-trivial task.}

Construction of the estimator

Consider the following parametrization of $\mathcal{G}$ local to the cutoff.

assumption[Semiparametric Specification] There is a one-to-one mapping from the class $\{ (t,e) \mapsto g(t,\bar{r},e): g \in \mathcal{G} \}$ of functions to a finite-dimensional parameter space $\Gamma \subset \mathbb{R}^{d_\Gamma}$. We write such parametrization as $\{g_\gamma(\cdot,\bar{r},\cdot): \gamma \in \Gamma\}$. Assume this parametric model is correctly specified, that is, there exists $\gamma^* \in \Gamma$ such that $g_{\gamma^*}(\cdot,\bar{r},\cdot) = g^*(\cdot,\bar{r},\cdot)$.

assumption[Normalization of $\Gamma$] For any $\gamma,\gamma' \in \Gamma$, if there exists a transformation $\lambda$ such that \begin{align*} g_\gamma(\cdot,\bar{r},\cdot) = g_{\gamma'}(\cdot,\bar{r},\lambda(\cdot)), \end{align*} then $\gamma = \gamma'$ and $\lambda$ is the identity transformation.

Assumption (ref) is a normalization condition that fixes the scale of the error term $\varepsilon$. An example is provided below to illustrate the parametrization of the structural function. One way to achieve such normalization is to have some treatment value $\tilde{t}$ such that $g_\gamma(\tilde{t},\bar{r},e) = e, \text{ for all } \gamma \in \Gamma.$

exampleLet $\tilde{t} \in [t_0',t_0''] \cap [t_1',t_1'']$. We can specify the structural function by \begin{align*} g_\gamma(T,\bar{r},\varepsilon) = \gamma_1 (T-\tilde{t}) + \gamma_2 (T-\tilde{t})^2 + \gamma_3 (T-\tilde{t})\varepsilon + \varepsilon. \end{align*} The parameter $\gamma = (\gamma_1,\gamma_2,\gamma_3)$ is three-dimensional. The function $g_\gamma(T,\bar{r},\varepsilon)$ is strictly increasing in $\varepsilon$ for all $\gamma$ satisfying $\mathbb{P}(\gamma_3(T-\tilde{t})+1 > 0)=1$. The parametrization satisfies Assumption (ref) because, by construction, $g_{\gamma}(\tilde{t},\bar{r},e) = e$ for all values of $\gamma$ and $e$.\footnote{This normalization strategy is presented in Equation (2.5) of matzkin2003nonparametric.} The model is quadratic in $T$ and nonseparable between $T$ and $\varepsilon$. The effect of $T$ on $Y$ is allowed to be nonlinear and contain unobserved heterogeneity. The distribution of $\varepsilon$ is not parametrized and thus can be very general.

The true parameter $\gamma^*$ in the normalized semiparametric model can be identified as follows. We use $h^* = (h_0^*,h_1^*)$ to signify the true conditional quantile functions and $h = (h_0,h_1)$ a generic pair of conditional quantile functions. Let $w(e,u)$ be a weighting function defined on $\mathbb{R} \times [0,1]$. Define the criterion function as

align*[align* omitted — 155 chars of source]

where $D_{\gamma,h}(e,u)$ is defined to be

align[align omitted — 230 chars of source]

This criterion function is based on Equation ((ref)), which by Lemma (ref) is a necessary characterization of the identified set. We take an integral form of Condition ((ref)) because it gives a faster convergence rate of the resulting estimator.

assumption[Weighting function $w$] The function $w$ is nonnegative, integrates to one, and is bounded on $\mathbb{R} \times [0,1]$. The support of $w$ contains $\mathcal{E} \times [0,1]$.

In practice, we may use a weighting function $w$ that is supported on the entire domain $\mathbb{R} \times [0,1]$ since $\mathcal{E}$ is unknown. The following corollary provides the semiparametric identification result, which is based on the nonparametric identification result in Section (ref). It shows that the criterion function, when evaluated at the true nuisance parameter value $h^*$, is uniquely minimized by the true $\gamma^*$.

corollary[Semiparametric Identification] Let Assumptions (ref) - (ref) hold. For any $\gamma \in \Gamma$ such that $\gamma \ne \gamma^*$, we have $\left\lVert D_{\gamma,h^*} \right\rVert_w > \left\lVert D_{\gamma^*,h^*} \right\rVert_w = 0.$ In other words, $\gamma^*$ is the unique minimizer of $\left\lVert D_{\gamma^*,h^*} \right\rVert_w$.

Assume there is an independent and identically distributed (iid) sample $(Y_i,T_i,R_i)_{i=1}^n$ available. We propose an estimation procedure based on the above semiparametric identification result. The idea is that we first estimate the nonparametric components $(h_0,h_1)$ and $(F^-_{Y|T,R},F^+_{Y|T,R})$ that appear in the criterion function. Then we construct an empirical version of the criterion function and take its minimizer to be the estimator.

The estimation procedure of $\gamma$ is more specifically divided into three steps. The first step is to estimate the conditional quantile functions $h_0$ and $h_1$. The second step uses local linear regression (LLR) to estimate the conditional distributions $F^-_{Y | T,R}$ and $F^+_{Y | T,R}$. It is standard to use local polynomials in the estimation of RD designs porter2003estimation,sun2005adaptive. The difference is that classical RD methods use local polynomial to estimate the conditional expectation function of $Y$ given $R$ while we estimate the conditional distribution of $Y$.\footnote{Local linear estimation of the conditional distribution function can be found in hansen2004nonparametric,xie2021uniform.} The third step constructs an estimate of the criterion function by replacing the nonparametric nuisance parameters in ((ref)) by their estimated counterparts and then finds the estimate of $\gamma^*$ by minimizing the estimated criterion function. We describe the detail of the estimation procedure as follows. Denote $\mathcal{Y}$ as the range of the outcome $Y$.

itemize• STEP 1. Choose estimators $\hat{h}_0(\bar{r},\cdot)$ and $\hat{h}_1(\bar{r},\cdot)$ of the corresponding conditional quantile processes, $h_0(\bar{r},\cdot)$ and $h_1(\bar{r},\cdot)$. Specific constructions are provided in Section (ref). • STEP 2. Choose two bandwidth sequences $b_1 = b_{1n}$ and $b_2 = b_{2n}$ and three kernel functions $k_Y$, $k_T$, and $k_R$. Define $K_Y(y) = \int_{-\infty}^y k_Y(\tilde{y}) d \tilde{y} $. For each $y \in \mathcal{Y}$ and $t \in [t_0',t_0'']$, solve the following minimization problem: \begin{align*} \min_{a^-,a_{T}^-,a_{R}^-} \sum_{i:R_i < \bar{r}} & \left( K_Y\left( \frac{y - Y_i}{b_2} \right) - a^- - a_{T}^-(T_i-t) - a_{R}^-(R_i - \bar{r}) \right)^2 \\ & \times k_T\left( \frac{T_i - t}{b_1} \right) k_R\left( \frac{R_i - \bar{r}}{b_1} \right). \end{align*} The minimizer $\hat{a}^-$ is the estimate $\hat{F}^-_{Y | T,R}(y | t,\bar{r})$. For each $y \in \mathcal{Y}$ and $t \in [t_1',t_1'']$, solve the following minimization problem: \begin{align*} \min_{a^+,a_{T}^+,a_{R}^+} \sum_{i:R_i \geq \bar{r}} & \left( K_Y\left( \frac{y - Y_i}{b_2} \right) - a^+ - a_{T}^+(T_i-t) - a_{R}^+(R_i - \bar{r}) \right)^2 \\ & \times k_T\left( \frac{T_i - t}{b_1} \right) k_R\left( \frac{R_i - \bar{r}}{b_1} \right). \end{align*} The minimizer $\hat{a}^+$ is the estimate $\hat{F}^+_{Y | T,R}(y | t,\bar{r})$. • STEP 3. Construct the empirical version of the criterion function: \begin{align*} \left\lVert \hat{D}_{\gamma,\hat{h}} \right\rVert_{w} = \left( \int_0^1 \int_{\mathbb{R}} | \hat{D}_{\gamma,\hat{h}}(e,u)|^2 w(e,u) dedu \right)^{1/2} , \end{align*} where $\hat{D}_{\gamma,\hat{h}}(e,u) $ is defined to be \begin{align*} \int_0^u \left( \hat{F}^-_{Y|T,R}( g_\gamma(\hat{h}_0(\bar{r},v),\bar{r},e) | \hat{h}_0(\bar{r},v),\bar{r} ) - \hat{F}^+_{Y|T,R}( g_\gamma(\hat{h}_1(\bar{r},v),\bar{r},e) | \hat{h}_1(\bar{r},v),\bar{r} ) \right) dv. \end{align*} The estimator $\hat{\gamma}$ is any parameter value in $\Gamma$ that satisfies \begin{align} \big\lVert \hat{D}_{\hat{\gamma},h} \big\rVert_{w} \leq \inf_{\gamma \in \Gamma} \big\lVert \hat{D}_{\gamma,h} \big\rVert_{w} + O_p\left(\alpha_n\right), \end{align} where $\alpha_n \rightarrow 0$ is specified by Equation ((ref)) in Appendix (ref).

In the third step, we construct the empirical version of the criterion function:

align*[align* omitted — 153 chars of source]

where

align*[align* omitted — 247 chars of source]

The estimator $\hat{\gamma}$ is defined as any parameter value that satisfies

align[align omitted — 215 chars of source]

\fi

Asymptotic normality

More regularity assumptions are imposed for the estimator to enjoy desirable statistical properties. Define

alignat*{3} f^-_{T | R}(t | r) & = \begin{cases} f_{T | R}(t | r), & if r < \bar{r}, \\ \lim_{r \uparrow \bar{r}} f_{T | R}(t | r), & if r = \bar{r}. \end{cases} \quad f^+_{T | R}(t | r) & = \begin{cases} f_{T | R}(t | r), & if r > \bar{r}, \\ \lim_{r \downarrow \bar{r}} f_{T | R}(t | r), & if r = \bar{r}. \end{cases}

The above left and right limits exist in view of Assumptions (ref) and (ref).

assumption[Distributions of $Y,T,$ and $R$] \ \begin{enumerate} [label = (\roman*)] • The support of $T$ does not vary with $R$ except when crossing the cutoff $\bar{r}$, i.e., Supp$(T | R=r) = [t_0',t_0'']$ for $r < \bar{r}$ and Supp$(T | R=r) = [t_1',t_1'']$ for $r > \bar{r}$. The density functions $f^-_{T,R}$ and $f^+_{T,R}$ are bound away from zero. • The density functions $f^-_{T,R}$ and $f^+_{T,R}$ are twice continuously differentiable, and $\frac{\partial^2}{\partial t^2} f^-_{T,R}(t,\bar{r})$ and $\frac{\partial^2}{\partial t^2} f^+_{T,R}(t,\bar{r})$ are Lipschitz continuous with respect to $t$. • The support of $Y$, $\mathcal{Y}$, is compact. The conditional distribution functions $F^-_{Y | T,R}$ and $F^+_{Y | T,R}$ are three-times continuously differentiable over $\mathcal{Y} \times [t_0',t_0''] \times [r_0,\bar{r}]$ and $\mathcal{Y} \times [t_1',t_1''] \times [\bar{r},r_1]$, respectively. \end{enumerate}

assumption[Complexity of the Parametric Model] The parametrization $\{g_\gamma(\cdot,\bar{r},\cdot):\gamma \in \Gamma\}$ satisfies the following conditions. \begin{enumerate} [label = (\roman*)] • The parameter space $\Gamma$ is compact. • The class of functions $\{ T \mapsto g_\gamma(T+v,\bar{r},e): \gamma \in \Gamma, v \in (-1,1), e \in \mathcal{E} \}$ is finite-dimensional. • The function $g_\gamma(t,\bar{r},e)$ is twice continuously differentiable over $\gamma \in \Gamma$, $t \in [t_0',t_0''] \cup [t_1',t_1'']$, and $ e \in \mathcal{E}$. • The gradient $\nabla_\gamma D_{\gamma^*,h^*}(e,u)$ is a vector of linearly independent functions of $(e,u)$. \end{enumerate}

assumption[Kernels] \ \begin{enumerate} [ label = (\roman*)] • The kernel functions $k_T$ and $k_R$ are (1) supported on $[-1,1]$, (2) strictly greater than zero in the interior of the support, (3) of bounded variation, (4) continuously differentiable on $\mathbb{R}$. • The kernel function $k_Y$ is (1) nonnegative and (2) integrable on $\mathbb{R}$ with $\int k_Y(y) dy =1$ and satisfies (3) $\int y k_Y(y) dy = 0$. \end{enumerate}

assumption[Bandwidth] The bandwidth $b_1$ and $b_2$ satisfy the following conditions: \begin{enumerate} [label = (\roman*)] • $b_1 \asymp b_2$.\footnote{The notation $b_1 \asymp b_2$ means that there exists $C>1$ such that $b_1/b_2 \in [1/C,C]$.} • $(n \log n) b_1^6 = o(1)$. • $n b_1^{\frac{13}{3} + \epsilon} \rightarrow \infty$, for some sufficiently small $\epsilon>0$. \end{enumerate}

assumption[First-step Conditional Quantile Estimators] The estimators $\hat{h}_0$ and $\hat{h}_1$ satisfy the following conditions. \begin{enumerate} [label = (\roman*)] • Monotonicity and smoothness: for every $n$ sufficiently large, there exist $C > 0$ and deterministic and finite partitions $\mathcal{P}_0^n$ and $\mathcal{P}_1^n$ on $(0,1)$ such that \begin{align*} \mathbb{P}\left(\hat{h}_0(\bar{r},\cdot) \notin \mathcal{H}_0(\mathcal{P}_0^n) \right), \mathbb{P}\left(\hat{h}_1(\bar{r},\cdot) \notin \mathcal{H}_1(\mathcal{P}_1^n) \right) = O \big( \sqrt{b_1}\big), \end{align*} where \begin{align*} \mathcal{H}_0(\mathcal{P}_0^n) = &\{ function h from [0,1] into [t_0',t_0”]: on each element of $\mathcal{P}_0^n$, h is strictly \\ & \quad increasing, its inverse $h^{-1}$ is three-times continuously differentiable, \\ & \quad \text{and $(h^{-1})^{(3)}$ is Lipschitz continuous } \}, \end{align*} and $\mathcal{H}_1(\mathcal{P}_1^n)$ is defined analogously by replacing $\mathcal{P}_0^n$ with $\mathcal{P}_1^n$. • Uniform Bahadur representation: \begin{align*} \hat{h}_0(\bar{r},u) - h_0^*(\bar{r},u) & = b_1^2 \nu_0(u) + O_p ( b_1^3)\\ & \quad + \frac{1}{nb_1} \sum_{i=1}^n q_{0}(T_i,R_i;u) k_{Q,0}\left( \frac{R_i - \bar{r}}{b_1} \right) \mathbf{1}\{R_i < \bar{r}\} + o_p \big( 1/\sqrt{nb_1} \big), \\ \hat{h}_1(\bar{r},u) - h_1^*(\bar{r},u) & = b_1^2 \nu_1(u) + O_p ( b_1^3)\\ & \quad + \frac{1}{nb_1} \sum_{i=1}^n q_{1}(T_i,R_i;u) k_{Q,1}\left( \frac{R_i - \bar{r}}{b_1} \right) \mathbf{1}\{R_i \geq \bar{r}\} + o_p \big(1/\sqrt{nb_1} \big), \end{align*} uniformly over $u \in (0,1)$. The functions $\nu_0$ and $\nu_1$ are bounded. The functions $q_0$ and $q_1$ are (1) bounded, (2) centered, that is, $\mathbb{E}[q_{0}(T,R;u)|T,R] = \mathbb{E}[q_{1}(T,R;u)|T,R] = 0$, and (3) does not vary with $n$. The functions $k_{Q,0}$ and $k_{Q,1}$ are bounded. • Uniform convergence rate: \begin{align*} \lVert \hat{h} - h^* \rVert_\infty & = \sup_{ u \in (0,1)} |\hat{h}_0(\bar{r},u) - h_0^*(\bar{r},u)| \vee |\hat{h}_1(\bar{r},u) - h_1^*(\bar{r},u)| \\ & = O_p \left( \sqrt{\log n / (n b_1)} + b_1^2 \right). \end{align*} \end{enumerate}

A brief discussion of the assumptions is in order. Assumption (ref) imposes smoothness restrictions on the joint distribution of $(Y,T,R)$. In the previous section, the identification result only requires continuity of the relevant functions. For estimation, we need higher-order smoothness regarding the distribution functions. Assumption (ref) imposes restrictions on the parametric model of the structural function. Part (ii) restricts the complexity of the model. Part (iii) imposes high-order smoothness on the structural function. Part (iv) is similar to Assumption D4 in TORGOVITSKY2017minimum and requires that $\nabla_\gamma D_{\gamma^*,h^*}$ to carry information about each component of the parameter.

Assumption (ref) imposes restrictions on the kernel functions $k_T$, $k_R$, and $k_Y$. The differentiability is needed to prove a stochastic equicontinuity condition. Assumption (ref) restricts that $b_1$ and $b_2$ are of the same asymptotic order, which is slightly faster than $n^{-1/6}$ and slightly slower than $n^{-3/13}$. This assumption is not restrictive and allows for the asymptotic mean squared error (AMSE) optimal bandwidth as well as undersmoothing.

Assumption (ref) imposes high-level restrictions on the first-stage nonparametric conditional quantile estimators. Part (i) assumes that the quantile estimators are piece-wise monotonic and smooth with a high probability. Part (ii) and (iii) give the uniform Bahadur representation and the uniform convergence rate, which are fairly standard in the quantile estimation literature. In Section (ref), we discuss a specific nonparametric quantile estimator that satisfies Assumption (ref).

theorem[Asymptotic Distribution of the Semiparametric Estimator] Let Assumptions (ref) - (ref) hold. Then $\left\lVert\hat{\gamma} - \gamma^*\right\rVert_2 = O_p(b_1^2 + 1/\sqrt{nb_1})$ and \begin{align*} \big( \sqrt{nb_1} (\Sigma_- + \Sigma_+)^{-1/2} \big) (\Delta (\hat{\gamma} - \gamma^*) - b_1^2( B_- - B_+)) \overset{d}{\rightarrow} N(0,\bm{I}_{d_\Gamma}), \end{align*} where $I_{d_\Gamma}$ is the $d_\Gamma$-dimensional identity matrix. The exact forms of $\Delta$, $B_-$, $B_+$, $\Sigma_-$ and $\Sigma_+$ are given in Equations ((ref)), ((ref)), ((ref)), ((ref)) and ((ref)) in Appendix (ref), respectively.

remarkThe convergence rate of $\hat{\gamma}$ is $b_1^2 + 1/\sqrt{nb_1}$, which is equal to $n^{-2/5}$ when $b_1 \asymp n^{-1/5}$. This rate is the same as the one obtained in the classical RD design with a binary treatment variable hahn2001identification. Having a continuous treatment does not slow down the convergence rate of the estimator in this case. This is due to the integral smoothing in the definition of the criterion function $D_{\gamma,h}$ in ((ref)).

remarkThe proof of Theorem (ref) follows the general steps of proving asymptotic normality of semiparametric estimators as in, for example, TORGOVITSKY2017minimum and chen2003estimation. The main difficulty is that the usual stochastic equicontinuity condition is not sharp because the criterion function is nonparametrically estimated. To overcome this issue, we first use the empirical process theory to derive a uniform convergence rate for the estimated criterion function, which gives an initial bound on the convergence rate of $\hat{\gamma}$. A sharper stochastic equicontinuity result is then derived based on this initial bound together with more applications of the empirical process theory. This sharper stochastic equicontinuity result helps demonstrate that the usual linearization of the criterion function is valid. See the proof in Appendix (ref) for details.

Once the asymptotic normal distribution of the estimator $\hat{\gamma}$ is established, we can conduct inference for $\gamma^*$. Theorem (ref) together with the undersmoothing condition that $nb_1^5 = o(1)$ gives that

align*[align* omitted — 139 chars of source]

A linear null hypothesis regarding $\gamma$ can be written as $H \gamma = \eta$, where $\eta \in \mathbb{R}^{d_\eta}$ and $H$ is a $d_{\eta} \times d_\gamma$ full-rank matrix. Consider the test statistic

align*[align* omitted — 177 chars of source]

where $\hat{\Delta}$, $\hat{\Sigma}_-$, and $\hat{\Sigma}_+$ are consistent estimators of $\Delta$, $\Sigma_-$, and $\Sigma_+$, respectively. By Slutsky's theorem, the above test statistic converges in distribution to the $\chi^2$ distribution with ${d_\eta}$ degrees of freedom. In Appendix (ref), we discuss how to construct consistent estimators for $\Delta$, $\Sigma_-$, and $\Sigma_+$.

First-step nonparametric quantile estimators

This section discusses how to construct nonparametric conditional quantile estimators that satisfy Assumption (ref). Consider the following two-step estimation procedure introduced by QU2015nonparametric. Define $\rho_u(t) = t(u - \mathbf{1}\{t < 0\})$.

itemize• STEP 1. Choose a bandwidth sequence $b_3 = b_{3n} = o(1)$ and a kernel function $k_{\textit{FS}}$. Partition the unit interval $(0,1)$ into a grid of equally spaced points $\{u_1,\cdots,u_{J_n}\}$, where $J_n/(nb_3)^{1/4} \rightarrow \infty$. Solve the following optimization problem: \begin{align} \min_{\{\mathrm{h}_j,\mathrm{h}'_j\}_{j=1}^{J_n}} \sum_{j=1}^{J_n} \sum_{i=1}^n \rho_{u_j} \left( T_i - \mathrm{h}_j - \mathrm{h}_j' (R_i - \bar{r}) \right) k_{FS} \left( \frac{R_i - \bar{r}}{b_3} \right) \mathbf{1}\{ R_i < \bar{r} \}. \end{align} Denote the minimizers by $(\hat{h}_0(\bar{r},u_1),\cdots,\hat{h}_0(\bar{r},u_{J_n}))$. • STEP 2. Let $u_0 = 0$ and $u_{J_{n+1}}=1$. Let $\hat{h}_0(\bar{r},u_0) = \min_{i:R_i < \bar{r}} T_i$ and $\hat{h}_0(\bar{r},u_{J_n + 1}) = \max_{i:R_i < \bar{r}} T_i$. Linearly interpolate between the estimates to obtain an estimate for the entire quantile process. That is, for any $u \in (u_j,u_{j+1})$, define \begin{align*} \hat{h}_0(\bar{r},u) = \frac{u_{j+1} - u}{u_{j+1} - u_j} \hat{h}_0(\bar{r},u_j) + \frac{u - u_j}{u_{j+1} - u_j} \hat{h}_0(\bar{r},u_{j+1}). \end{align*}

The estimator $\hat{h}_1(\bar{r},\cdot)$ can be analogously defined by using the data with $R_i \geq \bar{r}$.\footnote{There are three estimators of conditional quantile process in QU2015nonparametric. The estimator explained here is their second one, denoted by $\hat{\alpha}^*$ in that paper. Their third estimator imposes a monotonicity constrain to the minimization problem ((ref)).} We can verify Assumption (ref) for the estimator constructed above. Denote

align*[align* omitted — 195 chars of source]

propositionLet Assumptions (ref)(iv) and (ref)(i)-(ii) hold. Assume that the third-order derivatives $\frac{\partial^3}{\partial R^3} h_0^*(r,u)$ and $\frac{\partial^3}{\partial R^3} h_1^*(r,u)$ are Lipschitz continuous respectively on $[r_0,\bar{r}] \times [0,1]$ and $[\bar{r},r_1] \times [0,1]$. Assume that the bandwidth $b_3 = cb_1$ for some constant $c>0$. Assume that the kernel $k_{\textit{FS}}$ is nonnegative, of bounded variation, compactly supported, having finite first-order derivatives and satisfying \begin{align*} \int k_{FS}(x) dx= 1, \int x k_{FS}(x) dx= 0, \int x^2 k_{FS}(x) dx < \infty. \end{align*} Then the estimators $\hat{h}_0(\bar{r},\cdot)$ and $\hat{h}_1(\bar{r},\cdot)$ described above satisfy Assumption (ref). The specific forms of $\nu_0$ and $\nu_1$ are \begin{align*} \nu_0(u) & = \frac{c^2}{2} \frac{\partial^2}{\partial r^2} h_0^*(\bar{r},u) \iota' \Omega_{Q,0}^{-1} \int x^2 (1,x)' k_{FS}(x) \mathbf{1}\{x < 0\} dx, \\ \nu_1(u) & = \frac{c^2}{2} \frac{\partial^2}{\partial r^2} h_1^*(\bar{r},u) \iota' \Omega_{Q,1}^{-1} \int x^2 (1,x)' k_{FS}(x) \mathbf{1}\{x \geq 0\} dx. \end{align*} The specific forms of $k_{Q,0}$ and $k_{Q,1}$ are \begin{align*} k_{Q,j}(x) = \iota' \Omega_{Q,j}^{-1} (1,x/c)' k_{FS}(x/c)/c, j =0,1. \end{align*} The functions $q_0$ and $q_1$ are \begin{align*} q_0(T,R;u) = (u - \mathbf{1}\{ T \leq h_0^*(\bar{r},u) \}) / (f_R(\bar{r}) f^-_{T | R}(h_0^*(\bar{r},u) | \bar{r}) ), \\ q_1(T,R;u) = (u - \mathbf{1}\{ T \leq h_1^*(\bar{r},u) \}) / (f_R(\bar{r}) f^+_{T | R}(h_1^*(\bar{r},u) | \bar{r}) ), \end{align*} which take the form of an influence function for quantiles.

Other quantile estimation methods are also available. For example, one can consider the generic framework proposed by chernozhukov2010quantile for rearrangement. In particular, they show that the rearrangement of a preliminary estimated quantile process delivers a monotonic estimator that preserves the asymptotic properties. This result gives a different way to generate estimators that satisfy Assumption (ref). We can start with an estimator with desired asymptotic properties that give rise to Assumption (ref)(ii) and (iii), and then apply the rearrangement procedure. The resulting estimator would be monotonic on the entire domain, and partitioning is unnecessary.

Numerical results

This section presents the empirical application and the simulation studies. The empirical study shows that the semiparametric estimator is considerably better than the simple TSLS estimator in discovering quantitative information regarding the structural function. The simulation studies show that the semiparametric procedure can accurately estimate the parameters with a moderate sample size.\footnote{Replication files for the empirical and simulation studies are available from the author upon request.}

Empirical application

In the empirical study, we examine the causal effect of sleep time on health status by exploiting the discontinuity in the timing of natural light at time zone boundaries. The unit of observation is the individual in the American Time Use Survey (ATUS), the outcome $Y$ is the individual's health status measured by the body-mass index (BMI),\footnote{BMI is a person's weight in kilograms divided by the square of height in meters. The Centers for Disease Control and Prevention define overweight as BMI > 25 and obesity as BMI > 30.} the treatment $T$ is the sleep time, and the running variable $R$ is the longitudinal distance to the nearest time zone boundary, with cutoff $\bar{r} = 0$ denoting the time zone border. As explained in the introduction, the identification is based on the exogenous variation in the sleep time around the time zone boundary. This exogenous variation is due to the difference in the timing of natural light on each side of the time zone boundary.

Many studies in the medical literature examine the effect of sleep time on overweight issues. See beccuti2011sleep and the references therein. These studies typically use survey or laboratory data. This problem is first studied by using the RD design in Giuntella2019sunset.\footnote{Giuntella2019sunset study many health and economics-related issues. Here we only mention the relevant ones.} The relevant outcome variable they use is a binary indicator of the obesity (or overweight) status indicating whether the BMI is above some threshold. They use the TSLS procedure to estimate a linear structural function. We consider two improvements based on their work. First, we directly use BMI as the outcome variable, providing a more quantitative measure of the health status. Second, we use the proposed semiparametric estimator to estimate a nonlinear structural function. Previous medical studies have provided evidence of the nonlinearity of the structural function. For example, hairston2010sleep show that both undersleeping and oversleeping lead to an increase in BMI while sleeping around 8 hours leads to a more healthy BMI level.

The data for this empirical application is collected from IPUMS CPS ipums-cps and IPUMS ATUS ipums-atus during the periods 2006 - 2008 and 2014 - 2016. By linking these datasets, we can locate the county where the individual lives and then use the county's centroid as the location of the individual. We focus on counties near the time zone boundary between the Eastern and Central time zone. The counties are divided into two regions based on their latitude. We estimate the model separately for each region.

The estimated marginal effects of sleep on BMI from the semiparametric estimator and the TSLS estimator are shown in Figure (ref). Several interesting findings are observed based on the semiparametric estimates. First, the marginal effects are increasing and increase from negative to positive. This lends some support to the previous argument that the structural function is nonlinear and neither sleeping too little nor too much is preferable. Second, we can determine the optimal (in terms of BMI) sleep time by finding the zero of the marginal effect curve. In both cases, the optimal sleep time is between 7 and 8 hours, which also aligns with the findings in previous medical studies. Third, the results from the two regions are similar, meaning that the variation across different latitudes is small.

From Figure (ref), we can also see that the TSLS estimates are not capable of demonstrating the above results. First, the TSLS procedure only provides a constant estimate of the marginal effect across all levels of sleep time. This means an extra hour of sleep would lead to the same effect on health regardless of the person's current sleep time, which is inappropriate in this setting. Moreover, we cannot estimate the optimal sleep time based on the linear structural function. Second, the magnitude of the TSLS estimates is small. This is because the TSLS provides a weighted average of the marginal effects across the entire range of sleep time. By averaging the negative and positive effects, the TSLS delivers an estimate attenuated toward zero, which is not informative for the researcher.\footnote{Giuntella2019sunset find a more significant effect of sleep time on obesity. There are two possible reasons: they consider the binary indicator of obesity, and they include more control variables in the regression.}

figure[figure omitted — 848 chars of source]

Simulations

We use simulation studies to investigate the performance of the proposed semiparametric method and compare it with the performance of the TSLS estimator. The data generating process (DGP) for these simulations was chosen to roughly approximate the ATUS data used in the empirical application. Let marginal distributions of $U$ and $\varepsilon$ are given by $F_U = \textit{Unif}(0,1)$ and $F_\varepsilon = \textit{Beta}(2,2)$, respectively. Two marginal distributions for $R$ are considered, $\textit{Unif}(0,1)$ and $N(0,1)$. The joint distribution of $(R,\varepsilon,U)$ be characterized by the Gaussian copula with correlation structure corr$(R,U) = 0$, corr$(\varepsilon,R) = \rho_R$, and corr$(\varepsilon,U) = \rho_{U}$. The treatment choice model is given by $h^*_0(r,u) = r + 2\sin(\pi u /2)$ and $h^*_1(r,u) = r + 2u^3$. The structural function is given by

align*[align* omitted — 149 chars of source]

The true $\gamma^*$ is taken to be $(1,1,1)$.

Further implementation details are described below. Construct a kernel function $k$ that is an even function given by

align*[align* omitted — 169 chars of source]

We can verify that $k$ is continuously differentiable on the real line and compactly supported on $[0,1]$. Within the interior of its support, $k$ is strictly positive. We use this function $k$ to be the kernels $k_T$, $k_R$, $k_Y$, and $k_{\textit{FS}}$ in the estimation. The bandwidth is chosen to be $b_1 = b_2 = b_3 = 2n^{-1/5}$. The weighting function $w(e,u)$ is chosen to be constant in $u$ and equal to the standard normal density function with respect to $e$.

Table (ref) contains the simulation results of the performance of the semiparametric estimator for different choices of the marginal distribution of $R$, the correlation parameters $(\rho_U,\rho_R)$ and the sample size. In each case, the number of replications is set at 500. We can see that the estimator performs well with a moderate sample size ($n=1000$). The marginal distribution of $R$ does not have a large impact on the performance. When the sample size is small, larger values of $\rho_R$ or $\rho_U$ can lead to poorer performance of the estimator. The plausible reason is that larger values of the correlation parameters would lead to more severe endogeneity issues in finite samples. When the sample size becomes large, the performance of the estimator does not vary significantly with the choices of $(\rho_R,\rho_U)$.

table[table omitted — 3,971 chars of source]

It is also of interest to compare the semiparametric estimator with the TSLS estimator. Directly comparing the the two estimators can be difficult since they are of different dimensions and converge to different limits. Instead, we can compare their performance on estimating the marginal effect. For the structural function $g_\gamma(t,\bar{r},e) = \gamma_1 t + \gamma_2 t^2 + \gamma_3 te + e$, the marginal effect of the treatment on the outcome is $\frac{\partial}{\partial t} g_\gamma(t,\bar{r},e) = \gamma_1 + 2\gamma_2 t + \gamma_3 e$, which takes on different values for different treatment and outcome levels. For a given treatment level $t$, we can use the semiparametric estimator to obtain an estimate $\hat{\gamma}_1 + 2\hat{\gamma}_2 t + \hat{\gamma}_3 e$ of the marginal effect. However, a TSLS procedure would deliver a scalar estimate that is a mixture of marginal effects across different treatment and outcome levels. Figure (ref) shows that with a nonlinear specification the semiparametric estimator outperforms the TSLS estimator. Figure (ref) shows similar findings with a fully nonlinear and nonseparable specification.

Next, we compare the semiparametric estimator with the TSLS estimator when the structural function $g$ is linear. This is achieved by imposing $\gamma_2 = \gamma_3 = 0$. In this case, the TSLS estimator is consistent for the coefficient $\gamma_1$. However, the identification of the TSLS estimator is based solely on the difference between the two means. If the two distributions corresponding to $h_0$ and $h_1$ have the same mean, then the TSLS procedure suffers from weak identification issues. In contrast, the identification of the semiparametric estimator $\hat{\gamma}$ is based on the entire difference between $h_0$ and $h_1$. The semiparametric estimator continues to work even if the estimand of the TSLS estimator is weakly identified. For the simulation, we let $h_0$ be the quantile function of Beta$(0.1,0.1)$, and $h_1$ be the quantile function of Beta$(10,10)$. These two distributions are significantly different, but they have the same mean (0.5). Figure \ref*{fig:linear-sf} shows that in this case the semiparametric estimator outperforms the TSLS estimator even if the structural function is linear.

figure[figure omitted — 917 chars of source]

figure[figure omitted — 927 chars of source]

figure[figure omitted — 739 chars of source]

Conclusion

In this study, we have examined the identification and estimation of the structural function in an RD design with a continuous treatment variable. We have established the nonparametric identification result and proposed a semiparametric estimator for the possibly nonlinear and nonseparable structural function. The estimator is proven to be consistent and asymptotically normal. The empirical application and simulation studies demonstrate the advantage of the semiparametric estimator compared to the TSLS estimator.

There are two promising ways to extend the results in this paper in the future. First, we can consider extrapolating the identification result away from the cutoff. This can be done by identifying the derivative of the structural function with respect to the running variable at the cutoff as in dong2015identifying. Second, we can apply the methodology developed in this paper to the regression kink design model studied by card2015RKD,dong2018jump where the treatment choice function exhibits a kink instead of a discontinuity at the cutoff.

Extrapolation Away from the Cutoff

Thus far, we have been focusing on the structural function local to the cutoff. Next, we want to study how to extrapolate away from the cutoff. This is a difficult task since the exogenous variation in the treatment only occurs at the cutoff.

assumption[Local Independence] There exists $\delta >0$, such that $(\varepsilon,U)$ and $R$ are independent when $R$ is in the $\delta$-neighborhood of $\bar{r}$, i.e., \begin{align*} (\varepsilon,U) \perp R | R \in [\bar{r}-\delta,\bar{r}+\delta]. \end{align*}

The above local independence assumption is standard in the RD literature. For example, it is imposed in the paper by hahn2001identification.

theoremLet Assumptions (ref) - (ref) and (ref) hold. For any $g \in \mathcal{G}$ that satisfies Condition ((ref)), let $F^g_{\varepsilon | U,R}(\cdot | \cdot,\bar{r})$ be defined as in ((ref)) and $\lambda^g$ the transformation in Theorem (ref). Then for $t \in [t_0',t_0''] \cup [t_1',t_1''], r \in [\bar{r} - \delta,\bar{r}) \cup (\bar{r},\bar{r}+\delta] $ and $e \in \mathcal{E}$, \begin{align*} g^*(t,r,e) = F_{Y | T,R}^{-1}\big( F^g_{\varepsilon | U,R}(\lambda^g(e) | h^{-1}(r,t),\bar{r} ) | t,r \big) , \end{align*} where $F_{Y | T,R}^{-1}$ represents the inverse of $F_{Y | T,R}$ with respect to the first argument and is well-defined in view of Assumption (ref) and (ref). Furthermore, the CASF is identified by \begin{align*} \beta^*(t,r) = \int_0^1 \int_\mathcal{E} F_{Y | T,R}^{-1}\big( F^g_{\varepsilon | U,R}(e | h^{-1}(r,t),\bar{r} ) | t,r \big) dF^g_{\varepsilon | U,R}(e | h^{-1}(r,t),\bar{r} )du . \end{align*}

proof[Proof of Theorem (ref)] By Corollary (ref) and Assumption (ref), we have for $r$ in the $\delta$-neighborhood of $\bar{r}$, \begin{align*} F^g_{\varepsilon | U,R}(\lambda^g(e) | h^{-1}(r,t),\bar{r} ) & = F_{\varepsilon | U,R}(e | h^{-1}(r,t),\bar{r} ) = F_{\varepsilon | U,R}(e | h^{-1}(r,t),r ) \\ & = \mathbb{P}( \varepsilon \leq e | T = t, R = r ) = \mathbb{P}( Y \leq g^*(t,r,e) | T = t, R = r ) \\ & = F_{Y | T,R}(g^*(t,r,e) | t,r). \end{align*} The statement regarding the CASF follows from change of variables ($e$ to $\lambda^g(e)$).

The preceding theorem identifies the structural function as well as the CASF in the local neighborhood of $\bar{r}$ where the independence condition specified by Assumption (ref) holds. We can use Taylor's approximation to further extrapolate beyond this local neighborhood. For example, if $\beta^*$ is differentiable, then we can identify its derivatives at $\bar{r}$ by using Theorem (ref) and employ the following approximation:

align*[align* omitted — 161 chars of source]

Notice that the validity of the approximation above is based on the smoothness of $\beta^*$ and does not require the evaluation point $r$ to be within the $\delta$-neighborhood of $\bar{r}$. When the CASF $\beta^*$ is assumed to have higher-order derivatives, we can employ higher-order Taylor's approximations. This ability to identify derivatives (of any order) is, in fact, the addtional identification power of the independence condition (Assumption (ref)) as compared to the smoothness condition (Assumption (ref)(i)). Such a difference has been studied by dong2018alternative in the binary treatment case.

Regression Kink Designs

In the RK designs, the model ((ref)) - ((ref)) remains the same. However, there is a kink rather than a discontinuity in the first stage. That is,

align*[align* omitted — 178 chars of source]

Assumption (ref)(ii), the strong discontinuity condition, no longer holds in this case since the first stage is continuous in the running variable. Therefore, we cannot apply the previous identification results to the RK design.

To apply our identification results, we can reformulate the treatment variable and the first stage. Define

align[align omitted — 98 chars of source]

as the marginal effect of the treatment. Suppose the function $\frac{\partial}{\partial R} h(R,\cdot)$ is strictly increasing, then we can establish one-to-one correspondence between $(\tilde{T},R)$ and $(U,R)$ and further between $(\tilde{T},R)$ and $(T,R)$. Then by treating $\tilde{T}$ as the treatment variable and ((ref)) the first stage, we are back to the RD setup and can hence employ the previous identification results.

\fi