/*********************************************************
Simultaneous Equation Models   (Schoenberg July '95)

Written by Ronald Schoenberg, July 12, 1995.
rons@u.washington.edu

The programs below estimate parameters of linear and nonlinear
simultaneous equation systems.  There are two versions of each,
one requiring CML, Constrained Maximum Likelihood application
module, and another requiring MAXLIK 4.0, the unconstrained
Maximum Likelihood application module.  Any of these programs,
however, could be adapted quite easily for earlier versions of
MAXLIK.

The source code in this file is in the public domain and requires
GAUSS 3.2+.

The author and Aptech Systems cannot be held responsible for any
actual use of the software below.

Current versions of these programs are available by anonymous ftp from
ftp.u.washington.edu on /public/rons.
I have also submitted them to the GAUSS archive at American University.
They will also be made available in the CML and
MAXLIK 4.0 Applications modules as examples.

There are four programs:

CLEQ.GAU  -  constrained maximum likelihood estimation of linear
simultaneous equation models.

CNLEQ.GAU  -  constrained maximum likelihood estimation of nonlinear
simultaneous equation models.

LEQ.GAU  -  unconstrained maximum likelihood estimation of linear
simultaneous equation models.

NLEQ.GAU  -  unconstrained maximum likelihood estimation of nonlinear
simultaneous equation models.


CLEQ.GAU and CNLEQ.GAU require the CML applications module, and LEQ.GAU
and NLEQ.GAU require MAXLIK 4.0.  They all require GAUSS 3.2+, though it
should be straightforward to adapt these programs to earlier versions of
MAXLIK and GAUSS by changing global and procedure names, and deleting
references to procedures not available in those versions.


CLEQ.GAU illustrates constraining the absolute values of the eigenvalues
of the endogenous coefficient matrix.  If the system is to be "dynamically
stable", the absolute values of the eigenvalues of I - this matrix (where
I is a conformable identity matrix) must be on the unit interval.
CLEQ.GAU contains code to implement constraining the parameters to satisfy
this condition.  If there are explicit lagged variables in the model,
similar but more complicated constraints will need to be applied to the
parameters.  Green, in his Econometric Analysis book discussions these types
of conditions.

CNLEQ.GAU illustrates the same type of constraint in the context of the
nonlinear simultaneous equation model.  CNLEQ.GAU estimates the following
model:

    U = f(Y,X,B)

where U is a TxK matrix of residuals, Y a TxK matrix of endogenous
variables, X a TxL matrix of exogenous variables, and B an Mx1 vector of
coefficients.  Thus f() is a set of K simultaneous nonlinear equations.
The "dynamically stable" condition here is with respect to the Jacobian of
the residuals w.r.t. to the endogenous variables.  This has to be computed
by observation if the Jacobian is a function of the observations.
CNLEQ.GAU illustrates this type of constraint.  If there are explict
lagged variables in this model, more complicated constraints would have to
be considered.

CAVEAT:  the constraints in these programs are for illustration. I'm not
really sure that I have it right.   Greene is the only person I've found
who talks about it and even then only in the context of explicit lagged
variables.  I'd appreciate hearing from anyone who has knowledge about
any of this or who would like to correct my discussion above.


All of the programs refered to above also call procedures for generating
profile likelihood, bootstrap, and Bayesian statistical inference.
***********************************************************************\


/********************************************************************
   LEQ.GAU
********************************************************************/

/*
**   Linear Simultaneous Equations
**
**   The source code below, with CML, will estimate the coefficients
**   of a linear system of simultaneous equations:
**
**      Y * Gamma = X * Beta + U
**
**   where Y is an NxK matrix of endogenous variables, X an NxL matrix
**   of exogenous variables, Gamma a KxK coefficient matrix, Beta an LxK
**   coefficient matrix, and U an NxK matrix of residuals.
**  where U is a NxK matrix of residuals and f() is an arbitrary function
**  of Y, a NxK matrix of endogenous variables, X, an NxL matrix of
**  exogenous variables, and B, an Mx1 vector of coefficients.  Thus
**  there are N observations, and f() represents K simultaneous
**  equations.
**
**  The coefficient matrices are set up in GAUSS procs which return
**  the matrices given the vector of parameters.
**
*/

library maxlik;
#include maxlik.ext;
maxset;

_eq_Endo = { 1, 2 };   /* endogenous variables */
_eq_Exo  = { 3, 4 };   /* exogenous variables  */

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

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

proc _eq_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 = "leq";


/*
**  output file
*/

output file = leq.out reset;


/*
**  Likelihood function
*/

