# Maximum Likelihood Nonlinear Simultaneous Equation Model Estimation#

The source code below estimates the coefficients of a general nonlinear system of simultaneous equations:

$U = f(Y, X, B)$

where $$U$$ is a $$N \times K$$ matrix of residuals and $$f()$$ is an arbitrary function of $$Y$$, a $$N \times K$$ matrix of endogenous variables, $$X$$, an $$N \times L$$ matrix of exogenous variables, and $$B$$, an $$M \times 1$$ vector of coefficients. Thus there are $$N$$ observations, and $$f()$$ represents $$K$$ simultaneous equations.

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.

In this example, a procedure is provided for computing the Jacobian of the function with respect to the endogenous variables. When this proc, _nseq_EndoJcb, is provided, time of computation is reduced significantly.

## Code for estimation#

After loading the library, the first step is to set up the system of nonlinear simultaneous equations:

$Y_1 = b_1 Y_2 + c_1 Y_2 * X_1 + g_1 X_1 + d_1 + U_1$
$Y_2 = b_2 Y_1 + g_2 X_2 + d_2 + U_2$

where $$Y_1$$ and $$Y_2$$ are endogenous, $$X_1$$ and $$X_2$$ are exogenous, and $$U_1$$ and $$U_2$$ are residuals. The objective is to estimate the parameter vector $$\theta = {b_1, c_1, g_1, d_1, b_2, g_2, d_2}$$ using maximum likelihood.

new;
library cmlmt;

// Procedure for computing the residuals from
// the nonlinear simultaneous equations.
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;

// Procedure used in computing the Jacobians of the function
proc _nseq_EndoJcb(b, xx);
local j;
j = zeros(2, 2);

j[1, 1] = 1;
j[1, 2] = -b[1] - b[2] * xx[., 3];
j[2, 1] = -b[5];
j[2, 2] = 1;

retp(j);
endp;


Next, we’ll set up some basic controls for computation of the nonlinear system:

/*
**  Set endoIndices to the indices of the endogenous variables
**  in the dataset.
*/
endoIndices = { 1, 2 };

/*
**  If Proc _nseq_EndoJcb is NOT provided, set endoJcb to zero.
**  else if Proc _nseq_EndoJcb IS provided, set endoJcb to 1.
*/

endoJcb = 1;


Next, we’ll provide the likelihood function and the computational Jacobian function:

/*
**  Likelihood function
*/

proc lpr(b, z, endoIndices, endoJcb, ind);

local dev, psi_, g0, oldt, ldet0, ipsi, jb;

dev = _nseq_Fct(b, z);

psi_ = dev'dev/rows(z);

// Set trap to check for invertibility
oldt = trapchk(1);
trap 1, 1;
ipsi = invpd(psi_);
trap oldt, 1;

// Return error if non-invertible
if scalmiss(ipsi);
retp(error(0));
endif;

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

jb = _nseq_Jcb(b, z, endoIndices, endoJcb, 1); /* absolute value of log det of Jacobian */
/* of Fct w.r.t. endo. vars. */

