# A sequential test procedure for the choice of the number of regimes in
multivariate nonlinear models

Andrea Bucci  Department of Economics and Law, University of Macerata, Italy

###### Abstract

This paper proposes a sequential test procedure for determining the number of
regimes in nonlinear multivariate autoregressive models. The procedure relies
on linearity and no additional nonlinearity tests for both multivariate smooth
transition and threshold autoregressive models. We conduct a simulation study
to evaluate the finite-sample properties of the proposed test in small
samples. Our findings indicate that the test exhibits satisfactory size
properties, with the rescaled version of the Lagrange Multiplier test
statistics demonstrating the best performance in most simulation settings. The
sequential procedure is also applied to two empirical cases, the US monthly
interest rates and Icelandic river flows. In both cases, the detected number
of regimes aligns well with the existing literature.

JEL classification: C12; C32; C34; C52

Keywords: Nonlinear model; Regime Identification; Multivariate Time Series

##  1 Introduction

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 Sims, (1980). 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 (Anderson et al.,, 1999),
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 Hubrich and Teräsvirta, (2013) 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 (Davies,, 1987). In the univariate context, a linearity
test has been developed by Luukkonen et al., (1988), while Eitrheim and
Teräsvirta, (1996) 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, Camacho, (2004) and Teräsvirta and Yang, 2014a  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 Hansen, (1999), Gonzalo and Pitarakis, (2002), and
Strikholm and Teräsvirta, (2006). Inspired by the sequential test for
structural breaks in Bai and Perron, (1998), 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
Strikholm and Teräsvirta, (2006) and applies both to smooth and abrupt regime-
changing models. As for the univariate version in Strikholm and Teräsvirta,
(2006), 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 Saikkonen, (2008) and Kheifets and Saikkonen, (2020)
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 (Tsay,, 1998). 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 (Tong et al.,, 1985; Tsay,, 1998;
Teräsvirta and Yang, 2014b, ; Livingston Jr and Nur,, 2020). In both cases,
the number of regimes detected overlaps with what was found in the related
literature.

The paper is organized as follows. Section 2 describes the vector logistic
smooth transition autoregressive model. In Section 3, we define the linearity
test, while the sequential test procedure is introduced in Section 4. The
tests are then applied to simulated data in Section 5 to compute their
empirical sizes, empirical powers and selection frequencies, and to real data
in Section 6. Section 7 concludes.

##  2 The VLSTAR model

A specification for the general VLSTR model can be found in Teräsvirta and
Yang, 2014b . 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:

$\displaystyle\mathbf{y}_{t}$
$\displaystyle=\bm{\mu}_{0}+\sum_{j=1}^{p}\mathbf{\Phi}_{0,j}\mathbf{y}_{t-j}+\mathbf{G}_{t}^{(1)}\left(\mathbf{s}_{t};\bm{\gamma}_{1},\mathbf{c}_{1}\right)\left[\bm{\mu}_{1}+\sum_{j=1}^{p}\mathbf{\Phi}_{1,j}\mathbf{y}_{t-j}\right]+$
$\displaystyle\dots+\mathbf{G}_{t}^{(m-1)}\left(\mathbf{s}_{t};\bm{\gamma}_{m-1},\mathbf{c}_{m-1}\right)\left[\bm{\mu}_{m-1}+\sum_{j=1}^{p}\mathbf{\Phi}_{m-1,j}\mathbf{y}_{t-j}\right]+\bm{\varepsilon}_{t}$
(1)

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

$\mathbf{G}_{t}^{(d)}\left(\mathbf{s}_{t};\bm{\gamma}_{d},\mathbf{c}_{d}\right)=\text{diag}\left\\{g_{1,t}^{(d)}\left(s_{1,t};\gamma_{1d},c_{1d}\right),\ldots,g_{n,t}^{(d)}\left(s_{n,t};\gamma_{nd},c_{nd}\right)\right\\}$
(2)

where $s_{i,t}$, for $i=1,\ldots,n$, is a weakly stationary transition
variable for the $i$-th equation, while $\gamma_{id}$ 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 (Camacho,, 2004; Kheifets and
Saikkonen,, 2020), are usually chosen as $s_{i,t}$ for smooth transition
models. However, a stationary exogenous variable can also be used, and,
according to He et al., (2008), 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 (Hendry and Mizon,, 1998).

The elements of $\mathbf{G}_{t}^{(d)}$ in Eq. (2) are usually specified as
standard logistic functions111We use standard logistic functions because of
their simplicity, but a more general version of the logistic function can also
be used (He et al.,, 2008).

$g_{i,t}^{(d)}\left(s_{i,t};\gamma_{id},c_{id}\right)=\left[1+\exp\left\\{-\gamma_{id}\left(s_{i,t}-c_{id}\right)\right\\}\right]^{-1},\quad\gamma_{id}>0.$

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 Tsay, (1998), 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 4.1 for a discussion).

Model (2) can be reparametrized in the following form

$\mathbf{y}_{t}=\left\\{\sum_{d=1}^{m}\mathbf{G}_{t}^{(d-1)}\mathbf{B}_{r}^{\prime}\right\\}\mathbf{x}_{t}+\bm{\varepsilon}_{t}=\left[\mathbf{I}_{n}\
\mathbf{G}_{t}^{(1)}\ldots\
\mathbf{G}_{t}^{(m-1)}\right]\begin{bmatrix}\mathbf{B}_{1}\\\
\mathbf{B}_{2}\\\ \vdots\\\
\mathbf{B}_{m}\end{bmatrix}\mathbf{x}_{t}+\bm{\varepsilon}_{t}=\mathbf{\Psi}_{t}^{\prime}\mathbf{B}^{\prime}\mathbf{x}_{t}+\bm{\varepsilon}_{t}$
(3)

where
$\mathbf{\Psi}_{t}=\left(\mathbf{I}_{n},\mathbf{G}_{t}^{(1)},\ldots,\mathbf{G}_{t}^{(m-1)}\right)^{\prime}$
is a $mn\times n$ matrix, $\mathbf{I}_{n}$ is an $n\times n$ identity matrix,
$\mathbf{x}_{t}=\left[1,\mathbf{y}_{t-1}^{\prime},\mathbf{y}_{t-2}^{\prime},\ldots,\mathbf{y}_{t-p}^{\prime}\right]^{\prime}$
is a $(1+pn)\times 1$ vector and
$\mathbf{B}=\left(\mathbf{B}_{1},\mathbf{B}_{2},\ldots,\mathbf{B}_{m}\right)$
is a $(1+pn)\times mn$ matrix of parameters, where
$\mathbf{B}_{d}=\left(\bm{\mu}_{d}^{\prime},\mathbf{\Phi}_{d,1}^{\prime},\ldots,\mathbf{\Phi}_{d,p}^{\prime}\right)^{\prime}$.
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 (2)
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 (2) under the null
hypothesis. This problem, firstly studied by Davies, (1987) and Watson and
Engle, (1985), 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, Hansen, (1996) 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 Luukkonen et al., (1988),
Teräsvirta, (1994) for the univariate smooth transition model. Recently, Seong
et al., (2022) consider testing both the null hypotheses in a univariate
smooth transition model and combining the results in a single quasi-likelihood
ratio test statistic (Cho et al.,, 2011; White and Cho,, 2012). Following
Luukkonen et al., (1988) and Strikholm and Teräsvirta, (2006), 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 3.

##  3 Linearity test

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 Luukkonen et al., (1988). Otherwise, the joint
linearity test introduced in Teräsvirta and Yang, 2014a  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 Teräsvirta and Yang, 2014a .

###  3.1 Testing linearity with a common transition variable

The asymptotic normality of the score used to derive the test statistic is
guaranteed under the regularity conditions provided by Basawa et al., (1976)
and the Assumptions in the following Section.

By considering a 2-regime model (i.e., $m=2$), Eq. (3) becomes

$\mathbf{y}_{t}=\mathbf{B}_{1}^{\prime}\mathbf{x}_{t}+\mathbf{G}_{t}\mathbf{B}_{2}^{\prime}\mathbf{x}_{t}+\bm{\varepsilon}_{t}.$
(4)

Testing linearity in Eq. (4) 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. (4) 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

$g_{i,t}\left(s_{t}\mid\gamma_{i},c_{i}\right)\approx\sum_{l=0}^{L}\upsilon_{i,l}s_{t}^{l}+r_{i,t}$

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:

$\displaystyle\mathbf{G}_{t}$
$\displaystyle\approx\text{diag}\left\\{\sum_{l=0}^{L}\upsilon_{1,l}s_{t}^{l}+r_{1,t},\ldots,\sum_{l=0}^{L}\upsilon_{n,l}s_{t}^{l}+r_{n,t}\right\\}$
$\displaystyle\approx\sum_{l=0}^{L}\mathbf{\Upsilon}_{l}s_{t}^{l}+\mathbf{R}_{t}$
(5)

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. (3.1) in (4) yields:

$\displaystyle\mathbf{y}_{t}$
$\displaystyle=\mathbf{B}_{1}^{\prime}\mathbf{x}_{t}+\left(\sum_{l=0}^{L}\mathbf{\Upsilon}_{l}s_{t}^{l}+\mathbf{R}_{t}\right)\mathbf{B}_{2}^{\prime}\mathbf{x}_{t}+\bm{\varepsilon}_{t}$
$\displaystyle=\left(\mathbf{B}_{1}^{\prime}+\mathbf{\Upsilon}_{0}\mathbf{B}_{2}^{\prime}\right)\mathbf{x}_{t}+\sum_{l=1}^{L}\mathbf{\Upsilon}_{l}\mathbf{B}_{2}^{\prime}\mathbf{x}_{t}s_{t}^{l}+\mathbf{R}_{t}\mathbf{B}_{2}^{\prime}\mathbf{x}_{t}+\bm{\varepsilon}_{t}$
$\displaystyle=\mathbf{D}_{0}^{\prime}\mathbf{x}_{t}+\sum_{l=1}^{L}\mathbf{D}_{l}^{\prime}\mathbf{x}_{t}s_{t}^{l}+\bm{\varepsilon}_{t}^{*}$
(6)

where
$\mathbf{D}_{0}=\mathbf{B}_{1}+\mathbf{B}_{2}\mathbf{\Upsilon}_{0}^{\prime}$,
$\mathbf{D}_{l}=\mathbf{B}_{2}\mathbf{\Upsilon}_{l}^{\prime}$ and
$\bm{\varepsilon}_{t}^{*}=\mathbf{R}_{t}\mathbf{B}_{2}^{\prime}\mathbf{x}_{t}+\bm{\varepsilon}_{t}$.
In the auxiliary VAR in Eq. (3.1), 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)^{\prime}$,
$\mathbf{X}=\left(\mathbf{x}_{1},\ldots,\mathbf{x}_{T}\right)^{\prime}$,
$\mathbf{E}^{*}=\left(\bm{\varepsilon}_{1}^{*},\ldots,\bm{\varepsilon}_{T}^{*}\right)^{\prime}$,
$\mathbf{\tilde{D}}_{L}=\left(\mathbf{D}_{1}^{\prime},\ldots,\mathbf{D}_{L}^{\prime}\right)^{\prime}$,
and

$\mathbf{Z}_{L}=\begin{bmatrix}\mathbf{x}_{1}^{\prime}s_{1}&\mathbf{x}_{1}^{\prime}s_{1}^{2}&\ldots&\mathbf{x}_{1}^{\prime}s_{1}^{L}\\\
\mathbf{x}_{2}^{\prime}s_{2}&\mathbf{x}_{2}^{\prime}s_{2}^{2}&\ldots&\mathbf{x}_{2}^{\prime}s_{2}^{L}\\\
\vdots&\vdots&\ddots&\vdots\\\
\mathbf{x}_{T}^{\prime}s_{T}&\mathbf{x}_{T}^{\prime}s_{T}^{2}&\ldots&\mathbf{x}_{T}^{\prime}s_{T}^{L}\end{bmatrix},$
(7)

Eq. (3.1) can be written as

$\mathbf{Y}=\mathbf{X}\mathbf{D}_{0}+\mathbf{Z}_{L}\mathbf{\tilde{D}}_{L}+\mathbf{E}^{*}.$
(8)

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}^{\prime},\mathbf{d}_{1}^{\prime}\right)^{\prime}\in\Theta$
be the unknown parameters of the model (8) 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 (8) for $T$
observations must be specified as follows

$\ell_{T}(\bm{\theta})=\sum_{t=1}^{T}\ell_{t}(\bm{\theta})=k-\frac{1}{2}\sum_{t=1}^{T}\log|\mathbf{\Omega}_{t}|-\frac{1}{2}\sum_{t=1}^{T}\bm{\varepsilon}_{t}^{\prime}\mathbf{\Omega_{t}}^{-1}\bm{\varepsilon}_{t}$
(9)

where

$\bm{\varepsilon}_{t}=\mathbf{y}_{t}-\mathbf{D}_{0}^{\prime}\mathbf{x}_{t}-\sum_{l=1}^{L}\mathbf{D}_{l}^{\prime}\mathbf{x}_{t}s_{t}^{l}=\mathbf{y}_{t}-\mathbf{D}_{0}^{\prime}\mathbf{x}_{t}-\mathbf{\tilde{D}}_{L}^{\prime}\mathbf{z}_{t},$

with
$\mathbf{z}_{t}=\left(\mathbf{x}_{t}^{\prime}s_{t},\mathbf{x}_{t}^{\prime}s_{t}^{2},\ldots,\mathbf{x}_{t}^{\prime}s_{t}^{l}\right)^{\prime}$
and
$E\left\\{\bm{\varepsilon}_{t}\bm{\varepsilon}_{t}^{\prime}|\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}^{\prime}$
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 1.

The log-likelihood $\ell_{T}(\bm{\theta})$, defined as in Eq. (9), is twice
continuously differentiable with respect to $\bm{\theta}$ in an open
neighbourhood of $\mathbf{D}_{1}=\mathbf{0}$.

######  Assumption 2.

The maximum likelihood estimators of the parameters $\mathbf{D}_{0}$ are
consistent under the null hypothesis $\mathbf{D}_{1}=0$.

######  Assumption 3.

The stochastic sequence $\left\\{\bm{\varepsilon}_{t}\right\\}$ is a
martingale difference sequence with respect to an increasing sequence of
$\sigma$-fields, $\mathcal{F}_{t}$ with

$\underset{t}{\sup}\
E\left\\{|\varepsilon_{i,t}|^{2+\alpha}|\mathcal{F}_{t-1}\right\\}<\infty\qquad\text{a.s.}$

