multivariate GARCH model using SAS IML

[Hanyu]

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.

[Rick_SAS]

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.

[Hanyu]

Thank you Dr Wicklin. Can you tell me how to output the Hessian matrix from the subroutine NLPRNA?

[Rick_SAS]

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;

[Rick_SAS]

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)

[Hanyu]

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.

[Rick_SAS]

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 */

[Hanyu]

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.

[Hanyu]

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.