library maxlik;
#include maxlik.ext;
maxset;

/*
**  maxnleq.e  -  Nonlinear Simultaneous Equation Model
**
**  The source code below estimates the coefficients of a
**  general nonlinear system of simultaneous equations:
**
**      U = f(Y,X,B)
**
**  where U is an NxK matrix of residuals and f() is an arbitrary
**  function of Y, an NxK matrix of endogenous variables, X, an 
**  NxL matrix of exogenous variables, and B, and Mx1 vector of
**  coefficients.  Thus there are N observations, and f() represents
**  K simultaneous equations.
**
**  The Jacobian of the function with respect to the endogenous 
**  variables is computed numerically unless a proc to do that is
**  provided.  When this proc, _nseq_EndoJcb,is provided, time of
**  computation is reduced significantly.
**
**  If the Jacobian of U w.r.t. Y is an identity matrix, the system
**  is linear and use of the code below would work but would be 
**  significantly slower than a more appropriate method.  See
**  maxsimeq.e for an example of the estimation of a linear
**  simultaneous equation system.
**
**  Required:  a proc to compute the residuals given the coefficients
**             and the data.
**
**  Optional:  a proc to compute the Jacobian of the residuals w.r.t.
**             the endogenous variables, i.e., U w.r.t. Y.
**
**
**  In the model for this example, the endogenous variable Y2
**  interacts with the exogenous variable X1 in equation 1:
**
**     Y1 = b1 Y2 + c1 Y2 * X1 + g1 X1 + d1 + U1
**     Y2 = b2 Y1 + g2 X2 + d2 + U2
**
**  where Y1 and Y2 are endogenous, X1 and X2 are exogenous,
**  and U1 and U2 are residuals.
*/

library maxlik;
#include maxlik.ext;
maxset;



/*
**  Procedure for computing the function
**
**  In the system being estimated here, the endogenous variable Y2
**  interacts with the exogenous variable X1 in equation 1:
**
**     Y1 = b1 Y2 + b2 Y2 * X1 + b3 X1 + b4 + U1
**     Y2 = b5 Y1 + b6 X2 + b7 + U2
**
**  where Y1 and Y2 are endogenous, X1 and X2 are exogenous,
**  and U1 and U2 are residuals.
*/
proc _nseq_Fct(b,Z);
    local u;
    u = zeros(rows(z),2);
    u[.,1] = z[.,1] - b[1]*z[.,2] - b[2]*z[.,2].*z[.,3] - b[3]*z[.,3] - b[4];
    u[.,2] = z[.,2] - b[5]*z[.,1] - b[6]*z[.,4] - b[7];
    retp(u);
endp;



/*
**  Set _nseq_EndoInd to the indices of the endogenous variables
**  in the dataset.
*/

_nseq_EndoInd = { 1, 2 };



/*
** start values
*/

start = { .3, .1, .6, .8, .6, .5, .5 };
  

/*
**  dataset
*/

dataset = "maxnleq";


/*
**  output file
*/

output file = maxnleq.out reset;



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

_max_Options = { BFGS };  /* computing the Hessian is very time consuming */
                          /* so set to BFGS method for greater speed      */

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

__title = "Wald confidence limits";

cl1 = maxtlimits(b1,cov1);
call maxclprt(b1,f1,g1,cl1,ret1);


output off;


/********************
     Procedures
********************/

/*
**  Likelihood function
*/

proc lpr(b,x);
    local dev,psi,ibeta,g0,oldt,ldet0,ipsi,jb;

    dev = _nseq_Fct(b,x);

    psi = dev'dev/rows(x);
    oldt = trapchk(1);
    trap 1,1;
    ipsi = invpd(psi);
    trap oldt,1;
    if scalmiss(ipsi);
        retp(error(0));
    endif;
    ldet0 = ln(detl);
    dev = dev * chol(ipsi)';

    jb = _nseq_Jcb(b,x); /* Jacobian of Fct w.r.t. endo. vars. */

    retp(_ln2pi + jb - 0.5*(ldet0 + sumc((dev.*dev)')));
endp;


/*
**  Jacobian of function w.r.t. endogenous variables
*/

Proc _nseq_Jcb(b,x);
    local i,j,y,f0,dh,x0,x1,jb,abx,sgnx,arg,ret,en;

    ret = zeros(rows(x),1);

    i = 1;
    do until i > rows(x);

        if _nseq_IsJcbProc;
            ret[i] = ln(abs(det(_nseq_EndoJcb(b,x[i,.]))));

        else;
            jb = zeros(rows(_nseq_EndoInd),rows(_nseq_EndoInd));
            f0 = _nseq_Fct(b,x[i,.]);

            x0 = x[i,_nseq_EndoInd];
            x0 = x0 + 1e-2*(x0 == 0);
            abx = abs(x0);
            sgnx = x0./abx;

            dh = (1e-8)*maxc(abx|(1e-2)*ones(1,cols(x0))).*sgnx';

            j = 1;
            do until j > cols(x0);
                x1 = x[i,.];
                en = _nseq_EndoInd[j];
                x1[en] = x1[en] + dh[en];
                jb[j,.] = _nseq_Fct(b,x1);
                j = j + 1;
            endo;
            ret[i] = ln(abs(det((jb-f0)./dh)));

        endif;
        i = i + 1;
    endo;

    retp(ret);
endp;