for some $\alpha>0$ and $i=1,\ldots,n$, with
$E\left\\{\bm{\varepsilon}_{t}\bm{\varepsilon}_{t}^{\prime}|\mathcal{F}_{t-1}\right\\}=\mathbf{\Omega}_{t}$,
where $\mathbf{\Omega}_{t}$ is a positive definite matrix with the following
asymptotic limit

$\lim_{T\rightarrow\infty}(1/T)\sum_{t=1}^{T}\mathbf{\Omega}_{t}=\mathbf{\Omega}\qquad\text{a.s.}$

for some positive definite matrix $\mathbf{\Omega}$.

######  Assumption 4.

$\mathbf{X^{\prime}X}$ and
$\mathbf{Z}_{L}^{\prime}(\mathbf{I}-\mathbf{P}_{\mathbf{X}})\mathbf{Z}_{L}$,
where $\mathbf{P}_{X}$ is the limiting projection matrix of $\mathbf{X}$,
$\mathbf{P}_{\mathbf{X}}=\mathbf{X}(\mathbf{X^{\prime}X})^{-1}\mathbf{X}^{\prime}$,
are positive definite matrices.

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 (He et al.,, 2008) 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 (Wang et al.,, 2022). 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

$\displaystyle\frac{\partial\ell_{T}(\tilde{\bm{\theta}})}{\partial\mathbf{\tilde{D}}_{L}}$
$\displaystyle=-\frac{\partial}{\partial\mathbf{\tilde{D}}_{L}}(1/2)\sum_{t=1}^{T}\bm{\varepsilon}_{t}^{\prime}\bm{\Omega}^{-1}\bm{\varepsilon}_{t}=\frac{\partial}{\partial\mathbf{\tilde{D}}_{L}}\sum_{t=1}^{T}\mathbf{z}_{t}^{\prime}\mathbf{\tilde{D}}_{L}\bm{\Omega}^{-1}\bm{\varepsilon}_{t}$
$\displaystyle=\sum_{t=1}^{T}\mathbf{z}_{t}\bm{\varepsilon}_{t}^{\prime}\bm{\Omega}^{-1}=\mathbf{Z}_{L}^{\prime}\mathbf{E}\bm{\Omega}^{-1}$
(10)

see for example Lütkepohl, (1996) and Appendix A. Evaluated under $H_{0}$, the
score obtained in Eq. (10) becomes

$\frac{\partial\ell_{T}(\tilde{\bm{\theta}})}{\partial\mathbf{\tilde{D}}_{L}}\mid_{H_{0}}=\sum_{t=1}^{T}\mathbf{z}_{t}\hat{\bm{\varepsilon}}_{t}^{\prime}\hat{\bm{\Omega}}^{-1}=\mathbf{Z}_{L}^{\prime}\hat{\mathbf{E}}\hat{\bm{\Omega}}^{-1}$

where
$\hat{\mathbf{E}}=\left(\hat{\bm{\varepsilon}}_{1},\hat{\bm{\varepsilon}}_{2},\ldots,\hat{\bm{\varepsilon}}_{T}\right)^{\prime}$,
$\hat{\bm{\varepsilon}}_{t}=\mathbf{y}_{t}-\hat{\mathbf{D}}_{0}^{\prime}\mathbf{x}_{t}$,
and
$\hat{\bm{\Omega}}=(1/T)\sum_{t=1}^{T}\hat{\bm{\varepsilon}}_{t}\hat{\bm{\varepsilon}}_{t}^{\prime}$.
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 Kheifets
and Saikkonen, (2020).

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}$
(Breusch and Pagan,, 1980). As the score is normal and
$\mathbf{Z}_{L}^{\prime}(\mathbf{I}_{T}-\mathbf{P_{\mathbf{Z}}})\mathbf{Z}_{L}$
is positive definite, the vectorised LM test statistic

$LM_{L}=\text{vec}\left(\hat{\mathbf{E}}^{\prime}\mathbf{Z}_{L}\right)^{\prime}\left\\{\left(\mathbf{Z}_{L}^{\prime}(\mathbf{I}_{T}-\mathbf{P}_{\mathbf{X}})\mathbf{Z}_{L}\right)\otimes\hat{\bm{\Omega}}\right\\}^{-1}\text{vec}\left(\hat{\mathbf{E}}^{\prime}\mathbf{Z}_{L}\right)$
(11)

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 (11) can also be written as follows:

$\displaystyle LM_{L}$
$\displaystyle=\text{vec}\left(\hat{\mathbf{E}}^{\prime}\mathbf{Z}_{L}\right)^{\prime}\left\\{\left(\mathbf{Z}_{L}^{\prime}(\mathbf{I}_{T}-\mathbf{P}_{\mathbf{X}})\mathbf{Z}_{L}\right)\otimes\hat{\bm{\Omega}}\right\\}^{-1}\text{vec}\left(\hat{\mathbf{E}}^{\prime}\mathbf{Z}_{L}\right)$
$\displaystyle=\text{vec}\left(\hat{\mathbf{E}}^{\prime}\mathbf{Z}_{L}\right)^{\prime}\left\\{\left(\mathbf{Z}_{L}^{\prime}(\mathbf{I}_{T}-\mathbf{P}_{\mathbf{X}})\mathbf{Z}_{L}\right)^{-1}\otimes\hat{\bm{\Omega}}^{-1}\right\\}\text{vec}\left(\hat{\mathbf{E}}^{\prime}\mathbf{Z}_{L}\right)$
$\displaystyle=\text{vec}\left(\hat{\mathbf{E}}^{\prime}\mathbf{Z}_{L}\right)^{\prime}\text{vec}\left\\{\hat{\bm{\Omega}}^{-1}\hat{\mathbf{E}}^{\prime}\mathbf{Z}_{L}\left(\mathbf{Z}_{L}^{\prime}(\mathbf{I}_{T}-\mathbf{P}_{\mathbf{X}})\mathbf{Z}_{L}\right)^{-1}\right\\}$
$\displaystyle=\text{tr}\left\\{\mathbf{Z}_{L}^{\prime}\hat{\mathbf{E}}\hat{\bm{\Omega}}^{-1}\hat{\mathbf{E}}^{\prime}\mathbf{Z}_{L}\left(\mathbf{Z}_{L}^{\prime}(\mathbf{I}_{T}-\mathbf{P}_{\mathbf{X}})\mathbf{Z}_{L}\right)^{-1}\right\\}$
$\displaystyle=\text{tr}\left\\{\hat{\bm{\Omega}}^{-1}\hat{\mathbf{E}}^{\prime}\mathbf{Z}_{L}\left[\mathbf{Z}_{L}^{\prime}(\mathbf{I}_{T}-\mathbf{P}_{\mathbf{X}})\mathbf{Z}_{L}\right]^{-1}\mathbf{Z}_{L}^{\prime}\hat{\mathbf{E}}\right\\}.$
(12)

It should be noted that vectorisation and Kronecker products in Eq. (11) are
avoided in (12). Then, we have the following result:

######  Theorem 1.

The LM test statistic for the null hypothesis,
$\text{H}_{0}\colon\gamma_{i}=0$, $i=1,\ldots,n$ in Eq. (4), or
$\text{H}_{0}\colon\mathbf{D}_{L}=\mathbf{0}$ in Eq. (8), can be computed as
follows:

$LM_{L}=\text{tr}\left\\{\mathbf{\hat{\Omega}}^{-1}\left(\mathbf{Y}-\mathbf{X}\mathbf{\hat{D}}_{0}\right)^{\prime}\mathbf{Z}_{L}\left[\mathbf{Z}_{L}^{\prime}\left(\mathbf{I}_{T}-\mathbf{P}_{X}\right)\mathbf{Z}_{L}\right]^{-1}\mathbf{Z}_{L}^{\prime}\left(\mathbf{Y}-\mathbf{X}\mathbf{\hat{D}}_{0}\right)\right\\}$
(13)

where $\mathbf{\hat{D}_{0}}$ is the estimate of $\mathbf{D}_{0}$. Under the
null hypothesis the test statistic has a $\chi^{2}$-distribution with
$Ln\left(1+np\right)$ degrees of freedom.

###### Proof.

See Appendix B. ∎

In an asymptotically equivalent way, the test can be performed also in the
$TR^{2}$-form as follows

  1. 1.

Estimate the restricted model under the null hypothesis. Collect the residuals
$\mathbf{\hat{\bm{\varepsilon}}}_{t}=\mathbf{y}_{t}-\mathbf{X}\hat{\mathbf{D}}_{0}$.
Compute the matrix residual sum of squares
$\mathbf{\hat{E}}^{\prime}\mathbf{\hat{E}}$, where
$\mathbf{\hat{E}}=\left[\mathbf{\hat{\bm{\varepsilon}}}_{1},\ldots,\mathbf{\hat{\bm{\varepsilon}}}_{T}\right]^{\prime}$.

  2. 2.

Regress $\mathbf{\hat{E}}$ on $\mathbf{X}$ and $\mathbf{Z}_{L}$. Collect the
residuals, $\mathbf{\hat{\Xi}}$, and form the matrix residual sum of squares
$\mathbf{\hat{\Xi}}^{\prime}\mathbf{\hat{\Xi}}$.

  3. 3.

Compute the test statistic

$\displaystyle LM_{TR^{2}}$
$\displaystyle=T\cdot\text{tr}\left\\{\left(\mathbf{\hat{E}}^{\prime}\mathbf{\hat{E}}\right)^{-1}\left(\mathbf{\hat{E}}^{\prime}\mathbf{\hat{E}}-\mathbf{\hat{\Xi}}^{\prime}\mathbf{\hat{\Xi}}\right)\right\\}$
$\displaystyle=T\left(n-\text{tr}\left\\{\left(\mathbf{\hat{E}}^{\prime}\mathbf{\hat{E}}\right)^{-1}\mathbf{\hat{\Xi}}^{\prime}\mathbf{\hat{\Xi}}\right\\}\right).$
(14)

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
Luukkonen et al., (1988) 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

$g_{t}(s_{t};\gamma,c)=\sum_{l=0}^{L}\upsilon_{l}s_{t}^{l}+r_{t}$

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. (3.1),
the construction of the LM-type statistic remains the same as above.

##  4 Determining the number of regimes

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 Teräsvirta and Yang, 2014a  which can also be used in a
sequential procedure for the detection of the number of regimes. Following the
findings by Bai and Perron, (1998), Strikholm and Teräsvirta, (2006) 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
Strikholm and Teräsvirta, (2006). We further discuss in Section 4.1 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 (Kheifets and Saikkonen,, 2020),
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. (13) 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. (2) for $m=3$ regimes as follows

$\mathbf{y}_{t}=\mathbf{B}_{1}^{\prime}\mathbf{x}_{t}+\mathbf{G}_{t}^{(1)}\mathbf{B}_{2}^{\prime}\mathbf{x}_{t}+\mathbf{G}_{t}^{(2)}\mathbf{B}_{3}^{\prime}\mathbf{x}_{t}+\bm{\varepsilon}_{t},$
(15)

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 3, 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. (15)
becomes

$\mathbf{y}_{t}=\mathbf{B}_{1}^{\prime}\mathbf{x}_{t}+\mathbf{G}_{t}^{(1)}\mathbf{B}_{2}^{\prime}\mathbf{x}_{t}+\left(\sum_{l=0}^{L}\mathbf{\Upsilon}_{l}^{(2)}s_{t}^{l}+\mathbf{R}_{t}^{(2)}\right)\mathbf{B}^{\prime}_{3}\mathbf{x}_{t}+\bm{\varepsilon}_{t}$
(16)

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. (16) can be written as

$\displaystyle\mathbf{y}_{t}$
$\displaystyle=\mathbf{B}_{1}^{\prime}\mathbf{x}_{t}+\mathbf{G}_{t}^{(1)}\mathbf{B}_{2}^{\prime}\mathbf{x}_{t}+\mathbf{G}_{t}^{(2)}\mathbf{B}^{\prime}_{3}\mathbf{x}_{t}+\varepsilon_{t}$
$\displaystyle=\mathbf{B}_{1}^{\prime}\mathbf{x}_{t}+\mathbf{G}_{t}^{(1)}\mathbf{B}_{2}^{\prime}\mathbf{x}_{t}+\left(\sum_{l=0}^{L}\mathbf{\Upsilon}_{l}^{(2)}s_{t}^{l}+\mathbf{R}_{t}^{(2)}\right)\mathbf{B}^{\prime}_{3}\mathbf{x}_{t}+\bm{\varepsilon}_{t}$
$\displaystyle=\left(\mathbf{B}_{1}^{\prime}+\mathbf{G}_{t}^{(1)}\mathbf{B}_{2}^{\prime}+\mathbf{\Upsilon}_{0}^{(2)}\mathbf{B}^{\prime}_{3}\right)\mathbf{x}_{t}+\sum_{l=1}^{L}\mathbf{\Upsilon}_{l}^{(2)}\mathbf{B}^{\prime}_{3}\mathbf{x}_{t}s_{t}^{l}+\mathbf{R}_{t}^{(2)}\mathbf{B}^{\prime}_{3}\mathbf{x}_{t}+\bm{\varepsilon}_{t}$
$\displaystyle=\mathbf{\Psi}_{0}^{\prime}\mathbf{x}_{t}+\sum_{l=1}^{L}\mathbf{\Psi}_{l}^{\prime}\mathbf{x}_{t}s_{t}^{l}+\bm{\varepsilon}_{t}^{*}$
(17)

where
$\mathbf{\Psi}_{0}=\left(\mathbf{B}_{1}^{\prime}+\mathbf{G}_{t}^{(1)}\mathbf{B}_{2}^{\prime}+\mathbf{\Upsilon}_{0}^{(2)}\mathbf{B}^{\prime}_{3}\right)^{\prime}$,
$\mathbf{\Psi}_{l}=\left(\mathbf{\Upsilon}_{l}^{(2)}\mathbf{B}^{\prime}_{3}\right)^{\prime}$
and
$\bm{\varepsilon}_{t}^{*}=\mathbf{R}_{t}^{(2)}\mathbf{B}^{\prime}_{3}\mathbf{x}_{t}+\bm{\varepsilon}_{t}$.

The null hypothesis in the VAR in Eq. (4) is
$\text{H}_{0}\colon\mathbf{\Psi}_{1}=\ldots=\mathbf{\Psi}_{L}=0$. Let be
$\mathbf{Y}=\left(\mathbf{y}_{1}^{\prime},\ldots,\mathbf{y}_{T}^{\prime}\right)^{\prime}$,
$\mathbf{X}=\left(\mathbf{x}_{1}^{\prime},\ldots,\mathbf{x}_{T}^{\prime}\right)^{\prime}$,
$\mathbf{E}$ the $T\times n$ matrix of residuals from Eq. (4), and
$\mathbf{Z}_{L}$ as in Eq. (7), model (4) can be written as

