Modeling I(2) Processes Using Vector Autoregressions Where the Lag Length Increases with the Sample Size

Abstract

In this paper the theory on the estimation of vector autoregressive (VAR) models for I(2) processes is extended to the case of long VAR approximation of more general processes. Hereby the order of the autoregression is allowed to tend to infinity at a certain rate depending on the sample size. We deal with unrestricted OLS estimators (in the model formulated in levels as well as in vector error correction form) as well as with two stage estimation (2SI2) in the vector error correction model (VECM) formulation. Our main results are analogous to the I(1) case: We show that the long VAR approximation leads to consistent estimates of the long and short run dynamics. Furthermore, tests on the autoregressive coefficients follow standard asymptotics. The pseudo likelihood ratio tests on the cointegrating ranks (using the Gaussian likelihood) used in the 2SI2 algorithm show under the null hypothesis the same distributions as in the case of data generating processes following finite order VARs. The same holds true for the asymptotic distribution of the long run dynamics both in the unrestricted VECM estimation and the reduced rank regression in the 2SI2 algorithm. Building on these results we show that if the data is generated by an invertible VARMA process, the VAR approximation can be used in order to derive a consistent initial estimator for subsequent pseudo likelihood optimization in the VARMA model.

1. Introduction

Many macroeconomic variables have been found to exhibit trend-like behaviour that can be modelled by using vector autoregressions (VARs). Katarina Juselius (2006) states that empirical modelling led to the development of I(1) and I(2) models since certain features of the datasets considered required including first and second differences in order to obtain stationary time series. Additionally cointegrating relations were found in the corresponding analyses. Similar findings have reoccurred numerous times in the literature for example related to money demand Johansen (1992b); Juselius (1994), inflation Banerjee et al. (2001); Georgoutsos and Kouretas (2004), interest rates and real exchange rates Johansen et al. (2007); Juselius and Assenmacher (2017); Juselius and Stillwagon (2018); Stillwagon (2018) to mention only a few sources.

The predominant methodological approach to model integration and cointegration in the I(1) and the I(2) case in the vector autoregressive (VAR) framework has been established mainly by Søren Johansen and Katarina Juselius together with a number of coauthors (see the lists of references in Johansen (1995); Juselius (2006) for details) building on vector error correction models (see Engle and Granger (1987) for early comments on the history of using error correction models for co-integrated processes). Extending the main ideas for cointegration modeling for the I(1) setting Johansen (1997) see, e.g., Johansen (1992a) suggested a representation for the I(2) case. Johansen (1997) established asymptotic distributions for the suggested two step I(2) estimator (2SI2) as an approximation to pseudo maximum likelihood estimation involving numerical optimization. Asymptotics for the corresponding likelihood ratio tests has been developed in Paruolo (1994, 1996), its asymptotic equivalence to pseudo likelihood (using the Gaussian distribution) optimization (and hence in a certain sense statistical efficiency) is shown in Paruolo (2000). However, Nielsen and Rahbek (2007) shows that in finite samples the likelihood ratio test has size advantages. The testing of restrictions on the parameters has been investigated by Boswijk and Doornik (2004); Boswijk and Paruolo (2017); Johansen and Lütkepohl (2005). Due to the implicit vector error correction (VECM) modeling, deterministic terms in the VECM produce complex deterministic terms in the solutions processes. In the I(2) context Nielsen and Rahbek (2007); Paruolo (1994, 2006); Rahbek et al. (1999); Kurita et al. (2011) discuss the impacts of deterministic terms.

As the VECM representation includes the representation of reduced rank matrices by a product of two matrices, identification conditions are of particular importance, see Juselius (2006); Mosconi and Paruolo (2013, 2017). In this context also weak exogeneity has been studied Kurita (2012); Paruolo and Rahbek (1999).

The main idea underlying the VECM approach for estimating VAR models in the I(2) context is to reparameterize the problem such that integration and cointegration properties relate to the rank of two matrices. Assuming the data generating process to be a VAR of known finite order, the rank of matrices can be tested using (pseudo) likelihood ratio tests.

Sometimes the assumption of known order is not justified. For example it is known that a subset of variables that are generated using a finite order VAR cannot be described by a finite order VAR, but instead requires a vector autoregressive moving average (VARMA) model. However, the class of VARs provides flexibility in the sense that a VAR of infinite order can represent a large set of linear dynamical systems including all invertible VARMA systems. For stationary processes Berk (1974) and Lewis and Reinsel (1985) show that by letting the order of the VAR tend to infinity at a suitable function of the sample size, consistent estimation of the underlying transfer function can be achieved for data generating processes that can be described by a VAR(∞) subject to mild assumptions on the summability of the VAR coefficients. Additionally Lewis and Reinsel (1985) also establishes asymptotic normality (in a very specific sense) of linear combinations of the estimated autoregressive coefficients. Hannan and Deistler (1988) make the concepts operational by showing that in the case of a VARMA process generating the dataset the required rate of letting the order tend to infinity can be estimated using BIC model selection.

In the case of I(1) processes the estimation theory for long VAR approximations to VARMA processes has been extended based on the techniques in the stationary case of Lewis and Reinsel in a series of papers by Saikkonen and coauthors Saikkonen (1991, 1992); Lütkepohl and Saikkonen (1997); Saikkonen and Lütkepohl (1996); Saikkonen and Luukkonen (1997). Additionally also the Johansen framework of rank restricted estimation in the VECM model has been extended to the long VAR approximations by Saikkonen and Luukkonen (1997). Bauer and Wagner (2004) provide extensions to the multi frequency I(1) case where unit roots may occur at the seasonal frequencies.

For the I(2) case no such extensions are currently known. This is the research gap this paper tries to fill: First we establish consistency and asymptotic normality of estimated autoregressive coefficients (in the sense of Lewis and Reinsel) for unrestricted ordinary least squares (OLS) estimation in the VECM representation. This can be used in order to derive Wald type tests of linear restrictions on the autoregressive parameters. Secondly, we extend the rank restricted regression techniques in the I(2) case to the long VAR approximations showing that the asymptotics (for estimated cointegrating relations, likelihood ratio tests and the two step estimation procedures) are identical in the case of long VAR approximations and VARs of finite known order. Third, we show that if the data generating process is an invertible VARMA process the long VAR system estimator can be used in order to obtain consistent initial estimators for subsequent pseudo likelihood maximization in the VARMA model class. In all results we limit ourselves to the case of no deterministic terms being included in the VECM representation. The inclusion of deterministic terms requires changing the test distribution, compare the theory contained for example in Rahbek et al. (1999).

The paper is organized as follows: In the next section the data generating process and the main assumptions are described. Section 3 then provides the results for the unrestricted estimation. Section 4 deals with rank restricted regression in the 2SI2 procedure, while Section 5 investigates the initial guess in the VARMA setting for subsequent pseudo likelihood maximization. Finally Section 6 concludes the paper. Proofs are relegated to an appendix.

Throughout the paper we will use the notation introduced by Johansen (1997): For a matrix $C\in {\mathbb{R}}^{p\times s},s<p,$ of full column rank we use the notation $\overline{C}=C{\left(C^{\prime }C\right)}^{-1}$. Furthermore, $C_{\perp }$ denotes a full column rank matrix of dimension $p\times (p-s)$ such that $C_{\perp }^{\prime }C=0$. Whenever this notation is used the particular choice of $C_{\perp }$ is not of importance. For a matrix $C=\left(C_{i,j}\right)\in {\mathbb{R}}^{p\times s}$ we let $\parallel C\parallel$ denote the Frobenius norm $\parallel C\parallel =\sqrt{{\sum }_{i=1}^p{\sum }_{j=1}^s{C}_{i,j}^2}$.

2. Data Generating Process and Assumptions

In this paper we use the following assumptions on the data generating process:

Assumption 1 (DGP).

The process${\left(y_t\right)}_{t\in \mathbb{Z}},y_t\in {\mathbb{R}}^p,$is generated from the difference equation for$t\in \mathbb{Z}$:

$${\Delta }^2{y}_t=\alpha {\beta }^{\prime }{y}_{t-1}+\Gamma \Delta y_{t-1}+\sum _{j=1}^{\infty }{\Pi }_j{\Delta }^2{y}_{t-j}+{\epsilon }_t$$

(1)

where$\alpha ,\beta \in {\mathbb{R}}^{p\times r},0\le r<p$are full column rank matrices,$\Delta =(1-L)$with L denoting the backward shift operator such that$L{\left(y_t\right)}_{t\in \mathbb{Z}}={\left(y_{t-1}\right)}_{t\in \mathbb{Z}}$. The matrix function$A\left(z\right)={(1-z)}^2{I}_p-\alpha {\beta }^{\prime }z-\Gamma z(1-z)-{\sum }_{j=1}^{\infty }{\Pi }_j{(1-z)}^2{z}^j$fulfills the special marginal stability condition that

