/*
**   Confirmatory Factor Analysis
**
**   The source code below estimates the parameters of a confirmatory
**   factor analysis model.  Let
**
**      Y = Lambda * eta + E
**
**   where Y is an NxK matrix of observations on K variables, eta
**   is a NxL matrix of unobserved latent variables, Lambda a KxL
**   coefficient (or "loadings") matrix, and E an NxK matrix of
**   unobserved residuals (or "measurement error").  Then
**
**     Sigma = Lambda * Phi * Lambda' + Theta
**
**   where Sigma is the moment matrix of the variables in Y, Phi is
**   moment matrix of the latent variables, and Theta is the moment
**   matrix of the residuals or measurement errors.  If we assume
**   the variables are measured as deviations from their means (an
**   assumption without consequence to the estimation), the moment
**   matrices are all covariance matrices.
**
**   The elements of Lambda, Phi, and Theta may be estimated by
**   assuming multinormality of the variables in Y, eta, and E.
**   When this rather strong assumption seems untenable, the
**   PML or pseudo-maximum likelihood standard errors are recommended
**   (see Arminger and Schoenberg, 1989, "Pseudo maximum likelihood
**   estimation and a test for misspecification in mean- and covariance-
**   structure models", _Psychometrika_, 54: 409-425).
**
**
**   For particular applications, you must supply three procs, putPars,
**   putPhi, and putTheta which construct Lambda, Phi, and Theta given,
**   the vector of parameters, and a vector, numPars, with three elements
**   indicating the number of parameters to be estimated in each matrix.
**   MAXLIK does the rest.
**
*/


library maxlik;
#include maxlik.ext;
maxset;

output file = maxfact.out reset;


/*
**   PutLambda, PutPhi, PutTheta  -  enters the parameters in the
**   parameter vector into appropriate locations in the model matrices
*/

proc PutLambda(b,n);
    local lambda;
     lambda = zeros(4,2);

     if n;  /* Inserts a 1 for normalization for "reference indicators" */
            /* in function calculation.  Not inserted for derivatives.  */
        lambda[1,1] = 1;
        lambda[3,2] = 1;
     endif;

     lambda[2,1] = b[1];
     lambda[4,2] = b[2];

     retp(lambda);
endp;

proc PutPhi(b);
     local phi;
     phi = xpnd(b[3:5]);
     retp(phi);
endp;

proc PutTheta(b);
     local theta;
     theta = eye(4).*b[6:9];
     retp(theta);
endp;


/*
**  numPars - vector indicating the number of parameters to be
**            estimated in each model matrix
*/

numPars = { 2, 3, 4 };


/*
**  start values and labels for parameters
*/

start = { .6, .6, 1, .7, 1, .6, .6, .6, .6 };

_max_ParNames = { LA_21, LA_42, PH_11, PH_12, PH_22, TH_11, TH_22,
      TH_33, TH_44 };


__title = "Confirmatory Factor Analysis";

__output = 10;
_max_GradProc = &grd;
_max_HessProc = &hss;

{ b,f,g,cov,ret } = maxlik("maxfact",0,&lpr,start);

call maxprt(b,f,g,cov,ret);

output off;


/*
**   Procedures
*/

