A Sufficient Statistical Test for Dynamic Stability

Muhammad Ashfaq Ahmed, Nasreen Nawaz

Federal Board of Revenue, Islamabad, Pakistan.

DOI: 10.4236/jdaip.2023.112009                

Abstract

In the existing Statistics and Econometrics literature, there does not exist a statistical test which may test for all kinds of roots of the characteristic polynomial leading to an unstable dynamic response, i.e., positive and negative real unit roots, complex unit roots and the roots lying inside the unit circle. This paper develops a test which is sufficient to prove dynamic stability (in the context of roots of the characteristic polynomial) of a univariate as well as a multivariate time series without having a structural break. It covers all roots (positive and negative real unit roots, complex unit roots and the roots inside the unit circle whether single or multiple) which may lead to an unstable dynamic response. Furthermore, it also indicates the number of roots causing instability in the time series. The test is much simpler in its application as compared to the existing tests as the series is strictly stationary under the null (C01, C12).

Keywords

Dynamic Stability, Real and Complex Roots, Unit Circle

1. Introduction

A univariate time series that can be written in the form

$y_t=\mu +{\epsilon }_t+{\psi }_1{\epsilon }_{t-1}+{\psi }_2{\epsilon }_{t-2}+\cdots =\mu +\psi \left(L\right){\epsilon }_t,$

with ${\displaystyle {\sum }_{j=0}^{\infty }\left|{\psi }_j\right|}<\infty$ , roots of $\psi \left(\mathbb{Z}\right)=0$ outside the unit circle, and $\left\{{\epsilon }_t\right\}$ a white noise process with mean zero and variance ${\sigma }^2$ is a covariance stationary time series with 1) $E\left(y_t\right)=\mu$ and 2) $cov\left(y_t,y_{t-j}\right)={\gamma }_j$ , for all t where j is an integer. A process is said to be strictly stationary if the joint distribution of $\left(Y_t,Y_{t+j_1},Y_{t+j_2},\cdots ,Y_{t+j_n}\right)$ depends only on the intervals separating the dates, i.e., $j_1,j_2,\cdots ,j_n$ and not on the date itself, i.e., t. The most important feature of a stable dynamic response which distinguishes it from an unstable response is its convergence property, i.e., the future forecasts ${\hat{y}}_{t+s|t}\equiv \hat{E}\left(y_{t+s}|y_t,y_{t-1},\cdots \right)$ converge to the unconditional mean, i.e.,

$\underset{s\to \infty }{\lim }{\hat{y}}_{t+s|t}=\mu \text{.}$

Thus if someone is trying to forecast a series farther into the future, it becomes very important to know whether the series is dynamically stable or not. The most common approach in time series literature is testing for a real positive unit root. A number of tests have been proposed, including [1] - [27] and [28] .

The use of CUSUM for the stability (or stationarity) of regression equations was introduced by [29] . The CUSUM test which is based on the residuals from the recursive estimates is normally used to test the parameter change in time series model. A number of tests for testing the structural change are available in the existing literature, such as [30] [31] [32] , etc. There has been some work related to complex unit roots as well, e.g., [33] , etc.

However, if the null hypothesis of a unit root is rejected, the time series can still be unstable due to presence of other kinds of roots leading to instability as unit root is not the sole root as a cause of concern regarding stability. If the null hypothesis of unit root is rejected, forecast/prediction of time series can still be seriously flawed unless there is a test available to guarantee that the series is dynamically stable.

In the existing literature no test has been proposed which may test for all kinds of roots leading to an unstable response, i.e., real as well as complex unit roots along with the roots (of $\psi \left(\mathbb{Z}\right)=0$ ) which lie inside the unit circle.

This paper develops a test which is sufficient to prove the dynamic stability (in the context of roots of the characteristic polynomial) of a univariate as well as a multivariate time series without having a structural break. It covers all roots (positive and negative real unit roots, complex unit roots and roots inside the unit circle whether single or multiple) which may lead to an unstable dynamic response. Furthermore, it also indicates the number of roots causing instability in the time series. The test is much simpler in its application as compared to the existing tests as the series is strictly stationary under the null.

The remainder of this paper is organized as follows: Section 2 provides the background of the test explaining the Routh’s stability criterion, formation of the Routh array, the theorems of the Routh test and the bilinear transformation. Section 3 discusses the hypothesis testing. Section 4 explains the methodology regarding the constrained minimization with inequality constraints. Section 5 consists of the theorems regarding the distributions under the null for various kinds of stability tests. Section 6 provides the power and size performance of the test. Section 7 consists of a Monte Carlo study. Section 8 provides an empirical application of the test. Section 9 comprises of the conclusion and finally the proofs of the theorems and the derivation of the null hypothesis for a VAR is provided in the appendix.

2. Background

The main idea behind the test is to expolit the Routh-Hurwitz stability criterion which is a mathematical test that provides a necessary and sufficient condition for the stability of a linear time invariant (LTI) control system. The Routh test was proposed by an English mathematician Edward John Routh in 1876 which is an efficient recursive algorithm to determine whether all the roots of the characteristic polynomial of a linear system have negative real parts.

German mathematician Adolf Hurwitz arranged the coefficients of the polynomial into a square matrix in 1895, and showed that the polynomial is stable if and only if the sequence of determinants of its principal submatrices is all positive. These two procedures are exactly equivalent. Routh stability criterion provides a more efficient way to compute the Hurwitz determinants than computing them directly.

For discrete systems, the corresponding stability test can be handled through the bilinear transformation, the Jury test or the Bistritz test which are all equivalent, however the bilinear transformation is much simpler in its use.

[34] compares graphically, using the Arnold Tongues, some sufficient criteria for the stability of periodic differential equations.

For using these stability criteria, the parameters need to be estimated for which the procedure has been described in the later part.

2.1. Routh’s Stability Criterion

The Routh test is a purely algebraic method for determining how many roots of the characteristic equation have positive real parts; from this it can also be determined whether the system is stable, for if there are no roots with positive real parts, the system is stable in continuous time framework. The algorithm for applying Routh’s stability criterion requires the order of the polynomial (the characteristic equation) to be finite and is as follows:

Write the characteristic equation in the form

$a_0{w}^n+a_1{w}^{n-1}+a_2{w}^{n-2}+\cdots +a_n=0,$ (1)

where $a_0$ is positive (if $a_0$ is originally negative, both sides are multiplied by −1). In this form, it is necessary that all the coefficients

$a_0,a_1,a_2,\cdots ,a_{n-1},a_n$

be positive if all the roots are to lie in the left half plane. If any coefficient is negative, the system is definitely unstable, and the Routh test is not needed to answer the question of stability. However, in this case, the Routh test will tell us the number of roots in the right half plane. If all the coefficients are positive, the system may be stable or unstable. It is then necessary to apply the following procedure to determine stability.

2.1.1. Routh Array

Arrange the coefficients of Equation (1) into the first two rows of the Routh array (Table 1) as follows.

The array has been filled in for $n=7$ in order to simplify the discussion. For any other value of n, the array is prepared in the same manner. In general, there

Table 1. Routh array.

are $\left(n+1\right)$ rows. For n even, the first row has one more element than the second row. The elements in the remaining rows are found from the formulas

$\begin{array}{cccc}{b}_1=\frac{a_1{a}_2-a_0{a}_3}{a_1}& b_2=\frac{a_1{a}_4-a_0{a}_5}{a_1}& b_3=\frac{a_1{a}_6-a_0{a}_7}{a_1}& \cdots \\ c_1=\frac{b_1{a}_3-a_1{b}_2}{b_1}& c_2=\frac{b_1{a}_5-a_1{b}_3}{b_1}& c_3=\frac{b_1{a}_7-a_1{b}_4}{b_1}& \cdots \\ d_1=\frac{c_1{b}_2-b_1{c}_2}{c_1}& d_2=\frac{c_1{b}_3-b_1{c}_3}{c_1}& \cdots & \cdots \\ \cdots & \cdots & \cdots & \cdots \end{array}$

The elements for the other rows are found from formulas that correspond to those just given. The elements in any row are always derived from the elements of the two preceding rows. During the computation of the Routh array, any row can be divided by a positive constant without changing the results of the test. (The application of this rule often simplifies the arithmetic.)

Having obtained the Routh array, the following theorems are applied to determine stability.

2.1.2. Theorems of the Routh Test

Theorem 1 The necessary and sufficient condition for all the roots of the characteristic equation to have negative real parts (stable system) is that all elements of the first column of the Routh array be positive and non-zero.

Theorem 2 If some of the elements in the first column are negative, the number of roots with a positive real part (in the right half plane) is equal to the number of sign changes in the first column.

Theorem 3 If one pair of roots is on the imaginary axis, equidistant from the origin, and all other roots are in the left half plane, all the elements of the nth row will vanish and none of the elements of the preceding row will vanish. The location of the pair of imaginary roots can be found by solving the equation

$C{w}^2+D=0,$ (2)

where the coefficients C and D are the elements of the array in the (n − 1)th row as read from left to right, respectively.

The algebraic method for determining stability is limited in its usefulness in that all we can learn from it is whether a system is stable. It does not give us any idea of the degree of stability or the roots of the characteristic equation.

Example: Given the characteristic equation

$w^4+3w^3+5w^2+4w+2=0,$

let’s determine the stability by the Routh criterion as follows:

Since all the coefficients are positive, the system may be stable. To test this, form the following Routh array (Table 2).

Since there is no change in sign in the first column, there are no roots having positive real parts, and the system is stable.

2.2. Bilinear Transformation

The Routh test which is often used to examine the roots of the characteristic equation of a continuous system may also be used to examine the roots of the characteristic equation of a discrete data system. The Routh test detects the presence of roots in the right half of w-plane. Since the criterion of stability of a discrete data system requires that all roots fall within the unit circle of the z-plane (or outside the unit circle of the $\mathbb{Z}=L=z^{-1}$ plane), one must first apply a transformation that will map the inside of the unit circle of the z-plane into the left half of the w-plane. One can then apply the Routh test to discover roots in the right half of the w-plane, and if none are found, we know that the roots of the characteristic equation fall within the unit circle and that the discrete data system is stable.

A transformation that will map the inside of the unit circle of the z-plane into the left half of the w-plane is

$z=\frac{w+1}{w-1}\text{.}$ (3)

This transformation is called the bilinear-transformation. The regions involving the transformation are shown in Figure 1.

3. Hypothesis Testing

Let us consider the following AR(1) process:

$\left(1-{\varphi }_{1}L\right)y_t=c+{ϵ}_t,$ (4)

where ${ϵ}_t\sim N\left(0,1\right)\left[iid\right]$ .

The characteristic polynomial of the above expression is as follows: $1-{\varphi }_{1}L=0$ which can be written in the z-plane as: $1-{\varphi }_1{z}^{-1}=0$ , which implies that

$z-{\varphi }_1=0.$

Applying the bilinear transformation on the above expression gives:

$\frac{w+1}{w-1}-{\varphi }_1=0.$

Table 2. Example of routh array.

Figure 1. Bilinear transformation.

