Inverse Covariance Matrix of ARMA(p,q) Time Series Models

Introduction

The ARMA(p, q) model with mean zero may be written

(Subscript[z, t]) = (Subscript[ϕ, 1]) 
                  (Subscript[z, t - 1])+ ...   + (Subscript[ϕ, p]) 
              (Subscript[z, t - p])  + 
            (Subscript[a, t]) - (Subscript[θ, 1]) (Subscript[a, t - 1]) -  ...  - (Subscript[θ, q]) 
          (Subscript[a, t - q])

where it is assumed that 𝔈 {(Subscript[a, t])} = 0, 𝔈 {(Subscript[a, t]^2)} = 1 and 𝔈 {(Subscript[a, t]) (Subscript[a, s])} = 0 when t ≠ s, where 𝔈 {•} denotes mathematical expectation.  It is also assumed that the model is stationary anc casual so that the autocovariance function defined by \(TraditionalForm\`γ\_k = 𝔈 {z\_t\ z\_\(t - k\)}\) exists for all k.  A necessary and sufficient for stationarity is that all roots of the polynomial 1 - (Subscript[ϕ, 1]) z -  ...   - (Subscript[ϕ, p]) 
            (z^p)  = 0 lie outside the unit circle in the complex plane.

More conveniently, in operator form we have, ϕ(B) (Subscript[z, t]) = θ(
            B) (Subscript[a, t]), where B is the backshift operator and ϕ(
          B) = 1 - (Subscript[ϕ, 1]) -  ...   - (Subscript[ϕ, p]) (B^p) and θ(
          B)  =  1 - (Subscript[θ, 1]) -  ...   - (Subscript[θ, q]) (B^q).  The condition for a causal and stationary process is simply that all roots of ϕ(B) = 0 lie outside the unit circle.  Similarly the condition for an invertible process is that all roots of θ(B) = 0 lie outside the unit circle.  The model is identifiable if ϕ(B) and θ(B) have no roots in common.

Multiplying the above model equation for (Subscript[z, t])  by (Subscript[z, t - k]) and taking expectations,

\(TraditionalForm\`\(γ\_k - \(ϕ\_1\) γ\_\(k - 1\) -  ... \)\  - \(ϕ\_p\) γ\_\(k - p\) = \(\(γ\_\(z, a\)\)(
                k) - \(θ\_1\) \(\(γ\_\(z, a\)\)(
                  k - 1)\) - \  ... \)\  - \(θ\_q\) \(\(γ\_\(z, 
                  a\)\)(k - q)\)\),

where (Subscript[γ, z, a])(
          k) = 𝔈 {(Subscript[z, t - k]) (Subscript[a, t])} is the cross-covariance function between (Subscript[z, t]) and (Subscript[a, t]).  If we take k = 0, 1,   ...  , r where r ≥ max(p, q) and make use of the symmetry (Subscript[γ, k]) = (Subscript[γ, -k]) then the right-hand side contains unknowns \(TraditionalForm\`γ\_0, \ γ\_1, \  ... \ , γ\_r\) and if we consider (Subscript[ϕ, 1]),  (Subscript[ϕ, 2]),   ...  ,  (Subscript[ϕ, p]) as fixed then we could solve this linear system of equations given that we knew the right-hand side.  The only practical difficulty is that with numeric packages for solving systems of linear equations the equations must be in the form A x = b, where A is a specified (r + 1)-by-(r + 1) matrix.  The purpose of my previous algorithm (McLeod, 1975) was to give the form of the matrix A so that a standard numeric approach could be used.

Given a random sequence, (Subscript[z, 1]),   ...  ,  (Subscript[z, n]), of length n generated by an ARMA(p, q) model, the covariance matrix is given by

\(TraditionalForm\`Γ\  = \ \((\ γ\_\(i - j\)\ )\)\_\(n×n\)\)

We implement the algorithm of Zinde-Walsh (1988) for determining the exact inverse of (Γ^(-1)).   For symbolic computation it is simply not possible, in general, to use direct computation with Mathematica's Inverse unless n is very small.  For example in the ARMA(1,1) case, given in the first example below, with n = 100 our Mathematica implementation of the Zinde-Walsh algorithm took 18.5 seconds but the direct computation using Inverse was not possible on my computer.  Even n = 25 did not return any result after more than 10 hours.

The matrix (Γ^(-1)) plays an important role in many time series applications such as exact maximum likelihood estimation in ARMA models as well as in regression with ARMA disturbances.  For numerical computation the Trench algorithm (Trench, 1964; Golub and Van Loan, 1983, Algorithm 5.7-3) also provides an efficient method for computing (Γ^(-1)) and it has the advantage of being more general.  In our forthcoming article McLeod and Zinovi (2005) we provide a C code implementation of the Trench algorithm with a MathLink interface and we discuss in more detail the time series applications.  The Zinovi-Walsh exact algorithm was useful in our derivation of explicit formula for the power function in intervention analysis given in McLeod and Vingilis (2005).

References

Golub, G. and Van Loan (1983).  Matrix Computations.  John Hoptkins University Press: Baltimore.

McLeod, A.I. (1975),  Derivation of the theoretical autocorrelation function of autoregressive moving-average time series, Applied Statistics 24, 255-256.

McLeod, A.I. and Vingilis, E.R. (2005), Power Computations for Intervention Analysis, Technometrics. 47, 174-180.

McLeod, A.I. and Zinovi, K. (2005).  Trench algorithm with applications.  Journal of Statistical Software (submitted).

Zinde-Walsh, V. (1988).  Some exact formulae for autoregressive moving average processes.  Econometric Theory, 4, 384-402.

Trench, W.F. (1964), An algorithm for the inversion of finite Toeplitz matrices.  J. SIAM 12, 515-522.

Examples

If you put this package in your Mathematica Applications directory then you can load it like this,

<< InverseCovarianceMatrixARMA.m

Overview of the Package

({{InverseCovarianceMatrixARMA[ϕ, θ, 
            n], \(\(compute\)\(\ \)\(symbolic\)\(\ \)\(or\)\(\ \)\(exact\)\(\ \)\(inverse\)\(\ \)\(for\)\(\ \)\(fixed\)\(\ \)\(n\)\(\ \)\)}, {Toeplitz[x], Symmetric Toeplitz matrix from vector x}, {TacvfARMA[ϕ, θ, 
            m], (Subscript[γ, 0]), (Subscript[γ, 1]),  ... , (Subscript[γ, m])}})

?InverseCovarianceMatrixARMA

?TacvfARMA

ARMA(1,1) Example

The covariance matrix of 5 successive observations is given,

m = 
        Toeplitz[
          TacvfARMA[{(Subscript[ϕ, 1])}, {(Subscript[θ, 1])}, 4]];

Using the Zinde-Walsh algorithm we compute the inverse,

Timing[mInv = 
          InverseCovarianceMatrixARMA[{(Subscript[ϕ, 1])}, {(Subscript[θ, 1])}, 
            5];]

mInv = mInv // Simplify;

Short[mInv, 5]

As a check we compute the inverse using the Mathematica Inverse function and check that it agrees with the previous result

Timing[mi = Inverse[m];]

Simplify[mi - mInv]

Finally we compute the determinant of the covariance matrix as well,

InverseCovarianceMatrixARMA[{(Subscript[ϕ, 1])}, {(Subscript[θ, 1])}, 5,  
          Determinant -> True] // First // Simplify

ARMA(1,0) Example

m = Toeplitz[TacvfARMA[{(Subscript[ϕ, 1])}, {}, 4]];

mInv = 
        InverseCovarianceMatrixARMA[{(Subscript[ϕ, 1])}, {}, 5] // 
          Simplify;

Short[mInv]

Simplify[Inverse[m] - mInv]

InverseCovarianceMatrixARMA[{(Subscript[ϕ, 1])}, {}, 5,  
          Determinant -> True] // First // Simplify

ARMA(0,1) Example

m = Toeplitz[TacvfARMA[{}, {(Subscript[θ, 1])}, 4]];

mInv = 
        InverseCovarianceMatrixARMA[{}, {(Subscript[θ, 1])}, 5] // 
          Simplify;

Short[mInv, 3]

Simplify[Inverse[m] - mInv]

InverseCovarianceMatrixARMA[{}, {(Subscript[θ, 1])}, 5,  
          Determinant -> True] // First // Simplify

ARMA(2,2) Example

m = 
        Toeplitz[
          TacvfARMA[{(Subscript[ϕ, 1]), (Subscript[ϕ, 2])}, {(Subscript[θ, 1]), (Subscript[θ, 2])}, 
            4]];

mInv = 
        InverseCovarianceMatrixARMA[{(Subscript[ϕ, 1]), (Subscript[ϕ, 2])}, {(Subscript[θ, 1]), (Subscript[θ, 2])}, 5] // Simplify;

Short[mInv, 5]

Simplify[Inverse[m] - mInv]