options crt;  ? Multivariate ARCH via ML
?  Example by Clint Cummins 9/2/92
?  p=neq=4 equations
?  n=2 lags  (pure ARCH, not GARCH)
?
?  basic multivariate normal likelihood function;
?  define the components below
?
frml march logl = -.5*log(dets) -.5*eie/dets + con*neq;
set con = -.5*log(6.2831853);   ? constant term -- log(2*pi)
?
set neq = 4;
?
?  s.. = sigma matrix
?  a.. = constant part of sigma matrix
?  b.. = coefficients of first lag
?  c.. = coefficients of second lag
?  e12 = equation 1, lag 2 -- define this below
?  e21 = equation 2, lag 1, etc.
?
frml s11 a11 + b11*e11*e11 + c11*e12*e12;
frml s21 a21 + b21*e21*e11 + c21*e22*e12;
frml s31 a31 + b31*e31*e11 + c31*e32*e12;
frml s41 a41 + b41*e41*e11 + c41*e42*e12;
frml s22 a22 + b22*e21*e21 + c22*e22*e22;
frml s32 a32 + b32*e31*e21 + c32*e32*e22;
frml s42 a42 + b42*e41*e21 + c42*e42*e22;
frml s33 a33 + b33*e31*e31 + c33*e32*e32;
frml s43 a43 + b43*e41*e41 + c43*e42*e42;
frml s44 a44 + b44*e41*e41 + c44*e42*e42;
list se s11 s21 s31 s41 s22 s32 s42 s33 s43 s44;
list ie i11 i21 i31 i41 i22 i32 i42 i33 i43 i44;
list ap a11 a21 a31 a41 a22 a32 a42 a33 a43 a44;
list bp b11 b21 b31 b41 b22 b32 b42 b33 b43 b44;
list cp c11 c21 c31 c41 c22 c32 c42 c33 c43 c44;
list allp ap bp cp;
?
?  i.. = sigma-inverse (times dets) = cofactors
?
frml i11  s22*(s33*s44-s43*s43) - s32*(s32*s44-s43*s42) + s42*(s32*s43-s33*s42);
frml i21 -s21*(s33*s44-s43*s43) + s32*(s31*s44-s43*s41) - s42*(s31*s43-s33*s41);
frml i31  s21*(s32*s44-s43*s42) - s22*(s31*s44-s43*s41) + s42*(s31*s42-s32*s41);
frml i41 -s21*(s32*s43-s33*s42) + s22*(s31*s43-s33*s41) - s32*(s31*s42-s32*s41);
frml i22  s11*(s33*s44-s43*s43) - s31*(s31*s44-s43*s41) + s41*(s31*s43-s33*s41);
frml i32 -s11*(s32*s44-s43*s42) + s21*(s31*s44-s43*s41) - s41*(s31*s42-s32*s41);
frml i42  s11*(s32*s43-s33*s42) - s21*(s31*s43-s33*s41) + s31*(s31*s42-s32*s41);
frml i33  s11*(s22*s44-s42*s42) - s21*(s21*s44-s42*s41) + s41*(s21*s42-s22*s41);
frml i43 -s11*(s22*s43-s32*s42) + s21*(s21*s43-s32*s41) - s31*(s21*s42-s22*s41);
frml i44  s11*(s22*s33-s32*s32) - s21*(s21*s33-s32*s31) + s31*(s21*s32-s22*s31);
?
?  dets = determinant of sigma
?
frml dets s11*i11 + s21*i21 + s31*i31 + s41*i41;
?
?  eie = e'(sigma-inverse)*e (times dets)
?
frml eie e1*(e1*i11 + 2*e2*i21 + 2*e3*i31 + 2*e4*i41)
       + e2*(           e2*i22 + 2*e3*i32 + 2*e4*i42)
       + e3*(                      e3*i33 + 2*e4*i43)
       + e4*                                  e4*i44;
?
? simple example of structural equations and their lags
?  e12 = equation 1, lag 2 -- define this below
?
frml e1 y1-a1-b1*x;
frml e11 y1(-1)-a1-b1*x(-1);
frml e12 y1(-2)-a1-b1*x(-2);
frml e2 y2-a2;
frml e21 y2(-1)-a2;
frml e22 y2(-2)-a2;
frml e3 y3-a3;
frml e31 y3(-1)-a3;
frml e32 y3(-2)-a3;
frml e4 y4-a4;
frml e41 y4(-1)-a4;
frml e42 y4(-2)-a4;
list elag0 e1-e4;
list e1lags e11-e12;
list e2lags e21-e22;
list e3lags e31-e32;
list e4lags e41-e42;
param a1-a4 b1;
set T = 100;
smpl 1,T;   ? draw some data for example
random(seed=491479) x;
random y1; y1 = 1+x+y1;
random(mean=2) y2;
random(mean=3) y3;
random(mean=4) y4;
?
smpl 3,T;
lsq e1-e4;   ? estimate with LSQ to get Sigma starting values
?
param allp;  ? sets all covariance parameters to zero
eqsub march dets eie ie se elag0 e1lags-e4lags;
? initialize covariance parameters
set a11 = @covu(1,1);
set a21 = @covu(2,1);
set a31 = @covu(3,1);
set a41 = @covu(4,1);
set a22 = @covu(2,2);
set a32 = @covu(3,2);
set a42 = @covu(4,2);
set a33 = @covu(3,3);
set a43 = @covu(4,3);
set a44 = @covu(4,4);
options limnum=9999,limwarn=0;
const bp cp;  ? hold arch params at zero to start with
ml march;
param bp;  ? allow arch params on first lag to be non-zero
if (@ifconv); then; ml march;
param cp;  ? allow arch params on second lag to be non-zero
           ? (does not converge on this data, even in 100 iterations)
           ? (but note that the data-generating parameters are zero)
if (@ifconv); then; ml march;