1. INTRODUCTION
The concept of cointegration, due to Granger (1981), has proved to be an extremely useful tool in the analysis of many economic and financial time series, and it has given rise to an extensive literature. Lucid surveys of this development can be found in Banerjee, Dolado, Galbraith, and Hendry (1993) and Hatanaka (1996). Following the seminal papers by Engle and Granger (1987) and Johansen (1991), much of the empirical and theoretical work on cointegration has been conducted in the context of vector autoregressive (AR) processes, but more recently interest has been shown in extending the ideas to more general models. Yap and Reinsel (1995) and Lütkepohl and Claessen (1997), for example, consider estimation problems associated with cointegrated autoregressive moving-average, ARMA, processes, as does Dhrymes (1998). In such treatises it is commonly supposed that results on the identification of stationary vector time series models associated with the work of E.J. Hannan and M. Deistler (Hannan, 1974, 1976; Deistler, 1983, 1985; Hannan and Deistler, 1988) can be readily extended to the analysis of unit-root nonstationary and partially nonstationary (cointegrated) processes without modification. As the results presented subsequently will show, such an assumption significantly understates the theoretical issues associated with such an extension and is only partially correct.
This paper will provide a detailed discussion of the structure of nonstationary dynamic simultaneous equations systems of the form
In equation (1.1) the vector
denotes a v component observable output process, and
, if present, is a u component vector of observable exogenous input variables. The v × v and v × u matrix operators A(z) = A0 + A1 z1 + ··· + Ap zp and B(z) = B0 + B1 z1 + ··· + Bp zp in the unit-delay or lag operator
, namely,
, determine the basic evolutionary properties of yt, and the stochastic disturbance,
, which is unobserved, determines how chance or random influences enter the system. The endogenous process yt is assumed to evolve over the time period t = 1,…,T, according to the specification given in (1.1) starting from initial values given by yt and xt for t = 1 − p,…,0.
With economic and financial phenomena it will rarely if ever be appropriate to think of the process as having evolved unchanged from the infinite past. Conditioning on initial values, which is what the current paradigm implies and which corresponds to common current practice in the analysis of nonstationary time series, is therefore only natural. Thus we are faced with the task of analyzing a discrete time, time invariant, and causal dynamic system where time, following a finite sequence of initial values, is explicitly confined to the positive integers.
It will be assumed that ξt for t = 1,…,T behaves like a realization of a full rank, zero mean, stationary process ξt with covariance E [ξtξt+τ′] = Γξ(τ) = Γξ(−τ)′, τ = 0,±1,±2,…, Γξ(τ) = 0 for |τ| > p. It is well known (see Hannan, 1971, Thm. 10′ and the associated discussion) that for such a stationary process there exists a sequence of zero mean, uncorrelated random variables εt, defined on the same probability space as ξt, such that
where E [εtεt′] = Σε is positive definite and, without loss of generality, the v × v matrix operator M(z) = M0 + M1 z1 + ··· + Mp zp satisfies det M(z) ≠ 0, |z| < 1. If ξt is Gaussian then the autocovariances uniquely determine the structure of the process. Motivated by these features, and employing a slight abuse of notation, we will assume that there exists a sequence of zero mean, serially uncorrelated random variables εt, t = 1 − p,…,0,1,…,T with covariance E [εtεt′] = Σε such that
for t = 1,…,T. Expressed in the form
equation (1.1) gives us an autoregressive moving-average model with exogenous variables, commonly referred to as an ARMAX system.
An ARMAX process of the type described previously clearly violates the standard conditions for yt to be stationary. More importantly, we wish to explicitly allow for unit-root, partially nonstationary behavior. Suppose that the individual series yi,t, i = 1,…,v, are asymptotically stationary after first differencing, that is, [utri ]yt = yt − yt−1, t = 1,…,T, is I(0).1
A process is said to be asymptotically stationary, denoted I(0), if it admits a representation as in (1.1) or (1.2) where A(z) is stable (see Section 2) and the exogenous input is also I(0).
A specification often assumed in the analysis of ARMAX systems is the reduced form simply identified ARMAX(pa,pb,pm) structure in which the normalization A0 = M0 = Iv is imposed and the degrees pa = δ[A(z)], pb = δ[B(z)], and pm = δ[M(z)] are prescribed where δ[A(z)] equals the degree of [A(z)], and so on. The coefficient matrices Au, 1 ≤ u ≤ pa, Bu, 1 ≤ u ≤ pb, and Mu, 1 ≤ u ≤ pm, are then free to vary subject to the identifiability conditions that [A(z) : B(z) : M(z)] is left coprime and the matrix [Apa : Bpb : Mpm] has full row rank; see Hannan (1971, 1976). This is the structure considered in Dhrymes (1998, Sect. 1.4.2) and Hsiao (1997, Sect. 3), for example. Questions relating to the choice of model structure and identification are not straightforward, however, even in the stationary case. As pointed out by Hannan, the reduced form simply identified ARMAX(pa,pb,pm) specification is actually overidentifying in the sense that it excludes particular structures from consideration. In fact, from Theorem 5.1 of Gevers (1986) we know that a simply identified ARMAX model can only represent a system in which the McMillan degree is a multiple of v and it is not canonical. As pointed out by Lütkepohl and Poskitt (1996), such features can lead to serious practical problems when investigating observed time series. Reduced form simply identified ARMAX structures will generate similar problems for the current class of processes a fortiori.
In the analysis of cointegrated systems it has proved to be advantageous for both theoretical and practical purposes to separate the long-run behavior of the system from the more transient dynamics by using the error correction, EC, form of the model due to Engle and Granger (1987). It seems sensible, therefore, to contemplate combining the advantages of the EC specification with the merits of the echelon canonical form of ARMAX processes and to consider modeling an observed time series using an error correction—autoregressive moving-average with exogenous variables—expressed in echelon form, a specification that we will henceforth denote by the acronym ECARMAXE. ECARMAXE models offer a flexible class of specifications with which to capture the dynamics of a system under study, and recent work of Lütkepohl and Claessen (1997) and Poskitt (2003) indicates that in practice it is possible to construct a more parsimonious but equally adequate representation of an observed multiple time series using an ECARMAXE model rather than a more conventional ARX model while incurring little increase in either numerical or analytic complexity. This paper provides a rigorous foundation for the empirical implementation of these models.
Section 2 of this paper outlines the general framework of the discussion and presents additional definitions and notational conventions. In Section 2 we also establish some preliminary results for a broad class of nonstationary, possibly unstable, ARMAX systems that evolve according to the specification in (1.2). Section 3 extends the Kronecker index theory to such nonstationary systems. In Section 4 this extension is adapted and modified so as to cover cointegrated processes and the canonical form of ECARMAXE specifications is established.
As a prelude to the ensuing discussion note that approaches to the Kronecker index theory that parallel the development in Akaike (1974) and that depend on Hilbert space ideas involving prediction spaces stretching back into the infinite past are commonly used in the analysis of stationary processes, as are derivations that construct the indices from the stationary autocovariance sequence. See for example Hannan and Deistler (1988), or Reinsel (1993), and the references contained therein. Such approaches are obviously not available here. As will be seen, however, we can call upon results developed in the engineering literature in the context of what is known as partial realization and the finite identifiability problem; see Hanzon (1989, Ch. 3).
Motivated by results of Phillips (1991) indicating that the best way to proceed when analyzing cointegrated systems is via maximum likelihood incorporating all prior knowledge about the presence of unit roots and the short-run dynamics, an algorithm for evaluating the Gaussian likelihood of a nonstationary ARMAX model is given in Section 2.2. Section 4.3 presents the asymptotic distribution of the exact maximum likelihood estimator for ECARMAXE models constructed using the algorithm. The paper ends in Section 5 with some brief concluding remarks. Most proofs are provided in the Appendix.
2. THE MODEL, ASSUMPTIONS AND PRELIMINARY RESULTS
Isolating the input variables xt and εt on the right-hand side of (1.2) gives
where wt′ = (−xt′ : εt′)′ and N(z) = (B(z) : M(z)). For fixed values of v and u let [A : N] denote the set of pairs [A(z) : N(z)] such that M0 = A0 and det A(0) ≠ 0. Now set δ[A(z) : N(z)] equal to the degree of [A(z) : N(z)], defined as max1≤i≤v δi [A(z) : N(z)] where δi [A(z) : N(z)], i = 1,…,v denotes the polynomial degree of the ith row of [A(z) : N(z)]. Let {[A : N]}p denote the set {[A : N] : δ[A(z) : N(z)] = p}, with p finite. For [A(z) : N(z)] ∈ {[A : N]}p define the coefficient sequence {Φ0, Φ1, Φ2,…, ΦT+p−1} via the recursive relationships
Using standard nomenclature, {Φ0, Φ1, Φ2,…, ΦT+p−1} will be referred to as the impulse response sequence. Note that by construction ∥Φi∥ < ∞, i = 0,1,…,T + p − 1, where ∥·∥ denotes the Euclidean–Frobenius norm, and the condition that det A(0) ≠ 0 implies that the power series
will be convergent for |z| < c for some c > 0. If det A(z) ≠ 0, |z| ≤ 1, then ∥Φi∥ → 0 at an exponential rate as i → ∞, c = 1, and A(z) is said to be stable.
Assumption 2.1. The series yt admits an ARMAX representation as in (1.2), or equivalently (2.1), with [A(z) : N(z)] ∈ {[A : N]}p where (i) det M(z) ≠ 0, |z| < 1, (ii) the normalization M0 = A0, det A0 ≠ 0, is imposed, and (iii) Σε > 0.
The equations in (2.2) for T → ∞ define a mapping from [A(z) : N(z)] to Φ(z) that is sufficient to determine the first- and second-order characteristics of the data generating mechanism in the stationary case by virtue of the Wold representation theorem; see Lemma 1 of Deistler (1983). The properties of a process satisfying Assumption 2.1 also depend on the homogeneous solution to equation (2.1) when viewed as a stochastic difference equation, however, and it is this feature, among others, that serves to distinguish the current situation from the stationary case.
2.1. A Realization Theorem
For completeness let us briefly review the structure of the solutions to (2.1). It is well known that the solution to a difference equation can be expressed as the sum of a particular solution and a homogeneous solution, and this is reflected in the following theorem relating the specification in (2.1) to the representation of yt in input-output final form.
THEOREM 2.1. The process yt admits an ARMAX representation of the form
with [A(z) : N(z)] ∈ {[A : N]}p and initial conditions given by (yt′ : wt′)′, t = 1 − p,…,0, if and only if yt admits a linear input-output representation
in which the conditions
, are satisfied.
The final form in (2.3) expresses yt as a function of current and past values of the input, the initial values and the coefficients. The system is therefore described as being causal and time invariant, time invariance meaning that the coefficient values are held constant over time.
Reexpressing the final form using the original partition of wt into xt and εt and employing an obvious notation for the corresponding partition of the impulse response sequence we obtain the following mean value for yt conditional on the initial values and the exogenous inputs:
where mt is calculated deterministically from the initial values. The covariance function is
Because a Gaussian process with mean (2.4) and covariance (2.5) can be readily constructed we will, for ease of exposition, assume that yt is Gaussian.
2.2. The Partial Likelihood and Identification
Suppose then that εt is a zero mean Gaussian process with variance Σε for t = 1 − p,…,0,1,…,T and let Λξ = [Γξ(t − s)], t,s = 1,…,T, denote the Tv × Tv block Toeplitz covariance matrix of ξt, t = 1,…,T, where
, and is otherwise zero. Then the function
where y1T = (y1′,…,yT′)′, y1−p0 = (y1−p′,…,y0′)′, and so on, defines the partial likelihood for the parameter vector λ′ = (β′ : σ′) where β = vec[A0 : ··· : Ap : B0 : ··· : Bp : M0 : ··· : Mp] contains the structural coefficients and σ = vech[Σε] the scale parameters (Cox, 1975). Note that the density of the endogenous variable is conditional on both the initial values and the exogenous input.
Assumption 2.2. The statistic y1−p0 is, using statistical parlance, ancillary for λ, and (2.6) defines the partial likelihood for λ where the exogenous process satisfies Assumption 2.3. There are no restrictions on σ other than those that ensure Σε > 0, and there are no joint restrictions linking the elements of σ to those of β.
Assumption 2.3. Let
denote the projection of [utri ]xt onto the space spanned by [utri ]z2−pt−1 where [utri ]zt = ([utri ]yt′,[utri ]xt)′, t = 2 − p,…,T. Then there exist constants Φ and λ, 0 < Φ < ∞ and 0 ≤ λ < 1, such that ∥Ls∥ < Φλs as s → ∞. The exogenous disturbance [utri ]ηt = [utri ]xt − 〈[utri ]x〉t|t−1 is a Gaussian process that is independent of
where Γ[utri ]η(τ), τ = 0,1,2,… is a positive definite sequence.
Given that Assumptions 2.1–2.3 hold, identification can now be defined by the standard requirement that equality between the two likelihood values f (y1T|y1−p0,x1−pT;λ) and f (y1T|y1−p0,x1−pT;λ*) for all y1T, y1−p0 and x1−pT implies that λ = λ*; otherwise there would exist parameter values λ and λ* with λ ≠ λ* such that λ and λ* are observationally equivalent. Because there are no constraints linking the elements of σ to those of β, and M0 = A0, Σε is determined uniquely by σ without further ado. Further constraints must be placed upon β, however, to select a characteristic element from within each observational equivalence class.
Before constructing our identifiable parametrization note that the partial likelihood in (2.6) corresponds to the Gaussian likelihood of the MA(p) process
. The representation of an MA(p) process is uniquely determined by its first p + 1 autocovariances (see, e.g., Hanzon, 1989, Cor. 3.1.3.2-3), and for any identifiable parametrization Γξ(τ),τ = 0,1,…, p, also determine f (y1T|y1−p0,x1−pT;λ). It follows from Hanzon (1989, Thm. 3.2.2-8(i)) that T = p + 1 is the minimal number of observations necessary to ensure that for all λ ≠ λ* we have f (y1T|y1−p0,x1−pT;λ) ≠ f (y1T|y1−p0,x1−pT;λ*).
The partial likelihood can be readily computed in practice using the following algorithm.
PROPOSITION 2.1. Suppose that yt admits an ARMAX representation satisfying Assumptions 2.1–2.3. Then for t = 2,…,T and q = min{t − 1,p} set
and for j = q − 1,…,1,
where P〈t,j〉 = M〈t,j〉 Σt−j|t−j−11/2 and Σt|t−1 = Σt|t−11/2(Σt|t−11/2)′,
with initial value Σ1|01/2 = M0−1Γξ(0)1/2. Then the partial log-likelihood function
where
The recursive calculations given in (2.7)–(2.9) are derived from the finite span Wiener–Hopf equations due to Rissanen and Barbosa (1969) applied to the MA(p) process
. They are constructed from the Cholesky factorization of the covariance matrix of (ξ1′,…, ξT′)′ via a Gram–Schmidt orthonormalization based on ξ〈t|t−1〉, the projection of ξt onto the space spanned by ξt−1,…, ξ1. This gives a square-root, orthonormal version of what is often referred to as the innovations algorithm. Proposition 2.1 follows by noting that
for t = 1,…,T. The detailed steps of the argument, which can be deduced by consulting the manipulations presented in Rissanen and Barbosa (1969), are omitted.
2Observe that the possibility that M(z) is noninvertible is not ruled out. For such a process the minimum mean-squared error predictor of ξt based on ξt−1,…, ξ1, ξ〈t|t−1〉, is well defined and can be evaluated using the algorithm in Proposition 2.1. Because the algorithm is structured in terms of the covariances E [ξtξt+τ′] = Γξ(τ) it is independent of any assumptions concerning the invertibility of M(z). See Hannan (1974, Ch. III.2) for a discussion of the synthesis of ξ〈t|t−1〉 when det(M(z)) = 0 on the unit circle.
3. THE KRONECKER INDEX THEORY FOR NONSTATIONARY ARMAX PROCESSES
Consider recouching Theorem 2.1 in terms of a sequence of block Hankel matrices of finite dimension derived from the input-output representation. To this end, set Kτ = [Φτ : mτ−p+1]τ = 0,1,…,T + p − 1, and define HR,T to be the Rv × (T + p − R)(v + u + 1) block Hankel matrix with Ki+j−1 in the (i,j)th (v × (v + u + 1)) block, i = 1,…, R, j = 1,…,T + p − R. That is,
Let hR,T(i,r) denote row (i − 1)v + r of HR,T, i = 1,…, R, r = 1,…,v. From the Hankel structure of HR,T it follows that if hR,T(i,r) lies in the linear span of hR,T(i1,r1),…,hR,T(iL,rL) where ij < i, j = 1,…, L, then row hR+1,T(i + 1,r) of HR+1,T lies in the linear span of hR+1,T(i1 + 1,r1),…,hR+1,T(iL + 1,rL). Thus the block Hankel matrix sequence HR,T, R = 1,…,T + p − 1 exhibits similar linear dependence properties to those found in the infinite-dimensional block Hankel matrix conventionally analyzed in the stationary case (cf., e.g., Hannan and Deistler, 1988, exp. (2.3.5); Reinsel, 1993, exp. (3.2)). Arguments that parallel those employed in the stationary case can therefore be used to establish corresponding results.
3.1. The Echelon Form
Our interest centers on how the properties of HR,T, R = 1,…,T + p − 1, can be used to determine the structure of the input-output system. If we define the rank of the sequence as max1≤R≤T+p−1[ρ(HR,T)] where ρ(HR,T) denotes the rank of HR,T then the following corollary to Theorem 2.1 indicates that the rank provides a partial characterization of the system equivalent to the McMillan degree.
COROLLARY 3.1. Let v and p be fixed; then for each T > vp the following holds: The process yt admits an ARMAX representation of the form
with initial conditions given by (yt′ : wt′)′, t = 1 − p,…,0, if and only if max1≤R≤T+p−1[ρ(HR,T)] = ρ(Hp,T) ≤ vp.
The characterization is completed by selecting a basis for the row space of Hp+1,T, thereby obtaining a unique parameterization for the operator pair [A(z) : N(z)] in which
(i) arc,0 = nrc,0, r,c = 1,…,v,
(ii) arr(z) = 1 + arr,1 z + ··· + arr,nr znr, arc(z) = arc,nr−nrc+1 znr−nrc+1 + ··· + arc,nr znr, r,c = 1,…,v, and
(iii) nrc(z) = nrc,0 + nrc,1 z + ··· + nrc,nr znr, r = 1,…,v, c = 1,…,u + v,
where
A pair [A(z) : N(z)] satisfying (i)–(iii) is said to be in echelon form, and the nonnegative integers ni, i = 1,…,v are called the Kronecker indices. Such a pair defines a canonical structure called an echelon canonical form.
THEOREM 3.1. Assume that the process yt admits an ARMAX representation of the form
with initial conditions given by (yt′ : wt′)′, t = 1 − p,…,0, where [A(z) : N(z)] ∈ {[A : N]}p. Then for all T > vp there exists a unique, observationally equivalent, representation of yt in terms of an operator pair [A(z) : N(z)] ∈ {[A : N]}p where [A(z) : N(z)] is in echelon form with Kronecker indices ni ≤ p, i = 1,…,v. Moreover, under Assumptions 2.1–2.3 the echelon canonical form parametrization is identifiable via the partial likelihood.
Because by assumption p is finite ρ(Hp,T) = n1 + ··· + nv ≤ vp is bounded, and the echelon form depends on fixed, finite values of the Kronecker indices, as it does in the stationary case. It is clear, however, that the values nr, r = 1,…,v, are not invariant with respect to a reordering of the elements of yt, for if P denotes an arbitrary permutation matrix then
and the linear dependences in the Hankel matrix sequence that previously generated the Kronecker indices may no longer hold for the permuted sequence
R = 1,…,T + p − 1. To this extent the echelon canonical form is only unique modulo such rotations. Following an argument that exactly parallels the development of Gevers (1986, pp. 1750–1751), however, we can establish the following version of Gevers' Lemma 2.4, which yields a unique invariant form.
LEMMA 3.1. The variables in yt = (y1t,…, yvt)′ can be permuted such that the Kronecker indices of (yr(1)t,…, yr(v)t)′ are arranged in descending order, nr(1) ≥ nr(2) ≥ ··· ≥ nr(v), where r(j),j = 1,…,v, denotes a permutation of 1,…,v that induces the ordering. The nr(j),j = 1,…,v, are unique and are referred to as the Kronecker invariants.
Thus if P now denotes a permutation matrix such that P(1,…,v)′ = (r(1),…, r(v))′ then Pyt = (yr(1)t,…, yr(v)t)′ has an echelon form ARMAX representation with Kronecker indices (nr(1),…, nr(v)) equal to the Kronecker invariants. When expressed in terms of the Kronecker invariants the coefficient matrix A0 = M0 is lower triangular, the representation of the system is canonical, and the ordered variables yr(j)t, j = 1,…,v possess unique characterizations.
Example i
Suppose that yt is a v component process generated by the reduced form structure
with initial values y0 and ε0. Then Φ0 = I, Φ1 = M − A = D, Φj = (−A)j−1D, j = 2,3,…,T, and mt = (−A)td, t = 0,1,…,T, where d = y0 − ε0. Hence, from Theorem 2.1,
as can be verified by a direct sequence of successive substitutions, and
From HR,T, R = 1,…,T, it is clear that max1≤R≤T+p−1[ρ(HR,T)] = ρ(H1,T) ≤ v and the Kronecker indices nr ≤ 1, r = 1,…,v. Assume that ρ(H1,T) = k < v. Then the Kronecker invariants are nr(j) = 1, j = 1,…, k, nr(j) = 0, j = k + 1,…,v, and the echelon form of (yr(1)t,…, yr(k)t,yr(k+1)t,… yr(v)t)′ is given by
with M0 = A0 and M1 = D + A1 where * indicates a block of the (k : (v − k))′ × (k : (v − k)) partitioned matrix whose entries are not restricted by the canonical form.
3.2. Special Features
Although the current situation shares several features in common with the stationary case certain critical differences do arise. In particular, in the stationary case boundedness of the rank of the infinite block Hankel matrix [Φ(i − j + 1)]{i,j=1,…,∞} is associated with rationality and leads to the conclusion that the canonical representation will be coprime; see Hannan and Deistler (1988) for a self-contained discussion of rational transfer functions, coprimeness, and other related issues. For a nonstationary process of the type being considered here the Hankel matrices in the sequence HR,T, R = 1,…,T + p − 1, all have finite size and hence finite rank, but the echelon form need not be coprime.
Example i′
Observe that the canonical structure given in Example i is applicable whatever the values of A and M. Thus, if A = M = −I and d = y0 − ε0 ≠ 0 where, without loss of generality, d1 ≠ 0, then max1≤R≤T+p−1[ρ(HR,T)] = 1, and a simple calculation shows that
and A1 = M1 has first row (−1,0,…,0) and is otherwise zero. Expressed in terms of the individual components we have y1,t − y1,(t−1) = ε1,t − ε1,(t−1) with initial values y1,0 and ε1,0, and yi,t − εi,t = (di /d1)(y1,t − ε1,t) for i = 2,…,v, t = 1,…,T. Solving the echelon form leads to the representation yi,t − di = εi,t, i = 1,…,v, t = 1,…,T. Note that the constants di = yi,0 − εi,0, i = 1,…,v, are definitive, whereas the stationary solutions to yt − yt−1 = εt − εt−1, for
are given by yt = d* + εt where d* is arbitrary.
This example illustrates that in a nonstationary world with fixed starting points coprimeness is not a generic property of the canonical form. To examine the non-coprime situation in a little more depth and assess the practical implications, let us extend the preceding example.
Example ii
Presume that yt is observed for t = 0,1,…,T and it is known that nr = 1, r = 1,…,v, but the coefficients A and M in the echelon form
are unknown. Suppose also that there exist first-stage estimates t of εt for t = 0,1,…,T, such that
where HT is an increasing sequence of values such that HT /T → 0 as T → ∞. At least two versions of such estimates are available. The first is given by the residuals from an AR of order hT = [c(log T)a], c > 0, a > 1, fitted to yt, t = 1,…,T, supposing the initial values yt, t = −hT + 1,…,0, are available. Proposition 3.1 of Poskitt (2003) tells us that HT = O(log T) at most for this estimate. The second possibility is to construct residuals using instrumental variable estimates of A and M (see, among others, Yap and Reinsel, 1995). Substituting t−1 for εt−1 in yt = −Ayt−1 + Mεt−1 + εt it can be seen that second-stage least-squares estimates of A and M can be obtained by solving the normal equations
Now assume that A = M = C where the eigenvalues of C, λi(C), i = 1,…,v, say, belong to the closed unit disc
. Examining the components of (3.1) we find that
where
. The second and third terms are O(l2(T)) with probability one because the process ut = vec(εtd′(−C′)t) = ((−C)td × Iv)εt is a martingale difference sequence and therefore
. By the same arguments
Now, trivially,
The second term on the right-hand side is bounded in norm by
and the first term
Thus we can conclude that
. Similarly
and the equality
implies that
equals Σε + O(l2(T)) + O({HT /T}1/2) with probability one as T → ∞.
If maxi=1,…,v|λi(C)| < 1, then for any d and all t = 1,2,3,… the norm ∥(−C)td∥ < Cλt for some C < ∞ and λ < 1 and
. Thus RT → 0 as T → ∞, and the large-sample behavior of [T : T] will reflect that yt is asymptotically stationary. If C has at least one eigenvalue on the unit circle, however, then there is no guarantee that RT is not divergent. If, as T increases, the magnitude of RT dominates the O(l2(T)) and O({HT /T}1/2) remainder terms that appear in the sums of squares and cross-products that appear in (3.1), then from Lemma A.2 of Poskitt (2000) it follows that
Using standard formulas for partitioned inversion we can therefore conclude that [T : T] will converge to the value [C : C] as the time horizon increases. The consistency observed in this situation stems from the fact that mt = (−C)td plays a nonnegligible role in determining the evolution and structure of the process and feeds information about the parameters through to the observed statistics, information that is inevitably lost in the asymptotically stationary case because mt converges to zero at an exponential rate.
3.3. The Inverse Echelon Form
Consider interchanging the roles of yt and εt in (1.2). Making the replacements yt [map ] εt and A(z) [map ] M(z) on the left-hand side, and wt [map ] zt = (yt′ : xt′)′ and N(z) [map ] [A(z) : B(z)] on the right, we obtain
with initial values (εt′ : zt′)′, t = 1 − p,…,0. The arguments employed previously can now be repeated to produce the following parallel to Theorem 3.1.
THEOREM 3.2. Assume that the process yt admits an ARMAX representation of the form
with initial values (εt′ : yt′ : xt′)′, t = 1 − p,…,0 where [A(z) : B(z) : M(z)] satisfies the conditions
(i′) arc,0 = mrc,0, r,c = 1,…, v,
(ii′) mrr(z) = 1 + mrr,1 z + ··· + mrr,nr znr, mrc(z) = mrc,nr−nrc+1 znr−nrc+1 + ··· + mrc,nr znr r,c = 1,…,v,
(iii′) arc(z) = arc,0 + arc,1 z + ··· + arc,nr znr r,c = 1,…, v and brc(z) = brc,0 + brc,1 z + ··· + brc,nr znr r = 1,…,v, c = 1,…,u.
Then the parametrization is canonical, and under Assumptions 2.1–2.3 λ is identifiable via the partial likelihood.
Note that the canonical form in (i′)–(iii′) employs the same normalization as previously, namely, that A0 = M0 with unit diagonal elements, and δr [A(z) : B(z) : M(z)] = nr, r = 1,…, s, as before, but additional exclusion constraints are placed on the elements of M(z), rather than A(z) as in (i)–(iii). Those elements not so restricted are freely varying.3
The term freely varying is used here to indicate that the echelon form imposes no restrictions on the value of the parameter. Stability and minimum phase conditions will, of course, impose separate constraints that will limit the admissible parameter values.
4. COINTEGRATION AND THE ERROR CORRECTION ECHELON FORM
We now specialize the results of the previous sections to unit-root nonstationary cointegrated systems and investigate the consequences of the identification conditions presented earlier for the analysis of partially nonstationary ARMAX structures.
4.1. Error Correction and the Kronecker Indices
To facilitate discussion in a general framework let zt = (yt′ : xt′)′, as before, and consider embedding equation (1.2) in the ARMA system
with initial conditions (zt′,et′)′, t = 1 − p,…,0, where
, and et = (εt′ : ηt′)′ is an s = v + u component white noise process with covariance matrix Σ(ε,η). Assume that zt is an I(1) process and has cointegrating rank ϱamp; ≥ 1. To isolate the integrated components of the process we can apply the following variant of the Beveridge–Nelson decomposition.
PROPOSITION 4.1. Let K(z) = [sum ]j≥0 Kj zj. Then K(z) = K(1)z + (1 − z)L(z) where L(z) = [sum ]j≥0 Lj zj, L0 = K0, Lj = −[sum ]i≥j+1 Ki, j = 1, …. Furthermore, if [sum ]j≥0 jp∥Kj∥ ≤ ∞ then [sum ]j≥0 jp−1∥Lj∥ ≤ ∞.
The result is obtained by writing the identity K(z) ≡ K(1)z + [sum ]j≥0 Kj(zj − z) and then rearranging terms using the telescoping sum
. See Phillips and Solo (1992), who provide an interesting example of the use of this decomposition in a rather different context. Applied to (4.1) the proposition leads to the EC representation
with the same initial conditions where
with
. The only term that involves potentially integrated variables in levels is Πzt−1, and it is the coefficient matrix Π that summarizes the cointegrating relations.
Consider now the identification of (4.2). Suppose that the original ARMA system in (4.1) is expressed in inverse echelon canonical form and let αp = vec[Ψ0,…, Ψp, Θ0,…, Θp]. Then conditions (i′)–(iii′) of Theorem 3.2 can be expressed in the form of the imposition of linear constraints R{n1,…, ns}αp = r{n1,…, ns} where R{n1,…, ns}, d × s2p, and r{n1,…, ns}, d × 1, d = 2s2(p + 1) − ([sum ][sum ]i<j{min(ni,nj) + min(ni,nj + 1)} + (s + 1)[sum ]i ni), are known. In addition, the EC representation is obtained via the parametric transformation
where S is the s(p + 1) × s(p + 1) nonsingular matrix
This implies that
and hence that
where
Because the mapping from
is one to one it follows that the original ARMA system in (4.1) satisfies the constraints R{n1,…, ns}αp = r{n1,…, ns} if and only if the EC representation in (4.2) satisfies
, and the identification of one implies that of the other.
Examining these restrictions in detail we find that for any nonzero Kronecker index nr < p,r = 1,…, s, ψrc,j is freely varying for j = 1,…, nr but ψrc,j = 0,j = nr + 1,…, p, c = 1,…, s. It follows that
is freely varying for j = 1,…, nr − 1 but
. If nr = p then the ψrc,j,j = 1,…, nr,c = 1,…, s, are all freely varying, and the same is true for
. Thus the restrictions
incorporate the conditions that
. We also have that
is subject to the same restrictions in
as in R{n1,…, ns}αp = r{n1,…, ns}. Thus Θ(z) presents as the MA operator of an inverse echelon canonical form with Kronecker indices nr,r = 1,…, s in both (4.1) and (4.2).
There are s2 fewer parameters in
than in the original pair [Ψ(z) : Θ(z)], and the degrees of freedom so released are taken up by the elements of Π. Thus far Π remains unrestricted. By assumption, however, det Ψ(z) = ψ(z)(1 − z)ζ where ψ(z) ≠ 0, |z| ≤ 1, and 1 ≤ ζ < s, and it is well known (see, among others, Yap and Reinsel, 1995, Sect. 2) that the rank of Π = Ψ(1) equals s − ζ. Consequently, additional restrictions must be applied over and above those given by
if we require a particular cointegrating rank ϱamp; = s − ζ to hold.
To link the identification of Π to that of the original specification note that if rank(Π) = ϱamp; then Π = FG′, where F and G are (s × ϱamp;) matrices with full column rank. To ensure a one-to-one correspondence between Π and its reduced rank factorization we follow a standard procedure. Because the columns of G are linearly independent there exists a nonsingular (ϱamp; × ϱamp;) matrix E, constructed from a sequence of elementary column transformations, such that Γ = GE is column equivalent to G and in reduced column-echelon form. Postmultiplying F by (E′)−1 leads to an identified pair [ϒ : Γ] with ϒ = F(E′)−1. The condition that rank(Π) = ϱamp; is thereby obtained by imposing an additional ϱamp;2 constraints on Γ and leaving the elements of ϒ unconstrained.
Let (4.2) denote an EC representation in which the operator pair
are in inverse echelon canonical form but with the added restrictions
, where nr > 0, r = 1,…, s, imposed, and Π = ϒΓ′ where ϒ and Γ are (s × ϱamp;) matrices with full column rank and Γ is in reduced column-echelon form. Then the upshot of the preceding argument is that the structure in (4.2) is identified and is equivalent to an inverse echelon canonical form ARMA representation (4.1) in which the cointegrating rank ϱamp;, 1 ≤ ϱamp; ≤ s − 1, has been imposed.
If any of the Kronecker indices are zero, nq = 0, 1 ≤ q ≤ s, say, then the qth row of
and Π are equal, and hence the qth row of (4.1) is identical to that of (4.2). From the structure of the echelon form this implies that the variable zq,t can be expressed as a contemporaneous linear combination of the remaining variables in the system, and the innovations, and therefore it will inherit all its dynamics from these other variables. Hence zq,t must be either I(0) or cointegrated with some of the other variables so that ϱamp; ≥ 1. More generally, assume that the variables have been permuted such that the system is represented in terms of the Kronecker invariants. If nr(s) = ··· = nr(s−q+1) = 0 and nr(j) ≥ 1,j = 1,…, s − q, then some relatively straightforward manipulations indicate that
where Π11, Π10, and Π01 are ((s − q) × (s − q)),((s − q) × q), and (q × (s − q)) coefficient matrices, respectively. Obviously the rank of this matrix is at least q, so q ≤ ϱamp; = rank(Π) ≤ s − 1.
There are two basic conclusions to be drawn from the preceding analysis: First, the conditions for identifying the short-run dynamics in an EC inverse echelon canonical form, conditions (EC i′)–(EC iii′) of Theorem 4.1, which follows, can be separated from those that identify the long-run relationships, (EC iv′)–(EC v′) of Theorem 4.1. Moreover, the long-run relationships can be present among any of the variables in zt. Second, specifying an original dynamic system in which one or more of the Kronecker indices are zero amounts to a presumption that cointegration is present, the static equations corresponding to the zero indices representing the associated long-run equilibrium relationships assumed. On the whole it seems unlikely that preconditions of the latter type will be imposed a priori in a pure time series setting as the data are usually left to speak for themselves and such precise information is rarely available. On the other hand, for some relatively simple economic models the short-run dynamics and long-run relationships can be written down directly; see, for example, Wickens and Breusch (1988).
More generally, if a dynamic simultaneous equations perspective is taken then conclusions about the cointegrating structure are often implicit in the model formulation and the assumptions made. In the notation of the current section the dynamic structural equation (1.2) corresponds to the specialization Ψ11(z) = A(z), Ψ12(z) = B(z), Θ11(z) = M(z), and Θ12(z) = 0 where Ψij(z) and Θij(z), i,j = 1,2 denote partitions of the s × s operators Ψ(z) and Θ(z) into the first v and last u rows and columns. Adding the restrictions that Ψ21(z) = Θ21(z) = 0 and Σ(ε,η) = diag(Σε, Ση) amounts to the imposition of the assumption that xt is strictly exogenous; see Engle, Hendry, and Richard (1983). If we also suppose that det Ψ11(z) ≠ 0, |z| ≤ 1, and Ψ22(z) = D(z)(1 − z) where D(z) is stable then
and ϱamp; ≤ v. But following the argument employed by Hsiao (1997, p. 653) we find that the possibility that ϱamp; < v is ruled out, for otherwise there would exist a nonzero vector v such that v′[A(1) : B(1)] = 0, contradicting the assumption that Ψ11(z) = A(z) is stable. Thus ϱamp; = v,
and, following Wickens (1996), the variables in xt may be regarded as the common trends of Stock and Watson (1988).
The strength of the previous conclusion, which corresponds to that drawn by Hsiao (1997), obviously depends critically on the stringency of the assumptions. If the previous coefficient conditions are relaxed by replacing Ψ21(z) = 0 by Ψ21,0 = 0 so that Ψ21(z) = Ψ21,1 z1 + ··· + Ψ21,p zp, and the restrictions that Ψ11(z) is stable and Ψ22(z) = D(z)(1 − z), det D(z) ≠ 0, |z| ≤ 1, are removed, then a range of different cointegrating structures is possible. In particular, the dichotomy between the endogenous and exogenous variables and that between the short-run dynamics and long-run equilibrium relationships need no longer coincide. Such a system, in which xt exhibits feedback but is weakly exogenous, might be appropriate where xt contains policy instruments determined via a partial adjustment process or control variates determined via a linear-quadratic control rule, for example. For a discussion of causality and feedback in the context of nonstationary processes and exogeneity see Hosoya (1977) and Geweke (1984), respectively.
4.2. The Error Correction Echelon Form
In light of the previous discussion we will adopt the following assumption.
Assumption 4.1. The series yt is an I(1) process that admits an ARMAX representation as in (1.2), or equivalently (2.1), with [A(z) : N(z)] ∈ {[A : N]}p where det A(z) = (1 − z)ζ d(z), ζ ≤ v, and d(z) is stable.
The most general result concerning the canonical structure of our original partially nonstationary ARMAX model is then as follows.
THEOREM 4.1. Let
denote an EC representation in which
, and Π = D0 + D1 + ··· + Dp, where D(z) = [A(z) : B(z)] and the polynomial operators in
satisfy the conditions
where nr > 0, r = 1,…, v, and Π = ϒΓ′ where
Then under Assumptions 2.1–2.3 and 4.1 the structure in (4.3) is identified and is equivalent to an inverse echelon canonical form ARMAX representation in which the cointegrating rank ϱamp; ≤ v has been imposed.
It is the system in Theorem 4.1 that has previously been christened an ECARMAXE form.
To prevent a proliferation of notation the same symbolism is employed in Theorem 4.1 for the cointegrating relationships as was adopted earlier, namely, Π = FG′ = ϒΓ′, only now ϒ = FE−′ and Γ = GE are (v × ϱamp;) and ((v + u) × ϱamp;) matrices with full column rank. The reduced column-echelon form is a mathematical artifact often employed in matrix algebra that serves here to solve the statistical identification problem. In such a matrix the first nonzero entry in any column is unity and appears below the first nonzero entry in the preceding column. All other entries in the same row as the first nonzero entry in any column are zero. If, after suitable permutation denoted by RG, the first ϱamp; rows of G are linearly independent then the reduced column-echelon form becomes
This gives the triangular structure Γ′zt−1 = [Iϱamp; : Γϱamp;′]R′zt−1 introduced by Phillips (1991). For the reduced column-echelon form Γ any arrangement of the variables is permitted, and this allows the system to be ordered according to the permutation induced by the Kronecker invariants, which need not coincide with the reordering implicit in the triangular structure.
4.3. Parameter Estimation
Given that the ECARMAXE form is identified we are now interested in estimating the unknown parameters in λ. Suppose that ϱamp; and nr(i), i = 1,…,v, are known a priori and that the variables and equation system have been ordered according to the permutation r(1),…, r(v) induced by the Kronecker invariants to give the unique invariant form. Let {λ : ϱamp;,(nr(1),…, nr(v))} denote the set of parameter values λ = (β′ : σ′)′ such that β satisfies conditions (EC i′)–(EC v′). By definition the maximum likelihood estimator
is the value of λ that maximizes LT(λ) = log{f (y1T|y1−p0,x1−pT;λ)} over {λ : ϱamp;,(nr(1),…, nr(v))},
The determination of
will necessitate the use of numerical optimization techniques, and the function evaluations required to implement such numerical methods can be completed by using the innovations algorithm of Proposition 2.1 with
replaced by
, in the definitions of ε〈t|t−1〉, t = 1,…,T.
It is of interest to observe that if M(z) = M0 then the model in equation (4.3) reduces to
where ut = M0εt, a cointegrated ARX structure with row degrees nr(i), i = 1,…,v. In this case it is relatively straightforward to verify that the recursions in (2.7)–(2.9) of Proposition 2.1 yield M〈t,j〉 = 0, 1 ≤ j ≤ q, and Σt|t−1 = Σε for all t = 1,…,T. We are thereby led to the conclusion that ε〈t|t−1〉 = εt and that the partial log-likelihood function can be reexpressed as
where Σu = M0 ΣεM0′, because M0 is lower triangular with leading diagonal equal to the identity and hence det M0 = 1. Expressions (4.4) and (4.5) coincide with those commonly considered in the analysis of cointegrated autoregressive systems, following Johansen (1991).
To investigate the statistical properties of
consider the estimator
obtained by maximizing the following approximation to the partial log-likelihood function:
where εt for t = 1,…,T are calculated from the recursion
starting from the initial values εt = 0, t = 1 − p,…,0. Then
where
is given by the solution to the score equations
. These equations can be solved numerically, leading to the iterative Newton–Raphson approximation
where, by definition,
is the ith iterate. It is
that is commonly referred to as the Gaussian estimator, and it is this estimator that is considered by Yap and Reinsel (1995) and Dhrymes (1998). We will now specialize (4.7) to the ECARMAXE form.
Rearranging equation (4.3) so as to isolate the common parameters in D0 = [A0 : B0] and M0 = A0 gives
Following Poskitt (1992) let us now rewrite the left-hand side by vectorizing each term. The first two terms become (yt − εt) + ([(yt − εt)′ : xt′] [otimes ] Iv)ζ where ζ = vec([A0 − I : B0]), and the third is
where ξp−1(z)′ = (z,z2,…, zp−1) and δ = vec(D1 : ··· : Dp−1). The fourth term becomes −(zt−1′Γ [otimes ] Iv)υ where the parameter vector υ = vec(ϒ), and the fifth term gives
, where μ = vec(M1 : ··· : Mp). In the derivation of these expressions the well-known rule vec(ABC) = (C′ [otimes ] A)vec B has been employed. Set θ = (δ′ : ζ′ : μ′ : υ)′. The vector θ contains the parameters of D(z), M(z), and ϒ not restricted to be unity. For any ECARMAXE form with
the exclusion constraints implicit in (EC i′)–(EC iii′) are simply incorporated by deleting the corresponding elements of θ. To complete the parameterization of the model let γ = vec(Γϱamp;) denote the coefficients in (RΓ)′ = [Iϱamp; : Γϱamp;′] not restricted to be zero or one in the echelon form. This gives us a freely varying parameter vector β = (θ′, γ′)′ of dimension dβ = dθ + dγ, where
, such that λ = (β′, σ′)′ ∈ {λ : ϱamp;,(nr(1),…, nr(v))}.
Treating (4.8) as an implicit function of the stochastic disturbance, the data, and β and differentiating with respect to β we find that ∂εt /∂β′ = −Wt where Wt = [Wθt : Wγt] equals the v × dβ matrix obtained by selecting the appropriate columns of
where the v × v(v + u), v × v2, v × ϱamp;v, and v × ϱamp;(v + u − ϱamp;) derivative processes Wδt, Wμt, Wυt, and Wγt are generated from
H′ = [0 : I(v+u−ϱamp;)]R′. The Gaussian iterations of (4.7) can then be expressed as
Now let
denote the iterate obtained by substituting ε〈t|t−1〉 for εt, ∂ε〈t|t−1〉 /∂β′ = −W〈t|t−1〉 for −Wt and Σt|t−1 for T in (4.9). On convergence
will yield a critical point of the partial likelihood, but to ensure that the iterates
will converge to
when T is sufficiently large the iterations must be initiated using a consistent estimator. Suppose that the iterations commence at
where T denotes a preliminary estimator chosen such that NT(T − β) = Op(1), NT = diag[T1/2Idθ : TIdγ]. Then we have the following result.
THEOREM 4.2. If Assumptions 2.1–2.3 and 4.1 hold and, in addition, the operator M(z) is invertible, then the iterates
obtained using the initial value T are asymptotically equivalent in the sense that
for all i ≥ 1 as T → ∞. Moreover, both estimators converge in distribution to βT = (θT′, γT′)′ where
where
the asymptotic information matrix for θ; the vectors
and T(γT − γ) are asymptotically mutually uncorrelated; and T(γT − γ) = T vec(Γϱamp;,T − Γϱamp;) where the components of (Γϱamp;,T − Γϱamp;) satisfy the asymptotic mixed-normality result
where
When applying Theorem 4.2 Vθ and Vγ can be consistently estimated by replacing the unknown parameters in β by their Gaussian maximum likelihood values and substituting a consistent estimate for Σε, while dropping the expectation from Vθ. Two simple estimates of Σε that can be constructed as a by-product of the output from the algorithm in Proposition 2.1 are
.
COROLLARY 4.1. Both
are consistent estimators of Σε. Moreover,
have the same limiting distribution as
.
Finally, the initial estimator T must be chosen. One such estimate can be obtained by first using instrumental variables with instruments chosen from yt−τ, xt−τ, and zt−τ−1 plus [utri ]zt−τ−1,…,[utri ]zt−τ−p, for τ > p, to calculate T, T, T, and T, respectively. From these estimates T,t, to use an obvious notation, can be generated, and then T can be evaluated using the spectral factorization method of Tunnicliffe-Wilson (1972), as suggested by Yap and Reinsel (1995, Sect. 3.3) and Dhrymes (1998, pp. 325–326). Alternatively, T can be obtained by implementing the technique proposed by Brockwell and Davis (1988) following a preliminary pass through Proposition 2.1 with Γξ(τ) replaced by
.
Example ii′
To provide some indication of the possible impact of the results outlined previously a sequence of Monte-Carlo experiments has been conducted using the data generating mechanism yt + Ayt−1 = εt + Mεt−1, t = 1,…,T, where the pair [A : M] accord with the structure considered in Example ii, namely, A = M = C. The behavior of the maximum likelihood and Gaussian estimates
constructed from a single iteration of the Newton–Raphson (scoring) algorithm initiated at [T : T] has been monitored. The behavior of
was examined by calculating the observed mean-squared error, the results of which are summarized by presenting the empirical mean and the average squared Euclidean distance between the estimate and the true value [C : C], labeled M.S.E. in the tables that follow.
In all cases v = 2, n1 = n2 = 1, and the coefficient values are given by either
In the first case, process P1, the zeros of det(I + Cz) are 1.1 ± i0.8888, and in the second, process P2, ζ1 = −1.25 and ζ2 = 1.0. The covariance matrix
and the sample sizes considered were T = 50, 100, 150, and 250. All simulation results listed here were based on 1,000 replications.
The influence of the initial conditions was controlled by setting d = (1,1)′d where the scalar d was chosen so as to make ρy,m = det(RT)/det(RT + Σε) equal to 0.1QT, QT = (log(T))1/2l2(T), for P1, and 0.1 for P2. The rationale for this follows from observing that for both processes ρy,m is the (asymptotic) coefficient of vector correlation (Hotelling, 1936) between yt and mt. Process P1 is asymptotically stationary, and multiplication by the factor QT ensures that the influence of the initial conditions does not die away too quickly, so that RT → 0 as T → ∞ but RT /l2(T) = O((log(T))1/2). For P2, of course, RT is O(1), and such rescaling is not necessary. For process P1 ρy,m = 0.1QT takes the values 0.0327, 0.0265, 0.0232, and 0.0194 for T = 50, 100, 150, and 250, respectively, indicating that fluctuations in mt account for well below 5% of the observed variation in yt for both processes.
The values presented in Table 1 are typical of those obtained using different parameterizations of process P1. The figures indicate a somewhat superior performance for
, with the observed relative
decreasing from 0.3319 when T = 50 to 0.0849 when T = 250 whereas
decreases from 1.2133 to 0.255. Equivalent results for P2 are given in Table 2. A striking feature of the figures given in this second table is the relatively poor performance of
. This presumably reflects aspects of the Gaussian approximation that work less well in the presence of a unit root and further emphasizes the benefits of determining the maximum likelihood estimate.
Experimental outcomes for process 1
Experimental outcomes for process 2
Although both P1 and P2 are in a sense pathological processes this is not obvious from the figures reported in Tables 1 and 2. Mean-squared error does not capture the distributional properties of the estimators, however. Figure 1 presents kernel density estimates of the distribution of the Mahalanobis distance
compared to theoretical χ2 densities for process P1. Process P1 is asymptotically stationary, and the specification A = M means that it will also be asymptotically unidentified. Such a lack of identification is well known to manifest itself in
having less than full rank and QT having fewer degrees of freedom than might be anticipated on the basis of conventional asymptotic theory; see, for example, Poskitt and Tremayne (1982). Figure 1 lends clear testimony to this feature. For P2 both A and M are identified and the estimates
can be readily calculated, but Theorem 4.2 is not applicable because the invertibility condition is violated. It seems reasonable to conjecture that a similar theorem could be established that would allow for the singularity present under the less restrictive condition det M(z) ≠ 0, |z| < 1, of Assumption 2.1, but that avenue will not be pursued here. See Tanaka (1996) for a discussion of the issues associated with lack of invertibility.
Empirical distribution of Mahalanobis distance. Processes P1 with T = 150, ρε = 0.8, and ρy,m = 0.1QT. Kernel density estimate κ(QT) obtained using Gaussian kernel with optimal bandwidth.
5. CONCLUDING REMARKS
This paper has filled an important gap in the identification theory of nonstationarity vector ARMAX systems by showing that the Kronecker index theory can be extended to a broad class of nonstationary ARMAX systems and that ECARMAXE models provide a canonical form for partially nonstationary (cointegrated) ARMAX processes.
In the context of ECARMAXE models the asymptotic equivalence of the Gaussian estimator
and the maximum likelihood estimator constructed using the innovations algorithm
has been established. The large-sample distribution of the maximum likelihood estimator for such models was also presented. As pointed out by a referee, such results are contingent on the true values of the cointegrating rank and the Kronecker indices being known a priori. If ϱamp; and nr, r = 1,…,v, are chosen by the practitioner and they do not coincide with the true values then the likelihood function may not be well behaved. A maximum may not exist within the parameter space considered as a supremum to the likelihood could lie at the boundary of the parameter space. Even if a maximum exists it need not be unique; several parameter points could in principle correspond to the same maximum likelihood value. For a discussion of such issues see Deistler (1983, 1985) and Hannan and Deistler (1988, Ch. 4).
The normality assumption underlying the analysis conducted in this paper is commonly adopted in the literature on cointegration. Normality does not play a key role beyond motivating the estimators, however, and it seems likely that the asymptotic properties derived could be extended to more general processes under weaker regularity conditions.
In closing it is worth emphasizing that the identification conditions for an ECARMAXE model depend on fixed, finite values of the cointegrating rank ϱamp; and the Kronecker indices, nr, r = 1,…,v, and are applicable at any sample size T > vp. This latter point is of significance for any future development of exact finite-sample distribution theory for
or, perhaps more importantly, bootstrapping methodology.
APPENDIX: Proofs
Proof of Theorem 2.1. Let
, where the impulse response sequence is as given in (2.2). By construction
for τ = 1,2,…,T + p − 1, which when evaluated at τ = t + p gives
for t = 1,…,T. Hence yP,t provides a particular solution to (1.2).
Now let mt = yt − yP,t, t = 1 − p,…,0, and set
, t = 1,…,T. Then mt defines an appropriate complementary function and yP,t + mt gives the general solution to (1.2) because by construction yt = yP,t + mt, t = 1 − p,…,0, and
Conversely, suppose that the process yt admits an input-output representation as in (2.3) and that v × v coefficient values Aj, j = 0,…, p, exist such that the impulse response sequence satisfies the difference equation
and, similarly, mt solves
. Then for t = 1,…,T
Proof of Corollary 3.1. This necessarily follows by noting that the conditions
imply that
From the nonsingularity of A0 it follows that rows pv + r, r = 1,…,v, of Hp+1,T can be expressed as linear combinations of rows (i − 1)v + r, i = 1,…, p, r = 1,…,v and hence that ρ(Hp+1,T) ≤ vp. Theorem 3.1.2.3-29(iii) of Hanzon (1989) and the linear dependence properties of HR,T, R = 1,…,T + p − 1, now imply that ρ(HR,T) ≤ vp, R = p + 1,…,T + p − 1, and thus max1≤R≤T+p−1[ρ(HR,T)] = ρ(Hp,T) ≤ vp.
To establish sufficiency suppose that max1≤R≤T+p−1[ρ(HR,T)] = ρ(Hp,T) ≤ vp. Recognizing that Hp,T−1 can be obtained from Hp,T by deleting the last block of columns, and from Hp+1,T by deleting the last block of rows, we can use Hanzon (1989, Thm. 3.1.2.3-29(iii)), because T > vp, to conclude that ρ(Hp+1,T) = ρ(Hp,T−1) ≤ vp. Thus the rows of Hp+1,T are linearly dependent, and each of the last v rows can be expressed as a linear combination of the rows that precede it. These row combinations generate a sequence of v × v coefficient values Aj, j = 0,…, p, such that
where A0 is lower triangular and nonsingular, with leading diagonal equal to the identity, and all of A0,A1,…,Ap cannot be zero. That is,
. Appeal to Theorem 2.1 completes the proof. █
Proof of Theorem 3.1. Selecting the first basis rows of Hp+1,T in natural order produces v integers nr, r = 1,…,v, such that n1 + ··· + nv = ρ(Hp+1,T) and hp+1,T(1,1),…,hp+1,T(n1,1) through to hp+1,T(1,v),…,hp+1,T(nv,v) form a basis for the rows of Hp+1,T. Expressing row nr v + r of Hp+1,T as a linear combination of its linearly independent antecedents results in an equation system analogous to (A.1) that may be solved uniquely for the coefficient values arj,nr−nrj+1,…, arj,nr, r,j = 1,…,v. Equating arj,nr−s+1 to the (r,j)th element of Anr−s+1, s = 1,…, nrj, with arr,0 = 1, r = 1,…,v, and all other elements equal to zero, yields the autoregressive operator A(z) = A0 + A1 z1 + ··· + Ap zp where p = max1≤r≤v(nr). Given A(z) the coefficients of N(z) are obtained by evaluating
. By construction the row degrees of [A(z) : N(z)] equal the Kronecker indices, that is, δr [A(z) : N(z)] = nr, r = 1,…,v. The detailed steps in the preceding argument follow those used by Hannan and Deistler (1988) in establishing their Theorem 2.5.1.
Thus, for any given ARMAX system with [A(z) : N(z)] ∈ {[A : N]}p the sequence HR,T, R = 1,…,T + p − 1, can be readily determined, and from these ρ(Hp,T), the Kronecker indices ni, i = 1,…,v, and the observationally equivalent echelon form can subsequently be constructed. Conversely, every ARMAX system such that [A(z) : N(z)] satisfies (i)–(iii) obviously defines an echelon form representation with Kronecker indices nr, r = 1,…,v. Finally, by Hanzon (1989, Thm. 3.2.2-8(i)), for λ corresponding to an echelon canonical form, λ ≠ λ* implies f (y1T|y1−p0,x1−pT;λ) ≠ f (y1T|y1−p0,x1−pT;λ*) and the parametrization is identifiable. █
Proof of Theorem 4.2. To show that
we use the principle of induction. Assume that
are Op(1). Substituting
in (4.9) and using the stochastic equicontinuity results given in Saikkonen (1996) we obtain the asymptotic representation
A corresponding expression with
replaced by
, and W〈t|t−1〉 substituted for Wt, Σt|t−1 for T, and ε〈t|t−1〉 for εt, also obtains. It follows that the probability limit of
will be zero if
and
Suppressing explicit evaluation at the point β for notational convenience, both (A.3) and (A.4) are obtained by expressing the left-hand side in terms of the differences ∇εt = εt − ε〈t|t−1〉, ∂∇εt /∂β′ = W〈t|t−1〉 − Wt = −∇Wt = −[∇Wθt : ∇Wγt], and ∇ΣTt = T − Σt|t−1.
Consider first (A.3). Expanding each term in the summation as
we can see that the left-hand side of (A.3) can be decomposed into four series. The fourth of these,
, is bounded in norm by
But T = Σε + o(1) by ergodicity, and a direct application of Lemma 1 of Rissanen and Caines (1979) tells us that Σt|t−1 converges to Σε at a geometric rate. The latter implies that there exist constants Kσ and λσ with 0 ≤ Kσ < ∞ and 0 < λσ < 1, Kσ λσT > ∥T − Σε∥, such that ∥Σt|t−1−1∥ ≤ ∥Σε−1∥ + Kσ λσt and ∥∇ΣTt∥ ≤ ∥T − Σε∥ + ∥Σε − Σt|t−1∥ ≤ 2Kσ λσt. It follows that
and therefore A.5 is majorized by
Now let [ell ]T = −2 log(T)/log(λσ). Then
because
and [ell ]T /T → 0 as T → ∞ and λσ[ell ]T = T−2 by construction. Thus we can conclude that (A.5) is op(1) and hence
.
The first two series in the decomposition of (A.3) are bounded in norm by
and the norm of the third is bounded by
Thus these three series will converge to zero in probability if we can show that
To establish (A.6) note that simple manipulation of the equality between
gives
and therefore the proof given by Rissanen and Caines (1979, pp. 312–314) that
(cf. Rissanen and Caines, 1979, eqn. (A3.4)) is directly applicable. Differentiating (A.7) with respect to β and rearranging gives
By Lemma 1 of Rissanen and Caines (1979) there exist constants Kμ and λμ, 0 ≤ Kμ < ∞, 0 < λμ < 1, such that ∥M〈t,j〉 − Mj∥ < Kμ λμt uniformly in β, and therefore the same is true for ∥∂(M〈t,j〉 − Mj)/∂β′∥. A straightforward adaptation of the argument that gives (A.8) therefore leads to the conclusion that
as T → ∞, giving (A.6) as required. Equation (A.3) now follows.
Now consider equation (A.4). First we bound the norm of the left-hand side by
The fourth term
But
Now observe that for any λ, 0 ≤ λ < 1,
where [ell ]T = −2 log(T)/log(λ). An analogous derivation also gives
and
Substituting these bounds into the preceding inequalities gives
Equivalent bounds on the order of magnitude of the first three terms are derived using the fact that both ∥∇εt∥ and ∥∇Wt′∥ converge to zero at a geometric rate in t. Confirmation of this result is obtained by noting that the conditions ∥M〈t,j〉 − Mj∥ < Kμ λμt and det(M(z)) ≠ 0, |z| ≤ 1 when applied to (A.7) imply that
for some Kε, Kμ ≤ Kε < ∞, and λε, λμ ≤ λε2 < 1, where ∥ηε∥ < ∞ bounds the influence of the initial conditions. Hence ∥∇εt∥ ≤ λεtMε where Mε < ∞ but exceeds
Using the inequality ∥∇εt∥ ≤ λεtMε we also have
, leading to the conclusion that ∥∇Wt∥ ≤ λWtMW.
Applying these geometric bounds yields the inequalities
for the first term,
for the second, and, similarly,
for the third. But
by the Cauchy–Schwarz inequality and (A.10) and
because (A.9) and (A.11) imply that
. Finally,
Thus we have established (A.4).
The induction to show that
for all i ≥ 1, in addition to the convergence in distribution to βT, is now completed by verifying that the components of
converge in distribution as stated. The proof of the latter follows along lines that parallel the developments in Yap and Reinsel (1995) and Dhrymes (1998). Recall that for an ECARMAXE model the constraints of the echelon form are incorporated by simply deleting appropriate elements of θ. This implies that corresponding rows and columns in previous expressions involving Wt are similarly removed, and the same is true of (A.12) and (A.13), which follow. Hence the arguments of Yap and Reinsel (1995, cf. Sect. 4; see in particular the comment in Sect. 4.1) and Dhrymes (1998, cf. Sect. 6.4) can be applied at this stage with very little modification. The argument proceeds by showing that
equals
and that
can be rewritten as
Both (A.12) and (A.13) are established by applying Proposition 4.1 to M(z) to give
The sums of squares and cross-products involving Wγt will therefore be dominated by the components in (M(1)−1ϒ [otimes ] zt−1′H) because the process H′zt is integrated and [utri ]Wγt and Wθt are asymptotically stationary processes. Substituting (A.12) and (A.13) into (A.12) and (A.13) and (A.2) and applying standard central limit theorems to (A.13), following Yap and Reinsel (1995, p. 256) or Dhrymes (1998, pp. 326–333), yields the required result. █
Proof of Corollary 4.1. We have already shown that ∥∇εt∥ ≤ λεtMε where Mε < ∞ and 0 < λε < 1. From the inequality ∥εtεt′ − ε〈t|t−1〉ε〈t|t−1〉′∥ ≤ ∥∇εt∥2 + 2∥εt∥·∥∇εt∥ we therefore obtain the upper bound
and from (A.10) it follows that
. Similarly, the inequality ∥∇ΣTt∥ ≤ 2Kσ λσt implies that
and hence that
. █