$\mathbf{Y}=\mathbf{X}\mathbf{\Psi}_{0}+\mathbf{Z}_{L}\mathbf{\Psi}_{L}+\mathbf{E}.$
(18)

Let suppose that

$\mathbf{K}=\left[\text{vec}\left(\partial\mathbf{\Psi}_{1}^{\prime}\mathbf{B}^{\prime}\mathbf{x}_{1}/\partial\bm{\theta}\right)^{\prime}\ldots\
\text{vec}\left(\partial\mathbf{\Psi}_{T}^{\prime}\mathbf{B}^{\prime}\mathbf{x}_{T}/\partial\bm{\theta}\right)^{\prime}\right],$
(19)

and that
$\mathbf{P}_{\mathbf{K}}=\mathbf{K}(\mathbf{K^{\prime}K})^{-1}\mathbf{K}^{\prime}$,
the test statistic can be computed similarly to the one in Section 3,
therefore we can state the following result:

######  Theorem 2.

If the estimates of the parameters in Eq. (18) are consistent, under the null
$\text{H}_{0}\colon\mathbf{\Psi}_{L}=\mathbf{0}$, the LM test statistic for
non-additive nonlinearity

$LM_{L}=\text{tr}\left\\{\mathbf{\hat{\Omega}}^{-1}\mathbf{\hat{E}}^{\prime}\mathbf{Z}_{L}\left[\mathbf{Z}_{L}^{\prime}\left(\mathbf{I}_{T}-\mathbf{P}_{K}\right)\mathbf{Z}_{L}\right]^{-1}\mathbf{Z}_{L}^{\prime}\mathbf{\hat{E}}\right\\}$
(20)

has an asymptotic $\chi^{2}$-distribution with $Ln(1+np)$ degrees of freedom
under the Assumptions 1-3 from Section 3, and under the assumption that
$\mathbf{K}^{\prime}\mathbf{K}$ and
$\mathbf{Z}_{L}^{\prime}(\mathbf{I}_{T}-\mathbf{P}_{K})\mathbf{Z}_{L}$ are
positive definite matrices.

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 (Eitrheim and Teräsvirta,,
1996). 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. (19), the vector of first-order partial derivatives of
$\mathbf{\Psi}_{t}^{\prime}\mathbf{B}^{\prime}\mathbf{x}_{t}$ is necessary.
This can be found in Appendix C. As further discussed in Section 3, 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.

###  4.1 Applying the sequential procedure for a vector threshold
autoregressive model

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 Giannerini et
al., (2024). Let’s suppose a 2-regime VTAR model with a single transition
variable and a not-switching error term, defined as follows

$\mathbf{y}_{t}=\left(\bm{\mu}_{1}+\sum_{j=1}^{p}\mathbf{\Phi}_{1,j}\mathbf{y}_{t-j}\right)\mathbbm{1}\left(s_{t}\leq
c\right)+\left(\bm{\mu}_{2}+\sum_{j=1}^{p}\mathbf{\Phi}_{2,j}\mathbf{y}_{t-j}\right)\mathbbm{1}\left(s_{t}>c\right)+\bm{\varepsilon}_{t}$

which can alternatively be written as

$\mathbf{y}_{t}=\bm{\mu}_{1}^{*}+\sum_{j=1}^{p}\mathbf{\Phi}^{*}_{1,j}\mathbf{y}_{t-j}+\left(\mu_{2}^{*}+\sum_{j=1}^{p}\mathbf{\Phi}_{2,j}^{*}\mathbf{y}_{t-j}\right)\mathbbm{1}(s_{t}>c)+\bm{\varepsilon}_{t}.$
(21)

It follows that the indicator function $\mathbbm{1}(\cdot)$ in Eq. (21) 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 (see also Luukkonen et al.,, 1988). 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}$)

$\mathbf{y}_{t}=\left(\mathbf{\Phi}_{1}\mathbf{y}_{t-1}\right)\mathbbm{1}(s_{t}\leq
c_{1})+\left(\mathbf{\Phi}_{2}\mathbf{y}_{t-1}\right)\mathbbm{1}(c_{1}<s_{t}\leq
c_{2})+\left(\mathbf{\Phi}_{3}\mathbf{y}_{t-1}\right)\mathbbm{1}(s_{t}>c_{2})+\bm{\varepsilon}_{t}$

which can be reformulated as

$\mathbf{y}_{t}=\mathbf{\Phi}_{1}^{*}\mathbf{y}_{t-1}+\left(\mathbf{\Phi}_{2}^{*}\mathbf{y}_{t-1}\right)\mathbbm{1}(s_{t}>c_{1})+\left(\mathbf{\Phi}_{3}^{*}\mathbf{y}_{t-1}\right)\mathbbm{1}(s_{t}>c_{2})+\bm{\varepsilon}_{t}$
(22)

with $c_{1}<c_{2}$. The sequential test procedure can be summarized in the
following steps:

  1. 1.

Set $\mathbf{\Phi}_{3}^{*}=\mathbf{0}$ in Eq. (22) and approximate the
indicator function $\mathbbm{1}(s_{t}>c_{1})$ with
$\mathbf{G}_{t}^{(1)}(s_{t};\bm{\gamma}_{1},\mathbf{c}_{1})=g_{t}^{(1)}(s_{t};\gamma_{1},c_{1})\mathbf{I}_{n}$,
where $g_{t}^{(1)}$ is a standard logistic function, such that

$\mathbf{y}_{t}=\mathbf{\Phi}_{1}^{*}\mathbf{y}_{t-1}+\left(\mathbf{\Phi}_{2}^{*}\mathbf{y}_{t-1}\right)\mathbf{G}_{t}^{(1)}(s_{t};\bm{\gamma}_{1},\mathbf{c}_{1})+\bm{\varepsilon}_{t}.$

Then, apply the linearity test presented in Section 3 for a VLSTAR model,
imposing $H_{0}\colon\gamma_{1}=0$.

  2. 2.

If the null hypothesis is rejected at a given significance level $\alpha$,
estimate the coefficients in (22) model imposing
$\mathbf{\Phi}_{3}^{*}=\mathbf{0}$. As further detailed in Chan, (1993) and
Gonzalo and Pitarakis, (2002), the threshold estimator, $\hat{c}_{1}$, is
super consistent.

  3. 3.

Use the super consistent estimator, $\hat{c}_{1}$, in Eq. (22) and test the
linearity of the following model

$\mathbf{y}_{t}=\mathbf{\Phi}_{1}^{*}\mathbf{y}_{t-1}+\left(\mathbf{\Phi}_{2}^{*}\mathbf{y}_{t-1}\right)\mathbbm{1}(s_{t}>\hat{c}_{1})+\left(\mathbf{\Phi}_{3}^{*}\mathbf{y}_{t-1}\right)\mathbf{G}_{t}^{(2)}(s_{t};\bm{\gamma}_{2},\mathbf{c}_{2})+\bm{\varepsilon}_{t}$

where
$\mathbf{G}_{t}^{(2)}(s_{t};\bm{\gamma}_{2},\mathbf{c}_{2})=g_{t}^{(2)}(s_{t};\gamma_{2},c_{2})\mathbf{I}_{n}$.

  4. 4.

If the null hypothesis, $H_{0}\colon\gamma_{2}=0$, is rejected at a
significance level, estimate the parameters in Eq. (22) and use the super
consistent estimate, $\hat{c}_{2}$, to test the linearity of (22) against a
four-regime model.

  5. 5.

Continue the procedure until a non-rejection.

This, indeed, extends in the multivariate the sequential procedure for the
definition of the number of thresholds proposed by Strikholm and Teräsvirta,
(2006).

##  5 Simulation Study

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 Honda, (1988). As in Laitinen, (1978) and Meisner, (1979),
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 Bera et al., (1981), 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
Laitinen, (1978) and Meisner, (1979). The $F$-type LM test statistic, or
rescaled LM test statistic, can be computed as

$LM_{L}^{\text{resc}}=LM_{L}\cdot\frac{nT-S}{W\times nT}$ (23)

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 (Anderson,, 2003). In Appendix D, 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. (15). We do not include
any explanatory variable and we select a single lag for the simulation of
$n=3$ dependent variables, such that:

$\mathbf{y}_{t}=\mathbf{B}_{1}^{\prime}\mathbf{y}_{t-1}+\mathbf{G}_{t}^{(1)}\mathbf{B}_{2}^{\prime}\mathbf{y}_{t-1}+\mathbf{G}_{t}^{(2)}\mathbf{B}_{3}^{\prime}\mathbf{y}_{t-1}+\ldots+\mathbf{G}_{t}^{(m-1)}\mathbf{B}_{m1}^{\prime}\mathbf{y}_{t-1}+\bm{\varepsilon}_{t}.$
(24)

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. (24) with $m=2$ and we test the null
hypothesis of $m=2$ against $m=3$. The results are reported in Table 1. 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

$s_{t}=0.95s_{t-1}+\eta_{t}$ (25)

where $\eta_{t}\sim\mathcal{N}(0,1)$. As in Strikholm and Teräsvirta, (2006),
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 1, 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. (24) 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 2 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.

Table 1: Empirical size of additive nonlinearity test. The empirical size is in per cent based on 1000 replications from model (24) with $m=2$ regimes and $n=3$ | $LM_{3}$ | $LM_{3}^{\text{resc}}$ | Wilks
---|---|---|---
| 10% | 5% | 1% | 10% | 5% | 1% | 10% | 5% | 1%
$T=400$ |  |  |  |  |  |  |  |  |
$\rho_{i}\sim\mathcal{U}(0.3,0.5)$ | 14.1 | 7.6 | 1.6 | 9.9 | 5.1 | 0.7 | 14.0 | 7.8 | 1.6
$\rho_{i}\sim\mathcal{U}(0.5,0.8)$ | 15.1 | 8.6 | 1.7 | 11.2 | 5.9 | 1.2 | 14.9 | 8.8 | 1.8
$T=600$ |  |  |  |  |  |  |  |  |
$\rho_{i}\sim\mathcal{U}(0.3,0.5)$ | 10.4 | 6.3 | 1.6 | 8.8 | 5.0 | 1.1 | 10.5 | 6.3 | 1.6
$\rho_{i}\sim\mathcal{U}(0.5,0.8)$ | 12.4 | 7.1 | 2.4 | 10.2 | 5.4 | 2.0 | 12.3 | 6.9 | 2.4
$T=1000$ |  |  |  |  |  |  |  |  |
$\rho_{i}\sim\mathcal{U}(0.3,0.5)$ | 9.1 | 4.9 | 1.5 | 7.8 | 4.4 | 0.9 | 9.0 | 4.9 | 1.5
$\rho_{i}\sim\mathcal{U}(0.5,0.8)$ | 8.3 | 4.1 | 0.7 | 7.4 | 3.6 | 0.6 | 8.3 | 4.2 | 0.7
Note: Data are simulated from a model with $m=2$ regimes and we test the null
of $\text{H}_{0}\colon m=2$ against $H_{1}\colon m=3$. We simulate using
different values of $\mathbf{B}_{1}$ and for $\mathbf{B}_{2}=-\mathbf{B}_{1}$.
Table 2: Power of the additive nonlinearity test at different sample sizes based on 1000 replications from model (24) with $m=3$ regimes and $n=3$, testing the null $\text{H}_{0}\colon m=2$ |  $LM_{3}$ |  $LM_{3}^{\text{resc}}$ |  Wilks
---|---|---|---
$T=400$ |  |  |
$\rho_{i}\sim\mathcal{U}(0.3,0.5)$ |  85.7 |  83.0 |  84.8
$\rho_{i}\sim\mathcal{U}(0.5,0.8)$ |  89.5 |  86.2 |  89.0
$T=600$ |  |  |
$\rho_{i}\sim\mathcal{U}(0.3,0.5)$ |  89.4 |  87.5 |  88.6
$\rho_{i}\sim\mathcal{U}(0.5,0.8)$ |  92.7 |  91.8 |  92.7
$T=1000$ |  |  |
$\rho_{i}\sim\mathcal{U}(0.3,0.5)$ |  93.8 |  93.2 |  93.6
$\rho_{i}\sim\mathcal{U}(0.5,0.8)$ |  96.1 |  95.3 |  96.0
Note: Data are simulated from a model with $m=3$ regimes and we test the null
of $\text{H}_{0}\colon m=2$ against $H_{1}\colon m=3$. We simulate using
different values of $\mathbf{B}_{1}$, for $\mathbf{B}_{2}=-\mathbf{B}_{1}$ and
$\mathbf{B}_{3}=-0.7\cdot\mathbf{I}_{n}$, with $\gamma_{1}=\gamma_{2}=2$,
$c_{1}=2$ and $c_{2}=4$. The nominal significance level is set equal to 0.05.

In order to evaluate the regime choice of the procedures introduced in this
paper, we also report in Table 3 the selection frequencies when the data are
simulated, as before, from model (24) 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$.

