EconBase
← Back to paper

Inference on many jumps in nonparametric panel regression models

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.

67,907 characters

Inference on many jumps in nonparametric panel regression models




\def\spacingset#1{\renewcommand{1.0}
{#1}\small\normalsize} \spacingset{1}



\if11
{
  \title{\bf Inference on many jumps in nonparametric panel regression models}
  \author{Likai Chen\thanks{Chen's research is partially supported by NSF 23-512 and PD 18-1269. Su
thanks the National Natural Science Foundation of China (NSFC) for financial
support under the grant number 72133002. Wang's research is partially
supported by the ESRC (Grant Reference: ES/T01573X/1).}\hspace{.2cm}\\
    Washington University in St. Louis\\
    and \\
    Georg Keilbar \\
    Humboldt-Universit\"at zu Berlin \\
    and \\
    Liangjun Su \\
    Tsinghua University \\
    and \\
    Weining Wang \\
    University of Bristol}
    \date{}
  \maketitle
} \fi

\if01
{
  \bigskip
  \bigskip
  \bigskip
  \begin{center}
    {\LARGE\bf Inference on many jumps in nonparametric panel regression models}
\end{center}
  \medskip
} \fi


\begin{abstract}
We investigate the significance of change-points within fully nonparametric regression contexts, with a particular focus on panel data where data generation processes vary across units, and error terms may display complex dependency structures. In our setting the threshold effect depends on one specific covariate, and we permit the true nonparametric regression to vary based on additional (latent) variables. We propose two uniform testing procedures: one to assess the existence of change-points and another to evaluate the uniformity of such effects across units. Our approach involves deriving a straightforward analytical expression to approximate the variance-covariance structure of change-point effects under general dependency conditions. Notably, when Gaussian approximations are made to these test statistics, the intricate dependency structures within the data can be safely disregarded owing to the localized nature of the statistics. This finding bears significant implications for obtaining critical values. Through extensive simulations, we demonstrate that our tests exhibit excellent control over size and reasonable power performance in finite samples, irrespective of strong cross-sectional and weak serial dependency within the data. Furthermore, applying our tests to two datasets reveals the existence of significant nonsmooth effects in both cases.
\end{abstract}

\noindent
{\it Keywords:} high-dimensional time series; temporal and cross-sectional dependence; changepoint analysis; Gaussian
approximation;  simultaneous tests; threshold regression.
\vfill

\newpage
\spacingset{1.9}


\section{\protect\normalsize Introduction}

  Recently, there has been a notable surge in the study and application of changepoint analysis in high-dimensional time series data (see, e.g., \cite{cho2015multiple}, \cite{cao2018multi}, \cite{chen2022inference}, and \cite{li2024}, along with the references cited therein). Detecting changepoints in regression models is important and has wide applications in the field of biology \citep{pastor1998use}, agricultural sciences \citep{freeman1998credit}, econometrics \citep{li2024}, etc. The purpose of these studies is to estimate not only where the threshold effect exists but also the magnitude of the effects.

{\normalsize
\label{lin:test1}
There is extensive literature on parametric high-dimensional changepoint regression analysis. \cite{li2016panel}, \cite{manner2019testing} \cite{wang2021statistically}, \cite{rinaldo2021localizing}, and \cite{cho2024detection} investigate changepoint analysis of coefficients in high-dimensional linear regression models using penalization methods. \cite{lee2016lasso} examine a high-dimensional regression model with a changepoint due to a covariate threshold in a cross-sectional study. These studies majorly focus on parametric models with high-dimensional covariates. Similar to \cite{lee2016lasso}, this paper also considers changepoint regression due to a covariate threshold. However, our focus is on nonparametric regression with a high-dimensional changepoint effect driven by heterogeneity in the panel data structure.
Specifically, we {investigate the heterogeneous
nonparametric panel data with jumps in the conditional mean functions: }
\begin{align}
\label{modelstructural}
Y_{jt}=h_{j}(X_{jt},U_{jt})+\tau _{j}(X_{jt},U_{jt})\mathbf{1}_{\{X_{jt}\geq
c_{0j}\}}+\sigma _{j}(X_{jt},U_{jt})\epsilon _{jt},
\end{align}
where {$X_{jt}$ is a one-dimensional threshold variable,}  $Y_{jt}\in \mathbb{R}^1$ and {$U_{jt}\in \mathbb{R}^d$} are the outcome and covariate variables that can be observed or unobserved for individual $j\in \lbrack N]:=\{1,\ldots ,N\}$ at time point $t\in \lbrack T],$} $h_{j}(\cdot,\cdot )$, $\tau_j(\cdot,\cdot)$ and $\sigma _{j}(\cdot,\cdot)$ are continuous functions, $c_{0j}$ is a threshold value that
is either known or unknown and $\epsilon_{jt}$ is the error term. Our main focus is to test
whether changepoint effects related to $\tau _{j}(X_{jt},U_{jt})$ exist and to determine if the jump sizes are identical uniformly across individuals $j$ under general spatial temporal dependence assumptions.


\label{page:answerae1}As a motivating example for our model, we explore an application involving stock lagged returns and volatility using threshold regression. Threshold regression, dating back to \cite{chan1985multiple} and more recently studied in works such as \cite{duffy2023stationarity}, is mainly formulated within a linear parametric setting. Our approach naturally extends this framework to the non-parametric setting. Specifically, in our first application, revising the news impact curve of \cite{engle1993measuring},  $Y_{jt}$ is the volatility of firm $j$'s stock at time $t$ and $X_{jt}$ is the lagged stock return. A trader would like to discover whether there are asymmetric effects of $X_{jt}$ on $Y_{jt}$,  as well as whether and where a discontinuous relationship occurs, for trading and risk management purposes. If there are multiple stocks jump, are they all identical uniformly? In Section \ref{sec:application} in our empirical section, we identify significant threshold effects, demonstrating that our methodology is crucial for predicting stock volatility and effectively managing risk.
In another application, we study the effect of party incumbency as
determined by the 50\% threshold in the two-party vote share in the US House
Elections. Our findings confirm the presence of a significant positive
incumbency effect as studied by \citet{lee2008randomized}, and we detect
significant heterogeneity in the threshold effect across the US states.

The methodological contributions are twofold. First, we develop algorithms to test for overall threshold effects, examine the potential heterogeneity effects and derive the asymptotic properties of the resulting test statistics. Due to the use of local smoothing,
neither the weak dependence along the time dimension nor the weak or strong
dependence along the individual/group dimension enters the asymptotic
distribution of the final statistic. Interestingly, we find strong
cross-sectional dependence is allowed in our error terms. As a practical
recommendation, we find that under fairly weak assumptions, if there are
clustering dependency patterns in $e_{jt}$, the correlation in the estimated statistics can be safely ignored, and only adjustments for heteroskedasticity are needed.
Second, the standard central limit theorem is challenging to apply to high-dimensional time series when $N$ significantly exceeds the effective sample size. Previously, \citet{chernozhukov2015comparison, chernozhukov_central_2017}, \cite{zhang2017gaussian}, \cite{chernozhukov2019inference}, and others addressed this issue by introducing high-dimensional Gaussian approximation (GA) methods. However, the existing literature cannot address our problem. The key difficulty is to incorporate possibly non-identically distributed time series under general spatial and temporal dependency structures and more importantly latent variables. In this paper, we propose a new high dimensional GA for weighted statistics that address all these issues. These GA results are essential for determining the critical values of our supremum test.

\label{contribution}Our paper is related to three strands of literature in statistics. First, we add to the extensive body of work on threshold regression for high dimensional data. For example, \cite{miao2020panel_a} and
\cite{miao2020panel} consider the estimation of a panel threshold regression
model with interactive fixed effects (IFEs) and latent group structures, respectively; \cite
{massacci2022high} study inference for high-dimensional threshold
regressions with common stochastic trends. The aforementioned literature does not account for potential heterogeneous threshold effects. Ignoring these effects may bias estimators or lead to efficiency loss when estimating homogeneous threshold effects in finite samples. Second, our
paper contributes to the literature on testing homogeneity in a panel
data setting. For example, \citet{phillips2003dynamic},
\citet{pesaran2008testing} and \citet{su2013testing} propose tests for slope
coefficient heterogeneity in a linear regression model with or without
cross-sectional dependence. In the presence of slope heterogeneity, \cite
{su2016identifying} and \cite{su2018identifying} model the slope
heterogeneity via latent group structures. \cite{barassi2023threshold} study
a heterogeneous panel threshold regression model with IFEs. Above literature all focus on parametric setting. Our work complements these research streams by adopting nonparametric approaches, which address size and power distortions caused by model misspecifications.
Third, our nonparametric approach is related with literature in trends estimation. Inference on smooth trend functions in nonparametric models for time series data is a well-established field; see \cite{zhang2012inference}, \cite{karmakar2022simultaneous}, \cite{gao2024time}, \cite{chen2018testing}, among others.
Distinctly, while the above literature examines nonlinear dynamic models, they does not address discontinuity in trends, particularly heterogeneous jumps.
The nonparametric literature on jumps in trends mostly focuses on independent and identically distributed (i.i.d.) settings, without accounting for heterogeneity or high dimensionality; see, for example, \cite{qiu1998discontinuous}, \cite{muller1999discontinuous}, and \cite{spokoiny1998estimation}.


\label{anwer2ae:latent}Another distinct feature is that we explicitly handle the case with missing covariates or models with misspecified variables. Theoretically, validating such an inference procedure requires working with a reduced-form model that has specific variance-covariance structures. This aspect of theoretical validity has not been explored in the time-varying trend literature. Establishing the validity of uniform inference is a major theoretical contribution that sets our work apart in the study of varying trends.