$$\left|A\right(z\left)\right|=0implies that\left|z\right|>1orz=1.$$

(2)

Furthermore, there exists a real$\delta >0$such that the power series defining$A\left(z\right)$converges absolutely for$\left|z\right|<1+\delta$. Define${\beta }_2={\beta }_{\perp }{\eta }_{\perp },{\alpha }_2={\alpha }_{\perp }{\zeta }_{\perp }$where${\alpha }_{\perp }^{\prime }\Gamma {\beta }_{\perp }=\zeta {\eta }^{\prime },\eta ,\zeta \in {\mathbb{R}}^{(p-r)\times s}$are of full column rank$s<p-r$. Then it is assumed that the matrix

$${\alpha }_2^{\prime }(I_p+\Gamma \overline{\beta }{\overline{\alpha }}^{\prime }\Gamma -\sum _{j=1}^{\infty }{\Pi }_j){\beta }_2$$

(3)

is nonsingular.

Furthermore, the process${\left({\epsilon }_t\right)}_{t\in \mathbb{Z}}$denotes independent identically distributed (iid) white noise with mean zero and variance${\Sigma }_{ϵ}>0$.

It is well known that the conditions (2) and  (3) are necessary and sufficient for the existence of solutions to the difference equation that are I(2) processes, see for example Johansen (1992a). Moreover, note that the assumption of absolute convergence of $A\left(z\right)$ for $\left|z\right|<1+\delta$ implies that ${\sum }_{j=0}^{\infty }{j}^k\parallel {\Pi }_j\parallel <\infty$ for every $k\in \mathbb{N}$. In particular ${\sum }_{j=0}^{\infty }{j}^2\parallel {\Pi }_j\parallel <\infty$ follows as will be used frequently below.

Every vector autoregressive function $A\left(z\right)$ corresponding to the autoregression $A\left(L\right)y_t={\epsilon }_t$, that fulfills Assumption 1, allows a representation as $A\left(z\right)={(1-z)}^2{I}_p-\alpha {\beta }^{\prime }z-\Gamma z(1-z)-{\sum }_{j=1}^{\infty }{\Pi }_j{(1-z)}^2{z}^j=\tilde{g}\left(z\right)\tilde{B}\left(z\right),\tilde{B}\left(z\right)={(1-z)}^2{I}_p-\tilde{\Pi }z-\tilde{\Gamma }z(1-z),\tilde{g}\left(z\right)=I_p+{\sum }_{j=1}^{\infty }{G}_j{z}^j$. This can be seen as follows:

$$\begin{array}{ll}\hfill {\epsilon }_t& =A\left(L\right)y_t=(A\left(1\right)-\dot{A}\left(1\right)\Delta +A^{*}\left(L\right){\Delta }^2)y_t=(A\left(1\right)-\dot{A}\left(1\right)\Delta +A^{*}\left(L\right){\Delta }^2)\mathcal{B}{\mathcal{B}}^{\prime }{y}_t\hfill \\ \hfill & =\left([-\alpha ,0,0]+[\alpha ,0,0]\Delta -\Gamma \mathcal{B}\Delta +A^{*}\left(L\right)\mathcal{B}{\Delta }^2\right){\mathcal{B}}^{\prime }{y}_t\hfill \\ \hfill & =\left([-\alpha ,-\Gamma {\beta }_1,-\Gamma {\beta }_2]+[\alpha -\Gamma \beta ,A_1^{*}\left(L\right),A_2^{*}\left(L\right)]\Delta +\left[\begin{array}{c}{A}_0^{*}\left(L\right),0,0\end{array}\right]{\Delta }^2\right)\left(\begin{array}{c}{\beta }^{\prime }\\ {\beta }_1^{\prime }\Delta \\ {\beta }_2^{\prime }\Delta \end{array}\right)y_t\hfill \\ \hfill & =\left([-\alpha ,-\Gamma {\beta }_1,-\Gamma {\beta }_2+\alpha {\overline{\alpha }}^{\prime }\Gamma {\beta }_2]+[\alpha -\Gamma \beta ,A_1^{*}\left(L\right),{\tilde{A}}_2^{*}\left(L\right)]\Delta +\left[\begin{array}{c}{A}_0^{*}\left(L\right),0,-A_0^{*}\left(L\right){\overline{\alpha }}^{\prime }\Gamma {\beta }_2\end{array}\right]{\Delta }^2\right)\left(\begin{array}{c}{\beta }^{\prime }+{\overline{\alpha }}^{\prime }\Gamma {\beta }_2{\beta }_2^{\prime }\Delta \\ {\beta }_1^{\prime }\Delta \\ {\beta }_2^{\prime }\Delta \end{array}\right)y_t\hfill \\ \hfill & =\left([-\alpha ,-\Gamma {\beta }_1,{\tilde{A}}_2^{*}\left(L\right)]+[\alpha -\Gamma \beta ,A_1^{*}\left(L\right),-A_0^{*}\left(L\right){\overline{\alpha }}^{\prime }\Gamma {\beta }_2]\Delta +\left[\begin{array}{c}{A}_0^{*}\left(L\right),0,0\end{array}\right]{\Delta }^2\right)\left(\begin{array}{c}{\beta }^{\prime }+{\overline{\alpha }}^{\prime }\Gamma {\beta }_2{\beta }_2^{\prime }\Delta \\ {\beta }_1^{\prime }\Delta \\ {\beta }_2^{\prime }{\Delta }^2\end{array}\right)y_t\hfill \\ & =g\left(L\right)B\left(L\right)y_t\hfill \end{array}$$

where $\mathcal{B}=[\beta ,{\beta }_1,{\beta }_2],{\beta }_1={\beta }_{\perp }\eta ,$ is without restriction of generality assumed to be an orthonormal matrix, $A^{*}\left(L\right)\mathcal{B}=[A_0^{*}\left(L\right),A_1^{*}\left(L\right),A_2^{*}\left(L\right)],A\left(1\right)=-\alpha {\beta }^{\prime },\dot{A}\left(1\right)=-\alpha {\beta }^{\prime }+\Gamma$ and where we use that

$$\Gamma {\beta }_2-\alpha {\overline{\alpha }}^{\prime }\Gamma {\beta }_2=(I_p-\alpha {\overline{\alpha }}^{\prime })\Gamma {\beta }_2={\overline{\alpha }}_{\perp }{\alpha }_{\perp }^{\prime }\Gamma {\beta }_{\perp }{\eta }_{\perp }=0.$$

Here

$$B\left(L\right)=\left(\begin{array}{c}{\beta }^{\prime }+{\overline{\alpha }}^{\prime }\Gamma {\beta }_2{\beta }_2^{\prime }\Delta \\ {\beta }_1^{\prime }\Delta \\ {\beta }_2^{\prime }{\Delta }^2\end{array}\right).$$

In this representation

$$g\left(1\right)=\left[\begin{array}{c}-\alpha ,-\Gamma {\beta }_1,{\tilde{A}}_2^{*}\left(1\right)\end{array}\right]$$

is nonsingular due to assumption (3). Furthermore, $g\left(z\right)={\sum }_{j=0}^{\infty }{G}_j{z}^j$ is a transfer function with ${\sum }_{j=0}^{\infty }\parallel G_j\parallel j^2<\infty$ since ${\sum }_{j=1}^{\infty }\parallel {\Pi }_j\parallel j^2<\infty$ and thus the same holds for the power series coefficients $A^{*}\left(L\right)$. Since $\left|B\right(z\left)\right|\ne 0,z\ne 1$ it follows that $\left|g\right(z\left)\right|\ne 0,\left|z\right|\le 1$. Therefore

$$B\left(L\right)y_t=u_t,g\left(L\right)u_t={\epsilon }_t$$

(4)

is a VAR process. Note, however, that $g\left(0\right)=G_0\ne I_p$ in general. This constitutes a triangular representation of the process denoting $y_{1,t}={\beta }^{\prime }{y}_t\in {\mathbb{R}}^{p_1},y_{2,t}={\beta }_1^{\prime }{y}_t\in {\mathbb{R}}^{p_2},y_{3,t}={\beta }_2^{\prime }{y}_t\in {\mathbb{R}}^{p_3}$ such that