_ln2pi = -0.5 * rows(_eq_Endo) * ln(2*pi);

proc lpr(b,z);
    local u, psi, oldt, ldet0, ipsi, beta, gamm, cnst;
    gamm = _eq_Gamma(b);
    beta = _eq_Beta(b);
    cnst = _eq_B0(b);

    u = z[.,_eq_Endo] * gamm - z[.,_eq_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;


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

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


__title = "bootstrap confidence limits";
_max_BootFname = "boot";
 _max_NumSample = 100;
{ b2,f2,g2,cov2,ret2 } = maxboot(dataset,0,&lpr,start);
call maxclprt(b2,f2,g2,maxblimits("boot"),ret2);


__title = "Bayesian confidence limits";
_max_BootFname = "bayes";
_max_NumSample = 500;
_max_MaxTime = 600;

{ b3,f3,g3,cov3,ret3 } = maxbayes(dataset,0,&lpr,start);
call maxclprt(b3,f3,g3,maxblimits("bayes"),ret3);

output off;


/********************************************************************
   CLEQ.GAU
********************************************************************/

/*
**   Constrained Linear Simultaneous Equations
**
**   The source code below, with CML, will estimate the coefficients
**   of a linear system of simultaneous equations:
**
**      Y * Gamma = X * Beta + U
**
**   where Y is an NxK matrix of endogenous variables, X an NxL matrix
**   of exogenous variables, Gamma a KxK coefficient matrix, Beta an LxK
**   coefficient matrix, and U an NxK matrix of residuals.
**  where U is a NxK matrix of residuals and f() is an arbitrary function
**  of Y, a NxK matrix of endogenous variables, X, an NxL matrix of
**  exogenous variables, and B, an Mx1 vector of coefficients.  Thus
**  there are N observations, and f() represents K simultaneous
**  equations.
**
**  The coefficient matrices are set up in GAUSS procs which return
**  the matrices given the vector of parameters.
**
**  This example also illustrates constraining the Gamma coefficient
**  matrix.  For a model to be "dynamically stable" the absolute value
**  of the eigenvalues of I - Gamma must be on the unit interval.  The
**  example below provides such a constraint.
*/

library cml;
#include cml.ext;
cmlset;

_ceq_Endo = { 1, 2 };   /* endogenous variables */
_ceq_Exo  = { 3, 4 };   /* exogenous variables  */

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

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

proc _ceq_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 = "cleq";


proc ineqp(b);
    local ib;
    ib = _ceq_Gamma(b);
    retp(.99-maxc(abs(eig(eye(rows(ib))-ib))));
endp;
_cml_IneqProc = &ineqp;


/*
**  output file
*/

output file = cleq.out reset;


/*
**  Likelihood function
*/

_ln2pi = -0.5 * rows(_ceq_Endo) * ln(2*pi);

proc lpr(b,z);
    local u, psi, oldt, ldet0, ipsi, beta, gamm, cnst;
    gamm = _ceq_Gamma(b);
    beta = _ceq_Beta(b);
    cnst = _ceq_B0(b);

    u = z[.,_ceq_Endo] * gamm - z[.,_ceq_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;


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

__title = "Wald-type confidence limits";

/* cl1 = cmltlimits(b1,cov1);  if there are no inequality constraints */
cl1 = cmlclimits(b1,cov1);  /* if there are inequality constraints */

call cmlclprt(b1,f1,g1,cl1,ret1);


__title = "bootstrap confidence limits";
_cml_BootFname = "boot";
 _cml_NumSample = 100;
{ b2,f2,g2,cov2,ret2 } = cmlboot(dataset,0,&lpr,start);
call cmlclprt(b2,f2,g2,cmlblimits("boot"),ret2);


__title = "Bayesian confidence limits";
_cml_BootFname = "bayes";
_cml_NumSample = 500;
_cml_MaxTime = 600;

{ b3,f3,g3,cov3,ret3 } = cmlbayes(dataset,0,&lpr,start);
call cmlclprt(b3,f3,g3,cmlblimits("bayes"),ret3);

output off;



/********************************************************************
   CNLEQ.GAU
********************************************************************/

/*
**  Constrained Nonlinear Simultaneous Equations
**
**  The source code below, with CML, will estimate the coefficients
**  of a general nonlinear system of simultaneous equations:
**
**      U = f(Y,X,B)
**
**  where U is a NxK matrix of residuals and f() is an arbitrary function
**  of Y, a NxK matrix of endogenous variables, X, an NxL matrix of
**  exogenous variables, and B, an Mx1 vector of coefficients.  Thus
**  there are N observations, and f() represents K simultaneous
**  equations.
**
**  If the Jacobian of residuals w.r.t. the endogenous variables, i.e.
**  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 of solution.
**
**  Required:  a proc to compute the residuals given the coefficients
**             and the data.
**
**  Optional:  a proc to compute the Jacobians of the residuals w.r.t.
**             the endogenous variables, i.e., U w.r.t. Y.
**
**
**  This example also illustrates constraining the Jacobians.  For a
**  model to be "dynamically stable" the absolute value of the
**  eigenvalues of I - the Jacobian of the residuals  w.r.t. the endogenous
**  variables must be on the unit interval.  The example below
**  provides for such constraints.
*/

library cml;
#include cml.ext;
cmlset;


/*
**  Procedure for computing the function
**
**  The system being estimated is a set of two simultaneous equations
**  where in the first equation there is an interaction between the
**  second endogenous variable and the first exogenous variable:
**
**    Y1 = b1 * Y2 + b2 * Y2 * X1 + b3 * X1 + b4 + U1
**    Y2 = b5 * Y1 + b6 * X2 + b7 + U2
**
*/

proc _nl_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 _nl_EndoInd to the indices of the endogenous variables
**  in the dataset.
*/

_nl_EndoInd = { 1, 2 };


/*
**  If Proc _nl_EndoJcb is NOT provided, set _nl_IsJcbProc to zero.
**
**    _nl_IsJcbProc = 0;
*/


/*
**  If Proc _nl_EndoJcb IS provided, set _nl_IsJcbProc = 1
*/

_nl_IsJcbProc = 1;

proc _nl_EndoJcb(b,z);
    local j;
    j = zeros(2,2);
    j[1,1] = 1;
    j[1,2] = -b[1] - b[2] * z[.,3];
    j[2,1] = -b[5];
    j[2,2] = 1;
    retp(j);
endp;




/*
** start values
*/

start = { .2, .0, .2, .2, .2, .2, .2 };


/*
**  dataset
*/

_nl_data = loadd("cnleq");


/*
**  output file
*/

output file = cnleq.out reset;



/*
**  Constraining the eigenvalues of the Jacobians
**
**  The code below illustrates a method for constraining the
**  absolute value of the eigenvalues of I - the Jacobian of the
**  residuals w.r.t. the endogenous variables to a unit interval.
**  If the system of equations being estimated here is a cross-
**  sectional representation of a dynamic system, this constraint
**  ensures that it is a "stable" system, i.e., one in a
**  dynamic equilibrium.
**
**  If there are lagged variables in the system, you might want
**  to constrain instead the absolute values of the eigenvalues of
**  the product of a matrix of the coefficients of the endogenous
**  variables times the matrix of the coefficients of the
**  endogenous on the lagged variables.  See Greene, Econometric
**  Analysis, page 641ff, for details.
*/


proc neqp(b);
    retp(.99-_nl_Jcb(b,_nl_data,2));
endp;

_cml_IneqProc = &neqp;


_ln2pi = -0.5 * rows(_nl_EndoInd) * ln(2*pi);  /* constant needed by the  */
                                               /* log-likelihood          */

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

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


__title = "Wald-type confidence limits";

/* cl1 = cmltlimits(b1,cov1);  if there are no inequality constraints */
cl1 = cmlclimits(b1,cov1);  /* if there are inequality constraints */

call cmlclprt(b1,f1,g1,cl1,ret1);


__title = "profile likelihood confidence limits";

_cml_DirTol = .0001; /* no need for estimates to be accurate   */
                     /* because profile likelihood needs only  */
                     /* an accurate evaluation of the function */

__tol = 1e-4;        /* 4 places of accuracy in c.l.'s */

cl3 = cmlpflclimits(b1,f1,_nl_data,0,&lpr);
call cmlclprt(b1,f1,g1,cl3,ret1);




__title = "bootstrap confidence limits";
_cml_BootFname = "boot";
_cml_NumSample = 100;
_cml_DirTol = 1e-5;  /* we need more accuracy here */

{ b2,f2,g2,cov2,ret2 } = cmlboot(_nl_data,0,&lpr,b1);
call cmlclprt(b2,f2,g2,cmlblimits("boot"),ret2);





__title = "Bayesian confidence limits";
_cml_BootFname = "bayes";
_cml_NumSample = 500;

{ b3,f3,g3,cov3,ret3 } = cmlbayes(_nl_data,0,&lpr,b1);
call cmlclprt(b3,f3,g3,cmlblimits("bayes"),ret3);


output off;




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

/*
**  Likelihood function
*/

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

    dev = _nl_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 = _nl_Jcb(b,x,1); /* 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 _nl_Jcb(b,x,ind);
    local i,j,y,f0,dh,x0,x1,jb,abx,sgnx,arg,ret,en,dt;

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

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

        if _nl_IsJcbProc;

            jb = _nl_EndoJcb(b,x[i,.]);
            if ind == 1;
                ret[i] = ln(abs(det(jb)));
            else;
                ret[i] = maxc(abs(eig(eye(rows(jb))-jb)));
            endif;

        else;

            jb = zeros(rows(_nl_EndoInd),rows(_nl_EndoInd));
            f0 = _nl_Fct(b,x[i,.]);

            x0 = x[i,_nl_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 = _nl_EndoInd[j];
                x1[en] = x1[en] + dh[en];
                jb[j,.] = _nl_Fct(b,x1);
                j = j + 1;
            endo;
            if ind == 1;
                ret[i] = ln(abs(det((jb-f0)./dh)));
            else;
                ret[i] = maxc(abs(eig(eye(rows(jb))-(jb-f0)./dh)));
            endif;
        endif;
        i = i + 1;
    endo;

    retp(ret);
endp;



/********************************************************************
   NLEQ.GAU
********************************************************************/


/*
**  Nonlinear Simultaneous Equations
**
**  The source code below, with CML, will estimate the coefficients
**  of a general nonlinear system of simultaneous equations:
**
**      U = f(Y,X,B)
**
**  where U is a NxK matrix of residuals and f() is an arbitrary function
**  of Y, a NxK matrix of endogenous variables, X, an NxL matrix of
**  exogenous variables, and B, an Mx1 vector of coefficients.  Thus
**  there are N observations, and f() represents K simultaneous
**  equations.
**
**  If the Jacobian of residuals w.r.t. the endogenous variables, i.e.
**  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 of solution.
**
**  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.
**
*/

library maxlik;
#include maxlik.ext;
maxset;


/*
**  Procedure for computing the function
**
**  The system being estimated is a set of two simultaneous equations
**  where in the first equation there is an interaction between the
**  second endogenous variable and the first exogenous variable:
**
**    Y1 = b1 * Y2 + b2 * Y2 * X1 + b3 * X1 + b4 + U1
**    Y2 = b5 * Y1 + b6 * X2 + b7 + U2
**
*/

proc _ml_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 _ml_EndoInd to the indices of the endogenous variables
**  in the dataset.
*/

_ml_EndoInd = { 1, 2 };


/*
**  If Proc _ml_EndoJcb is NOT provided, set _ml_IsJcbProc to zero.
**
**    _ml_IsJcbProc = 0;
*/


/*
**  If Proc _ml_EndoJcb IS provided, set _ml_IsJcbProc = 1
*/

_ml_IsJcbProc = 1;

proc _ml_EndoJcb(b,z);
    local j;
    j = zeros(2,2);
    j[1,1] = 1;
    j[1,2] = -b[1] - b[2] * z[.,3];
    j[2,1] = -b[5];
    j[2,2] = 1;
    retp(j);
endp;




/*
** start values
*/

start = { .2, .0, .2, .2, .2, .2, .2 };


/*
**  dataset
*/

_ml_data = loadd("nleq");


/*
**  output file
*/

output file = nleq.out reset;


_ln2pi = -0.5 * rows(_ml_EndoInd) * ln(2*pi);  /* constant needed by the  */
                                               /* log-likelihood          */

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

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



__title = "Wald-type confidence limits";

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



__title = "profile likelihood confidence limits";

_max_DirTol = .0001; /* no need for estimates to be accurate   */
                     /* because profile likelihood needs only  */
                     /* an accurate evaluation of the function */

__tol = 1e-4;        /* 4 places of accuracy in c.l.'s */

cl3 = maxpflclimits(b1,f1,_ml_data,0,&lpr);
call maxclprt(b1,f1,g1,cl3,ret1);





__title = "bootstrap confidence limits";
_max_BootFname = "boot";
_max_NumSample = 100;
_max_DirTol = 1e-5;  /* we need more accuracy here */

{ b2,f2,g2,cov2,ret2 } = maxboot(_ml_data,0,&lpr,start);
call maxclprt(b2,f2,g2,maxblimits("boot"),ret2);





__title = "Bayesian confidence limits";
_max_BootFname = "bayes";
_max_NumSample = 500;

{ b3,f3,g3,cov3,ret3 } = maxbayes(_ml_data,0,&lpr,start);
call maxclprt(b3,f3,g3,maxblimits("bayes"),ret3);


output off;




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

/*
**  Likelihood function
*/

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

    dev = _ml_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 = _ml_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 _ml_Jcb(b,x);
    local i,j,y,f0,dh,x0,x1,jb,abx,sgnx,arg,ret,en,dt;

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

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

        if _ml_IsJcbProc;

            jb = _ml_EndoJcb(b,x[i,.]);
            ret[i] = ln(abs(det(jb)));

        else;

            jb = zeros(rows(_ml_EndoInd),rows(_ml_EndoInd));
            f0 = _ml_Fct(b,x[i,.]);

            x0 = x[i,_ml_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 = _ml_EndoInd[j];
                x1[en] = x1[en] + dh[en];
                jb[j,.] = _ml_Fct(b,x1);
                j = j + 1;
            endo;
            ret[i] = ln(abs(det((jb-f0)./dh)));
        endif;
        i = i + 1;
    endo;

    retp(ret);
endp;



/********************************************************************
   CNLEQ.SIM
********************************************************************/

/*
**  Simulate a nonlinear simultaneous equation system
*/

nobs = 100;
rndseed 64776897;


cr = { 1 .3, .3 1 };   /* correlation matrix of exogenous variables */



b = { .4, .1, .4 };   /* endogenous coefficients */

y = zeros(nobs,2);
X = rndn(nobs,2) * chol(cr)';     /* exogenous variables */

x[1:10,1] = 3*x[1:10,1];

g = { .2   0,        /* exogenous coefficient matrix */
       0  .2 };

c0 = { .2  .2 };     /* constants */

ib = eye(2);
i = 1;
do until i > nobs;

    ib[1,2] =  -b[1] - b[2]*x[i,1];
    ib[2,1] = -b[3];

    y[i,.] = (x[i,.] * g + c0 + rndn(1,2)) * inv(ib);

    i = i + 1;
endo;

dat = y~x;

lbl = { Y1, Y2, X1, X2 };

call saved(dat,"cnleq",lbl);



/********************************************************************
   LEQ.SIM
********************************************************************/

/*
**  Simulate a simultaneous equation system
*/

nobs = 500;
rndseed 64776897;

g0 = { 0 .2, .2 0 };  /* endogenous coefficient matrix */
b0 = { .3 0, 0 .3 };  /* exogenous coefficient matrix */
c0 = { 1 1 };         /* constants */

cr = {  1  .3,    /* correlation matrix of exogenous variables */
       .3   1 };

X = rndn(nobs,2) * chol(cr);

Y = ( X * b0 + c0  + rndn(nobs,2) ) * inv(eye(2) - g0);

dat = Y~X;

lbl = { Y1, Y2, X1, X2 };

call saved(dat,"leq",lbl);



/********************************************************************
   CLEQ.SIM
********************************************************************/


/*
**  Simulate a simultaneous equation system
*/

nobs = 100;
rndseed 64776897;

g0 = { 0 1.2, 1.2 0 };  /* endogenous coefficient matrix */

b0 = { .3 0, 0 .3 };  /* exogenous coefficient matrix */
c0 = { 1 1 };         /* constants */

cr = {  1  .3,    /* correlation matrix of exogenous variables */
       .3   1 };

X = rndn(nobs,2) * chol(cr);

Y = ( X * b0 + c0  + rndn(nobs,2) ) * inv(eye(2) - g0);

dat = Y~X;

lbl = { Y1, Y2, X1, X2 };

call saved(dat,"cleq",lbl);




/********************************************************************
   NLEQ.SIM
********************************************************************/

/*
**  Simulate a nonlinear simultaneous equation system
*/

nobs = 100;
rndseed 64776897;


cr = { 1 .6, .6 1 };   /* correlation matrix of exogenous variables */



b = { .3, .6, .3 };   /* endogenous coefficients */

y = zeros(nobs,2);
X = rndn(nobs,2) * chol(cr)';     /* exogenous variables */

g = { .5   0,        /* exogenous coefficient matrix */
       0  .5 };

c0 = { .5  .5 };     /* constants */

ib = eye(2);
i = 1;
do until i > nobs;

    ib[1,2] =  -b[1] - b[2]*x[i,2];
    ib[2,1] = -b[3];

    y[i,.] = (x[i,.] * g + c0 + rndn(1,2)) * inv(ib);

    i = i + 1;
endo;

dat = y~x;

lbl = { Y1, Y2, X1, X2 };

call saved(dat,"nleq",lbl);