{\normalsize
}



{\normalsize
}

{\normalsize
}

{\normalsize
}


 The rest of the paper is organized as follows. Section \ref
{sec:model} provides the model setup and the practical steps of the testing
procedure. Section \ref{sec:theorem} presents the assumptions and delivers
the main theorems. In Section \ref{sec:simulation} we provide simulation results.\footnote{The code for our method is available in the R package \texttt{hdthreshold} (\url{https://cran.r-project.org/web/packages/hdthreshold/index.html}), with an illustration accessible at \url{https://rpubs.com/Lk1110/1259698}.} Real data applications are discussed in Section \ref{sec:application}, and the paper concludes in Section \ref{sec:conclusion}. The Supplementary Material includes additional theorems, proofs of all results, extra simulation outcomes, and tables for the empirical applications.

 \textit{Notation.} For a vector $v=(v_{1},\ldots,v_{d})\in
\mathbb{R}^{d}$ and $q>0$, we denote $|v|_{q}=(
\sum_{i=1}^{d}|v_{i}|^{q})^{1/q}$ and $|v|_{\infty }=\max_{1\leq i\leq
d}|v_{i}|$. For a matrix $A=(a_{i,j})_{1\leq i\leq m,1\leq j\leq n}$, we
define the max norm $|A|_{\text{max}}=\max_{i,j}|a_{i,j}|$. For $s>0$ and a
random vector $X$, we say $X\in \mathcal{L}^{s}$ if $\lVert X\rVert _{s}=[
\mathbb{E}(|X|^{s})]^{1/s}<\infty $. For two positive number sequences $
(a_{n})$ and $(b_{n})$, we say $a_{n}=O(b_{n})$ or $a_{n}\lesssim b_{n}$
(resp. $a_{n}\asymp b_{n}$) if there exists $C>0$ such that $a_{n}/b_{n}\leq
C$ (resp. $1/C\leq a_{n}/b_{n}\leq C$) for all large $n$, and say $
a_{n}=o(b_{n})$ if $a_{n}/b_{n}\rightarrow 0$ as $n\rightarrow \infty $. We
set $(X_{n})$ and $(Y_{n})$ to be two sequences of random variables. Write $
X_{n}=O_{\mathbb{P}}(Y_{n})$ if for any $\epsilon >0$, there exists $C>0$
such that $\mathbb{P}(X_{n}/Y_{n}\leq C)>1-\epsilon $ for all large $n$, and
say $X_{n}=o_{\mathbb{P}}(Y_{n})$ if $X_{n}/Y_{n}\rightarrow 0$ in
probability as $n\rightarrow \infty $.


\section{\protect\normalsize Model Setup and Test Procedure}

{\normalsize \label{sec:model}
}

{\normalsize In this section, we present the model, the hypotheses and
estimators, and the test procedure. Subsections \ref{subsec:model} and \ref{subsec:hypotheses} are concerned with the known threshold case. In Subsection \ref{unknown_c1} we focus on the case in which the threshold locations $c_{0j}$ are unknown. The algorithms for our testing procedures are listed in Subsection \ref{eq:algo}.
}

\subsection{\protect\normalsize Model Setup}\label{subsec:model}
{
\label{modelsetup} In this subsection, we formulate our model.
Recall the full model in (\ref{modelstructural})
where $U_{jt} \in \mathbb{R}^{d}$ is a latent random vector that can be arbitrarily
correlated with $X_{jt},$ $\epsilon _{jt}$ is a mean zero term that is
independent of $\left( X_{jt},U_{jt}\right) $.


Without loss of generality set $c_{0j} = 0$ for the known case.
The full model in \eqref{modelstructural} can be rewritten into the following reduced form model:
\begin{align}
&Y_{jt}=\tilde{h}_{j}(X_{jt})+\tilde{\tau}_{j}(X_{jt})\mathbf{1}_{\{X_{jt}\geq c_{0j}\}}+e_{jt},
\label{mainmodel}
\end{align}
where $\tilde{\tau}_{j}(X_{jt})=\mathbb{E}[\tau _{j}(X_{jt},U_{jt})|X_{jt}]$,  $\tilde{h}_{j}(X_{jt})=\mathbb{E[}h_{j}(X_{jt},U_{jt})|X_{jt}]$, $\varepsilon _{j}(X_{jt},U_{jt})=h_{j}(X_{jt},U_{jt})-\tilde{h}
_{j}(X_{jt})+[\tau _{j}(X_{jt},U_{jt})-\tilde{\tau}_{j}(X_{jt})]\mathbf{1}
_{\{X_{jt}\geq c_{0j}\}}$ and the reduced form error
\begin{align}
\label{eq:deffjvarepj}
e_{jt}=\varepsilon _{j}(X_{jt},U_{jt})+\sigma _{j}(X_{jt},U_{jt})\epsilon
_{jt}.
\end{align}
The noise term $e_{jt}$ has a complex structure, influenced by both the latent variable and by the heterogeneous volatility of the observations. This intricate structure, combined with the temporal and cross-sectional dependencies of the noise, makes the inference process particularly challenging.

A few remarks are in order. First, we have $\mathbb{E}(e_{jt}| X_{jt})=0$ and
the conditional expectation of the observed outcome given $X_{jt}$ is $
\mathbb{E}(Y_{jt}|X_{jt})=\tilde{h}_{j}(X_{jt})+\tilde{\tau}_{j}(X_{jt})
\mathbf{1}_{\{X_{jt}\geq c_{0j}\}}.$
Second, in Section \ref{addtionalcovariate} we show the possibilities of extending our model to add more covariates and incoporating fixed effects.
Third, ($X_{jt},U_{jt}$) can be lagged variables. For example, in our application on modeling stock volatility (Section \ref{application1}), $X_{jt} = Y_{j(t-1)}$, functions $h_{j}(\cdot ,\cdot
) $, {$\tau _{j}(\cdot ,\cdot )$} and $\sigma _{j}(\cdot ,\cdot )$ represent baseline effect, the
jump effect and the volatility function, respectively. The $U_{jt}$ are unobserved risk factors or variables we exclude in the nonparametric regression.
Modeling $U_{jt}$ is important in this circumstance as the volatility
of a single stock is most likely affected by other risk factors such as news
events or unobserved market factors.


This paper aims to explore methods for conducting simultaneous inference on threshold effect estimators in a large $N$ and large $T$
setup. In particular, we allow $N\gg T.$ Our focus is on the threshold effect, starting with the case where the threshold value $c_{0j}$
is known and then extending the analysis to the unknown case in Subsections \ref
{unknown_c1} and \ref{unknown_c2}.

{\normalsize Assuming that $X_{jt}$ is a continuous random variable, the
jump effects can be identified as a \textquotedblleft gap" in the
conditional expectations and
their derivatives at the thresholds. Let $\mu _{j}(x):=\mathbb{E}(Y_{jt}|X_{jt}=x)$ and $\mu
_{j}(c_{0j}-)=\lim_{x\uparrow c_{0j}}\mu _{j}(x)$ and $\mu
_{j}(c_{0j}+)=\lim_{x\downarrow c_{0j}}\mu _{j}(x).$ Consider $\partial
_{+}\mu _{j}(\cdot )$ (resp. $\partial _{-}\mu _{j}(\cdot )$) as the right
(resp. left) derivative of function $\mu _{j}(\cdot ).$  The threshold effect
for the individual $j$ is defined as
\begin{eqnarray}
\gamma _{j}&=&\mu
_{j}(c_{0j}+)-
\mu
_{j}(c_{0j}-),\label{gammaj1}\\
\gamma _{j}^{[1]}&=&\partial _{+}\mu _{j}(c_{0j} )-
\partial _{-}\mu _{j}(c_{0j}). \label{gammaj}
\end{eqnarray}
}
Due to continuity of function $h_j(\cdot,\cdot)$ and $\tau_j(\cdot,\cdot)$, we have
$\gamma _{j}=\tilde{\tau}_{j}(c_{0j}),\gamma _{j}^{[1]}=\tilde{\tau}'_{j}(c_{0j}).$ ($\tilde{\tau}'_{j}(.)$ denotes the derivative of $\tilde{\tau}_{j}(.)$.)



\subsection{\protect\normalsize Hypotheses and Estimators}\label{subsec:hypotheses}
 In this section, we introduce the hypotheses and estimators. We
would like to conduct simultaneous inferences on $\gamma _{j}$ and $\gamma_j^{[1]}$ for $
j\in \left[ N\right] $. The importance of studying individual/group-specific
threshold effects is discussed in \citet{zimmert2019nonparametric}.



{\normalsize First, we are interested in testing the existence of
the threshold effects ($\gamma_j$ or $\gamma_j^{[1]}$).
The null and alternative hypotheses are defined as follows:
\begin{eqnarray}
&&H_{0}^{\left( 1\right) }:\gamma_{j} =0\ \forall \text{ }j\in \left[ N\right]
,\quad \text{and }H_{a}^{\left( 1\right) }:\gamma _{j}\neq 0\ \text{for some
}j\in \left[ N\right] .  \label{null1}
\end{eqnarray}

 The motivation for considering such a test is crucial. In Figure \ref{figure:comparison}, we present a simple comparison with the pooled threshold test in the literature, for example in \cite{fong2017chngpt}. It is clear that the conventional pooled tests exhibit very low power due to signal cancellation, with performance only comparable in cases of strictly positive signals.

\begin{figure}[tbp]
\centering
\includegraphics[width=.49\textwidth]{normal.pdf}\hfill
\includegraphics[width=.49\textwidth]{exponential.pdf}
\caption{Empirical power comparison between our uniform testing procedure
(black) for the existence of threshold effects and the test based on the pooled linear threshold model (red) at $5\%$ significance level for different proportions of non-zero coefficients ($\gamma_j$) on the x-axis. The non-zero $\gamma_j$s are drawn from a standard normal distribution (left panel) and an exponential distribution (right panel). The true process is linear, $f_j(x)=0.2x$ (for details see DGP7 in Section \ref{sec:simulation_appendix} of the supplement.)}
\label{figure:comparison}
\end{figure}

}

 When $N$ is large, assuming homogeneous threshold effects
$\gamma_j$ across all $j$ is restrictive, and inferences based on this assumption can be misleading if the threshold effects are heterogeneous.  Therefore, it is important to test for heterogeneous threshold effects. In this case, the null and alternative hypotheses are:
\begin{eqnarray*}
&&H_{0}^{\left( 2\right) }:\gamma _{j}=\gamma _{0}\ \text{for some }\gamma _{0},
\text{ }\forall \text{ }j\in \left[ N\right],\quad
\\
\text{and}\quad&&
H_{a}^{\left( 2\right) }:\text{there is no constant }
\gamma_0 \text{ such that }\gamma
_{j}=\gamma_0, \text{ }\forall \ j\in \left[ N\right] .  \label{test}
\end{eqnarray*}





{We can similarly define $H_0^{(1)[1]}$ and $H_0^{(2)[1]}$ for the existence and homogeneity tests of the derivative, respectively. }Rejecting of $H_{0}^{\left( 2\right) }$ suggests the
presence of heterogeneous threshold effects. Moreover, providing simultaneous confidence intervals for the threshold
effects is also tempting. If we would like to study the asymptotics of $\hat{
\gamma}=(\hat{\gamma}_{1},\ldots,\hat{\gamma}_{N})^{\top }$, the standard
central limit theorem may fail due to the high dimensionality $N$. One of the key theoretical contributions
is to provide a general framework for making uniform inferences on
heterogeneous threshold effects.


{\normalsize
}

{\normalsize



To estimate $\gamma
_{j} $, we need to estimate both $\mu _{j}(c_{0j}+)$ and $\mu _{j}(c_{0j}-).$
To this aim, it is natural to adopt local linear estimators, {see for example  \cite{fan2018local}}. Denote the
estimators:
\begin{align}
\big(\hat{\mu}_{j}(c_{0j}+),{\partial _{+}\hat{\mu}_{j}}(c_{0j})\big)=&
\mathrm{argmin}_{\beta _{j,0},\beta _{j,1}}\sum_{t=1}^{T}(Y_{jt}-\beta
_{j,0}-\beta _{j,1}X_{jt})^{2}K_{jt}\mathbf{1}_{\{X_{jt}\geq c_{0j}\}}, \\
\big(\hat{\mu}_{j}(c_{0j}-),{\partial _{-}\hat{\mu}_{j}}(c_{0j})\big)=&
\mathrm{argmin}_{\beta _{j,0},\beta _{j,1}}\sum_{t=1}^{T}(Y_{jt}-\beta
_{j,0}-\beta _{j,1}X_{jt})^{2}K_{jt}\mathbf{1}_{\{X_{jt}<c_{0j}\}},
\end{align}
where }$K_{jt}:=K((X_{jt}-c_{0j})/b_{j}),${\normalsize \ $b_{j}$ is a
bandwidth parameter, and $K(\cdot )$ is a kernel function.
For $l=0,1,2$,
define the weights as
\begin{align}
w_{j,b}^{+}(u,v)& =\frac{K((v-u)/b_{j})[S_{j2,b}^{+}(u)-(v-u)S_{j1,b}^{+}(u)]
}{S_{j2}^{+}(u)S_{j0}^{+}(u)-S_{j1,b}^{+2}(u)}\mathbf{1}_{\{v\geq u\}},  \label{eq:wjbuv} \\
w_{j,b}^{+[1]}(u,v)& =\frac{K((v-u)/b_{j})[S_{j1,b}^{+}(u)-(v-u)S_{j0,b}^{+}(u)]
}{-S_{j2}^{+}(u)S_{j0}^{+}(u)+S_{j1,b}^{+2}(u)}\mathbf{1}_{\{v\geq u\}}
\label{eq:wjbuv[1]}.
\end{align}
where
\begin{align}
S_{jl,b}^{+}(u)& =\sum_{t=1}^{T}(X_{jt}-u)^{l}K((X_{jt}-u)/b_{j})\mathbf{1}
_{\{X_{jt}\geq u\}},
\end{align}
and $w_{j,b}^{-}(u,v)$ (resp. $S_{jl,b}^{-}(u)$, $w_{j,b}^{-[1]}(u,v)$) is the same as $
w_{j,b}^{+}(u,v)$ (resp. $S_{jl,b}^{+}(u)$, $w_{j,b}^{+[1]}(u,v)$) with $+$ and $\geq $ therein
replaced by $-$ and $<,$ respectively.
{\color{red}


}




We can thus evaluate the threshold
effect near the cut-off point and define our local linear estimator of the
average threshold effect for individual/group $j$ as
\begin{align}
\hat{\mu}_{j}(c_{0j}+)
=\sum_{t=1}^{T}w_{jt,b}^{+}Y_{jt}
\quad\textrm{and}\quad {\partial _{+}\hat{\mu}_{j}}(c_{0j})=\sum_{t=1}^{T}w_{jt,b}^{+[1]}Y_{jt},
\end{align}
where the weights take the following form $
w_{jt,b}^{+}=w_{j,b}^{+}(c_{0j},X_{jt})$ and $
w_{jt,b}^{-}=w_{j,b}^{-}(c_{0j},X_{jt})$. }
{
Thus we have,
\begin{eqnarray*}
&\hat{\gamma}_{j}=
\sum_{t=1}^{T}w_{jt,b}^{+}Y_{jt}-\sum_{t=1}^{T}w_{jt,b}^{-}Y_{jt} \quad \mbox{and} \quad\hat{\gamma}_{j}^{[1]}=
\sum_{t=1}^{T}w_{jt,b}^{+[1]}Y_{jt}-\sum_{t=1}^{T}w_{jt,b}^{-[1]}Y_{jt}.
\label{estimate}
\end{eqnarray*}
We shall focus on $\hat{\gamma}_{j}$ from now on, and the algorithms and theorems concerning the statistical properties of $\hat{\gamma}_{j}^{[1]}$ are developed in Section \ref{derivatives} in the Appendix.

}
Under $H_{0}^{\left( 1\right) }:$ $\gamma _{j}=0$ for all $j\in
\left[ N\right] $, we are exclusively testing for the overall significance
of threshold effects. Since the functions $\tilde{h}_{j}(\cdot)$ and $\tilde{\tau}_{j}(\cdot)$ are continuous, by \eqref{mainmodel}, the difference between our estimator and the
parameter can be approximately expressed in terms of a weighted average of the
error sequence,
\begin{equation}
\hat{\gamma}_{j}\approx \gamma
_{j}+\sum_{t=1}^{T}(w_{jt,b}^{+}e_{jt}-w_{jt,b}^{-}e_{jt}),
\label{eq:hatgammajgammaj}
\end{equation}
where $e_{jt}$ is the composite error in (\ref{mainmodel}
). Under $H_{0}^{\left( 1\right) }$, one would expect $\left\vert \hat{
\gamma}_{j}\right\vert $ to be small for all $j\in \left[ N\right] $. {This
motivates us to consider the test statistic
\begin{equation}
\mathcal{I}=\max_{1\leq j\leq N}(Tb_{j})^{1/2}|\hat{\gamma}_{j}/v_{j}|,
\label{I1}
\end{equation}
where
\begin{equation}
v_{j}^{2}=v_{j}^2(c_{0j})=(Tb_{j})
\sum_{t=1}^{T}(w_{jt,b}^{+}-w_{jt,b}^{-})^{2}\mathrm{Var}\big(e_{jt}|X_{jt}\big)
\label{eq:def_sigma}
\end{equation}
}is the conditional variance {of }$(Tb_{j})^{1/2}\hat{\gamma}_{j}$ for
standardization purposes. For now, we assume $v_{j}^{2}$ is given, in
general, it is unknown and has to be replaced by its estimate. We discuss a
consistent estimator of $v_{j}^{2}$ in Subsection \ref{eq:algo}. Under the
assumptions in {Section \ref{sec:theorem}}, it can be shown that all
asymptotic results based on $v_{j}^{2}$ remain valid when it is replaced by
the consistent estimate $\hat{v}_{j}^{2}$.

Similarly, under $H_{0}^{\left( 2\right) }$, one would expect
$\left\vert \hat{\gamma}_{j}-\overline{\hat{\gamma}}\right\vert $ to be
small for all $j\in \left[ N\right] ,$ where $\overline{\hat{\gamma}}=\frac{1
}{N}\sum_{j=1}^{N}\hat{\gamma}_{j}$. This motivates us to consider the test
statistic
\begin{equation}
\mathcal{Q}=\max_{1\leq j\leq N}(Tb_{j})^{1/2}|\hat{\gamma}_{j}-\overline{
\hat{\gamma}}|/\tilde{v}_j,  \label{Q1}
\end{equation}
where $\tilde{v}_{j}^{2}=(1-1/N)^{2}v_{j}^{2}+\sum_{i\neq
j}^{N}v_{i}^{2}/N^{2}$. For practical implementation, we introduce a feasible version
of $\mathcal{Q}$ in Algorithm \ref{alg2} in Section \ref{sec:algo} of the Supplementary Material.


\subsection{{\protect\normalsize The Unknown Threshold Case \label{unknown_c1}}}

{\normalsize It shall be noted that in case the breakpoints $c_{0j}$'s are
not known, we shall develop an algorithm to estimate them. We therefore
generalize our methodology to estimate unknown breakpoints that vary over }$
{\normalsize j}$. To this end, we search over a grid $[c_1,c_{2},\ldots ,c_K]$ for each individual $j\in \lbrack N]$.
Let
\begin{align}
\label{eq:wjtbci}
w_{jt,b}(c_{i})=w_{j,b}^{+}(c_{i},X_{jt})-w_{j,b}^{-}(c_{i},X_{jt}) \quad \mbox{and}\quad
\hat{\gamma}_{j}(c_{i})=\sum_{t=1}^{T}w_{jt,b}(c_{i})Y_{jt}.
\end{align}
Further, let $v_{j}^{2}(c_{i})=(Tb_{j})\sum_{t=1}^{T}w_{jt,b}^{2}(c_{i})\mathrm{Var}
\big(e_{jt}|X_{jt}\big)$ and define
\begin{equation}
\mathcal{I}^{C}=\max_{1\leq i\leq K}\max_{1\leq j\leq N}(Tb_{j})^{1/2}|\hat{
\gamma}_{j}(c_{i})/v_{j}(c_{i})|.
\end{equation}
We propose to estimate $c_{0j}$ by $\hat{c}_{j}:=\arg \max_{1\leq i\leq
K}(Tb_{j})^{1/2}|\hat{\gamma}_{j}(c_{i})/\hat{v}_{j}(c_{i})|\ $for $j\in
\left[ N\right] ,${\normalsize \ where }$\hat{v}_{j}(c_{i})${\normalsize \
is a consistent estimator of }$v_{j}(c_{i})$ which will be introduced in Algorithm
 \ref{algunknown} in Section \ref{sec:algo} of the Supplementary Material.

\subsection{{\protect\normalsize Test Procedure \label{eq:algo} }}

In this subsection, we present the steps of the test procedures. To
proceed, we begin by discussing the estimation of the error variance $\sigma
_{e,j}^{2}$ of $e_{jt}$ when $c_{0j}$'s are known.

\begin{enumerate}
\item Recall the definition of $\hat{\gamma}_j$ as in Equation (\ref{estimate}). Regress $Y_{jt}-\hat{\gamma}_{j}\mathbf{1}_{\{X_{jt}\geq c_{0j}\}}$ on $X_{jt}$ using the
local linear method with kernel $K\left( \cdot \right) $ and bandwidth $
b_{j} $.

\item Calculate the regression residuals $\hat{e}_{jt}$.

\item  Estimate $\sigma_{e,j}^{2}$ by
\begin{align}
\label{eq:sigmaej2}
\hat{\sigma }_{e,j}^{2}=\frac{\sum_{t=1}^{T}
\mathbf{1}_{\{|X_{jt}-c_{0j}|\leq b_{j}\}}\hat{e}_{jt}^{2}}{\sum_{t=1}^{T}
\mathbf{1}_{\{|X_{jt}-c_{0j}|\leq b_{j}\}}} .
\end{align}
\end{enumerate}
The estimator $\hat{\sigma}_{e,j}^2$ is used to form $\hat{v}_{j}^{2}$ in the Algorithms \ref
{alg1} and \ref{alg2} in the case of known $c_{0j}$, which are
applied to conduct simultaneous tests for $H_{0}^{\left( 1\right) }$ and $
H_{0}^{\left( 2\right) }$. These algorithms are supported by
the theoretical GA results in Section \ref{sec:theorem} and Section \ref{sec:additionalthm} in Supplementary Material.
The following Algorithm \ref{alg1} is for testing $H_0^{(1)}$.
\begin{algorithm}
	\caption{Testing for existence of jumps in conditional means (in slope): i.e., $\gamma_j =0$ {( $\gamma_j^{[1]} =0$}) for all $j$.}\label{alg1}
	\begin{algorithmic}[1]
		\For {$j=1,\ldots,N$}
			\State Estimate $\gamma_j$ for each group/individual $j$ using equation (\ref{estimate}).
			\State Estimate the variance $v_j^2$ of $(Tb_j^3)^{1/2}\hat{\gamma}_j$ by $\hat{v}_j^2=Tb_j\sum_{t=1}^{T}\left(w_{jt,b}^{+}-w_{jt,b}^-\right)^2\hat{\sigma}_{e,j}^2$ {( $\hat{v}_j^{2[1]}=Tb_j^3\sum_{t=1}^{T}\left(w_{jt,b}^{+[1]}-w_{jt,b}^{-[1]}\right)^2\hat{\sigma}_{e,j}^2$}) , where $\hat{\sigma}_{e,j}^2$ is a consistent estimator for the error variance using the observations close to the threshold $c_{0j}$ as in \eqref{eq:sigmaej2}.
		\EndFor
    \State Calculate the feasible test statistic $\hat{\mathcal{I}}:=\mbox{max}_{1\leq j\leq N}(Tb_j)^{1/2}|\hat{\gamma}_j|/\hat{v}_j$
    {( $\hat{\mathcal{I}}^{[1]}:=\mbox{max}_{1\leq j\leq N}(Tb_j^3)^{1/2}|\hat{\gamma}_j^{[1]}|/\hat{v}_j^{[1]}$})
    .
    \State Reject the null hypothesis if $\hat{\mathcal{I}}>q_{\alpha}$ ( $\hat{\mathcal{I}}^{[1]}>q_{\alpha}$), where $q_{\alpha}$ is the $(1-\alpha)$ quantile of the Gaussian random variable $\max_{1\leq j\leq N}|Z_j|$ with $Z_j$'s being independent standard normal variables.
	\end{algorithmic}
\end{algorithm}

Furthermore, Algorithm \ref{alg2} in Section \ref{sec:algo} in the Supplementary Material is used for testing the homogeneity hypothesis $H_0^{(2)}$. For the case where the threshold is unknown, the corresponding algorithm is detailed in Algorithm \ref{algunknown} in the same Section.
Note that $\mathcal{\hat{I}}$ in Algorithm \ref{alg1} (resp. $\mathcal{\hat{Q}}$ in Algorithm \ref{alg2}, $\mathcal{\hat{I}}^{C}$ in Algorithm {\ref{algunknown}})  is a feasible version of $\mathcal{I}$ (resp. $\mathcal{Q}$, $\mathcal{{I}}^{C}.$). In Section \ref{sec:theorem}, we will study the asymptotic properties of
$\mathcal{\hat{I}}$, $\mathcal{\hat{Q}}$ and $
\mathcal{\hat{I}}^{C}.$



Given the well-documented cross-sectional and serial dependence in high dimensional data, it is crucial to develop tests that account for both types of dependence.
{\normalsize It is surprising that, according to the above algorithm, $q_{\alpha}$ can be obtained by simulating only i.i.d. Gaussian random variables, without considering the complicated dependency structure of $e_{jt}$.
Recall from (\ref{mainmodel}), we see that
the reduced-form error $e_{jt}$ inherits nontrivial dependency structure in $
\epsilon _{jt}$. {However, this dependence does not affect the
calculation of critical values, even in the case of fairly strong cross-sectional
dependence (e.g., a factor structure). This is mainly due to the localized nature of
our kernel-type estimators $\hat{\gamma }_{j}$'s in our test statistic.
In fact, we can show that under some
mild conditions, as long as the joint density of $X_{it}$ and $X_{jt}$ is not
degenerate for any $\left( i,j\right) $ pair (see Assumption \ref
{asmp:bddxj1j2density} in Section \ref{sec:theorem}), the covariance terms
between $\hat{\gamma}_{i}$ and $\hat{\gamma}_{j}$ will be of lower order compared to the variance terms, $O(1/T)$.
For a more detailed discussion on the key intuition on the impact of the dependence structure in a simplified setting, we refer to Section \ref{intuition} of the Supplementary Material.} For a discussion on the issue of bandwidth selection, we refer to
Remark \ref{remark_bandwith} in Section \ref{homo} in the Supplementary Material.}


\section{{\protect\normalsize Main Theorems \label{sec:theorem}}}

In this section, we present our main results, with a particular focus on Theorem \ref{thm:ga}, which provide Gaussian approximation results that determine the critical values for testing the existence of threshold effects.
Theorem \ref{thm:unknown} extends the result to unknown threshold case and Theorem \ref{consistency} shows the consistency of the estimated threshold locations.
Theorem \ref{thm:homotesting} shows GA for testing homogeneity.
Theorems \ref{thm:hatv} and \ref{thm:hatvap} derive the consistency of variance to ensure the validity of our test statistics. Theorem \ref{thm:power} analyzes the power of the tests for both the threshold effects and homogeneity. Results regarding the derivative cases are presented in Subsection \ref{derivatives}.



\subsection{\protect\normalsize Assumptions}
\label{assum}
In this subsection, we state the assumptions for the asymptotic
analyses in Section \ref{assum} and Section \ref{subsec:additionalassumption} in the Supplementary Material.  First of all, we will pose a general spatial and
temporal dependence assumption on the data generating processes. As for the covariate variable $X_{jt}$ and the
confounding variable $U_{jt} \in \mathbb{R}^d$, we assume that they take the form
\begin{equation}
X_{jt}=H_{j}(e_{t},e_{t-1},\ldots )\quad \text{and}\quad U_{jt}=\tilde{H}
_{j}(X_{jt},\tilde{e}_{t}),  \label{eq:xjtujtdef}
\end{equation}
where
$\tilde{H}_{j}=(\tilde{H}_{j1},\ldots, \tilde{H}_{jd})^\top$,  $\{\tilde{H}_{ji}\}_{j\in \lbrack N], i\in \lbrack d]}$ and
 $\left\{ H_{j}\right\} _{j\in \lbrack N]}$ are functions