$$\begin{array}{ll}\hfill y_{1,t}& =-{\overline{\alpha }}^{\prime }\Gamma {\beta }_2\Delta y_{3,t}+u_{1,t}=A\Delta y_{3,t}+u_{1,t}A:p_1\times p_3\hfill \\ \hfill \Delta y_{2,t}& =u_{2,t},\hfill \\ \hfill {\Delta }^2{y}_{3,t}& =u_{3,t}\hfill \end{array}$$

where $u_t={[u_{1,t}^{\prime },u_{2,t}^{\prime },u_{3,t}^{\prime }]}^{\prime }$ has a VAR(∞) representation. Furthermore, defining

$$\begin{array}{ll}\hfill \tilde{B}\left(L\right)& =\mathcal{B}\left(\begin{array}{ccc}{I}_{p_1}& 0& -{\overline{\alpha }}^{\prime }\Gamma {\beta }_2\\ 0& I_{p_2}& 0\\ 0& 0& I_{p_3}\end{array}\right)B\left(L\right)={\Delta }^2{I}_p+\beta {\beta }^{\prime }L+(\beta {\beta }^{\prime }+\beta {\overline{\alpha }}^{\prime }\Gamma {\beta }_2{\beta }_2^{\prime }+{\beta }_1{\beta }_1^{\prime })L\Delta ,\hfill \\ \hfill \tilde{g}\left(L\right)& =g\left(L\right){\left(\mathcal{B}\left(\begin{array}{ccc}{I}_{p_1}& 0& -{\overline{\alpha }}^{\prime }\Gamma {\beta }_2\\ 0& I_{p_2}& 0\\ 0& 0& I_{p_3}\end{array}\right)\right)}^{-1}\hfill \end{array}$$

we obtain $A\left(L\right)=g\left(L\right)B\left(L\right)=\tilde{g}\left(L\right)\tilde{B}\left(L\right)$ such that

$$\tilde{B}\left(L\right)y_t={\Delta }^2{y}_t+\tilde{\Pi }{y}_{t-1}+\tilde{\Gamma }\Delta y_{t-1}=v_t,\tilde{g}\left(L\right)v_t={\epsilon }_t$$

is another representation of the process ${\left(y_t\right)}_{t\in \mathbb{Z}}$ with $\tilde{B}\left(0\right)=I_p$. It follows that the triangular representation can be seen as a special case where one has partial information on the matrices $\beta ,{\beta }_1,{\beta }_2$. For estimation the VECM representation is approximated using a finite order h:

$${\Delta }^2{y}_t=\Phi y_{t-1}+\Psi \Delta y_{t-1}+\sum _{j=1}^{h-2}{\Pi }_j{\Delta }^2{y}_{t-j}+e_t$$

where $e_t={\epsilon }_t+e_{1t},e_{1t}={\sum }_{j=h-1}^{\infty }{\Pi }_j{\Delta }^2{y}_{t-j}$. As in the VECM representation the dimensions of $\beta ,{\beta }_1,{\beta }_2$ are linked to the rank of the matrices $\Phi$ and ${\alpha }_{\perp }^{\prime }\Psi {\beta }_{\perp }$. Restricting these matrices to be of particular rank is simpler than imposing the equivalent restrictions in the VAR(h) representation directly.

In the following we will first investigate the unrestricted ordinary least squares estimator in the VECM representation without taking rank restrictions into account. In the second step the 2SI2 procedure as presented in Paruolo (2000) for imposing the two rank restrictions in two steps is investigated.

For both procedures the selection of the order h is of importance. In this respect the following assumption will be used:

Assumption 2 (Lag order h).

The order h is chosen subject to the following restrictions:

• $h=o\left(T^{1/5}\right)$.
• $T^{1/2}{\sum }_{j=h+1}^{\infty }\parallel {\Pi }_j\parallel \to 0$ as $T,h\to \infty$.

This condition defines an upper bound for the order which is usually directly assured during order selection using for example information criteria. The upper bound is smaller than the usual rate $T^{1/3}$ for technical reasons. The stronger bound is not needed for all results. However, the implications for practical applications are minor as for example in the range $1\le T\le 950$ we have $2.5T^{1/5}>T^{1/3}$. The second condition of Assumption 2 implies a lower bound for the increase of h as a function of the sample size. Clearly ${\sum }_{j=h+1}^{\infty }\parallel {\Pi }_j\parallel \to 0$ for $h\to \infty$. The bound implies that for $h=h\left(T\right)$ this convergence needs to be fast enough such that $T^{1/2}{\sum }_{j=h\left(T\right)+1}^{\infty }\parallel {\Pi }_j\parallel$ still converges to zero. The lower bound depends on the underlying true parameters. For invertible VARMA processes – which can be seen as the leading case – $\parallel {\Pi }_j\parallel \le C{\rho }_0^j$ for some $0\le {\rho }_0<1$. Hannan and Deistler (1988) show that for an invertible stationary VARMA process the lower bound (in this case proportional to $logT$) can be achieved asymptotically by using BIC as the order selection procedure. Thus in this case also the stronger condition ($h=o\left(T^{1/5}\right)$) is satisfied. Bauer and Wagner (2004) extend this result to the multi frequency I(1) setting. For the I(2) case no analogous result is known, although the developments of Bauer and Wagner (2004) suggest that a similar result holds also there. This is left for future research.

Therefore the difference between the ’usual’ rates and the ones assumed above are deemed to be of minor practical consequences. Thus we are not explicit in the main text as to which results hold true under the less restrictive set of results and which do not. In the appendix, we will comment on this point, however.

3. Unrestricted Estimation

In this section the results of Lewis and Reinsel (1985) and Saikkonen and Lütkepohl (1996) are extended to the I(2) case. To simplify notation define $\langle a_t,b_t\rangle ={\sum }_{t=h+1}^T{a}_t{b}_t^{\prime }$ for sequences $a_t,b_t,t=1,\dots ,T$.1 Then the unrestricted least squares estimator in the finite VECM model uses the regressor vector $Z_{t,h}={[y_{t-1}^{\prime },\Delta y_{t-1}^{\prime },{\Delta }^2{y}_{t-1}^{\prime },\dots ,{\Delta }^2{y}_{t-h+2}^{\prime }]}^{\prime }\in {\mathbb{R}}^{ph}$. The corresponding ordinary least squares estimator is given as

$$\begin{array}{ll}\hfill \left[\widehat{\Phi },\widehat{\Psi },{\widehat{\Pi }}_1,\cdots ,{\widehat{\Pi }}_{h-2}\right]& =\left[\langle {\Delta }^2{y}_t,y_{t-1}\rangle ,\langle {\Delta }^2{y}_t,\Delta y_{t-1}\rangle ,\langle {\Delta }^2{y}_t,{\Delta }^2{y}_{t-1}\rangle ,\cdots ,\langle {\Delta }^2{y}_t,{\Delta }^2{y}_{t-h+2}\rangle \right]{\langle Z_{t,h},Z_{t,h}\rangle }^{-1}\hfill \\ \hfill & =\langle {\Delta }^2{y}_t,Z_{t,h}\rangle {\langle Z_{t,h},Z_{t,h}\rangle }^{-1}.\hfill \end{array}$$

The noise covariance is estimated from the residuals as usual as

$${\widehat{\Sigma }}_{ϵ}=N^{-1}\langle {\widehat{e}}_t,{\widehat{e}}_t\rangle ,{\widehat{e}}_t={\Delta }^2{y}_t-\widehat{\Phi }{y}_{t-1}-\widehat{\Psi }\Delta y_{t-1}-\sum _{j=1}^{h-2}{\widehat{\Pi }}_j{\Delta }^2{y}_{t-j}$$

(5)

where $N=T-h$ denotes the effective sample size.

3.1. Estimation in the Triangular VECM Representation

As typical for the cointegration framework, analysis is easier in the triangular representation which separates stationary components from I(1) and I(2) processes: Let $y_t={[y_{1,t}^{\prime },y_{2,t}^{\prime },y_{3,t}^{\prime }]}^{\prime }\in {\mathbb{R}}^p$ where $y_{i,t}\in {\mathbb{R}}^{p_i}$ is such that

$$\begin{array}{ll}\hfill y_{1,t}& =A\Delta y_{3,t}+u_{1,t},\hfill \\ \hfill \Delta y_{2,t}& =u_{2,t},\hfill \\ \hfill {\Delta }^2{y}_{3,t}& =u_{3,t}\hfill \end{array}$$

where $u_t={[u_{1,t}^{\prime },u_{2,t}^{\prime },u_{3,t}^{\prime }]}^{\prime }$ has a VAR(∞) representation $g\left(L\right)u_t={\epsilon }_t$ where

$$g\left(0\right)=\left(\begin{array}{ccc}I& 0& A\\ 0& I& 0\\ 0& 0& I\end{array}\right).$$

