CR Bound analysis of ARMA PDF - exact forms.

R=toeplitz(r(1:N)); lpx = -N/2*log(2*pi)-.5*log(det(R)) - .5*x'*inv(R)*x);We note that the derivative of

Ri=inv(R); lpx_theta = -.5*trace( Ri * D ) + .5 * x'*Ri*D*Ri*x;where

D=toeplitz(d(1:N)); lpx_theta = -.5*trace( Ri * D ) + .5 * x'*Ri*D*Ri*x;

We can write `r` in terms of the ARMA parameters in a simple form if we make use of the
fact that for a stable ARMA process, `r` must eventually decay to zero.
For large enough ,
for . Then, we can treat `r` as a circular ACF,
the inverse DFT of a circular power spectrum of size .
Thus, the values of the circular PSD of length
computed using (9.9) will be the same as the true PSD at discrete frequencies.
Therefore, the circular ACF will match the first values of the true ACF.
Let `rho` be the length- circular PSD. If,

r = real(ifft(rho));then the first values of

d = real(ifft(rho_theta));where

In summary,

d = real(ifft(rho_theta)); D=toeplitz(d(1:N)); Ri=inv(R); lpx_theta = -.5*trace( Ri * D ) + .5 * x'*Ri*D*Ri*x;The problem with this approach is the necessity of computing and storing the -by- matrices

y = toepsolve(r,x); % y = inv(toeplitz(r)) * x; d = real(ifft(rho_theta)); D=toeplitz(d(1:N)); H = toepsolve(r(1:N),D); lpx_theta = -.5*trace(H) + .5 * y'*D*y;where

function L=toep_special(B,r,y,N); % solves the problem: % % for i=1:size(B,2), % D=toeplitz(B(1:N,i)); % H = inv(R) * D; % L(i) = -.5*trace(H) + .5 * y'*D*y; % end; % where % R = toeplitz(r(1:N)); % y = inv(R) * x; % % for each column of B without having to store R, H, or D, % or take the inverse of a matrix!!! % If each column of B is the derivative of the power spectrum % with respect to a different parameter, and r is the ACF, % and x is the raw data, then the output are the derivatives of the log-PDF % lpx = -N/2*log(2*pi)-.5*log(det(R)) - .5*x'*inv(R)*x); % with respect to the parametersIn summary, we can compute the exact derivatives of (9.6) with respect to the ARMA parameters as follows:

A=fft([a(:); zeros(2*L-P-1,1)]); B=fft([b(:); zeros(2*L-Q-1,1)]); A2=msq(A); B2=msq(B); h = B2./A2; rho = sig2 * h; r = real(ifft(rho)); E=exp(-1i*2*pi/N/2*[0:2*N-1]' * [0:max(P,Q)]); [y,ldetR] = toepsolve(r(1:N),x(1:N),N); % y = R \ x; Rho_a = -2 * sig2 * real( repmat(h./A,1,P) .* E(:,2:P+1) ); Rho_b = 2 * sig2 * real( repmat(h./B,1,Q) .* E(:,2:Q+1) ); R_a = real(ifft(Rho_a)); R_b = real(ifft(Rho_b)); da = toep_special(R_a(1:N,:),r(1:N),y(1:N),N); db = toep_special(R_b(1:N,:),re(1:N),y(1:N),N);Note that

rho_sig = h; r_sig = real(ifft(rho_sig)); dsig = toep_special(r_sig(1:N),r(1:N),y(1:N),N);however, a more compact form can be shown to be

dsig = -.5*N/sig2 + .5 * y(1:N)'*x(1:N)/sig2;

The second derivatives may also be obtained and the Fisher's information computed. Let , be two general spectral parameters. Then

H_theta = toepsolve(r,D_theta);We can simplify matters if we notice that for two matrices

function Q=toep_fim(B,r,L); % solves the problem: % % R = toeplitz(r); % for i=1:size(B,2), % Di=toeplitz(B(:,i)); % for j=1:size(B,2), % Dj=toeplitz(B(:,j)); % Q(i,j) = .5*trace(inv(R)*Di*inv(R)*Dj); % end; % end; %See

Baggenstoss 2017-05-19