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.
82,711 characters · 7 sections · 51 citation commands
Estimating a Continuous Treatment Model with Spillovers: A Control Function Approach
Keywords: Continuous treatment, Endogeneity, Instrumental variables, Treatment spillover.
Treatment evaluation under cross-unit interference is currently one of the most intensively studied topics in the causal inference literature (see, e.g., halloran2016dependent, aronow2021spillover for recent reviews). These previous studies have highlighted the importance of accounting for treatment spillovers from other units via empirical applications in many fields, including political science, epidemiology, education, economics, etc. However, most of these studies have focused on the case of simple binary treatments, and estimating treatment spillover effects with continuous treatments has been rarely considered. Even in the context of binary treatments, they often assume that the treatments are (conditionally) randomly assigned with full compliance (with some exceptions including kang2016peer, imai2020causal, ditraglia2021identifying, hoshino2021causal, vazquez2021causal). To fill these gaps, this paper considers the estimation of a continuous treatment model with treatment spillovers while allowing the subjects to self-select their own treatment levels.
For conventional continuous treatment models without spillovers, several approaches have been proposed to estimate the average dose-response function -- the expected value of a potential outcome at particular treatment level. Among others, hirano2004propensity and imai2004causal propose the generalized propensity score (GPS) approach under an unconfoundedness assumption. A further extension of the GPS approach to the case of multiple correlated treatments is discussed in egger2013generalized. Then, it would be a natural idea to apply these methods to estimate continuous treatment and spillover effects jointly, as in del2020causal, in which the authors use a GPS method to adjust for confounding variables for both own and others' treatments.
However, the above-mentioned methods are not successful in the presence of unobservable confounding variables, which might lead to endogeneity issues not only with respect to own treatment but also to others' treatments. For example, suppose the outcome is academic performance for high school students, and the treatment of interest is hours spent on in-home private tutoring. Since close school mates are often taught by the same teachers and have similar socioeconomic backgrounds, they share a variety of unobservable attributes that affect the treatment and potential outcomes simultaneously. In this situation, not only one's own treatment level but also the treatments of his/her friends should be treated as endogenous variables. Other examples include the causal effect of minimum wage on the local employment rate, the effect of police force on the crime rate, and the effect of corporate tax subsidies on regional revenue. These are cases where “spatial” treatment spillovers are of concern and are typical empirical situations that can fit into the framework of this paper.
To deal with this complex endogeneity problem without relying on strong parametric distributional restrictions, we assume a specific type of triangular model in this study, in which the treatment equation has a certain form of separability and others' treatments influence one's outcome in the form of peers' weighted average treatment, possibly with some monotonic transformation. The latter assumption is similar to and more general than the mean interaction in manski2013identification, which assumes that the impacts from others can be summarized as the empirical mean of peers' treatments. Since causal inference is generally impossible if no assumptions are imposed on the interference structure (cf. imbens2015causal), assumptions similar to the above are widely adopted in the literature.
With these assumptions, we address the endogeneity by employing a control function approach (see, e.g., blundell2003endogeneity, florens2008identification, imbens2009identification), which introduces auxiliary regressors, the control variables, in the estimation of the outcome equation to eliminate the endogeneity bias. Under the availability of valid instrumental variables (IVs), the rank variable in the treatment equation can serve as the control variable for own treatment variable, which is a standard result in the literature. Meanwhile, to control the endogeneity for the treatment spillovers, a straightforward choice for the control variables would be to use the peers' rank variables by analogy. Although this approach is theoretically simple, it may not be practical unless the number of interacting partners is limited to one or two because of the curse of dimensionality. As a novel finding of this paper, we demonstrate that the rank variable of a “hypothetical” individual who is on average equivalent to real peers can be used as a valid control variable by utilizing the additive structure and interaction structure of our model. Since this rank variable is one dimensional, we can alleviate the dimensionality problem.
As is well-known, to achieve nonparametric identification of treatment parameters such as the average structural function based on a control function approach, we often require a strong support condition, such that the support of the control variable conditional on the treatment variable is equal to the marginal support of the control variable (see imbens2009identification). This is true for our model as well, but such condition may be rarely satisfied in reality. Thus, to improve the empirical tractability, we introduce additional functional form restrictions for estimation, similar to chernozhukov2020semiparametric and newey2021control. The resulting model takes the form of a multiplicative functional-coefficient regression model, which can be seen as a special case of the so-called varying-coefficient additive models (see, e.g., zhang2015varying,hu2019estimation). For this model, as a main parameter of interest, we focus on the estimation of the conditional average treatment response (CATR) -- the expectation of the potential outcome conditional on individual covariates.
For the estimation of the CATR parameter, we preliminarily need to estimate two control variables, one for own treatment effect and the other for the treatment spillover effect. This preliminary step can be easily implemented with a composite quantile regression (CQR) method (see, e.g., zou2008composite). Once the control variables are obtained, the CATR is estimated in three steps. In the first step, we globally estimate the model without considering the multiplicative structure using a penalized series regression approach. Next, we fine-tune the bias-correction function in the model by marginal integration to improve the efficiency. Given this result, we finally re-estimate the model using a local linear kernel regression. It is well-known that, to achieve pointwise asymptotic normality for a series estimator, we are often required to inefficiently undersmooth the estimator (cf. huang2003local). For this reason, the idea of first obtaining a globally consistent estimate using a series method and then locally re-adjusting the estimate using a kernel method is widely adopted in the literature (see, e.g., horowitz2004nonparametric, horowitz2005nonparametric, wang2007spline). With this multiple-step estimation procedure, we can show that the final CATR estimator is oracle, in that its asymptotic distribution is not affected by the estimation of the other nuisance functions.
To demonstrate the empirical usefulness of the proposed method, we conduct an empirical analysis on the causal effects of unemployment on crime. In this empirical analysis, using Japanese city-level data, we estimate the CATR parameter, with the crime rate as the outcome of interest and each city's own unemployment rate and the average unemployment rate of its neighboring cities as the treatments. To account for the endogeneity, we use the availability of childcare facilities as an instrument for the unemployment rate. As a result, we find that the CATR tends to increase as the unemployment rate increases, which is a common finding in the literature. Moreover, our model reveals new empirical evidence that the spillover effect of unemployment from the neighboring cities is important only for non-rural cities.
To summarize, the main contributions of this study are fourfold. First, to our knowledge, this study is the first to address a continuous treatment model in which both treatment spillover and endogeneity are present. Second, we propose a novel control function approach to establish the nonparametric identification of this model under some functional form restrictions. Third, considering the empirical feasibility, we propose a semiparametric multiplicative potential outcome model and develop a three-step estimation procedure for it, which is of independent interest in semiparametric estimation theory. Fourth, by applying the proposed method to Japanese city data, we provide new empirical evidence on the relationship between the local unemployment and crime rates.
\paragraph{Paper organization:} In Section (ref), we present our model and discuss nonparametric identification of the CATR parameter. In Section (ref), we introduce an empirically tractable semiparametric model and propose our three-step estimation procedure. The asymptotic properties of the estimator are also presented in this section. The empirical analysis on the Japanese crime data is presented in Section (ref), and Section (ref) concludes. The proofs of technical results, Monte Carlo simulation results, and supplementary information on the empirical analysis are all summarized in the online supplementary material.
\paragraph{Notation:} For natural numbers $a$ and $b$, $I_a$ denotes an $a \times a$ identity matrix, and $\mathbf{0}_{a \times b}$ denotes a matrix of zeros of dimension $a \times b$. For a matrix $A$, $|| A ||$ denotes the Frobenius norm. If $A$ is a square matrix, we use $\rho_{\max} (A)$ and $\rho_{\min} (A)$ to denote its largest and smallest eigenvalues, respectively. We use $\otimes$ and $\circ$ to denote the Kronecker and Hadamard (element-wise) product, respectively. For a set $\mathcal{A}$, $|\mathcal{A}|$ denotes its cardinality. We use $c$ (possibly with a subscript) to denote a generic positive constant whose value may vary in different contexts. For random variables $X$ and $Y$, $X \mathrel{\text{\scalebox{1.07}{$\perp\mkern-10mu\perp$}}} Y$ means that they are independent. Lastly, for positive sequences $a_n$ and $b_n$, $a_n \asymp b_n$ means that there exist $c_1$ and $c_2$ such that $0 < c_1 \le a_n/b_n \le c_2 < \infty$ for sufficiently large $n$.
Suppose that we have a sample of $N$ agents that form social networks whose connections are represented by an $N \times N$ adjacency matrix $\bm{A}_N = (A_{i,j})_{1 \le i,j \le N}$. These agents can be individuals, firms, or municipalities depending on the context. The networks can be directed, that is, regardless of the value of $A_{j,i}$, we may observe $A_{i,j} = 1$ if $j$ affects $i$ and $A_{i,j} = 0$ otherwise. The diagonal elements of $\bm{A}_N$ are all zero. Throughout the paper, we treat $\bm{A}_N$ as non-random. For each $i$, we denote $i$'s “reference group” (peers, colleagues, neighbors, etc) as $\mathcal{P}_i \coloneqq \{1 \le j \le N: A_{i,j} = 1\}$ and its size as $n_i = |\mathcal{P}_i|$. To simplify the discussion, we assume that $n_i > 0$ for all $i$. In addition, we write $\overline{\mathcal{P}}_i \coloneqq i \cup \mathcal{P}_i$. For a general variable $Q$, we denote $\bm{Q}_{\mathcal{P}_i} \coloneqq \{Q_j\}_{j \in \mathcal{P}_i}$. We define $\bm{Q}_{\overline{\mathcal{P}}_i}$ similarly.
Let $T \in \mathcal{T}$ denote the continuous treatment variable of interest. We assume that the support $\mathcal{T}$ is common across all individuals and is a real closed interval of $\mathbb{R}$. The outcome variable of interest is $Y \in \mathbb{R}$. In this study, we explicitly allow that the peers' treatments $\bm{T}_{\mathcal{P}_i}$ influence on $i$'s own outcome $Y_i$ via a known real-valued function $S_i$:
where $m_i$ is a strictly increasing continuous function which may depend on $i$, and $a_{i,j}$'s are weight terms satisfying $a_{i,j} > 0$ if and only if $j \in \mathcal{P}_i$ and $\sum_{j \in \mathcal{P}_i} a_{i,j} = 1$.\footnote{ Assuming that the functional form of $S_i$ is known a priori is a strong assumption. In the literature, there are several studies that investigate under what conditions we can make meaningful causal inference even when the “exposure” function (i.e., $S_i$) is unknown or mis-specified (e.g., hoshino2021causal, savje2021average, leung2022causal, VAZQUEZBARE2022). All these studies consider only binary treatment situations, and whether we can establish similar results for continuous treatment models would be an important open question. } For example, if $m_i$ is an identity function and $a_{i,j} = n_i^{-1}$, then $S_i$ simply returns the reference-group average. This form of treatment spillover is the most commonly used in empirical studies on peer effects. If $m_i(t) = n_i t$ instead, then $S_i$ in this case is the sum of peers' treatments. For another example, in the spatial statistics literature, researchers often assume that $a_{i,j}$ is inversely proportional to the geographical distance between $i$ and $j$. The support of $S_i$ is denoted as $\mathcal{S}_i$, which varies across individuals in general because of the transformation $m_i$ that may be specific to $i$. As a special case, if $m_i$ is an identity function, then we have the same support for all $i$ as $\mathcal{S}_i = \mathcal{T}$.
Now, letting $X_i = (X_{1,i}, \ldots, X_{dx, i})^\top$ be a $dx \times 1$ vector of observed individual covariates, we suppose that the outcome $Y_i$ is determined in the following equation:
where $y$ is an unknown function, $S_i = S_i(\bm{T}_{\mathcal{P}_i})$, and $\epsilon_i$ is an unobservable determinant of $Y_i$. Note that $\epsilon_i$ could be a vector of unknown dimension. The non-separable setting is essential for a treatment effect model since it permits general interactions between $(T_i, S_i)$ and $\epsilon_i$, allowing for treatment heterogeneity among observationally identical individuals. The potential outcome when $(T_i, S_i) = (t, s) \in \mathcal{T}\mathcal{S}_i$ is written as $Y_i(t, s) = y(t, s, X_i, \epsilon_i)$. The target parameter of interest in this study is the conditional average treatment response (CATR) function:\footnote{ Following the econometrics literature (e.g., blundell2003endogeneity), this parameter may be called the average structural function (conditional on $X_i = x$). In the causal inference literature, it is also known as the average dose response function. }
Next, suppose that we have a $dz \times 1$ vector of IVs, $Z_i = (Z_{1,i}, \ldots, Z_{dz,i})^\top$, which determines the value of $T_i$ through the following equation:
where $\pi$ and $\eta$ are unknown functions, and $U_i$ is a scalar unobservable variable. We allow that not only $U_i$ but also $\bm{U}_{\mathcal{P}_i}$ are potentially correlated with $\epsilon_i$, which are the sources of endogeneity for $T_i$ and $S_i$, respectively. To address the endogeneity issue, we make the following assumptions:
Assumption (ref)(i) is the exogeneity assumption that own and peers' IVs are independent of the unobservables given their covariates. This assumption would be particularly plausible in an experimental setup where $Z_i$ is a randomly assigned treatment inducement. Note that we do not preclude the cross-sectional correlation of the covariates, IVs, and unobservables, which should be an important feature of network data. For Assumption (ref)(ii), we emphasize that the functional form of $\eta$ is independent of $\bm{X}_{\overline{\mathcal{P}}_i}$. Note that although this independence restriction is not strictly necessary for a consistent estimation of the treatment equation (ref),\footnote{ When considering a more general treatment equation $T_i = \pi(X_i, Z_i) + \eta(U_i, \bm{X}_{\overline{\mathcal{P}}_i})$, the second part of Assumption (ref)(ii) is always satisfied using $\widetilde U_i$ and $\widetilde \eta$, where $\widetilde U_i = F_{U_i|\bm{X}_{\overline{\mathcal{P}}_i}}(U_i \mid \bm{X}_{\overline{\mathcal{P}}_i})$ and $\widetilde \eta(\widetilde U_i, \bm{X}_{\overline{\mathcal{P}}_i}) = \eta(F_{U_i|\bm{X}_{\overline{\mathcal{P}}_i}}^{-1}(\widetilde U_i \mid \bm{X}_{\overline{\mathcal{P}}_i}))$. Under Assumption (ref)(i), we have $\Pr(T_i \le \pi(X_i, Z_i) + \eta(u, \bm{X}_{\overline{\mathcal{P}}_i}) \mid Z_i, \bm{X}_{\overline{\mathcal{P}}_i}) = \Pr(U_i \le u \mid \bm{X}_{\overline{\mathcal{P}}_i}) = u$ for any $u \in (0,1)$, which implies that $\pi(X_i, Z_i) + \eta(u, \bm{X}_{\overline{\mathcal{P}}_i})$ can be estimated in a quantile regression framework. Furthermore, by checking whether the estimated conditional quantile function varies with the peers' covariates, we may be able to test the validity of assuming $\eta(U_i, \bm{X}_{\overline{\mathcal{P}}_i}) = \eta(U_i)$; however, directly implementing this would suffer from the curse of dimensionality. } we introduce it for constructing an effective control variable for the spillover effect.
Recall that the disturbance term in the outcome equation is allowed to be a vector and to enter the model in a non-additive way. Thus, we can see that the identification result heavily depends on the functional form of the treatment equation rather than that of the outcome equation, in line with the existing literature.
Below, we formally discuss the identification of the CATR parameter. Throughout this section, we use the term “identification” to indicate that the parameter of interest can be characterized through a moment of the observable random variables. Note that this does not automatically imply the estimability of the parameters because network data naturally follow a nonidentical and dependent data distribution. \footnote{ In this sense, the CATR parameter should be indexed by “$\:i\:$,” however, we suppress it for notational simplicity. } To estimate the CATR parameter, as discussed later, we require some additional stationarity assumptions.
As is well-known, we can deal with the endogeneity of $T_i$ by including $U_i$ as a control variable in the estimation of outcome equation (ref) (e.g., blundell2003endogeneity, florens2008identification, imbens2009identification). For the identification of $U_i$, in Lemma (ref), we prove that $\pi_i(X_i,Z_i)$ and $\eta(u)$ for any $u \in (0,1)$ are identifiable up to a location shift (see Remark (ref) as well). When these parameters are treated as known, $U_i$ can be obtained by solving $\min_{u \in [0,1]}|T_i - \pi(X_i,Z_i) - \eta(u)|$ under the monotonicity and continuity of $\eta$.
Similarly, to account for the endogeneity of $S_i$, we need to find a control variable(s) for it. One obviously valid candidate is to use $\bm{U}_{\mathcal{P}_i}$. Once the behaviors of $\bm{U}_{\mathcal{P}_i}$ are controlled, the variation of $S_i$ comes only from $(\bm{X}_{\mathcal{P}_i}, \bm{Z}_{\mathcal{P}_i})$, which are conditionally independent of the disturbances by Assumption (ref)(i). However, this approach is not practical unless $n_i$ is limited to one or two, because of the curse of dimensionality. In the next proposition, we show that under the additive separability in (ref) we can construct a one-dimensional control function for $S_i$.
Proposition (ref) implies that, for a general function $f$ of $\epsilon$, we have
where $c_1 = t- \eta(u)$ and $c_2 = m_i^{-1}(s) - \eta(v)$, by Assumption (ref)(i). The identification of $V_i$ can be achieved by solving $\min_{v \in [0,1]}|\sum_{j \in \mathcal{P}_i} a_{i,j} (T_j - \pi(X_j, Z_j)) - \eta(v)|$. Without the additive separability between $(X_i,Z_i)$ and $U_i$ as in (ref), $V_i$ generally depends on the instruments $\bm{Z}_{\mathcal{P}_i}$, and conditioning on $V_i$ also restricts their behavior (cf. kasy2011identification), leading to a failure of establishing (ref).
Now, we define the marginal treatment response (MTR) function:
From (ref), we have
Hence, we can identify $\mathrm{MTR}(t, s, u, v, \bm{x})$ as the left-hand side of (ref). Moreover, if $\mathrm{MTR}(t, s, u, v, \bm{x})$ is identifiable uniformly over $\mathrm{supp}( U_i, V_i, \bm{X}_{\mathcal{P}_i} \mid X_i = x)$, we can identify $\mathrm{CATR}(t, s, x)$ by
where the outer expectation is with respect to $\bm{X}_{\mathcal{P}_i}$ conditional on $X_i = x$, and $f_{U_i V_i|\bm{X}_{\overline{\mathcal{P}}_i}}$ is the conditional density function of $(U_i, V_i)$ given $\bm{X}_{\overline{\mathcal{P}}_i}$. For $\mathrm{MTR}(t, s, u, v, \bm{x})$ to be well-defined on the entire $\mathrm{supp}( U_i, V_i, \bm{X}_{\mathcal{P}_i} \mid X_i = x)$, we need the following additional condition:
From (ref), we can see that without any support restriction, $\mathrm{MTR}(t, s, u, v, \bm{x})$ can be obtained only on $\mathrm{supp}(U_i, V_i, \bm{X}_{\mathcal{P}_i} \mid T_i = t, S_i = s, X_i = x)$. Assumption (ref) requires that conditioning on $\{T_i = t, S_i = s\}$ should not affect the support of $(U_i, V_i, \bm{X}_{\mathcal{P}_i})$ conditional on $X_i = x$. An important implication of this condition is that the IV must have very rich variation so that the values of $T_i$ and $S_i$ do not restrict the range of values that $U_i$ and $V_i$ can take. We summarize the result obtained so far in the next theorem.
In order to construct an estimator for $\mathrm{CATR}(t, s, x)$ based on Theorem (ref) in practice, we need to cope with several obstacles. The first issue is the full-support condition in Assumption (ref). Finding IVs ensuring such a support condition would be quite difficult in reality. Even without this support condition, if there are informative upper and lower bounds of $Y_i$, it would be possible to partially identify $\mathrm{CATR}(t, s, x)$, as in Theorem 4 of imbens2009identification. For another more convenient approach to point-identify the CATR in the absence of Assumption (ref), in the next section and thereafter, we consider introducing additional functional form restrictions on the outcome equation.
Another issue is that, since the size of the reference group and the weights $a_{\mathcal{P}_i}$ are not necessarily identical among individuals, the support $\mathcal{S}_i$ of $S_i$ and the distribution of $V_i$ should vary with individuals in general. The dimension and distribution of $\bm{X}_{\mathcal{P}_i}$ are also obviously heterogeneous among individuals. These heterogeneities are problematic in the estimation stage. To circumvent this problem, we simply restrict our attention to a specific subsample $\mathcal{N}'$ such that those in this subsample have the same support $\mathcal{S}$ and are homogenous (in the sense of Assumption (ref)(ii) below). The simplest but empirically the most typical example would be $\mathcal{N}' = \{1 \le i \le N: m_i = m, \: n_i = c, \: a_{i,j} = c^{-1}\mathbf{1}\{j \in \mathcal{P}_i\}\}$ for some monotonically increasing $m$ and an integer $c \ge 1$. The size of this subsample is denoted as $n = |\mathcal{N}'|$. In the following, without loss of generality, we re-label the $N$ agents so that the first $n$ individuals belong to this subsample, that is, $\mathcal{N}' = \{1, \ldots, n\}$.
In this section, we propose our estimation procedure for $\mathrm{CATR}(t, s, x)$ in a specific semiparametric model, and investigate its asymptotic properties as $n$ increases to infinity. To increase $n$, we consider a sequence of networks $\{ \bm{A}_N \}$, where $\bm{A}_N$ is an $N \times N$ adjacency matrix. Since each $\bm{A}_N$ is assumed to be non-stochastic, our analysis can be interpreted as being conditional on the realization of $\bm{A}_N$. In this sense, the model and parameters presented in this study should be viewed as triangular arrays defined along the network sequence.
In newey2021control, they proposed the following type of potential outcome model as a baseline (in our notation): $Y_i(t,s) = g(t, s)^\top \epsilon_i$, where $g(t, s)$ is a vector of transformations of $(t,s)$; for example, $g(t,s) = (1,t,s)^\top$ (note that they do not consider models with two treatment variables). A similar parametric assumption was adopted in chernozhukov2020semiparametric as well. In this study, we generalize their models such that $g(t,s)$ is a general nonparametric function. In addition, we further extend their work by allowing the treatment effect to vary with the individual characteristics $X_i$. At the same time, to preserve empirical tractability, we assume that there are no interaction effects among $(X_i, \epsilon_i)$, and that $\epsilon_i$ is one-dimensional. Consequently, we focus on the following outcome model:
where $\beta(t,s) = (\beta_1(t,s), \ldots, \beta_{dx}(t,s))^\top$ and $g(t,s)$ are unknown functions to be estimated.\footnote{ If we add one more continuous treatment, then $\beta$'s and $g$ will become three-dimensional functions, which are practically difficult to estimate because of the curse of dimensionality. To estimate such a model, we would need to introduce additional functional form restrictions on (ref). } For normalization, we assume that $X_i$ includes a constant term so that $\operatorname*{\mathbb{E}}[\epsilon_i] = 0$ holds. We do not additionally restrict the treatment equation (ref), but we assume that the control variables $U_i$ and $V_i$ can be consistently estimated at a certain convergence rate (see Remark (ref) below). For further simplification, we strengthen Assumption (ref)(i) as follows:
\setcounter{assumption}{-1}
Recalling the definition of $V_i$, the first part of Assumption (ref)(ii) essentially requires that the joint distribution of $\bm{U}_{\overline{\mathcal{P}}_i}$ is stationary across all $i \in \mathcal{N}'$ and they adopt the same weighting scheme. For the second part, we do not require the treatment heterogeneity terms to be identically distributed. Since Assumption (ref) is fundamental to derive our estimation procedure, it will be assumed implicitly throughout the whole subsequent discussion.
Now, the CATR for this model at $(t,s,x)$ is simply given by
which in fact coincides with the “unconditional” expectation $\operatorname*{\mathbb{E}}[y(t,s,x,\epsilon_i)]$. Thus, the task of estimating $\textrm{CATR}(t,s,x)$ is greatly simplified to the estimation of $\beta$, and the support condition in Assumption (ref) is not necessary to recover the CATR. Further, an analogous argument to (ref) gives
This yields the following semiparametric multiplicative regression model:
By construction, we have $\operatorname*{\mathbb{E}}[\widetilde\epsilon_i \mid T_i, S_i, U_i, V_i, X_i] = 0$. It is clear that because of the multiplicative structure in (ref), in order to identify $g$ and $\lambda$ separately, we need some functional form normalizations. First, recall that the condition $\operatorname*{\mathbb{E}}[\epsilon_i] = 0$ must be maintained. This implies the following location normalization:
We further need some scale normalization for identification. To this end, we set
Regression models with a multiplicative structure as in (ref) can be found in the literature in different contexts (e.g., linton1995kernel, zhang2015varying, hu2019estimation, chen2020estimation). A typical approach to estimating this type of model is to apply some marginal integration (to the estimated nonparametric function or to the data itself in advance of the estimation). However, note that none of these prior studies considered a functional-coefficient specification as in ours. Thus, we need to develop a new estimation procedure for our model, which should be of independent interest in the semiparametric estimation literature. Specifically, we deal with the multiplicative structure by nonparametrically estimating the model ignoring the functional form restriction in the first step, and then splitting the estimated function into two multiplicative components by marginal integration in the second step in a similar manner to chen2020estimation.\footnote{ Although it is possible to estimate the model in a single step by estimating $g$ and $\lambda$ separately from the beginning, the resulting method requires solving a high-dimensional non-convex optimization problem, which is computationally challenging, as pointed out in zhang2020new. To circumvent this issue, zhang2020new proposed an iterative computational algorithm. } The first stage estimation involves a four-dimensional nonparametric regression on $(T_i,S_i, U_i, V_i)$, which often produces unstable estimates under a moderate sample size. Thus, we consider using a penalized regression in this step.
The whole estimation procedure is as follows. Before estimating the regression model in (ref), we need to consistently estimate the function $\pi$ and $\eta$ to obtain consistent estimates of $(U_i, V_i)$. Excluding this preliminary step, our procedure for estimating the CATR consists of three steps. In the first step, we estimate the model in (ref) ignoring the multiplicative structure using a penalized series regression approach. In the next step, we obtain estimates of $g$ and $\lambda$ by marginally integrating the nonparametric function obtained in the previous step. Finally, we re-estimate the coefficient functions $\beta$ using a local linear kernel regression. With these two additional steps, we can achieve an oracle property for the estimation of the CATR parameter (see Theorem (ref)). The proposed estimation procedure is new to the literature but may be viewed as a minor extension of a commonly used two-step series-kernel estimation method, in which chen2020estimation's marginal integration step is inserted between the series and kernel estimation steps.
\paragraph{Preliminary step: control variables estimation}
Under Assumptions (ref)(ii) and (ref)(i), it holds that $\Pr(T_i \le \pi(X_i, Z_i) + \eta(u) \mid X_i, Z_i) = u$ for any $u \in (0,1)$. This implies that we can estimate $\pi(X_i, Z_i)$ and $\eta(u)$ using the CQR method: for pre-specified $0 < u_1 < u_2 < \cdots < u_L < 1$,
subject to some location normalization (e.g., $\eta(0.5) = 0$), where $q_u(x) \coloneqq x(u - \mathbf{1}\{x < 0\})$. Then, we compute the residual, $res_i \coloneqq T_i - \widehat \pi_n(X_i, Z_i)$, for each $i$. Note that $res_i$ is an estimator of $\eta(U_i)$, whose convergence rate is governed by that for $\widehat \pi_n$. The estimator for $U_i$ can be obtained by $\widehat U_i \coloneqq \operatorname*{argmin}_{u \in [0,1]}|res_i - \widehat \eta_n(u)|$. Similarly, we can estimate $V_i$ by $\widehat V_i \coloneqq \operatorname*{argmin}_{v \in [0,1]}| \sum_{j \in \mathcal{P}_i} a_{i,j} res_j - \widehat \eta_n(v)|$. In practice, these minimization problems are solved by grid search with a sufficiently large $L$.
\paragraph{First step: penalized series estimation}
Let $\{(p_{T,j}(t), p_{S,j}(s), p_{U,j}(u), p_{V,j}(v)): j = 1, 2, \ldots \}$ be the basis functions on $\mathcal{T}$, $\mathcal{S}$, $[0,1]$, and $[0,1]$, respectively. We define $P_T(t) \coloneqq (p_{T,1}(t), \ldots , p_{T,K_T}(t))^\top$, $P_S(s) \coloneqq (p_{S,1}(s), \ldots , p_{S,K_S}(s))^\top$, and $P_{TS}(t, s) \coloneqq P_T(t) \otimes P_S(s)$. Then, we consider series approximating $\beta_l(\cdot, \cdot)$ by $\beta_l(t, s) \approx P_{TS}(t, s)^\top \theta_{\beta_l}$ for some $K_{TS} \times 1$ coefficient vector $\theta_{\beta_l}$, for each $l = 1, \ldots, dx$, where $K_{TS} \coloneqq K_T K_S$. Similarly, define $P_U(u) \coloneqq (p_{U,1}(u), \ldots , p_{U,K_U}(u))^\top$, $P_V(v) \coloneqq (p_{V,1}(v), \ldots , p_{V,K_V}(v))^\top$, $P_{UV}(u, v) \coloneqq P_U(u) \otimes P_V(v)$, and $K_{UV} \coloneqq K_U K_V$. We assume that $K_{TS}$ and $K_{UV}$ increase as $n$ increases at the same speed so that there exists an increasing sequence $\{\kappa_n\}$ satisfying $(K_{TS}, K_{UV}) \asymp \kappa_n$.
For the estimation of $g(t,s)\lambda(u,v)$, to impose the location condition (ref) in the estimation, we would like to normalize the basis function such that $\overline{P}_{UV}(u, v) \coloneqq P_{UV}(u, v) - \operatorname*{\mathbb{E}}[P_{UV}(U, V)]$ (recall that we have been focusing on a subsample in which $\{(U_i, V_i)\}$ are identically distributed). However, since both $\operatorname*{\mathbb{E}}[P_{UV}(U, V)]$ and $(U_i, V_i)$ are unknown, we instead use
Then, letting $\widehat P(t, s, u, v) \coloneqq P_{TS}(t, s) \otimes \widehat P_{UV}(u, v)$, we consider approximating $g(t,s)\lambda(u,v) \approx \widehat P(t, s, u, v)^\top \theta_{g\lambda}$ with a $K_{TSUV} \times 1$ coefficient vector $\theta_{g\lambda}$, where $K_{TSUV} \coloneqq K_{TS} K_{UV}$. Using these approximations, we have
where $P_X(X_i, T_i, S_i) \coloneqq X_i \otimes P_{TS}(T_i, S_i)$, and $\theta_\beta = (\theta_{\beta_1}^\top, \ldots, \theta_{\beta_{dx}}^\top)^\top$. Based on this model approximation, we estimate $\theta = (\theta_\beta^\top, \theta_{g\lambda}^\top)^\top$ by the penalized least squares method.
Let $\bm{D}$ be a $(dx K_{TS} + K_{TSUV})$-dimensional positive semidefinite symmetric matrix such that $\rho_{\max}(\bm{D}) = O(1)$. This matrix serves as a penalization matrix, in which typical elements of $\bm{D}$ are, for example, the integrated derivatives of the basis functions. Write $\widehat \Pi_i \coloneqq (P_X(X_i, T_i, S_i)^\top, \widehat P(T_i, S_i, \widehat U_i, \widehat V_i)^\top)^\top$, $\widehat{\bm{\Pi}}_n = (\widehat \Pi_1, \ldots, \widehat \Pi_n)^\top$, and $\bm{Y}_n = (Y_1, \ldots, Y_n)^\top$. Then, for a sequence of tuning parameters $\{\tau_n\}$ tending to zero at a certain rate as $n \to \infty$, the penalized least squares estimator of $\theta$ is given by
The first-stage estimator of $\beta_l(t, s)$ can be obtained by $\widehat \beta_{n,l}(t, s) \coloneqq P_{TS}(t, s)^\top \widehat \theta_{n, \beta_l}$ for $l = 1 , \ldots, dx$. Similarly, $g(t,s) \lambda(u, v)$ can be estimated as $\widehat P(t, s, u, v)^\top \widehat \theta_{n, g\lambda}$.
\paragraph{Second step: marginal integrations}
In the second step, we first recover the $\lambda$ function. Recalling the scale normalization in (ref), the estimator of $\lambda(u, v)$ can be naturally defined by
where $P^*_{TS} \coloneqq \int_{\mathcal{TS}}P_{TS}(t,s)\mathrm{d}t\mathrm{d}s$.
For the estimation of $g(t,s)$, note that simply dividing $\widehat P(t, s, u, v)^\top \widehat \theta_{n, g\lambda}$ by $\widehat \lambda_{n}(u, v)$ results in an inefficient estimator and is not well-defined if $\widehat \lambda_{n}(u, v)$ is close to zero. To obtain a more efficient estimator utilizing the information at all $(u,v)$ values, we apply a least-squares principle to $\widehat P(t, s, u, v)^\top \widehat \theta_{n, g\lambda} \approx \widehat \lambda_n(u, v) g(t,s)$; that is, we define our estimator $\widehat g_n(t,s)$ as the solution of $\min_{g^*} \int_0^1 \int_0^1 \left(\widehat P(t,s,u,v)^\top \widehat \theta_{n,g\lambda} - \widehat \lambda_n(u,v) g^* \right)^2 \mathrm{d}u\mathrm{d}v$. By simple calculation, we can find that the solution has a closed form expression:
where $\widehat P^*_{UV} \coloneqq \int_0^1\int_0^1 \widehat P_{UV}(u, v) \widehat \omega_n(u, v) \mathrm{d}u \mathrm{d}v$ and $\widehat \omega_n(u, v) \coloneqq \widehat\lambda_n(u, v)/\int_0^1\int_0^1 \widehat \lambda_n(u', v')^2 \mathrm{d}u'\mathrm{d}v'$.
\paragraph{Third step: kernel smoothing}
In the final step, we re-estimate $\beta$ using a local linear kernel regression. Let $W$ be a symmetric univariate kernel density function and $\{(h_T, h_S)\}$ be the sequence of bandwidths tending to zero as $n$ increases. The bandwidths may depend on the evaluation point $(t,s,x)$ of the CATR parameter in general (see Remark (ref)). Further, we define $X_i(t,s) \coloneqq (X_i^\top, X_i^\top(T_i - t)/h_T, X_i^\top(S_i - s)/h_S)^\top$, and $W_{i}(t, s) \coloneqq (h_T h_S)^{-1}W((T_i - t) / h_T)W( (S_i - s) /h_S)$ for a given $(t,s) \in \mathcal{TS}$. Then, our final estimator $\widetilde \beta_n(t, s)$ of $\beta(t,s)$ is obtained by
where $\mathbb{S}_{dx} \coloneqq [I_{dx}, \: \bm{0}_{dx \times 2dx}]$, and the estimator of $\mathrm{CATR}(t,s,x)$ is given by $\widehat{\mathrm{CATR}}_n(t,s,x) \coloneqq x^\top \widetilde \beta_n(t, s)$.
To investigate the asymptotic properties of the proposed estimators, we first introduce assumptions on the dependence structure underlying the data. Typical applications of treatment spillover models would be data from school classes, households, working places, neighborhoods, municipalities, etc., which are usually subject to some local dependence. To account for such dependency, we assume that all $N$ agents are located in a (latent) $d$-dimensional space $\mathcal{D} \subseteq \mathbb{R}^d$ for some $d < \infty$. For spatial data applications, the space $\mathcal{D}$ is typically defined by a geographical space with $d = 2$. For non-spatial data, it is possible that $\mathcal{D}$ is a space of general social and demographic characteristics (or a mixture of geographic space and such spaces), and in this case we should view it as an embedding of individuals rather than their actual locations. The set of observation locations, which we denote by $\mathcal{D}_N \subset \mathcal{D}$, may differ across different $N$. With a little abuse of notation, for each $i$, we use the same $i$ to denote his/her location. Let $\Delta(i,j)$ denote the Euclidean distance between $i$ and $j$.
Assumption (ref)(i) rules out the infill asymptotics (cressie1993statistics). Assumption (ref)(ii) is a “homophily” assumption in the space $\mathcal{D}$. This may not be too restrictive since the definition of $\mathcal{D}$ can be freely modified depending on the context.\footnote{ However, note that assuming that the space $\mathcal{D}$ is a subset of Euclidean space is also a restriction. For applications with a tree-like network structure (such as supply chain networks), $\mathcal{D}$ would more naturally fit in a hyperbolic space. For a comprehensive review of the implications of the choice of geometry in network data analysis, see smith2019geometry. } This framework would accommodate many empirically relevant situations. For example, clustered sample data with finite cluster size can be seen as a special case of it. An important implication from these two assumptions is that the size of reference group is uniformly bounded above by some constant of order $\overline{\Delta}^d$. Thus, we do not allow for the existence of “dominant units,” which may have increasingly many interacting partners as the network grows. Assuming $\overline{\Delta}$ to be an integer is only to simplify the proof.
For the precise definition of the $\alpha$-mixing random field and the mixing coefficient in this context, see Definition (ref) in Appendix (ref) (see also jenish2009central). Combined with Assumption (ref)(ii), Assumption (ref)(i) implies that $\{S_i\}$, $\{V_i\}$, and $\{\widetilde \epsilon_i\}$ are also $\alpha$-mixing processes. For Assumption (ref), the second assumption restricts the dependence structure in the treatment heterogeneity, and the third requires that all potential correlations among $\epsilon_i$'s are through the correlations of the other individual characteristics. Note that $\widetilde \epsilon_i$ can be written as $\widetilde \epsilon_i = g(T_i, S_i)[\epsilon_i - \lambda(U_i, V_i)]$. Thus, with this assumption we impose that the $\widetilde \epsilon_i$'s also have the fourth order moments and that they are independent of each other conditional on the individual characteristics.
Assumption (ref)(i) is introduced for analytical simplicity, which imposes some restrictions on the choice of basis functions. For example, B-spline basis and Fourier series satisfy this assumption. The uniform boundedness implies that $\sup_{t \in \mathcal{T}}||P_T(t)|| = O(\sqrt{K_T})$, and the similar result applies to $P_S$, $P_U$, and $P_V$. Assumptions (ref)(iii)--(v) are technical conditions to derive the uniform convergence rate for the second-stage estimators. Assumption (ref) is a standard non-singularity condition, which essentially serves as the identification condition for the CATR parameter. Note that this can hold even when the support of the IVs is small.
The constants $(\mu_\beta, \mu_g, \mu_\lambda)$ generally depends on the choice of the basis function and the dimension and the smoothness of the function to be approximated. For example, if the target function belongs to a $k$-dimensional H\"older class with smoothness $p$, it typically holds that $\mu = p/k$ for splines, wavelets, etc (chen2007large). Here, define $\theta^*_{g\lambda} \coloneqq \theta_g^* \otimes \theta_\lambda^*$ so that $\overline{P}(t,s,u,v)^\top \theta^*_{g\lambda} = (P_{TS}(t, s)^\top \theta_g^*) \cdot (\overline{P}_{UV}(u,v)^\top \theta_\lambda^*)$ holds, where $\overline{P}(t,s,u,v) \coloneqq P_{TS}(t, s) \otimes \overline{P}_{UV}(u, v)$. Then, we can easily see that
For Assumption (ref), if we adopt a parametric model specification for the treatment equation in (ref), then the assumption holds with $\nu = 1/2$. As discussed in Remark (ref), we require $\nu$ to be at least larger than $1/3$ under optimal bandwidths. Assumption (ref)(i) is used to establish a matrix law of large numbers. Recalling that $(K_{TS}, K_{UV}) \asymp \kappa_n$, the first part of the assumption requires that $K_{TS}$ and $K_{UV}$ must grow slower than $n^{1/4}$. It is easy to see that $\zeta_{\dagger}$ gives the order of $||\partial P_{UV}(u,v)/\partial u||$ and $||\partial P_{UV}(u,v)/\partial v||$. For example, when one uses a tensor product B-splines, it can be shown that $\zeta_{\dagger} = O(\kappa_n)$ (see, e.g., hoshino2021treatment). Condition (ii) implies the first part of (i). This assumption can be relaxed if we can strengthen Assumption (ref) so that the error terms $\{\widetilde \epsilon_i\}$ have the moments of order higher than four.
Now, in the following theorem, we derive the convergence rate for the first-stage series estimator.
The above results should be standard in the literature. In particular, result (i) indicates that we can estimate $\mathrm{CATR}(t,s,x)$ consistently by $x^\top\widehat \theta_{n, \beta}$, although less efficiently when compared with our final estimator. The magnitude of $\tau^*_n$ depends not only on the penalty matrix $\bm{D}$ but also on the choice of the basis function. Under $\rho_{\max}(\bm{D}) = O(1)$, it is clear that the upper bound for $\tau^*_n$ is $O(\tau_n\kappa_n)$. If the coefficients are decaying in the order of series, $\tau^*_n = O(\tau_n)$ would hold, as in han2020nonparametric.
For later use, we introduce the following miscellaneous assumptions.
In the next theorem, we establish the convergence rates for the second-stage estimators.
The above theorem clearly shows that the marginal integration procedure successfully resolves the slower convergence of the first-stage estimator. Note however that both estimators $\widehat \lambda_n$ and $\widehat g_n$ do not attain the optimal uniform convergence rate of stone1982optimal.\footnote{ Whether our marginal integration-based estimation method can achieve the optimal rate is left as an open question. For example, for $\widehat \lambda_n$, in order to achieve the optimal uniform convergence rate, the bias term should be of order just $O(K_{UV}^{-p/2})$, where $p$ is a smoothness parameter (see the discussion given after Assumption (ref)). For the attainability of the optimal uniform rate for standard linear series regression estimators, see huang2003local, belloni2015some, and chen2015optimal, among others. }
We next derive the limiting distribution of our final estimator for $\mathrm{CATR}(t,s,x)$, where $(t,s)$ is a given interior point of $\mathcal{TS}$. Let $f_i(t,s)$ be the joint density of $(T_i, S_i)$, and define
Assumption (ref) is fairly standard. The compact support assumption is used just for simplicity, and it can be dropped at the cost of lengthier proof. For Assumption (ref), it will be later shown that condition (i) is the optimal rate for the bandwidths. We use condition (ii) to ensure that the final CATR estimator becomes oracle efficient.
Let $\mathbf{1}_i(t,s) \coloneqq \mathbf{1}\left\{\frac{|T_i - t|}{h_T} \le c_W , \frac{|S_i - s|}{h_S} \le c_W \right\}$, $\bm{I}_n(t,s) \coloneqq \text{diag}\left(\frac{\mathbf{1}_1(t,s)}{h_T h_S}, \ldots, \frac{\mathbf{1}_n(t,s)}{h_T h_S}\right)$, $\bm{P}_{n,TS} \coloneqq (P_{TS}(T_1, S_1), \ldots, P_{TS}(T_n, S_n))^\top$, and $\bm{P}_{n,UV} \coloneqq (\overline{P}_{UV}( U_1, V_1), \ldots, \overline{P}_{UV}(U_n, V_n))^\top$.
Assumptions (ref), (ref), and (ref) are technical requirements to derive the asymptotic normality of our CATR estimator. In particular, Assumption (ref)(ii) and (iii) are introduced to utilize the central limit theorem for mixing random fields developed by jenish2009central in our context. Now we are ready to state our main theorem:
The proof of Theorem (ref) is straightforward from Lemma (ref), and thus is omitted. In Lemma (ref), we show that the estimation error for our CATR estimator caused by the estimations of $g$ and $\lambda$ are of order $o_P((n h_T h_S)^{-1/2})$. Thus, the asymptotic distribution presented in the theorem is in fact equivalent to that obtained when the estimators $\widehat g_n$ and $\widehat \lambda_n$ are replaced by their true counterparts; that is, our CATR estimator has an oracle property.
It is often considered that unemployment and crime are endogenously related because of their simultaneity (e.g., levitt2001alternative). In economic theory, criminal activity is typically characterized by the balance between the cost and benefit of committing illegal activities (becker1968crime). Thus, poor local labor market conditions result in a relative increase in the benefits of crime, increasing the number of crime incidents. On the other hand, crime drives away business owners and customers, exacerbating the working condition. In addition, from the viewpoint of criminals, if their communities are economically deteriorating with high unemployment, they might commit crimes in more “beneficial” neighborhoods. This would suggest the need to account for the spatial spillover effects of unemployment.
There are several prior studies that investigate the relationship between unemployment and crime using some IV-based methods (e.g., raphael2001identifying, lin2008does, altindag2012crime). For example, using U.S. state-level data, raphael2001identifying conducted two-stage least squares regression analysis with military spending and oil costs as the IVs for the unemployment. Their estimation results suggest that unemployment is indeed an important determinant of property crime rates. In this paper, we present an empirical analysis on the causal effect of local unemployment rate on crime based on Japanese city-level data. Our empirical study aims at extending the earlier works in two ways. Our model is based on a more flexible potential outcome framework and allows for the existence of spillover effects from neighboring regions.
The variables used and their definitions are summarized in Table (ref). The outcome variable of interest is the city-level crime rate, which is based on the number of crimes recorded in each city in 2006. The crime data cover all kinds of criminal offenses (the breakdown is unfortunately unavailable). The treatment variable is the regional unemployment rate as of 2005. As an IV for the unemployment rate, we employ the availability of child day-care facilities in 2006. It would be legitimate to assume that the availability of childcare facilities does not directly affect the crime rate. Meanwhile, in the labor economics literature, there is empirical evidence that expanding childcare services is effective in increasing female labor participation (e.g., bauernschuster2015public). Thus, the availability of childcare facilities can be viewed as an indicator of the local working environment, particularly for young females, and would contribute to reducing the total unemployment rate. For other control variables, we include population density, annual retail sales per capita, the ratio of elderly people, and the ratio of single households. The sales data are as of 2006, and the others are those in 2005.\footnote{ All this information is freely available from e-Stat (a portal site for Japanese Government Statistics). \url{https://www.e-stat.go.jp/en/regional-statistics/ssdsview/municipality} }
The network $\bm{A}_N$ is defined by whether the cities share a common boundary. Specifically, assuming that there is a limit on the number of cities each city can interact with, we set $A_{i,j} = 1$ if city $j$ is adjacent to $i$ and in the $k$-nearest neighbors of $i$. Below, we report the results when $k = 2$ for illustration. We also have tried several different specifications for $\bm{A}_N$, and confirmed that the results are overall similar (for more details, see Appendix (ref)). For all cases, the treatment spillover variable $S_i$ is defined by $i$'s reference-group mean, say $\overline{\textit{unempl}}$. Then, the model estimated is as follows:
Note that since the treatment equation does not involve any network interactions explicitly, the CQR estimation can be implemented using all data regardless of the specification of $\bm{A}_N$. For the estimation of the CATR parameter, we need to select a subsample to maintain distributional homogeneity of $(U,V)$. To this end we focus on the cities that have exactly two interacting partners (recalling that $k = 2$). After excluding observations with missing data, the estimation of CATR was performed on a sample of size 1773. The estimation procedure is the same as that in the Monte Carlo experiments in Appendix (ref). That is, we set the penalty parameter to $\tau_n = 5/n$ and use the B-spline basis with two internal knots.\footnote{ We have confirmed that the results reported here have a certain robustness to other choices of penalty parameters and basis orders. However, we have also observed that when the number of internal knots is three or higher, it is better to employ some regularized estimator when computing the second derivatives appearing in (ref) to stabilize the estimates. } The descriptive statistics of the data are summarized in Table (ref) in Appendix (ref).
The estimation results for the treatment equation are provided in Table (ref) in Appendix (ref), where we can find that childcare is significantly negatively related to the unempl variable, as expected. Figure (ref) presents the estimated $\mathrm{CATR}(t,s,x)$ for different values of $(t,s,x)$. In each panel, $t$ ranges over 0.1 to 0.9 empirical quantiles of $\{T_i\}$, and $s$ is either at 0.2 or 0.6 quantile of $\{S_i\}$. Because of the correlation between $T$ and $S$, the CATR estimates at more extreme $t$ cannot be estimate reliably and thus they are not reported (see Figure (ref) in Appendix (ref) for the joint density of $(T,S)$). For the value of $x$, we evaluate at the empirical median in the left panel. We consider two more cases for $x$: the median value for the cities in the bottom 20% level of population density among all cities (middle panel), and that for the cities in the top 20% level (right panel). The former would be considered as a typical city in a rural area, while the latter as a typical city in an urban area. In the figure, we also report the results obtained when $T$ and $S$ are treated as exogenous (this corresponds to $\widehat{\mathrm{CATR}}_n^\text{ex}$ estimator given in Appendix (ref)).
Now, we report our main empirical findings. First, we can observe that the CATR weakly increases in general as the unemployment rate increases. That is, as the number of unemployed people increases, the crime rate tends to increase, which is consistent with the findings in the prior studies. Second, in the left and right panels, the CATR with large $S$ tends to be significantly greater than the CATR with small $S$, indicating that the average unemployment rate of surrounding cities does affect the city's own crime rate in such cases. However, interestingly, the spillovers are no longer prominent when the city is in a rural area. One interpretation is that the cities in such areas are relatively independent from other cities and villages, and thus the spillover effects may be less impactful. A more comprehensive figure that summarizes the CATR estimates when $S$ is at 0.2, 0.3, \ldots, 0.8 quantiles is presented in Figure (ref) in Appendix (ref), where we can more clearly observe these tendencies. Lastly, we find that ignoring the potential endogeneity of $(\textit{unempl}, \overline{\textit{unempl}})$ leads to certain differences in the estimates, indicating the importance of accounting for the endogeneity.
In this paper, we considered a continuous treatment effect model that admits potential treatment spillovers through social networks and the endogeneity for both one's own treatment and that of others. We proved that the CATR, the conditional expectation of the potential outcome, can be nonparametrically identified under some functional form restrictions and the availability of appropriate IVs. We also considered a more empirically tractable semiparametric treatment model and proposed a three-step procedure for estimating the CATR. The consistency and asymptotic normality of the proposed estimator were established under certain regularity conditions. As an empirical illustration, using Japanese city-level data, we investigated the causal effects of unemployment rate on the number of crime incidents. As a result, we found that the unemployment rate indeed tends to increase the number of crimes and that, interestingly, the unemployment rates of neighboring cities matter more for non-rural cities than rural cities. These results illustrate the usefulness of the proposed method.
\setcounter{page}{1}