Note, however, that using the triangular representation implies that the matrix $B\left(L\right)$ is known up the value of the matrix A. For applications this is the case only seldom.

Thus letting $g\left(z\right)=g\left(1\right)+g^{*}\left(z\right)\Delta$ we obtain

$$\begin{array}{ll}\hfill {\epsilon }_t& =g\left(L\right)\left(\begin{array}{c}{y}_{1,t}-A\Delta y_{3,t}\\ \Delta y_{2,t}\\ {\Delta }^2{y}_{3,t}\end{array}\right)=g\left(L\right)\left(\begin{array}{c}{\Delta }^2{y}_{1,t}+\Delta y_{1,t-1}+y_{1,t-1}-A{\Delta }^2{y}_{3,t}-A\Delta y_{3,t-1}\\ {\Delta }^2{y}_{2,t}+\Delta y_{2,t-1}\\ {\Delta }^2{y}_{3,t}\end{array}\right)\hfill \\ \hfill & =g\left(L\right)\left(\begin{array}{ccc}I& 0& -A\\ 0& I& 0\\ 0& 0& I\end{array}\right){\Delta }^2{y}_t+g\left(L\right)\left(\begin{array}{c}{y}_{1,t-1}\\ 0\\ 0\end{array}\right)+g\left(L\right)\left(\begin{array}{c}\Delta y_{1,t-1}-A\Delta y_{3,t-1}\\ \Delta y_{2,t-1}\\ 0\end{array}\right)\hfill \\ \hfill & =\tilde{g}\left(L\right){\Delta }^2{y}_t+[g\left(1\right)+g^{*}\left(L\right)\Delta ]\left(\begin{array}{c}{y}_{1,t-1}\\ 0\\ 0\end{array}\right)+g\left(1\right)\left(\begin{array}{c}\Delta y_{1,t-1}-A\Delta y_{3,t-1}\\ \Delta y_{2,t-1}\\ 0\end{array}\right)\hfill \\ \hfill & =\pi \left(L\right){\Delta }^2{y}_t+g\left(1\right)\left(\begin{array}{c}{y}_{1,t-1}\\ 0\\ 0\end{array}\right)+\left[\begin{array}{ccc}{\mathcal{G}}_1+{\mathcal{G}}_1^{*}& {\mathcal{G}}_2& -{\mathcal{G}}_{1}A\end{array}\right]\left(\begin{array}{c}\Delta y_{1,t-1}\\ \Delta y_{2,t-1}\\ \Delta y_{3,t-1}\end{array}\right)\hfill \\ \hfill & =\pi \left(L\right){\Delta }^2{y}_t+\left[\begin{array}{ccc}{\mathcal{G}}_1& 0& 0\end{array}\right]y_{t-1}+\left[\begin{array}{ccc}{\mathcal{G}}_1+{\mathcal{G}}_1^{*}& {\mathcal{G}}_2& -{\mathcal{G}}_{1}A\end{array}\right]\Delta y_{t-1}\hfill \end{array}$$

with $\pi \left(L\right)=I_p-{\sum }_{j=1}^{\infty }{\Pi }_j{L}^j$ leads to the corresponding VECM representation:

$${\Delta }^2{y}_t=\Phi y_{t-1}+\Psi \Delta y_{t-1}+\sum _{j=1}^{\infty }{\Pi }_j{\Delta }^2{y}_{t-j}+{\epsilon }_t.$$

Here $\mathcal{G}:=g\left(1\right)={\sum }_{j=0}^{\infty }{G}_j=[{\mathcal{G}}_1,{\mathcal{G}}_2,{\mathcal{G}}_3]$, where ${\mathcal{G}}_i$ is $p\times p_i$ for $i=1,2,3$: Similarly, ${\mathcal{G}}^{*}:=g^{*}\left(1\right)=-{\sum }_{j=0}^{\infty }j{G}_j=[{\mathcal{G}}_1^{*},{\mathcal{G}}_2^{*},{\mathcal{G}}_3^{*}]$, where ${\mathcal{G}}_i^{*}$ is $p\times p_i$ for $i=1,2,3$. The sums exists since ${\sum }_{j=1}^{\infty }\parallel G_j\parallel j^2<\infty$ by assumption. Similarly, we partition $\Phi$, $\Psi$ and ${\Pi }_j$ into $[{\Phi }_1,{\Phi }_2,{\Phi }_3]$, $[{\Psi }_1,{\Psi }_2,{\Psi }_3]$ and $[{\Pi }_{j1},{\Pi }_{j2},{\Pi }_{j3}]$, respectively. The analogous partitioning is used for estimates.

Then $\Phi =-[{\mathcal{G}}_1,0,0],\Psi =[-{\mathcal{G}}_1^{*}-{\mathcal{G}}_1,-{\mathcal{G}}_2,{\mathcal{G}}_{1}A]$. Therefore ${\Psi }_3=-{\Phi }_{1}A$. Note that in this notation the I(2) components on the right hand side are $y_{t-1,3}$, the I(1) components are $y_{t-1,1},y_{t-1,2},\Delta y_{t-1,3}$, where $y_{t-1,1}-A\Delta y_{t-1,3}$ is stationary. Thus in order to separate regressors of different integration orders in the proof (as is usually done in the literature) we use a transformation using the unknown matrix A such that the regressor $y_{t-1,1}$ is replaced by $y_{t-1,1}-A\Delta y_{t-1,3}$. Consequently the estimate ${\widehat{\Psi }}_3$ of ${\Psi }_3$ is replaced by the estimate $\widehat{\Theta }={\widehat{\Psi }}_3+{\widehat{\Phi }}_{1}A$ of $\Theta ={\Psi }_3+{\Phi }_{1}A=0$.

Based on the estimates $\widehat{\Psi }$ and $\widehat{\Phi }$ then A can be estimated as

$$\widehat{A}=-{\left({\widehat{\Phi }}_1^{\prime }{\widehat{\Sigma }}_{ϵ}^{-1}{\widehat{\Phi }}_1\right)}^{-1}{\widehat{\Phi }}_1^{\prime }{\widehat{\Sigma }}_{ϵ}^{-1}{\widehat{\Psi }}_3.$$

(6)

Here the insertion of ${\widehat{\Sigma }}_{ϵ}^{-1}$ appears somewhat arbitrary. A motivation for this choice in the I(1) case can be found in Saikkonen (1992) equation (12). However, any other positive definite matrix could be used as well. Currently there is no knowledge on the optimality of the choice suggested above.

In the asymptotic distribution of the estimation error Brownian motions occur relating to the process ${\left(u_t\right)}_{t\in \mathbb{Z}}$: Under Assumption 1 we have

$$\frac{1}{\sqrt{T}}\sum _{t=1}^{\lfloor rT\rfloor }{u}_t\Rightarrow B\left(r\right)={[B_1{\left(r\right)}^{\prime },B_c{\left(r\right)}^{\prime }]}^{\prime }={[B_1{\left(r\right)}^{\prime },B_2{\left(r\right)}^{\prime },B_3\left(r\right)]}^{\prime }$$

where $B\left(r\right),0\le r\le 1,$ denotes a Brownian motion with corresponding variance

$$\Omega =\left[\begin{array}{cc}{\Omega }_{11}& {\Omega }_{1c}\\ {\Omega }_{c1}& {\Omega }_{cc}\end{array}\right]=\left[\begin{array}{ccc}{\Omega }_{11}& {\Omega }_{12}& {\Omega }_{13}\\ {\Omega }_{21}& {\Omega }_{22}& {\Omega }_{23}\\ {\Omega }_{31}& {\Omega }_{32}& {\Omega }_{33}\end{array}\right]=g{\left(1\right)}^{-1}{\Sigma }_{ϵ}{\left(g{\left(1\right)}^{\prime }\right)}^{-1},$$

where $B_{1.c}\left(r\right)=B_1\left(r\right)-{\Omega }_{1c}{\Omega }_{cc}^{-1}{B}_c\left(r\right)$ is a $p_1$-dimensional Brownian motion, which is independent of $B_c\left(r\right)$, with covariance

$${\Omega }_{1.c}={\Omega }_{11}-{\Omega }_{1c}{\Omega }_{cc}^{-1}{\Omega }_{c1}.$$

An estimator of ${\Omega }_{1.c}$ is given by2

$${\widehat{\Omega }}_{1.c}={\left({\widehat{\Phi }}_1^{\prime }{\widehat{\Sigma }}_{ϵ}^{-1}{\widehat{\Phi }}_1\right)}^{-1}.$$

(7)

With these definitions we can state our first result of the paper (which is proved in Appendix B):