so that $X_{jt}$ and $U_{jt}$ are well defined. For $t\in \mathbb{Z}$,   $e_{t}\in\mathbb{R}$ denote i.i.d. innovations, and $\tilde{e}
_{t}\in\mathbb{R}^{\tilde d}$ are i.i.d. random vectors with some constant $\tilde d>0$, independent of $\left\{
e_{t}\right\} _{t\in \mathbb{Z}}$. Denote
\begin{align}
\label{eq:FttildeFt}
\mathcal{F}_{t}=(e_{t},e_{t-1},
\ldots )\quad \textrm{and}\quad \tilde{\mathcal{F}}_{t}=(\tilde{e}_{t},\tilde{e}_{t-1},\ldots
).
\end{align}
For the innovation part $\epsilon _{t}=(\epsilon_{t1},\ldots,\epsilon_{tN})^\top$ as in (\ref{modelstructural}),
the dependence
among $\epsilon _{t}:=(\epsilon _{1t},\ldots ,\epsilon _{Nt})^{\top }$
is allowed to be strong (e.g., with a factor structure) along the cross-sectional dimension and it has to be weak along the time
dimension, as we only exploit the sample size in the time series dimension
for the estimation of heterogeneous threshold effects. Specifically, we
assume that $\epsilon _{t}$ follows an MA($\infty $) process as follows:
\begin{equation}
\epsilon _{t}=\sum_{k\geq 0}A_{k}\eta _{t-k},
\end{equation}
where $\eta _{t}=(\eta _{1t},\ldots ,\eta _{\tilde{N}t})^{\top }\in \mathbb{R
}^{\tilde{N}},$ $t\in \mathbb{Z},$ are i.i.d. random vectors with zero mean
and identity covariance matrix $I_{\tilde{N}}$, and $A_{k}\in \mathbb{R}
^{N\times \tilde{N}}$ are real-valued matrices with $\tilde{N}\leq cN$ for
some constant $c>0$. In particular, we allow $\tilde{N}=N$ so that the $A_{k}$'s
become square matrices. We assume that $\left\{ e_{t},\tilde{e}
_{t}\right\} _{t\in \mathbb{Z}}$ are independent of $\left\{ \eta
_{t}\right\} _{t\in \mathbb{Z}}.$