After rearranging, we get:

$\left(1-{\varphi }_1\right)w+\left(1+{\varphi }_1\right)=0.$ (5)

The Routh Array for the above expression is as follows:

The null hypothesis of stability of above AR(1) process is as follows:

$H_0:\left[\begin{array}{l}1-{\varphi }_1\hfill \\ 1+{\varphi }_1\hfill \end{array}\right]>0$ , or $H_0:-1<{\varphi }_1<1$ .

As another example, suppose there is an AR(2) process:

$\left(1-{\varphi }_{1}L-{\varphi }_2{L}^2\right)y_t=c+{ϵ}_t,$ (6)

where ${ϵ}_t~N\left(0,1\right)\left[iid\right]$ .

The characteristic polynomial of the above expression is as follows: $1-{\varphi }_{1}L-{\varphi }_2{L}^2=0$ which can be written in the z-plane as: $1-{\varphi }_1{z}^{-1}-{\varphi }_2{z}^{-2}=0$ , which implies that

$z^2-{\varphi }_{1}z-{\varphi }_2=0.$

Applying the bilinear transformation on the above expression gives:

${\left(\frac{w+1}{w-1}\right)}^2-{\varphi }_1\left(\frac{w+1}{w-1}\right)-{\varphi }_2=0.$

After rearranging, we get:

$\left(1-{\varphi }_1-{\varphi }_2\right)w^2+\left(2+2{\varphi }_2\right)w+\left(1+{\varphi }_1-{\varphi }_2\right)=0.$ (7)

The Routh Array for the above expression is as follows:

The null hypothesis of stability of above AR(2) process is as follows:

$H_0:\left[\begin{array}{l}1-{\varphi }_1-{\varphi }_2\hfill \\ 2+2{\varphi }_2\hfill \\ 1+{\varphi }_1-{\varphi }_2\hfill \end{array}\right]>0$ , or $H_0:\left[\begin{array}{l}-\left(1-{\varphi }_2\right)<{\varphi }_1<\left(1-{\varphi }_2\right)\hfill \\ -1<{\varphi }_2<1\hfill \end{array}\right]$ .

Now to express the hypothesis testing in a general form, we can write a univariate time series as follows:

$y_t={x^{\prime }}_t\beta +{ϵ}_t,$ (8)

where ${ϵ}_t~N\left(0,1\right)\left[iid\right]$ ; $y_t$ is a scalar time series variable, $x_t$ is a $\left(k\times 1\right)$ vector of regressors (each regressor represents some lagged value of $y_t$ ), $\beta$ is a $\left(k\times 1\right)$ vector of parameters of interest, and ${ϵ}_t$ represents the unexplained part. After estimating the above model, we are interested in knowing whether the roots of the characteristic equation associated with the above model satisfy the stability criterion or not. Our null hypothesis is as follows:

$H_0:a<R\beta <b,$

against the alternative

$H_1:b\le R\beta ,\text{or}R\beta \le a,$

where the matrix of constraints $R$ is a $\left(p\times k\right)$ matrix of rank p, where $p\le k$ . $a$ and $b$ are known $\left(p\times 1\right)$ vectors.

4. Minimization Problem with Inequality Constraints

The main challenge in the implementation of the above test is that in order to get the constrained estimates of $\beta$ ’s, we need to do the following constrained minimization:

$\underset{\beta }{\min }{\displaystyle {\sum }_{t=1}^T{\left(y_t-{x^{\prime }}_t\beta \right)}^2}$

subject to $a<R\beta <b$ .

There are various techniques in the current Mathematics literature that allow for constrained optimization with inequality constraints such as Linear programming (for linear objective function) [35] , Quadratic programming or characterizing the problem in terms of the Karush-Kuhn-Tucker conditions (for non-linear objective function). However, these techniques require the inequality constraints to have an equality sign as well. On account of relying on these techniques, the current Statistics and Econometrics literature only allows a null hypothesis with an inequality sign if the equality sign is also there. Furthermore, this type of inequality constraint is seldom found in hypothesis testing where the parameters are bounded by a lower as well as an upper limit. To avoid a strict inequality sign in the null hypothesis, usually the problem is reframed in such a way that the null hypothesis involves a less than or equal to, or greater than or equal to sign. However, in this situation, if we try to make the alternative hypothesis as the null, the series can be unstable for a variety of reasons making the distribution under the null nearly impossible to calculate.

In general, a constraint with an upper and lower bound on parameters is more convenient and practical to test as compared to an equality constraint since the exact value of a parameter is hardly known. If $a=0$ , and $b=\infty$ , then the above constraint is equivalent to the positiveness constraint of parameters. As the difference between $a$ and $b$ gets small, the inequality constraint approaches the equality constraint. In this regard, the above constraint is a general version of other type of constraints, and a methodology to handle this could have implications for hypothesis testing in general in econometrics, and solving the optimization problems with inequality constraints in Mathematics.

4.1. Methodology

4.1.1. $\beta$ a Scalar

The inequality constraint is $a<\beta <b$ . Let

$k=\ln \left[\frac{\beta -a}{b-\beta }\right]\text{.}$

Taking the partial differential of the above expression, we get:

$\delta k=\frac{b-a}{\left[\beta -a\right]\left[b-\beta \right]}\delta \beta ,$ (9)

$\delta \beta$ is calculated as follows:

$\delta \beta ={\left({\displaystyle {\sum }_{t=1}^T{x}_t{x^{\prime }}_t}\right)}^{-1}\left({\displaystyle {\sum }_{t=1}^T\left(x_t{y}_t-{x^{\prime }}_t{\beta }_0\right)}\right)\text{.}$ (10)

(See the detail of this step in the next section.)

Plugging this in Equation (9), we get $\delta k$ which we can use in the following expression:

$k=k_0+\delta k,$ (11)

$k_0$ is the initial value of k, i.e.,

$k_0=\ln \left[\frac{{\beta }_0-a}{b-{\beta }_0}\right]\text{.}$

From Equation (11), we have

$\ln \left[\frac{\beta -a}{b-\beta }\right]-\ln \left[\frac{{\beta }_0-a}{b-{\beta }_0}\right]=\delta k,$

$\Rightarrow \ln \left[\frac{\left(\beta -a\right)\left(b-{\beta }_0\right)}{\left(b-\beta \right)\left({\beta }_0-a\right)}\right]=\delta k,$

$\Rightarrow \left[\frac{\left(\beta -a\right)\left(b-{\beta }_0\right)}{\left(b-\beta \right)\left({\beta }_0-a\right)}\right]=\exp \left(\delta k\right),$

$\Rightarrow \left(\beta -a\right)\left(b-{\beta }_0\right)=\left(b-\beta \right)\left({\beta }_0-a\right)\exp \left(\delta k\right),$

$\Rightarrow \beta \left[\left(b-{\beta }_0\right)+\left({\beta }_0-a\right)\exp \left(\delta k\right)\right]=a\left(b-{\beta }_0\right)+b\left({\beta }_0-a\right)\exp \left(\delta k\right)\mathrm{,}$

$\Rightarrow {\beta }^{updated}=\frac{a\left(b-{\beta }_0\right)+b\left({\beta }_0-a\right)\exp \left(\delta k\right)}{\left(b-{\beta }_0\right)+\left({\beta }_0-a\right)\exp \left(\delta k\right)}\text{.}$

Now treat the updated $\beta$ as ${\beta }_0$ in Equation (10) and repeat the whole procedure until the estimate converges.

4.1.2. $\beta$ a Vector

Let $A\equiv \underset{\left(p\times p\right)}{\left[\begin{array}{cccc}{a}_1& 0& \cdots & 0\\ 0& a_2& \cdots & 0\\ ⋮& ⋮& \ddots & ⋮\\ 0& 0& \cdots & a_p\end{array}\right]}$ ; $B\equiv \underset{\left(p\times p\right)}{\left[\begin{array}{cccc}{b}_1& 0& \cdots & 0\\ 0& b_2& \cdots & 0\\ ⋮& ⋮& \ddots & ⋮\\ 0& 0& \cdots & b_p\end{array}\right]}$ ; $D\equiv \underset{\left(p\times p\right)}{\left[\begin{array}{cccc}{\left(Rß\right)}_1& 0& \cdots & 0\\ 0& {\left(Rß\right)}_2& \cdots & 0\\ ⋮& ⋮& \ddots & ⋮\\ 0& 0& \cdots & {\left(Rß\right)}_p\end{array}\right]}$ ; $e_p\equiv \underset{\left(p\times 1\right)}{\left[\begin{array}{c}1\\ 1\\ ⋮\\ 1\end{array}\right]}$ ; $a=A{e}_p=\underset{\left(p\times 1\right)}{\left[\begin{array}{c}{a}_1\\ a_2\\ ⋮\\ a_p\end{array}\right]}$ ; $b=B{e}_p=\underset{\left(p\times 1\right)}{\left[\begin{array}{c}{b}_1\\ b_2\\ ⋮\\ b_p\end{array}\right]}$ ; and $Rß=D{e}_p=\underset{\left(p\times 1\right)}{\left[\begin{array}{c}{\left(Rß\right)}_1\\ {\left(Rß\right)}_2\\ ⋮\\ {\left(Rß\right)}_p\end{array}\right]}$ .

The inequality constraint $a<Rß<b$ implies that $a_1<{\left(Rß\right)}_1<b_1$ , $a_2<{\left(Rß\right)}_2<b_2$ , $\cdots$ , $a_p<{\left(Rß\right)}_p<b_p$ .

Let $k_1=\ln \left[\frac{{\left(Rß\right)}_1-a_1}{b_1-{\left(Rß\right)}_1}\right]\mathrm{,}{k}_2=\ln \left[\frac{{\left(Rß\right)}_2-a_2}{b_2-{\left(Rß\right)}_2}\right]\mathrm{,}\cdots \mathrm{,}{k}_p=\ln \left[\frac{{\left(Rß\right)}_p-a_p}{b_p-{\left(Rß\right)}_p}\right]$ ,

$K\equiv \underset{\left(p\times p\right)}{\left[\begin{array}{cccc}{k}_1& 0& \cdots & 0\\ 0& k_2& \cdots & 0\\ ⋮& ⋮& \ddots & ⋮\\ 0& 0& \cdots & k_p\end{array}\right]}=\underset{\left(p\times p\right)}{\left[\begin{array}{cccc}\ln \left[\frac{{\left(Rß\right)}_1-a_1}{b_1-{\left(Rß\right)}_1}\right]& 0& \cdots & 0\\ 0& \ln \left[\frac{{\left(Rß\right)}_2-a_2}{b_2-{\left(Rß\right)}_2}\right]& \cdots & 0\\ ⋮& ⋮& \ddots & ⋮\\ 0& 0& \cdots & \ln \left[\frac{{\left(Rß\right)}_p-a_p}{b_p-{\left(Rß\right)}_p}\right]\end{array}\right]}$ or

$K=\ln d\left[{\left[B-D\right]}^{-1}\left[D-A\right]\right]\mathrm{,}$ (12)

