## AR Data PDF - exact

The exact PDF of an AR process is given by (9.6) with the ACF equal to the inverse FT of (9.40). The ACF may also be obtained (inefficiently) from the reverse Levinson algorithm
  r = rlevinson([a(:); zeros(N-P-1,1)],sig2);

In practice, we would obtain the length- ACF by the inverse DFT of the circular PSD (9.41) using a large enough transform of size that the ACF dies to zero for .

In the following analysis, we assume that the AR coefficients are known. The problem of computing (9.6) knowing just is inherently different because we must estimate the AR coefficients (See sections 9.4.3). Efficient means of computing the exact PDF for the more general ARMA model are provided in Sections 9.2.2 and 9.2.3. The special properties of the AR model can be used to advantage. First of all, , the impulse response length of the whitening filter is exactly . Thus, only the first samples of the whitened sequence need be corrected".

Also, a closed form expression for the inverse of the ACF matrix can be found. Suppose we obtain the AR coeffients of all orders up to and including . Let r(i) be the i-th ACF coefficient (lag i-1). Let Aa(p,:) be the -th order AR coefficient vector and e(p) the corresponding prediction error variance (note that ) :

  Aa = zeros(P+1,P+1);
for p=1:P+1,
[Aa(p,1:p),e(p)]=levinson(r(1:p));
end;

Now consider the operation :
  u(1) = (x(1)*Aa(1,1))/sqrt(e(1));
u(2) = (x(2)*Aa(2,1) + x(1)*Aa(2,2))/sqrt(e(2));
:            :
u(P) = (x(P)*Aa(P,1) + ...  + x(1)*Aa(P,P))/sqrt(e(P));
u(P+1) = (x(P+1)*Aa(P+1,1) + ...  + x(1)*Aa(P+1,P+1))/sqrt(e(P+1));
u(P+2) = (x(P+2)*Aa(P+1,1) + ...  + x(2)*Aa(P+1,P+1))/sqrt(e(P+1));
u(P+3) = (x(P+3)*Aa(P+1,1) + ...  + x(3)*Aa(P+1,P+1))/sqrt(e(P+1));
:            :

Notice that above index it is the same as u = filter(a,1,x);. Under the assumption that is a -th order AR process, this produces a sequence that is iid N(0,1) for i>P. The first equations can be written in matrix form as

 (9.42)

where and are the first elements of and , and where is the -by- diagonal matrix with with diagonal elements equal to the innovation variance from the levinson algorithm up to order : [e(1) e(2) ... e(P)]

 (9.43)

(shown for ). Note that this provides a way to write in terms of the AR coefficients since

It also gives us the determinant of since since . Also det(S) = prod(e(1:P)) . We can also deduce the determinant of the matrix because for i > P, e(i).

So, to obtain an efficient means of computing the exact PDF of for an AR process, which will have execution time of order , we can use:

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

Then, repace u(1:P) with the result of (9.42) on just the first samples. The log-determinant of is
   ldetR = sum(log(e(1:P))) + (N-P)*log(e(P+1));

For details, see software/pdf_ar_exact.m, or the more general function software/pdf_arma_exact.m.

Problem 2   In MATLAB, create a -by- Toeplitz ACF matrix and find its Cholesky factorization:
  Rp=toeplitz(rb(1:imax))
C = chol(Rp);

What is the relationship between C and A? How can the filtering operation above be written in terms of C?

Baggenstoss 2017-05-19