\begin{remark}(Dependence between $X_{jt}$ and $U_{jt}$.)
 Note that Equation (\ref{eq:xjtujtdef}) is a fairly general setup, where
almost all kinds of dependence structures between $X_{jt}$ and $U_{jt}$ can
be fulfilled for the following simple reasoning.
To see this, take $U_{jt}$ as a scalar random variable as an example.
For any two continuous
random variables $X$ and $U$ with the conditional cumulative distribution
function of $U$ given $X$ as $F_{U|X}(\cdot |X)$. Let $\tilde{e}
=F_{U|X}(U|X) $. Then by Lemma F.1 in Appendix F in \cite{rio2017asymptotic}, $\tilde{e}$ is a uniformly distributed random variable on $[0,1]$ and is
independent of $X$. Moreover, $U=F_{U|X}^{-1}(\tilde{e}|X).$
\end{remark}

Next, we impose some conditions on the elements of $\eta _{t}.$

\begin{assumption}
{\normalsize (Moment) \label{asmp:moment} $\eta _{jt}$'s are i.i.d. with
zero mean and unit variance across $j$ and $t$. Either one of the following
two conditions is satisfied: }

\begin{enumerate}
\item[(i)] {\normalsize (Finite moments) For some constant $q>2,$ $\eta
_{11} $ has finite $q$th moment. }