where $\ln d[.]\equiv$ log of diagonal of the matrix

(Note: This is just a defining expression. It does not say that we can take the log of only the diagonal of a matrix.)

$k=K{e}_p=\underset{\left(p\times 1\right)}{\left[\begin{array}{c}{k}_1\\ k_2\\ ⋮\\ k_p\end{array}\right]}=\underset{\left(p\times 1\right)}{\left[\begin{array}{c}\ln \left[\frac{{\left(R\text{ß}\right)}_1-a_1}{b_1-{\left(R\text{ß}\right)}_1}\right]\\ \ln \left[\frac{{\left(R\text{ß}\right)}_2-a_2}{b_2-{\left(R\text{ß}\right)}_2}\right]\\ ⋮\\ \ln \left[\frac{{\left(R\text{ß}\right)}_p-a_p}{b_p-{\left(R\text{ß}\right)}_p}\right]\end{array}\right]}=\ln d\left[{\left[B-D\right]}^{-1}\left[D-A\right]\right]e_p$ ,

${\left[B-D\right]}^{-1}=\left[\begin{array}{cccc}{\left[b_1-{\left(Rß\right)}_1\right]}^{-1}& 0& \cdots & 0\\ 0& {\left[b_2-{\left(Rß\right)}_2\right]}^{-1}& \cdots & 0\\ ⋮& ⋮& \ddots & ⋮\\ 0& 0& \cdots & {\left[b_p-{\left(Rß\right)}_p\right]}^{-1}\end{array}\right]$ ; and

$\left[D-A\right]=\left[\begin{array}{cccc}\left[{\left(Rß\right)}_1-a_1\right]& 0& \cdots & 0\\ 0& \left[{\left(Rß\right)}_2-a_2\right]& \cdots & 0\\ ⋮& ⋮& \ddots & ⋮\\ 0& 0& \cdots & \left[{\left(Rß\right)}_p-a_p\right]\end{array}\right]$ .

Taking the partial differential of the expression for $k_1$ , we get:

$\delta k_1=\frac{b_1-a_1}{\left[{\left(R\text{ß}\right)}_1-a_1\right]\left[b_1-{\left(R\text{ß}\right)}_1\right]}\delta {\left(R\text{ß}\right)}_1\mathrm{.}$

Similarly

$\delta k_2=\frac{b_2-a_2}{\left[{\left(R\text{ß}\right)}_2-a_2\right]\left[b_2-{\left(R\text{ß}\right)}_2\right]}\delta {\left(R\text{ß}\right)}_2\mathrm{,}$

$⋮$

$\delta k_p=\frac{b_p-a_p}{\left[{\left(R\text{ß}\right)}_p-a_p\right]\left[b_p-{\left(R\text{ß}\right)}_p\right]}\delta {\left(R\text{ß}\right)}_p\mathrm{.}$

$\delta k_1$ and $\delta {\left(Rß\right)}_1$ have the same sign as the multiplier of $\delta {\left(Rß\right)}_1$ is positive. In matrix notation, we can write:

$\underset{\left(p\times 1\right)}{\delta k}=\underset{\left(p\times p\right)}{{\left[D-A\right]}^{-1}}\underset{\left(p\times p\right)}{{\left[B-D\right]}^{-1}}\underset{\left(p\times p\right)}{\left[B-A\right]}\underset{\left(p\times 1\right)}{\delta \left(R\text{ß}\right)}$ , or

$\underset{\left(p\times 1\right)}{\delta k}=\left[\begin{array}{c}\delta k_1\\ \delta k_2\\ ⋮\\ \delta k_p\end{array}\right]=\left[\begin{array}{c}\frac{b_1-a_1}{\left[{\left(R\text{ß}\right)}_1-a_1\right]\left[b_1-{\left(R\text{ß}\right)}_1\right]}\delta {\left(R\text{ß}\right)}_1\\ \frac{b_2-a_2}{\left[{\left(R\text{ß}\right)}_2-a_2\right]\left[b_2-{\left(R\text{ß}\right)}_2\right]}\delta {\left(R\text{ß}\right)}_2\\ ⋮\\ \frac{b_p-a_p}{\left[{\left(R\text{ß}\right)}_p-a_p\right]\left[b_p-{\left(R\text{ß}\right)}_p\right]}\delta {\left(R\text{ß}\right)}_p\end{array}\right]\mathrm{,}$ (13)

${\left[D-A\right]}^{-1}=\left[\begin{array}{cccc}{\left[{\left(Rß\right)}_1-a_1\right]}^{-1}& 0& \cdots & 0\\ 0& {\left[{\left(Rß\right)}_2-a_2\right]}^{-1}& \cdots & 0\\ ⋮& ⋮& \ddots & ⋮\\ 0& 0& \cdots & {\left[{\left(Rß\right)}_p-a_p\right]}^{-1}\end{array}\right]$ ; and $\delta \left(Rß\right)=\left[\begin{array}{c}\delta {\left(R\text{ß}\right)}_1\\ \delta {\left(R\text{ß}\right)}_2\\ ⋮\\ \delta {\left(R\text{ß}\right)}_p\end{array}\right]$ .

$\delta \left(Rß\right)=R\cdot \delta ß$

Similarly

$\delta K=\underset{\left(p\times p\right)}{\left[\begin{array}{cccc}\delta k_1& 0& \cdots & 0\\ 0& \delta k_2& \cdots & 0\\ ⋮& ⋮& \ddots & ⋮\\ 0& 0& \cdots & \delta k_p\end{array}\right]}$ , or

$\delta K=\underset{\left(p\times p\right)}{\left[\begin{array}{cccc}\frac{\left(b_1-a_1\right)\cdot \delta {\left(R\text{ß}\right)}_1}{\left[{\left(R\text{ß}\right)}_1-a_1\right]\left[b_1-{\left(R\text{ß}\right)}_1\right]}& 0& \cdots & 0\\ 0& \frac{\left(b_2-a_2\right)\cdot \delta {\left(R\text{ß}\right)}_2}{\left[{\left(R\text{ß}\right)}_2-a_2\right]\left[b_p-{\left(R\text{ß}\right)}_2\right]}& \cdots & 0\\ ⋮& ⋮& \ddots & ⋮\\ 0& 0& \cdots & \frac{\left(b_p-a_p\right)\cdot \delta {\left(R\text{ß}\right)}_p}{\left[{\left(R\text{ß}\right)}_p-a_p\right]\left[b_p-{\left(R\text{ß}\right)}_p\right]}\end{array}\right]}$ (14)

$\delta ß$ is calculated as follows:

Let $f\left(x_t,y_t;ß\right)$ $={\displaystyle {\sum }_{t=1}^T{\left(y_t-{x^{\prime }}_t\text{ß}\right)}^2}$ , which we want to minimize with respect to ß, i.e., we want to find

$\frac{\partial f}{\partial \text{ß}}=0,$

4.1.3. Newton-Raphson Method

The equation of the tangent line to the curve

$Z=\frac{\partial f}{\partial \text{ß}},$

at point ß₀ is

(15)

Setting $Z\left(ß\right)=0$ , and $ß={ß}_{\text{1}}$ gives:

$\Rightarrow {\text{ß}}_1={\text{ß}}_0+{\left({\displaystyle {\sum }_{t=1}^T{x}_t{x^{\prime }}_t}\right)}^{-1}\left({\displaystyle {\sum }_{t=1}^T{x}_t}\left(y_t-{x^{\prime }}_t{\text{ß}}_0\right)\right)\mathrm{.}$

The above equation can also be written as

$\delta \text{ß}={\left({\displaystyle {\sum }_{t=1}^T{x}_t{x^{\prime }}_t}\right)}^{-1}\left({\displaystyle {\sum }_{t=1}^T{x}_t}\left(y_t-{x^{\prime }}_t{\text{ß}}_0\right)\right)\mathrm{.}$ (16)

Putting this in Equation (14), we get $\delta K$ which we can use in the following expression:

$K=K_0+\delta K,$ (17)

where $K_0$ is the initial value of $K$ , i.e.,

$K_0=\underset{\left(p\times p\right)}{\left[\begin{array}{cccc}\ln \left[\frac{{\left(R{\text{ß}}_0\right)}_1-a_1}{b_1-{\left(R{\text{ß}}_0\right)}_1}\right]& 0& \cdots & 0\\ 0& \ln \left[\frac{{\left(R{\text{ß}}_0\right)}_2-a_2}{b_2-{\left(R{\text{ß}}_0\right)}_2}\right]& \cdots & 0\\ ⋮& ⋮& \ddots & ⋮\\ 0& 0& \cdots & \ln \left[\frac{{\left(R{\text{ß}}_0\right)}_p-a_p}{b_p-{\left(R{\text{ß}}_0\right)}_p}\right]\end{array}\right]}\mathrm{.}$

Putting Equation (12) into Equation (17), we get:

$\ln d\left[{\left[B-D\right]}^{-1}\left[D-A\right]\right]-\ln d\left[{\left[B-D_0\right]}^{-1}\left[D_0-A\right]\right]=\delta K,$

$\Rightarrow \ln d\left[{\left[D_0-A\right]}^{-1}\left[B-D_0\right]{\left[B-D\right]}^{-1}\left[D-A\right]\right]=\delta K,$

$\Rightarrow {\left[D_0-A\right]}^{-1}\left[B-D_0\right]{\left[B-D\right]}^{-1}\left[D-A\right]=\exp d\left(\delta K\right),$

where $\exp d\left(\mathrm{.}\right)\equiv$ exponential of the diagonal of the matrix

(Note: This is just a defining expression. It does not say that we can take the exponential of only the diagonal of a matrix.)

$\Rightarrow \left[D-A\right]=\left[B-D\right]{\left[B-D_0\right]}^{-1}\left[D_0-A\right]\exp d\left(\delta K\right),$

$\Rightarrow D=A+\left[B-D\right]{\left[B-D_0\right]}^{-1}\left[D_0-A\right]\exp d\left(\delta K\right),$

$\Rightarrow D\left\{I+{\left[B-D_0\right]}^{-1}\left[D_0-A\right]\exp d\left(\delta K\right)\right\}=A+B{\left[B-D_0\right]}^{-1}\left[D_0-A\right]\exp d\left(\delta K\right),$

$\Rightarrow D\left\{\left[B-D_0\right]+\left[D_0-A\right]\exp d\left(\delta K\right)\right\}=\left[B-D_0\right]A+B\left[D_0-A\right]\exp d\left(\delta K\right),$

$\left(\text{as}B{\left[B-D_0\right]}^{-1}={\left[B-D_0\right]}^{-1}B\right),$

$\Rightarrow D={\left\{\left[B-D_0\right]+\left[D_0-A\right]\exp d\left(\delta K\right)\right\}}^{-1}\left\{\left[B-D_0\right]A+B\left[D_0-A\right]\exp d\left(\delta K\right)\right\}\mathrm{,}$

$\Rightarrow D{e}_p={\left\{\left[B-D_0\right]+\left[D_0-A\right]\exp d\left(\delta K\right)\right\}}^{-1}\left\{\left[B-D_0\right]A+\left[D_0-A\right]\exp d\left(\delta K\right)\cdot B\right\}{e}_p,$

