Maximum Likelihood Estimation with Nonlinear Equality Constraints¶
This GAUSS maximum likelihood example demonstrates the use of CMLMT to estimate parameters of a tobit model with nonlinear equality constraints.
Key example features¶
Usages of data from the file cmlmttobit.dat (included with cmlmt).
- User defined likelihood function,
lprwith four inputs: A PV structure storing parameters.
Additional X and y data matrices, which are passed to
cmlmt`()as optional arguments.The required ind input.
- User defined likelihood function,
The inclusion of analytic gradient computations, as specified in the
lprfunction.A user-defined function
eqpin combination with the c0.eqProc member of thecmlmtControlstructure to specify the equality constraints.
There are two equality constraints that are implemented in this example, one linear and one nonlinear:
where \(p[1], p[2], \ldots, p[5]\) are the model parameters to be estimated.
Code for estimation¶
/*
** Maximum likelihood tobit model
*/
new;
library cmlmt;
// Tobit likelihood function with 4 inputs
// i. p - The PV parameter structur
// ii-iii. x and y - Extra data needed by the objective procedure
// ii. ind - The indicator vector
proc lpr(struct PV p, x, y, ind);
local s2, b0, b, yh, u, res, g1, g2;
// Declare 'mm' to be a modelResults
// struct local to this procedure
struct modelResults mm;
// Unpack parameters from PV structure
b0 = pvUnpack(p, "b0");
b = pvUnpack(p, "b");
s2 = pvUnpack(p, "variance");
// Function computations
yh = b0 + x * b;
res = y - yh;
u = y[., 1] ./= 0;
// If first element of 'ind' is non-zero,
// compute function evaluation
if ind[1];
mm.function = u.*lnpdfmvn(res, s2) + (1 - u).*(ln(cdfnc(yh/sqrt(s2))));
endif;
// If second element of 'ind' is non-zero,
// compute function evaluation
if ind[2];
yh = yh/sqrt(s2);
g1 = ((res~x.*res)/s2) ~ ((res.*res/s2) - 1)/(2*s2);
g2 = ( -( ones(rows(x), 1) ~ x )/sqrt(s2) ) ~ (yh/(2*s2));
g2 = (pdfn(yh)./cdfnc(yh)).*g2;
mm.gradient = u.*g1 + (1 - u).*g2;
endif;
// Return modelResults struct
retp(mm);
endp;
// Pack parameters into PV structure
// note that first call to pvPack
struct PV p0;
p0 = pvPack(pvCreate, 1, "b0");
p0 = pvPack(p0, 1|1|1, "b");
p0 = pvPack(p0, 1, "variance");
// Load data
z = loadd(getGAUSSHome("pkgs/cmlmt/examples/cmlmttobit.dat"));
// Separate X and y
y = z[., 1];
x = z[., 2:4];
// Declare 'c0' to be a cmlmtControl struct
// and fill with default settings
struct cmlmtControl c0;
c0 = cmlmtControlCreate();
// Procedure to compute equality constraints
// this must specify the constraint such that
// eqp(x) = 0
proc eqp(p, x, y);
local c, b0, b;
// Extract parameters
b0 = pvUnpack(p, "b0");
b = pvUnpack(p, "b");
// This will be returned and
// it should be a vector of zeros
// with the same number of rows as constraints
c = zeros(2, 1);
// First constraint
c[1] = b0 - b[1];
// Second constraint
c[2] = b[2] * b[3] - 1;
retp(c);
endp;
// Assign pointer for equality procedure
c0.eqProc = &eqp;
// Declare 'out' to be a cmlmtResults
// struct to hold optimization results
struct cmlmtResults out;
out = cmlmtprt(cmlmt(&lpr, p0, x, y, c0));
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 -129.935
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
---------------------------------------------------------------------
b0[1,1] 0.7560 0.0862 8.775 0.0000 27.6776
b[1,1] 0.7560 0.0862 8.775 0.0000 -27.6779
b[2,1] 1.1077 0.1279 8.658 0.0000 -34.1711
b[3,1] 0.9028 0.1043 8.658 0.0000 -41.9260
variance[1,1] 1.2446 0.1883 6.610 0.0000 0.0085
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.
The results are consistent with our constraints: - b0 and b[1, 1] are equal (\(b0 - b[1, 1] = 0\)). - b[2, 1]*b[3, 1] - 1 = \(1.1077 * 0.9028 - 1 = 3.15e-05\).
The gradients are not equal to zero, which is indicative that the contraints are binding.
Parameter correlations¶
Correlation matrix of the parameters
1 1 -0.27931016 0.27931016 -0.0049885835
1 1 -0.27931016 0.27931016 -0.0049885835
-0.27931016 -0.27931016 1 -1 0.01958035
0.27931016 0.27931016 -1 1 -0.01958035
-0.0049885909 -0.0049885909 0.019580346 -0.019580346 1
Confidence intervals¶
Wald Confidence Limits
0.95 confidence limits
Parameters Estimates Lower Limit Upper Limit Gradient
----------------------------------------------------------------------
b0[1,1] 0.7560 0.5849 0.9270 27.6776
b[1,1] 0.7560 0.5849 0.9270 -27.6779
b[2,1] 1.1077 0.8537 1.3617 -34.1711
b[3,1] 0.9028 0.6958 1.1098 -41.9260
variance[1,1] 1.2446 0.8708 1.6184 0.0085
Number of iterations 10
Minutes to convergence 0.00012