\item[(ii)] {\normalsize (Sub-exponential) For some constant $\lambda _{0}>0$
, $\eta _{1}=(\eta _{11},\ldots ,\eta _{\tilde{N}1})^{\top }$ is
sub-exponential with $a_{0}:=\sup_{|v|_{2}\leq \lambda _{0}}\mathbb{E}
(e^{|v^{\top }\eta _{1}|})<\infty .$ }
\end{enumerate}
\end{assumption}


{\normalsize The following assumption imposes some conditions on the
temporal dependence for the processes $(X_{jt})_{t\in \mathbb{Z}}$ and $
(\epsilon _{jt})_{t\in \mathbb{Z}}.$ }
{Denote $\mathcal{C}_j$ as a set of search values for the time series $j$ containing at least the threshold values $c_{0j}$. If the true threshold is given, then $\mathcal{C}_j = \{c_{0j}\}$ is a single point in subsections \ref{known} and \ref{known2}. We will specify $\mathcal{C}_j$ later in the theorems for the unknown case.}
Let $g_{j}(\cdot )$ be the density function of $X_{jt}.$
Recall $\mathcal F_s$ in \eqref{eq:FttildeFt}, for $t\geq s+1,$ let $g_{j,t}(x|\mathcal{F}_{s}):=d\mathbb{P}(X_{jt}\leq x|
\mathcal{F}_{s})/dx.$

\begin{assumption}
{\normalsize (Dependence) \label{asmp:dependence}
(i) For any $x\in\mathcal C_j,$ assume $g_{j,t}(x|
\mathcal{F}_{t-1})$ has finite $p$th moment for some $p>2$.
Denote
\begin{equation*}
\theta _{k,p}=\Big\|\sup_{1\leq j\leq N,1\leq t\leq T,x\in\mathcal C_j}\big|g_{j,t}(x|
\mathcal{F}_{t-1})-g_{j,t}(x|\mathcal{F}_{t-1,\{t-1-k\}})\big|\Big\|
_{p}.
\end{equation*}
Assume $\sup_{m\geq 0}m^{\alpha }\sum_{k\geq m}\theta _{k,p}<\infty $ for
some $\alpha >0.$ }

{\normalsize (ii) For some $C>0$ and $\beta >0,$ we have $\max_{1\leq j\leq
N}\sum_{k\geq m}|A_{k,j,\cdot }|_{2}\leq C(m\vee 1)^{-\beta },\quad $where $
m\geq 0$ and $A_{k,j,\cdot }\in \mathbb{R}^{\tilde{N}}$ is the $j$th row of
matrix $A_{k}.$ }
\end{assumption}


 Assumption \ref{asmp:dependence} $(i)$-$(ii)$ essentially assumes
weak temporal dependence for the processes $(X_{jt})_{t\in \mathbb{Z}}$
and $(\epsilon _{jt})_{t\in \mathbb{Z}}$ for any $j\in \left[ N\right]$. We allow strong cross-sectional
dependence over the $j$ dimension. We shall also note that $\beta $
indicates the strength of the dependency structure of the error processes.


Furthermore, let $g_{j_{1},j_{2}}(\cdot ,\cdot |\mathcal{F}_{t-1})$ be
the conditional joint density of $X_{j_{1}t}$ and $X_{j_{2}t}$ given $\mathcal{F}_{t-1}.$ Finally, the following assumption imposes a condition on $
g_{j_{1},j_{2}}(\cdot ,\cdot |\mathcal{F}_{t-1}).$


\begin{assumption}
{\normalsize (Joint density) \label{asmp:bddxj1j2density} There is no
perfect or asymptotically perfect collinearity in $X_{j_{1}t}$ and $
X_{j_{2}t}$ for any $j_{1}\neq j_{2}.$ The conditional densities $
g_{j_{1},j_{2}}(x_{1},x_{2}|\mathcal{F}_{t-1})$ are uniformly upper bounded,
that is, $\max_{1\leq j_{1},j_{2}\leq N}\sup_{x_{1},x_{2}\in \mathbb{R}
}|g_{j_{1},j_{2}}(x_{1},x_{2}|\mathcal{F}_{t-1})|<\infty .$}
\end{assumption}

{\normalsize
}

{\normalsize
}

{\normalsize
}
Let
\begin{align}
\label{eq:defbarbunderb}
\bar{b}=\max_{1\leq j\leq N}b_{j}\quad\text{and}\quad \underline{b}
=\min_{1\leq j\leq N}b_{j}.
\end{align}
Assumptions \ref{boundedness}-\ref{kernel} in Section \ref{subsec:additionalassumption} in the Supplementary Material are concerning standard assumptions on boundedness, smoothness and kernel functions involved.

\subsection{\protect\normalsize Gaussian Approximation Results}

{\normalsize \label{known} }

{If the true threshold is given, then $\mathcal{C}_j = \{c_{0j}\}$ is a single point in Subsections \ref{known} and \ref{known2}.}
{We first consider the GA result for $\mathcal{I}$ defined in (
\ref{I1}). Define $d_{j}=(Tb_j)^{1/2}\gamma _{j}/v_{j}$ and $\underline{d}
=(d_{1},d_{2},\ldots ,d_{N})^{\top }$.
Further define $Z$ as a centered Gaussian random vector with identity covariance matrix.  Under the null }$H_{0}^{\left(
1\right) },$ the bias term $\underline{d}$
becomes a zero vector. The following theorem states that the infeasible statistic $\mathcal{I}$ can be approximated well by the maximum of Gaussian random
variables. Denote
\begin{equation*}
\Delta =\mathrm{log}^{7/6}(NT)(\underline{b}T)^{-1/6}+\mathrm{log}
^{1/2}(N)T^{1/2}\bar{b}^{5/2}+T^{-(\alpha p)\wedge (p/2-1)}+\bar{b}^{1/3}
\mathrm{log}^{2/3}(N).
\end{equation*}
Define }$\mathcal{R}_{NT}=\Delta +(\underline{b}T^{1-2/q})^{-1/3}\mathrm{log}
(NT)+\mathrm{log}(NT)N^{1/q}T^{-\beta }$ under Assumption 3.1$(i)$, and $
\mathcal{R}_{NT}=\Delta +\mathrm{log}(NT)T^{-\beta }$ under Assumption
3.1$(ii)$.


