Exact PDF of an ARMA process - Efficient method using Levinson algorithm

The exact PDF of an ARMA process can be written as in (10.6) where ${\bf R}_N$ is the symmetric Toeplitz $N\times N$ covariance matrix of an ARMA process formed from (10.15). Computing ${\bf R}_N$, and even just storing it, for large $N$ is inefficient. Using the Levinson algorithm, an efficient implementation of the exact PDF may be obtained.

When $N$ is large, the computing the term ${\bf x}^\prime {\bf R}_N^{-1} {\bf x}$ is prohibitive and when $N$ is very large, even storing ${\bf R}_N$ is inefficient. The problem is solved efficiently using the Levinson algorithm. Let x be an N-by-1 vector and let ${\bf R}_N$ be the symmetric N-by-N Toeplitz matrix formed from the N-by-1 ACF vector r. The function

   [y,ldetR] = toepsolve(r,x,N);

is the same as

   R = toeplitz(r(1:N));
   y = R \ x;
   ldetR = log(det(R));

We then have

$${\bf x}^\prime {\bf R}_N^{-1} {\bf x}= {\bf y}^\prime {\bf x}.$$

The execution time of software/toepsolve.m is much faster (see Figure 10.3) and includes as a by-product the determinant of R. As predicted by theory, the standard approach is order $N^3$, while Levinson is $N^2$.

Figure 10.3: Comparison of execution times of y = R b versus y = toepsolve(r,b,N) as a function of N.

The exact PDF can be therefore computed by

   [y,ldetR] = toepsolve(r,x,N);
   lpx = -N/2*log(2*pi) -.5*ldetR -.5*x'*y;