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.
86,667 characters
Factor-Driven Two-Regime Regression
\doparttoc
\faketableofcontents
\title{Factor-Driven Two-Regime Regression
\thanks{We would like to thank Don Andrews, Mehmet Caner, Greg Cox, Bruce Hansen, Zhongjun Qu and the seminar participants at BU, Emory, Michigan State, NYU, Wisconsin-Madison, Northwestern, Yale, and 2018 ASSA Winter Meeting for helpful comments. We would like to thank
the Ministry of Education of the Republic of Korea and the National Research Foundation of Korea (NRF-0405-20180026),
the Social Sciences and Humanities Research Council of Canada (SSHRC-435-2018-0275),
the European Research Council for financial support (ERC-2014-CoG-646917-ROMIA) and
the UK Economic and Social Research Council for research grant (ES/P008909/1) to the CeMMAP.}
}
\date{August 13, 2020}
\author{
Sokbae Lee\thanks{
Address: 420 West 118th Street, New York, NY 10027, USA. E-mail: \texttt{[email removed]}.} \\ \footnotesize Columbia University \and
Yuan Liao\thanks{Address: 75 Hamilton St., New Brunswick, NJ 08901, USA. Email:
\texttt{[email removed]}.}\\ \footnotesize Rutgers University
\and Myung Hwan Seo\thanks{
Address: 1 Gwanak-ro, Gwanak-gu, Seoul 08826, Korea. E-mail:
\texttt{[email removed]}.} \\ \footnotesize Seoul National University\and
Youngki Shin\thanks{Address: 1280 Main St.\ W.,\ Hamiloton, ON L8S 4L8, Canada. Email:
\texttt{[email removed]}.}\\ \footnotesize McMaster University
}
\maketitle
\begin{abstract}
We propose a novel two-regime regression model where regime switching is driven by a vector of possibly unobservable factors. When the factors are latent, we estimate them by the principal component analysis of a panel data set.
We show that the optimization problem can be reformulated as mixed integer optimization, and we present two alternative computational algorithms. We derive the asymptotic distribution of the resulting estimator under the scheme that the threshold effect shrinks to zero. In particular, we establish a phase transition that describes the effect of first-stage factor estimation as the cross-sectional dimension of panel data increases relative to the time-series dimension.
Moreover, we develop bootstrap inference and illustrate our methods via numerical studies.
\\
\\
Keywords: threshold regression, principal component analysis, mixed integer optimization, phase transition, oracle properties
\end{abstract}
\thispagestyle{empty}
\onehalfspacing
\newpage
\setcounter{page}{1}
\pagenumbering{arabic}
\section{Introduction}
Suppose that $y_t$ is generated from
\begin{align}
y_{t} & = x_{t}^{\prime }\beta _{0}+x_{t}^{\prime }\delta_{0}1\{f_t'\gamma_0>0\}+\varepsilon _{t}, \label{model1} \\
\mathbb{E}\left( \varepsilon _{t}|\mathcal{F}_{t-1}\right) & = 0, \; t =1, \ldots, T, \label{model2}
\end{align}
where $x_{t}\ $and $f_{t}$ are adapted to the filtration $\mathcal{F}_{t-1}$,
$(\beta_0, \delta_0, \gamma_0)$ is a vector of unknown parameters,
and the unobserved random variable $\varepsilon _{t}$ satisfies the conditional mean restriction
in \eqref{model2}. We interpret $f_t$ to be a vector of {factors} determining regime switching. When $f_t'\gamma_0 > 0$, the regression function becomes
$x_{t}^{\prime }(\beta _{0}+\delta_0)$; if $f_t'\gamma_0 \leq 0$, it reduces to $x_{t}^{\prime }\beta _{0}$.
We allow for either observable or unobservable factors.
For the latter, we assume that they can be recovered from a panel data set.
In light of this feature, we call the model in \eqref{model1} and \eqref{model2}
a \emph{factor-driven two-regime regression model}.
Our paper is closely related to the literature on threshold models with unknown change points
(see, e.g., \citep{chan1993consistency}, \citep{hansen2000sample}, \citep{ling_1999},
\citep{Seijo:Sen:11a}, \citep{Seo-Linton}, and \citep{tong1990non}, among many others).
In the conventional threshold regression model,
an intercept term and a scalar observed random variable constitute $f_t$.
For instance, Chan \citep{chan1993consistency} and Hansen \citep{hansen2000sample} studied the model in which $1\{f_t'\gamma_0>0\}$ in (\ref{model1}) is replaced by $1\{q_t> \widetilde{\gamma }_0 \}$ for some observable scalar variable $q_t$ with a scalar unknown parameter $\widetilde{\gamma}_0$.
In practice,
it might be controversial to choose which observed variable plays the role of $q_t$. For example, if the two different regimes represent
the status of two environments of the population, arguably it is difficult to assume that the change of the environment is governed by just a single variable.
On the contrary, our proposed model introduces a regime change due to a single index of factors that can be ``learned" from a potentially much larger dataset. Specifically, we consider the framework of latent approximate factor models in order to model a regime switch based on a potentially large number of covariates.
In view of the conditional mean restriction
in \eqref{model2}, a natural strategy to estimate $(\beta_0, \delta_0, \gamma_0)$ is to rely on least squares.
A least-squares estimator for our model brings new challenges in terms of both computation and asymptotic theory.
First of all, when the dimension of $f_t$ is larger than 2, it is computationally demanding to estimate $(\beta_0, \delta_0, \gamma_0)$. We overcome this difficulty by developing new computational algorithms based on the method of mixed integer optimization (MIO).
See, for example, section~2.1 in Bertsimas et al.\ \citep{bertsimas2016} for a discussion on computational advances in solving the MIO problems.
Second, we establish asymptotic properties of our proposed estimator by adopting a diminishing thresholding effect.
That is, we assume that $\delta_0=T^{-\varphi}d_0$ for some unknown $\varphi \in (0, 1/2)$ and unknown non-diminishing vector $d_0$.
The diminishing threshold has been one of the standard frameworks in the change point literature (e.g., \citep{bai1994least,hawkins1986simple,horvath1997effect}).
The unknown parameter $ \varphi $ reflects the difficulty of estimating $ \gamma_0 $ and affects the identification and estimation of the change-point $ \gamma_0 $.
Both the rate of convergence and the asymptotic distribution depend on $\varphi$.
This is a widely employed tool to allow for flexible signal strengths of the parameters in the nonlinear model.
For instance, McKeague and Sen \citep{mckeague2010fractals} studied a ``{point impact}" linear model, where the identification and estimation of $\gamma_0$ are affected by an unknown slope $\delta_0$.
While specifically assuming $\delta_0\neq0$, they encountered a similar parameter $\varphi$, reflecting the difficulty of estimating $\gamma_0$.
The asymptotic theory for the estimated $\delta_0$ under the diminishing jump setting is fundamentally different from the fixed jump setting: the former is determined by a Gaussian process (e.g., \citep{hansen2000sample}), and the latter by a compound poison process (e.g., \citep{chan1993consistency}). While both settings lead to important asymptotic implications, we focus on the diminishing setting because when the factors are estimated, there is a new and interesting \emph{phase transition} phenomenon that smoothly appears in the ``bias'' term of the Gaussian process. The phase transition characterizes the continuous change of the asymptotic distribution as the precision of the estimated factors increases relative to the size of the jump, which we shall detail below.
When the factor $f_t$ is latent, we estimate it using principal component analysis (PCA) from a potentially much larger dataset, whose dimension is $N$. It turns out that the asymptotic distribution for the estimator of $\alpha_0 \equiv (\beta_0',\delta_0')'$ is identical to that when $\gamma_0$ were known, regardless of whether factors are directly observable or not; therefore, the estimator of $\alpha_0$ enjoys an oracle property.
The issue is more sophisticated for the distribution of the estimator of $\gamma_0$. When factors are directly observable, we prove that
\begin{align*}
& T^{1-2\varphi }
\left( \widehat{\gamma }-\gamma _{0}\right)\overset{d}{\longrightarrow }
\operatorname*{argmin}_{g\in \mathcal{G}}B(g) +2W\left( g\right) ,
\end{align*}
where $B(g)$ represents a ``drift function" of the criterion function, which is linear with a kink at zero,
$W(g)$ is a mean-zero Gaussian process
and $\mathcal{G}$ is a rescaled parameter space. However, when factors are not directly observable, the estimation error from the PCA plays an essential role and may slow down the rates of convergence, depending on the relation between $N$ and $T$. Specifically, we show that
\begin{align*}
& \left( \left( NT^{1-2\varphi }\right) ^{1/3}\wedge T^{1-2\varphi }\right)
\left( \widehat{\gamma }-\gamma _{0}\right)\overset{d}{\longrightarrow }
\operatorname*{argmin}_{g\in \mathcal{G}}A\left( \omega, g\right) +2W\left( g\right),
\end{align*}
with a new drift function $A\left( \omega, g\right)$ that depends on $\omega=\lim \sqrt{N} T^{-(1-2\varphi)}\in[0,\infty]$. On one hand, when $\omega=\infty$, we find that $A(\omega, g)= B(g)$, so the limiting distribution becomes the same as if the factors were observable. This case
corresponds to
the \emph{super-consistency rate} (e.g., \citep{hansen2000sample}).
On the other hand, when $\omega=0$, it turns out that $A\left( \omega, g\right)$ is quadratic in $g$, corresponding to a \emph{cube root rate} similar to the maximum score estimator (e.g., \citep{kim1990, seo2018local}).
Furthermore, both the drift function and the resulting rates of convergence have continuous transitions as $\omega$ changes between
$0$ and $\infty$.
Therefore, one of our key findings for the estimator of $\gamma_0$ is the occurrence of a {phase transition} from a \emph{weak-oracle} limiting distribution to a \emph{semi-strong} oracle one, and then to a \emph{strong} oracle one as $\omega$ increases.
As the asymptotic distribution of $\widehat\gamma$ is non-pivotal, we propose a wild bootstrap for inference of $\gamma_0$. Importantly, we construct bootstrap confidence intervals for $\gamma_0$ that do not require knowledge of $\varphi$.
This facilitates applications in which the jump diminishing speed is not known in advance.
The remainder of the paper is organized as follows.
In Section \ref{sec:opt}, we propose the least-squares estimator and algorithms to compute the proposed estimator.
In Section \ref{sec:known factors}, we establish asymptotic theory when $f_t$ is directly observed.
In Section \ref{model:est:factors}, we consider estimation when $f_t$ is a vector of latent factors, we propose a two-step estimator via PCA, and we analyze asymptotic properties of our proposed estimator. In Section \ref{sec:inference}, we develop bootstrap inference, and in Section \ref{sec:MC} we give the results of Monte Carlo experiments. In Section \ref{sec:real-data-app}, we illustrate our methods by applying them to threshold autoregressive models of unemployment. We conclude in Section \ref{sec:conclusions}. The online appendices provides details that are omitted from the main text.
The notation used in the paper is as follows. The sample size is denoted by $T$ and the transpose of a matrix is denoted by a prime.
The true parameter is denoted by the subscript $0$, whereas a generic element has no subscript.
The Euclidean norm is denoted by $| \cdot |_2$,
the Frobenius norm of a matrix is $| \cdot |_F$, the spectral norm of a matrix is $|\cdot|_2$,
and the $\ell_0$-norm is $|\cdot |_0$.
For a generic random variable or vector $z_{t}$, let its density function be
denoted by $p_{z_{t}}$. Similarly, let
$p_{y_t|x_t}(y)$ denote the conditional density of $y_t$ given $x_t$ for the random vectors $y_t $ and $x_t$.
The abbreviation \emph{a.s.} means almost surely.
\section{Least-Squares Estimator via Mixed Integer Optimization}\label{sec:opt}
\subsection{Identifiability}
We use the convention that the constant $ 1 $ is the first element of $ x_t $ and $ -1 $ is the last element of $ f_t $. Define $\alpha :=(\beta ^{\prime },\delta ^{\prime
})^{\prime }$ and $Z_{t}(\gamma ):=(x_{t}^{\prime },x_{t}^{\prime }1\{f_{t}^{\prime
}\gamma >0\})^{\prime }$. Then, we can rewrite the model as
\begin{equation*}
y_{t}=Z_{t}\left( \gamma _{0}\right) ^{\prime }\alpha _{0}+\varepsilon _{t}.
\end{equation*}
Because only the sign of the index $f_t'\gamma_0$ determines the regime switching,
the scale of $\gamma_0$ is not identifiable.
We assume that the first element of $\gamma_0$ equals 1.
Let $d_x$ and $d_f$ denote the dimensions of $x_t$ and $f_t$, respectively.
\begin{assum}
\label{scale-normalization}
$\alpha_0 \in \mathbb{R}^{2d_x}$ and
$\gamma_0
\in \Gamma :=
\{ (1, \gamma_2')': \gamma_2 \in \Gamma_2 \}$, where
$\Gamma_2 \subset \mathbb{R}^{d_f-1}$ is a compact set.
\end{assum}
We decompose $f_{t}$ into a scalar random variable $f_{1t}$ and other variables $f_{2t}$, so that $f_{t}^{\prime}\gamma \equiv f_{1t} + f_{2t}^{\prime} \gamma_2$.
In view of the conditional mean zero restriction in \eqref{model2}, it is natural to impose
conditions under which both $\alpha _{0}$ and $\gamma _{0}$ are
identified by the $L_{2}$-loss. Introduce the excess loss
\begin{align}\label{eq:excess}
R(\alpha ,\gamma ) := \mathbb E(y_{t}-x_{t}^{\prime }\beta -x_{t}^{\prime }\delta
1\{f_{t}^{\prime }\gamma >0\})^{2}- \mathbb E( \varepsilon _{t}^{2} ).
\end{align}
In order to establish that
$R\left( \alpha,\gamma \right) >
R\left( \alpha _{0},\gamma _{0}\right) =0$ whenever
$(\alpha,\gamma) \neq (\alpha_0,\gamma_0)$,
we make the following regularity conditions.
\begin{assum}
\label{iden-assump}
For any $\varepsilon >0$,
$\left( \alpha_{0},\gamma _{0}\right) $ satisfies
\begin{equation*}
\inf_{\{(\alpha', \gamma')' \in \mathbb{R}^{2d_x} \times \Gamma: |(\alpha', \gamma') - (\alpha_0', \gamma_0')|_2 > \varepsilon \} }
R\left( \alpha,\gamma \right) >0.
\end{equation*}
\end{assum}
Online Appendix \ref{sec:identification} provides sufficient conditions for Assumption \ref{iden-assump}.
\subsection{Estimator}
We now propose the least-squares estimator and two alternative algorithms to compute the proposed estimator.
For computational purposes, we assume that $\alpha \in \mathcal{A} \subset \mathbb{R}^{2d_x}$ for some known compact set $\mathcal{A}$.
In practice, we can take a large $2 d_x$-dimensional hyper-rectangle so that the resulting estimator is not on the boundary of $\mathcal{A}$.
The unknown parameters can be estimated by least squares:
$\left( \widehat{\alpha},\widehat{\gamma}\right)$ solves
\begin{align}
\min_{(\alpha', \gamma')' \in \mathcal{A} \times \Gamma }
&\mathbb{S}_{T}\left( \alpha ,\gamma \right) \equiv \frac{1}{T}
\sum_{t=1}^{T}(y_{t}-x_{t}^{\prime }\beta -x_{t}^{\prime }\delta
1\{f_{t}^{\prime }\gamma >0\})^{2} \label{original-prob}
\\
&\text{subject to:}\quad \tau_1 \leq \frac{1}{T} \sum_{t=1}^T 1\{f_{t}^{\prime }\gamma >0\} \leq \tau_2.
\label{restrction:para:original-form}
\end{align}
We assume that the restriction (\ref{restrction:para:original-form}) is satisfied when $\gamma=\gamma_0$ \emph{a.s.}
Here, $0 < \tau_1 < \tau_2 < 1$ for some predetermined $\tau_1$ and $\tau_2$ (e.g., $\tau_1 = 0.05$ and $\tau_2 = 0.95$).
In the special case that $1\{f_{t}^{\prime }\gamma _{0}>0\} = 1\{q_t > \widetilde{\gamma}_0\}$ with a scalar variable $q_t$ and a parameter $\widetilde{\gamma}_0$, it is standard to assume that the parameter space for $\widetilde{\gamma}_0$ is between the $\tau$ and $(1-\tau)$ quantiles of $q_t$ for some known $0 < \tau < 1$. We can interpret \eqref{restrction:para:original-form} as a natural generalization of this restriction so that the proportion of one regime is never too close to 0 or 1.
When
$\gamma $ is of high dimension, the naive grid search would not work well.
Dynamic programming (e.g., \citep{Bai:Perron:2003})
or
smooth global optimization (e.g., \citep{Qu:Tkachenko:17})
might be considered but are
not readily available.
We overcome this computational difficulty by replacing the naive grid search with MIO.
We present two alternative algorithms based on MIO below.
\subsection{Mixed Integer Quadratic Programming}\label{sec:joint:estimation}
Our first algorithm is based on mixed integer quadratic programming (MIQP), which jointly estimates $(\alpha,\gamma )$. It is guaranteed to obtain a global solution once it is found. To write the original least-squares problem in MIQP, we introduce
$d_t:=1\{f_t'\gamma>0\} $ and $\ell_{t} := \delta d_t$ for $ t = 1,\ldots,T$.
Then, rewrite the objective function as
\begin{equation} \label{prob2}
\frac{1}{T}\sum_{t=1}^{T}
(y_{t}-x_{t}'{\beta}- x_t'\ell_{t})^{2},
\end{equation}
which is a quadratic function of $\beta$ and $\ell_t$.
The goal is to introduce only linear constraints with respect to variables of optimization, and to construct an MIQP that is equivalent to the original least-squares problem. Then, we can apply modern MIO packages (e.g., Gurobi) to solve MIQP.
The assumption
$\alpha \in \mathcal{A}$ implies that
there exist known upper and lower bounds for $\delta_{j}$: $L_j \leq \delta_{j} \leq U_j$,
where $\delta_{j}$ denotes the $j$th element of $\delta$ for $j =1,\ldots,d_x$.
In addition, to make sure that $\ell_{j,t} = \delta_{j} d_t$ for each $j$ and $t$, we impose two additional restrictions:
\begin{align}\label{eq:bilinear:d:L}
d_t L_j \leq \ell_{j,t} \leq d_t U_j
\ \ \text{ and } \ \
L_{j} (1-d_t) \leq \delta_{j} - \ell_{j,t} \leq U_{j} (1-d_t).
\end{align}
It is then straightforward to check that these constraints imply $\ell_{j, t} =\delta_j d_t$. To introduce another key constraint, we define
$
M_t \equiv \max_{\gamma \in \Gamma} | f_t' \gamma |
$
for each $t=1,\ldots,T$, where $\Gamma$ is the parameter space for $\gamma_0$.
We can compute $M_t$ easily for each $t$ using linear programming.
We store them as inputs to our algorithm. The following new constraints along with
\eqref{restrction:para:original-form} and
\eqref{eq:bilinear:d:L} ensure that
the reformulated problem \eqref{prob2} is the same as the original problem:
$$
(d_t - 1) (M_t + \epsilon) < f_t' \gamma \leq d_t M_t,
$$
where $\epsilon > 0$ is a small predetermined constant (e.g., $\epsilon = 10^{-6}$).
The following defines an algorithm for the MIQP algorithm.
\begin{algorithm}[h]
\KwInput{$\{(y_t, x_t, f_t, M_t): t=1,\ldots,T \}$}
\KwOutput{$\left( \widehat{\alpha},\widehat{\gamma}\right)$}
Let $\bm{d} = (d_1, \ldots, d_T)'$ and $\bm{\ell} = \{\ell_{j,t}: j =1,\ldots,d_x, t = 1,\ldots,T \}$, where
$\ell_{j,t}$ is a real-valued variable. Solve the following problem:
\begin{align}\label{prob3}
\min_{\beta, \delta, \gamma, \bm{d}, \bm{\ell} }
\mathbb{Q}_{T}\left( \beta ,\boldsymbol{\ell }\right) \equiv
\frac{1}{T}\sum_{t=1}^{T}
(y_{t}-x_{t}'{\beta}- \sum_{j=1}^{d_x} x_{j,t} \ell_{j,t} )^{2}
\end{align}
subject to
\begin{align}\label{main-constraints}
\begin{split}
& (\beta, \delta) \in \mathcal{A}, \; \gamma \in \Gamma, \; d_t \in \{0, 1\}, \; L_j \leq \delta_{j} \leq U_j, \\
& (d_t - 1) (M_t + \epsilon) < f_t' \gamma \leq d_t M_t, \\
& d_t L_j \leq \ell_{j,t} \leq d_t U_j, \\
& L_{j} (1-d_t) \leq \delta_{j} - \ell_{j,t} \leq U_{j} (1-d_t), \\
& \tau_1 \leq T^{-1} \sum_{t=1}^T d_t \leq \tau_2
\end{split}
\end{align}
for each $t=1,\ldots,T$ and each $j=1,\ldots,d_x$, where $0 < \tau_1 < \tau_2 < 1$.
\caption{Mixed Integer Quadratic Programming (MIQP)}
\end{algorithm}
Our proposed algorithm is mathematically equivalent to the original least-squares problem \eqref{original-prob} subject to \eqref{restrction:para:original-form} in terms of values of objective functions. Formally, we state it as the following theorem.
\begin{thm}\label{thm:computation:joint}
Let $(\bar{\alpha},\bar{\gamma})$ denote a solution using MIQP as described above.
Then, $\mathbb{S}_{T}\left( \widehat{\alpha},\widehat{\gamma}\right) =\mathbb{S}_{T}\left(
\bar{\alpha},\bar{\gamma}\right) $, where
$(\widehat{\alpha},\widehat{\gamma})$ is defined in \eqref{original-prob}.
\end{thm}
The proposed algorithm in Section \ref{sec:joint:estimation} may run slowly when the dimension $d_x$ of $x_t$ is large. To mitigate this problem, we
reformulate MIQP in Appendix \ref{sec:joint:estimation:appendix} and use the alternative formulation in our numerical work; however, we present a simpler form here to help readers follow our basic ideas more easily.
\subsection{Block Coordinate Descent}
\begin{algorithm}[h!tb]
\KwInput{$\{(y_t, x_t, f_t, M_t): t=1,\ldots,T \}$, \texttt{MaxTime\_1}, \texttt{MaxTime\_2}}
\KwOutput{$\left( \widehat{\alpha},\widehat{\gamma}\right)$}
Set $k = 1$\;
Step 1. Obtain an initial estimate
$\left( \widehat{\alpha}^{0},\widehat{\gamma}^{0}\right)$
using MIQP with the pre-specified time limit \texttt{MaxTime\_1}\;
\If{a solution is found before reaching \texttt{MaxTime\_1},}{
set the initial estimate as the final estimate and terminate\;
}
\While{elapsed time is no greater than \texttt{MaxTime\_2}}{
Step 2. For the given $\widehat{\alpha}^{k-1},$ obtain an estimate $\widehat{\gamma}
^{k}$ via MILP:
\begin{align}\label{prob3-local}
\min_{\gamma \in \Gamma, d_1,\ldots,d_T}\frac{1}{T}\sum_{t=1}^{T}
\left\{ (x_{t}'\widehat{\delta}^{k-1} )^2
- 2(y_{t}- x_{t}'\widehat{\beta}^{k-1}) x_t'\widehat{\delta}^{k-1} \right\} d_t
\end{align}
subject to
\begin{align}\label{main-constraints-iterative}
\begin{split}
& (d_{t}-1)(M_{t}+\epsilon )<f_{t}^{\prime }\gamma \leq d_{t}M_{t}, \\
& d_{t}\in \{0,1\} \ \ \text{for each $t =1,\ldots,T$}, \\
& \tau_1 \leq \frac{1}{T} \sum_{t=1}^T d_t \leq \tau_2\text{\;}
\end{split}
\end{align}
\If{$\mathbb{S}_{T}\left( \widehat{\alpha}^{k-1} ,\widehat{\gamma}^{k} \right)
\geq \mathbb{S}_{T}\left( \widehat{\alpha}^{k-1} ,\widehat{\gamma}^{k-1} \right)$,
}
{terminate\;}
Step 3. For the given $\widehat{\gamma}^{k},$ obtain
\begin{equation*}
\widehat{\alpha}^{k}= \left[\frac{1}{T}\sum_{t=1}^T Z_{t}\left( \widehat{\gamma}^{k} \right)Z_{t}\left( \widehat{\gamma}^{k} \right)' \right]^{-1} \frac{1}{T}\sum_{t=1}^T Z_{t}\left( \widehat{\gamma}^{k} \right) y_t\text{\;}
\end{equation*}
Let $k=k+1$\;
}
\caption{Block Coordinate Descent (BCD)}\label{algo:BCD}
\end{algorithm}
While the MIQP jointly estimates $(\alpha,\gamma )$ and aims at obtaining a global solution, it might not compute as fast as necessary in large-scale problems. To mitigate the issue of scalability, we introduce a faster alternative approach based on mixed integer linear programming (MILP), whose objective function is linear in $d_t$. The algorithm solves for $\alpha$ and $\gamma$ iteratively, which we call a block coordinate descent (BCD) algorithm, starting with an initial value that can be obtained through MIQP with an early stopping rule. At step $k$, given $\widehat\alpha^{k-1}$, which is obtained in the previous step, we estimate $\gamma$ by solving
\begin{equation}\label{eq3.6}
\min_{\gamma \in \Gamma, d_{1},\ldots ,d_{T}}\frac{1}{T}\sum_{t=1}^{T}\left(
y_{t}-x_{t}^{\prime }\widehat{\beta}^{k-1}-x_{t}^{\prime }\widehat{\delta}
^{k-1}d_{t}\right) ^{2}
\end{equation}
subject to similar constraints as in MIQP. Note that the least-squares problem (\ref{eq3.6}) is linear in $d_t$ as $d_t^2=d_t$. The BCD algorithm is defined in Algorithm \ref{algo:BCD}. Intuitively speaking, it runs the MIQP algorithm for the amount of time \texttt{MaxTime\_1}, then switches to the MILP for the amount of time \texttt{MaxTime\_2}.
The BCD approach is a descent algorithm in the sense that
the least-squares objective function is a non-increasing function of $k$.
In other words, BCD in Steps 2 and 3 can provide a higher-quality solution than MIQP with an early stopping rule \texttt{MaxTime\_1}.
The time limit \texttt{MaxTime\_2} in Step 2 can be smaller than \texttt{MaxTime\_1} as
it is easier to solve an MILP problem than to solve an MIQP problem.
Furthermore, the alternative minimization approach efficiently solves for $\widehat\alpha^k$ because it has
an explicit solution.
\begin{figure}[h]
\centering
\caption{Computation Example of MIQP and BCD}\label{fig-comp_ex}
\includegraphics[width=7cm, height=7cm]{plot/time-algorithms.pdf}
\end{figure}
Figure \ref{fig-comp_ex} illustrates the performance of MIQP and BCD in one simulation draw. After spending \texttt{MaxTime\_1} (600 seconds) in Step 1, BCD switches into Step 2 and it converges to the solution quickly just in one iteration. Meanwhile, MIQP achieves a similar objective function value after spending the whole time budget of 1800 seconds.
In Monte Carlo experiments, we compare MIQP with BCD more thoroughly, subject to the same total computing time restrictions, and we demonstrate the efficiency of BCD.
\section{Asymptotic Properties with Known Factors}\label{sec:known factors}
We split the asymptotic properties of the estimator into two cases: known and unknown factors.
In this section, we consider the former.
\begin{assum}
\label{A-mixing}
\begin{enumerate}[label=(\roman*)]
\item\label{A-mixing:itm1}
$\left\{ x_{t},f_{t},\varepsilon
_{t}\right\} $ is a sequence of strictly stationary, ergodic, and $\rho $
-mixing random vectors with $\sum_{m=1}^{\infty }\rho _{m}^{1/2}<\infty $, $ \mathbb E\left\vert
x_{t}\right\vert_2 ^{4}<\infty $, and there exists a constant $C < \infty$ such that $ \mathbb E (
\left\vert x_{t}\right\vert_2 ^{8} \big| f_{t}'\gamma=0 ) < C$ and $ \mathbb E (
\varepsilon_{t}^{8} \big| f_{t}'\gamma=0 ) <C$ for all $\gamma \in \Gamma$.
\item\label{A-mixing:itm2}
$\left\{ \varepsilon _{t}\right\}$ is a martingale
difference sequence, that is, $\mathbb{E}\left( \varepsilon _{t}|\mathcal{F}_{t-1}\right) = 0$,
where $x_{t}\ $and $f_{t}$ are adapted to the filtration $\mathcal{F}_{t-1}$.
\item\label{A-mixing:itm3}
The smallest eigenvalue of $ \mathbb E [ Z_{t}\left( \gamma \right) Z_{t}\left( \gamma \right)^{\prime } ]$ is bounded away from zero for all $\gamma \in \Gamma$.
\end{enumerate}
\end{assum}
We decompose $f_{t}$ into a scalar random variable $f_{1t}$ and the other variables $f_{2t}$, so that $f_{t}^{\prime}\gamma \equiv f_{1t} + f_{2t}^{\prime} \gamma_2$.
Define $u_t := f_t'\gamma_{0}$.
\begin{assum}
\label{A:diminishing dt}
\begin{enumerate}[label=(\roman*)]
\item \label{A:dim:itm1} For some $0<\varphi <1/2\ $and $d_{0}\neq 0,\ $
$\delta _{0}=d_{0}T^{-\varphi }$.
\item \label{A:dim:itm2} $p_{u_t|f_{2t}}(u) $,
$\mathbb E [ ( x_{t}^{\prime }d_{0})^{2}|f_{2t},u_t=u ]$ and
$ \mathbb E [ ( \varepsilon_{t}x_{t}^{\prime }d_{0} ) ^{2}|f_{2t},u_t=u ]$ are continuous and bounded away from zero
at $ u=0 $ a.s.
\item \label{A:dim:itm3} For some $M<\infty $,
$
\inf_{\left\vert r\right\vert_2 =1}\mathbb{E}\left( \left\vert f_{2t}^{\prime
}r\right\vert 1\left\{ \left\vert f_{2t}\right\vert_2 \leq M\right\}
\right) >0. $
\end{enumerate}
\end{assum}
Most of the conditions in Assumptions \ref{A-mixing} and \ref{A:diminishing dt} are a natural extension of the scalar case in the literature, when $ f_t = (q_t,-1)' $ for a scalar random variable (e.g., \citep{hansen2000sample}).
Assumption \ref{A:diminishing dt}\ref{A:dim:itm3} is a rank condition on $ f_{2t} $ due to the vector of threshold parameter to be estimated and it is in terms of the first moment because of the asymptotic linear approximation of criterion function near $ \gamma_0 $. It also allows for discrete variables in $ f_{2t} $.
Assumption \ref{A:diminishing dt}\ref{A:dim:itm2}
ensures the presence of a jump, not just a kink at the change point.
\begin{thm}\label{asdist-alpha-gamma}
Let
$\mathcal{G} := \{ g \in \mathbb{R}^{d_f}: g_1 = 0 \}$.
Let Assumptions \ref{scale-normalization}, \ref{iden-assump}, \ref{A-mixing}, and \ref{A:diminishing dt} hold.
Assume further that
$\alpha_0$ is in the interior of $\mathcal{A}$ and that
$\gamma_0$ is in the interior of $\Gamma$.
In addition, let $W$ denote a mean-zero Gaussian process whose covariance kernel is given by \begin{equation}\label{H-Gaussian-process-covariance-kernel}
H\left( s,g\right) :=\frac{1}{2} \mathbb E \left[ \left(
\varepsilon _{t}x_{t}^{\prime }d_{0}\right) ^{2}\left( \left\vert
f_{t}^{\prime }g\right\vert +\left\vert f_{t}^{\prime }s\right\vert
-\left\vert f_{t}^{\prime }\left( g-s\right) \right\vert \right) p_{u_t|f_{2t}}(0) \right].
\end{equation}
Then, as $T \rightarrow \infty$, we have
\begin{align*}
\sqrt{T}(\widehat\alpha-\alpha_0) &\overset{d}{\longrightarrow } \mathcal N(0, ( \mathbb EZ_t(\gamma_0)Z_t(\gamma_0)')^{-1}\mathrm{var}(Z_t(\gamma_0)\varepsilon_t) ( \mathbb EZ_t(\gamma_0)Z_t(\gamma_0)')^{-1} ), \\
T^{1-2\varphi }\left( \widehat{\gamma}-\gamma _{0}\right) &\overset{d}{
\longrightarrow }\operatorname*{argmin}_{g \in \mathcal{G}}
\left\{
\mathbb E\left[ \left( x_{t}^{\prime }d_{0}\right) ^{2}\left\vert
f_{t}^{\prime }g\right\vert p_{u_t|f_{2t}}(0) \right] +2W\left( g\right)
\right\},
\end{align*}
where $\sqrt{T}(\widehat\alpha-\alpha_0)$ and $T^{1-2\varphi }\left( \widehat{\gamma}-\gamma _{0}\right)$ are asymptotically independent.
\end{thm}
The normalization scheme is embedded in the asymptotic distribution.
Because $\gamma _{1}=1$, the minimum in the limit is taken after
fixing the first element of $g$ at zero (recall that $\mathcal{G} = \{ g \in \mathbb{R}^{d_f}: g_1 = 0 \}$). Also note that, in the scalar threshold case, $f_{t}=\left( q_{t},-1\right) ^{\prime }$ and $\gamma_0 = (1, \widetilde{\gamma}_0)'$,
$$
H(s,g)= \frac{1}{2} \mathbb E \left[ \left(
\varepsilon _{t}x_{t}^{\prime }d_{0}\right) ^{2}\left( 2\min\left(\left|g_{2}\right|,\left|s_{2}\right|\right)1\left\{ \mathrm{sgn}\left(g_{2}\right)=\mathrm{sgn}\left(s_{2}\right)\right\} \right) p_{u_t|f_{2t}}(0) \right],
$$
which becomes the two-sided
Brownian motion, as in Hansen \citep{hansen2000sample}.
\section{Estimation with Unobserved Factors}\label{model:est:factors}
In this section, we consider the case in which the factors are estimated.
\subsection{The Model}
Consider the following factor model,
\begin{align}\label{factor-model-reg}
\mathcal Y_t =\Lambda g_{1t}+ e_t, \; t=1,\ldots, T,
\end{align}
where $\mathcal Y_t$ is an $N \times 1$ vector of time series,
$\Lambda$ is an $N\times K$ matrix of factor loadings, $g_{1t}$ is a $K \times 1$ vector of common factors, and $e_t$ is an $N \times 1$ vector of idiosyncratic components.
Throughout this section,
we make it explicit that there is a constant term in the factors, and we
replace the regression model in \eqref{model1} with
\begin{align}\label{model1-est-factor}
y_t=x_t'\beta_0+x_t'\delta_01\{g_{t}'\phi_{0} > 0 \} + \varepsilon_t,
\end{align}
where $g_t = (g_{1t}', -1)'$ is a vector of unknown factors in \eqref{factor-model-reg} plus a constant term ($-1$), and $\phi_0$ is a vector of unknown parameters. In addition, we allow $g_{1t}$ to contain lagged (dynamic) factors, but we treat them as static factors and estimate them using the PCA without losing the validity of the estimated factors.
It is well known that $g_t$ is identifiable and estimable by the PCA up to an invertible
matrix transformation (i.e., $H_T'g_t$), whose exact form will be given in Section \ref{asymp:est:factors}. Therefore,
it is customary in the literature (see, e.g., \citep{bai03, BN06}) to treat $H_T' g_t$ as a centering object in the limiting distribution of estimated factors. Following this convention,
in this section, let
\begin{align}\label{rotation-convention}
f_t:=H_T'g_t\quad \text{and} \quad \gamma_0:=H_T^{-1}\phi_0.
\end{align}
Using the fact that $ g_t'\phi_0=f_t'\gamma_0$, we can rewrite (\ref{model1-est-factor}) as the original formulation in \eqref{model1}:
$$
y_t=x_t'\beta_0+x_t'\delta_01\{f_{t}'\gamma_{0} > 0 \} + \varepsilon_t.
$$
Hence, $\gamma_0$ depends on the sample in this section but we suppress dependence on $T$ for the sake of notational simplicity.
Our estimation procedure now consists of two steps. In the first step, a $(K+1) \times 1$ vector of estimated factors and the constant term (i.e., $\widetilde f_t := (\widetilde f_{1t}', 1)'$) are obtained by the method of principal components.
To describe estimated factors, let $\mathcal Y $ be the $T \times N$ matrix whose $t$-th row is $\mathcal Y_t' $. Let $ (\widetilde f_{11}, \ldots, \widetilde f_{1T})$ be the $ K\times T$ matrix, whose rows are $K$ eigenvectors (multiplied by $\sqrt{T}$) associated with the largest $K$ eigenvalues of $ \mathcal Y \mathcal Y '/{NT}$
in decreasing order.
In the second step, unknown parameters $(\alpha_0, \gamma_0)$ are estimated by the same algorithm in Section \ref{sec:opt} with $\widetilde f_t$ as inputs.
\subsection{Regularity Conditions}\label{subsection:asymp:est:factors}
We introduce assumptions needed for asymptotic results with estimated factors.
We first replace Assumptions \ref{scale-normalization}--\ref{A:diminishing dt} with the following assumption.
Define
\begin{align}\label{def-Phi-T}
\Phi_T := \{ \phi: \phi = H_T \gamma \; \text{ for some $\gamma \in \Gamma_\epsilon$} \},
\end{align}
where $\Gamma_\epsilon$ is an $\epsilon$-enlargement of $\Gamma$. Note that $ \phi $ cannot be a vector whose first $ K $ elements are zeros due to the normalization on $ \gamma $ and the block diagonal structure of $ H_T$ that will be defined in \eqref{Bai-type-expansion}.
The space $\Phi_T$ for $\phi$ is defined through $H_T$ and excludes the case that $g_t'\phi$ is degenerate.
The $\epsilon$-enlargement of $\Gamma$ is needed because the factors are latent.
\begin{assum} \label{as9}
\begin{enumerate}[label=(\roman*)]
\item\label{as9:itm1}
Assumptions \ref{scale-normalization}, \ref{iden-assump}, and
\ref{A:diminishing dt}\ref{A:dim:itm1} hold after replacing $f_t$ and $\gamma_0$ with $g_t$ and $\phi_0$, respectively.
\item\label{as9:itm11}
$\left\{ x_{t},g_{t}, e_{t}, \varepsilon
_{t}\right\} $ is a sequence of strictly stationary, ergodic, and $\rho $ -mixing random vectors with $\sum_{m=1}^{\infty }\rho _{m}^{1/2}<\infty $,
and there exists a constant $C < \infty$ such that
$\mathbb E(\left| x_t\right|_2^8 |g_t, e_t)<C$, $\mathbb E(\varepsilon_{t}^8 |g_t, e_t)<C$ a.s., and $ g_t'\phi$ has a density that is continuous and bounded by $C$ for all $\phi \in \Phi_T$.
\end{enumerate}
\end{assum}
Recall that by the normalization in Assumption \ref{scale-normalization}, the first element of $ \gamma $ is fixed at 1. One caveat of this normalization scheme is that the sign of the first element of $f_t$ might not be the same as that of the first element of $g_t$ due to random rotation $H_T$; however, if we assume that $\delta_0 \neq 0$ and we also know the sign of one of the non-zero coefficients of $\delta_0$, then we can determine the sign of the first element of $f_t$ after estimating the model. This is a ``labeling'' problem that is common in models with hidden regimes. For simplicity, we assume that the first element of $\gamma_0$ is 1.
The following assumption is standard in the literature. In particular, we allow weak serial correlation among $e_t$.
\begin{assum}\label{assmp:factor}
\begin{enumerate}[label=(\roman*)]
\item\label{assmp:factor:itm1}
$\lim_{N\to\infty}\frac{1}{N}\Lambda'\Lambda =\Sigma_{\Lambda}$ for some $K\times K$ matrix $\Sigma_{\Lambda}$, whose eigenvalues are bounded away from both zero and infinity.
\item\label{assmp:factor:itm2}
The eigenvalues of $\Sigma_{\Lambda}^{1/2} \mathbb E (g_{1t}g_{1t}') \Sigma_{\Lambda}^{1/2}$ are distinct.
\item\label{assmp:factor:itm3}
All the eigenvalues of the $N\times N$ covariance $\mathrm{var}(e_t)$ are bounded away from both zero and infinity.
\item\label{assmp:factor:itm4} For any $ t $, $\frac{1}{N}\sum_{s=1}^{T}\sum_{i=1}^{N}|\mathbb Ee_{it}e_{is}|<C$ for some $C>0.$
\end{enumerate}
\end{assum}
Define $\lambda_i'$ to be the $i$th row of $\Lambda$, so that $\Lambda=(\lambda_1, \ldots,\lambda_N)'$.
Further, let
\begin{align*}
\xi_{s,t}&:= N^{-1/2}\sum_{i=1}^N(e_{is}e_{it}- \mathbb Ee_{is}e_{it}),\quad \psi := (TN)^{-1/2}\sum_{t=1}^T\sum_{i=1}^N g_t e_{it} \lambda_i' ,\cr
\eta_t&:= (TN)^{-1/2}\sum_{s=1}^T \sum_{i=1}^Ng_{1s} (e_{is}e_{it}- \mathbb Ee_{is}e_{it}),\quad \zeta_t:= N^{-1/2}\sum_{i=1}^N \lambda_{it} e_{it}.
\end{align*}
We require the following additional exponential-tail conditions.
\begin{assum}
There exist finite, positive constants $C, C_1$ and $c_1$ such that for any $x>0$ and
for any $\varpi\in\Xi:=\{e_{it}, g_{1t}, \xi_{s,t}, \zeta_t, vec(\psi),\eta_{t}\}$,
\[\mathbb P(|\varpi|_2>x) \leq C \exp(-C_1x^{c_1}).\]
\end{assum}
These conditions impose exponential tail conditions on various terms. First, it requires weak cross-sectional correlations among $e_{it}.$ This assumption can be verified under some low-level conditions such as the $ \alpha$-mixing condition of the type of Merlev{\`e}de et al. \citep{MPR-2011} across both $(i,t)$ and individual exponential-tailed distributions on $\{e_{it}, g_t\}$.
While the quantities in $\Xi$ are often assumed to have finite moments in the high-dimensional factor model literature, these moment bounds would no longer be sufficient in the current context. Instead, exponential-type probability bounds are more useful for us to characterize the effect of the estimated factors. To see the point, note that we have the following asymptotic expansion:
\begin{align}\label{Bai-expansion}
\widetilde f_t = \widehat f_t + r_t,\quad \widehat f_t:=H_T'(g_t+N^{-1/2}h_t).
\end{align}
Here, $r_t$ is a remainder term,
\begin{align}\label{e:6.7:h}
H_T':= \begin{pmatrix}
\widetilde H_T'&0\\
0&1
\end{pmatrix},\quad
h_t:= \begin{pmatrix}
h_{1t} \\
0
\end{pmatrix}, \quad
h_{1t}:=
(\frac{1}{N}\Lambda'\Lambda)^{-1}\frac{1}{\sqrt{N}}\Lambda' e_t,
\end{align}
and the exact form of $\widetilde H_T$ is given in \eqref{Bai-type-expansion}.
The diagonality in $H_T$ and the zero element in $h_t$ reflect the inclusion of the constant in $g_t$.
We establish the following uniform approximation result:
uniformly for $\gamma$ over a compact set,
$$
\max_{t\leq T}\left| \mathbb P(\widetilde f_t'\gamma>0)- \mathbb P(\widehat f_t'\gamma>0)\right|\leq O \left(\frac{(\log T)^c}{T} \right) + \max_{t\leq T}\mathbb P\left( |r_t|>C\frac{(\log T)^c}{T} \right)
$$
for some constants $C, c>0$.
The above exponential-tail assumption then enables us to derive a sharp bound so that
$
\max_{t\leq T}\mathbb P( | r_t |>C(\log T)^c T^{-1})
$
is asymptotically negligible.
Next, we state important technical conditions to facilitate the local asymptotic expansion of the least-squares criterion function.
A technical challenge in the analysis is that even the expected criterion function is non-smooth with respect to the factors. As such, we introduce some conditional density conditions to study the effect of estimating factors $ H_T'h_t= \sqrt{N} (\widehat f_t- f_t) $.
\begin{assum}\label{as5}
\begin{enumerate}[label=(\roman*)]
\item\label{as5:itm1}
$
\sup_{x_t, g_t} \left| \mathbb P(h_t'\phi_0<0|x_t, g_t)- ({1}/{2}) \right|=O(N^{-1/2}).
$
\item\label{as5:itm2}
Let $\sigma^2_{h, x_t, g_t}:=\operatorname*{plim}_{N\to\infty} \mathbb E[(h_t'\phi_0)^2|x_t, g_t]$ and let $\mathcal Z_t$ be a sequence of Gaussian random variables whose conditional distribution, given $x_t$ and $g_t$, is $\mathcal N(0,\sigma^2_{h, x_t, g_t})$.
Then, there are positive constants $c$, $c_0$, and $C$ such that $
\sigma^2_{h,x_t, g_t} >c_0 $ a.s., $ \sup_{x_t, g_t} \sup_{|z|<c} p_{h_t^{\prime }\phi_0| g_t,x_t}(z) < C$, and
\begin{align*}
\sup_{x_t, g_t} \sup_{|z|<c} |p_{h_t'\phi_0| g_t,x_t}( z) -p_{\mathcal Z_t|g_t,x_t}(z) | =o(1).
\end{align*}
\end{enumerate}
\end{assum}
Assumption \ref{as5} is concerned with the asymptotic behavior of the distribution of $h_t$ as $N\to\infty.$
The rate $ N^{-1/2} $ in Assumption \ref{as5}\ref{as5:itm1} is a reminiscent of the Berry--Essen theorem.
The Edgeworth expansion of the sample means at zero implies that the approximation
error is $C N^{-1/2}$, where the universal constant $C$ depends on the moments of
the summand up to the third order \citep{hall1992bootstrap}. Thus, condition \ref{as5:itm1} holds
for a broad range of setups including heteroskedastic errors $e_{it}$. For
instance, if the idiosyncratic error has the form $e_{it}=\sigma \left( g_{t}\right) \xi
_{it}$, where $g_{t}$ and $\xi _{it}$ are two independent sequences and $
\left\{ \xi _{it}\right\} $ is an independent and identically distributed (i.i.d.) sequence across $i$, then the
condition is satisfied as long as both $\sigma \left( g_{t}\right) ^{3}$ and
$\mathbb{E} \left\vert \xi _{it}\right\vert ^{3}$ are bounded. Furthermore, it holds trivially if the conditional distribution of $h_t'\phi_0$ given $x_t$ and $g_t$ is symmetric around zero or more generally if its median is zero.
Assumption \ref{as5} ensures, among other things, that
for some function $\Psi(\cdot)$ such that $\mathbb E|\Psi(x_t, g_t)|<\infty$,
$$
\mathbb E\left[\Psi(x_t, g_t)\left(1\{h_t'\phi_0\leq 0\}-1\{\mathcal Z_t\leq 0\}\right)\bigg{|}x_t,g_t\right]=O(N^{-1/2}).
$$
Above all, because $h_t$ is a cross-sectional average multiplied by $\sqrt{N}$, this assumption can be verified by a cross-sectional central limit theorem (CLT), if $\{e_{it}: i\leq N\}$ satisfies some cross-sectional mixing condition.
In the next assumption, recall that, by the identification condition, we can write $\gamma=(1, \gamma_2)$, where $1$ is the first element of $\gamma$. Correspondingly, let $f_{2t}$ and $\widehat f_{2t}$ be the subvectors of $f_t$ and $\widehat f_t$, excluding their first elements. Also, let
$ u_t:=g_t'\phi_0 = f_t'\gamma_0 $ and $ \breve{g}_t:=g_t + h_t/\sqrt{N} $.
\begin{assum}\label{as8}
There exist positive constants $c$, $c_0$, $M_0$, and $M$ such that the following hold \emph{a.s.}.
\begin{enumerate}[label=(\roman*)]
\item \label{as8:itm1} $\inf_{|u|<c}p_{\widehat f_t'\gamma_0| \widehat f_{2t}, x_t } (u) \geq c_0$ and $ \sup_{|f|_2<M_0}p_{f_{2t}|h_t}(f)<M$.
\item \label{as8:itm2} $\inf_{|u|<c}p_{ u_t| f_{2t}, h_t, x_t } (u) \geq c_0 $.
For all $|u_1| < c, |u_2| < c$,
$$ |p_{u_t|h_t'\phi_0, f_{2t},x_{t} }(u_1)- p_{u_t| h_t'\phi_0, f_{2t},x_{t} }(u_2)|\leq M|u_1-u_2|.$$
\item\label{as8:itm3}
$\inf_{|r|_2=1} \mathbb E \left[ |f_{2t}'r|^k 1\{|f_{2t}|_2< M_0\} \right] \geq c_0$ for $ k=1,2. $
\item\label{as8:itm4}
$\sup_{|r|_2=1}\sup_{|u|<c}p_{g_t'r|h_t}(u) \leq M$.
\item\label{as8:itm6}
Each of $ \inf_{\phi \in \Phi_T} |
g_t'\phi|$, $\inf_{\phi \in \Phi_T}|\breve{g}_t'\phi|$, $ \sup_{\phi \in \Phi_T} |
h_t'\phi|$, and $ \breve{g}_t'\phi_0$ has a density function bounded and continuous at zero, with
$\Phi_T$ given in \eqref{def-Phi-T}.
\item \label{as8:itm7}
$ \mathbb E [ ( x_{t}^{\prime }d_{0} ) ^{2}|g_{t}, h_t ] $ is bounded above by $ M_0 $ and below by $ c_0 $.
\item \label{as8:itm_AD} For any $ s$ and $ w $ that are linearly independent of $ \phi_{0} $, $ p_{\breve{g}_{t}'\phi_0|\breve{g}_{t}'s, \breve{g}_t'w}(u) $
and $\mathbb E( (\varepsilon_t x_t'd_0)^{2}|\breve{g}_{t}'\phi_0=u,\breve{g}_{t}'s, \breve{g}_t'w ) $
are continuously differentiable at $u=0$ with bounded
derivatives. Furthermore, $\mathbb E( (\varepsilon_t x_t'd_0)^{4}\left\vert
\breve{g}_{t}\right\vert _{2}^{2}|\breve{g}_{t}'\phi_0 ) \leq M $.
\end{enumerate}
\end{assum}
These conditions control the local characteristics of the centered least-squares criterion function near the true parameter value. As the model is perturbed by the error in the estimated factors, the centered criterion is a drifting sequence $ \widehat{f}_t $.
Its leading term changes depending on whether $N=O(T^{2-4\varphi})$ or not. The lower bounds in the above assumption are part of rank conditions that ensure that the leading terms are well defined. As a result, it entails a phase transition on the distribution of $\widehat\gamma$. Because they are rather technical, we provide a more detailed discussion on Assumption \ref{as8} in Online Appendix \ref{sec:e:discuss}.
\subsection{Rates of Convergence}\label{rates:est:factors}
The following theorem presents the rates of convergence for the estimators.
\begin{thm} \label{thm:rate}
Let Assumptions \ref{as9}--\ref{as8} hold. Suppose $T=O(N)$. Then
\begin{align*}
|\widehat\alpha-\alpha_0|_2=O_P\left(\frac{1}{\sqrt{T}}\right)
\; \text{ and } \;
|\widehat\gamma-\gamma_0|_2 =O_P\left(\frac{1}{T^{1-2\varphi}}+\frac{1}{\left( NT^{1-2\varphi }\right) ^{1/3}}\right).
\end{align*}
\end{thm}
While the convergence rate for $\widehat\alpha$ is standard,
the convergence rate of $\widehat{\gamma}$ merits further explanation.
First of all, when $N$ is relatively large so that $T^{2-4\varphi }=o\left( N\right)$, $\widehat{\gamma}-\gamma_0$ converges at a super-consistent rate of ${T^{-(1-2\varphi)}}$.
Contrary to this case, when $N=o(T^{2-4\varphi})$, the estimated threshold parameter has a cube root rate, which is similar to that of the maximum score type estimators \citep{kim1990}.
Therefore, as $\sqrt{N}/ T^{1-2\varphi}$ varies in $[0,\infty]$, the rate of convergence varies between the super-consistency rate of the usual threshold models to the cube root rate of the maximum score type estimators.
The convergence rates exhibit a continuous transition from one to the other.
To explain this transition phenomenon, we can show that uniformly in $(\alpha,\gamma)$, the objective function has the following expansion: there are functions $ R_1(\cdot)$ and $ R_2(\cdot, \cdot)$ such that
$$
{\mathbb{S}}_{T}\left( \alpha , \gamma\right)- {\mathbb{S}}_{T}\left( \alpha_0 , \gamma_0 \right)
= R_1(\gamma) + R_2(\alpha, \gamma),
$$
where $\gamma \mapsto R_1(\gamma)$ is a non-stochastic function, representing the ``mean" of the loss function, but is also highly non-smooth with respect to $\gamma$,
and $R_2(\alpha, \gamma)$ is the remaining stochastic part.
A key step is to derive a sharp lower bound for $R_1(\gamma)$.
When $N$ is relatively large, the effect of estimating latent factors is negligible, and $R_1(\gamma)$ has a high degree of non-smoothness. Similar to the usual threshold model, we have
$$
R_1(\gamma)\geq CT^{-2\varphi} |\gamma-\gamma_0|_2 - O_P(T^{-1}).
$$
This lower bound leads to a super-consistency rate. On the other hand, when $N$ is relatively small, there are extra noises arising from the cross-sectional idiosyncratic errors when estimating the latent factors, which we call ``cross-sectional noises." A remarkable feature of our model is that the cross-sectional noises help {smooth} the objective function in this case. As a result,
the behavior of $R_1(\gamma) $ is similar to that of the maximum score type estimators, where a quadratic lower bound can be derived:
$$
R_1(\gamma)\geq CT^{-2\varphi} \sqrt{N} |\gamma-\gamma_0|_2^2
-O_P(T^{-2\varphi }N^{-5/6}) .
$$
The quadratic lower bound, together with a larger error rate, then leads to a cube root rate type of convergence.
See Online Appendix \ref{sec:roadmap} for a detailed description of the roadmap of the proof.
\subsection{Consistency of Regime-Classification}
We introduce an error rate in (in-sample)
regime-classification,
\[
\widehat{R}_{T}=\frac{1}{T}\sum_{t=1}^{T}\left\vert 1\left\{ \widetilde{f}
_{t}^{\prime }\widehat{\gamma}>0\right\} -1\left\{ f_{t}^{\prime }\gamma
_{0}>0\right\} \right\vert.
\]
The uncertainty about the regime classification comes from either $\widetilde f_t$ or $\widehat{\gamma}$ or both. We establish
its convergence rate in the following theorem.
\begin{thm} \label{thm:classification}
Let Assumptions \ref{as9}--\ref{as8} hold. Suppose $T=O(N)$. Then
\[
\widehat{R}_{T}=O_P\left( \left( NT^{1-2\varphi }\right)
^{-1/3}+T^{-1+2\varphi }+N^{-1/2}\right) .
\]
\end{thm}
This is a useful corollary of the derivation of the rates of convergence for the
threshold estimator. We expect a good performance of our regime classification rule even with a moderate size of $T$.
\subsection{Asymptotic Distribution}\label{asymp:est:factors}
To describe the asymptotic distribution, we introduce additional notation.
Let $V_T$ denote the $K \times K$ diagonal matrix whose elements are the $K$ largest eigenvalues of $ \mathcal Y \mathcal Y '/{NT}$.
Define
\begin{align}
\widetilde H_T':=V_T^{-1} \frac{1}{T}\sum_{t=1}^T\widetilde f_{1t} g_{1t}' \frac{1}{N}\Lambda'\Lambda,\quad
H_T:=\mathrm{diag}(\widetilde H_T, 1),
\label{Bai-type-expansion}
\end{align}
and
$ H:=\operatorname*{plim}_{T, N\rightarrow \infty }H_{T} $, which is well defined, following Bai \citep{bai03}.
Let
\begin{equation*}
\omega:=\lim_{N,T\rightarrow \infty }\frac{\sqrt{N}}{T^{1-2\varphi }}\in \lbrack
0,\infty ],\quad
\zeta_\omega :=\max\{\omega, \omega^{1/3}\},\quad \text{and} \quad
M_\omega:=\max\{1, \omega^{-1/3}\} .
\end{equation*}
Define, for $ u_t=f_t'\gamma_0$,
\begin{align*}
A(\omega,g)&:= M_\omega \mathbb E \left[(x_td_0)^2\left(\left| f_t'g + \zeta_\omega^{-1} \mathcal Z_t \right | -
\left | \zeta_\omega^{-1} \mathcal Z_t \right | \right) \bigg{|} u_t=0\right] p_{u_t}(0) \cr
\end{align*}
for $\omega\in \left( 0,\infty\right] $, with the convention that $1/\omega=0$ for $
\omega=\infty $, and
\begin{align*}
A(0, g)&:= \mathbb{E}\left[(x_t^{\prime }d_0)^2 ( f_t'g)^2\bigg{|}u_t=0, \mathcal Z_t=0\right]p_{u_t, \mathcal Z_t}(0,0)
\end{align*}
for $\omega=0 $.
Recall $ Z_{t}(\gamma) :=(x_{t}^{\prime },x_{t}^{\prime }1\{f_{t}^{\prime
}\gamma>0\})^{\prime }$.
\begin{thm}
\label{thm:AD with Estimated f}
Let Assumptions \ref{as9}--\ref{as8} hold. Suppose $T=O(N)$.
Let $\mathcal{G}:= \{0\} \times \mathbb{R}^{K} $.
In addition, let $W$ denote the same Gaussian process as in Theorem \ref{asdist-alpha-gamma}.
Then, as $N,T\rightarrow \infty $, we have
\begin{align*}
& \sqrt{T}(\widehat{\alpha }-\alpha _{0})\overset{d}{\longrightarrow }
\mathcal{N}\left( 0,\left( \mathbb{E}Z_t(\gamma_0)Z_t(\gamma_0)^{\prime }\right) ^{-1}\mathbb{E}\left( Z_t(\gamma_0)Z_t(\gamma_0)^{\prime }\varepsilon _{t}^{2}\right) \left( \mathbb{E}Z_t(\gamma_0)Z_t(\gamma_0)^{\prime }\right) ^{-1}\right) , \\
& \left( \left( NT^{1-2\varphi }\right) ^{1/3}\wedge T^{1-2\varphi }\right)
\left( \widehat{\gamma }-\gamma _{0}\right)\overset{d}{\longrightarrow }
\operatorname*{argmin}_{g\in \mathcal{G}}A\left( \omega, g\right) +2W\left( g\right) ,
\end{align*}
and $\sqrt{T}(\widehat{\alpha }-\alpha _{0})$ and $( (
NT^{1-2\varphi }) ^{1/3}\wedge T^{1-2\varphi }) ( \widehat{\gamma }-\gamma _{0}) $ are asymptotically independent. Moreover,
$A(0, g)=\lim_{w\to0}A( w, g).$
\end{thm}
It is worth noting that
$A\left( \omega, g\right) $ is continuous
everywhere,
which implies
that the distribution of the argmin of the limit processes $A\left(
\omega,g\right) +2W\left( g\right) $ is also continuous in $\omega$ in virtue of
the argmax continuous mapping theorem [see e.g.,\citep{VW}].
Furthermore,
the asymptotic distribution of $\widehat{\gamma}$ is well defined for any $
\omega $ due to Lemma 2.6 of Kim and Pollard \citep{kim1990}. Specifically, the argmin of the limit Gaussian process is $O_P\left(
1\right) $ since $A\left( \omega, g\right) $ is a deterministic function of order
at least $\left\vert g\right\vert $ for any $\omega$ while the variance of $
W\left( g\right) $ grows at the rate of $\left\vert g\right\vert $ as $
g\rightarrow \infty $. It also possesses a unique minimizer almost
surely.
In the literature, Bai and Ng \citep{BN06, BN08} have shown that the oracle property (with regard to the estimation of the factors) holds for the linear regression if $T^{1/2}=o\left( N\right) $ and for the extremum estimation if $T^{5/8}=o\left( N\right) $, in the presence of estimated factors. Thus, it appears that the oracle property demands a larger $N$ as the nonlinearity of the estimating equation rises. In view of this, we regard our condition, $T=O(N)$, as not too stringent because we need to deal with estimated factors inside the indicator functions.
\subsection{Phase Transition}
To demonstrate that our asymptotic results are sharp, we consider a special case that
$N = T^\kappa$ for $\kappa \geq 1$.
In this case, the asymptotic results can be depicted on the $(\kappa, \varphi)$-space.
We categorize the results of Theorem \ref{thm:AD with Estimated f} into three groups.
In all three cases, the estimators enjoy certain oracle properties.
\begin{itemize}
\item Strong oracle: $T^{2-4\varphi }=o\left( N\right)
$ or $ \omega = \infty $.
This is equivalent to $\kappa> 2-4\varphi$.
The drift function $A(\infty, g)$ has a kink at $g=0$. Intuitively, a bigger $N$ makes the
estimated factors more precise. This yields the oracle result for both $\widehat{\gamma }$ and $ \widehat{\alpha} $, and the same asymptotic distribution as in the known factor case.
\item Weak oracle: $N=o\left(
T^{2-4\varphi }\right) $ or $ \omega = 0 $. This is equivalent to $\kappa< 2-4\varphi$. The drift function $A(0, g)$
is approximately quadratic in $g$ near the origin.
Because it is harder to identify the minimum when the function is
smooth than when it has a
kink at the minimum, this results in a non-oracle asymptotic distribution as well as a slower rate of convergence for $\widehat\gamma$ to $\left( NT^{1-2\varphi }\right) ^{-1/3}$. However, the asymptotic distribution for $\widehat \alpha$ are still the same as those when the unknown factors are observed. So the oracle property for $ \widehat{\alpha} $ is preserved.
\item Semi-strong oracle: $N\asymp T^{2-4\varphi }$ or $\omega\in (0,\infty)$. This is equivalent to $\kappa=2-4\varphi$. In this case, $A(\omega, g)$ has a continuous transition between the two polar cases discussed above. The effect of estimating factors is non-negligible for $ \widehat{\gamma} $ and yet the estimator enjoys the same rate of convergence. The estimator $ \widehat{\alpha} $ continues to achieve the oracle efficiency.
\end{itemize}
The phase transition occurs when $\kappa=2-4\varphi$, which is the \emph{semi-strong oracle case} and the \emph{critical boundary} of the phase transition. Changes in the convergence rates and asymptotic distributions are continuous along the critical boundary.
\begin{figure}[htb]
\caption{Phase Diagram}
\label{fig:phase}
\begin{center}
\begin{tikzpicture}[scale = 1.1]
\draw [green!20!white, fill=green!20!white] (0,0) -- (4,0) -- (0,1) -- cycle;
\draw [blue!20!white, fill=blue!20!white] (0,1) -- (0,2) -- (6,2) -- (6,0) -- (4,0);
\draw[thick, ->] (0,0) -- (6,0) node[below] {$\kappa$};
\draw[thick, ->] (0,0) -- (0,2.4) node[above] {$\varphi$};
\draw[blue,line width=3pt,densely dotted] (0,1) -- (4,0);
\draw[color=black, fill=white] (0,2) circle (.1);
\draw[color=black, fill=white] (0,1) circle (.1);
\draw[color=black, fill=white] (4,0) circle (.1);
\node[black] at (-0.3,2) {$\frac{1}{2}$};
\node[black] at (-0.3,1) {$\frac{1}{4}$};
\node[black] at (0.1,-0.1) {$1$};
\node[black] at (-0.1,0.1) {$0$};
\node[black] at (4,-0.3) {$2$};
\node[black] at (4,1.2) {Strong Oracle};
\node[black] at (1,0.3) {Weak Oracle};
\end{tikzpicture}
\end{center}
\end{figure}
Figure \ref{fig:phase} depicts a phase transition from the strong oracle phase to the weak oracle phase.
The critical boundary $\kappa=2-4\varphi$ is shown by closely dotted points in the figure.
On one hand, as $\varphi$ moves from 0 to $1/2$, the strong oracle region for $\kappa$ increases. That is, as the convergence rate for $\widehat \gamma$ becomes slower, the requirement for the minimal sample size $N$ for factor estimation becomes less stringent.
On the other hand, as $\kappa$ becomes larger, the strong oracle region for $\varphi$ increases. In other words, as $N$ becomes larger, the range of attainable oracle rates of convergence for $\widehat \gamma$ becomes wider. In this way, we provide a thorough characterization of the effect of estimated factors.
\subsection{Graphical Representation of $A\left( \omega, g\right)$}
\begin{figure}[htbp]
\caption{An Example of $A\left( \omega, g\right)$}
\label{fig-Akg}
\begin{center}
\graphicspath{ {plot/} }
\includegraphics[scale=0.25]{Fig-A-3d.pdf}
\includegraphics[scale=0.25]{Fig-A-2d1.pdf}
\includegraphics[scale=0.25]{Fig-A-2d2.pdf}
\end{center}
\end{figure}
To plot $A\left( \omega, g\right)$, we consider the simple case that
$g_{t}=\left( q_{t},-1\right) ^{\prime }$, $g=\left( 0,g_{2}\right) ^{\prime },$ $x_{t}=1$, $d_{0}=1$,
and $h_{t}$ and $q_{t}$ are independent of each other. We write $g_2=g$ for simplicity.
The left panel of Figure \ref{fig-Akg} shows the three-dimensional graph of $A\left( \omega, g\right)$,
the middle panel depicts the profile of $A\left( \omega, g\right)$ as a function of $\omega$ for several values of $g$,
and the right panel exhibits that of $A\left( \omega, g\right)$ as a function of $g$ for given values of $\omega$.
First of all, it can be seen that $A\left( \omega, g\right)$ is continuous everywhere but has a kink at $\omega=1$.
As $\omega$ approaches zero, the shape of $A\left( \omega, g\right)$ is clearly quadratic in $g$; whereas, as $\omega$ becomes larger,
it becomes almost linear in $g$. Also, $A\left( \omega, g\right)$ is quite flat around its minimum at $g=0$ when $\omega$ is close to zero; however, $A\left( \omega, g\right)$ has a sharp minimum at zero for a larger value of $\omega$. This reflects the fact that the rate of convergence increases as $\omega$ becomes larger.
\section{Inference}\label{sec:inference}
In this section, we consider inference.
Regarding $\alpha_0$,
Theorems \ref{asdist-alpha-gamma} and \ref{thm:AD with Estimated f} imply that inference for $\alpha_0$ can be carried out as if $\gamma_0$ were known. Therefore, the standard inference method based on the asymptotic normality can be carried out for $\alpha_0$ for both observed and estimated $f_t$.
We now focus on the inference issue regarding $\gamma_0$.
Let $ \theta_0 = h(\gamma_0)$ denote the parameter of interest for some known linear transformation $h(\cdot)$. For instance, this can be a particular element of $\gamma_0$ or a linear combination of the elements of $\gamma_0$. We use a quasi-likelihood ratio statistic:
\begin{align*}
LR(\theta) &:=\frac{\mathbb{
S}_{T}\left( \widehat\alpha_h ,\widehat\gamma_h \right) - {\mathbb{S}}_{T}(\widehat\alpha,\widehat\gamma)}{ {\mathbb{S}}_{T}(\widehat\alpha,\widehat\gamma)},\cr
( \widehat\alpha_h ,\widehat\gamma_h)&:=\arg\min_{\alpha ,h\left( \gamma \right) = \theta}\mathbb{
S}_{T}\left( \alpha ,\gamma \right),\quad
( \widehat\alpha ,\widehat\gamma):=\arg\min_{\alpha ,\gamma}\mathbb{
S}_{T}\left( \alpha ,\gamma \right),
\end{align*}
where $\mathbb S_T$ denotes the least-squares loss function, using $f_t$ when factors are observable, and $\widetilde f_t$ when factors are estimated.
Then, the $ 100(1-a) \% $-level confidence set for $ \theta_0 $ is $ \{\theta : LR(\theta) \leq \texttt{cv}_a \} $, where $ \texttt{cv}_a $ denotes a critical value.
As Theorem \ref{t6.1} shows, the asymptotic distribution is non-pivotal, so the critical value is computed based on the bootstrap.
\subsection{The Bootstrap with Estimated Factors}
We focus on the case of estimated factors, where we use
$\widetilde f_t$ as the ``true" factors, and denote by $f_t^*$ as the \emph{estimated factors}
in the bootstrap world.
To preserve the phase transition brought by the effect of PCA factor estimators,
$f_t^*$ should be a ``perturbed" version of $\widetilde f_t$. Specifically,
let $f_t^*$ be re-estimated factors in the bootstrap sample via PCA. This is given by Gon{\c{c}}alves and Perron \citep{gonccalves2018bootstrapping}. To maintain the cross-sectional dependence among the idiosyncratic components in the bootstrap factor models, we generate bootstrap data by
$$
\mathcal Y_t^*:= \widehat\Lambda \widetilde f_{t} + \widehat \mathrm{var}(e_t)^{1/2}\mathcal W_t^*,
$$
where $\{\mathcal W_t^*:t\leq T\}$ is a sequence of independent $N\times 1$ multivariate standard normal random vectors and $ \widehat \mathrm{var}(e_t)$ is the estimated covariance matrix of $e_t$. If the covariance is a sparse matrix,
we apply the thresholding covariance estimator of Fan, Liao, and Mincheva \citep{POET}.
Then, we apply PCA to estimate factors to obtain $\widetilde F_{t}^*$.
However,
$\widetilde F_{t}^*$ estimates $\widetilde f_t$, the ``true factors" in the bootstrap sample, up to a new rotation matrix $H_T^*$. Fortunately, such a rotation indeterminacy can be removed because $H_T^*$ is known in the bootstrap world. Following Gon{\c{c}}alves and Perron \citep{gonccalves2014bootstrapping, gonccalves2018bootstrapping}, we define
\begin{align}\label{f-star-bootstrap}
f_t^*:= H_T^{*'-1} \widetilde F_{t}^*
\end{align}
as the final ``estimated factors" in the bootstrap sample.
The bootstrap distribution of $f_t^*- \widetilde f_t$ mimics well the asymptotic sampling distribution of $\widetilde f_t -H_T'g_t$, that is $\mathcal N(0, \Sigma_h)$.
We give more details of this method, the definition of $H_T^*$,
and an alternative method based on Gaussian perturbation in Online Appendix \ref{sec:appendix:est:factors}.
\subsection{The k-Step Bootstrap Algorithm}
We now describe the bootstrap algorithm in detail.
Define
\begin{align}\label{Z-star:bootstrap}
\widetilde Z_t(\gamma):= (x_t', x_t' 1\{\widetilde f_t'\gamma>0\})'
\; \text{ and } \;
Z^*_t(\gamma):= (x_t', x_t'1\{f_t^{*'}\gamma>0\})'.
\end{align}
For each $t=1,\ldots, T$, construct $\left\{ y_{t}^{\ast }\right\}_{t\leq T} $ by
\begin{align}\label{y-star:bootstrap}
y_{t}^{\ast }:=
\widetilde Z_t \left( \widehat \gamma \right)^{\prime }\widehat{\alpha}+\eta
_{t}\widehat{\varepsilon}_{t}
\; \text{ with } \;
\widehat\varepsilon_t:=y_t-\widetilde Z_{t}\left( \widehat{\gamma}\right) ^{\prime }\widehat{\alpha},
\end{align}
where $\eta_t$ is an i.i.d.\ sequence whose mean is
zero and whose variance is one.
For example, $\eta_t \sim \mathcal{N} (0,1)$ or it can be simulated from a discrete distribution
(e.g., the Rademacher distribution).
The bootstrap least-squares loss is given by
\begin{equation}\label{eq7.1}
\mathbb S_T^*(\alpha,\gamma):=\frac{1}{T}\sum_{t=1}^T[y_t^*-Z^*_{t}\left( {\gamma}\right)'\alpha]^2.
\end{equation}
In principle, the bootstrap analog of the original constraint
is $h(\gamma)=h(\widehat\gamma)$ and the bootstrap analogous $ LR $ is defined as
$$
\widetilde{LR}^*:=\frac{\min_{\alpha ,h\left(\gamma \right) =h(\widehat\gamma)}\mathbb{
S}_{T}^*\left( \alpha ,\gamma \right) -\min_{\alpha, \gamma} {\mathbb{S}}^*_{T}( \alpha, \gamma)}{ \min_{\alpha, \gamma} {\mathbb{S}}^*_{T}( \alpha, \gamma)}.
$$
A potential computational problem for $ \widetilde{LR}^*$ is that it is necessary to fully solve two joint MIO problems: $ \min_{\alpha, \gamma} {\mathbb{S}}^*_{T}( \alpha, \gamma) $ and $ \min_{\alpha, h(\gamma)=h(\widehat\gamma)} {\mathbb{S}}^*_{T}( \alpha, \gamma)$ in each of the bootstrap repetitions. To circumvent this problem,
we adopt the approach of Andrews \citep{andrews2002higher}. Because a solution based on the original data should be close to a solution based on the bootstrapped data, within each bootstrap replication, we can employ the MILP algorithm, with
$(\widehat\alpha,\widehat\gamma)$ as the initial value, and iteratively update the algorithm for $k$ steps rather than computing the full bootstrap solutions.
A computationally convenient $k$-step LR statistic (${LR}_k^*$) and
its computational details are given in Algorithm \ref{algo:bootstrap}.
\begin{algorithm}[h!tb]
\KwInput{$\{(y_t, x_t, \widetilde f_{t}, M_t, \widehat{\varepsilon}_t): t=1,\ldots,T \}$, $\widehat \mathrm{var}(e_t)$, $\widehat\Lambda$, $ \widehat{\alpha} $, $\widehat\gamma$, $\widehat\gamma_h$, $ B $}
\KwOutput{bootstrap critical value $\texttt{cv}_a^*$}
Set $ b=1 $\;
\While{$ b \leq B $}{
Generate an i.i.d. sequence $\left\{ \eta _{t}\right\}_{t\leq T} $ whose mean is
zero and variance is one and an i.i.d. sequence of multivariate vectors $\left\{ \mathcal W^* _{t}\right\}_{t\leq T} $ from $\mathcal N(0, I)$\;
Generate $
\mathcal Y_t^*= \widehat\Lambda \widetilde f_{t} + \widehat \mathrm{var}(e_t)^{1/2}\mathcal W_t^*, t=1,...,T
$\;
Apply PCA to $ \{\mathcal Y_t^* \}$ and obtain $\widetilde F_{t}^*$ as the PCA factor estimates\;
Compute $ H_T^* $ and $ f_t^* = H_T^{*'-1} \widetilde{F}_t $, $ t=1,...,T$\;
Construct $
y_{t}^{\ast }=
\widetilde Z_t \left( \widehat \gamma \right)^{\prime }\widehat{\alpha}+\eta
_{t}\widehat{\varepsilon}_{t}$, $ t=1,...,T $, where $ \widetilde Z_t(\gamma)= (x_t', x_t' 1\{\widetilde f_t'\gamma>0\})'
$\;
Initialize at $\widehat\gamma^{* ,0}=\widehat\gamma, $ $\widehat\gamma_h^{* , 0}=\widehat\gamma_h $\;
Set $ l=1$;
\While{$ l \leq k $}{
Compute $ \widehat{\alpha}^{*, l}= \alpha^*(\widehat{\gamma}^{*, l-1} )$ and $\widehat{\alpha}_h^{*, l}=\alpha^*(\widehat{\gamma}_h^{*, l-1})$, where
\begin{eqnarray*}
\alpha^*(\gamma)&=& \left[\frac{1}{T}\sum_{t=1}^T Z^*_{t}\left(\gamma\right)Z^*_{t}\left( \gamma \right)' \right]^{-1} \frac{1}{T}\sum_{t=1}^T Z^*_{t}\left( \gamma \right) y_t^*\text{\;} \end{eqnarray*}
For the given $(\widehat{\alpha}^{*l},\widehat{\alpha}^{*l}_h)$, compute the following by MILP:
\begin{eqnarray*}
\widehat\gamma^{*,l} &=&\arg\min_{\gamma}\mathbb S_T^*(\widehat{\alpha}^{*, l}, \gamma ),\cr
\widehat\gamma^{*,l}_h &=&\arg\min_{h(\gamma)=h(\widehat\gamma)}\mathbb S_T^*(\widehat{\alpha}_h^{*, l}, \gamma )\text{\;}
\end{eqnarray*}
Let $l=l+1$\;
}
Compute
\[ {LR}_k^*:=\frac{ \mathbb{
S}_{T}^*\left( \widehat \alpha_h^{*, k} ,\widehat \gamma_h^{*, k} \right)- \mathbb{
S}_{T}^*\left( \widehat \alpha^* ,\widehat \gamma^* \right) }{ \mathbb{
S}_{T}^*\left( \widehat \alpha^* ,\widehat \gamma^* \right) }\text{\;}
\]
Let $b=b+1$\;
}
Obtain $\texttt{cv}_a^*$ by the $(1-a)$ th quantile of the empirical distribution of $LR^*_k$.
\caption{Bootstrap for Estimated Factors }\label{algo:bootstrap}
\end{algorithm}
\subsection{Asymptotic Distribution}
To describe the asymptotic distribution of the quasi-likelihood ratio statistic, let
$\sigma_\varepsilon^{2}$ be the variance of $\varepsilon_t$. In addition, recall the asymptotic distributions of $\widehat\gamma$, the minimizer of
$$
\mathbb Q(\omega, g) := A(\omega, g) +2W\left( g\right),
$$
and,
as we discussed for Theorem \ref{thm:AD with Estimated f}, $\omega=\infty$ also corresponds to the case of known factors.
Note that $A(\omega, g)$ depends on the true value $\phi_0$, the rotation matrix $H$, and the covariance matrix $\Sigma_h$. For the bootstrap sampling distribution, we consider drifting sequences around these values. For this, define
\begin{align*}
&\mathbb A(\omega, g, \Sigma, \bar H,\phi)\cr
&:=M_\omega \mathbb E \left[(x_td_0)^2\left(\left| g_t'Hg + \zeta_\omega^{-1} \mathcal W_t^{*'}\Sigma^{1/2}\bar H^{-1} \phi\right | -
\left | \zeta_\omega^{-1} \mathcal W_t^{*'}\Sigma^{1/2}\bar H^{-1} \phi \right | \right) \bigg{|} g_t'\phi=0\right] p_{g_t'\phi}(0)
\end{align*}
for $\omega\in(0,\infty]$, and
\begin{align*}
&\mathbb A( 0, g, \Sigma, \bar H,\phi)\cr
&:= \mathbb{E}\left[(x_t^{\prime }d_0)^2 ( g_t'Hg)^2\bigg{|}g_t'\phi=0, \mathcal W_t^{*'}\Sigma^{1/2}\bar H^{-1} \phi =0\right]p_{g_t'\phi, \mathcal W_t^{*'}\Sigma^{1/2}\bar H^{-1} \phi }(0,0) .
\end{align*}
Note that $A(\omega, g)= \mathbb A(\omega, g, H'\Sigma_h H, H,\phi_0)$.
\begin{assum} \label{a9.1}
(i) Uniformly for $\phi$ inside a neighborhood of $\phi_0$,
$
\sup_{x_t, f_{2t}}| p_{\breve g_t'\phi|x_t, f_{2t}}(0)-p_{g_t'\phi_1|x_t, f_{2t}}(0)|=o(1).
$
(ii) For each fixed $\omega\in[0,\infty]$ and $g$, $\mathbb A(\omega, g, S)$ is continuous with respect to $S=(\Sigma, \bar H, \phi)$.
(iii) The factor idiosyncratic component $e_t$ is independent of $(x_t, g_t)$, and
$|\widehat \mathrm{var}(e_t)- \mathrm{var}(e_t)|_2=o_P(1)$ under the matrix spectral norm.
(iv) $\inf_{\gamma}|\widehat f_t^{*'}\gamma|$ has a density (jointly with respect to $(e_t, g_t, \mathcal W_t^*)$) bounded and continuous at zero, where $\widehat f_t^{*}=\widehat f_t+ N^{-1/2} \widehat\Sigma_h^{1/2}\mathcal W_t^*$.
\end{assum}
Fan, Liao, and Mincheva \citep{POET} showed that under mild sparsity assumptions, for the matrix spectral norm, $|\widehat \mathrm{var}(e_t)- \mathrm{var}(e_t)|_2=o_P(1)$, given that $\log N$ does not grow too fast relative to $T$. The following theorem presents the asymptotic distribution of $LR$, and the validity of the $k$-step bootstrap procedure.
\begin{thm}\label{t6.1}
Suppose that assumptions of Theorem \ref{asdist-alpha-gamma} (for the known factor case) or assumptions of Theorem \ref{thm:AD with Estimated f} (for the estimated factor case)
and Assumption \ref{a9.1} hold.
Let $h(\cdot)$ be a $\mathbb R^m$-valued linear function with a fixed $m$ and let $
r_{NT}:=\left( NT^{1-2\varphi }\right) ^{1/3}\wedge T^{1-2\varphi } $, where we set $N = T^2 $ in case of the known factor.
Then, under $\mathcal H_0: h(\gamma_0)=\theta$, we have
$$
\sqrt{r_{NT}T^{1+2\varphi }} \cdot LR\to^d \sigma_{\varepsilon}^{-2}\min_{ g_h'\nabla h=0} \mathbb Q( \omega,g_h) - \sigma_{\varepsilon}^{-2}\min_{ g} \mathbb Q( \omega,g),
$$
and for any $k\geq 1$ as the number of iterations in the $k$-step bootstrap,
$$
\sqrt{r_{NT}T^{1+2\varphi }} \cdot LR_k^*\to^{d^*} \sigma_{\varepsilon}^{-2}\min_{ g_h'\nabla h=0} \mathbb Q( \omega,g_h) - \sigma_{\varepsilon}^{-2}\min_{ g} \mathbb Q( \omega,g).
$$
In the above, $\to^{d^*}$ represents the convergence in distribution with respect to the conditional distribution of $\left\{ \eta _{t},\mathcal W_t^*\right\}_{t\leq T} $ given the original data. Also, $\nabla h$ denotes the gradient of $h(\cdot)$, which is independent of $\gamma_0$ as $h$ is linear.
\end{thm}
\section{Monte Carlo Experiments}\label{sec:MC}
In this section, we study the finite sample properties of the proposed method via Monte Carlo experiments. The data are generated from the following design:
\begin{align*}
y_{t} & = x_{t}^{\prime }\beta _{0}+x_{t}^{\prime }\delta
_{0}1\{g_{t}^{\prime }\phi _{0}>0\}+\varepsilon _{t} ~~\text{for}~~t=1,\ldots,T,
\end{align*}
where $\ensuremath{\varepsilon}_t \sim N(0,0.5^2)$, $x_t \equiv (1,x_{2,t}')'$, and $g_t \equiv (g_{1,t}',-1)'$. Both $x_{2,t}$ and $g_{1,t}$ follow the vector autoregressive model of order 1:
$
x_{2,t} = \rho_x x_{2, t-1} + \nu_t,
g_{1,t} = \rho_g g_{1, t-1} + u_t,
$
where $\nu_t \sim N(0, I_{d_x-1})$ and $u_t \sim N(0,I_{K})$. When the factor $g_t$ is not observable, we instead observe $\mathcal{Y}_t$ that is generated from
$\mathcal{Y}_t = \Lambda g_{1,t} + \sqrt{K}e_t, e_t = \rho_e e_{t-1} + \omega_t$,
where $\mathcal{Y}_t$ is an $N\times 1$ vector and $\omega_t$ is an i.i.d.\ innovation generated from $N(0,I_N)$.
The terms $\ensuremath{\varepsilon}_t$, $\nu_t$, $u_t$, and $\omega_t$ are mutually independent.
In the baseline model, we set $T=N=200$, $d_x=2$, and $K=3$, and apply the MIQP algorithm.
The additional parameter values are set as follows:
$\beta_0=\delta_0=(1,1)$;
$\phi_0=(1,2/3,0,2/3)$;
$\rho_x = \mathrm{diag}(0.5,\ldots,0.5)$;
$\rho_g = \mathrm{diag}(\rho_{g,1},\ldots,\rho_{g,K})$, where $\rho_{g,k} \sim U(0.2,0.8)$ for $k=1,\ldots,K$,
the $i$th row of $\Lambda$, $\lambda_i' \sim N(0', K \cdot I_K)$; and
$\rho_e = \mathrm{diag}(\rho_{e,1},\ldots,\rho_{e,N})$, where $\rho_{e,i} \sim U(0.3,0.5)$ for $i=1,\ldots,N$.
The values of $\rho_g$ and $\rho_{e}$ are drawn only once and kept for the whole replications. The factor model design is similar to Bai and Ng \citep{bai2009boosting} and Cheng and Hansen \citep{cheng2015forecasting}.
All simulation results are based on 1,000 replications unless otherwise mentioned. We use a desktop computer equipped with an AMD RYZEN Threadripper 1950X CPU (16 cores with 3.4 GHz) and 64 GB RAM. The replication R codes for both the Monte Carlo experiments and empirical applications are available at \url{https://github.com/yshin12/fadtwo}. Also, the full simulation results can be found in Tables \ref{tb1:base}--\ref{tb:com-large-1000} in Online Appendix \ref{sec:add:sim:appendix}.
\begin{figure}[thb]
\centering
\caption{Simulation Results: Baseline Model}\label{fig-base-model}
\begin{tabular}[t]{c c}
\includegraphics[width=6cm, height=4cm]{plot/base-rmse.pdf} & \includegraphics[width=6cm, height=4cm]{plot/base-coverage.pdf}
\end{tabular}
\end{figure}
First, we study the baseline model under four scenarios: (i) when we know the correct regime, i.e.\ $\phi_0$, (Oracle); (ii) when we observe $g_t$ and know that the third factor is irrelevant (Observed Factors/No Selection); (iii) when we observe $g_t$ and have to select the relevant factors (Observed Factors/Selection); and (iv) when we do not observe $g_t$ but estimate factors from $\mathcal{Y}_t$ by PCA. We set the dimension of $\gamma$ to be 4 in (iv). Figure \ref{fig-base-model} reports the relative size of the root-mean-square errors (RMSEs) for $\beta$, $\delta$ as well as the coverage rate for the 95\% confidence intervals. As predicted by the asymptotic theory in the previous sections, the relative RMSEs over Oracle are close to 1 in all scenarios. The coverage rates for the 95\% confidence intervals are also close to the nominal value. Not surprisingly, these results on $\alpha$ are based on the good estimation performance of $\phi$ (or $\gamma$).
\begin{figure}[h]
\centering
\caption{Unobserved Factors with Different $N$}\label{fig-sim2}
\begin{tabular}[t]{c c c}
\includegraphics[scale=0.25]{plot/sim2-mean-bias.pdf} & \includegraphics[scale=0.25]{plot/sim2-rmse.pdf} & \includegraphics[scale=0.25]{plot/sim2-predict.pdf}
\end{tabular}
\par
\parbox{5in}{\footnotesize Note. The whisker plot in the panel on the right denotes one standard deviation computed over replication draws.}
\end{figure}
\begin{table}[h]
\caption{Size of Bootstrap Test}\label{tb:bootstrap}
\centering
\begin{tabular}{llcc}
\hline\hline
\multirow{2}{*}{Null hypothesis} & \multirow{2}{*}{Scenarios} & \multicolumn{2}{c}{Significance level} \\
\cline{3-4}
& & \multicolumn{1}{c}{5\%} & \multicolumn{1}{c}{1\%} \\
\hline
$H_0: \gamma_{02}=0$ & Estimated factor & 3.8\% & 0.7\% \\
$H_0: \phi_{02}=0$ & Known factor & 4.3\% & 0.5\% \\
$H_0: \phi_{02}=0 \mbox{ and } \phi_{03}=0$& Many known factors & 7.5\% & 1.1\% \\
\hline
\end{tabular}
\end{table}
\begin{figure}[h]
\centering
\caption{Computation Time over $T$, $d_x$, and $d_g$}\label{fig-time1}
\begin{tabular}[t]{c c c}
\includegraphics[width=5cm, height=4cm]{plot/time-T.pdf} &
\includegraphics[width=5cm, height=4cm]{plot/time-dx.pdf} &
\includegraphics[width=5cm, height=4cm]{plot/time-dg.pdf}
\end{tabular}
\end{figure}
\begin{figure}[h]
\centering
\caption{Large Dimensional Models ($T=500$)}\label{fig-large}
\begin{tabular}[t]{c c}
$d_x = 6$ & $d_x=10$ \\
\includegraphics[width=5cm, height=5cm]{plot/time-dg-T500-dx6.pdf} &
\includegraphics[width=5cm, height=5cm]{plot/time-dg-T500-dx10.pdf}
\\
\includegraphics[width=5cm, height=5cm]{plot/time-dg-T500-dx6-obj.pdf} &
\includegraphics[width=5cm, height=5cm]{plot/time-dg-T500-dx10-obj.pdf} \\
\end{tabular}
\par
\parbox{5.5in}{\footnotesize Note. The relative measures are calculated by dividing the outcome of BCD by that of MIQP.}
\end{figure}
Second, we focus on the unobserved factor model and investigate the performance as $N$ increases. For each simulated sample of $\{y_t, x_t, g_t\}$, we generate $\mathcal{Y}_t$ with $N=100, 200, 400, 1600$. We use the same baseline design with $T=200$, $d_x=2$, but $K=1$ to speed up computations. Figure \ref{fig-sim2} summarizes the results. The regimes are predicted more precisely as $N$ increases and the performance of the estimator improves. We observe relatively more improvements in $\gamma$ rather than $\alpha$. This is because $\widehat \alpha$ already enjoys the oracle property, provided that $T= O(N)$.
Third, we investigate the performance of the bootstrap test under three scenarios: (i) an estimated factor; (ii) a known factor; (iii) many known factors. The parameters are set as follows: $T=200$, $N=400$, $B=499$, $\ensuremath{\varepsilon}_t\sim N(0,1)$, $\eta_t \sim N(0,1)$, $\beta_0=(1,1)$, $\delta_0=(0.5,0.5)$, $\gamma_0 (\mbox{or }\phi_0) = (1,0)$ in (i) and (ii), and $\phi_0 = (1,0,0,0)$ in (iii). We test a simple null hypothesis of $H_0: \gamma_{02}(\mbox{or }\phi_{02})=0$ in (i) and (ii) and a joint hypothesis of $H_0: \phi_{02}=\phi_{03}=0$ in (iii). There is no serial correlation in the model ($\rho_x=\rho_g=\rho_e=0$). Table \ref{tb:bootstrap} reports the size of the bootstrap test in each scenario and it is satisfactory but we observe over-rejection in the joint hypothesis case.
We next investigate the computation time. We start from a set of simple models and extend to large dimensional models.
We simplify the baseline model by considering scenario (ii) (i.e., Observed/No Selection), and by setting $\rho_x=\rho_g=0$. The results are based on 100 replications. We set $T=200$, $d_x=1$, and $d_g=2$, initially and increase each dimension as follows: $T=\{200, 300, 400, 500\}$, $d_x=\{1, 2, 3, 4\}$ while keeping $T=200$ and $d_g=2$; $d_g=\{2,3,4,5\}$ while keeping $T=200$ and $d_x=1$.
Figure \ref{fig-time1} reports the computation time of MIQP.
The results indicate that the computation time stays in a reasonable bound and increases linearly as $T$ and $d_x$ increase. However, it increases exponentially as $d_g$ increases.
\begin{figure}[h]
\centering
\caption{Larger Dimensional Models using BCD ($T=1000$ and $d_x = 6$)}\label{fig-very-large}
\begin{tabular}[t]{c c}
\includegraphics[width=5cm, height=5cm]{plot/time-dg-T1000.pdf} &
\includegraphics[width=5cm, height=5cm]{plot/time-dg-T1000-obj.pdf} \\
\end{tabular}
\end{figure}
We now consider large dimensional models and handle the computational challenge by implementing the BCD algorithm in addition to MIQP. We extend the dimension of the models as $T=\{500, 1000\}$, $d_x=\{6, 8, 10\}$, and $d_g = \{6, 8, 10\}$. Note that $d_g=10$ would be quite challenging and the standard grid search method would be infeasible in practice with $T=$ 1,000. The results are based on 10 iterations of each model. We set the total time budget as 1,800 seconds for both MIQP and BCD so that each estimation terminates after that even if it does not converge. In BCD, we set \texttt{MaxTime\_1}=600 (seconds) and \texttt{MaxTime\_2}=60 (seconds).
Figure \ref{fig-large} reports the ratio of the median computation time and median objective function values between BCD and MIQP when $T=500$. BCD spends a third of the computation time, whereas MIQP spends the total time budget. BCD achieves better objective function values in all cases and the performance of MIQP deteriorates quickly as $d_g$ increases when $d_x = 10$. Figure \ref{fig-very-large} reports the summary statistics of computation time and the median objective function values of BCD when $T=$ 1,000. As the computation is more challenging, we observe that the maximum computation time is higher for all $d_g$. However, the median computation time is still around 600 seconds and the achieved objective function values are quite stable.
Based on our simulation studies, we propose to use the BCD algorithm by assigning $1/3$ of the total time budget into the maximum time (\texttt{MaxTime\_1}) for Step 1 (MIQP). When the global solution is not attainable within \texttt{MaxTime\_1}, the BCD algorithm would switch into Steps 2--3 (MILP) automatically. We recommend assigning 1/30 of the total time budget into the maximum time (\texttt{MaxTime\_2}) for each cycle of Step 2.
In summary, the simulation studies reveal that the proposed method achieves the properties predicted by the asymptotic theory, especially the oracle property of $\alpha$ and the inference based on the bootstrap method. The BCD algorithm also shows quite satisfactory results in a large dimensional change-point model whose computation is infeasible with grid search.
\section{Classifying the Regimes of US Unemployment}\label{sec:real-data-app}
We revisit the empirical application of Hansen \citep{Hansen:97}, who
considered threshold autoregressive models for the US unemployment rate.
Specifically, Hansen \citep{Hansen:97} used monthly unemployment rates (i.e., $u_t$) for males age 20 and over, and set $y_t = \Delta u_t$ in \eqref{model1}. The lag length in the autoregressive model was $p=12$ and the preferred threshold variable was
$q_{t-1} = u_{t-1} - u_{t-12}$. In this section, we investigate the usefulness of using unknown but estimated factors. We use the first principal component (i.e., $F_t$) of Ludvigson and Ng \citep{Ludvigson:Ng:09} that is estimated from 132 macroeconomic variables.
This factor not only explains the largest fraction of the total variation in their panel data set but also loads heavily on employment, production, and so on. Ludvigson and Ng call it a \emph{real factor} and thus it is a legitimate candidate for explaining the unemployment rate.
We consider three different specifications for $f_t$: (1) $f_{1t} = (q_{t-1}, -1)$, (2) $f_{2t} = (F_{t-1}, -1)$, and (3) $f_{3t} = (q_{t-1}, F_{t-1}, -1)$.
We combined the updated estimates of the real factor, which are available on Ludvigson's web page at \url{https://www.sydneyludvigson.com},
with Hansen's data, yielding a monthly sample from March 1960 to July 1996.
\begin{table}[h]
\caption{Estimation Results}
\begin{center}
\begin{tabular}{lrrrrrr}
\hline\hline
Specification & \multicolumn{2}{c}{(1) $f_{1t} = (q_{t-1}, -1)$} & \multicolumn{2}{c}{(2) $f_{2t} = (F_{t-1}, -1)$} & \multicolumn{2}{c}{(3) $f_{3t} = (q_{t-1}, F_{t-1}, -1)$} \vspace*{1ex} \\
\hline
Regime 1 (``Expansion'') & \multicolumn{2}{c}{$q_{t-1} \leq 0.302$} & \multicolumn{2}{c}{$F_{t-1} \leq -0.28$} & \multicolumn{2}{c}{$q_{t-1} + 3.55 F_{t-1} \leq -1.60$} \vspace*{1ex} \\
Prediction error & \multicolumn{2}{c}{0.0264} & \multicolumn{2}{c}{0.0272} & \multicolumn{2}{c}{0.0252} \vspace*{1ex} \\
Classification error & \multicolumn{2}{c}{0.193} & \multicolumn{2}{c}{0.106} & \multicolumn{2}{c}{0.104} \vspace*{1ex} \\
\hline
\end{tabular}
\end{center}
\label{table-hansen97-short}
\parbox{6.5in}{\footnotesize Note. See Table \ref{table-hansen97} in the Online Appendix for estimated coefficients and their heteroskedasticity-robust standard errors.
Regime 2 (``Contraction'') is the complement of regime 1.
``Prediction Error'' refers to the average of squared residuals $(T^{-1} \sum_{i=1}^T \widehat{\varepsilon}_t^2)$.
``Classification Error'' corresponds to the proportion of misclassification defined in \eqref{def:correct_matches}.
}
\end{table}
\begin{figure}[h]
\caption{Regime Classification}
\label{fig-hansen97}
\begin{center}
\graphicspath{ {plot/} }
\includegraphics[scale=0.3]{UR_nber.pdf}
\includegraphics[scale=0.3]{UR_hansen.pdf}
\includegraphics[scale=0.3]{UR_F1.pdf}
\includegraphics[scale=0.3]{UR_all.pdf}
\end{center}
\parbox{5in}{Note. The leftmost panel shows NBER recession dates in the shaded area, and
the other three panels display those with specifications (1), (2) and (3), respectively.}
\end{figure}
Table \ref{table-hansen97-short} reports estimation results that are obtained by the MIQP algorithm. We
show the goodness of fit by reporting the average squared residuals and also the results of regime misclassification relative to the NBER business cycle dates.
The latter is obtained by
\begin{align}\label{def:correct_matches}
\frac{1}{T}\sum_{t=1}^{T}\Big\vert 1\left\{ {f}_{jt}^{\prime }\widehat{\gamma}_j>0\right\} - 1_{\text{NBER}, t} \Big\vert \; \text{ for each $j=1,2,3$},
\end{align}
where $1_{\text{NBER}, t}$ is the indicator function that has value 1 if and only if the economy is in contraction according to the NBER dates.
Accordingly, we label regime 1 ``expansion'' and regime 2 ``contraction'', respectively.
Figure \ref{fig-hansen97} gives the graphical representation of regime classification.
Specification (1) suffers from the highest level of misclassification and tends to classify recessions more often than NBER.
Specification (2) mitigates the misclassification risk but at the expense of a worse goodness of fit. On one hand,
the threshold autoregressive model solely by $q_{t-1}$ fittingly explains the unemployment rate but is short of classifying the overall economic conditions satisfactorily. On the other hand, the model based only on $F_{t-1}$ is adequate at describing the underlying overall economy but does not explain the unemployment rate well.
It turns out that specification (3) has the lowest misclassification error and best explains unemployment. Thus, we have shown the real benefits of using a vector of possibly unobserved factors to explain the unemployment dynamics.
As an additional check, we tested the null hypothesis of no threshold effect. The resulting $p$-value is 0.002 based on 500 bootstrap replications, thus providing strong evidence for the existence of two regimes. See Table \ref{table-hansen97} in the Online Appendix for details and additional results.
\section{Conclusions}\label{sec:conclusions}
We have proposed a new method for estimating a two-regime regression model where
regime switching is driven by a vector of possibly unobservable factors. We show that our optimization problem can be reformulated as MIO and have presented two alternative computational algorithms.
We have also derived the asymptotic distribution of the resulting estimator under the scheme that the threshold effect shrinks to zero as the sample size tends to infinity.
As a possible interesting extension, we can consider nonparametric regime switching, where the switching indicator is replaced by $1\{F(w_t)>0\}$ with a vector of observables $w_t$ and a nonparametric function $F(\cdot)$. We intend to study this in the future.
\newpage