\begin{theorem}[Gaussian Approximation 1]
\label{thm:ga} Let Assumptions \ref{asmp:moment}-{\ref
{asmp:bddxj1j2density}} and \ref{boundedness}-\ref{kernel} hold. Then we have
\begin{equation}
\sup_{u\in \mathbb{R}}\big|\mathbb{P}(\mathcal{I}\leq u)-\mathbb{P}(|Z+
\underline{d}|_{\infty }\leq u)\big|\lesssim \mathcal{R}_{NT}.
\label{eq:ga1}
\end{equation}
\end{theorem}
Now we list the related rate requirement on $T$, $\bar b$ and $N$ under different moment conditions.
Under Assumption \ref{asmp:moment}$(i)$,
if $T\bar{b}^{5}\mathrm{log}(N)\rightarrow 0,$ $\underline{b}
^{-1}T^{2/q-1}\mathrm{log}^{3}(NT)\rightarrow 0$ and $\mathrm{log}
(NT)N^{1/q}T^{-\beta }\rightarrow 0,$ then $\Delta \rightarrow 0$ and $\mathcal{R}_{NT}\to 0$.
Under Assumption \ref{asmp:moment}$(ii)$, if $T\bar{b}^{5}\mathrm{log}(N)\rightarrow 0,$ $(T\underline{b})^{-1}
\mathrm{log}^{7}(NT)\rightarrow 0,$ $\bar b\mathrm {log}^3(N)\rightarrow0$ and $\mathrm{log}(NT)T^{-\beta
}\rightarrow 0,$ then $\Delta \rightarrow 0$ and $\mathcal{R}_{NT}\to 0$.

{\normalsize
}

{\normalsize
}

{\normalsize Theorem \ref{thm:ga} lays down the foundation for testing the
null hypothesis $H_{0}^{\left( 1\right) }.$ The proof of Theorem \ref{thm:ga}
is technically challenging and involved for two reasons. First, we allow for
both temporal and cross-sectional dependence among $\left\{ \left(
X_{jt},\epsilon _{jt}\right) \right\} .$ Second, the latent variable \(U_{jt}\) makes the error term \(e_{jt}\) in (\ref{mainmodel}) highly complex. Such complicated structure makes
it impossible to apply any existing GA results or their proof strategies
directly. In Section \ref{SecA.1} of the online
supplement, we outline the main idea used in the proof of the above theorem, highlighting that the key step is proving the GA result in Lemma \ref{lem:gapprox}. We note that the term $\mathrm{log}^{7/6}(NT)(\underline{b} T)^{-1/6} + \mathrm{log}^{1/2}(N)T^{1/2}\bar{b}^{5/2}$ arises from extending the original GA results in \cite{chernozhukov_central_2017}. The term $\mathrm{log}(NT)N^{1/q}T^{-\beta }$ is due to the temporal dependence and the moment condition of the {innovation term $\eta _{jt}$}. The term $\bar{b}^{1/3} \mathrm{log}^{2/3}(N)$ results from the covariance approximation and Gaussian comparison.
 } {Here we adopt maximum based test statistic, extending to other type of statistic can also be considered, see for example the statistics in   \cite{fan2015power} and \cite{giessing2020bootstrapping}.}



Next, we consider the GA result for {$\mathcal{Q}$} defined in (\ref{Q1}). Note that
\begin{align}
\mathcal{Q}& =\max_{1\leq j\leq N}(Tb_{j})^{1/2}\Big|
\sum_{t=1}^{T}\Big((w_{jt,b}^{+}-w_{jt,b}^{-})Y_{jt}-\sum_{j=1}^N(w_{jt,b}^{+}-w_{jt,b}^{-})Y_{jt}/N\Big)\Big|/
\tilde{v}_{j}\nonumber\\
&\approx \max_{1\leq j\leq N}\Big|(Tb_{j})^{1/2}
\sum_{t=1}^{T}\Big((w_{jt,b}^{+}-w_{jt,b}^{-})e_{jt}-\sum_{j=1}^N(w_{jt,b}^{+}-w_{jt,b}^{-})e_{jt}/N\Big)/
\tilde{v}_{j}+\tilde d_j\Big|,\label{Q2}
\end{align}
where $\tilde{d}_{j}=(Tb_{j})^{1/2}(\gamma _{j}-\frac{1}{N}\sum_{i=1}^{N}\gamma _{i})/
\tilde{v}_{j}.$ The approximation is valid due to the smoothness of $\tilde h_j(\cdot)$ and $\tilde \tau_j(\cdot)$, the remainder is addressed in Subsection \ref{SecB.3}.
Let $\tilde{\underline{d}}=(\tilde{d}_{1},\tilde{d}
_{2},\ldots ,\tilde{d}_{N})^{\top }.$
We deliver the GA results for $\mathcal{Q}$ in Supplementary Material \ref{homo} and the proofs in Section \ref{proofhomo}. Furthermore, the theoretical results regarding the derivative are given in Section \ref{derivatives} and the proofs in Section \ref{proofderi}. In addition,  in Section \ref{known2}, we establish the validity of the test statistics when variance estimators are utilized.


The most interesting finding is that in the GA result \eqref{eq:ga1}, neither temporal nor cross-sectional dependencies appear in $|Z+\underline{d}|_\infty$, despite allowing  strong cross-sectional dependence for the underlying data generating process $X_{jt}, U_{jt}, \epsilon_{jt}$ within $\mathcal{I}$. The theoretical intuition behind this result is that local smoothing effectively whitens the noise. In practice, this observation suggests that, when the Assumption \ref{asmp:bddxj1j2density} hold for cross-sectional dependence, the correlation in the estimated statistics can be safely ignored, requiring only adjustments for heteroskedasticity.


\subsection{Theoretical Results for the Unknown
Threshold Case\label{unknown_c2}}
\label{unknown} In this section, we present theoretical results in the case of unknown thresholds.
For simplicity of notation, we define $\mathcal{C}_j$ as the unified interval $\mathcal{C}_j = [c_{\min}, c_{\max}]$, where $c_{\min}$ and $c_{\max}$ are constants in $\mathbb{R}$. The theorem still holds even when $\mathcal C_j$ differs for different $j$. Recall that $[c_1,\ldots, c_K]$ is the grid we search over. Let $c_1=c_{\min}$ and  $c_{K}=c_{\max}.$ As our statistics aggregate over different threshold grids, the problem becomes more demanding.
Let $\Delta_{c,\min }=\min_{1\leq i\leq {K-1}}|c_{i+1}-c_{i}|$ and $\Delta_{c,\max}=\max_{1\leq i\leq {K-1}}|c_{i+1}-c_{i}|$.

Recall $w_{jt,b}(c_{i})$ in \eqref{eq:wjtbci} and define ${\gamma }_{j}(c_{i})=\sum_{t=1}^{T}w_{jt,b}(c_{i})
\gamma _{j}$. Let $$d_{jK}=(\gamma _{j}(c_{1})/v_{j}(c_{1}),\ldots ,\gamma
_{j}(c_{K})/v_{j}(c_{K}))^{\top }$$ and $\underline{d}
^{C}=(d_{1K}^{\top },d_{2K}^{\top },\ldots ,d_{NK}^{\top })^{\top }.$
Define $Z^{C}=\left( {{Z_{1}^{C\top },...,Z_{N}^{C\top }}}\right) ^{\top }$
as an $NK$-vector of mean zero Gaussian distributed random vector with
covariance matrix $\Sigma^C$
specified below. For each $j\in \lbrack N]$ and $1\leq i_{1},i_{2}\leq K$,
the covariance between the
$i_{1}$- and $i_{2}$-th elements of $Z_{j}^C$  {{is
given by{
\begin{equation*}
\mathrm{Cov}(Z_{j,i_1}^C, Z_{j,i_2}^C)=(v_{j}(c_{i_{1}})v_{j}(c_{i_{2}}))^{-1}(Tb_{j})
\sum_{t=1}^{T}w_{jt,b}(c_{i_{1}})w_{jt,b}(c_{i_{2}})\mathrm{Var}\big(
e_{jt}|X_{jt}\big).
\end{equation*}
}}}${\normalsize {{Z_{j_1}^C}}}$ and $Z_{j_2}^C$ are
independent for $ j_1\neq j_2$,  so that the matrix $\Sigma^C$ is block diagonal. Note that $
\underline{d}^{C}$ denotes a bias term which is a zero vector under the null hypothesis. In addition, when $2\max_{1\leq j\leq
N}b_{j}<\Delta_{c,\min },$ one can ensure the matrix to be diagonal. The following theorem establishes the
asymptotic property of our test statistics with unknown thresholds.