struct modelResults mm;
mm.function = _ln2pi + jb - 0.5*(ldet0 + sumc((dev.*dev)'));

retp(mm);
endp;

/*
**  Jacobian of function w.r.t. endogenous variables
**
**  ind = 1, compute log of the absolute value of the determinant
**  ind = 2, compute maximum absolute eigenvalue
*/

proc _nseq_Jcb(b,  z, endoIndices, endoJcb, ind);
local i, j, f0, dh, x0, jb, abx, sgnx, ret, en, xx;

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

i = 1;
do until i > rows(z);
xx = z[i,.];
if endoJcb == 1;
if ind == 1;
ret[i] = ln(abs(det(_nseq_EndoJcb(b, xx))));
else;
jb = _nseq_EndoJcb(b,xx);
ret[i] = maxc(abs(eig(eye(rows(jb)) - jb)));
endif;
else;
jb = zeros(rows(endoIndices), rows(endoIndices));
f0 = _nseq_Fct(b, xx);

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

endif;
i = i + 1;
endo;

retp(ret);
endp;


Now, we’re ready to set up the cmlmt() optimization:

// Load the data

// Vector of starting values for our 7 parameters
b0 = { .2, .0, .2, .2, .2, .2, .2 };

// Declare 'c0' to be a cmlmtControl struct
// and fill with default settings
struct cmlmtControl c0;
c0 = cmlmtControlCreate();

/*
**   Convergence can sometimes fail because a feasible line search can't
**   be found.  Thus the line search feasibility test is turned off.
*/
c0.feasibleTest = 0;

/*
**  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 will need
**  a more complicated constraint.  See Greene, Econometric
**  Analysis, page 641ff, for details.
*/
proc neqp(b, z, endoIndices, endoJcb);
retp(.99-_nseq_Jcb(b, z, endoIndices, endoJcb, 2));
endp;

c0.ineqProc = &neqp;

// Constant used in log-likelihood procedure
// pulled out of procedure for computational efficiency
_ln2pi = -0.5 * rows(endoIndices) * ln(2*pi);

// Declare 'out' to be a cmlmtResults
// struct to hold optimization results
struct cmlmtResults out;
out = cmlmt(&lpr, b0, z, endoIndices, endoJcb, c0);

// Print results
call cmlmtprt(out);


## Results#

The cmlmtprt() procedure prints three output tables:

• Estimation results.

• Correlation matrix of parameters.

• Wald confidence limits.

### Estimation results#

===============================================================================
CMLMT Version 3.0.0
===============================================================================

return code =    0
normal convergence

Log-likelihood        -265.265
Number of cases     100

Covariance of the parameters computed by the following method:
ML covariance matrix
Parameters    Estimates     Std. err.  Est./s.e.  Prob.    Gradient
---------------------------------------------------------------------
x[1,1]   -0.1224        0.4956      -0.247   0.8050      0.0000
x[2,1]    0.1176        0.0863       1.362   0.1731      0.0000
x[3,1]    0.6036        0.3093       1.951   0.0510      0.0000
x[4,1]    0.6122        0.3734       1.640   0.1011      0.0000
x[5,1]    0.7219        0.1851       3.899   0.0001      0.0000
x[6,1]    0.2951        0.1097       2.690   0.0071     -0.0011
x[7,1]    0.3362        0.1348       2.495   0.0126     -0.0001


The estimation results reports:

• That the model has converged normally with a return code of 0. Any return code other than 0, indicates an issue with convergence. The cmlmt() documentation provides details on how to interpret non-zero return codes.

• The log-likelihood value and number of cases.

• Parameter estimates, standard errors, t-statistics and associated p-values, and gradients.

### Parameter correlations#

Correlation matrix of the parameters
1      -0.50689926      -0.92087721      -0.95115909       0.22738692      -0.33543423       -0.1759722
-0.50689926                1       0.32167663       0.38002565     -0.079802709       0.13764327       0.06266076
-0.92087721       0.32167663                1       0.90688866       -0.3135579       0.31601986       0.23706339
-0.95115909       0.38002565       0.90688866                1      -0.22985309       0.32358823        0.1165813
0.22738692     -0.079802709       -0.3135579      -0.22985309                1      -0.59469942      -0.73925046
-0.33543423       0.13764327       0.31601986       0.32358823      -0.59469942                1       0.46284909
-0.1759722       0.06266076       0.23706339        0.1165813      -0.73925046       0.46284909                1


### Confidence intervals#

Wald Confidence Limits

0.95 confidence limits
Parameters    Estimates     Lower Limit   Upper Limit    Gradient
----------------------------------------------------------------------
x[1,1]          -0.1224       -1.1065        0.8618        0.0000
x[2,1]           0.1176       -0.0538        0.2889        0.0000
x[3,1]           0.6036       -0.0106        1.2178        0.0000
x[4,1]           0.6122       -0.1293        1.3537        0.0000
x[5,1]           0.7219        0.3542        1.0895        0.0000
x[6,1]           0.2951        0.0773        0.5130       -0.0011
x[7,1]           0.3362        0.0686        0.6038       -0.0001

Number of iterations    9
Minutes to convergence     0.00292