Theorem 1.

Under Assumptions 1 and 2 for the triangular VECM representation we have:

(A) Consistency:

$$\left(i\right)\widehat{\Phi }\stackrel{p}{\to }\Phi ;\left(ii\right){\widehat{\Sigma }}_{ϵ}\stackrel{p}{\to }{\Sigma }_{ϵ};\left(iii\right){\widehat{\Omega }}_{1.c}\stackrel{p}{\to }{\Omega }_{1.c};\left(iv\right)\widehat{\Psi }\stackrel{p}{\to }\Psi ;\left(v\right)\widehat{\Theta }\stackrel{p}{\to }0;\left(vi\right)\widehat{A}\stackrel{p}{\to }A.$$

(B) Asymptotic distribution of coefficients to nonstationary regressors: Under Assumptions 1 and 2 we have ($N=T-h$):

$$\begin{array}{ll}\hfill \left(i\right)[N{\widehat{\Phi }}_2,N\widehat{\Theta },N^2{\widehat{\Phi }}_3]\stackrel{d}{\to }g\left(1\right){\int }_0^{1}dB{F}^{\prime }{\left({\int }_0^{1}F{F}^{\prime }\right)}^{-1},& \left(ii\right)N(\widehat{A}-A)\stackrel{d}{\to }{\int }_0^{1}d{B}_{1.c}{L}^{\prime }{\left({\int }_0^{1}L{L}^{\prime }\right)}^{-1}\hfill \end{array}$$

(8)

where$F\left(u\right)=\left[\begin{array}{c}{B}_c\left(u\right)\\ {\int }_0^u{B}_3\left(v\right)dv\end{array}\right],F_a\left(u\right)=\left[\begin{array}{c}{B}_2\left(u\right)\\ {\int }_0^u{B}_3\left(v\right)dv\end{array}\right]$and$L\left(u\right)=B_3\left(u\right)-{\int }_0^1{B}_3{F}_a^{\prime }{\left({\int }_0^1{F}_a{F}_a^{\prime }\right)}^{-1}{F}_a\left(u\right)$.

(C) Asymptotic distribution of coefficients to stationary regressors: Let$L_h$be a sequence of$(p^2(h-2)+p(2p_1+p_2))\times J$matrices such that$L_h^{\prime }({\Gamma }_{ECM}^{-1}\otimes {\Sigma }_{ϵ})L_h\to M>0$where${\Gamma }_{ECM}=\mathbb{E}\left(X_t{X}_t^{\prime }\right)$with$X_t:={\left[u_{1,t-1}^{\prime },\Delta y_{1,t-1}^{\prime },\Delta y_{2,t-1}^{\prime },{\Delta }^2{y}_{t-1}^{\prime },\dots ,{\Delta }^2{y}_{t-h+2}^{\prime }\right]}^{\prime }$.

Let

$$\underline{\Pi }=\left[\begin{array}{cccccc}{\Phi }_1& {\Psi }_1& {\Psi }_2& {\Pi }_1& \dots & {\Pi }_{h-2}\end{array}\right].$$

Then

$$N^{\frac{1}{2}}{L}_h^{\prime }vec(\widehat{\Pi }-\underline{\Pi })\stackrel{d}{\to }N(0,M).$$

(D) Asymptotic distribution on Wald type tests: Finally letting

$${\widehat{\Gamma }}_{ECM}=N^{-1}(\langle {\tilde{X}}_t,{\tilde{X}}_t\rangle -\langle {\tilde{X}}_t,\Delta y_{3,t-1}\rangle {\langle \Delta y_{3,t-1},\Delta y_{3,t-1}\rangle }^{-1}\langle \Delta y_{3,t-1},{\tilde{X}}_t\rangle )$$

where${\tilde{X}}_t={\left[y_{1,t-1}^{\prime },\Delta y_{1,t-1}^{\prime },\Delta y_{2,t-1}^{\prime },{\Delta }^2{y}_{t-1}^{\prime },\dots ,{\Delta }^2{y}_{t-h+2}^{\prime }\right]}^{\prime }$, the Wald test for the null hypothesis$H_0:L_h^{\prime }vec\left(\underline{\Pi }\right)=l_h$is given by

$${\widehat{\lambda }}_{Wald}=N{\left(L_h^{\prime }vec\left(\widehat{\Pi }\right)-l_h\right)}^{\prime }{\left(L_h^{\prime }({\widehat{\Gamma }}_{ECM}^{-1}\otimes {\widehat{\Sigma }}_{ϵ})L_h\right)}^{-1}\left(L_h^{\prime }vec\left(\widehat{\Pi }\right)-l_h\right).$$

Then if$L_h$is such that$L_h^{\prime }({\Gamma }_{ECM}^{-1}\otimes {\Sigma }_{ϵ})L_h\to M>0$, under the null hypothesis${\widehat{\lambda }}_{Wald}\stackrel{d}{\to }{\chi }^2\left(J\right)$.

The theorem provides the asymptotic distributions of the OLS estimates in the triangular system. Note that in this somewhat special case the properties of the regressor components (stationary or not) are known such that for each entry the convergence speed is known. Correspondingly the definition of the regressor vector ${\tilde{X}}_t$ involves only lags of $y_t$ but omits all nonstationary regressors except the ones cointegrated with $\Delta y_{3,t-1}$.

The assumptions on $L_h$ are more restrictive than needed. Lewis and Reinsel (1985) and Saikkonen and Lütkepohl (1996) only require that $L_h$ has full column rank when deriving the normalized convergence to normal distribution with unit variance as the limit for

$$N^{\frac{1}{2}}{\left(L_h^{\prime }({\Gamma }_{ECM}^{-1}\otimes {\Sigma }_{ϵ})L_h\right)}^{-1/2}{L}_h^{\prime }vec(\widehat{\Pi }-\underline{\Pi }).$$

Similar arguments could be used here.

3.2. Estimation in the General VECM Representation

The previous section dealt with the special case that a triangular representation is used and hence knowledge on the matrices $[\beta ,{\beta }_1,{\beta }_2]$ is given. This section provides a result for the general case, which, however, is limited to the coefficients to the stationary components. Since a general process generated according to Assumption 1 can be rewritten into a triangular representation using the knowledge of $[\beta ,{\beta }_1,{\beta }_2]$, some asymptotic properties of the unrestricted OLS estimators can be derived from Theorem 1 for the general case (which is proved in Appendix C):

Theorem 2.

Let the regressor vector$Z_{t,h}$$=[y_{t-1}^{\prime },\Delta y_{t-1}^{\prime },{\Delta }^2{y}_{t-1}^{\prime },\dots ,{\Delta }^2{y}_{t-h+2}^{\prime }]$and define

$$\begin{array}{ll}\hfill \underline{\mathsf{\Lambda }}& =\left[\begin{array}{ccccc}\Phi & \Psi & {\Pi }_1& \dots & {\Pi }_{h-2}\end{array}\right],\tilde{\mathsf{\Lambda }}=\langle {\Delta }^2{y}_t,Z_{t,h}\rangle {\langle Z_{t,h},Z_{t,h}\rangle }^{-1},{\tilde{\Gamma }}_{ECM}=N^{-1}\langle Z_{t,h},Z_{t,h}\rangle .\hfill \end{array}$$

Then under Assumptions 1 and 2 it follows that$\tilde{\mathsf{\Lambda }}-\underline{\mathsf{\Lambda }}=o_P\left(1\right)$.

Furthermore, let$L_h\in {\mathbb{R}}^{p^2(h+2)\times J}$be such that$L_h^{\prime }({\tilde{\Gamma }}_{ECM}^{-1}\otimes {\Sigma }_{ϵ})L_h\to M>0$. Then

$$N^{\frac{1}{2}}{L}_h^{\prime }vec(\tilde{\mathsf{\Lambda }}-\underline{\mathsf{\Lambda }})\stackrel{d}{\to }N(0,M).$$

Beside consistency the theorem implies that linear combination of OLS estimators show asymptotic normality and hence standard inference, if the asymptotic variance is nonsingular. One application of such results consists in the so called ‘surplus lag’ formulation in the context of Granger causality testing, see Bauer and Maynard (2012); Dolado and Lütkepohl (1996).

Finally note that this section does not contain results with regard to the cointegrating rank or the cointegrating space. The theorem above merely allows to test coefficients corresponding to stationary regressors. Therefore the usage is limited to somewhat special situations like the surplus-lag causality tests. However, it is also relevant for impulse response analysis, compare Inoue and Kilian (2020).

4. Rank Restricted Regression