Table 3: Selection frequencies for a VLSTAR model. The frequency is in per cent based on 1000 replications from model (24) with $m=2$ regimes and $n=3$ | $LM_{3}$ | $LM_{3}^{\text{resc}}$ | Wilks
---|---|---|---
| $\hat{m}=1$ | $\hat{m}=2$ | $\hat{m}\geq 3$ | $\hat{m}=1$ | $\hat{m}=2$ | $\hat{m}\geq 3$ | $\hat{m}=1$ | $\hat{m}=2$ | $\hat{m}\geq 3$
$T=400,\rho_{i}\sim\mathcal{U}(0.3,0.5)$ |  |  |  |  |  |  |
$\alpha=0.10$ | 0.0 | 85.9 | 14.1 | 0.0 | 90.1 | 9.9 | 11.7 | 74.8 | 13.5
$\alpha=0.05$ | 0.1 | 92.3 | 7.6 | 0.1 | 94.8 | 5.1 | 17.8 | 74.7 | 7.5
$\alpha=0.01$ | 0.4 | 98.0 | 1.6 | 0.7 | 98.6 | 0.7 | 34.9 | 63.8 | 1.3
$T=400,\rho_{i}\sim\mathcal{U}(0.5,0.8)$ |  |  |  |  |  |  |  |  |
$\alpha=0.10$ | 0.0 | 84.9 | 15.1 | 0.0 | 88.8 | 11.2 | 0.5 | 84.6 | 14.9
$\alpha=0.05$ | 0.0 | 91.4 | 8.6 | 0.0 | 94.1 | 5.9 | 1.1 | 90.1 | 8.8
$\alpha=0.01$ | 0.0 | 98.3 | 1.7 | 0.0 | 98.8 | 1.2 | 1.9 | 96.3 | 1.8
$T=600,\rho_{i}\sim\mathcal{U}(0.3,0.5)$ |  |  |  |  |  |  |  |  |
$\alpha=0.10$ | 0.0 | 89.6 | 10.4 | 0.0 | 91.2 | 8.8 | 3.0 | 86.7 | 10.3
$\alpha=0.05$ | 0.0 | 93.7 | 6.3 | 0.0 | 95.0 | 5.0 | 4.7 | 89.0 | 6.3
$\alpha=0.01$ | 0.0 | 98.4 | 1.6 | 0.0 | 98.9 | 1.1 | 10.1 | 88.3 | 1.6
$T=600,\rho_{i}\sim\mathcal{U}(0.5,0.8)$ |  |  |  |  |  |  |  |  |
$\alpha=0.10$ | 0.0 | 87.6 | 12.4 | 0.0 | 89.8 | 10.2 | 0.1 | 87.6 | 12.3
$\alpha=0.05$ | 0.0 | 92.9 | 7.1 | 0.0 | 94.6 | 5.4 | 0.1 | 93.0 | 6.9
$\alpha=0.01$ | 0.0 | 97.6 | 2.4 | 0.0 | 98.0 | 2.0 | 0.5 | 97.1 | 2.4
$T=1000,\rho_{i}\sim\mathcal{U}(0.3,0.5)$ |  |  |  |  |  |  |  |  |
$\alpha=0.10$ | 0.0 | 90.9 | 9.1 | 0.0 | 92.2 | 7.8 | 0.0 | 91.0 | 9.0
$\alpha=0.05$ | 0.0 | 95.1 | 4.9 | 0.0 | 95.6 | 4.4 | 0.0 | 95.1 | 4.9
$\alpha=0.01$ | 0.0 | 98.5 | 1.5 | 0.0 | 99.1 | 0.9 | 0.2 | 98.3 | 1.5
$T=1000,\rho_{i}\sim\mathcal{U}(0.5,0.8)$ |  |  |  |  |  |  |  |  |
$\alpha=0.10$ | 0.0 | 91.7 | 8.3 | 0.0 | 92.6 | 7.4 | 0.0 | 91.7 | 8.3
$\alpha=0.05$ | 0.0 | 95.9 | 4.1 | 0.0 | 96.4 | 3.6 | 0.0 | 95.8 | 4.2
$\alpha=0.01$ | 0.0 | 99.3 | 0.7 | 0.0 | 99.4 | 0.6 | 0.0 | 99.3 | 0.7
Note: Data are simulated from a model with $m=2$ regimes. We simulate using
different values of $\mathbf{B}_{1}$ and for $\mathbf{B}_{2}=-\mathbf{B}_{1}$.
Underlined values denote the identified number of regimes.

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

$\mathbf{y}_{t}=\left(\mathbf{\Phi}_{1}\mathbf{y}_{t-1}+\bm{\varepsilon}_{1,t}\right)\mathbbm{1}(s_{t}<c)+\left(\mathbf{\Phi}_{2}\mathbf{y}_{t-1}+\bm{\varepsilon}_{2,t}\right)\mathbbm{1}(s_{t}\geq
c)$ (26)

where $s_{t}$ is generated from an autoregressive model as in Eq. (25). 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 4, we evaluate the empirical sizes of the test applied to the data
simulated from a VTAR specified as in Eq. (26). 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 5 (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 6 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.

Table 4: Empirical size of additive nonlinearity test for a VTAR model. The empirical size is in per cent based on 1000 replications from model (26) with $m=2$ regimes and $n=3$ | $LM_{3}$ | $LM_{3}^{\text{resc}}$ | Wilks
---|---|---|---
| 10% | 5% | 1% | 10% | 5% | 1% | 10% | 5% | 1%
$T=400$ |  |  |  |  |  |  |  |  |
$\rho_{i}\sim\mathcal{U}(0.3,0.5)$ | 12.6 | 6.9 | 0.7 | 9.2 | 4.3 | 0.3 | 12.2 | 6.8 | 0.9
$\rho_{i}\sim\mathcal{U}(0.5,0.8)$ | 15.6 | 8.5 | 2.6 | 10.4 | 5.8 | 1.6 | 15.1 | 8.6 | 2.6
$T=600$ |  |  |  |  |  |  |  |  |
$\rho_{i}\sim\mathcal{U}(0.3,0.5)$ | 10.6 | 5.6 | 0.7 | 8.2 | 3.3 | 0.6 | 10.5 | 5.5 | 0.7
$\rho_{i}\sim\mathcal{U}(0.5,0.8)$ | 13.5 | 7.6 | 1.8 | 11.1 | 5.4 | 1.4 | 13.4 | 7.5 | 1.9
$T=1000$ |  |  |  |  |  |  |  |  |
$\rho_{i}\sim\mathcal{U}(0.3,0.5)$ | 11.2 | 5.9 | 1.6 | 9.4 | 4.9 | 1.4 | 11.1 | 5.9 | 1.6
$\rho_{i}\sim\mathcal{U}(0.5,0.8)$ | 13.1 | 6.2 | 1.5 | 10.6 | 5.2 | 1.1 | 13.0 | 6.1 | 1.5
Note: Data are simulated from a VTAR model with $m=2$ regimes and we test the
null of $\text{H}_{0}\colon m=2$ against $H_{1}\colon m=3$. We simulate using
different values of $\bm{\Phi}_{1}$ and for $\bm{\Phi}_{2}=-\bm{\Phi}_{1}$.
Table 5: Power of the additive nonlinearity test at different sample sizes based on 1000 replications from model (26) with $m=3$ regimes and $n=3$, testing the null $\text{H}_{0}\colon m=2$ |  $LM_{3}$ |  $LM_{3}^{\text{resc}}$ |  Wilks
---|---|---|---
$T=400$ |  |  |
$\rho_{i}\sim\mathcal{U}(0.3,0.5)$ |  93.2 |  91.6 |  76.6
$\rho_{i}\sim\mathcal{U}(0.5,0.8)$ |  94.6 |  92.6 |  80.3
$T=600$ |  |  |
$\rho_{i}\sim\mathcal{U}(0.3,0.5)$ |  94.7 |  94.0 |  80.6
$\rho_{i}\sim\mathcal{U}(0.5,0.8)$ |  95.7 |  95.3 |  84.2
$T=1000$ |  |  |
$\rho_{i}\sim\mathcal{U}(0.3,0.5)$ |  95.5 |  95.3 |  86.7
$\rho_{i}\sim\mathcal{U}(0.5,0.8)$ |  97.2 |  96.9 |  85.1
Note: Data are simulated from a model with $m=3$ regimes and we test the null
of $\text{H}_{0}\colon m=2$ against $H_{1}\colon m=3$. We simulate using
different values of $\bm{\Phi}_{1}$, for $\bm{\Phi}_{2}=-\bm{\Phi}_{1}$ and
$\bm{\Phi}_{3}=-0.7\cdot\mathbf{I}_{n}$, with $c_{1}=2$ and $c_{2}=4$. The
nominal significance level is set equal to 0.05.
Table 6: Selection frequencies for a VTAR model. The frequency is in per cent based on 1000 replications from model (26) with $m=2$ regimes and $n=3$ | $LM_{3}$ | $LM_{3}^{\text{resc}}$ | Wilks
---|---|---|---
| $\hat{m}=1$ | $\hat{m}=2$ | $\hat{m}\geq 3$ | $\hat{m}=1$ | $\hat{m}=2$ | $\hat{m}\geq 3$ | $\hat{m}=1$ | $\hat{m}=2$ | $\hat{m}\geq 3$
$T=400,\rho_{i}\sim\mathcal{U}(0.3,0.5)$ |  |  |  |  |  |  |
$\alpha=0.10$ | 0.0 | 87.4 | 12.6 | 0.0 | 90.8 | 9.2 | 0.0 | 87.8 | 12.2
$\alpha=0.05$ | 0.0 | 93.1 | 6.9 | 0.0 | 95.7 | 4.3 | 0.0 | 93.2 | 6.8
$\alpha=0.01$ | 0.0 | 99.3 | 0.7 | 0.0 | 99.7 | 0.3 | 0.0 | 99.1 | 0.9
$T=400,\rho_{i}\sim\mathcal{U}(0.8)$ |  |  |  |  |  |  |  |  |
$\alpha=0.10$ | 0.0 | 84.4 | 15.6 | 0.0 | 89.6 | 10.4 | 0.0 | 84.9 | 15.1
$\alpha=0.05$ | 0.0 | 91.5 | 8.5 | 0.0 | 94.2 | 5.8 | 0.0 | 91.4 | 8.6
$\alpha=0.01$ | 0.0 | 97.4 | 2.6 | 0.0 | 98.4 | 1.6 | 0.0 | 97.4 | 2.6
$T=600,\rho_{i}\sim\mathcal{U}(0.3,0.5)$ |  |  |  |  |  |  |  |  |
$\alpha=0.10$ | 0.0 | 89.4 | 10.6 | 0.0 | 91.8 | 8.2 | 0.0 | 89.5 | 10.5
$\alpha=0.05$ | 0.0 | 94.4 | 5.6 | 0.0 | 96.7 | 3.3 | 0.0 | 94.5 | 5.5
$\alpha=0.01$ | 0.0 | 99.3 | 0.7 | 0.0 | 99.4 | 0.6 | 0.1 | 99.2 | 0.7
$T=600,\rho_{i}\sim\mathcal{U}(0.5,0.8)$ |  |  |  |  |  |  |  |  |
$\alpha=0.10$ | 0.0 | 86.5 | 13.5 | 0.0 | 88.9 | 11.1 | 0.0 | 86.6 | 13.4
$\alpha=0.05$ | 0.0 | 92.4 | 7.6 | 0.0 | 94.6 | 5.4 | 0.0 | 92.5 | 7.5
$\alpha=0.01$ | 0.0 | 98.2 | 1.8 | 0.0 | 98.6 | 1.4 | 0.0 | 98.1 | 1.9
$T=1000,\rho_{i}\sim\mathcal{U}(0.3,0.5)$ |  |  |  |  |  |  |  |  |
$\alpha=0.10$ | 0.0 | 88.8 | 11.2 | 0.0 | 90.6 | 9.4 | 0.0 | 88.9 | 11.1
$\alpha=0.05$ | 0.0 | 94.1 | 5.9 | 0.0 | 95.1 | 4.9 | 0.0 | 94.1 | 5.9
$\alpha=0.01$ | 0.0 | 98.4 | 1.6 | 0.0 | 98.6 | 1.4 | 0.0 | 98.4 | 1.6
$T=1000,\rho_{i}\sim\mathcal{U}(0.5,0.8)$ |  |  |  |  |  |  |  |  |
$\alpha=0.10$ | 0.0 | 86.9 | 13.1 | 0.0 | 89.4 | 10.6 | 0.0 | 87.0 | 13.0
$\alpha=0.05$ | 0.0 | 93.8 | 6.2 | 0.0 | 94.8 | 5.2 | 0.0 | 93.9 | 6.1
$\alpha=0.01$ | 0.0 | 98.5 | 1.5 | 0.0 | 98.9 | 1.1 | 0.0 | 98.5 | 1.5
Note: Data are simulated from a model with $m=2$ regimes. We simulate using
different values of $\bm{\Phi}_{1}$ and for $\bm{\Phi}_{2}=-\bm{\Phi}_{1}$.
Underlined values denote the identified number of regimes.

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)$ (stationarity was assessed through the
method proposed in Kheifets and Saikkonen,, 2020).

Table 7 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 8. 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.