$\left(\text{as}B\left[D_0-A\right]\exp d\left(\delta K\right)=\left[D_0-A\right]\exp d\left(\delta K\right)\cdot B\right).$

This implies that

$R{\text{ß}}^{updated}={\left\{\left[B-D_0\right]+\left[D_0-A\right]\exp d\left(\delta K\right)\right\}}^{-1}\left\{\left[B-D_0\right]a+\left[D_0-A\right]\exp d\left(\delta K\right)\cdot b\right\}.$ (18)

In expanded form, the above expression can be written as follows:

$R{\text{ß}}^{updated}=\underset{\left(p\times 1\right)}{\left[\begin{array}{c}{\left(R{\text{ß}}^{updated}\right)}_1\\ {\left(R{\text{ß}}^{updated}\right)}_2\\ ⋮\\ {\left(R{\text{ß}}^{updated}\right)}_p\end{array}\right]}=\underset{\left(p\times 1\right)}{\left[\begin{array}{c}\frac{a_1\left[b_1-{\left(R{\text{ß}}_0\right)}_1\right]+b_1\left[{\left(R{\text{ß}}_0\right)}_1-a_1\right]\exp \left(\delta k_1\right)}{\left[b_1-{\left(R{\text{ß}}_0\right)}_1\right]+\left[{\left(R{\text{ß}}_0\right)}_1-a_1\right]\exp \left(\delta k_1\right)}\\ \frac{a_2\left[b_2-{\left(R{\text{ß}}_0\right)}_2\right]+b_2\left[{\left(R{\text{ß}}_0\right)}_2-a_2\right]\exp \left(\delta k_2\right)}{\left[b_2-{\left(R{\text{ß}}_0\right)}_2\right]+\left[{\left(R{\text{ß}}_0\right)}_2-a_2\right]\exp \left(\delta k_2\right)}\\ ⋮\\ \frac{a_p\left[b_p-{\left(R{\text{ß}}_0\right)}_p\right]+b_p\left[{\left(R{\text{ß}}_0\right)}_p-a_p\right]\exp \left(\delta k_p\right)}{\left[b_p-{\left(R{\text{ß}}_0\right)}_p\right]+\left[{\left(R{\text{ß}}_0\right)}_p-a_p\right]\exp \left(\delta k_p\right)}\end{array}\right]}\mathrm{.}$ (19)

We need to keep the restricted estimate as far from the boundaries as if the unrestricted estimate was equal to a boundary value, i.e., either $a$ or $b$ , we were able to reject the null. Therefore, we can allow the restricted estimate to take a value at the most at that point, e.g., if ß is a scalar, and the restricted estimate is converging toward a value close to b, we need a t-statistic (at 95% confidence level and 100 degrees of free1dom) as follows:

$t=\frac{b-\tilde{\beta }}{se\left(\hat{\beta }\right)}\geqq \mathrm{1.984.}$

This implies that

$\tilde{\beta }\le b-1.984\ast se\left(\hat{\beta }\right)\mathrm{.}$

This leads to the intended performance of the test, i.e., if the unrestricted estimate is less than b, we shall be able “not to reject” the null, whereas if the unrestricted estimate is greater than or equal to b, we will be able to reject the null.

4.1.4. Summary

Now let us summarize the step-wise methodology as follows:

1) Assume some initial values of ß, i.e., ß₀ in the interval $a<Rß<b$ .

2) Calculate $\delta ß$ from Equation (16) using the assumed initial values ß₀.

3) Plug this value of $\delta ß$ along with the initial values ß₀ in Equation (14) to get $\delta K$ .

4) Plug the value of $\delta K$ in Equation (18) to get updated ß.

5) Now treat the updated ß as ß₀ and repeat steps 2 to 5.

6) Stop the algorithm (if necessary) at the point discussed above.

5. Stability Test

Consider an ARMA model of the form

$Y_t=c+{\varphi }_1{Y}_{t-1}+{\varphi }_2{Y}_{t-2}+\cdots +{\varphi }_p{Y}_{t-p}+{ϵ}_t+{\theta }_1{ϵ}_{t-1}+{\theta }_2{ϵ}_{t-2}+\cdots +{\theta }_q{ϵ}_{t-q}\mathrm{,}$

with ${ϵ}_t$ as white noise:

$E\left({ϵ}_t\right)=\mathrm{0,}$

