/*
**  maxsimeq.e  -  Simultaneous Equation Models
**
**  Supply three procs, one that returns _seq_Gamma, the matrix
**  of coefficients of the regression of the endogenous
**  variables on themselves given the vector of parameters,
**  _seq_Beta, the matrix of the coefficients of the regression
**  of the endogenous variables on the exogenous variables,
**  and _seq_B0, the vector of constants.  Thus
**
**       Y * _seq_Gamma = X * _seq_Beta + _seq_B0 + U
**
**  where Y is the NxL matrix of endogenous variables, X is
**  the NxK matrix of exogenous variables, and U is a residual
**  matrix.
**
**  Also required are two index vectors, _seq_Endo, and _seq_Exo,
**  that specify which columns of the dataset are the endogenous
**  variables in Y and exogenous variables in X respectively.
**
**
**  For this example, the model is
**
**     Y1 = B1 Y2 + g1 X1 + c1 + U1
**     Y2 = B2 Y1 + g2 X2 + c2 + U2
**
**  where Y1 and Y2 are endogenous, X1 and X2 are exogenous,
**  and Z1 and Z2 are the residuals.
*/


library maxlik;
#include maxlik.ext;
maxset;


_seq_Endo = { 1, 2 };   /* endogenous variables */
_seq_Exo  = { 3, 4 };   /* exogenous variables  */

proc _seq_Gamma(b);
    local g;
    g = eye(2);
    g[1,2] = -b[1];
    g[2,1] = -b[2];
    retp(g);
endp;

proc _seq_Beta(b);
    local be;
    be = zeros(2,2);
    be[1,1] = b[3];
    be[2,2] = b[4];
    retp(be);
endp;

proc _seq_B0(b);
     local c;
     c = zeros(2,1);
     c[1] = b[5];
     c[2] = b[6];
     retp(c);
endp;


/*
** start values
*/

start = .5 * ones(6,1);


/*
** dataset
*/

dataset = "maxsimeq";


/*
**  output file
*/

output file = maxsimeq.out reset;


_ln2pi = -0.5 * rows(_seq_Endo) * ln(2*pi);  /* constant needed for the */
                                             /* log-likelihood function */

{ b1,f1,g1,cov1,ret1 } = maxlik(dataset,0,&lpr,start);

__title = "Wald-type confidence limits";
cl2 = maxtlimits(b1,cov1);
call maxclprt(b1,f1,g1,cl2,ret1);


output off;



/*
**  Likelihood function
*/


proc lpr(b,z);
    local u, psi, oldt, ldet0, ipsi, beta, gamm, cnst;
    gamm = _seq_Gamma(b);
    beta = _seq_Beta(b);
    cnst = _seq_B0(b);

    u = z[.,_seq_Endo] * gamm - z[.,_seq_Exo] * beta - cnst';

    psi = u'u/rows(u);
    oldt = trapchk(1);
    trap 1,1;
    ipsi = invpd(psi);
    trap oldt,1;
    if scalmiss(ipsi);
        retp(error(0));
    endif;

    ldet0 = ln(detl);
    u = u * chol(ipsi)';

    retp(_ln2pi + ln(abs(det(gamm))) - 0.5*(ldet0 + sumc((u.*u)')));
endp;