Table 7: Empirical size of additive nonlinearity test. The empirical size is in per cent based on 1000 replications from a model with $m=2$ regimes and $n=5$ | $LM_{3}$ | $LM_{3}^{\text{resc}}$ | Wilks
---|---|---|---
| 10% | 5% | 1% | 10% | 5% | 1% | 10% | 5% | 1%
| Panel A: VLSTAR
$T=400$ |  |  |  |  |  |  |  |  |
$\rho_{i}\sim\mathcal{U}(0.3,0.5)$ | 7.2 | 3.7 | 1.4 | 3.3 | 1.6 | 0.2 | 7.0 | 3.7 | 1.4
$T=600$ |  |  |  |  |  |  |  |  |
$\rho_{i}\sim\mathcal{U}(0.3,0.5)$ | 7.8 | 5.0 | 1.8 | 5.0 | 3.0 | 1.4 | 7.6 | 5.0 | 1.6
$T=1000$ |  |  |  |  |  |  |  |  |
$\rho_{i}\sim\mathcal{U}(0.3,0.5)$ | 3.2 | 1.4 | 0.4 | 2.4 | 1.0 | 0.2 | 3.2 | 1.2 | 0.4
| Panel B: VTAR
$T=400$ |  |  |  |  |  |  |  |  |
$\rho_{i}\sim\mathcal{U}(0.3,0.5)$ | 10.3 | 7.4 | 3.7 | 5.9 | 3.9 | 1.2 | 10.3 | 7.2 | 3.7
$T=600$ |  |  |  |  |  |  |  |  |
$\rho_{i}\sim\mathcal{U}(0.3,0.5)$ | 16.3 | 9.8 | 3.8 | 16.9 | 9.7 | 2.9 | 16.7 | 10.7 | 2.7
$T=1000$ |  |  |  |  |  |  |  |  |
$\rho_{i}\sim\mathcal{U}(0.3,0.5)$ | 15.6 | 9.5 | 3.4 | 12.0 | 7.1 | 2.1 | 15.1 | 9.4 | 3.5
Note: Data are simulated from a model with $m=2$ regimes and we test the null
of $\text{H}_{0}\colon m=2$ against $H_{1}\colon m=3$. We simulate using
different values of $\mathbf{B}_{1}$ ($\bm{\Phi}_{1}$) and for
$\mathbf{B}_{2}=-\mathbf{B}_{1}$ ($\bm{\Phi}_{2}=-\bm{\Phi}_{1}$ in the VTAR).
Table 8: Power of the additive nonlinearity test at different sample sizes based on 1000 replications from a model with $m=3$ regimes and $n=5$, testing the null $\text{H}_{0}\colon m=2$ |  $LM_{3}$ |  $LM_{3}^{\text{resc}}$ |  Wilks
---|---|---|---
| Panel A: VLSTAR
$T=400$ |  |  |
$\rho_{i}\sim\mathcal{U}(0.3,0.5)$ |  95.2 |  91.6 |  94.0
$T=600$ |  |  |
$\rho_{i}\sim\mathcal{U}(0.3,0.5)$ |  92.4 |  90.8 |  92.2
$T=1000$ |  |  |
$\rho_{i}\sim\mathcal{U}(0.3,0.5)$ |  99.2 |  98.6 |  99.2
| Panel B: VTAR
$T=400$ |  |  |
$\rho_{i}\sim\mathcal{U}(0.3,0.5)$ |  93.2 |  91.6 |  76.6
$T=600$ |  |  |
$\rho_{i}\sim\mathcal{U}(0.3,0.5)$ |  98.9 |  98.3 |  99.8
$T=1000$ |  |  |
$\rho_{i}\sim\mathcal{U}(0.3,0.5)$ |  99.9 |  100.0 |  100.0
Note: Data are simulated from a model with $m=3$ regimes and we test the null
of $\text{H}_{0}\colon m=2$ against $H_{1}\colon m=3$. We simulate using
different values of $\mathbf{B}_{1}$ ($\bm{\Phi}_{1}$), for
$\mathbf{B}_{2}=-\mathbf{B}_{1}$ ($\bm{\Phi}_{2}=-\bm{\Phi}_{1}$ in the VTAR)
and $\mathbf{B}_{3}=-0.7\cdot\mathbf{I}_{n}$
($\bm{\Phi}_{3}=-0.7\cdot\mathbf{I}_{n}$ in the VTAR), with $c_{1}=2$ and
$c_{2}=4$ in both the models, and $\gamma_{1}=\gamma_{2}=2$ in the VLSTAR
model. The nominal significance level is set equal to 0.05.
Table 9: Selection frequencies. The frequency is in per cent based on 1000 replications from a model with $m=2$ regimes and $n=3$ | $LM_{3}$ | $LM_{3}^{\text{resc}}$ | Wilks
---|---|---|---
| $\hat{m}=1$ | $\hat{m}=2$ | $\hat{m}\geq 3$ | $\hat{m}=1$ | $\hat{m}=2$ | $\hat{m}\geq 3$ | $\hat{m}=1$ | $\hat{m}=2$ | $\hat{m}\geq 3$
| Panel A: VLSTAR
$T=400,\rho_{i}\sim\mathcal{U}(0.3,0.5)$ |  |  |  |  |  |  |
$\alpha=0.10$ | 0.0 | 92.8 | 7.2 | 0.0 | 96.7 | 3.3 | 3.7 | 89.3 | 7.0
$\alpha=0.05$ | 0.0 | 96.3 | 3.7 | 0.1 | 98.4 | 1.6 | 5.7 | 90.8 | 3.5
$\alpha=0.01$ | 0.0 | 98.6 | 1.4 | 0.7 | 99.8 | 0.2 | 12.9 | 85.9 | 1.2
$T=600,\rho_{i}\sim\mathcal{U}(0.3,0.5)$ |  |  |  |  |  |  |  |  |
$\alpha=0.10$ | 0.0 | 92.2 | 7.8 | 0.0 | 95.0 | 5.0 | 0.6 | 91.8 | 7.6
$\alpha=0.05$ | 0.0 | 95.0 | 5.0 | 0.0 | 97.0 | 3.0 | 0.8 | 94.2 | 5.0
$\alpha=0.01$ | 0.0 | 98.2 | 1.8 | 0.0 | 98.6 | 1.4 | 2.2 | 96.2 | 1.6
$T=1000,\rho_{i}\sim\mathcal{U}(0.3,0.5)$ |  |  |  |  |  |  |  |  |
$\alpha=0.10$ | 0.0 | 96.8 | 3.2 | 0.0 | 97.6 | 2.4 | 0.0 | 96.8 | 3.2
$\alpha=0.05$ | 0.0 | 98.6 | 1.4 | 0.0 | 99.1 | 1.0 | 0.0 | 98.8 | 1.2
$\alpha=0.01$ | 0.0 | 99.6 | 0.4 | 0.0 | 99.8 | 0.2 | 0.2 | 99.6 | 0.4
| Panel B: VTAR
$T=400,\rho_{i}\sim\mathcal{U}(0.3,0.5)$ |  |  |  |  |  |  |
$\alpha=0.10$ | 0.0 | 89.7 | 10.3 | 0.0 | 94.1 | 5.9 | 0.0 | 89.7 | 10.3
$\alpha=0.05$ | 0.0 | 92.6 | 7.4 | 0.0 | 96.1 | 3.9 | 0.0 | 92.8 | 7.2
$\alpha=0.01$ | 0.0 | 96.3 | 3.7 | 0.0 | 98.8 | 1.2 | 0.0 | 96.3 | 3.7
$T=600,\rho_{i}\sim\mathcal{U}(0.3,0.5)$ |  |  |  |  |  |  |  |  |
$\alpha=0.10$ | 0.0 | 76.7 | 23.3 | 0.0 | 83.1 | 16.9 | 0.0 | 77.3 | 22.7
$\alpha=0.05$ | 0.0 | 83.2 | 16.8 | 0.0 | 88.3 | 11.7 | 0.0 | 83.3 | 16.7
$\alpha=0.01$ | 0.0 | 92.2 | 7.8 | 0.0 | 94.7 | 5.3 | 0.0 | 92.3 | 7.7
$T=1000,\rho_{i}\sim\mathcal{U}(0.3,0.5)$ |  |  |  |  |  |  |  |  |
$\alpha=0.10$ | 0.0 | 84.4 | 15.6 | 0.0 | 88.0 | 12.0 | 0.0 | 84.9 | 15.1
$\alpha=0.05$ | 0.0 | 90.5 | 9.5 | 0.0 | 92.9 | 7.1 | 0.0 | 90.6 | 9.4
$\alpha=0.01$ | 0.0 | 96.6 | 3.4 | 0.0 | 97.9 | 2.1 | 0.2 | 96.5 | 3.5
Note: Data are simulated from a model with $m=2$ regimes. We simulate using
different values of $\mathbf{B}_{1}$ ($\bm{\Phi}_{1}$) and for
$\mathbf{B}_{2}=-\mathbf{B}_{1}$ ($\bm{\Phi}_{2}=-\bm{\Phi}_{1}$ in the VTAR).
Underlined values denote the identified number of regimes.

##  6 Empirical applications

###  6.1 Interest rate term structure

We first apply the sequential test procedure to the U.S. monthly interest
rates data already used in Tsay, (1998). 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)^{\prime}$,
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. 1.

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 (Harvey,, 1988), this would reflect the business cycle of the U.S.
economy. Following Tsay, (1998), 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. 1). 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 Camacho, (2004).

The results of the tests are reported in Table 10. 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
Camacho, (2004) which rejects at 5 and 10%. This means that, in line with what
was supposed and proved by Tsay, (1998), 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 Camacho, (2004).
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 Tsay, (1998).

Table 10: Tests for linearity and additive nonlinearity in the empirical application with interest rates |  $H_{0}\colon m=1$ |  $H_{0}\colon m=2$ |  $H_{0}\colon m=3$
---|---|---|---
|  Test statistic ($p$-value) |  Test statistic ($p$-value) |  Test statistic ($p$-value)
$LM_{3}$ |  $\underset{(1.11e-15)}{98.399}$ |  $\underset{(2.08e-18)}{107.819}$ |  $\underset{(0.770)}{8.185}$
$LM_{3}^{\text{resc}}$ |  $\underset{(4.89e-15)}{8.141}$ |  $\underset{(1.11e-16)}{8.920}$ |  $\underset{(0.775)}{0.677}$
Wilks |  $\underset{(3.33e-16)}{101.087}$ |  $\underset{(3.54e-09)}{64.381}$ |  $\underset{(0.875)}{6.766}$
Camacho |  $\underset{(0.011)}{25.805}$ |  $\underset{(0.966)}{0.265}$ |  -
Note: In the first column the linearity test is reported, while columns two
and three belong to the sequential procedure for the identification of the
number of regimes. $p$-values are reported between parentheses.
Figure 1: Time series plots of U.S. interest rates. Note: The time series of
the monthly 3-bill rate growth (blue line) and the 3-year Treasury rate growth
(orange line) highlight the possible existence of several regimes. The 3-month
moving average of the maturity spread (green line) is used as a transition
variable.

###  6.2 River flows data

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ökulsá and the Vatndalsá, as well
as the temperature and the precipitation, see Figure 2. River flow data has
been shown to be nonlinear in several former applications. For instance, Tong
et al., (1985) use a univariate threshold model to estimate their relationship
with temperature and precipitation, while Tsay, (1998), Teräsvirta and Yang,
2014b , and Livingston Jr and Nur, (2020) apply a multivariate nonlinear
model.

Figure 2: Time series plots of Icelandic river flow data. Note: The time
series of the Jökulsá (blue line) and the Vatndalsá (orange line) river flows
highlight a nonlinear path for both rivers. The temperature (green line) and
the precipitation level (red line) can be used as a transition variable.

Following Tsay, (1998) and Teräsvirta and Yang, 2014b , we first select the
lagged temperature as the transition variable for both flow equations (Panel A
of Table 11) 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
Camacho, (2004) 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 Teräsvirta and Yang,
2014b .

We also apply the tests using the lagged precipitation as a transition
variable (Panel B of Table 11). 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.

Table 11: Tests for linearity and additive nonlinearity in the empirical application with river flows | Panel A: Lagged temperature as the transition variable
---|---
|  $H_{0}\colon m=1$ |  $H_{0}\colon m=2$ |  $H_{0}\colon m=3$
|  Test statistic ($p$-value) |  Test statistic ($p$-value) |  Test statistic ($p$-value)
$LM_{3}$ |  $\underset{(2.22e-15)}{89.301}$ |  $\underset{(0.021)}{23.939}$ |  $\underset{(0.143)}{17.174}$
$LM_{3}^{\text{resc}}$ |  $\underset{(2.31e-15)}{21.098}$ |  $\underset{(0.022)}{1.984}$ |  $\underset{(0.148)}{1.423}$
Wilks |  $\underset{(4.44e-16)}{100.663}$ |  $\underset{(0.004)}{29.333}$ |  $\underset{(0.065)}{20.095}$
Camacho |  $\underset{(1.88e-10)}{71.231}$ |  $\underset{(0.289)}{14.184}$ |  -
| Panel B: Lagged precipitation as the transition variable
|  $H_{0}\colon m=1$ |  $H_{0}\colon m=2$ |  $H_{0}\colon m=3$
|  Test statistic ($p$-value) |  Test statistic ($p$-value) |  Test statistic ($p$-value)
$LM_{3}$ |  $\underset{(4.14e-14)}{90.393}$ |  $\underset{(0.719)}{8.810}$ |  -
$LM_{3}^{\text{resc}}$ |  $\underset{(1.06e-13)}{7.491}$ |  $\underset{(0.722)}{0.730}$ |  -
Wilks |  $\underset{(2.21e-14)}{91.807}$ |  $\underset{(0.670)}{9.377}$ |  -
Camacho |  $\underset{(0.007)}{27.222}$ |  $\underset{(0.478)}{11.606}$ |  -
Note: In the first column the linearity test is reported, while columns two
and three belong to the sequential procedure for the identification of the
number of regimes. $p$-values are reported between parentheses.

##  7 Conclusions

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 Teräsvirta and
Yang, 2014a .

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 Camacho, (2004). 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
Camacho, (2004) 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.

## Acknowledgements

I offer my sincerest gratitude to Timo Të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.

