## Exact PDF of an ARMA process - Very efficient method using filtering

In spite of the speedup of the Levinson algorithm, notice that for , 100 milliseconds are needed, and at the rate of , 10 seconds would be needed for . In some applications, that may be excessive. An even more efficient method exists that based on the idea of a whitening filter. The basic idea is the following. Let be any lower-triangular (LT) matrix such that

Thus,

Now consider the vector . The covariance of is the identity matrix:

So, is said to be whitened, because it consists of iid samples of Gaussian noise with mean 0 and variance 1, denoted by N(0,1), having a flat power spectrum. The log PDF of can be written in terms of u.
   lpx = -N/2*log(2*pi) + ldetC -.5*sum(u.^2);

where ldetC is the log-determinant of . We may, if we want, break down into an LT part with diagonal elements equal to 1, and a diagonal matrix ,

Then, we have

Thus,

Because we imposed the constraint that is lower-triangular with diagonal elements of 1, has a determinant of 1. The determinant of R is therefore the same as the product of the diagonal elements of . The whitened vector u is not unique because any matrix which accomplishes the whitening, in fact any linear transformation which accomplishes this is suitable. Note that we can obtain by Cholesky factorization. By convention, the Cholesky factor is upper-triangular, therefore, in MATLAB, we need to twist it around:
   C = fliplr(flipud(chol(inv(R))));


For an ARMA process, it is only necessary to reverse the AR and MA coefficients to obtain an inverse filter that recovers the innovation sequence and whitens the data:

   u = filter(a,b,x)/sqrt(sig2);

If we ignore filter startup errors, the data PDF could be written:
>> u=filter(a,b,x) / sqrt(sig2);
>> lpx = -N/2*log(2*pi) +ldetC -.5*sum(u.^2);

Note that by comparing the filtering approach with the above matrix approach, we can deduce that the rows of the matrix converge to the impulse response of the whitening filter and that the diagonal elements of converge to . Thus, ldetC = -N/2*log(sig2).

But, not knowing the correct filter initial values causes a startup transient that dies away at a rate depending on the numerator polynomial . After some amount of time, say samples, the samples of u will agree closely with e. There is a requirement that the MA polynomial

has roots inside the unit circle. We can determine how long we have to wait by solving for the autocorrelation function of the inverse of the ARMA process.
>> ri = real(ifft( 1./ rho));
>> K = max(find(abs(rbi(1:floor(N/2))) > ri(1)/30));

Note that due to the symmetry of the ACF, the time-reversed whitener will also work, but will produce a sequence that does not converge to .
   v = filter(a,b,flipud(x))/sqrt(sig2);

Our approach will be to use the exact approach for the first K values and the filter method for indexes [K+1:N]. As increases, the overhead of computing the first K values will be constant and we will have an algorithm with order execution time. We end up with a hybrid filtering/exact approach:
>> %% PART 1 - EXACT on first K samples
>> [y,ldetR] = toepsolve(r(1:K),x(1:K),K);
>> lpx1 = -K/2*log(2*pi) -.5*ldetR -.5*x(1:K)'*y;
>> %% PART 2 - FILTERING APPROACH
>> u=filter(a,b,x) / sqrt(sig2); % filtering approach
>> lpx2 = -(N-K)/2*log(2*pi)  -(N-K)/2*log(sig2) -.5*sum(u(K+1:end).^2);
>> lpx = lpx1 + lpx2; % combine

In PART 2, we have made use of the equivalence of the filtering approach and the matrix approach to deduce the determinant term
>>   -(N-K)/2*log(sig2)


Note that above, the first elements of u are not white". If we really need the first samples of u to be iid, we can correct" them:

>> Rb=toeplitz(r(1:K));
>> C = fliplr(flipud(chol(inv(Rb))));
>> u(1:K) = C * x(1:K); % exact "fix" of first K samples

however, it would be only for academic purposes - in general we don't need u(1:K). In Figure 9.4, we compare u generated from u=filter(a,b,x)/sqrt(sig2); with e and with the corrected" u. See software/pdf_arma_exact.m An example is run by script software/test_arma.m.
For more details see software/pdf_arma_exact.m, which is used in software/module_acfx.m and software/test_arma.m. There is more on this in Sections 9.4.10 and 9.4.10.

Baggenstoss 2017-05-19