\begin{theorem}\label{thm:unknown}
{\normalsize (Gaussian Approximation) Suppose that Assumptions \ref
{asmp:moment}-\ref{asmp:bddxj1j2density} and \ref{boundedness}-\ref{kernel} hold. Suppose that }$\Delta_{c,\max }\rightarrow 0.$ Then
\begin{equation}
\sup_{u\in \mathbb{R}}\big|\mathbb{P}(\mathcal{I}^{C}\leq u)-\mathbb{P}
(|Z^C+\underline{d}^C|_{\infty }\leq u)\big|\lesssim \mathcal{R}_{(NK)T},
\label{eq:ga}
\end{equation}
where $\mathcal{R}_{(NK)T}$ is $\mathcal{R}_{NT} $ in Theorem
3.1 with $N$ replaced by $NK$.
\end{theorem}

The proof of the above theorem is a generalization
of Theorem \ref{thm:ga} with a different variance-covariance structure. It should be noted that the conclusion of the above theorem follows the same pattern as Theorem \ref{thm:ga}, except that $N$ is replaced with $NK$.
The following theorem provides the theoretical support for the consistency
of the threshold estimator. Note that we only focus on the location $j$ where $\gamma_j$ does not equal to $0$, and we require those breaks to be significant.

\begin{theorem}[Consistency of the estimation of thresholds]
{\normalsize \label{consistency} Suppose that Assumptions \ref{asmp:moment}-
\ref{asmp:bddxj1j2density} and \ref{boundedness}-\ref{kernel}  hold. Assume $\min_{1 \leq j\leq N: \gamma_j \neq 0}(Tb_{j})^{1/2}|\gamma _{j}|\gg 1\ $ and }$\Delta_{c,\max }\rightarrow 0.$
Suppose $c_{0j}\in\{c_i, 1\leq i\leq K\}$, then
\begin{equation}
\max_{1\leq j\leq N,\gamma_j\neq 0} \gamma _{j}^{2}|\hat{c}_{j}-c_{0j}|=O_{\mathbb P}((\mathrm {log} N)/T).
\end{equation}

\end{theorem}

If the true break $c_{0j}$ does not fall in the grid, then the difference is controlled by $\max\{\Delta_{c,\max},  \mathrm {log}(N)/(T\gamma_j^2)\}.$ Under the assumption of a minimum break signal, we can achieve
consistency in detecting the break. The rate is expected to depend
on both the signal strength $\gamma_j$ and the available sample size $T$.
Note that the precision of threshold estimation is not affected by the bandwidth $b_j$. This rate is in line with the high-dimensional changepoint literature, see for example Theorem 2 in \cite{li2024}.

{
To further evaluate the feasibility of the test statistic $\mathcal{\hat{I}}^{C}$, we will also examine the consistency of the variance estimator $\hat{\sigma}_{e,j}^2(c_i)$ defined in Section \ref{subsec:unknown} of the Supplementary Material. Estimating the variance in the presence of an unknown threshold requires additional caution, particularly when accounting for jumps, see Section \ref{sec:algo}. We will consider a truncation operation to eliminate the contamination from jumps. The consistency of the estimator will be established in a similar manner.  For completeness, we provide Theorems \ref{thm:hatvap}  and \ref{thm:powerap} in Section \ref{subsec:unknown} of the Supplementary Material as counterparts to Theorems \ref{thm:hatv} and \ref{thm:power}, specifically addressing the case with an unknown threshold.
}

{\normalsize
}

{\normalsize
}

{\normalsize
}



\section{{\protect\normalsize Monte Carlo Simulations \label{sec:simulation}}
}

In this section, we analyze the finite sample performance of
our GA results. In particular, we study the size and
power of our uniform testing procedures in different settings. The data generating processes(DGPs) 1-7 are presented in Section \ref{sec:simulation_appendix} in Supplementary Materials.


\subsection{\protect\normalsize Simulation Results}
Throughout the study, we use a local linear estimator with a uniform kernel and select the bandwidth according to the procedure described in \citet{calonico2014robust}. We consider three
significance levels, namely, $\alpha =0.1,$ $0.05,$ and $0.01$. The results
are based on $1000$ Monte Carlo iterations. All tables and figures are presented in Section \ref{sec:simulation_appendix} in Supplementary Materials.

{\normalsize
The results of simultaneous jump effects testing for DGP1--DGP3 are displayed in Table \ref{table:sim1}.
First, the empirical size closely matches the nominal size across all $N$ values and for large $T$ in all DGPs. For DGP3, we observe slightly
oversize when $N=10\ $ which may be caused by the presence of
heteroskedasticity in comparison with results with DGP2. Overall, these
simulation results confirm our theoretical findings in Theorem \ref{thm:ga}.
In particular, it shows that the test procedure, which is based on
critical values neglecting the dependence in the covariance structure, has fairly accurate size control
even in the case of strong cross-sectional dependence in the innovations. Second, the local
power properties of the test are reasonably well for all DGPs under investigation, and
the power increases as either $N$ or $T$ increases. This result confirms that our
testing procedure is well-suited to detect sparse alternatives.


In addition, we achieve comparable simulation performance in terms of size and power results for the test of \emph{homogeneity} of the threshold effects in Table \ref{table:sim2}. The results confirm our Gaussian approximation result in Theorem \ref{thm:homotesting}.
For robustness, we also consider three additional settings with even stronger cross-sectional dependence for DGP4--6 in Table \ref{table:sim3}. The results show that size and power are indeed robust to strong cross-sectional dependence. Furthermore, Table \ref{table:sim4} demonstrates that detecting the existence of jump effects simultaneously across cross-sectional dimensions in the case of unknown thresholds has size and power properties as good as in the known threshold case, confirming our theoretical results from Theorem \ref{thm:unknown}. We additionally show the good finite sample performance for the estimation accuracy of the threshold location estimator in Table \ref{table:sim5}, providing Monte Carlo evidence supporting the consistency result of the threshold in Theorem \ref{consistency}.}


\label{compare}For comparison, one direct method is the sup-Wald or sup-likelihood ratio test from the pooled threshold regression literature, neither of which accounts for nonlinearity or heterogeneity (see, for example, \cite{andrews1993tests}, \cite{hansen2000sample} and \cite{fong2017chngpt}). Notably, recent studies, such as \cite{barassi2023threshold}, propose linear models with heterogeneous threshold effects in panel settings. However, they do not provide uniform testing methods.
To facilitate comparison, we firstly restrict ourself to linear threshold model DGP7.
And under the alternative we sample a fraction of coefficients $\gamma_j$ from a standard normal distribution. As shown in Table \ref{table:sim6}, our method has similar performance as the sup-Wald type statistic under the null hypothesis, while it has much better power under the alternative. See also Figure \ref{figure:comparison} in Section \ref{sec:model}. The reasons for this are twofold. First, our uniform testing procedure is robust towards settings with sparse signal, in which a pooled test suffers from dilution of signals. And second, our uniform testing procedure does not have the problem of signal cancellation in case the signals for different $j$'s have opposite signs and cancel each other out. As a second comparison, we consider the nonparametric pooled test statistic under DGP3. Again, we draw the non-zero threshold coefficients from a standard normal. Such a setup without accounting for parameter heterogeneity is covered such as in \citet{spokoiny1998estimation}, \citet{yang2014jump} and \cite{chiou2018nonparametric}. As shown in Table \ref{table:sim7}, similar to the linear case, our method works with better power. The above superior performance is mainly due to the sparsity of the signal under the alternative hypothesis. In addition, our method has better power even for dense alternative comparing to the pooled type test statistic when there might exist signal cancellation as illustrated in Figure \ref{figure:comparison}.

We also consider the test of existence of a threshold effect in the derivative under DGP1--3, to achieve similar size and power performance, oversmoothing is needed, as shown in Table \ref{table:sim8}.




\section{Empirical Applications \label{sec:application}}

In this section, we apply our testing procedure to two empirical applications. First, we study possible jumps in the news impact curve for constituents of the S\&P 500. As a second application, we study the impact of the incumbency effect in U.S. House elections.


\subsection{Application I: Jumps in the News Impact Curve}
\label{application1}
In this application, we revisit the news impact curve introduced by \citet{engle1993measuring}. In particular, we want to identify possible jumps in the relationship between news events in the form of past returns on the volatility of stocks. The application serves as an illustration of our uniform testing procedure when the threshold $c_{0j}$ is unknown. Originally proposed in the context of autoregressive volatility models (e.g., ARCH), \citet{engle1993measuring} and \citet{linton2005estimating} consider partially nonparametric specifications for the news impact curve which are able to capture possible asymmetries. A limitation of these specifications is the restriction to continuous functions. However, \cite{zakoian1994threshold} and \cite{fornari1997sign} advocate for (parametric) GARCH specifications that allow for discontinuities in the way past returns impact volatility. We therefore want to investigate this issue by testing for the presence of jumps.