## References

  * Anderson et al., (1999) Anderson, H. M., Nam, K., and Vahid, F. (1999).  Asymmetric Nonlinear Smooth Transition GARCH Models, pages 191–207.  Springer US, Boston, MA.
  * Anderson, (2003) Anderson, T. W. (2003).  An Introduction to Multivariate Statistical Analysis.  Wiley, New York, 3rd edition.
  * Bai and Perron, (1998) Bai, J. and Perron, P. (1998).  Estimating and Testing Linear Models With Multiple Structural Changes.  Econometrica, 66:47–78.
  * Bartlett, (1954) Bartlett, M. S. (1954).  A Note on the Multiplying Factors for Various $\chi^{2}$ Approximations.  Journal of the Royal Statistical Society. Series B (Methodological), 16(2):296–298.  Full publication date: 1954.
  * Basawa et al., (1976) Basawa, I., Feigin, P., and Heyde, C. (1976).  Asymptotic properties of maximum likelihood estimators for stochastic processes.  Sankhyā, A, 38(3):259–270.
  * Bera et al., (1981) Bera, A. K., Byron, R. P., and Jarque, C. M. (1981).  Further evidence on asymptotic tests for homogeneity and symmetry in large demand systems.  Economics Letters, 8(2):101–105.
  * Breusch and Pagan, (1980) Breusch, T. S. and Pagan, A. R. (1980).  The Lagrange Multiplier Test and its Applications to Model Specification in Econometrics.  The Review of Economic Studies, 47(1):239–253.
  * Camacho, (2004) Camacho, M. (2004).  Vector smooth transition regression models for US GDP and the composite index of leading indicators.  Journal of Forecasting, 23:173–196.
  * Chan, (1993) Chan, K. S. (1993).  Consistency and limiting distribution of the least squares estimator of a threshold autoregressive model.  Annals of Statistics, 21:520–533.
  * Cho et al., (2011) Cho, J. S., Ishida, I., and White, H. (2011).  Revisiting Tests for Neglected Nonlinearity Using Artificial Neural Networks.  Neural Computation, 23(5):1133–1186.
  * Davies, (1987) Davies, R. B. (1987).  Hypothesis Testing when a Nuisance Parameter is Present Only Under the Alternatives.  Biometrika, 74(1):33–43.
  * Eitrheim and Teräsvirta, (1996) Eitrheim, Ø. and Teräsvirta, T. (1996).  Testing the adequacy of smooth transition autoregressive models.  Journal of Econometrics, 74(1):59–75.
  * Giannerini et al., (2024) Giannerini, S., Goracci, G., and Rahbek, A. (2024).  The validity of bootstrap testing for threshold autoregression.  Journal of Econometrics, 239(1):105379.
  * Gonzalo and Pitarakis, (2002) Gonzalo, J. and Pitarakis, J.-Y. (2002).  Estimation and model selection based inference in single and multiple threshold models.  Journal of Econometrics, 110(2):319–352.
  * Hansen, (1999) Hansen, B. (1999).  Testing for linearity.  Journal of Economic Surveys, 13(5):551–576.
  * Hansen, (1996) Hansen, B. E. (1996).  Inference When a Nuisance Parameter Is Not Identified Under the Null Hypothesis.  Econometrica, 64(2):413–430.  Full publication date: Mar., 1996.
  * Harvey, (1988) Harvey, C. R. (1988).  The real term structure and consumption growth.  Journal of Financial Economics, 22(2):305–333.
  * He et al., (2008) He, C., Teräsvirta, T., and González, A. (2008).  Testing Parameter Constancy in Stationary Vector Autoregressive Models Against Continuous Change.  Econometric Reviews, 28(1-3):225–245.
  * Hendry and Mizon, (1998) Hendry, D. F. and Mizon, G. E. (1998).  Exogeneity, causality, and co-breaking in economic policy analysis of a small econometric model of money in the UK.  Empirical Economics, 23(3):267–294.
  * Honda, (1988) Honda, Y. (1988).  A size correction to the Lagrange multiplier test for heteroskedasticity.  Journal of Econometrics, 38(3):375–386.
  * Hubrich and Teräsvirta, (2013) Hubrich, K. and Teräsvirta, T. (2013).  Thresholds and Smooth Transitions in Vector Autoregressive Models, volume 32.
  * Kheifets and Saikkonen, (2020) Kheifets, I. L. and Saikkonen, P. J. (2020).  Stationarity and ergodicity of vector STAR models.  Econometric Reviews, 39(4):407–414.
  * Laitinen, (1978) Laitinen, K. (1978).  Why is demand homogeneity so often rejected?  Economics Letters, 1(3):187–191.
  * Livingston Jr and Nur, (2020) Livingston Jr, G. and Nur, D. (2020).  Bayesian estimation and model selection of a multivariate smooth transition autoregressive model.  Environmetrics, 31(6):e2615.
  * Lütkepohl, (1996) Lütkepohl, H. (1996).  Handbook of Matrices.  Wiley, New York.
  * Luukkonen et al., (1988) Luukkonen, R., Saikkonen, P., and Teräsvirta, T. (1988).  Testing Linearity Against Smooth Transition Autoregressive Models.  Biometrika, 75:491–499.
  * Mardia et al., (1979) Mardia, K., Kent, J., and Bibby, J. (1979).  Multivariate analysis.  Acad. Press.
  * Meisner, (1979) Meisner, J. F. (1979).  The sad fate of the asymptotic Slutsky symmetry test for large systems.  Economics Letters, 2:231–233.
  * Saikkonen, (2008) Saikkonen, P. (2008).  Stability of Regime Switching Error Correction Models under Linear Cointegration.  Econometric Theory, 24(1):294–318.
  * Seong et al., (2022) Seong, D., Cho, J. S., and Teräsvirta, T. (2022).  Comprehensively testing linearity hypothesis using the smooth transition autoregressive model.  Econometric Reviews, 41(8):966–984.
  * Sims, (1980) Sims, C. A. (1980).  Macroeconomics and reality.  Econometrica, 48(1):1–48.  Full publication date: Jan., 1980.
  * Strikholm and Teräsvirta, (2006) Strikholm, B. and Teräsvirta, T. (2006).  A sequential procedure for determining the number of regimes in a threshold autoregressive model.  The Econometrics Journal, 9(3):472–491.  Full publication date: 2006.
  * Teräsvirta, (1994) Teräsvirta, T. (1994).  Specification, estimation and evaluation of smooth transition autoregressive models.  Journal of the American Statistical Association, 89:208–218.
  * (34) Teräsvirta, T. and Yang, Y. (2014a).  Linearity and Misspecification Tests for Vector Smooth Transition Regression Models.  CREATES Research Paper, Aarhus University.
  * (35) Teräsvirta, T. and Yang, Y. (2014b).  Specification, estimation and evaluation of vector smooth transition autoregressive models with applications.  Technical report, Department of Economics and Business Economics, Aarhus University.  CREATES Research Paper, Aarhus University.
  * Tong et al., (1985) Tong, H., Thanoon, B., and Gudmundsson, G. (1985).  Threshold time series modeling of two Icelandic riverflow systems.  JAWRA Journal of the American Water Resources Association, 21(4):651–662.
  * Tsay, (1998) Tsay, R. S. (1998).  Testing and Modeling Multivariate Threshold Models.  Journal of the American Statistical Association, 93(443):1188–1202.  Full publication date: Sep., 1998.
  * Wang et al., (2022) Wang, G., Zhu, K., and Shao, X. (2022).  Testing for the Martingale Difference Hypothesis in Multivariate Time Series Models.  Journal of Business & Economic Statistics, 40(3):980–994.
  * Watson and Engle, (1985) Watson, M. W. and Engle, R. F. (1985).  Testing for Regression Coefficient Stability with a Stationary AR(1) Alternative.  The Review of Economics and Statistics, 67:341–346.
  * White and Cho, (2012) White, H. and Cho, J. S. (2012).  Higher-Order Approximations for Testing Neglected Nonlinearity.  Neural Computation, 24(1):273–287.

##  Appendix A The score vector

Let’s assume a 3-regime model with a single lag

$\mathbf{y}_{t}=\bm{\mu}_{0}+\mathbf{\Phi}_{0}\mathbf{y}_{t-1}+\mathbf{G}_{t}^{(1)}\left(\bm{\mu}_{1}+\mathbf{\Phi}_{1}\mathbf{y}_{t-1}\right)+\mathbf{G}_{t}^{(2)}\left(\bm{\mu}_{2}+\mathbf{\Phi}_{2}\mathbf{y}_{t-1}\right)+\bm{\varepsilon}_{t}$
(27)

which can be written as

$\mathbf{y}_{t}=\mathbf{\Psi}_{t}^{\prime}\mathbf{B}^{\prime}\mathbf{x}_{t}+\mathbf{G}_{t}^{(2)}\mathbf{B}_{2}\mathbf{x}_{t}+\bm{\varepsilon}_{t}$
(28)

where
$\mathbf{\Psi}=\left(\mathbf{I}_{n},\mathbf{G}_{t}^{(1)}\right)^{\prime}$ is a
$2n\times n$ matrix, $\mathbf{B}=\left(\mathbf{B}_{0},\mathbf{B}_{1}\right)$
is a $(1+n)\times 2n$ matrix of parameters,
$\mathbf{x}_{t}=\left[1,\mathbf{y}_{t-1}^{\prime}\right]^{\prime}$ is a
$(n+1)\times 1$ vector, and
$\mathbf{B}_{j}=\left(\bm{\mu}_{j}^{\prime},\mathbf{\Phi}_{j}^{\prime}\right)^{\prime}$,
for $j=0,1,2$.

To derive the score and the relevant derivatives, we first assume a single
transition model, i.e. excluding
$\mathbf{G}_{t}^{(2)}\left(\bm{\mu}_{2}+\mathbf{\Phi}_{2}\mathbf{y}_{t-1}\right)$
from model (27). The score of the log-likelihood for model (27) is

$\mathbf{s}_{t}(\bm{\theta})=\frac{\partial\ell_{t}(\bm{\theta})}{\partial\bm{\theta}}=-\frac{\partial\bm{\varepsilon}_{t}^{\prime}}{\partial\bm{\theta}}\mathbf{\Omega}^{-1}\bm{\varepsilon}_{t}.$

Now, we have to compute
$\partial\bm{\varepsilon}_{t}^{\prime}/\partial\bm{\theta}$, where

$\bm{\varepsilon}_{t}=\mathbf{y}_{t}-\bm{\mu}_{0}-\mathbf{\Phi}_{0}\mathbf{y}_{t-1}-\mathbf{G}_{t}(\bm{\mu}_{1}+\mathbf{\Phi}_{1}\mathbf{y}_{t-1}).$

The $j$-th equation (out of $n$) is

$\displaystyle\bm{e}_{j}\bm{\varepsilon}_{t}$
$\displaystyle=\bm{e}_{j}^{\prime}\left(\mathbf{y}_{t}-\bm{\mu}_{0}-\mathbf{\Phi}_{0}\mathbf{y}_{t-1}-\mathbf{G}_{t}(\bm{\mu}_{1}+\mathbf{\Phi}_{1}\mathbf{y}_{t-1})\right)$
$\displaystyle=\bm{e}_{j}^{\prime}\left(\mathbf{y}_{t}-\bm{\mu}_{0}-\sum_{i=1}^{n}\phi_{0i}y_{i,t-1}-\mathbf{G}_{t}(\bm{\mu}_{1}+\sum_{i=1}^{n}\phi_{1i}y_{i,t-1})\right)$

where $\mathbf{e}_{j}$ is the $j$-th column of the $n\times n$ identity
matrix. Denote the parameters of the $i$-th equation as $\bm{\theta}_{i}$.
Then,
$\partial\mathbf{e}_{j}^{\prime}\bm{\varepsilon}_{t}/\partial\bm{\theta}_{i}=\mathbf{0}$
is the $2n(n+1)$-dimensional null vector. For $i=j$, one obtains, component by
component,

$\frac{\partial\mathbf{e}_{j}^{\prime}\bm{\varepsilon}_{t}}{\partial\mu_{0j}}=1$

so $\partial\bm{\varepsilon}_{t}/\partial\mu_{0j}=\mathbf{e}_{j}$ and
$\partial\bm{\varepsilon}_{t}^{\prime}/\partial\bm{\mu}_{0}=\mathbf{I}_{n}$.
Similarly, it holds that
$\partial\bm{\varepsilon}_{t}^{\prime}/\partial\bm{\mu}_{1}=\mathbf{G}_{t}$.
Further,
$\partial\bm{\varepsilon}_{t}/\partial\phi_{0ij}=-\partial\mathbf{\Phi}_{0}\mathbf{y}_{t-1}/\partial\phi_{0ij}=-\mathbf{e}_{i}\mathbf{e}_{j}^{\prime}\mathbf{y}_{t-1}$,
and with a slight abuse of notation it follows that

$\displaystyle\frac{\partial\bm{\varepsilon}_{t}^{\prime}}{\partial\mathbf{\Phi}_{0}}$
$\displaystyle=-\begin{bmatrix}\mathbf{e}_{1}\mathbf{e}_{1}^{\prime}&\dots&\mathbf{e}_{1}\mathbf{e}_{n}^{\prime}\\\
\vdots&\vdots&\vdots\\\
\mathbf{e}_{n}\mathbf{e}_{1}^{\prime}&\dots&\mathbf{e}_{n}\mathbf{e}_{n}^{\prime}\end{bmatrix}\otimes\mathbf{y}_{t-1}$
$\displaystyle=-\begin{bmatrix}y_{1,t-1}\mathbf{e}_{1}\mathbf{e}_{1}^{\prime}&\dots&y_{1,t-1}\mathbf{e}_{1}\mathbf{e}_{n}^{\prime}\\\
\vdots&\vdots&\vdots\\\
y_{1,t-1}\mathbf{e}_{n}\mathbf{e}_{1}^{\prime}&\dots&y_{1,t-1}\mathbf{e}_{n}\mathbf{e}_{n}^{\prime}\\\
\vdots&\vdots&\vdots\\\
y_{n,t-1}\mathbf{e}_{1}\mathbf{e}_{1}^{\prime}&\dots&y_{1,t-1}\mathbf{e}_{1}\mathbf{e}_{n}^{\prime}\\\
\vdots&\vdots&\vdots\\\
y_{n,t-1}\mathbf{e}_{n}\mathbf{e}_{1}^{\prime}&\dots&y_{n,t-1}\mathbf{e}_{n}\mathbf{e}_{n}^{\prime}\end{bmatrix}$
(29)

and

$\displaystyle\frac{\partial\bm{\varepsilon}_{t}^{\prime}}{\partial\mathbf{\Phi}_{1}}$
$\displaystyle=-\begin{bmatrix}\mathbf{G}_{t}\mathbf{e}_{1}\mathbf{e}_{1}^{\prime}&\dots&\mathbf{G}_{t}\mathbf{e}_{1}\mathbf{e}_{n}^{\prime}\\\
\vdots&\vdots&\vdots\\\
\mathbf{G}_{t}\mathbf{e}_{n}\mathbf{e}_{1}^{\prime}&\dots&\mathbf{G}_{t}\mathbf{e}_{n}\mathbf{e}_{n}^{\prime}\end{bmatrix}\otimes\mathbf{y}_{t-1}$
$\displaystyle=-\begin{bmatrix}g_{1t}y_{1,t-1}\mathbf{e}_{1}\mathbf{e}_{1}^{\prime}&\dots&g_{1t}y_{1,t-1}\mathbf{e}_{1}\mathbf{e}_{n}^{\prime}\\\
\vdots&\vdots&\vdots\\\
g_{nt}y_{1,t-1}\mathbf{e}_{n}\mathbf{e}_{1}^{\prime}&\dots&g_{nt}y_{1,t-1}\mathbf{e}_{n}\mathbf{e}_{n}^{\prime}\\\
\vdots&\vdots&\vdots\\\
g_{1t}y_{n,t-1}\mathbf{e}_{1}\mathbf{e}_{1}^{\prime}&\dots&g_{1t}y_{1,t-1}\mathbf{e}_{1}\mathbf{e}_{n}^{\prime}\\\
\vdots&\vdots&\vdots\\\
g_{nt}y_{n,t-1}\mathbf{e}_{n}\mathbf{e}_{1}^{\prime}&\dots&g_{nt}y_{n,t-1}\mathbf{e}_{n}\mathbf{e}_{n}^{\prime}\end{bmatrix}.$
(30)

The matrices $\mathbf{e}_{1}\mathbf{e}_{1}^{\prime}$ have only one nonzero
element that equals one, so for the score, the corresponding $n^{2}\times n$
matrix becomes $\mathbf{y}_{t-1}\otimes\mathbf{I}_{n}$. Terms for the average
score from (29) are as a vector

$\mathbf{s}_{T}(\mathbf{\Phi}_{0})=\frac{1}{T}\sum_{t=1}^{T}(\mathbf{y}_{t-1}\otimes\mathbf{I}_{n})$

which is an $n^{2}\times 2$ matrix, and from (30) it follows that

$\mathbf{s}_{T}(\mathbf{\Phi}_{1})=\frac{1}{T}\sum_{t=1}^{T}(\mathbf{G}_{t}\mathbf{y}_{t-1}\otimes\mathbf{I}_{n}).$

Finally,

