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.
83,918 characters · 12 sections · 82 citation commands
A sequential test procedure for the choice of the number of regimes in multivariate nonlinear models
JEL classification: C12; C32; C34; C52
Keywords: Nonlinear model; Regime Identification; Multivariate Time Series
Linear vector autoregressive models (VAR) have been a cornerstone in the analysis of multivariate time series for over four decades since the seminal paper by Sims1980. Despite their widespread use, linear models often fail to capture the complexity of real-world data, particularly when relationships exhibit nonlinear dynamics. For example, financial asset prices respond asymmetrically to unexpected macroeconomic news Anderson1999, hence requiring models that can accommodate such nonlinear behaviors.
Advances in computational power have facilitated the development of more complex models, such as the vector logistic smooth transition regression (VLSTR) and the vector threshold regression (VTR). These models offer greater flexibility by allowing for regime changes based on the value of a transition variable, see hute13 for a comprehensive review. Despite the increasing popularity of these models, their empirical application on real problems is yet limited, partly due to the challenges in model specification and the lack of robust tests for linearity and misspecification.
Therefore, proper specification tests can be crucial for these models, which are not identified if the linear or a lower-regime model is the data-generating process Davies1987. In the univariate context, a linearity test has been developed by lusate88, while Eitrheim1996 construct misspecification tests for smooth transition autoregressive (STAR) models, including an error autocorrelation test, a test of no additional nonlinearity and a test against parameter non-constancy. For multivariate models, camacho04 and yate14 have extended these tests. The former has developed a modelling strategy for a bivariate VLSTAR model along with several misspecification tests employing an equation-by-equation approach. The latter builds upon Camacho's approach by extending the linearity and misspecification tests to a system-based approach and generalizing the tests beyond two time series. This provides additional flexibility in modelling complex multivariate nonlinear dependencies and permits capturing the possible nonlinear interactions between the variables, which may be missed in an equation-by-equation approach.
In this context, determining the number of regimes is not straightforward. Several attempts to identify the number of regimes have been made in the univariate framework in Hansen1999, Gonzalo2002, and Strikholm2006. Inspired by the sequential test for structural breaks in baipe98, these last two approaches suggest choosing the number of regimes starting from a linear model, i.e. a single regime model, and testing iteratively between $m$ and $m+1$ regimes until rejection of the null hypothesis of $m$ regimes. This paper aims to fill the gap in the existing literature, proposing a easy-to-implement sequential procedure for the selection of the number of regimes in multivariate problems. The sequential procedure is a mere extension of the approach proposed in Strikholm2006 and applies both to smooth and abrupt regime-changing models. As for the univariate version in Strikholm2006, the practitioner has full control of the asymptotic significance level of the test at each step. We demonstrate that the finite-sample properties of the test procedure are satisfying either if the data are generated from a VLSTAR or a Threshold Vector Autoregressive model (TVAR).
One of the possible challenges of a system-based approach for the tests is the existence of stationarity and ergodicity conditions for the VLSTAR model. Although the papers from Saikkonen2008 and Kheifets2020 provide the conditions for stationarity and ergodicity in particular cases, the conditions for the general model are not available, therefore the test procedure proposed here works asymptotically properly only for a single-transition null hypothesis.
To validate our approach, we apply the sequential procedure to two empirical problems. On the one hand, we try to detect the number of regimes in US monthly interest rates Tsay1998. On the other hand, the sequential procedure is applied to daily Icelandic river flow data, which have shown to be nonlinear in several former applications Tong1985, Tsay1998, teya14, LivingstonJr2020. In both cases, the number of regimes detected overlaps with what was found in the related literature.
The paper is organized as follows. Section (ref) describes the vector logistic smooth transition autoregressive model. In Section (ref), we define the linearity test, while the sequential test procedure is introduced in Section (ref). The tests are then applied to simulated data in Section (ref) to compute their empirical sizes, empirical powers and selection frequencies, and to real data in Section (ref). Section (ref) concludes.
A specification for the general VLSTR model can be found in teya14. For ease of notation, in this study we do not include exogenous variables in the model, this means that we are analysing a vector logistic smooth transition autoregressive (VLSTAR) model. Let $\mathbf{y}_t$ be an $n \times 1$ vector of dependent variables, the VLSTAR model with $m$ regimes can be defined as follows:
where $\bm{\mu}_d$ is an $n \times 1$ vector of intercepts, for $d = 0, \ldots, m-1$, $\mathbf{\Phi}_{d,j}$ is an $n \times n$ matrix of parameters for the $j$-th lag and $\mathbf{G}_t^{(d)}\left(\mathbf{s}_t; \bm{\gamma}_{d}, \mathbf{c}_{d}\right)$ is a diagonal matrix of transition functions such that
where $s_{i,t}$, for $i = 1, \ldots, n$, is a weakly stationary transition variable for the $i$-th equation, while $\gamma_{i d}$ and $c_{id}$ are respectively the slope parameter and the location parameter where the transitions occur for the $d$-th regime. The lagged values of $\mathbf{y}_t$, or a combination of them camacho04, Kheifets2020, are usually chosen as $s_{i,t}$ for smooth transition models. However, a stationary exogenous variable can also be used, and, according to He2008, a temporal trend such as $s_{i,t} = t/T$ can be employed as well without violating the asymptotic theory. In this case, the VLSTAR model can be considered a special case of a time-varying autoregressive (TV-VAR) model and, for $\gamma_{i,d} \rightarrow \infty$, the changes of regimes identify structural breaks in the model. This could provide a good alternative to already existing methods for the identification of co-shifting in multivariate time series Hendry1998.
The elements of $\mathbf{G}_t^{(d)}$ in Eq. (ref) are usually specified as standard logistic functions\footnote{We use standard logistic functions because of their simplicity, but a more general version of the logistic function can also be used He2008.}
This specification is extremely flexible, since for $\gamma_{d} \rightarrow \infty$, $\forall d$, the diagonal elements of $\mathbf{G}^{(d)}_t\left(\mathbf{s}_t; \bm{\gamma}_{d}, \mathbf{c}_{d}\right)$ (for ease of notation we will refer to this function as $\mathbf{G}^{(d)}_t$) approach the indicator function, $\mathbbm{1}(s_{i,t}> c_{id})$, thus the model becomes a vector threshold autoregressive (VTAR) model as the one introduced by Tsay1998, while for $\gamma_{d} \rightarrow 0$, the model becomes a simple VAR. This means that the approach proposed in this study based on a smooth transition model can also be implemented for the selection of the number of regimes in a VTAR, for $\gamma_{d}$ sufficiently large (see Section (ref) for a discussion).
Model (ref) can be reparametrized in the following form
where $\mathbf{\Psi}_t = \left(\mathbf{I}_n, \mathbf{G}_{t}^{(1)}, \ldots, \mathbf{G}_t^{(m-1)}\right)'$ is a $m n \times n$ matrix, $\mathbf{I}_n$ is an $n \times n$ identity matrix, $\mathbf{x}_t = \left[1, \mathbf{y}_{t-1}', \mathbf{y}_{t-2}', \ldots, \mathbf{y}_{t-p}'\right]'$ is a $(1+p n) \times 1$ vector and $\mathbf{B} = \left(\mathbf{B}_1, \mathbf{B}_2, \ldots, \mathbf{B}_m\right)$ is a $(1 + p n) \times m n$ matrix of parameters, where $\mathbf{B}_d = \left(\bm{\mu}_d', \mathbf{\Phi}_{d,1}', \ldots, \mathbf{\Phi}_{d,p}'\right)'$. Setting $\mathbf{G}_t^{(0)} = \mathbf{I}_n$ indicates that no transitions are allowed before the first change of regime. The set of parameters to be estimated is $\bm{\theta} = \left\{\mathbf{B}, \bm{\Gamma}, \mathbf{C}\right\}$, where $\bm{\Gamma}$ and $\mathbf{C}$ are $n \times m$ matrices of parameters of the transition functions.
The linearity and additive nonlinearity testing problems in the model (ref) concern testing the additive ($m-1$)-th component, therefore the null hypothesis is that $\bm{\mu}_{m-1} = \mathbf{0}$, $\mathbf{\Phi}_{m-1,j} = \mathbf{0}$, $j = 1, \ldots, p$, in which case $\mathbf{G}_t^{(m-1)}$ is not identified since it contains unidentified nuisance parameters. Equivalently, the null hypothesis can be specified as $H_0 \colon \gamma_{i, m-1} = 0$, consequently $\mathbf{G}_t^{(m-1)} = (1/2)\mathbf{I}_n$, where $\mathbf{I}_n$ is an $n \times n$ identity matrix. This implies that the model is not identified because the linear component contains too many parameters that cannot be estimated consistently. The fact that a null hypothesis can be specified in different ways indicates a lack of identification of model (ref) under the null hypothesis. This problem, firstly studied by Davies1987 and Watson1985, has the direct consequence that the standard asymptotic inference does not hold as the asymptotic distribution of the test is not known under the null. To overcome it, Hansen1996 has provided an empirical null distribution by simulation and has given the asymptotic theory for inference. Nevertheless, this method is computationally demanding and applies only in the case of a common transition function among all the equations, $\mathbf{G}_t = g(s_t|\gamma,c)\mathbf{I}_n$, so the number of nuisance parameters is restricted to two. Alternatively, the use of a Taylor series approximation around the null and a Lagrange multiplier (LM) test has been used to circumvent the identification problem, see lusate88, tera94 for the univariate smooth transition model. Recently, Seong2022 consider testing both the null hypotheses in a univariate smooth transition model and combining the results in a single quasi-likelihood ratio test statistic Cho2011a, White2012. Following lusate88 and Strikholm2006, we propose to approximate the logistic function in the alternative hypothesis through a $L$-order Taylor approximation around $\gamma_i = 0$, as further discussed in Section (ref).
In our sequential procedure the linearity test is the first step, since the smooth transition model is not identified if the linear model is the true data-generating process. When the system foresees a different transition variable for each equation, linearity can be tested equation-by-equation through the test introduced by lusate88. Otherwise, the joint linearity test introduced in yate14 can be performed when a single transition variable is used. In the next sections, we deepen the theory behind the linearity and no additional nonlinearity tests already proposed in yate14.
The asymptotic normality of the score used to derive the test statistic is guaranteed under the regularity conditions provided by Basawa1976 and the Assumptions in the following Section.
By considering a 2-regime model (i.e., $m = 2$), Eq. (ref) becomes
Testing linearity in Eq. (ref) equals testing the null hypothesis $\text{H}_0: \gamma_i = 0$, $i = 1, \ldots, n$. Under the null, we have that $\mathbf{G}_t = \left(1/2\right)\mathbf{I}_n$ and that Eq. (ref) is linear, meaning that the null hypothesis creates an identification problem for the parameters in the linear combination $\mathbf{B}_1 + (1/2) \mathbf{B}_2$ and for the location parameter, $c_i$. As already pointed out above, this identification problem can be overcome by approximating the logistic function through an $L$-order Taylor approximation around $\gamma_i = 0$, such that
where $\upsilon_{i,0}, \ldots, \upsilon_{i,L}$ are the coefficients and $r_{i,t}$ is the reminder term. This means that $\mathbf{G}_t$ can be written as follows:
where $\mathbf{\Upsilon}_{l} = \text{diag}\left(\upsilon_{1,l}, \ldots, \upsilon_{n, l}\right)$ and $\mathbf{R}_t = \text{diag}\left(r_{1,t}, \ldots, r_{n,t}\right)$. Inserting Eq. (ref) in (ref) yields:
where $\mathbf{D}_0 = \mathbf{B}_1 + \mathbf{B}_2 \mathbf{\Upsilon}_0'$, $\mathbf{D}_l = \mathbf{B}_2 \mathbf{\Upsilon}_l'$ and $\bm{\varepsilon}_t^* = \mathbf{R}_t \mathbf{B}_2'\mathbf{x}_t + \bm{\varepsilon}_t$. In the auxiliary VAR in Eq. (ref), testing linearity is equal to testing the null hypothesis $\text{H}_0\colon \mathbf{D}_1 = \dots = \mathbf{D}_L = \mathbf{0}$. Under the null hypothesis $\mathbf{R}_t = \mathbf{0}$, therefore the error term is $\bm{\varepsilon}_t^* = \bm{\varepsilon}_t$, so that the distributional properties of the error process are not affected by the Taylor approximation under the null hypothesis.
Denoting $\mathbf{Y} = \left(\mathbf{y}_1, \ldots, \mathbf{y}_T\right)'$, $\mathbf{X} = \left(\mathbf{x}_1, \ldots, \mathbf{x}_T\right)'$, $\mathbf{E}^* = \left(\bm{\varepsilon}_1^*, \ldots, \bm{\varepsilon}_T^*\right)'$, $\mathbf{\tilde{D}}_L = \left(\mathbf{D}_1', \ldots, \mathbf{D}_L'\right)'$, and
Eq. (ref) can be written as
The null hypothesis is $\mathbf{\tilde{D}}_L = \mathbf{0}$, while the subscript in $\mathbf{Z}$ and $\mathbf{\tilde{D}}$ indicates the order of the Taylor expansion.
Let $\bm{\theta} = \left(\mathbf{d}_0',\mathbf{d}_1'\right)' \in \Theta$ be the unknown parameters of the model (ref) with the true values $\bm{\theta}_0$, where $\mathbf{d}_0 =\text{vec}(\mathbf{D}_0)$, $\mathbf{d}_1 = \text{vec}(\mathbf{\tilde{D}}_L)$, $\Theta = \Theta_{\mathbf{d}_0} \times \Theta_{\mathbf{d}_1}$ is the parametric space with $\Theta_0 \in \mathbbm{R}^{\tau_0}$ and $\Theta_1 \in \mathbbm{R}^{\tau_1}$, with $\tau_0 = (1+pn)n$ and $\tau_1 = (1+pn)n + 2n$. Below, we assume that $\Theta_0$ and $\Theta_1$ are compact and $\bm{\theta}_0$ is an interior point of $\Theta$. To compute a test for the null hypothesis, the log-likelihood of model (ref) for $T$ observations must be specified as follows
where
with $\mathbf{z}_t = \left(\mathbf{x}_t's_t, \mathbf{x}_t's_t^2, \ldots, \mathbf{x}_t's_t^l\right)'$ and $E\left\{\bm{\varepsilon}_t\bm{\varepsilon}_t' | \mathcal{F}_{t-1}\right\} = \mathbf{\Omega}_t$ is a positive definite covariance matrix, with $\lim_{T \rightarrow \infty}(1/T)\sum_{t=1}^{T}\mathbf{\Omega}_t = \mathbf{\Omega}$, see the following Assumption 3 for further details. Consequently, the limiting covariance matrix can be estimated from $(1/T)\sum_{t=1}^{T}\hat{\bm{\varepsilon}}_t\hat{\bm{\varepsilon}}_t'$ and can be used in the construction of the test statistic.
We need to specify the following assumptions in order to define an LM test.
Assumption 2 is a high-level assumption, while Assumption 3 guarantees the existence of the second moments for $\mathbf{y}_t$ and the convergence of the sample moments to their true values He2008 and permits the use of asymptotic theory for a martingale difference sequence (MDS), even when the assumption of i.i.d. errors is not valid, e.g., in the case of conditionally heteroskedastic errors Wang2022. Assumption 4 is a moment condition: for instance, if the model is a VLSTAR and $s_t = y_{i,t-d}$, $d >0$, this implies that $\mathbf{y}_t$ has a finite $2(L+1)$-th moment.
The block of the score vector involving the parameters under test, $\tilde{\bm{\theta}}$, can be written as follows
see for example Lutkepohl1996 and Appendix (ref). Evaluated under $H_0$, the score obtained in Eq. (ref) becomes
where $\hat{\mathbf{E}} = \left(\hat{\bm{\varepsilon}}_1, \hat{\bm{\varepsilon}}_2, \ldots, \hat{\bm{\varepsilon}}_T\right)'$, $\hat{\bm{\varepsilon}}_t = \mathbf{y}_t - \hat{\mathbf{D}}_0'\mathbf{x}_t$, and $\hat{\bm{\Omega}} = (1/T)\sum_{t=1}^{T}\hat{\bm{\varepsilon}}_t\hat{\bm{\varepsilon}}_t'$. The matrix $\hat{\mathbf{D}}_0$ is the maximum likelihood (ML) estimator of $\mathbf{D}_0$ under the null hypothesis. The consistency of the ML estimator is guaranteed under the stationarity conditions provided by Kheifets2020.
Under Assumptions 1-4, the score vector is asymptotically normally distributed with $n \cdot \text{cd}(\mathbf{Z}_L)$ degrees of freedom, where $\text{cd}(\mathbf{Z}_L)$ is the column dimension of $\mathbf{Z}_L$ Breusch1980. As the score is normal and $\mathbf{Z}_L'(\mathbf{I}_T - \mathbf{P_{\mathbf{Z}}})\mathbf{Z}_L$ is positive definite, the vectorised LM test statistic
has an asymptotic $\chi^2$-distribution with $n \cdot \text{cd}(\mathbf{Z}_L)$ degrees of freedom when the null hypothesis holds.
The statistics in (ref) can also be written as follows:
It should be noted that vectorisation and Kronecker products in Eq. (ref) are avoided in (ref). Then, we have the following result:
In an asymptotically equivalent way, the test can be performed also in the $TR^2$-form as follows
In this setting, the choice of $L$ is somewhat arbitrary. A higher order will increase the column dimension of $\mathbf{Z}_L$, but rejecting the null hypothesis would become easier, since a higher order often increases the the test. On the other hand, a lower order may lead to a test with better size properties, because it uses fewer parameters. As further discussed in lusate88 for the univariate case, choosing $L=1$ is not a good choice when $s_t = y_{t-d, i}$ for some $1 \leq d \leq p$ and for $i = 1, \ldots, n$, since the LM statistic has only trivial power against this alternative. The problem is typically solved by choosing a third-order Taylor expansion.
A special case of a VLSTAR with a common transition variable is the one with a transition function that is common to all the equations, i.e., $\mathbf{G}_t = g_t(s_t; \gamma, c) \mathbf{I}_n$ with $g(s_t; \gamma, c)$ being a scalar. In such a case, we have that
which leads to $\mathbf{\Upsilon}_l = \upsilon_l \mathbf{I}_n$ and $\mathbf{R}_t = r_t \mathbf{I}_n$. Inserting these elements in Eq. (ref), the construction of the LM-type statistic remains the same as above.
Once rejected the null of linearity, the practitioner should account for the possible presence of some nonlinearity not gathered from a 2-regime model. This means that there may exist an additional nonlinear component that enters the model additively. In this paper, we build upon the additive nonlinearity test introduced by yate14 which can also be used in a sequential procedure for the detection of the number of regimes. Following the findings by baipe98, Strikholm2006 suggest the use of a sequential testing procedure for additive nonlinearity in the univariate framework. We here extend such a procedure in the multivariate framework to both test additional nonlinearity and specify the number of regimes, $m$.
If the equations do not share the same transition variable, identifying the number of regimes is not straightforward and a suitable choice would be to select the minimum number of regimes identified in an equation-by-equation test. When a common transition variable is assumed throughout the system, the number of regimes can be identified from the following procedure which mainly extends in the multivariate framework the sequential test proposed by Strikholm2006. We further discuss in Section (ref) how this procedure can be applied also in the case of a VTAR model as the data-generating process. It should be noticed that the stationarity conditions are available only for a two-regime VLSTAR model Kheifets2020, this means that the consistency of the ML estimator, and the stability of the LM test results are guaranteed only for $H_0\colon m = 2$. Consequently, the tests can be only used to suggest the presence of at least three regimes.
The sequential testing procedure can start directly from the case of linearity testing against a 2-regime model. Hence, the first step of the procedure foresees the implementation of the linearity test shown in Eq. (ref) to test the null hypothesis of $m = 1$ against $m = 2$. If $\text{H}_0$ is rejected at a given level, $\alpha$, there could exist additive nonlinearity in the model. Therefore, the purpose of the practitioner may be sequentially testing for $m-1$ versus $m$ regimes until a non-rejection.
If we write Eq. (ref) for $m = 3$ regimes as follows
testing for non-additive nonlinearity equals to test $\text{H}_0 \colon \gamma_{2,i} = 0$, $i = 1, \ldots, n$, against the alternative $H_1 \colon \exists \gamma_{2,i} > 0$. Clearly, the test can be extended to a generic number of $m$ regimes.
As for the linearity test in Section (ref), the alternative model is not identified under the null hypothesis. Once again, Taylor's approximation of $\mathbf{G}_t^{(2)}$ allows us to overcome this problem and obtain a feasible test statistic. Using an $L$-order Taylor approximation, Eq. (ref) becomes
where $\mathbf{\Upsilon}_{l}^{(2)}$ is the diagonal matrix of coefficients of the $L$-order Taylor expansion of $g_{i,t}^{(2)}$. As for the linearity test, the null hypothesis implies $\mathbf{\Upsilon}_l^{(2)} = \mathbf{0}$ for $l = 1, \ldots, L$. By reparametrizing, Eq. (ref) can be written as
where $\mathbf{\Psi}_0 = \left(\mathbf{B}_1' + \mathbf{G}_t^{(1)}\mathbf{B}_2' + \mathbf{\Upsilon}_0^{(2)}\mathbf{B}'_{3}\right)'$, $\mathbf{\Psi}_l = \left(\mathbf{\Upsilon}_{l}^{(2)}\mathbf{B}'_{3}\right)'$ and $\bm{\varepsilon}_t^* = \mathbf{R}_t^{(2)}\mathbf{B}'_{3}\mathbf{x}_t + \bm{\varepsilon}_t$.
The null hypothesis in the VAR in Eq. (ref) is $\text{H}_0 \colon \mathbf{\Psi}_1 = \ldots = \mathbf{\Psi}_L = 0$. Let be $\mathbf{Y} = \left(\mathbf{y}_1', \ldots, \mathbf{y}_T'\right)'$, $\mathbf{X} = \left(\mathbf{x}_1', \ldots, \mathbf{x}_T'\right)'$, $\mathbf{E}$ the $T \times n$ matrix of residuals from Eq. (ref), and $\mathbf{Z}_L$ as in Eq. (ref), model (ref) can be written as
Let suppose that
and that $\mathbf{P}_{\mathbf{K}} = \mathbf{K}(\mathbf{K'K})^{-1}\mathbf{K}'$, the test statistic can be computed similarly to the one in Section (ref), therefore we can state the following result:
The asymptotic distribution of the LM statistic has the desired null distribution only when $m = 2$ in testing $\mathbf{G}_t^{(m-1)} = (1/2)\mathbf{I}_n$. There are moment conditions for the asymptotic distribution theory to be valid Eitrheim1996. In the univariate case, a STAR model with logistic-type transition functions must satisfy the condition $E(\varepsilon_t^8) < \infty$. A sufficient condition in the multivariate case is $E(\varepsilon_{i,t}^8) < \infty$, for $i = 1, \ldots, n$. To compute $K$ as in Eq. (ref), the vector of first-order partial derivatives of $\mathbf{\Psi}_t' \mathbf{B}' \mathbf{x}_t$ is necessary. This can be found in Appendix (ref). As further discussed in Section (ref), the test can also be performed using the $TR^2$-form in a multi-step regression problem.
The LM test statistics can be, hence, used in a top-down sequential testing procedure that foresees testing the null hypothesis of $\gamma_{m-1,i} = 0$ for a growing number of regimes until non-rejection. Therefore, the number of regimes to be included in the model is the minimum for which the null hypothesis of no-additive nonlinearity cannot be rejected.
The tests presented in this study are also valid if the practitioner is analysing a VTAR model. Since the VLSTAR nests the VTAR for $\gamma_{d}$ sufficiently large, the idea is to apply the tests directly on a VLSTAR approximation of the VTAR. This should also solve the drawback of finding the first derivative of the indicator function in the VTAR model or using a bootstrap procedure as proposed for the univariate framework in Giannerini2024. Let's suppose a 2-regime VTAR model with a single transition variable and a not-switching error term, defined as follows
which can alternatively be written as
It follows that the indicator function $\mathbbm{1}(\cdot)$ in Eq. (ref) can be approximated by a logistic function where the slope parameter $\gamma$ is fixed and equal to a sufficiently large value. Consequently, the estimated parameters of the approximation are consistent under the same assumptions of the VLSTAR model lusate88. This means that the aforementioned tests for a VLSTAR model can also be applied when the data-generating process is a VTAR. Moreover, in the sequential procedure for choosing the number of regimes, the estimates of the threshold parameters are super-consistent and assuming them known makes it easier to test $m = m_0-1$ against $m = m_0$, where $m_0 \geq 2$.
Consider a single lag three-regime VTAR model (where, for simplicity, we imposed $\bm{\mu}_1 = \bm{\mu}_2 = \bm{\mu}_3 = \mathbf{0}$)
which can be reformulated as
with $c_1 <c_2$. The sequential test procedure can be summarized in the following steps:
This, indeed, extends in the multivariate the sequential procedure for the definition of the number of thresholds proposed by Strikholm2006.
In vector models, the standard LM-type tests can be strongly oversized when the null hypothesis foresees the estimation of a large set of parameters and when the size of the sample is not large. In practice, the nominal size of the test tends to overestimate the true probability of type I error in finite samples, see also Honda1988. As in Laitinen1978 and Meisner1979, to overcome this limitation, we use a Bartlett-type correction that allows to rescale the degrees of freedom of the test and apply an $F$-statistic. In a Monte Carlo simulation study conducted by Bera1981, the authors show that this correction is able to correct the oversize of LM.
The Laitinen-Meisner correction consists of a degree of freedom rescaling of the form $(nT - S)/(W \times nT)$, where $n$ and $T$ are defined as before, $S$ is the number of parameters, and $W$ is the number of restrictions, see Laitinen1978 and Meisner1979. The $F$-type LM test statistic, or rescaled LM test statistic, can be computed as
and follows an $F\left(W, nT - S\right)$ distribution.
We carry out some simulation experiments for the finite-sample performance of the test procedure. Specifically, we investigate the empirical size and the power of the test, and we also report the selection frequencies for the sequential procedure. We also consider a Wilks' lambda test statistic based on Wilks' $\Lambda$-distribution anderson2003. In Appendix (ref), we show that Wilks' $\Lambda$ is applicable in our testing situation and how the test is performed in our framework.
To compute the empirical sizes of the sequential test procedure, we generate $1000$ replications from the model specified in Eq. (ref). We do not include any explanatory variable and we select a single lag for the simulation of $n = 3$ dependent variables, such that:
where $\bm{\varepsilon}_t \sim \mathcal{N}\left(0, \mathbf{I}_n\right)$, thus we are supposing uncorrelated errors.
For each realization, we estimate the VLSTAR model and compute the residuals matrix. Since the VLSTAR model is estimated numerically, relatively large samples are required for a reasonable estimation accuracy, therefore we choose $T = 400, 600, 1000$.
In a first attempt, we assess the empirical size of the test procedure by simulating $\mathbf{y}_t$ from Eq. (ref) with $m = 2$ and we test the null hypothesis of $m = 2$ against $m = 3$. The results are reported in Table (ref). The DGP relies on a parameter matrix $\mathbf{B}_1$ with 0.1 entries and diagonal values $\rho_i$, for $i = 1, \ldots, n$, sampled either from a uniform distribution $\mathcal{U}(0.3, 0.5)$ or from $\mathcal{U}(0.5,0.8)$, we also set $\mathbf{B}_2 = -\mathbf{B}_1$ and we use $c = \gamma = 2$. Finally, the common transition variable is generated by an exogenous first-order AR process, such that
where $\eta_t \sim \mathcal{N}(0, 1)$. As in Strikholm2006, we use three different nominal sizes, $\alpha = 0.10, 0.05, 0.01$. To compute the $TR^2$-form of the test statistics, we use a third-order Taylor expansion, therefore $L = 3$.
Analysing the empirical sizes of the test in Table (ref), it can be noticed that these are close to the nominal values in all of the three tests. When simulating $T = 1000$ observations, the empirical size of the LM test statistics (first three columns) and the Wilks' statistics (last three columns) almost coincides with the nominal size for $\alpha = 0.10$ and $\alpha = 0.05$, while the rescaled LM has an empirical size close to the real one for $T = 400$ and $T=600$. The persistence level of the dependent variables seems to be relevant, since a diagonal value $\rho_i \sim \mathcal{U}(0.5, 0.8)$ leads to slightly divergent results in terms of empirical size.
We further evaluate the empirical power of the tests in finite-size samples by applying it to simulated data from Eq. (ref) with $m=3$ regimes (with $\mathbf{B}_3 = -0.7\cdot \mathbf{I}_n$, $\gamma_2= 2$ and $c_2 = 4$). The power of the test is then calculated by testing the null hypothesis of $m = 2$ regimes. The results in Table (ref) suggest that the LM test has a good empirical power and that this increases with sample size. This is not entirely surprising since the VLSTAR model, especially with higher-order regimes, requires the estimation of a large set of parameters to obtain reasonable estimation accuracy.
In order to evaluate the regime choice of the procedures introduced in this paper, we also report in Table (ref) the selection frequencies when the data are simulated, as before, from model (ref) with $m = 2$, with the same characteristic presented above. The procedures start with a linearity test against a two-regime model, then foresee testing $m=2$ vs $m=3$ regimes and continue until a non-rejection.
For any sample size and persistence level, the selection frequencies indicate a high tendency to identify two regimes correctly. As expected, increasing the sample size enhances the accuracy of the tests in identifying the correct number of regimes. In fact, for both $\rho_i \sim \mathcal{U}(0.3, 0.5)$ and $\rho_i \sim \mathcal{U}(0.5, 0.8)$, the selection frequencies of $\hat{m} = 2$ are all above 97% for all the tests when the significance level is $\alpha = 0.01$ and $T=1000$.
We then compute the empirical size of the test and the selection frequencies when data are simulated from a VTAR model with $2$ regimes, specified as follows
where $s_t$ is generated from an autoregressive model as in Eq. (ref). We let $\mathbf{\Phi}_1$ vary as $\mathbf{B}_1$ before, we also set $\mathbf{\Phi}_2 = -\mathbf{\Phi}_1$, the threshold for $s_t$ is $c = 2$, while $\bm{\varepsilon}_{1,t}, \bm{\varepsilon}_{2,t} \sim \mathcal{N}(0, \mathbf{I}_n)$.
In Table (ref), we evaluate the empirical sizes of the test applied to the data simulated from a VTAR specified as in Eq. (ref). With the exception of the rescaled LM, all the tests tend to slightly over-reject the null at any significance level. Nevertheless, the difference between the empirical and the nominal size is negligible when the persistence changes.
When analysing the empirical power in Table (ref) (simulating from a three-regime model with $\bm{\Phi}_3 = -0.7\cdot \mathbf{I}_n$ and $c_2 = 4$), it can be observed that the power of the test is generally high and that improves with larger sample sizes, underlining, again, the importance of sample size in detecting additional regimes. Although the empirical power is above 75%, Wilks' $\Lambda$ tends to perform worse than the other two specifications.
Table (ref) complements these findings by showing the selection frequencies for a VTAR model when the real DGP has $m=2$. The results are consistent with those of the VLSTAR model, confirming that the sequential procedure is equally applicable and reliable for VTAR models. As with the VLSTAR model, the accuracy of the test procedure increases with sample size and significance levels.
As a robustness check, we also observe what happens when the number of dependent variables increases (we use $n=5$). We simulate only with $\rho_i \sim \mathcal{U}(0.3, 0.5)$, because the DGPs are not stationary under $\rho_i \sim \mathcal{U}(0.5, 0.8)$ Kheifets2020.
Table (ref) assesses the empirical size of the additive nonlinearity test in models with $m=2$ regimes and $n=5$ variables. The results slightly diverge from what observed with $n=3$, since the test statistics generally exhibit sizes lower to the nominal levels for the VLSTAR model and higher for the VTAR model. Nevertheless, the empirical sizes for both VLSTAR and VTAR models are overall not too far from the expected values.
The findings are more encouraging when one analyses the empirical powers in Table (ref). The power of all three tests is always higher than 0.90 and improves with larger sample sizes. For instance, with $T=1000$, the empirical power is nearly perfect, reflecting the tests' ability to correctly identify additional regimes in large samples.
As for the case of $n=3$, when the number of dependent variables is equal to 5, the procedure is capable of correctly identifying the real number of regimes for any model, sample size and level of persistence. In fact, the number of selected regimes is always equal to 2, with percentage values ranging from 76.7 to 99.8.
We first apply the sequential test procedure to the U.S. monthly interest rates data already used in Tsay1998. The dataset contains the 3-month treasury bill rates ($Y_{1,t}$) and 3-year treasury notes ($Y_{2,t}$) for the period from June 1953 to September 2022 ($T = 832$). These represent the short-term and intermediate-term series in the term structure of the interest rates. To obtain weakly stationary time series, the data have been considered as growth rates, i.e., $\mathbf{y}_t = \left(y_{1,t}, y_{2,t}\right)'$, with $y_{i,t} = \ln(Y_{i,t}) - \ln(Y_{i,t-1})$ for $i = 1,2$. The plots of the time series are shown in Fig. (ref).
As a candidate transition variable, we select the maturity spread computed as $x_t = \ln(Y_{1,t}) - \ln(Y_{2,t})$ since, according to the inverted yield curve theory Harvey1988, this would reflect the business cycle of the U.S. economy. Following Tsay1998, to avoid random fluctuations in the interest rates term structure, we use the 3-month moving average of $x_t$ as a transition variable, therefore $s_t = (x_t + x_{t-1} + x_{t-2})/3$ (see the green line in Fig. (ref)). To select the VLSTAR lag length, we use AIC and BIC criteria which suggest a single lag specification, i.e., $p = 1$. To understand if the results of our methodology are similar to the ones obtained with an alternative method, we compare them with the equation-by-equation approach proposed by camacho04.
The results of the tests are reported in Table (ref). It can be noticed in the first column that all the tests strongly reject the null hypothesis of linearity at any significance level, except for the test implemented in camacho04 which rejects at 5 and 10%. This means that, in line with what was supposed and proved by Tsay1998, a nonlinearity is present in the dynamics of the interest rates and the top-down procedure for the selection of the number of regimes presented in this paper can be applied. The procedure foresees testing the null hypothesis of $m = 2$ regimes against the alternative of $m = 3$. Once again, the null hypothesis is rejected in all the test statistics (second column of the table) introduced in this article, while the null cannot be rejected with the equation-by-equation test of camacho04. According to the results from the system-based test statistics, a 2-regime model is not enough to consider all the nonlinearity in the model. The third column of the table reports the test statistics for the null of $m = 3$ regimes. None of the tests rejects the null hypothesis. Consequently, it can be deduced that the optimal number of regimes for these time series is three, which is also what was originally supposed in Tsay1998.
As a second empirical example, the linearity and no remaining nonlinearity tests have been applied to the daily Icelandic river flow data for the period from 1972 to 1974. The time series in this dataset include river flows in cubic meters per second for two rivers, the J\"{o}kuls\'{a} and the Vatndals\'{a}, as well as the temperature and the precipitation, see Figure (ref). River flow data has been shown to be nonlinear in several former applications. For instance, Tong1985 use a univariate threshold model to estimate their relationship with temperature and precipitation, while Tsay1998, teya14, and LivingstonJr2020 apply a multivariate nonlinear model.
Following Tsay1998 and teya14, we first select the lagged temperature as the transition variable for both flow equations (Panel A of Table (ref)) and we compute both the linearity test, and the sequential procedure to identify the number of regimes in these series. According to the linearity test statistics, the null hypothesis of linearity is strongly rejected. The procedure introduced in this paper foresees to use the additive nonlinearity test on an increasing number of regimes when the null of linearity is rejected. Therefore, we perform the additive nonlinearity test when the null hypothesis is $\text{H}_0 \colon m = 2$. While the test by camacho04 does not reject the null hypothesis, the null of $m = 2$ is rejected at some significance level (i.e., $\alpha = 0.10, 0.05, 0.01$) in all the other tests, meaning that a residual nonlinearity is still present. The procedure iterates until a non-rejection is obtained at all the significance levels. We iterate again the procedure by testing the null hypothesis of $\text{H}_0 \colon m = 3$ regimes. In this case, the null hypothesis cannot be rejected for the standard and rescaled LM-type tests at any standard significance level, while Wilks' $\Lambda$ still rejects the null hypothesis if $\alpha = 0.10$. Based on these results, our top-down sequential procedure points to $m=3$ regimes as the optimal number of regimes for the Icelandic rivers' time series, which is in line with what found in teya14.
We also apply the tests using the lagged precipitation as a transition variable (Panel B of Table (ref)). As for the case of the lagged temperature as a transition variable, the null hypothesis of linearity is rejected for all the tests. Once rejected the hypothesis of a linear model, we conduct the top-down procedure for the selection of the number of regimes. In this case, the additive nonlinearity tests are not able to reject the null hypothesis of $m = 2$ regimes. This means that, with the lagged precipitation as a transition variable, the optimal model is a 2-regime model.
In this paper, we developed a simple method for selecting the number of regimes in multivariate nonlinear models with no restrictions on the number of dependent variables and transitions, also giving a more formal context for the linearity and no additional nonlinearity tests introduced in yate14.
The results on small-sample properties of the tests are of interest because they highlight that the empirical sizes are affected by the dimension of the model, the size of the sample and the persistence of the time series. We find that the standard LM tests tend to be size-distorted when the time series are almost non-stationary. We also show that Wilks' $\Lambda$ statistic has satisfying size properties, and is recommended for empirical use. Nevertheless, the size of the LM test can be adjusted using a proper bootstrapped version, although this has not been addressed in this work. Not surprisingly, the power experiments demonstrate that the joint test is more powerful in finite samples when the number of temporal observations is large. The selection frequencies of the number of regimes reflect what is already observed through empirical sizes and powers. The sequential procedure is capable of correctly identifying the number of regimes for any sample size, and for both smooth and abrupt regime changes. Finally, our simulation study underlines that the sequential procedure introduced in this paper can be applied either to detect the number of regimes in smoothly changing time series and abrupt regime-changing time series.
When we apply the sequential test procedure to real data, we can observe that the tests introduced in this paper lead to more rejections with respect to the test introduced by camacho04. On the one hand, our approach foresees that the VLSTAR model is correctly specified, and that deviations from linearity are due to remaining nonlinearity. In contrast, the test proposed by camacho04 may be less sensitive to misspecification, since each equation is tested separately and may have different sources of nonlinearity. On the other hand, more rejections may indicate that the tests on the whole system have greater power to detect nonlinearity, since they account for the joint behaviour of all the equations in the system. As a result, the system-based test may be more likely to identify nonlinear relationships that are present across multiple equations in the system.
A possible implementation for future research could foresee the use of the procedure to detect the number of structural breaks in a multivariate linear model. In fact, if the transition variable is a temporal trend, the VLSTAR model becomes a time-varying parameter model and the changes in regimes coincide with smooth structural breaks.
I offer my sincerest gratitude to Timo T\"{e}rasvirta who has contributed to the early version of this paper, and whose comments have helped me improve the theoretical background behind the sequential procedure. I retain the responsibility for any errors and shortcomings in this work.