In our model setup, the dependent variable $Y_{jt}$ is the Garman-Klass volatility of stock $j$ at time $t$, and $X_{jt}$ is the lagged stock return.
We want to test the null hypothesis $H_0^{(1)}:\gamma_i=0$ for all $i\in[N]$ and all $c_{0j}$, versus the alternative $H_a^{(1)}:\gamma_i\neq0$ for some $i$ and some $c_{0j}$.  Our testing procedure is robust to the presence of unobserved risk factors, $U_{jt}$, which are more than likely to occur in our application. Our data includes the daily stock returns and volatility data of the S\&P 500 constituents from January 2010 to May 2024, which we obtained from Kaggle \footnote{https://www.kaggle.com/datasets/andrewmvd/sp-500-stocks.}. Therefore we have a fairly high-dimensional setting with $N=500$ and $T_j$ is the number of days of stock $j$ being a member of the stock index. See Figure \ref{figure:return_vola} for a visualization of the time series of returns and volatility for the stock of AT\&T. As candidates for the threshold location, we consider $c_{i}\in\{-0.01,-0.005,0,0.005,0.01\}$.

\begin{figure}[tbp]
\centering
\includegraphics[width=.49\textwidth]{return.pdf}\hfill
\includegraphics[width=.49\textwidth]{volatility.pdf}
\caption{Time series for daily stock returns (left panel) and daily volatility (right panel) for AT\&T. The dashed lines in the left panel indicate the $10\%$ and $90\%$ quantiles of the return distribution.}
\label{figure:return_vola}
\end{figure}

\begin{figure}[tbp]
\centering
\includegraphics[width=.49\textwidth]{significant.pdf}\hfill
\includegraphics[width=.49\textwidth]{insignificant.pdf}
\caption{Local linear fit of a significant threshold effect for AT\&T (left panel) and an insignificant effect for Akamai (right panel).}
\label{figure:comparison_application}
\end{figure}



The main findings are summarized in Table \ref{table:stocks}.
We find significant results for $11$ out of the $500$ stocks we consider. For stocks with significant threshold effects, the locations of the estimated thresholds vary, and none of the significant threshold effects occur precisely at zero. This contradicts the rationale of the sign-switching ARCH model of \cite{fornari1997sign}, which models a jump effect based on the sign of the lagged returns, and thus expects a significant threshold effect at zero. For $c_{i}=-0.01$, we get significant negative effects at the $1\%$ level for five stocks, including Allstate Corporation and American Express. On the other hand, we obtain significant positive effects for the $c_{i}=0.01$ threshold for four stocks, e.g.  AT\&T and McKesson. These results reaffirm the findings of \cite{chen2011news}, who find that both
bad and good news have an increasing effect on volatility. Additionally, we get another negatively significant coefficient at the $-0.005$ threshold location and another positively significant coefficient at $c_{i}=0.005$. Overall, we get evidence for the presence of jump effects in the news impact curve. Figure \ref{figure:comparison_application} shows a visualization of the significant positive threshold effect for AT\&T at $c_i=0.01$, and as a comparison we show the fit for Atamai for which we did not find a significant effect at any threshold location.

Furthermore, we are interested in studying the predictive advantage of our method compared to a baseline model without jumps. The baseline model can be written as $Y_{jt}=\tilde{h}_j(X_{jt})+e_{jt}$, where $\tilde h_j(\cdot)$ and $e_{jt}$ are defined in \eqref{mainmodel}.
For this purpose, we split the data into two parts, with training data from 2010--2019 and test data from 2020--2024. We restrict our analysis to the stocks identified in Table \ref{table:stocks} with significant threshold effects at $1\%$ level. The mean squared error of prediction is lower for all but one of the 11 stocks we consider. The same is true if we restrict our analysis to these observations close to the respective threshold, $c_{0j}$. We conclude that including threshold effects can slightly improve the predictive accuracy with a median gain of $1.6\%$ in predictive power.

\subsection{Application II: Party Incumbency Effects on U.S. House Elections}\label{application2}

In the second application, we are interested in estimating the advantage a
candidate has in the U.S. House elections when the seat is currently
occupied by the party of the candidate \footnote{We use the data provided by the \citet{electiondata}. The data include the
vote shares of all candidates in the elections for the US House of
Representatives from 1976--2020. Due to redistricting at the start of each
decade, we have to exclude years ending with a `0' or `2'.}. This is called the party incumbency
effect. We follow the research design of \citet{lee2008randomized}. See
\citet{caughey2011elections} for an analysis of more recent elections. The
dependent variable $Y_{jt}$ is the democratic two-party vote share in year $
j $ in district $t$. The covariate $X_{jt}$ is the difference in the
two-party vote share in the previous election. Since the winner of the
election is determined by the 0 threshold, our model setup is appropriate.
We want to test the one-sided null hypothesis $H_{0}^{(1)}:\gamma _{i}<0$ for all $
i\in \lbrack N]$, versus the alternative hypothesis $H_{a}^{(1)} $: $\gamma
_{i}\geq 0$ for some $i$. Unlike \citet{lee2008randomized} who pools the
data together to estimate a single incumbency effect, we conduct a separate
analysis for different election years and also a separate analysis for
different states by interchanging the roles of $j$ and $t$ and allowing for
possible heterogeneity of incumbency effects along either the $j$ or $t$
dimension. In addition, it is easy to see our theories extend to the
unbalanced data in which case we can use $T_{j}$ to denote the number of
observations associated with individual/group $j$.

First, we are interested in the election year-specific effects. The election
year-specific results are reported in Table \ref{table:vote_new} in Appendix \ref{application2_appendix}. Our
analysis covers 13 election years, so $N=13$ and $T_{j}$ is equal to the
number of districts included in the sample of the respective election
year/group. For each $j,$ the number of effective observations is given by
the number of non-zero weights in the local linear estimation, which depends
on the bandwidth parameters $b_{j}$'s which are chosen as in the simulation
section. Table \ref{table:vote_new} reports the results for the effect
estimate ($\hat{\gamma}_{j}$)$,$ the standard
error ($(T_{j}b_{j})^{-1/2}\hat{v}_{j}$)$,$ the individual test statistic ($(T_{j}b_{j})^{1/2}
\hat{\gamma}_{j}/\hat{v}_{j}$)$,$ the number of observations $T_{j}$ (Obs)
and the number of effective observations (Eobs) as well. We find five years
with significant effects at the $1\%$ level and one year which is significant at the
$5\%$ level by using the simulated critical value for $\mathcal{\hat{I}}$.
Since we are interested in the maximum of the test statistics in our uniform
testing procedure, the null hypothesis of no effect ($H_{0}^{\left( 1\right)
}$) can be rejected at the $1\%$ level. We can observe that the estimated
incumbency effect is stronger in the period from 1978 to 1998, while in the
most recent elections the effect is insignificant. However, the
sign of the estimated coefficients is positive for all election years. We
also run the test for the homogeneity of incumbency effects ($H_{0}^{\left(
2\right) }$) over the 13 election years, see Table \ref
{table:vote_new_hetero}. We rely on the median ($\hat{\gamma}_{0.5}$) as a
more robust estimator of the average effect across election years. The test statistic is given by $\hat{\mathcal{
Q}}=3.011$, which implies that we can reject the null hypothesis of
homogeneous effects at a 5\% level.


Now, we are interested in the state-specific effects. We restrict our
attention to states with a sufficient number of elections held in the
sample. We choose a threshold of 50 elections. This gives us $N=29$ states
in our analysis, and $T_{j}$ now denotes the number of elections included
for the $j$th state. We argue that this setting constitutes a fairly high-dimensional setting since the effective number of observations is smaller
than $N$ for a fairly large proportion of the states. The results are
reported in Table \ref{table:vote_new_states}. Consistent with the previous
results, almost all of the estimated coefficients are positive. By using the
critical values for the simultaneous tests, we find significant effects for
six states, for Missouri at $1\%$ level, and for
California, Iowa, Virginia, Alabama and South Carolina at $5\% $. Again, we
can reject the null hypothesis of no jump effect for our uniform testing
procedure at a $1\%$ level. The significant states tend to belong to those
with the largest number of observations. However, for the other two largest
states, New York and Texas, the estimated effect is rather small and far
from being significant. The incumbency effect seems to be stronger in some
states compared to others. See Figure \ref{figure:application_election} in the Supplement. for a visualization of the estimated threshold effects for the states of California and Pennsylvania. We also run the test for the homogeneity of
incumbency effects ($H_{0}^{\left( 2\right) }$) over the 29 states. The
results are displayed in Table \ref{table:vote_new_states_bw_hetero} and the
test statistic is given by $\hat{\mathcal{Q}}=3.011$, which implies that we
can reject the null hypothesis of homogeneous effects at the $5\%$ level. We
find heterogeneity in the effects across election periods and states.


\section{Conclusion}
\label{sec:conclusion}
This paper focuses on the estimation and inference of heterogeneous threshold effects in nonparametric high-dimensional data, accounting for both cross-sectional and time dependencies in the covariates and error terms. We propose a test to assess the significance of the threshold effect and another to evaluate whether the threshold effects are homogeneous across individuals or groups. Additionally, we introduce a consistent method for estimating unknown threshold locations. Our tests are based on high-dimensional Gaussian approximation results. Despite the complexity of the underlying dependency structures, we show that the variance-covariance structure of the threshold effect estimators has a simple analytical expression under general dependence conditions. Unlike changepoint analysis in time series, our approach considers a setup involving latent variables, leading to a significantly different development of the GA results. Simulations demonstrate significant power enhancement compared to the existing literature under various settings.

{\small
\bibliographystyle{apalike}
\bibliography{literature.bib}
}

{\normalsize
}

{\normalsize
}

{\normalsize
}

{\normalsize
}

{\normalsize
}

{\normalsize
}

{\normalsize
}

{\normalsize
}

{\normalsize
}

{\normalsize
}

{\normalsize
}

{\normalsize
}

{\normalsize
}

{\normalsize
}

{\normalsize
}

{\normalsize
}

{\normalsize
}

{\normalsize
}

{\normalsize
}

{\normalsize
}

{\normalsize
}

{\normalsize
}

{\normalsize
}

{\normalsize
}

{\normalsize
}

{\normalsize
}
\newpage

\setlength{\baselineskip}{12pt}
\spacingset{1.5}

\bigskip
\begin{center}
{\large\bf SUPPLEMENTARY MATERIALS}
\end{center}

\tableofcontents