$\displaystyle\frac{\partial\mathbf{G}_{t}}{\partial\gamma_{i}}$
$\displaystyle=\frac{\partial\
\text{diag}(g_{1t},\ldots,g_{nt})}{\gamma_{i}}=\text{diag}\left(0,\ldots,\frac{\partial
g_{it}}{\partial\gamma_{i}},\ldots,0\right)$ $\displaystyle=\frac{\partial
g_{it}}{\partial\gamma_{i}}\text{diag}(0,\ldots,1,\ldots)=\frac{\partial
g_{it}}{\gamma_{i}}\mathbf{e}_{i}\mathbf{e}_{i}^{\prime}$

where $\partial g_{it}/\partial\gamma_{i}=g_{it}(1-g_{it})(s_{i,t}-c_{i})$.
This implies

$\frac{\partial\bm{\varepsilon}_{t}}{\partial\gamma_{i}}=-\frac{\partial
g_{it}}{\partial\gamma_{i}}\mathbf{e}_{i}\mathbf{e}_{i}^{\prime}(\bm{\mu}_{1}+\mathbf{\Phi}_{1}\mathbf{y}_{t-1}).$

For the $i$-th error, the derivative w.r.t. $\gamma_{i}$ equals

$\frac{\partial\varepsilon_{it}}{\partial\gamma_{i}}=-\frac{\partial
g_{it}}{\partial\gamma_{i}}(\bm{\mu}_{1}+\mathbf{\Phi}_{1}\mathbf{y}_{t-1}).$

Denoting $\mathbf{g}_{\bm{\gamma},t}=\left(\partial
g_{1t}/\partial\gamma_{1},\ldots,\partial
g_{nt}/\partial\gamma_{n}\right)^{\prime}$ and
$\bm{\gamma}=\left(\gamma_{1},\ldots,\gamma_{n}\right)^{\prime}$, one can
write, after removing the zero elements and using the notation
$\partial\bm{\varepsilon}_{t}/\partial\bm{\gamma}=\left(\partial\varepsilon_{1t}/\partial\gamma_{1},\ldots,\partial\varepsilon_{nt}/\partial\gamma_{n}\right)^{\prime}$,

$\frac{\partial\bm{\varepsilon}_{t}^{\prime}}{\partial\bm{\gamma}}=-\mathbf{g}_{\bm{\gamma},t}^{\prime}\otimes(\bm{\mu}_{1}+\mathbf{\Phi}_{1}\mathbf{y}_{t-1}),$

which is an $n\times n$ matrix. Now, $\partial g_{1t}/\partial
c_{i}=-\gamma_{i}g_{it}(1-g_{it})$. Analogously to the previous notation,
setting $\mathbf{c}=\left(c_{1},\ldots,c_{n}\right)^{\prime}$ and
$\mathbf{g}_{\mathbf{c},t}=\left(\partial g_{1t}/\partial
c_{1},\ldots,\partial g_{nt}/\partial c_{n}\right)^{\prime}$, one obtains

$\frac{\partial\bm{\varepsilon}_{t}^{\prime}}{\partial\mathbf{c}}=-\mathbf{g}_{\mathbf{c},t}\otimes(\bm{\mu}_{1}+\mathbf{\Phi}_{1}\mathbf{y}_{t-1}).$

Drawing things together, the score vector equals (assuming $\mathbf{y}_{0}$
known)

$\begin{bmatrix}\mathbf{s}_{T}(\mathbf{\Phi}_{0})\\\
\mathbf{s}_{T}(\mathbf{\Phi}_{1})\\\ \mathbf{s}_{T}(\bm{\gamma})\\\
\mathbf{s}_{T}(\mathbf{c})\end{bmatrix}=-\frac{1}{T}\sum_{t=1}^{T}\begin{bmatrix}\left(\mathbf{I}_{n}\otimes\mathbf{y}_{t-1}\right)\mathbf{\Omega}_{t}^{-1}\bm{\varepsilon}_{t}\\\
\left(\mathbf{I}_{n}\otimes\mathbf{G}_{t}\mathbf{y}_{t-1}\right)\mathbf{\Omega}_{t}^{-1}\bm{\varepsilon}_{t}\\\
\left\\{\mathbf{g}_{\bm{\gamma},t}^{\prime}\otimes(\bm{\mu}_{1}+\mathbf{\Phi}_{1}\mathbf{y}_{t-1})\right\\}\mathbf{\Omega}^{-1}_{t}\bm{\varepsilon}_{t}\\\
\left\\{\mathbf{g}_{\mathbf{c},t}^{\prime}\otimes(\bm{\mu}_{1}+\mathbf{\Phi}_{1}\mathbf{y}_{t-1})\right\\}\mathbf{\Omega}^{-1}_{t}\bm{\varepsilon}_{t}\end{bmatrix}.$
(31)

When the model is the one in Eq. (27), the purpose could be testing one
transition against two. In this case, the term
$\mathbf{G}_{t}^{(2)}\left(\bm{\mu}_{2}+\mathbf{\Phi}_{2}\mathbf{y}_{t-1}\right)$
can be approximated through a third-order Taylor expansion, so that the error
vector becomes

$\bm{\varepsilon}_{t}=\mathbf{y}_{t}-\bm{\mu}_{0}-\mathbf{\Phi}_{0}\mathbf{y}_{t-1}-\mathbf{G}_{t}^{(1)}-\sum_{l=1}^{L=3}\mathbf{\Upsilon}_{l}(\mathbf{y}_{t-1}\odot\mathbf{s}_{t}^{l})$

where
$\mathbf{\Upsilon}_{l}=\text{diag}\left(\upsilon_{l1},\ldots,\upsilon_{ln}\right)$
and $\mathbf{s}_{t}^{l}=\left(s_{1t}^{l},\ldots,s_{nt}^{l}\right)^{\prime}$,
$l=1,2,3$. Then,
$\partial\mathbf{\Upsilon}_{l}/\partial\upsilon_{lj}=\mathbf{e}_{j}\mathbf{e}_{j}^{\prime}$
and
$\partial\bm{\varepsilon}_{t}/\partial\upsilon_{lj}=\mathbf{e}_{j}^{\prime}(\mathbf{y}_{t-1}\odot\mathbf{s}_{t}^{i})$.
Finally, defining

$\frac{\partial\bm{\varepsilon}_{t}}{\partial\bm{\upsilon}_{l}}=\left(\frac{\partial\bm{\varepsilon}_{t}}{\partial\upsilon_{l1}},\ldots,\frac{\partial\bm{\varepsilon}_{t}}{\partial\upsilon_{ln}}\right)^{\prime},\quad\text{for
}l=1,2,3$

the corresponding components for the average score are

$\mathbf{s}_{T}(\bm{\upsilon}_{l})=\frac{1}{T}\sum_{t=1}^{T}\left\\{\left(\mathbf{y}_{t-1}\odot\mathbf{s}_{t}^{l}\right)\otimes\mathbf{I}_{n}\right\\}\mathbf{\Omega}_{t}^{-1}\bm{\varepsilon}_{t}.$

Equivalently, the average score for the special case of testing linearity
against a single transition can be obtained by suppressing the redundant rows
from (31).

##  Appendix B Proof of Theorem 1

###### Proof.

Assuming a first-order Taylor expansion, the Lagrange multiplier test under
the null is derived from the score matrix

$\frac{\partial\ell_{T}(\hat{\bm{\theta}})}{\partial\mathbf{D}_{1}}=\sum_{t=1}^{T}\left\\{\mathbf{x}_{t}\mathbf{s}_{t}\left(\mathbf{y}_{t}-\mathbf{\hat{D}}_{0}^{\prime}\mathbf{x}_{t}\right)^{\prime}\hat{\mathbf{\Omega}}^{-1}\right\\}.$
(32)

Setting $\mathbf{Y}=\left[\mathbf{y}_{1}^{\prime}\ \mathbf{y}_{2}^{\prime}\
\ldots\ \mathbf{y}_{T}^{\prime}\right]^{\prime}$,
$\mathbf{Z}=\left[x_{1}^{\prime}s_{1}\ x_{2}^{\prime}s_{2}\ \ldots\
x_{T}^{\prime}s_{T}\right]$, we have that

$\frac{\partial\ell_{T}(\hat{\bm{\theta}})}{\partial\mathbf{D}_{1}}=\mathbf{Z}^{\prime}\left(\mathbf{Y}-\mathbf{X}\mathbf{\hat{D}}_{0}\right)\hat{\mathbf{\Omega}}^{-1},$
(33)

where $\hat{\mathbf{D}}_{0}$ and $\hat{\mathbf{\Omega}}$ are estimates under
the null hypothesis. Under regularity conditions, the score converges in
probability to a matrix-normal distribution with zero mean and variance
$\mathbf{Z}^{\prime}(\mathbf{I}-\mathbf{P})Z\otimes\mathbf{\Omega}^{-1}$
conditional on $\mathbf{X}$ and $\mathbf{Z}$, where
$\mathbf{P}_{X}\equiv\mathbf{X}(\mathbf{X}^{\prime}\mathbf{X})^{-1}\mathbf{X}$
is the projection matrix of $\mathbf{X}$.

To see this, we write (33) as follows

$\displaystyle\mathbf{Q}$
$\displaystyle\equiv\frac{\partial\ell_{T}(\hat{\bm{\theta}})}{\partial\mathbf{D}_{1}}=\mathbf{Z}^{\prime}\left(\mathbf{Y}-\mathbf{X}\hat{\mathbf{D}}_{0}\right)\hat{\mathbf{\Omega}}^{-1}$
$\displaystyle=\mathbf{Z}^{\prime}\left(\mathbf{Y}-\mathbf{X}(\mathbf{X}^{\prime}\mathbf{X})\mathbf{X}^{\prime}\mathbf{Y}\right)\hat{\mathbf{\Omega}}^{-1}$
$\displaystyle=\mathbf{Z}^{\prime}\left(\mathbf{I}-\mathbf{P}_{X}\right)\left(\mathbf{X}\mathbf{D}_{0}+\mathbf{E}\right)\hat{\mathbf{\Omega}}^{-1}$
$\displaystyle=\mathbf{Z}^{\prime}\left(\mathbf{I}-\mathbf{P}_{X}\right)\mathbf{E}\hat{\mathbf{\Omega}}^{-1}.$

Under the null hypothesis, $\mathbf{Y}=\mathbf{X}\mathbf{D}_{0}+\mathbf{E}$,
where $\mathbf{E}=\left[\bm{\varepsilon}_{1}^{\prime}\ \ldots\
\bm{\varepsilon}_{T}^{\prime}\right]^{\prime}$ and
$\text{vec}(\mathbf{E}^{\prime})$ follows a
$\mathcal{N}\left(\mathbf{0},\mathbf{I}_{T}\otimes\mathbf{\Omega}\right)$
distribution. Under the null hypothesis, $\hat{\mathbf{\Omega}}$ will converge
in probability to $\mathbf{\Omega}$. Set

$\mathbf{S}=\left(\mathbf{Z}^{\prime}(\mathbf{I}-\mathbf{P}_{X})\mathbf{Z}\right)^{-1/2}\mathbf{Q}\hat{\mathbf{\Omega}}^{1/2}$

which will asymptotically converge to a zero-mean matrix-normal distribution
with variance $\mathbf{I}\otimes\mathbf{I}$. Thus, we have the chi-square
version LM test statistic

$LM=\text{tr}\left\\{\mathbf{S}^{\prime}\mathbf{S}\right\\}=\text{tr}\left\\{\hat{\mathbf{\Omega}}^{-1}(\mathbf{Y}-\mathbf{X}\hat{\mathbf{D}})^{\prime}\mathbf{Z}[\mathbf{Z}^{\prime}(\mathbf{I}_{T}-\mathbf{P}_{X})\mathbf{Z}]^{-1}\mathbf{Z}^{\prime}(\mathbf{Y}-\mathbf{X}\hat{\mathbf{D}}_{0})\right\\}$

which converges to the $\chi^{2}_{n(1+np)}$-distribution when the null
hypothesis is valid. ∎

##  Appendix C First-order partial derivatives of
$\mathbf{\Psi}_{t}^{\prime}\mathbf{B}^{\prime}\mathbf{x}_{t}$

The vectorised first-order derivative of
$\mathbf{\Psi}_{t}^{\prime}\mathbf{B}^{\prime}\mathbf{x}_{t}$ with respect to
parameters $\bm{\theta}$ can be easily found in both univariate and
multivariate cases, see Eitrheim and Teräsvirta, (1996) and Teräsvirta and
Yang, 2014a . The set of parameters $\bm{\theta}$ consists of
$\text{vec}(\mathbf{B})$, $\bm{\gamma}$ and $\mathbf{C}$, where
$\mathbf{B}=\left[b_{ij}\right]$, $\bm{\gamma}=\left[\gamma_{ij}\right]$ and
$\mathbf{C}=\left[c_{ij}\right]$.

For the $ij$-th parameter of $\mathbf{B}$, $b_{ij}$, we have

$\frac{\partial\mathbf{\Psi}_{t}^{\prime}\mathbf{B}^{\prime}\mathbf{x}_{t}}{\partial
b_{ij}}=\mathbf{\Psi}_{t}^{\prime}\mathbf{H}_{ij}\mathbf{x}_{t}$ (34)

where $\mathbf{H}_{ij}=\left[h_{kl}\right]$ is a matrix in which $h_{ij}=1$
and $h_{kl}=0$ for $k\neq i$ and $l\neq i$. The first derivative in (34) is
the directional derivative of
$\mathbf{\Psi}_{t}^{\prime}\mathbf{B}^{\prime}\mathbf{x}_{t}$ with respect to
the unit length matrix $\mathbf{H}_{ij}$.

For the parameter matrices $\bm{\gamma}=\left[\gamma_{ij}\right]$ and
$\mathbf{C}=\left[c_{ij}\right]$, letting $\beta_{ij}=\gamma_{ij},c_{ij}$, we
have

$\frac{\partial\mathbf{\Psi}_{t}^{\prime}\mathbf{B}^{\prime}\mathbf{x}_{t}}{\partial\beta_{ij}}=\left(0_{n},\ldots,\frac{\partial\mathbf{G}_{t}^{(d)}}{\partial\beta_{ij}},\ldots,0_{n}\right)\mathbf{B}^{\prime}\mathbf{x}_{t}=\frac{\partial\mathbf{G}_{t}^{(d)}}{\partial\beta_{ij}}\mathbf{B}_{d+1}^{\prime}\mathbf{x}_{t}$
(35)

for $d=1,\ldots,m-1$, where

$\frac{\partial\mathbf{G}_{t}^{(d)}}{\beta_{ij}}=\text{diag}\left\\{0,\ldots,\frac{\partial
g_{t}^{ij}}{\beta_{ij}},\ldots,0\right\\}$ (36)

