The exact contents of citations.db main_text.text for this paper — one flattened LaTeX string, title through conclusion, appendix excluded, unmodified except for removing email addresses. This is what our citation measures are computed over.
91,094 characters
Threshold Regression with Nonparametric Sample Splitting
\title{Threshold Regression with Nonparametric Sample Splitting\thanks{
We are grateful to Xiaohong Chen, Jonathan Dingel, Bo
Honor\'{e}, Sokbae Lee, Yuan Liao, Francesca Molinari, Ingmar Prucha, Myung
Seo, Ping Yu, and participants at numerous seminar/conference presentations
for very helpful comments. Financial supports from the Appleby-Mosher grant
and the CUSE grant are highly appreciated.}}
\author{\textsc{Yoonseok Lee}\thanks{\textit{Address}: Department of
Economics and Center for Policy Research, Syracuse University, 426 Eggers
Hall, Syracuse, NY 13244. \textit{E-mail}: \texttt{[email removed]}}
\\
Syracuse University \and \textsc{Yulong Wang}\thanks{\textit{Address}:
Department of Economics and Center for Policy Research, Syracuse University,
127 Eggers Hall, Syracuse, NY 13244. \textit{E-mail}: \texttt{
[email removed]}} \\
Syracuse University}
\date{January 2021}
\maketitle
\begin{abstract}
\noindent This paper develops a threshold regression model where an unknown
relationship between two variables nonparametrically determines the
threshold. We allow the observations to be cross-sectionally dependent so
that the model can be applied to determine an unknown spatial border for
sample splitting over a random field. We derive the uniform rate of
convergence and the nonstandard limiting distribution of the nonparametric
threshold estimator. We also obtain the root-n consistency and the
asymptotic normality of the regression coefficient estimator. Our model has
broad empirical relevance as illustrated by estimating the tipping point in
social segregation problems as a function of demographic characteristics;
and determining metropolitan area boundaries using nighttime light intensity
collected from satellite imagery. We find that the new empirical results are
substantially different from those in the existing studies.
\vspace{0.15in}
\noindent \noindent \textit{\noindent Keywords}: threshold regression,
sample splitting, nonparametric, random field, tipping point, metropolitan
area boundary.
\noindent \textit{\noindent JEL Classifications}: C14, C21, C24, R1
\end{abstract}
\thispagestyle{empty}\setcounter{page}{0}\newpage
\setcounter{page}{1}
\section{Introduction}
Sample splitting and threshold regression models have spawned a vast
literature in econometrics and statistics. Existing studies typically
specify the sample splitting criteria in a parametric way as whether a
single random variable or a linear combination of variables crosses some
unknown threshold. See, for example, \cite{Hansen00a}, \cite{Caner04}, \cite
{SeoLinton07}, \cite{LeeSeoShin11}, \cite{LiLing12}, \cite{Yu12}, \cite
{LLSS18}, \cite{Hidalgo19}, and \cite{Yu19}. In this paper, we study a novel
extension to consider a \textit{nonparametric} sample splitting model. Such
an extension leads to new theoretical results and substantially generalizes
the empirical applicability of threshold models.
Specifically, we consider a model given by
\begin{equation}
y_{i}=x_{i}^{\top }{\Greekmath 010C} _{0}+x_{i}^{\top }{\Greekmath 010E} _{0}\mathbf{1}\left[
q_{i}\leq {\Greekmath 010D} _{0}\left( s_{i}\right) \right] +u_{i} \label{model}
\end{equation}
for $i=1,\ldots ,n$, where $\mathbf{1}\left[ \cdot \right] $ is the binary
indicator. In this model, the marginal effect of $x_{i}$ to $y_{i}$ can be
different across $i$ as $({\Greekmath 010C} _{0}+{\Greekmath 010E} _{0})$ or ${\Greekmath 010C} _{0}$ depending
on whether $q_{i}\leq {\Greekmath 010D} _{0}\left( s_{i}\right) $ or not. The threshold
function ${\Greekmath 010D} _{0}(\cdot )$ is unknown, and the main parameters of
interest are ${\Greekmath 010C} _{0}$, ${\Greekmath 010E} _{0}$, and ${\Greekmath 010D} _{0}(\cdot )$. The
novel feature of this model is that the sample splitting is determined by an
unknown relationship between two variables $q_{i}$ and $s_{i}$, and their
relationship is characterized by the nonparametric threshold function $
{\Greekmath 010D} _{0}(\cdot )$. In contrast, the classical threshold regression models
assume ${\Greekmath 010D} _{0}\left( \cdot \right) $ to be a constant or a linear
index. Our new specification can cover interesting cases that have not been
studied. For example, we can consider the threshold to be heterogeneous and
specific to each observation $i$ if we see ${\Greekmath 010D} _{0}\left( s_{i}\right)
={\Greekmath 010D} _{0i}$; or the threshold to be determined by the direction of some
moment condition ${\Greekmath 010D} _{0}(s_{i})=\mathbb{E}[q_{i}|s_{i}]$. Apparently,
when ${\Greekmath 010D} _{0}(s)={\Greekmath 010D} _{0}$ or ${\Greekmath 010D} _{0}(s)={\Greekmath 010D} _{0}s$ for some
parameter ${\Greekmath 010D} _{0}$ and $s\neq 0$, it reduces to the standard threshold
regression model.
The new model is motivated by the following two applications: estimating
potentially heterogeneous thresholds in public economics and determining
spatial sample splitting in urban economics. The first one is about the
tipping point model proposed by \cite{Schelling71}, who analyzes the
phenomenon that a neighborhood's white population substantially decreases
once the minority share exceeds a certain threshold, called the tipping
point. \cite{Card08} empirically estimate the tipping point model by
considering the constant threshold regression, $y_{i}={\Greekmath 010C} _{10}+{\Greekmath 010E}
_{10}\mathbf{1}\left[ q_{i}\leq {\Greekmath 010D} _{0}\right] +x_{2i}^{\top }{\Greekmath 010C}
_{20}+u_{i}$, where $y_{i}$ is the white population change in a decade and $
q_{i}$ is the initial minority share in the $i$th tract. The parameters $
{\Greekmath 010E} _{10}$ and ${\Greekmath 010D} _{0}$ denote the change size and the threshold,
respectively. In Section VII of \cite{Card08}, however, they find that the
tipping point ${\Greekmath 010D} _{0}$ varies depending on the attitudes of white
residents toward the minority. This finding raises the concern on the
constant threshold model and motivates us to study the more general model (
\ref{model}) by specifying the tipping point ${\Greekmath 010D} _{0}$ as a
nonparametric function of local demographic characteristics. We estimate
such a tipping function in Section \ref{Section tipping}.
For the second application, we use the model (\ref{model}) to define
metropolitan area boundaries, which is a fundamental problem in urban
economics. Recently, many studies propose to use nighttime light intensity
collected from satellite imagery to define the metropolitan area. They set
an \textit{ad hoc} level of light intensity as a threshold and categorize a
pixel in the satellite imagery as a part of the metropolitan area if the
light intensity of that pixel is higher than the threshold. See, for
example, \cite{Rozenfeld11}, \cite{Henderson12}, \cite{Dingel19}, and \cite
{Vogel19}. In contrast, the model (\ref{model}) can provide a data-driven
guidance of choosing the intensity threshold from the econometric
perspective, if we let $y_{i}$ as the light intensity in the $i$th pixel and
$(q_{i},s_{i})$ as the location information of that pixel (more precisely,
the coordinate of a point on a rotated map as described in Section \ref
{Section contour}). In Section \ref{Section boundary}, we estimate the
metropolitan area of Dallas, Texas, especially its development from 1995 to
2010, and find substantially different results from the conventional
approaches. To the best of our knowledge, this is the first study to
nonparametrically determine the metropolitan area using a threshold model.
We develop a two-step estimation procedure of (\ref{model}), where we
estimate ${\Greekmath 010D} _{0}\left( \cdot \right) $ by the local constant least
squares. Under the shrinking threshold asymptotics as in \cite{Bai97b}, \cite
{Bai98}, and \cite{Hansen00a}, we show that the nonparametric estimator $
\widehat{{\Greekmath 010D} }(\cdot )$ is uniformly consistent and has a highly
nonstandard limiting distribution. Based on such distribution, we develop a
pointwise specification test of ${\Greekmath 010D} _{0}(s)$ for any given $s$, which
enables us to construct a confidence interval by inverting the test.
Besides, the parametric part $(\widehat{{\Greekmath 010C} }^{\top },\widehat{{\Greekmath 010E} }
^{\top })^{\top }$ is shown to satisfy the root-$n$ asymptotic normality.
We highlight some novel technical features of the new estimator as follows.
First, since the nonparametric function ${\Greekmath 010D} _{0}\left( \cdot \right) $
is inside the indicator function, technical proofs of the asymptotic results
are non-standard.\ In particular, we establish the uniform rate of
convergence of $\widehat{{\Greekmath 010D} }\left( \cdot \right) $, which involves
substantially more complicated derivations than the standard (constant)
threshold regression model. Second, we find that, unlike the standard kernel
estimator, $\widehat{{\Greekmath 010D} }(\cdot )$ is asymptotically unbiased even if
the optimal bandwidth is used. Also, when the change size ${\Greekmath 010E} _{0}$
shrinks very slowly, the optimal rate of convergence of $\widehat{{\Greekmath 010D} }
(\cdot )$ becomes close to the root-$n$ rate. In the standard kernel
regression, such a fast rate of convergence can be obtained when the unknown
function is infinitely differentiable, while we only require the
second-order differentiability of ${\Greekmath 010D} _{0}\left( \cdot \right) $. Third,
to limit the effect of estimating ${\Greekmath 010D} _{0}\left( \cdot \right) $ to $(
\widehat{{\Greekmath 010C} }^{\top },\widehat{{\Greekmath 010E} }^{\top })^{\top }$, we propose to
use the observations that are sufficiently away from the estimated threshold
in the second-step parametric estimation. The choice of this distance is
obtained by the uniform convergence rate of $\widehat{{\Greekmath 010D} }\left( \cdot
\right) $. Fourth, we let the variables be cross-sectionally dependent by
considering the strong-mixing random field as in \cite{Conley99} and \cite
{Conley07}. This generalization allows us to study nonparametric sample
splitting of spatial observations. For instance, if we let $(q_{i},s_{i})$
correspond to the geographical location (i.e., latitude and longitude on the
map), then the threshold $\mathbf{1}\left[ q_{i}\leq {\Greekmath 010D} _{0}\left(
s_{i}\right) \right] $ identifies the unknown border yielding a
two-dimensional sample splitting. In more general contexts, the model can be
applied to identify social or economic segregation over interacting agents.
Finally, noting that $\mathbf{1}\left[ q_{i}\leq {\Greekmath 010D} _{0}\left(
s_{i}\right) \right] $ can be considered as the special case of $\mathbf{1}
\left[ g_{0}\left( q_{i},s_{i}\right) \leq 0\right] $ when $g_{0}$ is
monotonically increasing in $q_{i}$, we discuss how to extend the proposed
method to such a more general case that leads to a threshold contour model.
The rest of the paper is organized as follows. Section \ref{Section
estimation} sets up the model, establishes the identification, and defines
the estimator. Section \ref{Section asymptotics} derives the asymptotic
properties of the estimators and develops a likelihood ratio test of the
threshold function. Section \ref{Section contour} describes how to extend
the main model to estimate a threshold contour. Section \ref{Section
simulation} studies small sample properties of the proposed statistics by
Monte Carlo simulations. Section \ref{Section empirics} applies the new
method to estimate the tipping point function and to determine metropolitan
areas. Section \ref{Section conclusion} concludes this paper with some
remarks. The main proofs are in the Appendix, and all the omitted proofs are
collected in the supplementary material.
We use the following notations. Let $\rightarrow _{p}$ denote convergence in
probability, $\rightarrow _{d}$ convergence in distribution, and $
\Rightarrow $ weak convergence of the underlying probability measure as $
n\rightarrow \infty $. Let $\left\lfloor r\right\rfloor $ denote the biggest
integer smaller than or equal to $r$, $\mathbf{1}[E]$ the indicator function
of a generic event $E$, and $\left\Vert A\right\Vert $ the Euclidean norm of
a vector or matrix $A$. For any set $B$, let $|B|$ as the cardinality of $B$.
\section{Model Setup\label{Section estimation}}
We assume spatial processes located on an evenly spaced lattice $\Lambda
\subset
\mathbb{R}
^{2}$, following \cite{Conley99}, \cite{Conley07}, and \cite
{CarbonFrancqTran2007}.\footnote{
It can be extended to an unevenly spaced lattice as in \cite{Bolthausen82}
and \cite{Jenish09} with substantially more complicated notations.\ (cf.\
footnote 9 in \cite{Conley99}).} We consider the threshold regression model
given by (\ref{model}), which is
\begin{equation*}
y_{i}=x_{i}^{\top }{\Greekmath 010C} _{0}+x_{i}^{\top }{\Greekmath 010E} _{0}\mathbf{1}\left[
q_{i}\leq {\Greekmath 010D} _{0}\left( s_{i}\right) \right] +u_{i}\relax\protect\ifmmode\expandafter\text@\else\expandafter\mbox\fi{,}
\end{equation*}
where the observations $\{(y_{i},x_{i}^{\top },q_{i},s_{i})^{\top }\in
\mathbb{R}^{1+\dim (x)+1+1};i\in \Lambda _{n}\}$ are a triangular array of
real random variables defined on some probability space with $\Lambda _{n}$
being a fixed sequence of finite subsets of $\Lambda $. In this setup, the
cardinality of $\Lambda _{n}$, $n=|\Lambda _{n}|$, is the sample size and $
\sum_{i\in \Lambda _{n}}$ denotes the summation of all observations. For
readability, we postpone the regularity conditions on $\Lambda _{n}$ in
Assumption A later. The threshold function ${\Greekmath 010D} _{0}:\mathbb{R\rightarrow
R}$ as well as the regression coefficients ${\Greekmath 0112} _{0}=({\Greekmath 010C} _{0}^{\top
},{\Greekmath 010E} _{0}^{\top })^{\top }\in \mathbb{R}^{2\dim (x)}$ are unknown, and
they are the parameters of interest.\footnote{
The main results of this paper can be extended to consider multi-dimensional
$s_{i}$ using multivariate kernels. However, we only consider the scalar
case for the expositional simplicity. Furthermore, the results are readily
generalized to the case where only a subset of parameters differ between
regimes.} Since we consider a shrinking threshold effect, the parameter $
{\Greekmath 010E} _{0}$ is to depend on the sample size $n$ as in Assumption A-(ii)
below; hence ${\Greekmath 010E} _{0}$ and ${\Greekmath 0112} _{0}$ should be written as ${\Greekmath 010E}
_{n0}$ and ${\Greekmath 0112} _{n0}$, respectively. However, we write ${\Greekmath 010E} _{0}$ and
${\Greekmath 0112} _{0}$ for simplicity. We let $\mathcal{Q}\subset
\mathbb{R}
$ and $\mathcal{S}\subset
\mathbb{R}
$ denote the supports of $q_{i}$ and $s_{i}$, respectively. Suppose the
space of ${\Greekmath 010D} _{0}\left( s\right) $ for any $s$ is a compact set $\Gamma
\subset
\mathbb{R}
$.
First, we establish the identification, which requires the following
conditions.
\paragraph{Assumption ID}
\begin{description}
\item \textit{(i) }$\mathbb{E}\left[ u_{i}|x_{i},q_{i},s_{i}\right] =0$
\textit{\ almost surely.}
\item \textit{(ii) }$\mathbb{E}\left[ x_{i}x_{i}^{\top }\right] >\mathbb{E}
\left[ x_{i}x_{i}^{\top }\mathbf{1}\left[ q_{i}\leq {\Greekmath 010D} \right] \right]
>0 $\textit{\ for any }${\Greekmath 010D} \in \Gamma $.
\item \textit{(iii) For any }$s\in \mathcal{S}$, \textit{there exists }$
{\Greekmath 0122} (s)>0$\textit{\ such that }${\Greekmath 0122} (s)<\mathbb{P}\left(
q_{i}\leq {\Greekmath 010D} _{0}(s_{i})|s_{i}=s\right) <1-{\Greekmath 0122} (s)$\textit{\ and
}${\Greekmath 010E} _{0}^{\top }\mathbb{E}\left[ x_{i}x_{i}^{\top }|q_{i}=q,s_{i}=s
\right] {\Greekmath 010E} _{0}>0$ \textit{for all} $\left( q,s\right) \in \mathcal{Q}
\times \mathcal{S}$.
\item \textit{(iv) }$q_{i}$\textit{\ is continuously distributed with its
conditional density }$f(q|s)$\textit{\ satisfying that }$
0<C_{1}<f(q|s)<C_{2}<\infty $\textit{\ for all} $\left( q,s\right) \in
\Gamma \times \mathcal{S}$\textit{\ and some constants }$C_{1}$ and $C_{2}$.
\end{description}
\medskip
Assumption ID is mild. The condition (i) excludes endogeneity, and (ii) is
the full rank condition to identify the global parameters ${\Greekmath 010C} _{0}$ and $
{\Greekmath 010E} _{0}$. The conditions (ii) and (iii) require that the location of the
threshold is not on the boundary of the support of $q_{i}$ for any $s\in
\mathcal{S}$, which is inevitable for identification and has been commonly
assumed in the existing threshold literature (e.g., \cite{Hansen00a}). If $
{\Greekmath 010D} _{0}(s)$ reaches the boundary of $q_{i}$ for some $s\in \mathcal{S}$,
then no observation can be generated from one side of the threshold function
at this $s$, and identification is failed. The second condition in (iii)
assumes the coefficient change exists (i.e., ${\Greekmath 010E} _{0}\neq 0$). Note that
it does not require $\mathbb{E}\left[ x_{i}x_{i}^{\top }|q_{i}=q,s_{i}=s
\right] $ to be of full rank, and hence $q_{i}$ or $s_{i}$ can be one of the
elements of $x_{i}$ (e.g., the threshold autoregressive model by \cite
{Tong83}) or a linear combination of $x_{i}$. The condition (iv) requires
the conditional density of $q_{i}$ given any $s_{i}$ is positive and bounded
in $\Gamma $.
Under Assumption ID, the following theorem establishes the identification of
the semiparametric threshold regression model (\ref{model}).
\begin{theorem}
\label{Thm id}Under Assumption ID, the parameters $({\Greekmath 010C} _{0}^{\top
},{\Greekmath 010E} _{0}^{\top })^{\top }$ are the unique minimizer of $\mathbb{E}
[(y_{i}-x_{i}^{\top }{\Greekmath 010C} -x_{i}^{\top }{\Greekmath 010E} \mathbf{1}\left[ q_{i}\leq
{\Greekmath 010D} \right] )^{2}]$ for any ${\Greekmath 010D} \in \Gamma $, and the threshold
function ${\Greekmath 010D} _{0}\left( s\right) $ is the unique minimizer of $\mathbb{E}
[(y_{i}-x_{i}^{\top }{\Greekmath 010C} _{0}-x_{i}^{\top }{\Greekmath 010E} _{0}\mathbf{1}\left[
q_{i}\leq {\Greekmath 010D} (s_{i})\right] )^{2}|s_{i}=s]$ for each given $s\in
\mathcal{S}$.
\end{theorem}
Given identification, we proceed to estimate this semiparametric model in
two steps. First, for given $s\in \mathcal{S}$, we fix ${\Greekmath 010D}
_{0}(s)={\Greekmath 010D} $ and obtain $\widehat{{\Greekmath 010C} }\left( {\Greekmath 010D} ;s\right) $ and $
\widehat{{\Greekmath 010E} }\left( {\Greekmath 010D} ;s\right) $ by the local constant least
squares conditional on ${\Greekmath 010D} $:
\begin{equation}
(\widehat{{\Greekmath 010C} }\left( {\Greekmath 010D} ;s\right) ^{\intercal },\widehat{{\Greekmath 010E} }
\left( {\Greekmath 010D} ;s\right) ^{\intercal })^{\intercal }=\arg \min_{{\Greekmath 010C} ,{\Greekmath 010E}
}Q_{n}\left( {\Greekmath 010C} ,{\Greekmath 010E} ,{\Greekmath 010D} ;s\right) \relax\protect\ifmmode\expandafter\text@\else\expandafter\mbox\fi{,} \label{first}
\end{equation}
where
\begin{equation}
Q_{n}\left( {\Greekmath 010C} ,{\Greekmath 010E} ,{\Greekmath 010D} ;s\right) =\mathop{\displaystyle \sum }_{i\in \Lambda
_{n}}K\left( \frac{s_{i}-s}{b_{n}}\right) \left( y_{i}-x_{i}^{\top }{\Greekmath 010C}
-x_{i}^{\top }{\Greekmath 010E} \mathbf{1}\left[ q_{i}\leq {\Greekmath 010D} \right] \right) ^{2}
\label{sse}
\end{equation}
for some kernel function $K\left( \cdot \right) $ and a bandwidth parameter $
b_{n}$. Then ${\Greekmath 010D} _{0}(s)$ is estimated by
\begin{equation}
\widehat{{\Greekmath 010D} }\left( s\right) =\arg \min_{{\Greekmath 010D} \in \Gamma
_{n}}Q_{n}\left( {\Greekmath 010D} ;s\right) \label{g-hat0}
\end{equation}
for given $s$, where $\Gamma _{n}=\Gamma \cap \{q_{1},\ldots ,q_{n}\}$ and $
Q_{n}\left( {\Greekmath 010D} ;s\right) $ is the concentrated sum of squares defined as
\begin{equation}
Q_{n}\left( {\Greekmath 010D} ;s\right) =Q_{n}\left( \widehat{{\Greekmath 010C} }\left( {\Greekmath 010D}
;s\right) ,\widehat{{\Greekmath 010E} }\left( {\Greekmath 010D} ;s\right) ,{\Greekmath 010D} ;s\right) \relax\protect\ifmmode\expandafter\text@\else\expandafter\mbox\fi{.
} \label{reg}
\end{equation}
To avoid any additional technical complexity, we focus on estimation of $
{\Greekmath 010D} _{0}(s)$ at $s\in \mathcal{S}_{0}\subset \mathcal{S}$ for some
compact interior subset $\mathcal{S}_{0}$ of the support, say the middle
70\% quantiles. Note that, given $s$, the nonparametric estimator $\widehat{
{\Greekmath 010D} }\left( s\right) $ can be seen as a local version of the standard
(constant) threshold regression estimator. Therefore, the computation of (
\ref{g-hat0}) requires one-dimensional grid search of the threshold for only
$n$ times as in the standard threshold regression estimation.
In the second step, we estimate the parametric components ${\Greekmath 010C} _{0}$ and $
{\Greekmath 010E} _{0}$. To minimize any potential effects from the first step
estimation, we estimate ${\Greekmath 010C} _{0}$ and ${\Greekmath 010E} _{0}^{\ast }={\Greekmath 010C}
_{0}+{\Greekmath 010E} _{0}$ using the observations that are sufficiently away from the
estimated threshold. This is implemented by considering
\begin{eqnarray}
\widehat{{\Greekmath 010C} } &=&\arg \min_{{\Greekmath 010C} }\mathop{\displaystyle \sum }_{i\in \Lambda _{n}}\left(
y_{i}-x_{i}^{\top }{\Greekmath 010C} \right) ^{2}\mathbf{1}\left[ q_{i}>\widehat{{\Greekmath 010D} }
\left( s_{i}\right) +{\Greekmath 0119} _{n}\right] \mathbf{1}[s_{i}\in \mathcal{S}_{0}]
\relax\protect\ifmmode\expandafter\text@\else\expandafter\mbox\fi{,} \label{para-b} \\
\widehat{{\Greekmath 010E} }^{\ast } &=&\arg \min_{{\Greekmath 010E} ^{\ast }}\mathop{\displaystyle \sum }_{i\in \Lambda
_{n}}\left( y_{i}-x_{i}^{\top }{\Greekmath 010E} ^{\ast }\right) ^{2}\mathbf{1}\left[
q_{i}<\widehat{{\Greekmath 010D} }\left( s_{i}\right) -{\Greekmath 0119} _{n}\right] \mathbf{1}
[s_{i}\in \mathcal{S}_{0}] \label{para-d}
\end{eqnarray}
for some constant ${\Greekmath 0119} _{n}>0$ satisfying ${\Greekmath 0119} _{n}\rightarrow 0$ as $
n\rightarrow \infty $, which is defined later. The change size ${\Greekmath 010E} $ can
be estimated as $\widehat{{\Greekmath 010E} }=\widehat{{\Greekmath 010E} }^{\ast }-\widehat{{\Greekmath 010C} }
$.
For the asymptotic behavior of the threshold estimator, the existing
literature typically assumes martingale difference arrays (e.g., \cite
{Hansen00a} and \cite{LLSS18}) or random samples (e.g., \cite{Yu12} and \cite
{Yu19}). In this paper, we allow for cross-sectional dependence by
considering spatial ${\Greekmath 010B} $-mixing processes as in \cite{Bolthausen82} and
\cite{Conley99}. More precisely, for any indices (or locations) $i,j\in
\Lambda $, we define the metric ${\Greekmath 0115} \left( i,j\right) =\max_{1\leq \ell
\leq \dim \left( \Lambda \right) }\left\vert i_{\ell }-j_{\ell }\right\vert $
and the corresponding norm $\max_{1\leq \ell \leq \dim \left( \Lambda
\right) }\left\vert i_{\ell }\right\vert $, where $i_{\ell }$ denotes the $
\ell $th component of $i$. The distance of any two subsets $\Lambda
_{1},\Lambda _{2}\subset \Lambda $ is defined as ${\Greekmath 0115} (\Lambda
_{1},\Lambda _{2})=\inf \{{\Greekmath 0115} (i,j):i\in \Lambda _{1},j\in \Lambda
_{2}\} $. We let $\mathcal{F}_{\Lambda }$ be the ${\Greekmath 011B} $-algebra generated
by a random sequence $\left( x_{i}^{\top },q_{i},s_{i},u_{i}\right) ^{\top }$
for $i\in \Lambda $ and define the spatial ${\Greekmath 010B} $-mixing coefficient as
\begin{equation}
{\Greekmath 010B} _{k,l}(m)=\sup \left\{ \left\vert \mathbb{P}\left( A\cap B\right) -
\mathbb{P}\left( A\right) \mathbb{P}\left( B\right) \right\vert :A\in
\mathcal{F}_{\Lambda _{1}}\relax\protect\ifmmode\expandafter\text@\else\expandafter\mbox\fi{, }B\in \mathcal{F}_{\Lambda _{2}}\relax\protect\ifmmode\expandafter\text@\else\expandafter\mbox\fi{, }
{\Greekmath 0115} \left( \Lambda _{1},\Lambda _{2}\right) \geq m\right\} \relax\protect\ifmmode\expandafter\text@\else\expandafter\mbox\fi{,}
\label{mix}
\end{equation}
where $\left\vert \Lambda _{1}\right\vert \leq k$ and $\left\vert \Lambda
_{2}\right\vert \leq l$. Without loss of generality, we assume ${\Greekmath 010B}
_{k,l}(0)=1$ and ${\Greekmath 010B} _{k,l}(m)$ is monotonically decreasing in $m$ for
all $k$ and $l$.
The following conditions are imposed for deriving the asymptotic properties
of our two-step estimator. Let $f\left( q,s\right) $\ be the joint density
function of $(q_{i},s_{i})$ and
\begin{eqnarray}
D\left( q,s\right) &=&\mathbb{E}[x_{i}x_{i}^{\top }|\left(
q_{i},s_{i}\right) =\left( q,s\right) ]\relax\protect\ifmmode\expandafter\text@\else\expandafter\mbox\fi{,} \label{D} \\
V\left( q,s\right) &=&\mathbb{E}[x_{i}x_{i}^{\top }u_{i}^{2}|\left(
q_{i},s_{i}\right) =\left( q,s\right) ]\relax\protect\ifmmode\expandafter\text@\else\expandafter\mbox\fi{.} \label{V}
\end{eqnarray}
\paragraph{Assumption A}
\begin{description}
\item \textit{(i) The lattice }$\Lambda _{n}\subset
\mathbb{R}
^{2}$\textit{\ is infinitely countable; all the elements in }$\Lambda _{n}$
\textit{\ are located at distances at least }${\Greekmath 0115} _{0}>1$\textit{\ from
each other (i.e., for any }$i,j\in \Lambda _{n}$, ${\Greekmath 0115} \left( i,j\right)
\geq {\Greekmath 0115} _{0}$\textit{); and }$\lim_{n\rightarrow \infty }\left\vert
\partial \Lambda _{n}\right\vert /n=0$\textit{, where }$\partial \Lambda
_{n}=\{i\in \Lambda _{n}:\exists j\not\in \Lambda _{n}$\textit{\ with }$
{\Greekmath 0115} (i,j)=1\}$\textit{.}
\item \textit{(ii) }${\Greekmath 010E} _{0}=c_{0}n^{-{\Greekmath 010F} }$\textit{\ for some }$
c_{0}\neq 0$\textit{\ and }${\Greekmath 010F} \in (0,1/2)$; $\left( c_{0}^{\top
},{\Greekmath 010C} _{0}^{\top }\right) ^{\top }$\textit{\ belongs to some compact
subset of }$\mathbb{R}^{2\dim (x)}$\textit{.}
\item \textit{(iii) }$\left( x_{i}^{\top },q_{i},s_{i},u_{i}\right) ^{\top }$
\textit{\ is strictly stationary and }${\Greekmath 010B} $\textit{-mixing with the
mixing coefficient }${\Greekmath 010B} _{k,l}(m)$\textit{\ defined in (\ref{mix}),
which satisfies that for all }$k$ and $l$,\textit{\ }${\Greekmath 010B} _{k,l}(m)\leq
C_{1}\exp (-C_{2}m)$\textit{\ for some positive constants }$C_{1}$\textit{\
and }$C_{2}$\textit{.}
\item \textit{(iv) }$0<\mathbb{E}\left[ u_{i}^{2}|x_{i},q_{i},s_{i}\right]
<\infty $\textit{\ almost surely.}
\item \textit{(v) Uniformly in }$\left( q,s\right) $\textit{, there exist
some finite constants }${\Greekmath 0127} >0$ and $C>0$\textit{\ such that }$\mathbb{E}
[||x_{i}x_{i}^{\top }||^{4+{\Greekmath 0127} }|(q_{i},s_{i})=(q,s)]<C$\textit{\ and }$
\mathbb{E}[||x_{i}u_{i}||^{4+{\Greekmath 0127} }|(q_{i},s_{i})=(q,s)]<C$.
\item \textit{(vi) }${\Greekmath 010D} _{0}:\mathcal{S}\mapsto \Gamma $\textit{\ is a
twice continuously differentiable function with bounded derivatives.}
\item \textit{(vii) }$D\left( q,s\right) $\textit{, }$V\left( q,s\right) $
\textit{, and }$f\left( q,s\right) $\textit{\ are uniformly bounded in }$
(q,s)$\textit{, continuous in }$q$\textit{, and twice continuously
differentiable in }$s$ \textit{with bounded derivatives. For any }$i,j\in
\Lambda _{n}$,\textit{\ the joint density of }$(q_{i},q_{j},s_{i},s_{j})^{
\intercal }$\textit{\ is uniformly bounded above by some constant }$C<\infty
$\textit{\ and continuously differentiable in all components.}
\item \textit{(viii) }$c_{0}^{\top }D\left( {\Greekmath 010D} _{0}(s),s\right) c_{0}>0$
\textit{, }$c_{0}^{\top }V\left( {\Greekmath 010D} _{0}(s),s\right) c_{0}>0$\textit{,
and }$f\left( {\Greekmath 010D} _{0}(s),s\right) >0$ \textit{for all }$s\in \mathcal{S}$
\textit{.}
\item \textit{(ix) As }$n\rightarrow \infty $\textit{, }$b_{n}\rightarrow 0$
\textit{, }$n^{1-2{\Greekmath 010F} }b_{n}/\log n\rightarrow \infty $\textit{, }$\log
n/(nb_{n}^{2})\rightarrow 0$, \textit{and }$nb_{n}^{(2+2{\Greekmath 0127} )/(2+{\Greekmath 0127}
)}\rightarrow \infty $\textit{\ for some }${\Greekmath 0127} >0$ \textit{given in (v).}
\item \textit{(x) }$K\left( \cdot \right) $\textit{\ is a positive
second-order kernel, which is Lipschitz, symmetric around zero, and
nonincreasing on }$
\mathbb{R}
^{+}$ \textit{and satisfies }$\int K\left( v\right) dv=1$\textit{, }$\int
v^{2}K\left( v\right) dv<\infty $.
\end{description}
\medskip
We provide some discussions about Assumption A. First, we assume that $q_{i}$
and $s_{i}$ are continuous random variables as in the example in Section \ref
{Section tipping}. This setup also covers the two-dimensional
\textquotedblleft spatial structural break\textquotedblright\ model as a
special case as in the example in Section \ref{Section boundary}. For the
latter case, we denote $n_{1}$ and $n_{2}$ as the numbers of rows
(latitudes) and columns (longitudes) in the grid of pixels, and normalize $q$
and $s$ in the way that $q\in \{1/n_{1},2/n_{1},\ldots ,1\}$ and $s\in
\{1/n_{2},2/n_{2},\ldots ,1\}$. Under similar (and even weaker) regularity
conditions as Assumption A, we can show that the asymptotic results in the
following sections extend to this case once we treat $(q_{i},s_{i})^{\top }$
as independently and uniformly distributed random variables over $\left[ 0,1
\right] ^{2}$. Such similarity is also found in the standard structural
break and the threshold regression models (e.g.,\ Proposition 5 in \cite
{Bai98} and Theorem 1 in \cite{Hansen00a}).
Second, Assumption A is mild and common in the existing literature. In
particular, the condition (i) is the same as in \cite{Bolthausen82} to
define the latent random field. Note that ${\Greekmath 0115} _{0}$ can be any strictly
positive value, and hence we can impose ${\Greekmath 0115} _{0}>1$ without loss of
generality. The condition (ii) adopts the widely used shrinking change size
setup as in \cite{Bai97b}, \cite{Bai98}, and \cite{Hansen00a} to obtain a
simple limiting distribution. In contrast, a constant change size (when $
{\Greekmath 010F} =0$) leads to a complicated asymptotic distribution of the
threshold estimator, which depends on nuisance parameters (e.g., \cite
{Chan93}). The condition (iii) is required to establish the maximal
inequality and uniform convergence in a spatially dependent random field. We
impose a stronger condition than \cite{Jenish09} to obtain the maximal
inequality uniformly over ${\Greekmath 010D} $ and $s$. We could weaken this condition
such that ${\Greekmath 010B} _{k,l}(m)$ decays at a polynomial rate (e.g., ${\Greekmath 010B}
_{k,l}(m)\leq Cm^{-r}$ for some $r>8$ and a constant $C$ as in \cite
{CarbonFrancqTran2007}) if we impose higher moment restrictions in the
condition (v). However, this exponential decay rate simplifies the technical
proofs. The conditions (iv) to (viii) are similar to Assumption 1 of \cite
{Hansen00a}. The condition (ix) imposes restrictions on the bandwidth $b_{n}$
, which now depends on ${\Greekmath 010F} $ and ${\Greekmath 0127} $. The condition (x) holds
for many commonly used kernel functions including the Gaussian and the
uniform kernels.
Third, we assume ${\Greekmath 010D} _{0}\left( \cdot \right) $ to be a function from $
\mathcal{S}$ to $\Gamma $ in Assumption A-(vi), which is not necessarily
one-to-one. For this reason, sample splitting based on $\mathbf{1}\left[
q_{i}\leq {\Greekmath 010D} _{0}\left( s_{i}\right) \right] $ can be different from
that based on $\mathbf{1}\left[ s_{i}\geq \breve{{\Greekmath 010D}}_{0}\left(
q_{i}\right) \right] $ for some function $\breve{{\Greekmath 010D}}_{0}\left( \cdot
\right) $. Instead of restricting ${\Greekmath 010D} _{0}\left( \cdot \right) $ to be
one-to-one in this paper, we presume that one knows which variables should
be respectively assigned as $q_{i}$ and $s_{i}$ from the context.
Alternatively, we can consider a function $g_{0}\left( q,s\right) $ such
that $g_{0}$ is monotonically increasing in $q$ for any $s$. Then, $\mathbf{1
}\left[ q_{i}\leq {\Greekmath 010D} _{0}\left( s_{i}\right) \right] $ is viewed as a
special case of $\mathbf{1}\left[ g_{0}\left( q_{i},s_{i}\right) \leq 0
\right] $ by inverting $g_{0}\left( \cdot ,s\right) =q^{\ast }$, where $
g_{0}\left( q^{\ast },s\right) =0$. We discuss such extension to identify a
threshold contour in Section \ref{Section contour}.
\section{Asymptotic Results\label{Section asymptotics}}
We first obtain the asymptotic properties of $\widehat{{\Greekmath 010D} }\left(
s\right) $. The following theorem derives the pointwise consistency and the
pointwise rate of convergence of $\widehat{{\Greekmath 010D} }\left( s\right) $ at the
interior points of $\mathcal{S}$.
\begin{theorem}
\label{p-roc}For a given $s\in \mathcal{S}_{0}$, under Assumptions ID and A,
$\widehat{{\Greekmath 010D} }\left( s\right) \rightarrow _{p}{\Greekmath 010D} _{0}\left( s\right)
$ as $n\rightarrow \infty $. Furthermore,
\begin{equation}
\widehat{{\Greekmath 010D} }\left( s\right) -{\Greekmath 010D} _{0}\left( s\right) =O_{p}\left(
\frac{1}{n^{1-2{\Greekmath 010F} }b_{n}}\right) \label{proc}
\end{equation}
provided that $n^{1-2{\Greekmath 010F} }b_{n}^{2}$ does not diverge.
\end{theorem}
The pointwise rate of convergence of $\widehat{{\Greekmath 010D} }\left( s\right) $
depends on two parameters, ${\Greekmath 010F} $ and $b_{n}$. It is decreasing in $
{\Greekmath 010F} $ like the parametric (constant) threshold case: a larger ${\Greekmath 010F}
$ reduces the threshold effect ${\Greekmath 010E} _{0}=c_{0}n^{-{\Greekmath 010F} }$ and hence
decreases the effective sampling information on the threshold. Since we
estimate ${\Greekmath 010D} _{0}(\cdot )$ using the kernel estimation method, the rate
of convergence depends on the bandwidth $b_{n}$ as well. As in the standard
kernel estimator case, a smaller bandwidth decreases the effective local
sample size, which reduces the precision of the estimator $\widehat{{\Greekmath 010D} }
\left( s\right) $. Therefore, in order to have a sufficiently fast rate of
convergence, we need to choose $b_{n}$ large enough when the threshold
effect ${\Greekmath 010E} _{0}$ is expected to be small (i.e., when ${\Greekmath 010F} $ is
close to $1/2$).
Unlike the standard kernel estimator, there seems no bias-variance trade-off
in $\widehat{{\Greekmath 010D} }\left( s\right) $ in (\ref{proc}), implying that we
could improve the rate of convergence by choosing a larger bandwidth $b_{n}$
. However, as we can find in Theorem \ref{g-an} below, $b_{n}$ cannot be
chosen too large to result in $n^{1-2{\Greekmath 010F} }b_{n}^{2}\rightarrow \infty $
, under which $n^{1-2{\Greekmath 010F} }b_{n}(\widehat{{\Greekmath 010D} }\left( s\right)
-{\Greekmath 010D} _{0}\left( s\right) )$ is no longer $O_{p}(1)$. Therefore, we can
obtain the optimal bandwidth using the restriction that $n^{1-2{\Greekmath 010F}
}b_{n}^{2}$ does not diverge.
Under this restriction, we find the optimal bandwidth as $b_{n}^{\ast
}=n^{-(1-2{\Greekmath 010F} )/2}c^{\ast }$ for some constant $0<c^{\ast }<\infty $,
which yields the optimal pointwise rate of convergence of $\widehat{{\Greekmath 010D} }
\left( s\right) $ as $n^{-(1-2{\Greekmath 010F} )/2}$. However, such a bandwidth
choice is not feasible because of the unknown constant $c^{\ast }$ and the
nuisance parameter ${\Greekmath 010F} $ that are not estimable. In practice, we
suggest cross validation as we implement in Section \ref{Section empirics},
although its statistical properties need to be studied further. Note that,
when the change size ${\Greekmath 010E} _{0}$ shrinks very slowly with $n$ (i.e., $
{\Greekmath 010F} $ is close to $0$), the optimal rate of convergence of $\widehat{
{\Greekmath 010D} }(\cdot )$ is close to $n^{-1/2}$. This $\sqrt{n}$-rate is obtained
in the standard kernel regression if the unknown function is infinitely
differentiable, while we only require the second-order differentiability of $
{\Greekmath 010D} _{0}\left( \cdot \right) $.
The next theorem derives the limiting distribution of $\widehat{{\Greekmath 010D} }
\left( s\right) $. We let $W(\cdot )$ be a two-sided Brownian motion defined
as in \cite{Hansen00a}:
\begin{equation}
W(r)=W_{1}(-r)\mathbf{1}\left[ r<0\right] +W_{2}(r)\mathbf{1}\left[ r>0
\right] \relax\protect\ifmmode\expandafter\text@\else\expandafter\mbox\fi{,} \label{2BM}
\end{equation}
where $W_{1}(\cdot )$ and $W_{2}(\cdot )$ are independent standard Brownian
motions on $[0,\infty )$.
\begin{theorem}
\label{g-an}Under Assumptions ID and A, for a given $s\in \mathcal{S}_{0}$,
if $n^{1-2{\Greekmath 010F} }b_{n}^{2}\rightarrow {\Greekmath 0125} \in (0,\infty )$,
\begin{equation}
n^{1-2{\Greekmath 010F} }b_{n}\left( \widehat{{\Greekmath 010D} }\left( s\right) -{\Greekmath 010D}
_{0}\left( s\right) \right) \rightarrow _{d}{\Greekmath 0118} \left( s\right) \arg
\max_{r\in
\mathbb{R}
}\left( W\left( r\right) +{\Greekmath 0116} \left( r,{\Greekmath 0125} ;s\right) \right)
\label{a-normal}
\end{equation}
as $n\rightarrow \infty $, where
\begin{eqnarray*}
{\Greekmath 0116} \left( r,{\Greekmath 0125} ;s\right) &=&-\left\vert r\right\vert {\Greekmath 0120} _{0}\left(
r,{\Greekmath 0125} ;s\right) +\frac{{\Greekmath 0125} |\dot{{\Greekmath 010D}}_{0}(s)|}{{\Greekmath 0118} (s)}{\Greekmath 0120}
_{1}\left( r,{\Greekmath 0125} ;s\right) \relax\protect\ifmmode\expandafter\text@\else\expandafter\mbox\fi{,} \\
{\Greekmath 0120} _{j}\left( r,{\Greekmath 0125} ;s\right) &=&\int_{0}^{|r|{\Greekmath 0118} (s)/\left( {\Greekmath 0125} |
\dot{{\Greekmath 010D}}_{0}(s)|\right) }t^{j}K\left( t\right) dt\ \ \ \relax\protect\ifmmode\expandafter\text@\else\expandafter\mbox\fi{for }j=0,1
\relax\protect\ifmmode\expandafter\text@\else\expandafter\mbox\fi{,} \\
{\Greekmath 0118} \left( s\right) &=&\frac{{\Greekmath 0114} _{2}c_{0}^{\top }V\left( {\Greekmath 010D}
_{0}\left( s\right) ,s\right) c_{0}}{\left( c_{0}^{\top }D\left( {\Greekmath 010D}
_{0}\left( s\right) ,s\right) c_{0}\right) ^{2}f\left( {\Greekmath 010D} _{0}\left(
s\right) ,s\right) }
\end{eqnarray*}
with ${\Greekmath 0114} _{2}=\int K(v)^{2}dv$ and $\dot{{\Greekmath 010D}}_{0}\left( s\right) $ is
the first derivative of ${\Greekmath 010D} _{0}$ at $s$. Furthermore,
\begin{equation*}
\mathbb{E}\left[ \arg \max_{r\in
\mathbb{R}
}\left( W\left( r\right) +{\Greekmath 0116} \left( r,{\Greekmath 0125} ;s\right) \right) \right] =0.
\end{equation*}
\end{theorem}
The drift term ${\Greekmath 0116} \left( r,{\Greekmath 0125} ;s\right) $ in (\ref{a-normal}) depends
on the constant $0<{\Greekmath 0125} <\infty $, which is the limit of $n^{1-2{\Greekmath 010F}
}b_{n}^{2}=(n^{1-2{\Greekmath 010F} }b_{n})b_{n}$, and $|\dot{{\Greekmath 010D}}_{0}(s)|$, the
steepness of ${\Greekmath 010D} _{0}(\cdot )$ at $s$. Interestingly, it resembles the
typical $O(b_{n})$ boundary bias of the standard local constant estimator.
However, this non-zero drift term is not because of the typical boundary
effect but because of the inequality restriction inside the indicator
function, $\mathbf{1}\left[ q_{i}\leq {\Greekmath 010D} _{0}\left( s_{i}\right) \right]
$, which characterizes the sample splitting.
It is important to note that having this non-zero drift term in the limiting
expression does not mean that the limiting distribution of $\widehat{{\Greekmath 010D} }
\left( s\right) $ has a non-zero mean, even when we use the optimal
bandwidth $b_{n}^{\ast }=O(n^{-(1-2{\Greekmath 010F} )/2})$ satisfying $
n^{1-2{\Greekmath 010F} }b_{n}^{\ast 2}\rightarrow {\Greekmath 0125} \in (0,\infty )$. This is
mainly because the drift function ${\Greekmath 0116} \left( r,{\Greekmath 0125} ;s\right) $ is
symmetric about zero and hence the limiting random variable $\arg \max_{r\in
\mathbb{R}
}\left( W\left( r\right) +{\Greekmath 0116} \left( r,{\Greekmath 0125} ;s\right) \right) $ is mean
zero. In general, we can show that the random variable $\arg \max_{r\in
\mathbb{R}
}\left( W\left( r\right) +{\Greekmath 0116} \left( r,{\Greekmath 0125} ;s\right) \right) $ always
has mean zero if ${\Greekmath 0116} \left( r,{\Greekmath 0125} ;s\right) $ is a non-random function
that is symmetric about zero and monotonically decreasing fast enough. This
result might be of independent research interest and is summarized in Lemma
\ref{lemma-drift} in the Appendix. Figure \ref{fig drift} depicts the drift
function ${\Greekmath 0116} \left( r,{\Greekmath 0125} ;s\right) $ for various kernels when ${\Greekmath 0118}
(s)/\left( {\Greekmath 0125} |\dot{{\Greekmath 010D}}_{0}(s)|\right) =1$.
\begin{figure}[h]
\centering
\caption{Drift function ${\Greekmath 0116} \left( r,{\Greekmath 0125} ;s\right) $ for different kernels (color online)}
\includegraphics[width=0.6\textwidth]{F10523.jpg}
\label{fig drift}
\end{figure}
Since the limiting distribution in (\ref{a-normal}) depends on unknown
components, like ${\Greekmath 0125} $ and $\dot{{\Greekmath 010D}}_{0}(s)$, it is hard to use
this result for further inference. We instead suggest undersmoothing for
practical use. More precisely, if we suppose $n^{1-2{\Greekmath 010F}
}b_{n}^{2}\rightarrow 0$ as $n\rightarrow \infty $, then the limiting
distribution in (\ref{a-normal}) simplifies to\footnote{
We let ${\Greekmath 0120} _{0}\left( r,0;s\right) =\int_{0}^{\infty }K\left( t\right)
dt=1/2$ and ${\Greekmath 0120} _{1}\left( r,0;s\right) =\int_{0}^{\infty }tK\left(
t\right) dt<\infty $.}
\begin{equation}
n^{1-2{\Greekmath 010F} }b_{n}\left( \widehat{{\Greekmath 010D} }\left( s\right) -{\Greekmath 010D}
_{0}\left( s\right) \right) \rightarrow _{d}{\Greekmath 0118} \left( s\right) \arg
\max_{r\in
\mathbb{R}
}\left( W\left( r\right) -\frac{\left\vert r\right\vert }{2}\right)
\label{an0}
\end{equation}
as $n\rightarrow \infty $, which appears the same as in the parametric case
in \cite{Hansen00a} except for the scaling factor $n^{1-2{\Greekmath 010F} }b_{n}$.
The distribution of $\arg \max_{r\in
\mathbb{R}
}\left( W\left( r\right) -\left\vert r\right\vert /2\right) $ is known
(e.g., \cite{Bhattacharya76} and \cite{Bai97b}), which is also described in
Hansen (2000, p.581). The ${\Greekmath 0118} \left( s\right) $ term determines the scale
of the distribution at given $s$ in the way that it increases in the
conditional variance $\mathbb{E}[u_{i}^{2}|x_{i},q_{i},s_{i}]$ and decreases
in the size of the threshold constant $c_{0}$ and the density of $
(q_{i},s_{i})$ near the threshold.
Even when $n^{1-2{\Greekmath 010F} }b_{n}^{2}\rightarrow 0$ as $n\rightarrow \infty $
, the asymptotic distribution in (\ref{an0}) still depends on the unknown
parameter ${\Greekmath 010F} $ (or equivalently $c_{0}$) in ${\Greekmath 0118} \left( s\right) $
that is not estimable. Thus, this result cannot be directly used for
inference of ${\Greekmath 010D} _{0}\left( s\right) $. Alternatively, given any $s\in
\mathcal{S}_{0}$, we can consider a pointwise likelihood ratio test
statistic for
\begin{equation}
H_{0}:{\Greekmath 010D} _{0}\left( s\right) ={\Greekmath 010D} _{\ast }\left( s\right) \relax\protect\ifmmode\expandafter\text@\else\expandafter\mbox\fi{ \ \
against \ \ }H_{1}:{\Greekmath 010D} _{0}\left( s\right) \neq {\Greekmath 010D} _{\ast }\left(
s\right) , \label{hypo}
\end{equation}
\ which is given as
\begin{equation}
LR_{n}(s)=\mathop{\displaystyle \sum }_{i\in \Lambda _{n}}K\left( \frac{s_{i}-s}{b_{n}}\right)
\frac{Q_{n}\left( {\Greekmath 010D} _{\ast }\left( s\right) ,s\right) -Q_{n}\left(
\widehat{{\Greekmath 010D} }\left( s\right) ,s\right) }{Q_{n}\left( \widehat{{\Greekmath 010D} }
\left( s\right) ,s\right) }\relax\protect\ifmmode\expandafter\text@\else\expandafter\mbox\fi{.} \label{LR}
\end{equation}
The following corollary obtains the limiting null distribution of this test
statistic that is free of nuisance parameters. By inverting the likelihood
ratio statistic, we can form a pointwise confidence interval for ${\Greekmath 010D}
_{0}\left( s\right) $.
\begin{corollary}
\label{g-test}Suppose $n^{1-2{\Greekmath 010F} }b_{n}^{2}\rightarrow 0$ as $
n\rightarrow \infty $. Under the same condition in Theorem \ref{g-an}, for
any fixed $s\in \mathcal{S}_{0}$, the test statistic in (\ref{LR}) satisfies
\begin{equation}
LR_{n}(s)\rightarrow _{d}{\Greekmath 0118} _{LR}\left( s\right) \max_{r\in
\mathbb{R}
}\left( 2W\left( r\right) -\left\vert r\right\vert \right) \label{limLR}
\end{equation}
as $n\rightarrow \infty $ under the null hypothesis (\ref{hypo}), where
\begin{equation*}
{\Greekmath 0118} _{LR}\left( s\right) =\frac{{\Greekmath 0114} _{2}c_{0}^{\top }V\left( {\Greekmath 010D}
_{0}\left( s\right) ,s\right) c_{0}}{{\Greekmath 011B} ^{2}(s)c_{0}^{\top }D\left(
{\Greekmath 010D} _{0}\left( s\right) ,s\right) c_{0}}
\end{equation*}
with ${\Greekmath 011B} ^{2}(s)=\mathbb{E}\left[ u_{i}^{2}|s_{i}=s\right] $ and ${\Greekmath 0114}
_{2}=\int K(v)^{2}dv$.
\end{corollary}
When $\mathbb{E}[u_{i}^{2}|x_{i},q_{i},s_{i}]=\mathbb{E}[u_{i}^{2}|s_{i}]$,
which is the case of local conditional homoskedasticity, the scale parameter
${\Greekmath 0118} _{LR}\left( s\right) $ is simplified as ${\Greekmath 0114} _{2}$, and hence the
limiting null distribution of $LR_{n}(s)$ becomes free of nuisance
parameters and the same for all $s\in \mathcal{S}_{0}$. Though this limiting
distribution is still nonstandard, the critical values in this case can be
simulated using the same method as Hansen (2000, p.582) with the scale
adjusted by ${\Greekmath 0114} _{2}$. More precisely, since the distribution function
of ${\Greekmath 0110} =\max_{r\in
\mathbb{R}
}\left( 2W\left( r\right) -\left\vert r\right\vert \right) $ is given as $
\mathbb{P}({\Greekmath 0110} \leq z)=(1-\exp (-z/2))^{2}\mathbf{1}\left[ z\geq 0\right] $
, the distribution function of ${\Greekmath 0110} ^{\ast }={\Greekmath 0114} _{2}{\Greekmath 0110} $ is $
\mathbb{P}({\Greekmath 0110} ^{\ast }\leq z)=(1-\exp (-z/2{\Greekmath 0114} _{2}))^{2}\mathbf{1}
\left[ z\geq 0\right] $, where ${\Greekmath 0110} ^{\ast }$ is the limiting random
variable of $LR_{n}(s)$ given in (\ref{limLR}) under the local conditional
homoskedasticity. By inverting it, we can obtain the critical values given a
choice of $K(\cdot )$. For instance, the critical values for the Gaussian
kernel is reported in Table \ref{tbl cv}, where ${\Greekmath 0114} _{2}=(2\sqrt{{\Greekmath 0119} }
)^{-1}\simeq 0.2821$ in this case.
\begin{table}[tbp]
\centering
\caption{Simulated Critical Values of the LR Test (Gaussian Kernel)}\label
{tbl cv}
\bigskip
\begin{tabular}{crrrrrrrr}
\hline\hline
$\mathbb{P}({\Greekmath 0110} ^{\ast }>cv)$ & & {\small 0.800} & {\small 0.850} &
{\small 0.900} & {\small 0.925} & {\small 0.950} & {\small 0.975} & {\small
0.990} \\ \hline
$cv$ & & {\small 1.268} & {\small 1.439} & {\small 1.675} & {\small 1.842}
& {\small 2.074} & {\small 2.469} & {\small 2.988} \\ \hline
\end{tabular}
\bigskip
\raggedright{\footnotesize Note: ${\Greekmath 0110} ^{\ast }$ is the limiting
distribution of $LR_{n}(s)$ under the local conditional homoskedasticity.
The Gaussian kernel is used.}\bigskip
\end{table}
In general, we can estimate ${\Greekmath 0118} _{LR}\left( s\right) $ by
\begin{equation*}
\widehat{{\Greekmath 0118} }_{LR}\left( s\right) =\frac{{\Greekmath 0114} _{2}\widehat{{\Greekmath 010E} }^{\top
}\widehat{V}\left( \widehat{{\Greekmath 010D} }\left( s\right) ,s\right) \widehat{
{\Greekmath 010E} }}{\widehat{{\Greekmath 011B} }^{2}(s)\widehat{{\Greekmath 010E} }^{\top }\widehat{D}\left(
\widehat{{\Greekmath 010D} }\left( s\right) ,s\right) \widehat{{\Greekmath 010E} }},
\end{equation*}
where $\widehat{{\Greekmath 010E} }$ is from (\ref{para-b}) and (\ref{para-d}), and $
\widehat{{\Greekmath 011B} }^{2}(s)$, $\widehat{D}\left( \widehat{{\Greekmath 010D} }\left(
s\right) ,s\right) $, and $\widehat{V}\left( \widehat{{\Greekmath 010D} }\left(
s\right) ,s\right) $ are the standard Nadaraya-Watson estimators. In
particular, we let $\widehat{{\Greekmath 011B} }^{2}(s)=\sum_{i\in \Lambda _{n}}{\Greekmath 0121}
_{1i}(s)\widehat{u}_{i}^{2}$ with $\widehat{u}_{i}=y_{i}-x_{i}^{\top }
\widehat{{\Greekmath 010C} }-x_{i}^{\top }\widehat{{\Greekmath 010E} }\mathbf{1}\left[ q_{i}\leq
\widehat{{\Greekmath 010D} }\left( s_{i}\right) \right] $,
\begin{equation*}
\widehat{D}\left( \widehat{{\Greekmath 010D} }\left( s\right) ,s\right) =\mathop{\displaystyle \sum }_{i\in
\Lambda _{n}}{\Greekmath 0121} _{2i}(s)x_{i}x_{i}^{\top }\relax\protect\ifmmode\expandafter\text@\else\expandafter\mbox\fi{, \ and }\ \widehat{V}
\left( \widehat{{\Greekmath 010D} }\left( s\right) ,s\right) =\mathop{\displaystyle \sum }_{i\in \Lambda
_{n}}{\Greekmath 0121} _{2i}(s)x_{i}x_{i}^{\top }\widehat{u}_{i}^{2}\relax\protect\ifmmode\expandafter\text@\else\expandafter\mbox\fi{,}
\end{equation*}
where
\begin{equation*}
{\Greekmath 0121} _{1i}(s)=\frac{K\left( (s_{i}-s)/b_{n}\right) }{\sum_{j\in \Lambda
_{n}}K\left( (s_{j}-s)/b_{n}\right) }\relax\protect\ifmmode\expandafter\text@\else\expandafter\mbox\fi{ \ and \ }{\Greekmath 0121} _{2i}(s)=\frac{
\mathbb{K}\left( (q_{i}-\widehat{{\Greekmath 010D} }\left( s\right) )/b_{n}^{\prime
},(s_{i}-s)/b_{n}^{\prime \prime }\right) }{\sum_{j\in \Lambda _{n}}\mathbb{K
}\left( (q_{j}-\widehat{{\Greekmath 010D} }\left( s\right) )/b_{n}^{\prime
},(s_{j}-s)/b_{n}^{\prime \prime }\right) }
\end{equation*}
for some bivariate kernel function $\mathbb{K}(\cdot ,\cdot )$ and bandwidth
parameters $(b_{n}^{\prime },b_{n}^{\prime \prime })$.
Finally, we show the $\sqrt{n}$-consistency of the semiparametric estimators
$\widehat{{\Greekmath 010C} }$ and $\widehat{{\Greekmath 010E} }^{\ast }$ in (\ref{para-b}) and (
\ref{para-d}). For this purpose, we first obtain the uniform rate of
convergence of $\widehat{{\Greekmath 010D} }\left( s\right) $.
\begin{theorem}
\label{u-roc}Under Assumptions ID and A,
\begin{equation*}
\sup_{s\in \mathcal{S}_{0}}\left\vert \widehat{{\Greekmath 010D} }\left( s\right)
-{\Greekmath 010D} _{0}\left( s\right) \right\vert =O_{p}\left( \frac{\log n}{
n^{1-2{\Greekmath 010F} }b_{n}}\right)
\end{equation*}
provided that $n^{1-2{\Greekmath 010F} }b_{n}^{2}$ does not diverge.
\end{theorem}
\noindent Apparently, the uniform consistency of $\widehat{{\Greekmath 010D} }\left(
s\right) $ follows when $\log n/(n^{1-2{\Greekmath 010F} }b_{n})\rightarrow 0$ as $
n\rightarrow \infty $. Based on this uniform convergence, the following
theorem derives the joint limiting distribution of $\widehat{{\Greekmath 010C} }$ and $
\widehat{{\Greekmath 010E} }^{\ast }$. We let $\widehat{{\Greekmath 0112} }^{\ast }=(\widehat{
{\Greekmath 010C} }^{\top },\widehat{{\Greekmath 010E} }^{\ast \top })^{\top }$ and ${\Greekmath 0112}
_{0}^{\ast }=({\Greekmath 010C} _{0}^{\top },{\Greekmath 010E} _{0}^{\ast \top })^{\top }$.
\begin{theorem}
\label{bd}Suppose the conditions in Theorem \ref{u-roc} hold. If we let ${\Greekmath 0119}
_{n}>0$ such that ${\Greekmath 0119} _{n}\rightarrow 0$ and $\{\log n/(n^{1-2{\Greekmath 010F}
}b_{n})\}/{\Greekmath 0119} _{n}\rightarrow 0$ as $n\rightarrow \infty $, we have
\begin{equation}
\sqrt{n}\left( \widehat{{\Greekmath 0112} }^{\ast }-{\Greekmath 0112} _{0}^{\ast }\right)
\rightarrow _{d}\mathcal{N}\left( 0,\Sigma _{X}^{\ast -1}\Omega ^{\ast
}\Sigma _{X}^{\ast -1}\right) \label{th-an*}
\end{equation}
as $n\rightarrow \infty $, where
\begin{equation*}
\Sigma _{X}^{\ast }=\left[
\begin{array}{cc}
\mathbb{E}\left[ x_{i}x_{i}^{\top }\mathbf{1}_{i}^{+}\right] & 0 \\
0 & \mathbb{E}\left[ x_{i}x_{i}^{\top }\mathbf{1}_{i}^{-}\right]
\end{array}
\right] \relax\protect\ifmmode\expandafter\text@\else\expandafter\mbox\fi{ \ and \ }\Omega ^{\ast }=\lim_{n\rightarrow \infty }\frac{1}{n
}Var\left[
\begin{array}{c}
\sum_{i\in \Lambda _{n}}x_{i}u_{i}\mathbf{1}_{i}^{+} \\
\sum_{i\in \Lambda _{n}}x_{i}u_{i}\mathbf{1}_{i}^{-}
\end{array}
\right]
\end{equation*}
with $\mathbf{1}_{i}^{+}=\mathbf{1[}q_{i}>{\Greekmath 010D} _{0}(s_{i})]\mathbf{1}
[s_{i}\in \mathcal{S}_{0}]$ and $\mathbf{1}_{i}^{-}=\mathbf{1[}q_{i}<{\Greekmath 010D}
_{0}(s_{i})]\mathbf{1}[s_{i}\in \mathcal{S}_{0}]$.
\end{theorem}
For the second-step estimator $\widehat{{\Greekmath 0112} }^{\ast }$, we use (\ref
{para-b}) and (\ref{para-d}), instead of the conventional plug-in
estimation, say $\arg \min_{{\Greekmath 010C} ,{\Greekmath 010E} }\sum_{i\in \Lambda
_{n}}(y_{i}-x_{i}^{\top }{\Greekmath 010C} -x_{i}^{\top }{\Greekmath 010E} \mathbf{1}[q_{i}\leq
\widehat{{\Greekmath 010D} }\left( s_{i}\right) ])^{2}\mathbf{1}[s_{i}\in \mathcal{S}
_{0}]$. The reason is that the first-step nonparametric estimator $\widehat{
{\Greekmath 010D} }(\cdot )$ may not be asymptotically orthogonal to the second step.
Unlike the standard semiparametric literature (e.g.,\ Assumption N(c) in
\cite{Andrews94a}), the asymptotic effect of $\widehat{{\Greekmath 010D} }\left(
s\right) $ to the second-step estimation is not easily derived due to the
discontinuity. The new estimation idea above, however, only uses the
observations that are little affected by the estimation error in the first
step to achieve asymptotic orthogonality. As we verify in Lemma \ref{bias1}
in the Appendix, this is done by choosing a large enough ${\Greekmath 0119} _{n}$ in (\ref
{para-b}) and (\ref{para-d}) such that the observations that are included in
the second step are outside the uniform convergence bound of $\left\vert
\widehat{{\Greekmath 010D} }\left( s\right) -{\Greekmath 010D} _{0}\left( s\right) \right\vert $.
Thanks to the threshold regression structure, we can estimate the parameters
on each side of the threshold even using these subsamples. Meanwhile, we
also want ${\Greekmath 0119} _{n}\rightarrow 0$ fast enough to include more observations.
By doing so, though we lose some efficiency in finite samples, we can derive
the asymptotic normality of $\widehat{{\Greekmath 0112} }=(\widehat{{\Greekmath 010C} }^{\top },
\widehat{{\Greekmath 010E} }^{\top })^{\top }$ that has zero mean and achieves the same
asymptotic variance as if ${\Greekmath 010D} _{0}(\cdot )$ was known.
By the delta method, Theorem \ref{bd} readily yields the limiting
distribution of $\widehat{{\Greekmath 0112} }=(\widehat{{\Greekmath 010C} }^{\top },\widehat{{\Greekmath 010E}
}^{\top })^{\top }$ as
\begin{equation}
\sqrt{n}\left( \widehat{{\Greekmath 0112} }-{\Greekmath 0112} _{0}\right) \rightarrow _{d}\mathcal{
N}\left( 0,\Sigma _{X}^{-1}\Omega \Sigma _{X}^{-1}\right) \relax\protect\ifmmode\expandafter\text@\else\expandafter\mbox\fi{ \ as }
n\rightarrow \infty \relax\protect\ifmmode\expandafter\text@\else\expandafter\mbox\fi{,} \label{th-an}
\end{equation}
where
\begin{equation*}
\Sigma _{X}=\mathbb{E}\left[ z_{i}z_{i}^{\top }\mathbf{1}\left[ s_{i}\in
\mathcal{S}_{0}\right] \right] \relax\protect\ifmmode\expandafter\text@\else\expandafter\mbox\fi{ \ and \ }\Omega =\lim_{n\rightarrow
\infty }\frac{1}{n}Var\left[ \sum_{i\in \Lambda _{n}}z_{i}u_{i}\mathbf{1}
\left[ s_{i}\in \mathcal{S}_{0}\right] \right]
\end{equation*}
with $z_{i}=[x_{i}^{\top },x_{i}^{\top }\mathbf{1}\left[ q_{i}\leq {\Greekmath 010D}
_{0}\left( s_{i}\right) \right] ]^{\top }$. The asymptotic variance
expressions in (\ref{th-an*}) and (\ref{th-an}) allow for cross-sectional
dependence as they have the long-run variance (LRV) forms $\Omega ^{\ast }$
and $\Omega $. They can be consistently estimated by the robust estimator
developed by \cite{Conley07} using $\widehat{u}_{i}=(y_{i}-x_{i}^{\top }
\widehat{{\Greekmath 010C} }-x_{i}^{\top }\widehat{{\Greekmath 010E} }\mathbf{1}[q_{i}\leq \widehat{
{\Greekmath 010D} }\left( s_{i}\right) ])\mathbf{1}[s_{i}\in \mathcal{S}_{0}]$. The
terms $\Sigma _{X}^{\ast }$ and $\Sigma _{X}$ can be estimated by their
sample analogues.
\section{Threshold Contour\label{Section contour}}
When we consider sample splitting over a two-dimensional space (i.e., $q_{i}$
and $s_{i}$ respectively correspond to the latitude and longitude on the
map), the threshold model (\ref{model}) can be generalized to estimate a
nonparametric contour threshold model:
\begin{equation}
y_{i}=x_{i}^{\top }{\Greekmath 010C} _{0}+x_{i}^{\top }{\Greekmath 010E} _{0}\mathbf{1}\left[
g_{0}\left( q_{i},s_{i}\right) \leq 0\right] +u_{i}\relax\protect\ifmmode\expandafter\text@\else\expandafter\mbox\fi{,} \label{model2}
\end{equation}
where the unknown function\ $g_{0}:\mathcal{Q}\times \mathcal{S}\mapsto
\mathbb{R}$ determines the threshold contour on a random field that yields
sample splitting. An interesting example includes identifying an unknown
closed boundary over the map, such as a city boundary, and an area of a
disease outbreak or airborne pollution. In social science, it can identify a
group boundary or a region in which the agents share common demographic,
political, or economic characteristics.
To relate this generalized form to the original threshold model (\ref{model}
), we suppose there exists a known center at $\left( q_{i}^{\ast
},s_{i}^{\ast }\right) $ such that $g_{0}\left( q_{i}^{\ast },s_{i}^{\ast
}\right) <0$. Without loss of generality, we can normalize $\left(
q_{i}^{\ast },s_{i}^{\ast }\right) $ to be $\left( 0,0\right) $ and
re-center the original location variables $(q_{i},s_{i})$ accordingly. In
addition, we define the radius distance $l_{i}$ and angle $a_{i}^{\circ }$
of the $i$th observation relative to the origin as
\begin{eqnarray*}
l_{i} &=&(q_{i}^{2}+s_{i}^{2})^{1/2}\relax\protect\ifmmode\expandafter\text@\else\expandafter\mbox\fi{,} \\
a_{i}^{\circ } &=&\bar{a}_{i}^{\circ }\mathbf{I}_{i}+\left( 180^{\circ }-
\bar{a}_{i}^{\circ }\right) \mathbf{II}_{i}+\left( 180^{\circ }+\bar{a}
_{i}^{\circ }\right) \mathbf{III}_{i}+\left( 360^{\circ }-\bar{a}_{i}^{\circ
}\right) \mathbf{IV}_{i}\relax\protect\ifmmode\expandafter\text@\else\expandafter\mbox\fi{,}
\end{eqnarray*}
where $\bar{a}_{i}^{\circ }=\arctan \left( \left\vert q_{i}/s_{i}\right\vert
\right) $, and each of $(\mathbf{I}_{i},\mathbf{II}_{i},\mathbf{III}_{i},
\mathbf{IV}_{i})$ respectively denotes the indicator that the $i$th
observation locates in the first, second, third, and forth quadrant.
We suppose that there is only one threshold at any angle and the threshold
contour is star-shaped. For each chosen angle $a^{\circ }\in \lbrack
0^{\circ },360^{\circ })$, we rotate the original coordinate
counterclockwise and implement the least squares estimation (\ref{reg}) only
using the observations in the first two quadrants after rotation. It will
ensure that the threshold mapping after rotation is a well-defined function.
\begin{figure}[tbp]
\centering
\caption{Illustration of rotation (color online)}\label{fig contour}
\includegraphics[width=1\textwidth]{F20523.jpg}
\end{figure}
In particular, the angle relative to the origin is $a_{i}^{\circ }-a^{\circ
} $ after rotating the coordinate by $a^{\circ }$ degrees counterclockwise,
and the new location (after the rotation) is given as $(q_{i}\left( a^{\circ
}\right) ,s_{i}\left( a^{\circ }\right) )$, where
\begin{equation*}
\left(
\begin{array}{c}
q_{i}\left( a^{\circ }\right) \\
s_{i}\left( a^{\circ }\right)
\end{array}
\right) =\left(
\begin{array}{c}
q_{i}\cos \left( a^{\circ }\right) -s_{i}\sin \left( a^{\circ }\right) \\
s_{i}\cos \left( a^{\circ }\right) +q_{i}\sin \left( a^{\circ }\right)
\end{array}
\right) \relax\protect\ifmmode\expandafter\text@\else\expandafter\mbox\fi{.}
\end{equation*}
After this rotation, we estimate the following nonparametric threshold model:
\begin{equation}
y_{i}=x_{i}^{\top }{\Greekmath 010C} _{0}+x_{i}^{\top }{\Greekmath 010E} _{0}\mathbf{1}\left[
q_{i}\left( a^{\circ }\right) \leq {\Greekmath 010D} _{a^{\circ }}\left( s_{i}\left(
a^{\circ }\right) \right) \right] +u_{i} \label{model3}
\end{equation}
using only the observations $i$ satisfying $q_{i}\left( a^{\circ }\right)
\geq 0$ and in the neighborhood of $s_{i}\left( a^{\circ }\right) =0$, where
${\Greekmath 010D} _{a^{\circ }}\left( \cdot \right) $ is the unknown threshold curve
as in the original model (\ref{model}) on the $a^{\circ }$-degree-rotated
coordinate plane. Such reparametrization guarantees that ${\Greekmath 010D} _{a^{\circ
}}\left( \cdot \right) $ is always positive and it is estimated at the
origin. Figure \ref{fig contour} illustrates the idea of such rotation and
pointwise estimation over a bounded support so that only the red cross
points are included for estimation at different angles. Thus, the estimation
and inference procedures developed in the previous sections are directly
applicable, though we expect some efficiency loss as we only use the
subsample with $q_{i}\left( a^{\circ }\right) \geq 0$ at each $a^{\circ }$.
\section{Monte Carlo Experiments\label{Section simulation}}
We examine the small sample performance of the semiparametric threshold
regression estimator by Monte Carlo simulations. We generate $n$ draws from
\begin{equation}
y_{i}=x_{i}^{\top }{\Greekmath 010C} _{0}+x_{i}^{\top }{\Greekmath 010E} _{0}\mathbf{1}\left[
q_{i}\leq {\Greekmath 010D} _{0}\left( s_{i}\right) \right] +u_{i}\relax\protect\ifmmode\expandafter\text@\else\expandafter\mbox\fi{,}
\label{sim-model}
\end{equation}
where $x_{i}=(1,x_{2i})^{\top }$ and $x_{2i}\in \mathbb{R}$. We let ${\Greekmath 010C}
_{0}=({\Greekmath 010C} _{10},{\Greekmath 010C} _{20})^{\top }=0{\Greekmath 0113} _{2}$ and consider three
different values of ${\Greekmath 010E} _{0}=({\Greekmath 010E} _{10},{\Greekmath 010E} _{20})^{\top }={\Greekmath 010E}
{\Greekmath 0113} _{2}$ with ${\Greekmath 010E} =1,2,3,4$, where ${\Greekmath 0113} _{2}=(1,1)^{\top }$. For
the threshold function, we let ${\Greekmath 010D} _{0}\left( s\right) =\sin (s)/2$. We
consider the cross-sectional dependence structure in $\left(
x_{2i},q_{i},s_{i},u_{i}\right) ^{\top }$ as follows:
\begin{equation}
\left\{
\begin{array}{l}
\left( q_{i},s_{i}\right) ^{\top }\sim iid\mathcal{N}\left( 0,I_{2}\right)
\relax\protect\ifmmode\expandafter\text@\else\expandafter\mbox\fi{;} \\
x_{2i}|\left( q_{i},s_{i}\right) \sim iid\mathcal{N}\left( 0,(1+{\Greekmath 011A} \left(
s_{i}^{2}+q_{i}^{2}\right) )^{-1}\right) \relax\protect\ifmmode\expandafter\text@\else\expandafter\mbox\fi{;} \\
\underline{\mathbf{u}}|\{(x_{i},q_{i},s_{i})\}_{i=1}^{n}\sim \mathcal{N}
\left( 0,\Sigma \right) \relax\protect\ifmmode\expandafter\text@\else\expandafter\mbox\fi{,}
\end{array}
\right. \label{DGP corr}
\end{equation}
where $\underline{\mathbf{u}}=(u_{1},\ldots ,u_{n})^{\top }$. The $(i,j)$th
element of $\Sigma $ is $\Sigma _{ij}={\Greekmath 011A} ^{\lfloor \ell _{ij}n\rfloor }
\mathbf{1}[\ell _{ij}<m/n]$, where $\ell _{ij}=\{\left( s_{i}-s_{j}\right)
^{2}+\left( q_{i}-q_{j}\right) ^{2}\}^{1/2}$ is the $L^{2}$-distance between
the $i$th and $j$th observations. The diagonal elements of $\Sigma $ are
normalized as $\Sigma _{ii}=1$. This $m$-dependent setup follows from the
Monte Carlo experiment in \cite{Conley07} in the sense that each unit can be
cross-sectionally correlated with at most $2m^{2}$ observations. Within the $
m$ distance, the dependence decays at a polynomial rate as indicated by $
{\Greekmath 011A} ^{\lfloor \ell _{ij}n\rfloor }$. The parameter ${\Greekmath 011A} $ describes the
strength of cross-sectional dependence in the way that a larger ${\Greekmath 011A} $
leads to stronger dependence relative to the unit standard deviation. In
particular, we consider the cases with ${\Greekmath 011A} =0$ (i.e., i.i.d.
observations), $0.5$, and $1$. We consider the sample size $n=100$, $200$,
and $500$, and set $\mathcal{S}_{0}$ to include the middle 70\% observations
of $s_{i}$.
\begin{table}[tbp]
\centering
\caption{Rej. Prob. of the LR Test with i.i.d. Data}\label{tbl r1}
\bigskip
\begin{tabular}{cccccccccccccccc}
\hline\hline
& & \multicolumn{4}{c}{${\small s=0.0}$} & & \multicolumn{4}{c}{${\small
s=0.5}$} & & \multicolumn{4}{c}{${\small s=1.0}$} \\
\cline{3-6}\cline{8-11}\cline{13-16}
${\small n}$ & \multicolumn{1}{l}{${\small {\Greekmath 010E} =}$} & {\small 1} &
{\small 2} & {\small 3} & {\small 4} & & {\small 1} & {\small 2} & {\small 3
} & {\small 4} & & {\small 1} & {\small 2} & {\small 3} & {\small 4} \\
\hline
\multicolumn{1}{r}{\small 100} & & {\small 0.16} & {\small 0.09} & {\small
0.06} & {\small 0.08} & & {\small 0.19} & {\small 0.10} & {\small 0.09} &
{\small 0.08} & & {\small 0.25} & {\small 0.18} & {\small 0.16} & {\small
0.12} \\
\multicolumn{1}{r}{\small 200} & & {\small 0.09} & {\small 0.06} & {\small
0.06} & {\small 0.07} & & {\small 0.12} & {\small 0.06} & {\small 0.04} &
{\small 0.06} & & {\small 0.18} & {\small 0.09} & {\small 0.08} & {\small
0.06} \\
\multicolumn{1}{r}{\small 500} & & {\small 0.08} & {\small 0.05} & {\small
0.05} & {\small 0.06} & & {\small 0.08} & {\small 0.04} & {\small 0.04} &
{\small 0.05} & & {\small 0.09} & {\small 0.04} & {\small 0.04} & {\small
0.03} \\ \hline
\end{tabular}
\bigskip
\raggedright {\footnotesize Note: Entries are rejection probabilities of the
LR test (\ref{LR}) when data are generated from (\ref{sim-model}) with $
{\Greekmath 010D} _{0}\left( s\right) =\sin (s)/2$. The dependence structure is given
in (\ref{DGP corr})\ with ${\Greekmath 011A} =0$.\ The significance level\ is $5\%$\ and
the results are based on 1000 simulations.}
\end{table}
\begin{table}[tbp]
\centering
\caption{Rej. Prob. of the LR Test with Cross-sectionally Correlated Data}
\label{tbl r3}
\bigskip
\begin{tabular}{cccccccccccccccc}
\hline\hline
& & \multicolumn{4}{c}{${\small s=0.0}$} & & \multicolumn{4}{c}{${\small
s=0.5}$} & & \multicolumn{4}{c}{${\small s=1.0}$} \\
\cline{3-6}\cline{8-11}\cline{13-16}
${\small n}$ & \multicolumn{1}{l}{${\small {\Greekmath 010E} =}$} & {\small 1} &
{\small 2} & {\small 3} & {\small 4} & & {\small 1} & {\small 2} & {\small 3
} & {\small 4} & & {\small 1} & {\small 2} & {\small 3} & {\small 4} \\
\hline
\multicolumn{1}{r}{\small 100} & & {\small 0.18} & {\small 0.09} & {\small
0.07} & {\small 0.08} & & {\small 0.21} & {\small 0.11} & {\small 0.10} &
{\small 0.06} & & {\small 0.28} & {\small 0.20} & {\small 0.17} & {\small
0.13} \\
\multicolumn{1}{r}{\small 200} & & {\small 0.12} & {\small 0.06} & {\small
0.06} & {\small 0.06} & & {\small 0.13} & {\small 0.08} & {\small 0.06} &
{\small 0.05} & & {\small 0.20} & {\small 0.37} & {\small 0.09} & {\small
0.06} \\
\multicolumn{1}{r}{\small 500} & & {\small 0.08} & {\small 0.04} & {\small
0.06} & {\small 0.06} & & {\small 0.06} & {\small 0.06} & {\small 0.04} &
{\small 0.05} & & {\small 0.13} & {\small 0.08} & {\small 0.05} & {\small
0.02} \\ \hline
\end{tabular}
\bigskip
\raggedright {\footnotesize Note: Entries are rejection probabilities of the
LR test (\ref{LR}) when data are generated from (\ref{sim-model}) with $
{\Greekmath 010D} _{0}\left( s\right) =\sin (s)/2$. The dependence structure is given
in (\ref{DGP corr})\ with ${\Greekmath 011A} =1$\ and $m=10$.\ The significance level\ is
$5\%$\ and the results are based on $1000$ simulations.}
\end{table}
First, Tables \ref{tbl r1} and \ref{tbl r3} report the small sample
rejection probabilities of the LR test in (\ref{LR}) for $H_{0}:{\Greekmath 010D}
_{0}(s)=\sin (s)/2$ against $H_{1}:{\Greekmath 010D} _{0}(s)\neq \sin (s)/2$ at the 5\%
nominal level at three different locations $s=0$, $0.5$, and $1$. In
particular, Table \ref{tbl r1} examines the case with no cross-sectional
dependence (${\Greekmath 011A} =0$), while Table \ref{tbl r3} examines the case with
cross-sectional dependence whose dependence decays slowly with ${\Greekmath 011A} =1$ and
$m=10$. For the bandwidth parameter, we normalize $s_{i}$ and $q_{i}$ to
have zero mean and unit standard deviation, and choose $b_{n}=0.5n^{-1/2}$
in the main regression. This choice is for undersmoothing so that $
n^{1-2{\Greekmath 010F} }b_{n}^{2}=n^{-2{\Greekmath 010F} }\rightarrow 0$. To estimate $
D\left( {\Greekmath 010D} _{0}\left( s\right) ,s\right) $ and $V\left( {\Greekmath 010D}
_{0}\left( s\right) ,s\right) $, we use the rule-of-thumb bandwidths from
the standard kernel regression satisfying $b_{n}^{\prime }=O(n^{-1/5})$ and $
b_{n}^{\prime \prime }=O(n^{-1/6})$. All the results are based on $1000$
simulations. In general, the test for ${\Greekmath 010D} _{0}$ performs better as (i)
the sample size gets larger; (ii) the coefficient change gets more
significant; (iii) the cross-sectional dependence gets weaker; and (iv) the
target gets closer to the mid-support of $s$. When ${\Greekmath 010E} _{0}$ and $n$ are
large, the LR test is conservative, which is also found in the classical
threshold regression (e.g., \cite{Hansen00a}).
\begin{table}[tbp]
\centering
\caption{Coverage Prob. of the Plug-in Confidence Interval}\label{tbl bd
plugin}
\bigskip
\begin{tabular}{lccccccccccccccc}
\hline\hline
& & \multicolumn{4}{c}{${\Greekmath 010C} _{20}$} & & \multicolumn{4}{c}{${\Greekmath 010C}
_{20}+{\Greekmath 010E} _{20}$} & & \multicolumn{4}{c}{${\Greekmath 010E} _{20}$} \\
\cline{3-6}\cline{8-11}\cline{13-16}
\multicolumn{1}{c}{${\small n}$} & \multicolumn{1}{l}{${\small {\Greekmath 010E} =}$} &
{\small 1} & {\small 2} & {\small 3} & {\small 4} & & {\small 1} & {\small 2
} & {\small 3} & {\small 4} & & {\small 1} & {\small 2} & {\small 3} &
{\small 4} \\ \hline
\multicolumn{1}{r}{\small 100} & & {\small 0.85} & {\small 0.87} & {\small
0.90} & {\small 0.89} & & {\small 0.82} & {\small 0.89} & {\small 0.88} &
{\small 0.88} & & {\small 0.83} & {\small 0.88} & {\small 0.89} & {\small
0.90} \\
\multicolumn{1}{r}{\small 200} & & {\small 0.87} & {\small 0.91} & {\small
0.91} & {\small 0.91} & & {\small 0.87} & \multicolumn{1}{l}{\small 0.90} &
{\small 0.93} & {\small 0.93} & & {\small 0.86} & {\small 0.90} & {\small
0.92} & {\small 0.94} \\
\multicolumn{1}{r}{\small 500} & & {\small 0.89} & {\small 0.92} & {\small
0.95} & {\small 0.94} & & {\small 0.87} & \multicolumn{1}{l}{\small 0.93} &
{\small 0.93} & {\small 0.94} & & {\small 0.85} & {\small 0.92} & {\small
0.95} & {\small 0.92} \\ \hline
\end{tabular}
\bigskip
\raggedright {\footnotesize Note: Entries are coverage probabilities of 95\%
confidence intervals for ${\Greekmath 010C} _{20}$, ${\Greekmath 010C} _{20}+{\Greekmath 010E} _{20},$\ and $
{\Greekmath 010E} _{20}$. Data are generated from (\ref{sim-model})\ with ${\Greekmath 010D}
_{0}\left( s\right) =\sin (s)/2$, where the dependence structure is given in
(\ref{DGP corr})\ with ${\Greekmath 011A} =0.5$\ and $m=3$.\ The results are based on
1000 simulations.}
\end{table}
\begin{table}[tbp]
\centering
\caption{Coverage Prob. of the Plug-in Confidence Interval (w/ LRV adj.)}
\label{tbl bd lrv}
\bigskip
\begin{tabular}{lccccccccccccccc}
\hline\hline
& & \multicolumn{4}{c}{${\Greekmath 010C} _{20}$} & & \multicolumn{4}{c}{${\Greekmath 010C}
_{20}+{\Greekmath 010E} _{20}$} & & \multicolumn{4}{c}{${\Greekmath 010E} _{20}$} \\
\cline{3-6}\cline{8-11}\cline{13-16}
\multicolumn{1}{c}{${\small n}$} & \multicolumn{1}{l}{${\small {\Greekmath 010E} =}$} &
{\small 1} & {\small 2} & {\small 3} & {\small 4} & & {\small 1} & {\small 2
} & {\small 3} & {\small 4} & & {\small 1} & {\small 2} & {\small 3} &
{\small 4} \\ \hline
\multicolumn{1}{r}{\small 100} & & {\small 0.94} & {\small 0.94} & {\small
0.94} & {\small 0.94} & & {\small 0.91} & {\small 0.95} & {\small 0.95} &
{\small 0.95} & & {\small 0.92} & {\small 0.94} & {\small 0.96} & {\small
0.96} \\
\multicolumn{1}{r}{\small 200} & & {\small 0.94} & {\small 0.95} & {\small
0.96} & {\small 0.96} & & {\small 0.93} & \multicolumn{1}{l}{\small 0.94} &
{\small 0.96} & {\small 0.96} & & {\small 0.92} & {\small 0.96} & {\small
0.97} & {\small 0.97} \\
\multicolumn{1}{r}{\small 500} & & {\small 0.94} & {\small 0.95} & {\small
0.98} & {\small 0.97} & & {\small 0.92} & \multicolumn{1}{l}{\small 0.96} &
{\small 0.97} & {\small 0.96} & & {\small 0.91} & {\small 0.97} & {\small
0.97} & {\small 0.96} \\ \hline
\end{tabular}
\bigskip
\raggedright {\footnotesize Note: Entries are coverage probabilities of 95\%
confidence intervals for ${\Greekmath 010C} _{20}$, ${\Greekmath 010C} _{20}$$+{\Greekmath 010E} _{20},$\ and $
{\Greekmath 010E} _{20} $\ with a small sample adjustment of the LRV estimator. Data
are generated from (\ref{sim-model})\ with ${\Greekmath 010D} _{0}\left( s\right)
\left. =\right. \sin (s)/2$, where the dependence structure is given in (\ref
{DGP corr})\ with ${\Greekmath 011A} =0.5$\ and $m=3$.\ The results are based on 1000
simulations.}
\end{table}
Second, Table \ref{tbl bd plugin} shows the finite sample coverage
properties of the 95\% confidence intervals for the parametric components $
{\Greekmath 010C} _{20}$, ${\Greekmath 010E} _{20}^{\ast }={\Greekmath 010C} _{20}+{\Greekmath 010E} _{20}$, and ${\Greekmath 010E}
_{20}$. The results are based on the same simulation design as above with $
{\Greekmath 011A} =0.5$ and $m=3$. Regarding the tuning parameters, we use the same
bandwidth choice $b_{n}=0.5n^{-1/2}$ as before and set the truncation
parameter ${\Greekmath 0119} _{n}=\left( nb_{n}\right) ^{-1/2}$. Unreported results
suggest that choice of the constant in the bandwidth matters particularly
with small samples like $n=100$, but such effect quickly decays as the
sample size gets larger. For the estimator of the LRV, we use the spatial
lag order of $5$ following \cite{Conley07}. Results with other lag choices
are similar and hence omitted. The result suggests that the asymptotic
normality is better approximated with larger samples and larger change
sizes. Table \ref{tbl bd lrv} shows the same results with\ a small sample
adjustment of the LRV estimator for $\Omega ^{\ast }$ by dividing it by the
sample truncation fraction $\sum_{i\in \Lambda _{n}}(\mathbf{1}[q_{i}>
\widehat{{\Greekmath 010D} }(s_{i})+{\Greekmath 0119} _{n}]+\mathbf{1}[q_{i}<\widehat{{\Greekmath 010D} }
(s_{i})-{\Greekmath 0119} _{n}])\mathbf{1}[s_{i}\in \mathcal{S}_{0}]/\sum_{i\in \Lambda
_{n}}\mathbf{1}[s_{i}\in \mathcal{S}_{0}]$. This ratio enlarges the LRV
estimator and hence the coverage probabilities, especially when the change
size is small. It only affects the finite sample performance as it
approaches one in probability as $n\rightarrow \infty $.
\section{Applications\label{Section empirics}}
\subsection{Tipping point and social segregation\label{Section tipping}}
The first application is about the tipping point problem in social
segregation, which stimulates a vast literature in labor, public, and
political economics. \cite{Schelling71} initially proposes the tipping point
model to study the fact that the white population decreases substantially
once the minority share exceeds a certain tipping point. \cite{Card08}
empirically estimate this model and find strong evidence for such a tipping
point phenomenon. In particular, they specify the threshold regression model
as
\begin{equation*}
y_{i}={\Greekmath 010C} _{10}+{\Greekmath 010E} _{10}\mathbf{1}\left[ q_{i}\leq {\Greekmath 010D} _{0}\right]
+x_{2i}^{\top }{\Greekmath 010C} _{20}+u_{i}\relax\protect\ifmmode\expandafter\text@\else\expandafter\mbox\fi{,}
\end{equation*}
where for tract $i$ in a certain city, $q_{i}$ is the minority share in
percentage at the beginning of a certain decade, $y_{i}$ is the normalized
white population change in percentage within this decade, and $x_{2i}$ is a
vector of control variables. They apply the least squares method to estimate
the tipping point ${\Greekmath 010D} _{0}$. For most cities and for the periods
1970-80, 1980-90, and 1990-2000, they find that white population flows
exhibit the tipping-like behavior, with the estimated tipping points ranging
approximately from 5\% to 20\% across cities.
In Section VII of \cite{Card08}, they also find that the location of the
tipping point substantially depends on white residents' attitudes toward the
minority. Specifically, they first construct a city-level index that
measures white attitudes and regress the estimated tipping point from each
city on this index. The regression coefficient is significantly different
from zero, suggesting that the tipping point is heterogeneous across cities.
We go one step further by considering a more flexible model in the tract
level given as
\begin{equation*}
y_{i}={\Greekmath 010C} _{10}+{\Greekmath 010E} _{10}\mathbf{1}\left[ q_{i}\leq {\Greekmath 010D} _{0}(s_{i})
\right] +x_{2i}^{\top }{\Greekmath 010C} _{20}+u_{i}\relax\protect\ifmmode\expandafter\text@\else\expandafter\mbox\fi{,}
\end{equation*}
where ${\Greekmath 010D} _{0}(\cdot )$ denotes an unknown tipping point function and $
s_{i}$ denotes the attitude index. The nonparametric function ${\Greekmath 010D}
_{0}(\cdot )$ here allows for heterogeneous tipping points across tracts
depending on the level of the attitude index $s_{i}$ in tract $i$.
Unfortunately, the attitude index by \cite{Card08} is only available at the
aggregated city-level, and hence we cannot use it to analyze the census
tract-level observations. For this reason, we instead use the tract-level
unemployment rate as $s_{i}$ to illustrate the nonparametric threshold
function, which is readily available in the original dataset. Such a
compromise is far from being perfect but can be partially justified since
race discrimination has been widely documented to be correlated with
employment (e.g., \cite{DarityMason1998}).
\begin{figure}[tbp]
\begin{center}
\caption{Estimate of the tipping point as a function of the unemployment
rate}\label{fig tipping}
\vspace{+4ex}
{\footnotesize Panel A: Estimated tipping point function in Chicago 1980-90}
\includegraphics[width=6.06in,height=2.03in]{fig_Chicago8090_20210119_choose_c05_55.jpg}
{\footnotesize Panel B: Estimated tipping point function in Los Angeles 1980-90}
\includegraphics[width=6.06in,height=2.03in]{fig_LA8090_20210119_choose_c05_55.jpg}
{\footnotesize Panel C: Estimated tipping point function in New York City 1980-90}
\includegraphics[width=6.06in,height=2.03in]{fig_NewYork8090_20210119_choose_c05_55.jpg}
\end{center}
\begin{small}
Note: The figure depicts the point estimate of the tipping points as a
function of the unemployment rate, using the data in Chicago, Los Angeles, and
New York City in 1980-1990. The vertical axis is the estimated tipping point
in percentage, and the horizontal axis is the tract-level unemployment
normalized into quantile level. Data are available from \cite{Card08}.
\end{small}
\end{figure}
We use the data provided by \cite{Card08} and estimate the tipping point
function ${\Greekmath 010D} _{0}(\cdot )$ over census tracts by the method introduced
in Section \ref{Section estimation}. As in their work, we drop the tracts
where the minority shares are above 60 percentage points and use five
control variables as $x_{2i}$, including the logarithm of mean family
income, the fractions of single-unit, vacant, and renter-occupied housing
units, and the fraction of workers who use public transport to travel to
work. The bandwidth is set as $b_{n}=cn^{-1/2}$ for some $c>0$, so that it
satisfies the technical conditions in the previous sections, where the
constant $c$ is chosen by the leave-one-out cross validation. In particular,
we first construct the leave-one-out estimate, $\widehat{{\Greekmath 010D} }_{-i}\left(
s_{i}\right) $, of ${\Greekmath 010D} _{0}\left( s_{i}\right) $ as in (\ref{g-hat0})
without using the $i$th observation. Then, leaving the $i$th observation
out, we construct $\widehat{{\Greekmath 010C} }_{-i}$ and $\widehat{{\Greekmath 010E} }_{-i}$ as in
(\ref{para-b}) and (\ref{para-d}) with ${\Greekmath 0119} _{n}=\left( nb_{n}\right)
^{-1/2} $ using the bandwidth $b_{n}$ chosen in the previous step. We choose
the bandwidth that minimizes $\sum_{i\in \Lambda _{n}}(y_{i}-\widehat{{\Greekmath 010C} }
_{1,-i}-\widehat{{\Greekmath 010E} }_{1,-i}\mathbf{1}\left[ q_{i}\leq \widehat{{\Greekmath 010D} }
_{-i}(s_{i})\right] -x_{2i}^{\top }\widehat{{\Greekmath 010C} }_{2,-i})^{2}\mathbf{1}
\left[ s_{i}\in \mathcal{S}_{0}\right] $, where $\mathcal{S}_{0}$ again
includes the middle 70\% quantiles of $s_{i}$.
Figure \ref{fig tipping} depicts the estimated tipping points and the 95\%
pointwise confidence intervals by inverting the likelihood ratio test
statistic (\ref{LR}) in the years 1980-90 in Chicago, Los Angeles, and New
York City, whose sample sizes are relatively large. For each city, the
constant $c$ of the bandwidth $b_{n}=cn^{-1/2}$ chosen by the aforementioned
cross validation is $3.20$, $4.87$, and $3.42$, respectively. We make the
following comments. First, the estimates of the tipping points vary
substantially in the unemployment rate within all three cities. Therefore,
the standard constant tipping point model is insufficient to characterize
the segregation fully. Second, the tipping points as functions of the
unemployment rate do not exhibit the same pattern across cities, reinforcing
the heterogeneous tipping points in the city-level as found in \cite{Card08}
. Finally, the estimated tipping point $\widehat{{\Greekmath 010D} }\left( s\right) $
as a function of $s$ can be discontinuous, which does not contrast with
Assumption A-(vi), that is, the true function ${\Greekmath 010D} _{0}\left( \cdot
\right) $ is smooth. The discontinuity comes from the fact that $\widehat{
{\Greekmath 010D} }\left( s\right) $ is obtained by grid search and can only take
values among the discrete points $\{q_{1},...,q_{n}\}$ in finite samples.
\subsection{Metropolitan area determination\label{Section boundary}}
The second application is about determining the boundary of a metropolitan
area, which is a fundamental question in urban economics. Recently,
researchers propose to use nighttime light intensity obtained by satellite
imagery to define metropolitan areas. The intuition is straightforward:
metropolitan areas are bright at night while rural areas are dark.
\begin{figure}[tbp]
\begin{center}
\caption{Nighttime light intensity in Dallas, Texas, in 2010}\label{fig raw}
\vspace{+4ex}
\includegraphics[width=3.646in,height=2.508in]{rawdata_dallas_stablelight_2010.jpg}
\end{center}
\begin{small}
Note: The figure depicts the intensity of the stable nighttime light in
Dallas, TX 2010. Data are available from https://www.ncei.noaa.gov/.\
\end{small}
\end{figure}
Specifically, the National Oceanic and Atmospheric Administration (NOAA)
collects satellite imagery of nighttime lights at approximately 1-kilometer
resolution since 1992. NOAA further constructs several indices measuring the
annual light intensity. Following the literature (e.g., \cite{Dingel19}), we
choose the \textquotedblleft average visible, stable
lights\textquotedblright\ index that ranges from 0 (dark) to 63 (bright).
For illustration, we focus on Dallas, Texas and use the data from the years
1995, 2000, 2005, and 2010. In each year, the data are recorded as a 240$
\times $360 grid that covers the latitudes from 32$^{\circ }$N to 34$^{\circ
}$N and the longitudes from 98.5$^{\circ }$W to 95.5$^{\circ }$W. The total
sample size is 240$\times $360=86400 each year. These data are available at
NOAA's website and also provided on the authors' website. Figure \ref{fig
raw} depicts the intensity of the stable nighttime light of the Dallas area
in 2010 as an example.
Let $y_{i}$ be the level of nighttime light intensity and $\left(
q_{i},s_{i}\right) $ be the latitude and longitude of the $i$th pixel, which
is normalized into the equally-spaced grids on $[0,1]^{2}$. To define the
metropolitan area, existing literature in urban economics first chooses an
\textit{ad hoc} intensity threshold, say 95\% quantile of $y_{i}$, and
categorizes the $i$th pixel as a part of the metropolitan area if $y_{i}$ is
larger than the threshold. See \cite{Dingel19}, \cite{Vogel19}, and
references therein. In particular, on p.3 in \cite{Dingel19}, they note that
\textquotedblleft \lbrack ...] the choice of the light-intensity threshold,
which governs the definitions of the resulting metropolitan areas, is not
pinned down by economic theory or prior empirical
research.\textquotedblright\ Our new approach can provide a data-driven
guidance of choosing the intensity threshold from the econometric
perspective.
\begin{figure}[tbp]
\begin{center}
\caption{Kernel density estimate of nighttime light intensity, Dallas 2010}
\label{fig ksdensity}
\includegraphics[width=4.1632in,height=2.9334in]{ksdensity_dallas_stablelight_2010.jpg}
\end{center}
\begin{small}
Note: The figure depicts the kernel density estimate of the strength of the
stable nighttime light in Dallas, TX 2010. Data are available from
https://www.ncei.noaa.gov/.\
\end{small}
\end{figure}
\begin{figure}[tbp]
\begin{center}
\caption{Metropolitan area determination in Dallas (color online)}
\label{fig cityboundary}
\includegraphics[width=3.0in,height=2.3609in]{dallas_stablelight_1995.jpg}
\includegraphics[width=3.0in,height=2.3609in]{dallas_stablelight_2000.jpg}
\includegraphics[width=3.0in,height=2.3609in]{dallas_stablelight_2005.jpg}
\includegraphics[width=3.0in,height=2.3609in]{dallas_stablelight_2010.jpg}
\end{center}
\begin{small}
Note: The figure depicts the city boundary determined by either the new
method or by taking the 0.95 quantile of nighttime light strength as the
threshold, using the satellite imagery data for Dallas, TX in the years
1995, 2000, 2005, and 2010. Data are available from
https://www.ncei.noaa.gov/.\
\end{small}
\end{figure}
To this end, we first examine whether the light intensity data exhibits a
clear threshold pattern. We plot the kernel density estimates of $y_{i}$ in
the year 2010 in Figure \ref{fig ksdensity}. The bandwidth is the standard
rule-of-thumb one. The estimated density exhibits three peaks at around the
intensity levels 0, 8, and 63. They respectively correspond to the rural
area, small towns, and the central metropolitan area. It shows that the
threshold model is appropriate in characterizing such a mean-shift pattern.
Now we implement the rotation and estimation method described in Section \ref
{Section contour}. In particular, we pick the center point in the bright
middle area as the Dallas metropolitan center, which corresponds to the
pixel point in the 181st column from the left and the 100th row from the
bottom. Then for each $a^{\circ }$ over the 500 equally-spaced grid on $
[0^{\circ },360^{\circ }]$, we rotate the data by $a^{\circ }$ degrees
counterclockwise and estimate the model (\ref{model3}) with $x_{i}=1$. The
bandwidth is chosen as $cn^{-1/2}$ with $c=1$. Other choices of $c$ lead to almost identical results,
given the large sample size. Figure \ref{fig cityboundary} presents the
estimated metropolitan areas using our nonparametric approach (red) and the
area determined by the \textit{ad hoc} threshold of the 95\% quantile of $
y_{i}$ (black) in the years 1995, 2000, 2005, and 2010. It clearly shows the
expansion of the Dallas metropolitan area over the 15 years of the sample
period.
Several interesting findings are summarized as follows. First, the estimated
boundary is highly nonlinear as a function of the angle. Therefore, any
parametric threshold model could lead to a substantially misleading result.
Second, our estimated area is larger than that determined by the \textit{ad
hoc} threshold, by 80.31\%, 81.56\%, 106.46\%, and 102.09\% in the years
1995, 2000, 2005, and 2010, respectively. In particular, our nonparametric
estimates tend to include some suburban areas that exhibit strong light
intensity and that are geographically close to the city center. For example,
the very left stretch-out area in the estimated boundary corresponds to Fort
Worth, which is 30 miles from downtown Dallas. Residents can easily commute
by train or driving on the interstate highway 30. It is then reasonable to
include Fort Worth as a part of the metropolitan Dallas area for economic
analysis. Third, given the large sample size, the 95\% confidence intervals
of the boundary are too narrow to be distinguished from the estimates and
therefore omitted from the figure. Such narrow intervals apparently exclude
the boundary determined by the \textit{ad hoc }method. Finally, the
estimated value of ${\Greekmath 010C} _{0}+{\Greekmath 010E} _{0}$ is approximately $53$ in these
sample periods, which corresponds to the 89\% quantile of $y_{i}$ in the
sample. This suggests that a more proper choice of the level of light
intensity threshold is the 89\% quantile of $y_{i}$, instead of the 95\%
quantile, if one needs to choose the light-intensity threshold to determine
the Dallas metropolitan area.
\section{Concluding Remarks\label{Section conclusion}}
This paper proposes a novel approach to conduct sample splitting. In
particular, we develop a nonparametric threshold regression model where two
variables can jointly determine the unknown threshold boundary. Our approach
can be easily generalized so that the sample splitting depends on more
numbers of variables, though such an extension is subject to the curse of
dimensionality, as usually observed in the kernel regression literature. The
main interest is in identifying the threshold function that determines how
to split the sample. Thus our model should be distinguished from the
smoothed threshold regression model or the random coefficient regression
model. It instead could be seen as an unsupervised learning tool for
clustering.
This new approach is empirically relevant in broad areas studying sample
splitting (e.g., segregation and group-formation) and heterogeneous effects
over different subsamples. We illustrate some of them with the tipping point
problem in social segregation and metropolitan area determination using
satellite imagery datasets. Though we omit in this paper, we also estimate
the economic border between Brooklyn and Queens boroughs in New York City
using housing prices.\footnote{
The result is available upon request.} The estimated border is substantially
different from the existing administrative border, which was determined in
1931 and cannot reflect the dramatic city development. Interestingly, the
estimated border coincides with the Jackson Robinson Parkway and the Long
Island Railroad. This finding provides new evidence that local
transportation corridors could increase community segregation (cf.\ \cite
{Ananat11} and \cite{Heilmann18}).
We list some related works, which could motivate potential theoretical
extensions. First, while we focus on the local constant estimation in this
paper, one could consider the local linear estimation using the threshold
indicator $\mathbf{1}\left[ q_{i}\leq {\Greekmath 010D} _{1}+{\Greekmath 010D} _{2}(s_{i}-s)\right]
$ in (\ref{sse}). Although grid search is very difficult in determining the
two threshold parameters (${\Greekmath 010D} _{1}$ and ${\Greekmath 010D} _{2}$), we could use the
MCMC algorithm developed by \cite{Yu19} and the mixed integer optimization
(MIO) algorithms developed by \cite{LLSS18}. Besides the computational
challenge, however, the asymptotic derivation is more involved since we need
to consider higher-order expansions of the objective function. Second, while
our nonparametric setup is on the threshold function ${\Greekmath 010D} _{0}(\cdot )$,
some recent literature studies the nonparametric regression model with a
parametric threshold, such as $y_{i}=m_{1}(x_{i})+m_{2}(x_{i})\mathbf{1}
[q_{i}\leq {\Greekmath 010D} _{0}]+u_{i}$, where $m_{1}\left( \cdot \right) $ and $
m_{2}\left( \cdot \right) $ are different nonparametric functions. See, for
example, \cite{Henderson17}, \cite{Chiou18}, \cite{Yu18}, \cite
{YuLiaoPhillips19}, and \cite{DelgadoHidalgo2000}.
{\footnotesize \newpage }
\setcounter{equation}{0}\scalefont{0.96}\baselineskip=14pt