The previous sections show that for the estimators discussed in that sections full inference on all coefficients is only possible when information on the matrices $\beta ,{\beta }_1$ and ${\beta }_2$ exists. The dimensions of the matrices relate to the ranks of the matrices $\Phi =\alpha {\beta }^{\prime }$ and, conditional on this, to the rank of ${\overline{\alpha }}_{\perp }^{\prime }\Psi {\overline{\beta }}_{\perp }$. The two rank restrictions make estimation and specification more complex than in the I(1) case.

Johansen (1995) provides the two-step approach 2SI2 that can be used for estimation and specification of the two integer valued parameters $p_1$ and $p_2$. Paruolo and Rahbek (1999) extend the 2SI2 procedure suggested in section 8 of Johansen (1997). Paruolo (2000) shows that this 2SI2 procedure achieves the same asymptotic distribution as pseudo maximum likelihood estimation which could be performed subsequent to 2SI2 estimation. This makes the procedure attractive from a practical point of view. In this section we show that these approaches extend naturally to the long VAR case. The main focus here lies on the derivation of the asymptotic properties of the rank tests.

Recall the long VAR approximation given as

$${\Delta }^2{y}_t=\Phi y_{t-1}+\Psi \Delta y_{t-1}+\sum _{j=1}^{h-2}{\Pi }_j{\Delta }^2{y}_{t-j}+e_t$$

(9)

where $\Phi =\alpha {\beta }^{\prime }$ has reduced rank $r<p$ and ${\overline{\alpha }}_{\perp }^{\prime }\Psi {\overline{\beta }}_{\perp }=\zeta {\eta }^{\prime }$ has reduced rank $s<p-r$. In this notation the 2SI2 procedure works as follows: In the first step the rank constraint on ${\overline{\alpha }}_{\perp }^{\prime }\Psi {\overline{\beta }}_{\perp }$ is neglected estimating $\alpha$ and $\beta$ by using reduced-rank regression (RRR). Then in the second step the reduced rank of ${\overline{\alpha }}_{\perp }^{\prime }\Psi {\overline{\beta }}_{\perp }$ is imposed using RRR in a transformed equation.

In more detail using the Johansen notation we denote with $R_{0t}$, $R_{1t}$ and $R_{2t}$ the residuals of regressing ${\Delta }^2{y}_t$, $\Delta y_{t-1}$ and $y_{t-1}$ on ${\Delta }^2{y}_{t-1},\dots ,{\Delta }^2{y}_{t-h+2}$, respectively; then we can rewrite (9) as

$$R_{0t}=\alpha {\beta }^{\prime }{R}_{2t}+\Psi R_{1t}+{\tilde{e}}_t.$$

(10)

Concentrating out $R_{1t}$ and denoting the residuals as $R_{0.1t}$ and $R_{2.1t}$ we obtain with $S_{ij.1}=\langle R_{it},R_{jt}\rangle -\langle R_{it},R_{1t}\rangle {\langle R_{1t},R_{1t}\rangle }^{-1}\langle R_{1t},R_{jt}\rangle$ the solution to the RRR problem from solving the eigenvalue problem

$$|\lambda S_{22.1}-S_{20.1}{S}_{00.1}^{-1}{S}_{02.1}|=0,$$

(11)

with solutions $1>{\widehat{\lambda }}_1\ge \cdots \ge {\widehat{\lambda }}_p>0$ ordered with decreasing size and corresponding vectors $V=(v_1,\cdots ,v_p)$. Then as usual the trace statistic of testing the model $H_r$ with $\mathrm{rank}(\Phi )\le r$, $r<p$, in the model $H_p$ with $\mathrm{rank}(\Phi )\le p$, is given as

$$Q_r=-2logQ\left(H_r|H_p\right)=-T\sum _{i=r+1}^{p}log(1-{\widehat{\lambda }}_i).$$

(12)

The optimizers for $\alpha ,\beta$ are given by

$$\widehat{\beta }=(v_1,\cdots ,v_r),\widehat{\alpha }=S_{02.1}\widehat{\beta },{\widehat{\Sigma }}_{ϵ}=S_{00.1}-\widehat{\alpha }{\widehat{\alpha }}^{\prime }.$$

(13)

In the second step, given $\alpha$ and $\beta$ known, we can obtain by multiplying (10) by ${\overline{\alpha }}_{\perp }^{\prime }$ that

$${\overline{\alpha }}_{\perp }^{\prime }{R}_{0t}={\overline{\alpha }}_{\perp }^{\prime }\Psi ({\overline{\beta }}_{\perp }{\beta }_{\perp }^{\prime }+\overline{\beta }{\beta }^{\prime })R_{1t}+{\overline{\alpha }}_{\perp }^{\prime }{\tilde{e}}_t=\zeta {\eta }^{\prime }\left({\beta }_{\perp }^{\prime }{R}_{1t}\right)+C\left({\beta }^{\prime }{R}_{1t}\right)+{\overline{\alpha }}_{\perp }^{\prime }{\tilde{e}}_t.$$

(14)

Note that ${\beta }^{\prime }{R}_{1t}$ is stationary. Thus concentrating out C and denoting the residuals as $R_{{\overline{\alpha }}_{\perp }.\beta ,t}$ and $R_{{\beta }_{\perp }.\beta ,t}$, respectively, we can define $S_{ab.\beta }:=\langle R_{a.\beta ,t},R_{b.\beta ,t}\rangle$, for $a,b={\overline{\alpha }}_{\perp }$ or ${\beta }_{\perp }$. Then the likelihood ratio test of the model $H_{r,s}$ with $\mathrm{rank}\left(\zeta {\eta }^{\prime }\right)\le s$, $s<p-r$ in the model $H_r^0$ with $\mathrm{rank}({\overline{\alpha }}_{\perp }^{\prime }\Psi {\overline{\beta }}_{\perp })=p-r$ is given by

$$Q_{r,s}=-2logQ\left(H_{r,s}|H_r^0\right)=-T\sum _{i=s+1}^{p-r}log(1-{\widehat{\rho }}_i).$$

(15)

where $1>{\widehat{\rho }}_1\ge \cdots \ge {\widehat{\rho }}_{p-r}>0$ are the solutions of the eigenvalue problem

$$|\rho S_{{\beta }_{\perp }{\beta }_{\perp }.\beta }-S_{{\beta }_{\perp }{\overline{\underline{\alpha }}}_{\perp }.\beta }{S}_{{\overline{\underline{\alpha }}}_{\perp }{\overline{\alpha }}_{\perp }.\beta }^{-1}{S}_{{\overline{\alpha }}_{\perp }{\beta }_{\perp }.\beta }|=0,$$

(16)

and the corresponding eigenvectors are $W=(w_1,\cdots ,w_{p-r})$. Estimators of $\zeta$ and $\eta$ are given by

$$\widehat{\eta }=(w_1,\cdots ,w_s),\widehat{\zeta }=S_{{\overline{\alpha }}_{\perp }{\beta }_{\perp }.\beta }\widehat{\eta }.$$

(17)

For the 2SI2 procedure in this second step the first step estimates $\widehat{\alpha }$ and $\widehat{\beta }$ are used in place of the unknown true quantities. Then we obtain the following analogon to the results in the finite order VAR framework (the proof is given in Appendix D):

Theorem 3.

Let the data be generated according to Assumption 1 and let the VAR order fulfill Assumption 2. Then the following asymptotic results hold:

(A) The asymptotic distribution of the likelihood ratio statistic$Q_r$under the null hypothesis$H_r$is given by

$$Q_r\stackrel{d}{\to }tr\left\{{\int }_0^{1}d{W}_{†}{F}_{†}^{\prime }{\left({\int }_0^1{F}_{†}{F}_{†}^{\prime }du\right)}^{-1}{\int }_0^1{F}_{†}d{W}_{†}^{\prime }\right\}.$$

(18)

where $W_{†}={\left({\alpha }_{\perp }^{\prime }{\Sigma }_{ϵ}{\alpha }_{\perp }\right)}^{-1/2}{\alpha }_{\perp }^{\prime }W,F_a\left(u\right)=\left[\begin{array}{c}{B}_2\left(u\right)\\ {\int }_0^u{B}_3\left(v\right)dv\end{array}\right]$ and $F_{†}\left(u\right)=F_a\left(u\right)-{\int }_0^1{F}_a{B}_3^{\prime }{\left({\int }_0^1{B}_3{B}_3^{\prime }\right)}^{-1}{B}_3\left(u\right)$. This is identical to the distribution achieved in the finite VAR case.

(B) The asymptotic distribution of the likelihood ratio statistic $Q_{r,s}$ under the null hypothesis $H_{r,s}$ is given by