for $j=1,\ldots,n$. When $\beta_{ij}=\gamma_{ij}$,

$\frac{\partial
g^{ij}_{t}}{\gamma_{ij}}=\left(g_{t}^{ij}\right)^{2}\exp\left\\{-\gamma_{ij}\left(s_{t}-c_{ij}\right)\right\\}\left(s_{t}-c_{ij}\right)=\left(s_{t}-c_{ij}\right)g_{t}^{ij}\left(1-g_{t}^{ij}\right)$
(37)

and when $\beta_{ij}=c_{ij}$,

$\frac{\partial g_{t}^{ij}}{\partial
c_{ij}}=-\left(g_{t}^{ij}\right)^{2}\text{exp}\left\\{-\gamma_{ij}\left(s_{t}-c_{ij}\right)\right\\}\gamma_{ij}=-\gamma_{ij}g_{t}^{ij}\left(1-g_{t}^{ij}\right).$
(38)

The dimension of the first-order derivative of
$\mathbf{\Psi}_{t}^{\prime}\mathbf{B}^{\prime}\mathbf{x}_{t}$ with respect to
$\bm{\theta}$ is $n\times\left[(pn+1)mn+2(m-1)n\right]$.

##  Appendix D Wilks’ extension of the test

As in Teräsvirta and Yang, 2014a , we use an improvement to the rescaled LM-
type test presented in Eq. (23) based on the so-called Wilks’
$\Lambda$-distribution.

######  Theorem 3.

Let $\mathbf{\text{RSS}}_{0}$ and $\mathbf{\text{RSS}}_{1}$ be the $n\times n$
residual sum of squares from respectively the restricted and the auxiliary
regression, and let
$\mathbf{W}_{1}=\mathbf{\text{RSS}}_{0}-\mathbf{\text{RSS}}_{1}$ and
$\mathbf{W}_{2}=\mathbf{\text{RSS}}_{1}$. In the linearity test, under the
null, $\mathbf{W}_{1}$ and $\mathbf{W}_{2}$ are two independent Wishart-
distributed random matrices

$\mathbf{W}_{1}\sim\mathcal{W}_{n}\left(\mathbf{\Omega},\text{cd}(\mathbf{Z})\right),\qquad\mathbf{W}_{2}\sim\mathcal{W}_{n}\left(\mathbf{\Omega},T-\text{cd}(\mathbf{X})-\text{cd}(\mathbf{Z})\right)$
(39)

where $\text{cd}(\cdot)$ is the column dimension of a matrix.

###### Proof.

The score matrix evaluated under the null hypothesis has the general form

$\frac{\partial\ell_{T}(\hat{\bm{\theta}})}{\partial\mathbf{D}_{1}}=\mathbf{Z}^{\prime}\left(\mathbf{Y}-\mathbf{X}\hat{\mathbf{D}}_{0}\right)\hat{\mathbf{\Omega}}^{-1}.$

Using the auxiliary approach for computing the test statistic produces two
residual sum of squares, $\mathbf{\text{RSS}}_{0}$ and
$\mathbf{\text{RSS}}_{1}$. The former is the residual sum of the squares
matrix from the restricted regression, i.e.,
$\mathbf{\text{RSS}}_{0}=\hat{\mathbf{E}}^{\prime}\hat{\mathbf{E}}$,
$\hat{\mathbf{E}}=(\mathbf{I}-\mathbf{P}_{X})\mathbf{Y}$, where
$\mathbf{P}_{X}$ is the projection matrix of $\mathbf{X}$. Notice that under
the null hypothesis, $\mathbf{Y}=\mathbf{X}\mathbf{D}_{0}+\mathbf{E}$, where
$\mathbf{E}=\left[\bm{\varepsilon}_{1}^{\prime}\ \ldots\
\bm{\varepsilon}_{T}^{\prime}\right]^{\prime}$, and
$\text{vec}(\mathbf{E}^{\prime})\sim\mathcal{N}(\mathbf{0},\mathbf{I}_{T}\otimes\mathbf{\Omega})$.
Consequently, we have that
$\hat{\mathbf{E}}=(\mathbf{I}-\mathbf{P}_{X})\mathbf{Y}=(\mathbf{I}-\mathbf{P}_{X})\mathbf{E}$.
$\mathbf{RSS}_{1}$ is the residual sum of the squares matrix from the
auxiliary regression, i.e.,
$\mathbf{RSS}_{1}=\hat{\mathbf{\Xi}}^{\prime}\hat{\mathbf{\Xi}}$, where
$\hat{\mathbf{\Xi}}=(\mathbf{I}-\mathbf{P}_{XZ})\hat{\mathbf{E}}$, where
$\mathbf{P}_{XZ}$ is the projection matrix of the matrix
$\left[\mathbf{X},\mathbf{Z}\right]$, such that

$\mathbf{P}_{XZ}=\left[\mathbf{X}\
\mathbf{Z}\right]=\begin{bmatrix}\mathbf{X^{\prime}X}&\mathbf{X^{\prime}Z}\\\
\mathbf{Z^{\prime}X}&\mathbf{Z^{\prime}Z}\end{bmatrix}^{-1}\begin{bmatrix}\mathbf{X^{\prime}}\\\
\mathbf{Z^{\prime}}\end{bmatrix}.$

Therefore, $\mathbf{W}_{1}$ can be reformulated as follows

$\displaystyle\mathbf{W}_{1}$
$\displaystyle=\mathbf{\text{RSS}}_{0}-\mathbf{\text{RSS}}_{1}=\hat{\mathbf{E}}^{\prime}\hat{\mathbf{E}}-\hat{\mathbf{\Xi}}^{\prime}\hat{\mathbf{\Xi}}$
$\displaystyle=\hat{\mathbf{E}}^{\prime}\mathbf{P}_{XZ}\hat{\mathbf{E}}=\hat{\mathbf{E}}^{\prime}\mathbf{Z}(\mathbf{Z}^{\prime}(\mathbf{I}_{T}-\mathbf{P}_{X})\mathbf{Z})^{-1}\mathbf{Z}^{\prime}\hat{\mathbf{E}}$
$\displaystyle=\mathbf{E}^{\prime}(\mathbf{I}_{T}-\mathbf{P}_{X})\mathbf{Z}(\mathbf{Z}^{\prime}(\mathbf{I}_{T}-\mathbf{P}_{X})\mathbf{Z})^{-1}\mathbf{Z}^{\prime}(\mathbf{I}_{T}-\mathbf{P}_{X})\mathbf{E}.$

Let $\mathbf{I}_{T}-\mathbf{P}_{X}=\mathbf{R}^{\prime}\mathbf{R}$, where
$\mathbf{R}$ is orthogonal to $\mathbf{X}$ and
$\mathbf{R}^{\prime}\mathbf{R}=\mathbf{I}_{T-\text{cd}(X)}$, then

$\mathbf{W}_{1}=\mathbf{E}^{\prime}\mathbf{R}^{\prime}\mathbf{R}\mathbf{Z}(\mathbf{Z}^{\prime}\mathbf{R}\mathbf{R}^{\prime}\mathbf{Z})^{-1}\mathbf{Z}^{\prime}\mathbf{R}\mathbf{R}^{\prime}\mathbf{E}.$

Setting
$\mathbf{V}_{1}=\mathbf{Z}^{\prime}\mathbf{R}\mathbf{R}^{\prime}\mathbf{E}$,
we have that
$\mathbf{V}_{1}\sim\mathcal{N}(\mathbf{0},\mathbf{Z^{\prime}RR^{\prime}Z}\otimes\mathbf{\Omega})$,
therefore $\mathbf{W}_{1}$ follows a Wishart distribution generated by
$\mathbf{V}_{1}$:

$\mathbf{W}_{1}=\mathbf{V}_{1}^{\prime}(\mathbf{Z^{\prime}RR^{\prime}Z})^{-1}\mathbf{V}_{1}\sim\mathcal{W}_{n}(\mathbf{\Omega},\text{cd}(\mathbf{Z})).$

For $\mathbf{W}_{2}$, we obtain

$\displaystyle\mathbf{W}_{2}$
$\displaystyle=\mathbf{RSS}_{1}=\hat{\mathbf{\Xi}}^{\prime}\hat{\mathbf{\Xi}}=\hat{\mathbf{E}}^{\prime}(\mathbf{I}-\mathbf{P}_{XZ})\hat{\mathbf{E}}=\hat{\mathbf{E}}^{\prime}\hat{\mathbf{E}}-\hat{\mathbf{E}}^{\prime}\mathbf{P}_{XZ}\hat{\mathbf{E}}$
$\displaystyle=\hat{\mathbf{E}}^{\prime}\hat{\mathbf{E}}-\hat{\mathbf{E}}^{\prime}\mathbf{Z}(\mathbf{Z}^{\prime}(\mathbf{I}-\mathbf{P}_{X})\mathbf{Z})^{-1}\mathbf{Z}^{\prime}\hat{\mathbf{E}}$
$\displaystyle=\mathbf{E}^{\prime}\mathbf{RR^{\prime}(I-Z(Z^{\prime}RR^{\prime}Z)^{-1}}\mathbf{Z^{\prime}})\mathbf{RR^{\prime}E}$
$\displaystyle=\mathbf{E^{\prime}R}(\mathbf{I}-\mathbf{R^{\prime}Z}(\mathbf{Z^{\prime}RR^{\prime}Z})^{-1}\mathbf{Z^{\prime}R})\mathbf{R^{\prime}E}.$

By imposing
$\mathbf{I}_{T-\text{cd}(X)}-\mathbf{R^{\prime}Z}(\mathbf{Z^{\prime}RR^{\prime}Z})^{-1}\mathbf{Z^{\prime}R}=\mathbf{QQ^{\prime}}$,
where $\mathbf{Q}\bot\mathbf{R^{\prime}Z}$ (with $\bot$ indicating
orthogonality) and
$\mathbf{QQ^{\prime}}=\mathbf{I}_{T-\text{cd}(X)-\text{cd}(Z)}$, we have that
$\mathbf{W}_{2}=\mathbf{E^{\prime}RQQ^{\prime}R^{\prime}E}$. Setting
$\mathbf{V}_{2}=\mathbf{Q^{\prime}R^{\prime}E}$, we have that
$\mathbf{V}_{2}\sim\mathcal{N}(\mathbf{0},\mathbf{I}\otimes\mathbf{\Omega})$
and that $\mathbf{W}_{2}$ follows a Wishart distribution generated by
$\mathbf{V}_{2}$, such that

$\mathbf{W}_{2}=\mathbf{V}_{2}^{\prime}\mathbf{V}_{2}\sim\mathcal{W}_{n}(\mathbf{\Omega},T-\text{cd}(X)-\text{cd}(Z)).$

Stacking the columns of $\mathbf{V}_{1}$ and $\mathbf{V}_{2}$ yields the
random matrix

$\mathbf{U}=\begin{pmatrix}\mathbf{V}_{1}\\\
\mathbf{V}_{2}\end{pmatrix}=\begin{pmatrix}\mathbf{Z^{\prime}R}\\\
\mathbf{Q^{\prime}}\end{pmatrix}\mathbf{R^{\prime}E}.$

It follows that
$\mathbf{U}\sim\mathcal{N}(\mathbf{0},\mathbf{\Sigma}\otimes\mathbf{\Omega})$,
where the row covariance matrix is

$\mathbf{\Sigma}=\begin{pmatrix}\mathbf{Z^{\prime}R}\\\
\mathbf{Q^{\prime}}\end{pmatrix}\mathbf{R^{\prime}R}(\mathbf{R^{\prime}ZQ})=\begin{pmatrix}\mathbf{Z^{\prime}RR^{\prime}Z}&\mathbf{Z^{\prime}RQ}\\\
\mathbf{Q^{\prime}R^{\prime}Z}&\mathbf{Q^{\prime}Q}\end{pmatrix}=\begin{pmatrix}\mathbf{Z^{\prime}RR^{\prime}Z}&\mathbf{0}\\\
\mathbf{0}&\mathbf{I}\end{pmatrix}$

given that $\mathbf{Q}\bot\mathbf{R^{\prime}Z}$. We can conclude that
$\mathbf{V}_{1}$ and $\mathbf{V}_{2}$ are uncorrelated and independent due to
normality, consequently $\mathbf{W}_{1}$ and $\mathbf{W}_{2}$ are independent
as desired. ∎

When $\mathbf{A}\sim\mathcal{W}_{n}(\mathbf{\Sigma},m)$ and
$\mathbf{B}\sim\mathcal{W}_{n}(\mathbf{\Sigma},d)$ are independent,
$\mathbf{\Sigma}$ is an $n\times n$ positive definite matrix, $m\geq n$, the
Wilks’ statistic is

$\Lambda=\frac{|\mathbf{A}|}{|\mathbf{A}+\mathbf{B}|}=|\mathbf{I}_{n}+\mathbf{A}^{-1}+\mathbf{B}|^{-1}\sim\mathcal{W}(n,m,d)$
(40)

and has a Wilks’ $\Lambda$-distribution with parameters $n$, $m$ and $d$, see
Mardia et al., (1979) and Anderson, (2003) for a discussion on the Wilks’
$\Lambda$ distribution. Setting $\mathbf{A}=\mathbf{W}_{2}$ and
$\mathbf{B}=\mathbf{W}_{2}$, the test statistic can be written as

$\Lambda=\frac{|\mathbf{W}_{2}|}{|\mathbf{W}_{2}+\mathbf{W}_{1}|}=\frac{|\mathbf{RSS}_{1}|}{|\mathbf{RSS}_{0}|}$
(41)

which, under linearity, follows a Wilks’ $\Lambda$-distribution
$\mathcal{L}(n,T-\text{cd}(\mathbf{X})-\text{cd}(\mathbf{Z}),\text{cd}(\mathbf{Z}))$.
If $T$ is large, the Bartlett’s approximation can be used

$\lambda=\left(\frac{1}{2}(n+\text{cd}(\mathbf{Z})+1)+\text{cd}(\mathbf{X})-T\right)\log\Lambda\sim\chi^{2}_{\text{cd}(\mathbf{Z})n},$

see Bartlett, (1954) and Anderson, (2003).

The test statistic can be carried out after performing steps 1 and 2 in the
algorithm in Section 3. Since matrix $\mathbf{Z}$ can be set to any
$\mathbf{Z}$ in the sequential test procedure in Section 4, the improved test
statistic can be computed for both the tests used in this paper for the
sequential procedure.