$E\left({ϵ}_t{ϵ}_{\tau }\right)=\{\begin{array}{l}{\sigma }^2\text{for}t=\tau \mathrm{,}\\ 0\text{otherwise}\text{.}\end{array}$

Let $\beta \equiv {\left(c,{\varphi }_1\mathrm{,}\cdots \mathrm{,}{\varphi }_p\mathrm{,}{\theta }_1\mathrm{,}\cdots \mathrm{,}{\theta }_q\mathrm{,}{\sigma }^2\right)}^{\prime }$ denote the $k\times 1$ vector of population parameters.

Theorem 4 Suppose that $y_1\mathrm{,}{y}_2\mathrm{,}\cdots \mathrm{,}{y}_T$ have the joint probability density

$f_{Y_T\mathrm{,}{Y}_{T-1}\mathrm{,}\cdots \mathrm{,}{Y}_1}\left(y_T\mathrm{,}{y}_{T-1}\mathrm{,}\cdots \mathrm{,}{y}_1\mathrm{<mo>\; ;\; </mo>}\beta \right)\mathrm{,}$

$\beta \in \Theta$ , and $f_{Y_T\mathrm{,}{Y}_{T-1}\mathrm{,}\cdots \mathrm{,}{Y}_1}\left(y_T\mathrm{,}{y}_{T-1}\mathrm{,}\cdots \mathrm{,}{y}_1\mathrm{<mo>\; ;\; </mo>}\beta \right)$ satisfies the following assumptions:

Assumption 1: Identifiability; ${\beta }_1\ne {\beta }_2$ implies that

$F_{Y_T\mathrm{,}{Y}_{T-1}\mathrm{,}\cdots \mathrm{,}{Y}_1}\left(y_T\mathrm{,}{y}_{T-1}\mathrm{,}\cdots \mathrm{,}{y}_1\mathrm{<mo>\; ;\; </mo>}{\beta }_1\right)\ne F_{Y_T\mathrm{,}{Y}_{T-1}\mathrm{,}\cdots \mathrm{,}{Y}_1}\left(y_T\mathrm{,}{y}_{T-1}\mathrm{,}\cdots \mathrm{,}{y}_1\mathrm{<mo>\; ;\; </mo>}{\beta }_2\right)\mathrm{.}$

Assumption 2: For each $\beta \in \Theta$ , $F_{Y_T\mathrm{,}{Y}_{T-1}\mathrm{,}\cdots \mathrm{,}{Y}_1}\left(y_T\mathrm{,}{y}_{T-1}\mathrm{,}\cdots \mathrm{,}{y}_1\mathrm{<mo>\; ;\; </mo>}\beta \right)$ has the same support not depending on $\beta$ .

Assumption 3: For each $\beta \in \Theta$ , the first three partial derivatives of

$\left(\log f_{Y_T\mathrm{,}{Y}_{T-1}\mathrm{,}\cdots \mathrm{,}{Y}_1}\left(y_T\mathrm{,}{y}_{T-1}\mathrm{,}\cdots \mathrm{,}{y}_1\mathrm{<mo>\; ;\; </mo>}\beta \right)\right)\mathrm{,}$

with respect to $\beta$ exist for $y_T\mathrm{,}{y}_{T-1}\mathrm{,}\cdots \mathrm{,}{y}_1$ in the support of

$F_{Y_T\mathrm{,}{Y}_{T-1}\mathrm{,}\cdots \mathrm{,}{Y}_1}\left(y_T\mathrm{,}{y}_{T-1}\mathrm{,}\cdots \mathrm{,}{y}_1\mathrm{<mo>\; ;\; </mo>}\beta \right).$

Assumption 4: For each $\beta \in \Theta$ , there exists a function $g\left(y_T\mathrm{,}{y}_{T-1}\mathrm{,}\cdots \mathrm{,}{y}_1\right)$ (possibly depending on $\beta$ ), such that in a neighborhood of the given $\beta$ and for all $l\mathrm{,}m\mathrm{,}n\in \left\{\mathrm{1,}\cdots \mathrm{,}k\right\}$ ,

$\left|\frac{{\partial }^3}{\partial {\beta }_l\partial {\beta }_m\partial {\beta }_n}\log f_{Y_T\mathrm{,}{Y}_{T-1}\mathrm{,}\cdots \mathrm{,}{Y}_1}\left(y_T\mathrm{,}{y}_{T-1}\mathrm{,}\cdots \mathrm{,}{y}_1\mathrm{<mo>\; ;\; </mo>}\beta \right)\right|\le g\left(y_T\mathrm{,}{y}_{T-1}\mathrm{,}\cdots \mathrm{,}{y}_1\right)\mathrm{,}$

where ${\displaystyle \int g}\left(y_T\mathrm{,}{y}_{T-1}\mathrm{,}\cdots \mathrm{,}{y}_1\right)\text{d}{F}_{Y_T\mathrm{,}{Y}_{T-1}\mathrm{,}\cdots \mathrm{,}{Y}_1}\left(y_T\mathrm{,}{y}_{T-1}\mathrm{,}\cdots \mathrm{,}{y}_1\mathrm{<mo>\; ;\; </mo>}\beta \right)<\infty$ .

Assumption 5: For each $\beta \in \Theta$ , $E\left[\partial \log f\left(y_T\mathrm{,}{y}_{T-1}\mathrm{,}\cdots \mathrm{,}{y}_1\mathrm{<mo>\; ;\; </mo>}\beta \right)/\partial \beta \right]=0$ ,

$\begin{array}{c}I\left(\beta \right)=E\left[\frac{\partial }{\partial \beta }\log f\left(y_T\mathrm{,}{y}_{T-1}\mathrm{,}\cdots \mathrm{,}{y}_1\mathrm{<mo>\; ;\; </mo>}\beta \right)\frac{\partial }{\partial {\beta }^{\prime }}\log f\left(y_T\mathrm{,}{y}_{T-1}\mathrm{,}\cdots \mathrm{,}{y}_1\mathrm{<mo>\; ;\; </mo>}\beta \right)\right]\\ =E\left[-\frac{{\partial }^2}{\partial \beta \partial {\beta }^{\prime }}\log f\left(y_T\mathrm{,}{y}_{T-1}\mathrm{,}\cdots \mathrm{,}{y}_1\mathrm{<mo>\; ;\; </mo>}\beta \right)\right]\mathrm{,}\end{array}$

and $I\left(\beta \right)$ is nonsingular.

$\hat{\beta }$ satisfies $S\left(\hat{\beta }\right)={\partial L\left(\beta \right)/\partial \beta |}_{\beta =\hat{\beta }}=0$ where

$L\left(\beta \right)=\log f_{Y_T\mathrm{,}{Y}_{T-1}\mathrm{,}\cdots \mathrm{,}{Y}_1}\left(y_T\mathrm{,}{y}_{T-1}\mathrm{,}\cdots \mathrm{,}{y}_1\mathrm{<mo>\; ;\; </mo>}\beta \right)\mathrm{,}$

and $\hat{\beta }\stackrel{p}{\to }\beta$ as $T\to \infty$ .

Theorem 5 Suppose that $y_1,y_2,\cdots ,y_T$ have the joint probability density

$f_{Y_T\mathrm{,}{Y}_{T-1}\mathrm{,}\cdots \mathrm{,}{Y}_1}\left(y_T\mathrm{,}{y}_{T-1}\mathrm{,}\cdots \mathrm{,}{y}_1\mathrm{<mo>\; ;\; </mo>}\beta \right),$

$\beta \in \Theta$ , and all the assumptions of theorem 4 hold, then under the null

$H_0\mathrm{<mo>\; :\; </mo>}\beta =\tilde{\beta }\mathrm{,}$

where $\tilde{\beta }$ is the restricted estimator,

$LR=-2\log \left(L\left(\tilde{\beta }\right)/L\left(\hat{\beta }\right)\right)=2\left(\log L\left(\hat{\beta }\right)-\log L\left(\tilde{\beta }\right)\right)\stackrel{d}{\to }{\chi }_p^2\left(\text{under}{H}_0\right)\mathrm{,}$

where $L\left(\tilde{\beta }\right)$ and $L\left(\hat{\beta }\right)$ denote the values of the log likelihood function at the restricted ( $\tilde{\beta }$ ) and the unrestricted ( $\hat{\beta }$ ) estimates respectively.

Theorem 6 If all the assumptions of theorem 4 hold, then under the null hypothesis that the restrictions are true

$LM=T^{-1}S{\left(\tilde{\beta }\right)}^{\prime }I{\left(\tilde{\beta }\right)}^{-1}S\left(\tilde{\beta }\right)\stackrel{d}{\to }{\chi }_p^2\left(\text{under}{H}_0\right)\mathrm{.}$

Theorem 7 If all the assumptions of theorem 4 hold and the null hypothesis is true, then

$Wald=T{\left(\hat{\beta }-\tilde{\beta }\right)}^{\prime }{\left[I{\left(\hat{\beta }\right)}^{-1}\right]}^{-1}\left(\hat{\beta }-\tilde{\beta }\right)\stackrel{d}{\to }{\chi }_p^2\left(\text{under}{H}_0\right)\mathrm{.}$

For Equation (8), and the null hypothesis of the form $R\beta =R\tilde{\beta }$ , the above Wald statistic can be written as follows:

$\begin{array}{c}W=\sqrt{T}{\left(R\hat{\beta }-R\tilde{\beta }\right)}^{\prime }{\left[{\hat{\sigma }}^{2}R{\left(T^{-1}{\displaystyle {\sum }_{t=1}^T{x}_t{x^{\prime }}_t}\right)}^{-1}{R}^{\prime }\right]}^{-1}\sqrt{T}\left(R\hat{\beta }-R\tilde{\beta }\right)\\ \stackrel{d}{\to }{\chi }_p^2\left(\text{under}{H}_0\right)\mathrm{,}\end{array}$

where $\hat{\beta }$ and $\tilde{\beta }$ are the unrestricted and the restricted estimators respectively.

Theorem 8 If an ARMA (p, q) process is misspecified and there is some serial correlation in the noise term, we need to use the robust inference given in [36] .

Theorem 9 Let us consider a VAR of the form

$Y_t=c+{\Phi }_1{Y}_{t-1}+{\Phi }_2{Y}_{t-2}+\cdots +{\Phi }_p{Y}_{t-p}+{ϵ}_t\mathrm{,}$

where $Y_t$ denote an $\left(n\times 1\right)$ vector containing the values that n variables take at date t.

${ϵ}_t~i\mathrm{.}i\mathrm{.}d\mathrm{.}\left(0\mathrm{,}\Omega \right)\mathrm{.}$

Let ${\Pi }^{\prime }=\left[\begin{array}{ccccc}c& {\Phi }_1& {\Phi }_2& \cdots & {\Phi }_p\end{array}\right]$ which is $n\times \left(np+1\right)$ matrix. $\beta \equiv Vec\Pi$ denote the $\left(\left(n^{2}p+n\right)\times 1\right)$ vector of population parameters. Suppose that we have observed each of n variables for $\left(T+p\right)$ time periods and $y_1,y_2,\cdots ,y_T$ have the conditional joint probability density

$f_{Y_T,Y_{T-1},\cdots ,Y_1|Y_0,Y_{-1},\cdots ,Y_{-p+1}}\left(y_T,y_{T-1},\cdots ,y_1|y_0,y_{-1},\cdots ,y_{-p+1}\mathrm{<mo>\; ;\; </mo>}\beta \right)\mathrm{,}$

$\beta \in \Theta$ , and the above conditional joint density satisfies the following assumptions:

Assumption 6: Identifiability; ${\beta }_1\ne {\beta }_2$ implies that

$\begin{array}{l}{F}_{Y_T,Y_{T-1},\cdots ,Y_1|Y_0,Y_{-1},\cdots ,Y_{-p+1}}\left(y_T,y_{T-1},\cdots ,y_1|y_0,y_{-1},\cdots ,y_{-p+1}\mathrm{<mo>\; ;\; </mo>}{\beta }_1\right)\\ \ne F_{Y_T,Y_{T-1},\cdots ,Y_1|Y_0,Y_{-1},\cdots ,Y_{-p+1}}\left(y_T,y_{T-1},\cdots ,y_1|y_0,y_{-1},\cdots ,y_{-p+1}\mathrm{<mo>\; ;\; </mo>}{\beta }_2\right)\mathrm{.}\end{array}$

Assumption 7: For each $\beta \in \Theta$ ,

$F_{Y_T,Y_{T-1},\cdots ,Y_1|Y_0,Y_{-1},\cdots ,Y_{-p+1}}\left(y_T,y_{T-1},\cdots ,y_1|y_0,y_{-1},\cdots ,y_{-p+1}\mathrm{<mo>\; ;\; </mo>}\beta \right)$

has the same support not depending on $\beta$ .

Assumption 8: For each $\beta \in \Theta$ , the first three partial derivatives of

$\log f_{Y_T,Y_{T-1},\cdots ,Y_1|Y_0,Y_{-1},\cdots ,Y_{-p+1}}\left(y_T,y_{T-1},\cdots ,y_1|y_0,y_{-1},\cdots ,y_{-p+1}\mathrm{<mo>\; ;\; </mo>}\beta \right),$

with respect to $\beta$ exist for $y_T,y_{T-1},\cdots ,y_1$ in the support of

$F_{Y_T,Y_{T-1},\cdots ,Y_1|Y_0,Y_{-1},\cdots ,Y_{-p+1}}\left(y_T,y_{T-1},\cdots ,y_1|y_0,y_{-1},\cdots ,y_{-p+1}\mathrm{<mo>\; ;\; </mo>}\beta \right).$

Assumption 9: For each $\beta \in \Theta$ , there exists a function $g\left(y_T,y_{T-1},\cdots ,y_1\right)$ (possibly depending on $\beta$ ), such that in a neighborhood of the given $\beta$ and for all $l\mathrm{,}m\mathrm{,}n\in \left\{\mathrm{1,}\cdots \mathrm{,}\left(n^{2}p+n\right)\right\}$ ,

$\left|\frac{{\partial }^3}{\partial {\beta }_l\partial {\beta }_m\partial {\beta }_n}\log f\left(y_T,y_{T-1},\cdots ,y_1|y_0,y_{-1},\cdots ,y_{-p+1}\mathrm{<mo>\; ;\; </mo>}\beta \right)\right|\le g\left(y_T,y_{T-1},\cdots ,y_1\right)\mathrm{,}$

where

${\displaystyle \int g}\left(y_T,y_{T-1},\cdots ,y_1\right)\text{d}F\left(y_T,y_{T-1},\cdots ,y_1|y_0,y_{-1},\cdots ,y_{-p+1}\mathrm{<mo>\; ;\; </mo>}\beta \right)<\infty \mathrm{.}$

Assumption 10: For each $\beta \in \Theta$ ,

$E\left[\partial \log f\left(y_T,y_{T-1},\cdots ,y_1|y_0,y_{-1},\cdots ,y_{-p+1}\mathrm{<mo>\; ;\; </mo>}\beta \right)/\partial \beta \right]=\mathrm{0,}$

$\begin{array}{c}I\left(\beta \right)=E[\frac{\partial }{\partial \beta }\log f\left(y_T,y_{T-1},\cdots ,y_1|y_0,y_{-1},\cdots ,y_{-p+1}\mathrm{<mo>\; ;\; </mo>}\beta \right)\\ \cdot \frac{\partial }{\partial {\beta }^{\prime }}f\left(y_T,y_{T-1},\cdots ,y_1|y_0,y_{-1},\cdots ,y_{-p+1}\mathrm{<mo>\; ;\; </mo>}\beta \right)]\\ =E\left[-\frac{{\partial }^2}{\partial \beta \partial {\beta }^{\prime }}\log f\left(y_T,y_{T-1},\cdots ,y_1|y_0,y_{-1},\cdots ,y_{-p+1}\mathrm{<mo>\; ;\; </mo>}\beta \right)\right]\mathrm{,}\end{array}$

and $I\left(\beta \right)$ is nonsingular.

$\hat{\beta }$ satisfies $S\left(\hat{\beta }\right)={\partial L\left(\beta \right)/\partial \beta |}_{\beta =\hat{\beta }}=0$ where

$L\left(\beta \right)=\log f_{Y_T,Y_{T-1},\cdots ,Y_1|Y_0,Y_{-1},\cdots ,Y_{-p+1}}\left(y_T,y_{T-1},\cdots ,y_1|y_0,y_{-1},\cdots ,y_{-p+1}\mathrm{<mo>\; ;\; </mo>}\beta \right)\mathrm{,}$

and $\hat{\beta }\stackrel{p}{\to }\beta$ as $T\to \infty$ .

Theorem 10 Suppose that $y_1,y_2,\cdots ,y_T$ have the conditional joint probability density

$f_{Y_T,Y_{T-1},\cdots ,Y_1|Y_0,Y_{-1},\cdots ,Y_{-p+1}}\left(y_T,y_{T-1},\cdots ,y_1|y_0,y_{-1},\cdots ,y_{-p+1}\right)$

$\beta \in \Theta$ , and all the assumptions of theorem 9 hold, then under the null

$H_0\mathrm{<mo>\; :\; </mo>}\beta =\tilde{\beta }\mathrm{,}$

where $\tilde{\beta }$ is the restricted estimator

$LR=-2\log \left(L\left(\tilde{\beta }\right)/L\left(\hat{\beta }\right)\right)=2\left(\log L\left(\hat{\beta }\right)-\log L\left(\tilde{\beta }\right)\right)\stackrel{d}{\to }{\chi }_p^2\left(\text{under}{H}_0\right)\mathrm{,}$

where $L\left(\tilde{\beta }\right)$ and $L\left(\hat{\beta }\right)$ denote the values of the log likelihood function at the restricted ( $\tilde{\beta }$ ) and the unrestricted ( $\hat{\beta }$ ) estimates respectively.

Theorem 11 If all the assumptions of theorem 9 hold, then under the null hypothesis that the restrictions are true

$LM=T^{-1}S{\left(\tilde{\beta }\right)}^{\prime }I{\left(\tilde{\beta }\right)}^{-1}S\left(\tilde{\beta }\right)\stackrel{d}{\to }{\chi }_p^2\left(\text{under}{H}_0\right)\mathrm{.}$

Theorem 12 If all the conditions of theorem 9 hold and the null hypothesis is true, then

$Wald=T{\left(\hat{\beta }-\tilde{\beta }\right)}^{\prime }{\left[I{\left(\hat{\beta }\right)}^{-1}\right]}^{-1}\left(\hat{\beta }-\tilde{\beta }\right)\stackrel{d}{\to }{\chi }_p^2\left(\text{under}{H}_0\right)\mathrm{.}$

If each equation of the VAR mentioned in theorem 6 is expressed in the form of Equation (8), and the null hypothesis is of the form $R\beta =R\tilde{\beta }$ , the above Wald statistic can be written as follows:

$\begin{array}{c}W=\sqrt{T}{\left(R\hat{\beta }-R\tilde{\beta }\right)}^{\prime }{\left[R\left[{\hat{\Omega }}_T\otimes {\left(T^{-1}{\displaystyle {\sum }_{t=1}^T{x}_t{x^{\prime }}_t}\right)}^{-1}\right]R^{\prime }\right]}^{-1}\sqrt{T}\left(R\hat{\beta }-R\tilde{\beta }\right)\\ \stackrel{d}{\to }{\chi }_p^2\left(\text{under}{H}_0\right)\mathrm{,}\end{array}$

where $\hat{\beta }$ and $\tilde{\beta }$ are the unrestricted and the restricted estimators respectively.

6. Power and Size of the Test

Since the distribution under the alternative is unknown, therefore the power of ST has been estimated through Monte Carlo simulations as follows:

$\widehat{Pow}=\frac{1}{N}{\displaystyle {\sum }_{i=1}^{N}I}\left(S{T}_i>c_{\alpha }\right)\mathrm{,}$

where $S{T}_i$ is the test statistic for the ith Monte Carlo sample, $c_{\alpha }$ is a given critical value, and I is the indicator function having value 1 if $S{T}_i>c_{\alpha }$ and otherwise. The following data generating process has been used:

$y_t=\mu +{\varphi }_1{y}_{t-1}+{ϵ}_t\mathrm{,}$

$\mu =\mathrm{0,}$

${ϵ}_t~N\left(\mathrm{0,1}\right)\mathrm{.}$

The values of parameters have been estimated through least squares with the help of methodology described earlier. A power plot for a range of alternatives starting from a unit root, i.e., ${\varphi }_1=1$ , to ${\varphi }_1=1.08$ for $N=1000$ and $T=100$ is shown in Figure 2. The Sufficient test has the minimum power for a unit root. The power increases as we move farther from a unit root. The null rejection probabilities are reported in Table 3. For the true values satisfying the null hypothesis, the test has the correct size.

7. Monte Carlo Simulations

Some specific examples are listed below in order to highlight the performance of the test in a variety of situations, i.e., real, complex, single and multiple roots. The examples also illustrate how the test identifies the number of roots causing instability in the dynamic response. A comparison with the Augmented Dickey Fuller and ADF-GLS test is also provided for the unit root cases.

7.1. AR(1)

Suppose that the model is correctly specified, i.e., we are estimating an AR(1) process. The null hypothesis for stability of an AR(1) process is as follows:

$H_0:\left[\begin{array}{c}1-{\varphi }_1\\ 1+{\varphi }_1\end{array}\right]>0$ , or

$H_0:-1<{\varphi }_1<1$ .

Under the null hypothesis, the dynamic response is stable and the t-statistic follows an asymptotic normal distribution. The Sufficient test statistic is defined as follows:

$ST=\frac{\sqrt{T}\left(\hat{\beta }-\tilde{\beta }\right)}{\sqrt{{\hat{\sigma }}^2{\left(T^{-1}{\displaystyle {\sum }_{t=1}^T{x}_t{x^{\prime }}_t}\right)}^{-1}}}\mathrm{.}$

The following data generating process has been used for generating data:

$y_t=\mu +{\varphi }_1{y}_{t-1}+{ϵ}_t\mathrm{,}$

${ϵ}_t~N\left(\mathrm{0,1}\right)\mathrm{.}$

7.1.1. Case (a)

$\mu =0,{\varphi }_1=1.2$

For 1000 replications ( $T=100$ ), the number of rejections has been recorded, and the result is as follows:

The null hypothesis is rejected for all runs in the simulation, and the estimated response is dynamically unstable.

7.1.2. Case (b)

$\mu =0,{\varphi }_1=1$

Figure 2. Power plot for the sufficient test for an AR(1) process.

Table 3. Empirical null rejection probabilities 5% level.

For 1000 replications ( $T=100$ ), the number of rejections by the ST and the number of “No rejections” by the ADF test have been recorded, and the result is as follows:

7.1.3. Case (c)

$\mu =0,{\varphi }_1=0.99$

For 1000 replications ( $T=100$ ), the number of “No rejections” by the ST and the number of rejections by the ADF test (which is optimal in the absence of an intercept and time trend in terms of power) have been recorded, and the result is as follows:

Although we cannot directly compare the power of ST with that of the ADF test as the null hypotheses are different, however, cases (c) and (d) somehow give a reflection of the power of ST as compared to that of the Augmented Dickey Fuller test. The ST has the minimum power in case (c) (a unit root case), i.e., 0.089. As the roots move farther from a value of one, the power of ST increases. In contrast, in case (d), the ADF test has only a power of 0.047.

7.2. AR (2)

Suppose the model is correctly specified, i.e., we are estimating an AR(2) process. The null hypothesis for stability of an AR(2) process is as follows:

$H_0:\left[\begin{array}{l}1-{\varphi }_1-{\varphi }_2\hfill \\ 2+2{\varphi }_2\hfill \\ 1+{\varphi }_1-{\varphi }_2\hfill \end{array}\right]>0$ , or

$H_0:\left[\begin{array}{l}-\left(1-{\varphi }_2\right)<{\varphi }_1<\left(1-{\varphi }_2\right)\hfill \\ -1<{\varphi }_2<1\hfill \end{array}\right]$ .

Under the null hypothesis, the dynamic response is stable and the Wald statistic follows a Chi squared distribution with two degress of freedom. The Sufficient test statistic is defined as follows:

$\begin{array}{c}ST=\sqrt{T}{\left(R\hat{\beta }-R\tilde{\beta }\right)}^{\prime }{\left[{\hat{\sigma }}^{2}R{\left(T^{-1}{\displaystyle {\sum }_{t=1}^T{x}_t{x^{\prime }}_t}\right)}^{-1}{R}^{\prime }\right]}^{-1}\sqrt{T}\left(R\hat{\beta }-R\tilde{\beta }\right)\\ \stackrel{d}{\to }{\chi }_2^2\left(\text{under}{H}_0\right)\mathrm{.}\end{array}$

The following data generating process has been used for generating data:

$y_t=\mu +{\varphi }_1{y}_{t-1}+{\varphi }_2{y}_{t-2}+{ϵ}_t\mathrm{,}$

${ϵ}_t~N\left(\mathrm{0,1}\right)\mathrm{.}$

7.2.1. Case (a)

$\mu =0,{\varphi }_1=0.85,{\varphi }_2=\mathrm{0.3.}$

For 1000 replications ( $T=100$ ), the number of rejections has been recorded, and the result is as follows:

The null hypothesis is rejected for all runs in the simulation, and the estimated response is dynamically unstable. In order to see how many roots are casuing instability, a t-test is performed for all the three elements of the null hypothesis vector as follows:

$H_0:\left[\begin{array}{l}1-{\varphi }_1-{\varphi }_2\hfill \\ 2+2{\varphi }_2\hfill \\ 1+{\varphi }_1-{\varphi }_2\hfill \end{array}\right]>0,$

$H_{01}:1-{\varphi }_1-{\varphi }_2>0,$

$H_{02}:2+2{\varphi }_2>0,$

$H_{03}:1+{\varphi }_1-{\varphi }_2>0.$

For 1000 replications ( $T=100$ ), the number of rejections has been recorded, and the result is as follows:

The t-statistic for mainly the first element of null hypothesis vector is rejected (which is equivalent to one sign change in the first column of the Routh array), which implies that only one root is causing instability in the AR(2) process. As the two roots are 1.118271 and −0.268271, therefore the test correctly determines the number of roots causing instability in the response.

7.2.2. Case (b)

$\mu =0,{\varphi }_1=1.4,{\varphi }_2=-\mathrm{0.4.}$

For 1000 replications ( $T=100$ ), the number of rejections by the ST and the number of “No rejections” by the Augmented Dickey Fuller test have been recorded, and the result is as follows:

7.2.3. Case (c)

$\mu =0,{\varphi }_1=1.2,{\varphi }_2=-\mathrm{0.21.}$

For 1000 replications ( $T=100$ ), the number of “No rejections” by the ST and the number of rejections by the Augmented Dickey Fuller test have been recorded, and the result is as follows:

Again the cases (c) and (d) provide an indirect comparison of the ST and the ADF test. The ST does not perform as badly for a unit root as the ADF test does for a root close to one.

7.2.4. Case (d)

$\mu =1,{\varphi }_1=0.8,{\varphi }_2=-\mathrm{1.2.}$

This is a case of complex roots causing instability. For 1000 replications ( $T=100$ ), the number of rejections has been recorded, and the result is as follows:

The null hypothesis is rejected for all runs in the simulation, and the estimated response is dynamically unstable. In order to see how many roots are casuing instability, a t-test is performed for all the three elements of the null hypothesis vector as follows:

$H_0:\left[\begin{array}{l}1-{\varphi }_1-{\varphi }_2\hfill \\ 2+2{\varphi }_2\hfill \\ 1+{\varphi }_1-{\varphi }_2\hfill \end{array}\right]>0,$

$H_{01}:1-{\varphi }_1-{\varphi }_2>0,$

$H_{02}:2+2{\varphi }_2>0,$

$H_{03}:1+{\varphi }_1-{\varphi }_2>0.$

For 1000 replications ( $T=100$ ), the number of rejections has been recorded, and the result is as follows:

The t-statistic for the second element of null hypothesis vector is rejected (which is equivalent to two sign changes in the first column of the Routh array), which implies that two roots are causing instability in the AR(2) process. As the two roots are $0.4+1.019804i$ and $0.4-1.019804i$ , therefore the test correctly determines the number of roots causing instability in the response.

7.2.5. Case (e)

$\mu =1,{\varphi }_1=0.8,{\varphi }_2=-\mathrm{1.02.}$

This is a case of complex unit roots causing instability. For 1000 replications ( $T=100$ ), the number of rejections has been recorded, and the result is as follows:

The null hypothesis is rejected for 64.7 percent of the times, and the estimated response is dynamically unstable. In order to see how many roots are casuing instability, a t-test is performed for all the three elements of the null hypothesis vector as follows:

$H_0:\left[\begin{array}{l}1-{\varphi }_1-{\varphi }_2\hfill \\ 2+2{\varphi }_2\hfill \\ 1+{\varphi }_1-{\varphi }_2\hfill \end{array}\right]>0,$

$H_{01}:1-{\varphi }_1-{\varphi }_2>0,$

$H_{02}:2+2{\varphi }_2>0,$

$H_{03}:1+{\varphi }_1-{\varphi }_2>0.$

For 1000 replications ( $T=100$ ), the number of rejections has been recorded, and the result is as follows:

The t-statistic for the second element of null hypothesis vector is rejected. For a pair of unit roots in an AR(2) process, the third theorem of the Routh test is applicable, which implies that there are two roots which are causing instability in the AR(2) process. As the two roots are $0.4+0.927362i$ and $0.4-0.927362i$ , therefore the test correctly determines the number of roots causing instability in the response.

7.2.6. Case (f)

$\mu =1,{\varphi }_1=2,{\varphi }_2=-1.$

This is a case of multiple unit roots. For 1000 replications ( $T=100$ ), the number of rejections by the ST and the number of “No rejections” by ADF test have been recorded, and the result is as follows:

It is evident from the above table that ST has a higher power in case of multiple unit roots. In order to see how many roots are casuing instability, a t-test is performed for all the three elements of the null hypothesis vector as follows:

$H_0:\left[\begin{array}{l}1-{\varphi }_1-{\varphi }_2\hfill \\ 2+2{\varphi }_2\hfill \\ 1+{\varphi }_1-{\varphi }_2\hfill \end{array}\right]>0,$

$H_{01}:1-{\varphi }_1-{\varphi }_2>0,$

$H_{02}:2+2{\varphi }_2>0,$

$H_{03}:1+{\varphi }_1-{\varphi }_2>0.$

For 1000 replications ( $T=100$ ), the number of rejections has been recorded, and the result is as follows:

The t-statistic for the first two elements of the null hypothesis vector is rejected, which implies that there are two roots which are causing instability in the AR(2) process. As both the roots are equal to 1, therefore the test correctly determines the number of roots causing instability in the response.

7.3. VAR

Suppose we want to test the stability of the following VAR:

$y_{1t}={\mu }_1+{\varphi }_{11}{y}_{\mathrm{1,}t-1}+{\varphi }_{12}{y}_{\mathrm{2,}t-1}+{\varphi }_{13}{y}_{\mathrm{3,}t-1}+{ϵ}_{1t}\mathrm{,}$

$y_{2t}={\mu }_2+{\varphi }_{21}{y}_{\mathrm{1,}t-1}+{\varphi }_{22}{y}_{\mathrm{2,}t-1}+{\varphi }_{23}{y}_{\mathrm{3,}t-1}+{ϵ}_{2t}\mathrm{,}$

$y_{3t}={\mu }_3+{\varphi }_{31}{y}_{\mathrm{1,}t-1}+{\varphi }_{32}{y}_{\mathrm{2,}t-1}+{\varphi }_{33}{y}_{\mathrm{3,}t-1}+{ϵ}_{3t}\mathrm{,}$

${\epsilon }_{1t}\mathrm{,}{\epsilon }_{2t}\text{and}{\epsilon }_{3t}~N\left(\mathrm{0,1}\right)\mathrm{.}$

The null hypothesis for stability of the above VAR process is as follows (see appendix for detailed derivation of the null hypothesis):

$H_0:\left[\begin{array}{l}1+A_0+A_1-A_2\hfill \\ 3-3A_0-A_1-A_2\hfill \\ 1-A_1-A_0^2-A_0{A}_2\hfill \\ 1-A_0+A_1+A_2\hfill \end{array}\right]>0$ ; r

$H_0:\left[\begin{array}{l}{a}_{11}<{\varphi }_{11}<b_{11}\hfill \\ a_{12}<{\varphi }_{12}<b_{12}\hfill \\ a_{13}<{\varphi }_{13}<b_{13}\hfill \\ a_{21}<{\varphi }_{21}<b_{21}\hfill \\ a_{22}<{\varphi }_{22}<b_{22}\hfill \\ a_{23}<{\varphi }_{23}<b_{23}\hfill \\ a_{31}<{\varphi }_{31}<b_{31}\hfill \\ a_{32}<{\varphi }_{32}<b_{32}\hfill \\ a_{33}<{\varphi }_{33}<b_{33}\hfill \end{array}\right],$

or $H_0\mathrm{<mo>\; :\; </mo>}a<vec\left({\Phi }^{\text{T}}\right)<b$ ; where

$a=\left[\begin{array}{c}{a}_{11}\\ a_{12}\\ a_{13}\\ a_{21}\\ a_{22}\\ a_{23}\\ a_{31}\\ a_{32}\\ a_{33}\end{array}\right]\text{;}\text{and}b=\left[\begin{array}{c}{b}_{11}\\ b_{12}\\ b_{13}\\ b_{21}\\ b_{22}\\ b_{23}\\ b_{31}\\ b_{32}\\ b_{33}\end{array}\right]\mathrm{.}$

Under the null hypothesis, the dynamic response is stable and the Wald statistic follows a Chi squared distribution with nine degress of freedom. The Sufficient test statistic is defined as follows:

$\begin{array}{c}ST=\sqrt{T}{\left(R\hat{\beta }-R\tilde{\beta }\right)}^{\prime }{\left[R\left[{\hat{\Omega }}_T\otimes {\left(T^{-1}{\displaystyle {\sum }_{t=1}^T{x}_t{x^{\prime }}_t}\right)}^{-1}\right]R^{\prime }\right]}^{-1}\sqrt{T}\left(R\hat{\beta }-R\tilde{\beta }\right)\\ \stackrel{d}{\to }{\chi }_9^2\left(\text{under}{H}_0\right)\mathrm{,}\end{array}$

${\hat{\Omega }}_T=T^{-1}{\displaystyle {\sum }_{t=1}^T{\widehat{ϵ}}_t{{\widehat{ϵ}}^{\prime }}_t}\mathrm{,}$

${{\widehat{ϵ}}^{\prime }}_t=\left[\begin{array}{ccc}{\widehat{ϵ}}_{1t}& {\widehat{ϵ}}_{2t}& {\widehat{ϵ}}_{3t}\end{array}\right]\mathrm{.}$

7.3.1. Case (a)

$\begin{array}{l}\mu =1,{\varphi }_{11}=0.5,{\varphi }_{12}=0.4,{\varphi }_{13}=0.3,{\varphi }_{21}=0.3,\\ {\varphi }_{22}=0.4,{\varphi }_{23}=0.2,{\varphi }_{31}=0.1,{\varphi }_{32}=0.1,{\varphi }_{33}=\mathrm{0.99.}\end{array}$

For 1000 replications ( $T=100$ ), the number of rejections by ST has been recorded, and the result is as follows:

The null hypothesis is rejected for all 1000 runs in the simulation, and the estimated response is dynamically unstable. In order to see how many roots are casuing instability, a Wald test is performed for all the four elements of the null hypothesis vector as follows:

$H_0:\left[\begin{array}{l}1+A_0+A_1-A_2\hfill \\ 3-3A_0-A_1-A_2\hfill \\ 1-A_1-A_0^2-A_0{A}_2\hfill \\ 1-A_0+A_1+A_2\hfill \end{array}\right]>0,$

$H_{01}:1+A_0+A_1-A_2>0,$

$H_{02}:3-3A_0-A_1-A_2>0,$

$H_{03}:1-A_1-A_0^2-A_0{A}_2>0,$

$H_{04}:1-A_0+A_1+A_2>0.$

For 1000 replications ( $T=100$ ), the number of rejections has been recorded, and the result is as follows:

The Wald statistic for only the first element of null hypothesis vector is rejected (which is equivalent to one sign change in the first column of the Routh array), which implies that only one root is causing instability in the VAR. As the three roots are 1.137951, 0.652049 and 0.1, therefore the test correctly determines the number of roots causing instability in the response.

7.3.2. Case (b)

$\begin{array}{l}\mu =1,{\varphi }_{11}=0.1,{\varphi }_{12}=0.1,{\varphi }_{13}=0.1,{\varphi }_{21}=0.1,\\ {\varphi }_{22}=0.1,{\varphi }_{23}=0.1,{\varphi }_{31}=0.1,{\varphi }_{32}=0.1,{\varphi }_{33}=\mathrm{0.98.}\end{array}$

This is a case of a real unit root. For 1000 replications ( $T=100$ ), the number of rejections by ST has been recorded, and the result is as follows:

The null hypothesis is rejected for 702 runs in the simulation, and the estimated response is dynamically unstable. In order to see how many roots are casuing instability, a Wald test is performed for all the four elements of the null hypothesis vector as follows:

$H_0:\left[\begin{array}{l}1+A_0+A_1-A_2\hfill \\ 3-3A_0-A_1-A_2\hfill \\ 1-A_1-A_0^2-A_0{A}_2\hfill \\ 1-A_0+A_1+A_2\hfill \end{array}\right]>0,$

$H_{01}:1+A_0+A_1-A_2>0,$

$H_{02}:3-3A_0-A_1-A_2>0,$

$H_{03}:1-A_1-A_0^2-A_0{A}_2>0,$

$H_{04}:1-A_0+A_1+A_2>0.$

For 1000 replications ( $T=100$ ), the number of rejections has been recorded, and the result is as follows:

The Wald statistic for only the first element of null hypothesis vector is rejected (which is equivalent to one sign change in the first column of the Routh array), which implies that only one root is causing instability in the VAR. As the three roots are 10.175 and 0, therefore the test correctly determines the number of roots causing instability in the response.

7.3.3. Case (c)

$\begin{array}{l}\mu =1,{\varphi }_{11}=0.7,{\varphi }_{12}=-0.5,{\varphi }_{13}=-0.6,{\varphi }_{21}=0.4,\\ {\varphi }_{22}=0.1,{\varphi }_{23}=-0.2,{\varphi }_{31}=0.7,{\varphi }_{32}=0.8,{\varphi }_{33}=\mathrm{0.35.}\end{array}$

This is a case of complex unit roots. For 1000 replications ( $T=100$ ), the number of rejections by ST has been recorded, and the result is as follows:

The null hypothesis is rejected for 622 runs in the simulation, and the estimated response is dynamically unstable. In order to see how many roots are casuing instability, a Wald test is performed for all the four elements of the null hypothesis vector as follows:

$H_0:\left[\begin{array}{l}1+A_0+A_1-A_2\hfill \\ 3-3A_0-A_1-A_2\hfill \\ 1-A_1-A_0^2-A_0{A}_2\hfill \\ 1-A_0+A_1+A_2\hfill \end{array}\right]>0,$

$H_{01}:1+A_0+A_1-A_2>0,$

$H_{02}:3-3A_0-A_1-A_2>0,$

$H_{03}:1-A_1-A_0^2-A_0{A}_2>0,$

$H_{04}:1-A_0+A_1+A_2>0.$

For 1000 replications ( $T=100$ ), the number of rejections has been recorded, and the result is as follows:

The Wald statistic for only the third element of null hypothesis vector is rejected and the third theorem of the Routh test is applicable, which implies that there are two roots which are causing instability in the VAR. As the three roots are $0.126425,0.511788+0.859458i$ and $0.511788-0.859458i$ , therefore the test correctly determines the number of roots causing instability in the response.

8. Empirical Application

This empirical application is quite similar to the one in Marco Del Negro, and Frank Schorfheide (2004). The three variables used in the trivariate VAR (with a lag length of one quarter) are the US interest rates (effective federal funds rate for the first month of the quarter) and the quarterly data of percentage changes in real output growth and inflation for the period 1955:I to 2013:II. The data for real output growth come from the Bureau of Economic Analysis; the data for inflation come from the Bureau of Labor Statistics and the EFFR data source (not seasonally adjusted) is Board of Governors of the Federal Reserve System.

The null hypothesis is not rejected and the response of the VAR is dynamically stable.

9. Conclusions

In this paper, a sufficient test for dynamic stability (in the context of the roots of the characteristic polynomial) of a univariate as well as a multivariate time series has been proposed, which may test for all kinds of roots (positive and negative real unit roots, complex unit roots and roots inside the unit circle whether single or multiple) causing instability in the dynamic response. The test is much simpler in its application as the response is dynamically stable under the null. The test also indicates the number of roots causing instability in the dynamic response. In order to formulate the null hypothesis, Routh Hurwitz stability criterion (a mathematical test) is exploited which provides a necessary and sufficient condition for the stability of a dynamic response. To use the Routh stability test in the discrete data framework, bilinear transformation has been used which maps the inside of the unit circle of the z-plane into the left half of the w-plane. In order to find the restricted estimators which satisfy the Routh Hurwitz stability criterion (given the data), an algorithm for minimization of the regression objective function subject to the inequality constraints has been devised.

For the sufficient test, a, t, LR, LM and a Wald statistic are used. The t statistic follows an asymptotic normal distribution and LR, LM and Wald follow an asymptotic chi squared distribution under the null with degrees of freedom equal to the number of restrictions, when the model is correctly specified. In case of serial correlation, robust test in Nawaz (2020) has been proposed

Conflicts of Interest

The authors declare no conflicts of interest regarding the publication of this paper.

References

[1]	Dickey, D.A. and Fuller, W.A. (1979) Distribution of the Estimators for Autoregressive Time Series with a Unit Root. Journal of the American Statistical Association, 74, 427-431. https://doi.org/10.1080/01621459.1979.10482531
[2]	Bhargava, A. (1986) On the Theory of Testing for Unit Roots in Observed Time Series. Review of Economic Studies, 53, 369-384. https://doi.org/10.2307/2297634
[3]	Park, J.Y. and Choi, B. (1988) A New Approach to Testing for a Unit Root. Working Paper 88. Center for Analytic Economics, Cornell University, New York, 23.
[4]	Guilkey, D.K. and Schmidt, P. (1989) Extended Tabulations for Dickey-Fuller Tests. Economics Letters, 31, 355-357. https://doi.org/10.1016/0165-1765(89)90028-1
[5]	Kim, K. and Schmidt, P. (1990) Some Evidence on the Accuracy of Phillips-Perron Tests Using Alternative Estimates of Nuisance Parameters. Economics Letters, 34, 345-350. https://doi.org/10.1016/0165-1765(90)90143-O
[6]	Kwiatkowski, D. and Schmidt, P. (1990) Dickey-Fuller Tests with Trend. Communications in Statistics-Theory and Methods, 19, 3645-3656. https://doi.org/10.1080/03610929008830402
[7]	Schmidt, P. (1990) Dickey-Fuller Tests with Drift. In: Fomby, T.B. and Rhodes, J., Eds., Advances in Econometrics: A Research Annual. Volume 8: Cointegration, Spurious Regression, and Unit Roots, JAI Press, Greenwich, 161-200.
[8]	Guilkey, D.K. and Schmidt, P. (1991) An Interpolation Test for a Unit Root in the Presence of Drift. Biometrika, 78, 809-815. https://doi.org/10.1093/biomet/78.4.809
[9]	Schmidt, P. and Lee, J. (1991) A Modification of the Schmidt-Phillips Unit Root Test. Economics Letters, 36, 285-289. https://doi.org/10.1016/0165-1765(91)90034-I
[10]	Stock, J.H. (1991) Confidence Intervals for the Largest Autoregressive Root in U.S. Macroeconomic Time Series. Journal of Monetary Economics, 28, 435-459. https://doi.org/10.1016/0304-3932(91)90034-L
[11]	Kwiatkowski, D., Phillips, P., Schmidt, P. and Shin, Y. (1992) Testing the Null Hypothesis of Stationarity Against the Alternative of a Unit Root: How Sure Are We That Economic Time Series Have a Unit Root? Journal of Econometrics, 54, 154-178. https://doi.org/10.1016/0304-4076(92)90104-Y
[12]	Perron, P. and Vogelsang, T.J. (1992) Testing for a Unit Root in a Time Series with a Changing Mean: Corrections and Extensions. Journal of Business and Economic Statistics, 10, 467-470. https://doi.org/10.1080/07350015.1992.10509923
[13]	Schmidt, P. and Phillips, P.C. (1992) LM Tests for a Unit Root in the Presence of Deterministic Trends. Oxford Bulletin of Economics and Statistics, 54, 257-287. https://doi.org/10.1111/j.1468-0084.1992.tb00002.x
[14]	Shin, Y. and Schmidt, P. (1992) The KPSS Stationarity Test as a Unit Root Test. Economics Letters, 38, 387-392. https://doi.org/10.1016/0165-1765(92)90023-R
[15]	Hwang, J. and Schmidt, P. (1993) On the Power of Point Optimal Tests of the Trend Stationarity Hypothesis. Economics Letters, 43, 143-147. https://doi.org/10.1016/0165-1765(93)90027-A
[16]	Kim, K. and Schmidt, P. (1993) Unit Root Tests with Conditional Heteroskedasticity. Journal of Econometrics, 59, 287-300. https://doi.org/10.1016/0304-4076(93)90027-3
[17]	Schmidt, P. (1993) Some Results on Testing for Stationarity Using Data Detrended in Differences. Economics Letters, 41, 1-6. https://doi.org/10.1016/0165-1765(93)90101-H
[18]	Lee, J. and Schmidt, P. (1994) Unit Root Tests Based on Instrumental Variables Estimation. International Economic Review, 35, 449-462. https://doi.org/10.2307/2527063
[19]	Hwang, J. and Schmidt, P. (1996) Alternative Methods of Detrending and the Power of Unit Root Tests. Journal of Econometrics, 71, 227-248. https://doi.org/10.1016/0304-4076(94)01702-6
[20]	Lee, D. and Schmidt, P. (1996) On the Power of the KPSS Test of Stationarity against Fractionally Integrated Alternatives. Journal of Econometrics, 73, 285-302. https://doi.org/10.1016/0304-4076(95)01741-0
[21]	Vogelsang, T.J. and Perron, P. (1998) Additional Tests for a Unit Root Allowing for a Break in the Trend Function at an Unknown Time. International Economic Review, 39, 1073-1100. https://doi.org/10.2307/2527353
[22]	Xiao, Z. and Phillips, P.C. (1998) An ADF Coefficient Test for a Unit Root in ARMA Models of Unknown Order with Empirical Applications to the US Economy. The Econometrics Journal, 1, 27-43. https://doi.org/10.1111/1368-423X.12016
[23]	Vogelsang, T.J. (1999) Two Simple Procedures for Testing for a Unit Root When There Are Additive Outliers. Journal of Time Series Analysis, 20, 237-252. https://doi.org/10.1111/1467-9892.00135
[24]	De Jong, R.M., Amsler, C. and Schmidt, P. (2007) A Robust Version of the KPSS Test Based on Indicators. Journal of Econometrics, 137, 311-333. https://doi.org/10.1016/j.jeconom.2006.01.001
[25]	Amsler, C., Schmidt, P. and Vogelsang, T.J. (2009) The KPSS Test Using Fixed-b Critical Values: Size and Power in Highly Autocorrelated Time Series. Journal of Time Series Econometrics, 1. https://doi.org/10.2202/1941-1928.1027
[26]	Su, J.-J., Amsler, C. and Schmidt, P. (2012) A Note on the Size of the KPSS Unit Root Test. Economics Letters, 117, 697-699. https://doi.org/10.1016/j.econlet.2012.08.019
[27]	Vogelsang, T.J. and Wagner, M. (2013) A Fixed-b Perspective on the Phillips-Perron Unit Roottests. Econometric Theory, 29, 609-628. https://doi.org/10.1017/S0266466612000485
[28]	Zeren, F. and Kızılkaya, F. (2021) A New Combination of Fourier Unit Root Tests: A PPP Application for Fragile Economies. Applied Economics Letters, 28, 1707-1711. https://doi.org/10.1080/13504851.2020.1851647
[29]	Brown, R.L., Durbin, J. and Evans, J.M. (1975) Techniques for Testing the Constancy of Regression Relationships over Time. Journal of the Royal Statistical Society B, 37, 149-163. https://doi.org/10.1111/j.2517-6161.1975.tb01532.x
[30]	Chow, G.C. (1960) Tests of Equality between Sets of Coefficients in Two Linear Regressions. Econometrica: Journal of the Econometric Society, 28, 591-605. https://doi.org/10.2307/1910133
[31]	Hansen, B.E. (1992) Testing for Parameter Instability in Linear Models. Journal of Policy Modeling, 14, 517-533. https://doi.org/10.1016/0161-8938(92)90019-9
[32]	Andrews, D.W.K. and Ploberger, W. (1994) Optimal Tests When a Nuisance Parameter Is Present Only under the Alternative. Econometrica, 62, 1383-1414. https://doi.org/10.2307/2951753
[33]	Bierens, H.J. (2001) Complex Unit Roots and Business Cycles: Are They Real? Econometric Theory, 17, 962-983. https://doi.org/10.1017/S0266466601175055
[34]	Franco, C.A. and Collado, J. (2017) Comparison on Sufficient Conditions for the Stability of Hill Equation: An Arnold’s Tongues Approach. Applied Mathematics, 8, 1481-1514. https://doi.org/10.4236/am.2017.810109
[35]	Schrijver, A. (1998) Theory of Linear and Integer Programming. John Wiley & Sons, Hoboken.
[36]	Nawaz, N. (2020) Robust Inference by Sub-Sampling. Journal of Quantitative Economics, 18, 657-681. https://doi.org/10.1007/s40953-020-00207-x