$$Q_{r,s}\stackrel{d}{\to }tr\left\{{\int }_0^{1}d{W}_2^{\prime }{B}_3^{\prime }{\left({\int }_0^1{B}_3{B}_3^{\prime }du\right)}^{-1}{\int }_0^1{B}_{3}d{W}_2^{\prime }\right\}.$$

(19)

where $W_2\left(u\right)={\left({\alpha }_2^{\prime }{\Sigma }_{ϵ}{\alpha }_2\right)}^{-1/2}{\alpha }_2^{\prime }W\left(u\right)$.

(C) The asymptotic distribution of the test statistic $S_{r,s}=Q_r+Q_{r,s}$ under the null hypothesis $H_{r,s}$ is given by

$$S_{r,s}\stackrel{d}{\to }tr\left\{{\int }_0^{1}d{W}_{†}{F}_{†}^{\prime }{\left({\int }_0^1{F}_{†}{F}_{†}^{\prime }du\right)}^{-1}{\int }_0^1{F}_{†}d{W}_{†}^{\prime }\right\}+tr\left\{{\int }_0^{1}d{W}_2{B}_3^{\prime }{\left({\int }_0^1{B}_3{B}_3^{\prime }du\right)}^{-1}{\int }_0^1{B}_{3}d{W}_2^{\prime }\right\}.$$

(20)

(D) Using suitable normalizations all estimators are consistent: $\widehat{\alpha }{\left(c_{\alpha }^{\prime }\widehat{\alpha }\right)}^{-1}\stackrel{p}{\to }\alpha ,\widehat{\beta }{\left(c_{\beta }^{\prime }\widehat{\beta }\right)}^{-1}\stackrel{p}{\to }\beta ,\widehat{\zeta }{\left(c_{\zeta }^{\prime }\widehat{\zeta }\right)}^{-1}\stackrel{p}{\to }\zeta ,\widehat{\eta }{\left(c_{\eta }^{\prime }\widehat{\eta }\right)}^{-1}\stackrel{p}{\to }\eta ,\widehat{\Psi }\stackrel{p}{\to }\Psi ,{\widehat{\Pi }}_j\stackrel{p}{\to }{\Pi }_j$ where for example $c_{\alpha }^{\prime }\alpha =I_r$.

(E) The asymptotic distributions of the coefficients to the nonstationary regressors are identical to the ones in the finite order VAR case stated in Paruolo (2000). The asymptotic distribution of the coefficients ${\widehat{\Pi }}_j$ are identical to the ones in Theorem 1.

The main message of the theorem is that the 2SI2 procedure shows the same asymptotic properties including the rank tests as in the finite order VAR case. As usual also restricting the coefficients for the non-stationary regressors does not influence the asymptotics for the coefficients corresponding to the stationary regressors.

Note that Paruolo (2000) shows that in the finite VAR case 2SI2 estimates have the same asymptotic distribution as pseudo maximum likelihood (pML) estimates maximizing the Gaussian likelihood. The first order conditions for the pML estimates of the coefficients to the non-stationary regressors provided in the first display on p. 548 in Paruolo (2000) depend on the data only via the matrices $S_{ij}$ defined above. These matrices depend on the lag length of the VECM only via the concentration step. The proof of our Theorem 3 shows that these terms have the same asymptotic distributions for the finite order VAR and the long VAR. Theorem 4.3 of Paruolo (2000) shows that the asymptotic distribution of the coefficients due to stationary regressors does not depend on the distribution of the coefficients corresponding to the non-stationary regressors as long as they are estimated super-consistently. Thus our results imply that also in the long VAR case the asymptotic distribution of all estimates for the 2SI2 and the pML approach is identical.

5. Initial Guess for VARMA Estimation

One usage of long VAR approximations is as preliminary estimate for VARMA model estimation. Hannan and Kavalieris (1986) provide properties of such an approach in the stationary case, Lütkepohl and Claessen (1997) extend the procedure to the I(1) case. Here we extend this idea to the I(2) case.

The goal is to provide a consistent initial guess for the estimation of a VARMA model for I(2) processes. In this respect we assume the following data generating process:

Assumption 3 (VARMA dgp).

The process ${\left(y_t\right)}_{t\in \mathbb{Z}}$ is generated as the solution to the state space equations

$$y_t=C{x}_t+{\epsilon }_t,x_{t+1}=A{x}_t+B{\epsilon }_t$$

(21)

where ${\left({\epsilon }_t\right)}_{t\in \mathbb{Z}}$ denotes white noise subject to the same assumptions as in Assumption 1.

Here $x_t\in {\mathbb{R}}^n$ is the unobserved state process. The system $(A,B,C)$ is assumed to be minimal and in the canonical form of Bauer and Wagner (2012), that is

$$A=\left[\begin{array}{cccc}{I}_c& I_c& 0& 0\\ 0& I_c& 0& 0\\ 0& 0& I_d& 0\\ 0& 0& 0& A_{•}\end{array}\right],B=\left[\begin{array}{c}{B}_1\\ B_2\\ B_3\\ B_{•}\end{array}\right],C=\left[\begin{array}{cccc}{C}_1& C_2& C_3& C_{•}\end{array}\right],$$

where $|{\lambda }_{max}\left(A_{•}\right)|<1$ (the matrix $A_{•}$ is stable), $C_1^{\prime }{C}_1=I_c,C_3^{\prime }{C}_3=I_d,C_1^{\prime }{C}_3=0,C_1^{\prime }{C}_2=0,C_2^{\prime }{C}_3=0$. Furthermore, the system is strictly minimum-phase, that is ${\rho }_0=\left|{\lambda }_{max}(A-BC)\right|<1$. Finally the matrix $\overline{A}=A-BC$ is nonsingular.

At time $t=0$ the state $x_0={[x_{0,u}^{\prime },x_{•}^{\prime }]}^{\prime },x_{0,u}\in {\mathbb{R}}^{2c+d},$ is such that $x_{0,u}$ is deterministic and $x_{0,•}={\sum }_{j=1}^{\infty }{A}_{•}^{j-1}{B}_{•}{\epsilon }_{-j}$ denotes the stationary solution to the stable part of the system.

In this situation it follows that ${\left(y_t\right)}_{t\in \mathbb{Z}}$ is an I(2) process in the definition of Bauer and Wagner (2012), that is its second difference is a stationary VARMA process. The integers c and d are connected to the integers $p_1,p_2,p_3$ via $c=p_3,d=p_2$ such that $p_1=p-c-d$. It can furthermore be shown that a process generated using Assumption 3 possesses a VAR(h) approximation:

$$y_t+\sum _{j=1}^h{A}_j{y}_{t-j}={\epsilon }_t+C{(A-BC)}^h{x}_{t-h}$$

where $A_j=-C{(A-BC)}^{j-1}B,\parallel A_j\parallel \le \mu {\rho }^j$ ($0\le {\rho }_0<\rho <1$) converges to zero exponentially fast for $j\to \infty$ due to the strict minimum-phase condition. Letting $h\to \infty$ then implies the existence of a VAR(∞) representation. It follows that for such systems $A\left(z\right)$ converges absolutely for $\left|z\right|<{\rho }^{-1}$ where $1<{\rho }^{-1}$.

From the autoregressive representation the VECM representation can be obtained:

$$a\left(z\right)=I_p+\sum _{j=1}^{\infty }{A}_j{z}^j=I_p-\sum _{j=1}^{\infty }C{\overline{A}}^{j-1}B{z}^j={(1-z)}^2{I}_p-\Phi z-\Psi z(1-z)-{(1-z)}^2\sum _{j=1}^{\infty }{\Pi }_j{z}^j$$

where $\overline{A}=A-BC$ such that

$${\Delta }^2{y}_t=\Phi y_{t-1}+\Psi \Delta y_{t-1}+\sum _{j=1}^{\infty }{\Pi }_j{\Delta }^2{y}_{t-j}+{\epsilon }_t.$$

A comparison of power series coefficients provides the identities:

$$\begin{array}{ll}\hfill \Phi & =-I_p+C{(I-\overline{A})}^{-1}B,\hfill \\ \hfill \Psi & =-I_p-C{(I-\overline{A})}^{-2}\overline{A}B,\hfill \\ \hfill {\Pi }_j& =\left[C{\overline{A}}^2{(I-\overline{A})}^{-2}\right]{\overline{A}}^{j-1}B=D{\overline{A}}^{j-1}B,j=1,2,\dots \hfill \end{array}$$

It follows that the coefficients ${\Pi }_j,j=1,2,\dots$ form the impulse response of a rational transfer function of order smaller or equal to n. If $\overline{A}$ is nonsingular then the order equals n and the system $(\overline{A},B,D)$ is minimal. Furthermore, it follows that for arbitrary $\Phi$ and $\Psi$ the transfer function

