proc iml; start Loglik_org(parm) global(eps,R_bar,T,N,Q,R,D,h,y,s); a=parm[1];b=parm[2]; T=nrow(eps); N=ncol(eps); R_bar=1/T#eps`*eps; Q=J(N*T,N); R=J(N*T,N); D=J(N*T,N); h=J(T,1); y=J(T,1); s=J(T,1); Q[1:N,]=R_bar; D[1:N,]=sqrt(inv(diag(Q[1:N,]))); R[1:N,]=D[1:N,]*Q[1:N,]*D[1:N,]; G=root(R[1:N,]); det=2*sum(log(vecdiag(G))); square=eps[1,]*eps[1,]`; h[1]=det; y[1]=eps[1,]*inv(R[1:N,])*eps[1,]`; s[1]=square; do i=2 to T; w=N*(i-1)+1; u=N*(i-2)+1; m=N*i; v=N*(i-1); Q[w:m,]=R_bar*(1-a-b)+a*eps[i-1,]`*eps[i-1,]+b*Q[u:v,]; D[w:m,]=sqrt(inv(diag(Q[w:m,]))); R[w:m,]=D[w:m,]*Q[w:m,]*D[w:m,]; G=root(R[w:m,]); det=2*sum(log(vecdiag(G))); square=eps[i,]*eps[i,]`; h=det; y=eps[i,]*inv(R[w:m,])*eps[i,]`; s=square; end; f=-0.5*sum(h+y-s); return(f); finish; store module=Loglik_org; quit; according to this formulation, the matrix R[w:m,] which is one block of matrix R, is no way to be negative definite. I do not understand how this is possible. You see in the definition, the block of R evolves as a weighted average of a positive definite matrix eps[i-1,]`*eps[i-1,] and a semi-definite matrix Q[w:m,] which is initialized by R_bar. R_bar is simply the sample covariance matrix of eps.
... View more
