I am trying to estimate a multivariate garch model and I want to use sas iml becuase the proc varmax only provide one simple BEKK model. . I have read the article Maximum likelihood estimation in SAS/IML - The DO Loop But the likelihood function of a multivariate garch model contain lag of matrix. Could any one give me a hint? Thank you in advance. The specification of a multivariate GARCH could be found in sas ets manual.
The SAS/IML language supports the LAG and DIF function, which you might find useful.
Here's an article that demonstrates the LAG function: http://blogs.sas.com/content/iml/2012/04/30/the-lag-function/
You might also want to read this article about vectorizing time series computations: http://blogs.sas.com/content/iml/2014/01/13/how-to-vectorize-time-series-computations/
If you want help with your code, supply statements that compute the likelihood function in PROC IML. Even if you aren't sure of every step, having a template to work from will be useful for those who want to help.
Thank you Dr Wicklin. Can you tell me how to output the Hessian matrix from the subroutine NLPRNA?
Do you mean "evaluate the Hessian at the optimal solution and print it out"? If so, then just pass the optimal solution to a function that evaluates the Hessian. If the Hessian is too complex to program (or there is not analytical expression for it), then you can use the NLPFDD subroutine to compute finite difference derivatives (FDDs). For example, in the MLE blog post that you reference, I could evaluate the objective function, gradient, and Hessian matrix at the optimal solution by using the followiong statement:
call nlpfdd(f,gradient,hessian,"LogLik",result);
print f,gradient, hessian;
Can you provide a driver program that calls LogLik and that defines eps as a small matrix, maybe 4 rows and 3 columns? (T and N small)
2764
2765 proc iml;
NOTE: IML Ready
2766 use intra_garch_norm;
2767 read all var{at be fr de it nl es} into eps [colname=country];
2768 close intra_garch_norm;
2769
2770 con={0 0,
2771 1 1};
2772 opt={1 4};
2773 p={0.4 0.4};
2774 tc=1000;
2775 load module=Loglik_org;
NOTE: Opening storage library WORK.IMLSTOR
2776 call nlpnra(rc,result,"Loglik_org",p,opt,con);
ERROR: (execution) Matrix should be positive definite.
operation : ROOT at offset 30 column 7
operands : _TEM1002
_TEM1002 7 rows 7 cols (numeric)
1 0.0429409 0.0924234 -0.516006 0.3325295 -0.254749 0.2209257
0.0429409 1 0.0785693 -0.595686 0.3531359 -0.271494 0.3314657
0.0924234 0.0785693 1 -0.71971 0.5530654 -0.340426 0.4299508
-0.516006 -0.595686 -0.71971 1 -0.670053 0.4233309 -0.556491
0.3325295 0.3531359 0.5530654 -0.670053 1 -0.24287 0.6664104
-0.254749 -0.271494 -0.340426 0.4233309 -0.24287 1 -0.126188
0.2209257 0.3314657 0.4299508 -0.556491 0.6664104 -0.126188 1
statement : ASSIGN at offset 30 column 1
traceback : module LOGLIK_ORG at offset 30 column 1
NOTE: PROCEDURE IML used (Total process time):
real time 6:34.24
cpu time 6:30.61
Hi Dr. Ricklin, I got this error from sas. Does this mean that the second parameter is the reason that the matrix is not positive definite.
The error message says that in your module called LOGLIK_ORG you are calling the ROOT function. At some point in the optimization, the matrix that you are attempting to decompose is not positive definite. Therefore the ROOT function is gving an error.
If you copy/paste that matrix into a matrix called X, then you can see that the matrix is not positive definite:
v = eigval(X);
print v; /* one eigenvalue is negative */
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.
I try to locate this negative definite matrix in the matrix R using loc function. I cannot find any matrix has these values. I cannot understand how is this possible.
SAS Innovate 2025 is scheduled for May 6-9 in Orlando, FL. Sign up to be first to learn about the agenda and registration!
Learn how to run multiple linear regression models with and without interactions, presented by SAS user Alex Chaplin.
Find more tutorials on the SAS Users YouTube channel.