$$a\left(z\right)={(1-z)}^2{I}_p-\Phi z-\Psi z(1-z)-{(1-z)}^{2}zD{(I-z\overline{A})}^{-1}B$$

is a rational transfer function with the additional property that

$$a\left(1\right)=-\Phi =-\alpha {\beta }^{\prime },{\overline{\alpha }}_{\perp }^{\prime }\dot{a}\left(1\right){\overline{\beta }}_{\perp }={\overline{\alpha }}_{\perp }^{\prime }(-\Phi +\Psi ){\overline{\beta }}_{\perp }=-{\overline{\alpha }}_{\perp }^{\prime }C{(I-\overline{A})}^{-2}B{\overline{\beta }}_{\perp }=\zeta {\eta }^{\prime }.$$

Consequently $\Phi$ and $\Psi$ determine the integration properties of processes generated using $a\left(z\right)$.

Conversely whenever the constraints

$$-I_p+C{(I-\overline{A})}^{-1}B=\alpha {\beta }^{\prime },-{\overline{\alpha }}_{\perp }^{\prime }C{(I-\overline{A})}^{-2}B{\overline{\beta }}_{\perp }=\zeta {\eta }^{\prime }$$

hold the corresponding triple $(A,B,C)$ corresponds to an I(2) process (if the eigenvalues of A are in the closed unit disc). Defining $C_{*}={\overline{\alpha }}^{\prime }C,C_{†}={\overline{\alpha }}_{\perp }^{\prime }C$ we obtain

$$-{\overline{\alpha }}^{\prime }+C_{*}{(I-\overline{A})}^{-1}B={\beta }^{\prime },-{\overline{\alpha }}_{\perp }^{\prime }+C_{†}{(I-\overline{A})}^{-1}B=0,-C_{†}{(I-\overline{A})}^{-2}B{\overline{\beta }}_{\perp }=\zeta {\eta }^{\prime }.$$

(22)

The third equation does not have a solution for fixed $B{\overline{\beta }}_{\perp },\zeta ,\eta$, if the row space of $B{\overline{\beta }}_{\perp }$ does not contain the space spanned by the rows of ${\eta }^{\prime }$. In this case row-wise projection of ${\eta }^{\prime }$ onto the space spanned by the rows of $B{\overline{\beta }}_{\perp }$ allows for (not necessarily unique) solutions in $C_{†}$. In the limit no projection is needed. Consequently for large enough T the projected matrix will have full row rank. The second equation then determines ${\overline{\alpha }}_{\perp }$ which in turn determines $\overline{\alpha }$ up to the choice of the basis such that ${\overline{\alpha }}^{\prime }=T_C{\overline{{\alpha }_o}}^{\prime }$ for some full row rank matrix ${\overline{{\alpha }_o}}^{\prime }\in {\mathbb{R}}^{r\times p},{\overline{{\alpha }_o}}^{\prime }{\overline{\alpha }}_{\perp }=0$. The first equation then can be rewritten as

$$\left[T_C,C_{*}\right]\underset{R_1}{\underbrace{\left[\begin{array}{c}-{\overline{{\alpha }_o}}^{\prime }\\ {(I-\overline{A})}^{-1}B\end{array}\right]}}={\beta }^{\prime }.$$

The second equation shows that the row space of ${(I-\overline{A})}^{-1}B$ contains the row space of ${\overline{\alpha }}_{\perp }^{\prime }$. Thus the matrix $R_1$ has full row rank. It follows that this equation has solutions.

Having obtained a solution for $C_{*},C_{†},\overline{\alpha },{\overline{\alpha }}_{\perp }$ then C is obtained from

$$C=\left[\begin{array}{cc}\alpha & {\alpha }_{\perp }\end{array}\right]\left[\begin{array}{c}{C}_{*}\\ C_{†}\end{array}\right].$$

A unique solution then can be obtained from adding the restrictions ${\Pi }_j=C{(I-\overline{A})}^{-2}{\overline{A}}^{j+1}B,j=1,2,\dots ,2n$ which for the estimates are to be solved in a least squares sense among all solutions to equations (22).

It then follows that for the true matrices $\Phi ,\Psi ,{\Pi }_j$ the only solution for given $\overline{A},B$ consists in the corresponding true C. These facts therefore can be used in order to develop an initial guess for subsequent pseudo likelihood maximization using the parameterization of I(2) processes in state space representation: Given the integer valued parameters $n,c$ and d:

• Obtain a long VAR approximation $\widehat{\Phi },\widehat{\Psi },{\widehat{\Pi }}_j,j=1,2,\dots$, including $\widehat{\Phi }=\widehat{\alpha }{\widehat{\beta }}^{\prime }$ and $\widehat{\zeta }{\widehat{\eta }}^{\prime }={\widehat{{\overline{\alpha }}_{\perp }}}^{\prime }\widehat{\Psi }\widehat{{\overline{\beta }}_{\perp }}$ using the 2SI2 approach.
• Choose the integer $f\ge n$. Use the algorithm described in Appendix F to obtain estimates $(\widehat{\overline{A}},\widehat{B},\widehat{D})$ realizing the impulse response ${\widehat{\Pi }}_j,j=1,\dots ,2f$ from the Hankel matrix with f block columns and f block rows.
• Project rows of ${\widehat{\eta }}^{\prime }$ onto the space spanned by the rows of $\widehat{B}\widehat{{\overline{\beta }}_{\perp }}$ to obtain ${\tilde{\eta }}^{\prime }$.
• Obtain a unique solution $\widehat{C}$ solving (22) such that the matrices ${\tilde{\Pi }}_j=\widehat{C}{(I-\widehat{\overline{A}})}^{-2}{\widehat{\overline{A}}}^{j+1}\widehat{B},j=1,2,\dots ,2n$ have minimal Euclidean distance to ${\widehat{\Pi }}_j,j=1,2,\dots ,2n$.
• Transform the corresponding system $(\widehat{\overline{A}}+\widehat{B}\widehat{C},\widehat{B},\widehat{C})$ to the canonical form of Bauer and Wagner (2012) to obtain the estimate $(\tilde{A},\tilde{B},\tilde{C})$.

The algorithm obtains a minimal state space system of order n in the canonical form for I(2) processes given in Bauer and Wagner (2012) and hence can be used as an initial guess for subsequent pseudo-likelihood optimization in the set $M_n(r,s)$ of all order n rational transfer functions corresponding to I(2) processes with state space unit root structure $\left(\right(0,(c,c+d)\left)\right)$.

Theorem 4 (Consistent initial guess).

Let ${\left(y_t\right)}_{t\in \mathbb{Z}}$ denote a process generated using the system $(A_0,B_0,C_0)$ according to Assumption 3 and let the system $(\tilde{A},\tilde{B},\tilde{C})$ be estimated based on the long VAR approximation with lag order chosen according to Assumption 2. Then $(\tilde{A},\tilde{B},\tilde{C})$ is a weakly consistent estimator of the data generating system $(A_0,B_0,C_0)$ in the sense that $\tilde{C}{\tilde{A}}^j\tilde{B}\stackrel{p}{\to }{C}_0{A}_0^j{B}_0,j=0,1,\dots$ and hence the corresponding transfer functions converge in pointwise topology.

The proof of this theorem can be found in Appendix E.

6. Conclusions

In this paper the theory on long VAR approximation of general linear dynamical processes is extended to the case of I(2) processes. We find that we need slightly narrower upper and lower bounds in the approximations. The tighter bounds are not needed for all results and appear not very restrictive for applications.

The main results are completely analogous to the I(1) case: The asymptotics in many respects is identical to the finite order VAR case. Asymptotic distributions for the coefficients to non- stationary variables are the same as in the finite order VAR case. This holds true both for unrestricted OLS estimates as well as the 2SI2 approach in the Johansen framework. Tests on cointegrating ranks show identical asymptotic distributions under the null as in the finite order VAR case and hence do not require other tables. In this respect the main conclusion is that the usual procedure of estimating the lag order in the first step and then applying the Johansen procedure for estimated lag order is justified also for processes generated from a VAR(∞) that is approximated with a choice of the lag order lying within the prescribed bounds.

Additionally in the VARMA case the long VAR approximation can be used in order to derive consistent initial guesses that can be used in subsequent pseudo likelihood estimation.

Thus the paper provides both a full extension of results that have been achieved in the I(1) case as well as a useful starting point for subsequent VARMA modeling which might be preferable in situations which require a high VAR order or show a large number of variables to be modeled, a situation where VARMA models can be more parsimonious than VAR models.