options crt double;
name mgarch 'Multivariate GARCH with BEKK parameterization';
? Following Engle and Kroner "Multivariate Simultaneous Generalized ARCH",
? Econometric Theory, March 1995
? by Clint Cummins 8/23/96
const nob,100;
smpl 1,nob;
supres smpl;
? generate artificial data from true model
random(seedin=90841);
?  2 equations, GARCH(1,1), and k=1 here for exposition
const neq,2 maxlag,1;

?mat ht = c0'c0 + a1'(epslag*epslag')*a1 + g1'htlag*g1;
?  for k=2, we would have a2 and g2 matrices (with zero restriction in them)

read(type=tri,nrow=neq) c0;
1 1.5
  2;
list c0ps c0_11 c0_12 c0_22;
list c00ps c0_11 0 c0_12 c0_22;  ? elements by column w/ zeros
param c0ps;
?mmake ttmp c0ps;
?mat c0 = tri(ttmp);  ? this method does not work at present
?mmake vtmp c00ps;
?mform(type=tri,nrow=neq,ncol=neq) c0 = vtmp;

read(type=gen,nrow=neq,ncol=neq) a1;
.2   .25
-.15 .3;
list a1ps a1_11 a1_21 a1_12 a1_22;  ? list elements by column for mmake/mform
param a1ps;
?mmake vtmp a1ps;
?mform(type=gen,nrow=neq,ncol=neq) a1 = vtmp;

read(type=gen,nrow=neq,ncol=neq) g1;
.3  -.25
.15 .2;
list g1ps g1_11 g1_21 g1_12 g1_22;  ? list elements by column for mmake/mform
param g1ps;
?mmake vtmp g1ps;
?mform(type=gen,nrow=neq,ncol=neq) g1 = vtmp;

?mat ht = c0'c0 + a1'(epslag*epslag')*a1 + g1'htlag*g1;
? check covariance stationarity
mat agk = (a1#a1)' + (g1#g1)';
mat eigr1 = eigval(agk);
mat maxmod = max( sqrt(eigr1%eigr1 + @eigvali%@eigvali) );
print maxmod;  ? must be less than 1
set neq2 = neq*neq;
mat Ineq2 = ident(neq2);
mat v2tmp = (Ineq2 - agk)"vec(c0'c0);          ? 2.6b
mform(type=sym,nrow=neq,ncol=neq) h0 = v2tmp;  ? unconditional variance
  print v2tmp h0;
mat ht = h0;

list(last=neq,prefix=e) es;  ? e1-e2
list(last=neq,prefix=eps) epss;  ? eps1-eps2
list(last=neq) nums;  ? 1-2
random es x1 x2;
frml xb1 b1_0 + b1_1*x1;   ? covariates for y1,y2
frml xb2 b2_0 + b2_2*x2;
param b1_0,1 b1_1,1 b2_0,1.5 b2_2,2;
list betas b1_0 b1_1 b2_0 b2_2;
genr xb1 yhat1;  ? all that remains is to add the residual eps1
genr xb2 yhat2;

mat c0pc0 = c0'c0;
trend t;
do it=1,nob;
  select it=t;
  mmake e es;   ? e1,e2 for current observation
  mat eps = (e*chol(ht))' ;
  unmake eps epss;
  dot nums;
    y. = yhat. + eps.;
  enddot;
  mat ht = sym(   ? the sym() function should not be needed here....
    c0pc0 + (a1'eps)*(a1'eps)' + g1'ht*g1 );  ? update ht using BEKK eqn
enddo;
select 1;

? Estimate the parameters for the above artificial data

? Use true values as starting values.  This involves copying the
? matrix values to the param names.  The betas are already params.
mat ttmp = vech(c0);
unmake ttmp c0ps;
mat vtmp = vec(a1);
unmake vtmp a1ps;
mat vtmp = vec(g1);
unmake vtmp g1ps;

? Last-minute initializations for PROC
set pi = 4*atan(1);  ? tan(pi/4) = 1
set loglcons = neq*log(2*pi);
logl = 0;

? performance with default HITER=BFGS is not so good
?ML mgarch betas c0ps a1ps g1ps;

? try with HITER=B  (feature implemented 11/96)
ML(hiter=b,hcov=b) mgarch betas c0ps a1ps g1ps;

? Proc mgarch evaluates the log likelihood for each obs.  It's pretty
? similar to the data generation code above.
?===============================================

Proc mgarch;
? Check parameter constraints
? 1. (1,1) elements of BEKK matrices positive, just to insure uniqueness.
if (c0_11 <= 0) .or.
   (a1_11 <= 0) .or.
   (g1_11 <= 0) ; then; goto 90;

? construct BEKK matrices from params
mmake vtmp c00ps;
mform(type=tri,nrow=neq,ncol=neq) c0 = vtmp;
mmake vtmp a1ps;
mform(type=gen,nrow=neq,ncol=neq) a1 = vtmp;
mmake vtmp g1ps;
mform(type=gen,nrow=neq,ncol=neq) g1 = vtmp;

?2. check covariance stationarity
mat agk = (a1#a1)' + (g1#g1)';
mat eigr1 = eigval(agk);
mat maxmod = max( sqrt(eigr1%eigr1 + @eigvali%@eigvali) );
if (maxmod >= 1); then; goto 90;

mat c0pc0 = c0'c0;
mat v2tmp = (Ineq2 - agk)"vec(c0pc0);          ? 2.6b
mform(type=sym,nrow=neq,ncol=neq) ht = v2tmp;  ? unconditional variance

dot nums;
  genr xb. yhat.;
  eps. = y. - yhat.;
enddot;
? Make these initializations above the ML command, so they don't have
? to be repeatedly executed each time the PROC is called
?set pi = 4*atan(1);  ? tan(pi/4) = 1
?set loglcons = neq*log(2*pi);
do it=1,nob;
  select it=t;
  mmake eps epss;
  mat eps = eps';
  mat ehe = eps'ht"eps;
  mat @det = det(ht);
  ?if @det <= 0; then; goto 90;  ? we don't have to check this because BEKK
  ? guarantees positive definiteness
  logl = -(loglcons + log(@det) + ehe)/2;
  mat ht = c0pc0 + (a1'eps)*(a1'eps)' + g1'ht*g1;  ? update ht using BEKK eqn
enddo;
select 1;
mat @logl = sum(logl);
goto 100;

? error return
90 set @logl = @miss;

100 nodbug;
endproc;