options crt;  ? FIML via ML
name fiml4 'FIML for 4 Equation Model Using ML Procedure' ;
?
?    TSP 4.2A (Sept 1992)
?    Author: Clint Cummins
?
?  instead of -.5*log(det(sigma)), make simplification:
?   det(sigma) = 1/det(sigma-inverse) = 1/DL1*DL1*DL2*DL2*...
?   log(det(sigma)) = -[2*log(DL1) + 2*log(DL2) + ...]
?                   = -2*ldetd     (and -.5*(-2) = 1)
?                  or could use log(DL1*DL1) if there are negative problems
?
set con = -.5*log(6.2831853);
frml fiml logl = ldetd + ldetjac -.5*epsie + con*neq;
?
?  parameterize sigma inverse in terms of its LDL' factorization,
?  so that its determinant is simplified
?  s.. = sigma inverse element (lower triangle)
?  d.  = D element (diagonal matrix)
?  DL. = square root of D element
?  L.. = L element (unit triangular:  ones on diagonal, zeros above diagonal)
?  a.. = A element (A = LD)
?
?  2 x 2 case example:
?
?  sigma inverse = s11 s21 = [  1  0 ] * [ DL1*DL1    0    ] * [ 1 L21 ]
?                  s21 s22 = [ L21 1 ] * [    0    DL2*DL2 ] * [ 0  1  ]
?
?  These elements of sigma inverse are the same for higher order models;
?  just add new rows to sigma inverse.
?
frml d1 DL1*DL1;
frml s11 d1;
list es1 s11 d1;  ? list of substitution equations
list fp1 DL1 1;   ? list of covariance parameters
frml det1 ldetd = log(DL1);
frml ese1 epsie = e1*e1*s11;
?
frml s21 s11*L21;
frml d2 DL2*DL2;
frml s22 s21*L21 + d2;
list es2 s22 d2 s21 es1;
list fp2 fp1 L21 DL2 1;
frml det2 ldetd = log(DL1) + log(DL2);
frml ese2 epsie = e1*(e1*s11 + 2*e2*s21) + e2*e2*s22;
?
frml s31 s11*L31;
frml a32 d2*L32;
frml s32 s31*L21 + a32;
frml d3 DL3*DL3;
frml s33 s31*L31 + a32*L32 + d3;
list es3 s33 d3 s32 a32 s31 es2;
list fp3 fp2 L31 L32 DL3 1;
frml det3 ldetd = log(DL1) + log(DL2) + log(DL3);
frml ese3 epsie = e1*(e1*s11 + 2*e2*s21 + 2*e3*s31) +
                  e2*(e2*s22 + 2*e3*s32) + e3*e3*s33;
?
frml s41 s11*L41;
frml a42 d2*L42;
frml s42 s41*L21 + a42;
frml a43 d3*L43;
frml s43 s41*L31 + a42*L32 + a43;
frml d4 DL4*DL4;
frml s44 s41*L41 + a42*L42 + a43*L43 + d4;
list es4 s44 d4 s43 a43 s42 a42 s41 es3;
list fp4 fp3 L41 L42 L43 DL4 1;
frml det4 ldetd = log(DL1) + log(DL2) + log(DL3) + log(DL4);
frml ese4 epsie = e1*(e1*s11 + 2*e2*s21 + 2*e3*s31 + 2*e4*s41) +
                  e2*(e2*s22 + 2*e3*s32 + 2*e4*s42) +
                  e3*(e3*s33 + 2*e4*s43) + e4*e4*s44;
?
frml e1 y1-a1-b1*x;   ? simple example of structural equations
frml e2 y2-a2;
frml e3 y3-a3;
frml e4 y4-a4;
param a1-a4 b1;
smpl 1,50;            ? draw some data for example
random x; random y1; y1 = 1+x+y1;
random(mean=2) y2;
random(mean=3) y3;
random(mean=4) y4;
?
? Estimate with regular FIML for comparison and to obtain 
? Sigma start values.
?
FIML(endog=(y1-y4)) e1-e4;   
?
frml Lj4 ldetjac=0;          ? for identity Jacobian
set neq=4;
param fp4;
eqsub fiml det4 ese4 es4 e1-e4 Lj4;
?
? Initialize covariance parameters
?
inv @covu si;  yldfac si dsig tsig;
set DL1 = sqrt(dsig(1,1));
set DL2 = sqrt(dsig(2,2));
set DL3 = sqrt(dsig(3,3));
set DL4 = sqrt(dsig(4,4));
set L21 = tsig(1,2);
set L31 = tsig(1,3);
set L32 = tsig(2,3);
set L41 = tsig(1,4);
set L42 = tsig(2,4);
set L43 = tsig(3,4);
ml(hiter=n,hcov=nb) fiml;
?
? Display estimated covariance matrix in matrix form.
?
smpl 1 4 ;    ? Temporary for construction of 4 by 4 matrix.
dot 1-4 ; set DS. = DL.*DL. ; enddot ;
mmake dmat DS1-DS4 ;
mat dmat = diag(dmat) ;
mmake col1 1   0   0   0   ; unmake col1 col1 ;
mmake col2 L21 1   0   0   ; unmake col2 col2 ;
mmake col3 L31 L32 1   0   ; unmake col3 col3 ;
mmake col4 L41 L42 L43 1   ; unmake col4 col4 ;
mmake lmat col1-col4 ;
mat sigma = (lmat'dmat*lmat)" ;    ? Sigma = L'DL inverse
print sigma ;
smpl 1,50 ;   ? Restore SMPL if you plan to do more estimation.
stop ; end ;