Extracted main text — title through conclusion, appendix excluded. This is what our citation measures are computed over, published so the extraction can be checked by eye.
86,685 characters · 22 sections · 37 citation commands
Factor-Driven Two-Regime Regression
\doparttoc \faketableofcontents
\thispagestyle{empty}
\onehalfspacing
\setcounter{page}{1} \pagenumbering{arabic}
Suppose that $y_t$ is generated from
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 (ref). 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 (ref) and (ref) a factor-driven two-regime regression model.
Our paper is closely related to the literature on threshold models with unknown change points (see, e.g., chan1993consistency, hansen2000sample, ling_1999, Seijo:Sen:11a, Seo-Linton, and 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 chan1993consistency and Hansen hansen2000sample studied the model in which $1\{f_t'\gamma_0>0\}$ in ((ref)) 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 (ref), 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.\ 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., 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 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., hansen2000sample), and the latter by a compound poison process (e.g., 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 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
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
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 super-consistency rate (e.g., hansen2000sample). On the other hand, when $\omega=0$, it turns out that $A\left( \omega, g\right)$ is quadratic in $g$, corresponding to a cube root rate similar to the maximum score estimator (e.g., 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 weak-oracle limiting distribution to a semi-strong oracle one, and then to a 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), we propose the least-squares estimator and algorithms to compute the proposed estimator. In Section (ref), we establish asymptotic theory when $f_t$ is directly observed. In Section (ref), 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), we develop bootstrap inference, and in Section (ref) we give the results of Monte Carlo experiments. In Section (ref), we illustrate our methods by applying them to threshold autoregressive models of unemployment. We conclude in Section (ref). 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 a.s. means almost surely.
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
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.
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 (ref), 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
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.
Online Appendix (ref) provides sufficient conditions for Assumption (ref).
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
We assume that the restriction ((ref)) is satisfied when $\gamma=\gamma_0$ 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 (ref) 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., Bai:Perron:2003) or smooth global optimization (e.g., 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.
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
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:
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 (ref) and (ref) ensure that the reformulated problem (ref) 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.
Our proposed algorithm is mathematically equivalent to the original least-squares problem (ref) subject to (ref) in terms of values of objective functions. Formally, we state it as the following theorem.
The proposed algorithm in Section (ref) may run slowly when the dimension $d_x$ of $x_t$ is large. To mitigate this problem, we reformulate MIQP in Appendix (ref) 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.
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
subject to similar constraints as in MIQP. Note that the least-squares problem ((ref)) is linear in $d_t$ as $d_t^2=d_t$. The BCD algorithm is defined in Algorithm (ref). Intuitively speaking, it runs the MIQP algorithm for the amount of time MaxTime_1, then switches to the MILP for the amount of time 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 MaxTime_1. The time limit MaxTime_2 in Step 2 can be smaller than 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.
Figure (ref) illustrates the performance of MIQP and BCD in one simulation draw. After spending 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.
We split the asymptotic properties of the estimator into two cases: known and unknown factors. In this section, we consider the former.
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}$.
Most of the conditions in Assumptions (ref) and (ref) are a natural extension of the scalar case in the literature, when $ f_t = (q_t,-1)' $ for a scalar random variable (e.g., hansen2000sample). Assumption (ref)(ref) 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)(ref) ensures the presence of a jump, not just a kink at the change point.
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 hansen2000sample.
In this section, we consider the case in which the factors are estimated.
Consider the following factor model,
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 (ref) with
where $g_t = (g_{1t}', -1)'$ is a vector of unknown factors in (ref) 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). Therefore, it is customary in the literature (see, e.g., 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
Using the fact that $ g_t'\phi_0=f_t'\gamma_0$, we can rewrite ((ref)) as the original formulation in (ref): $$ 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) with $\widetilde f_t$ as inputs.
We introduce assumptions needed for asymptotic results with estimated factors. We first replace Assumptions (ref)--(ref) with the following assumption. Define
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 (ref). 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.
Recall that by the normalization in Assumption (ref), 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$.
Define $\lambda_i'$ to be the $i$th row of $\Lambda$, so that $\Lambda=(\lambda_1, \ldots,\lambda_N)'$. Further, let
We require the following additional exponential-tail conditions.
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. 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:
Here, $r_t$ is a remainder term,
and the exact form of $\widetilde H_T$ is given in (ref). 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) $.
Assumption (ref) is concerned with the asymptotic behavior of the distribution of $h_t$ as $N\to\infty.$ The rate $ N^{-1/2} $ in Assumption (ref)(ref) 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 hall1992bootstrap. Thus, condition (ref) 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) 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} $.
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) in Online Appendix (ref).
The following theorem presents the rates of convergence for the estimators.
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 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) for a detailed description of the roadmap of the proof.
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.
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$.
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
and $ H:=\operatorname*{plim}_{T, N\rightarrow \infty }H_{T} $, which is well defined, following Bai bai03. Let
Define, for $ u_t=f_t'\gamma_0$,
for $\omega\in \left( 0,\infty\right] $, with the convention that $1/\omega=0$ for $ \omega=\infty $, and
for $\omega=0 $. Recall $ Z_{t}(\gamma) :=(x_{t}^{\prime },x_{t}^{\prime }1\{f_{t}^{\prime }\gamma>0\})^{\prime }$.
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.,VW]. Furthermore, the asymptotic distribution of $\widehat{\gamma}$ is well defined for any $ \omega $ due to Lemma 2.6 of Kim and Pollard 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 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.
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) into three groups. In all three cases, the estimators enjoy certain oracle properties.
The phase transition occurs when $\kappa=2-4\varphi$, which is the semi-strong oracle case and the critical boundary of the phase transition. Changes in the convergence rates and asymptotic distributions are continuous along the critical boundary.
Figure (ref) 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.
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) 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.
In this section, we consider inference. Regarding $\alpha_0$, Theorems (ref) and (ref) 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:
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) shows, the asymptotic distribution is non-pivotal, so the critical value is computed based on the bootstrap.
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 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 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 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 gonccalves2014bootstrapping, gonccalves2018bootstrapping, we define
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).
We now describe the bootstrap algorithm in detail. Define
For each $t=1,\ldots, T$, construct $\left\{ y_{t}^{\ast }\right\}_{t\leq T} $ by
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
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 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).
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), $\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
for $\omega\in(0,\infty]$, and
Note that $A(\omega, g)= \mathbb A(\omega, g, H'\Sigma_h H, H,\phi_0)$.
Fan, Liao, and Mincheva 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.
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:
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 bai2009boosting and Cheng and Hansen 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)--(ref) in Online Appendix (ref).
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) 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$).
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) 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) 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) 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.
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 MaxTime_1=600 (seconds) and MaxTime_2=60 (seconds). Figure (ref) 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) 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 (MaxTime_1) for Step 1 (MIQP). When the global solution is not attainable within 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 (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.
We revisit the empirical application of Hansen Hansen:97, who considered threshold autoregressive models for the US unemployment rate. Specifically, Hansen Hansen:97 used monthly unemployment rates (i.e., $u_t$) for males age 20 and over, and set $y_t = \Delta u_t$ in (ref). 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 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 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.
Table (ref) 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
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) 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) in the Online Appendix for details and additional results.
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.