proc lpr(b,x);
    local lambda,phi,theta,oldt,logl;

    lambda = putLambda(b,1);
    phi = putPhi(b);
    theta = putTheta(b);

    oldt = trapchk(0);
    trap 1;
    logl = lnpdfn2(x,lambda*phi*lambda'+theta);
    call trapchk(oldt);
    retp(logl);
endp;


proc grd(b,x);
    local lambda,phi,theta,sigma,oldt,si;
    local dsigma,dlambda, n, i, k, d, x1, g;

    n = rows(b);
    g = zeros(rows(x),n);
    lambda = putLambda(b,1);
    phi = putPhi(b);
    theta = putTheta(b);

    oldt = trapchk(1);
    trap 1,1;
    si = invpd(lambda*phi*lambda' + theta);
    trap oldt,1;
    if scalmiss(si);
        retp(error(0));
    endif;
    k = 1;
    i = 1;
    do until i > numPars[1];
        dlambda = putLambda(mu(k,n),0);
        dsigma = si*(dlambda*phi*lambda' + lambda*phi*dlambda');
        g[.,k] = sumc( ((x*dsigma*si).*x)' ) 
             - sumc(diag(dsigma));
        i = i + 1;
        k = k + 1;
    endo;

    i = 1;
    do until i > numPars[2];
        dsigma = si*(lambda*putPhi(mu(k,n))*lambda');
        g[.,k] = sumc( ((x*dsigma*si).*x)' ) 
             - sumc(diag(dsigma));
        i = i + 1;
        k = k + 1;
    endo;

    i = 1;
    do until i > numPars[3];
        dsigma = si*putTheta(mu(k,n));
        g[.,k] = sumc( ((x*dsigma*si).*x)' ) 
             - sumc(diag(dsigma));
        i = i + 1;
        k = k + 1;
    endo;
    retp(0.5*g);
endp;
  
  
  
proc hss(b,x);   /* HESSIAN */
    local lambda,phi,theta,sigma,dsigma,oldt,si,h,i,j,n,cnum;
    local dsigmai, dsigmaj, W;

    n = rows(b);
    h = zeros(n,n);
    lambda = putLambda(b,1);
    phi = putPhi(b);
    theta = putTheta(b);

    oldt = trapchk(1);
    trap 1,1;
    si = invpd(lambda*phi*lambda' + theta);
    trap oldt,1;
    if scalmiss(si);
        retp(error(0));
    endif;

    cnum = cumsumc(numPars);

/* lambda, lambda */
    i = 1;
    do until i > cnum[1];
        W = putLambda(mu(i,n),0)*phi*lambda';
        dsigmai = si*(W + W');           
        j = 1;
        do until j > i;
             W = putLambda(mu(j,n),0)*phi*lambda';
            dsigmaj = si*(W + W');           
            h[i,j] = sumc(diag(dsigmai*dsigmaj));
            if i /= j;
               h[j,i] = h[i,j];
            endif;
            j = j + 1;
        endo;
        i = i + 1;
    endo;

/*  phi, lambda */
    i = cnum[1] + 1;
    do until i > cnum[2];
        dsigmai = si*(lambda*putPhi(mu(i,n))*lambda');
        j = 1;
        do until j > cnum[1];
            W = putLambda(mu(j,n),0)*phi*lambda';
            dsigmaj = si*(W + W');           
            h[i,j] = sumc(diag(dsigmai*dsigmaj));
            if i /= j;
               h[j,i] = h[i,j];
            endif;
            j = j + 1;
        endo;
        i = i + 1;
    endo;

/* phi, phi */
    i = cnum[1] + 1;
    do until i > cnum[2];
        dsigmai = si*lambda*putPhi(mu(i,n))*lambda';
        j = cnum[1] + 1;
        do until j > i;
            dsigmaj = si*lambda*putPhi(mu(j,n))*lambda';
            h[i,j] = sumc(diag(dsigmai*dsigmaj));
            if i /= j;
               h[j,i] = h[i,j];
            endif;
            j = j + 1;
       endo;
        i = i + 1;
    endo;

/* psi, lambda */
    i = cnum[2] + 1;
    do until i > cnum[3];
        dsigmai = si*putTheta(mu(i,n));
        j = 1;
        do until j > cnum[1];
            W = putLambda(mu(j,n),0)*phi*lambda';
            dsigmaj = si*(W + W');           
            h[i,j] = sumc(diag(dsigmai*dsigmaj));
            if i /= j;
               h[j,i] = h[i,j];
            endif;
            j = j + 1;
        endo;
        i = i + 1;
    endo;

/* psi, phi */
    i = cnum[2] + 1;
    do until i > cnum[3];
        dsigmai = si*putTheta(mu(i,n));
        j = cnum[1] + 1;
        do until j > cnum[2];
            dsigmaj = si*lambda*putPhi(mu(j,n))*lambda';
            h[i,j] = sumc(diag(dsigmai*dsigmaj));
            if i /= j;
               h[j,i] = h[i,j];
            endif;
            j = j + 1;
        endo;
        i = i + 1;
    endo;

/* psi, psi */
    i = cnum[2] + 1;
    do until i > cnum[3];
        dsigmai = si*putTheta(mu(i,n));
        j = cnum[2] + 1;
        do until j > i;
            dsigmaj = si*putTheta(mu(j,n));
            h[i,j] = sumc(diag(dsigmai*dsigmaj));
            if i /= j;
               h[j,i] = h[i,j];
            endif;
            j = j + 1;
        endo;
        i = i + 1;
    endo;

    retp(-0.5*rows(x)*h);
endp;
  
  

proc mu(i,k);
   retp(seqa(1,1,k).==i);
endp;