A note on calculating autocovariances of periodic ARMA models

An analytically simple and tractable formula for the start-up autocovariances of periodic ARMA (PARMA) models is provided.

💡 Research Summary

The paper addresses the problem of computing the start‑up autocovariances of periodic autoregressive moving‑average (PARMA) models, a task that traditionally relies on the periodic Yule‑Walker equations. Earlier works (Li and Hui, 1988; Shao and Lund, 2004) formulate the problem as a linear system A γ = y or Γ U γ = κ, where the matrices involved are either not given in closed form or become increasingly cumbersome as the model order or period grows. Consequently, practical implementation can be difficult and computationally expensive.

In this note the author proposes a more transparent and analytically tractable approach based on a lattice representation. Consider a causal PARMA(p,q) model with period S:

  Xₙ₊ₛ = Σ_{j=0}^{p} φ₍ᵥ₎ⱼ X_{n‑j+S} + Σ_{j=0}^{q} θ₍ᵥ₎ⱼ ε_{n‑j+S}, 1 ≤ v ≤ S,

where the coefficients φ₍ᵥ₎ⱼ and θ₍ᵥ₎ⱼ may vary with the season v, and εₜ is a periodic white‑noise sequence with variance σ²ᵥ. The autocovariance function γ₍ᵥ₎(h) satisfies the periodic Yule‑Walker difference equation (2) and depends on the normalized cross‑autocovariances ψ₍ᵥ₎(k) defined by the recursive relation (3). To start the recursion one needs the (p + 1) initial values γ₍ᵥ₎(h) for h = 0,…,p and for every season v.

The author stacks these initial values into a single vector

 γ = (γ₍1₎(0),…,γ₍1₎(p), γ₍2₎(0),…,γ₍S₎(p))′

of dimension S(p + 1), and similarly defines a right‑hand‑side vector ζ whose entries are the known terms of (2). The key contribution is the explicit construction of a block‑circulant matrix Φ such that

 Φ γ = ζ.

Φ is built from smaller blocks ϕ₍ᵥ₎(h) (equation 4) and Φ₍ᵥ₎(k) (equation 5). The matrix ϕ₍ᵥ₎(h) has a triangular shape: for h ≤ p it contains the AR coefficients φ₍ᵥ₎ⱼ placed in appropriate positions, while for h ≥ p + 1 it is a zero matrix. Consequently, each Φ₍ᵥ₎(k) is a finite sum of at most p non‑zero ϕ‑blocks, despite the formal infinite summation in (5). The global matrix Φ arranges the Φ₍ᵥ₎(k) blocks in a circular pattern (see the displayed block matrix), reflecting the periodic nature of the model.

Because of the circular structure, only the first set of blocks Φ₍1₎(k) needs to be computed directly; the remaining blocks are obtained by substituting the season‑specific AR coefficients φ₍v₎ⱼ for φ₍1₎ⱼ. This observation can reduce the effort of forming Φ, especially when the period S is large but the AR order p is modest.

Solving the linear system Φ γ = ζ yields the unique start‑up autocovariances provided the underlying PARMA model is causal. The method is analytically simple: all matrices are given in closed form, and the circular pattern aligns naturally with the periodicity of the process.

However, the author also points out a significant limitation. Since the system involves the full S(p + 1) × S(p + 1) matrix, a direct solution (e.g., Gaussian elimination) requires O(