new;
/* C17LIML.gss (Example 17.3.2)
**
** To demonstrate an application of the LIML estimator, we provide the GAUSS
** program C17LIML.gss in the examples manual. The method of estimation
** utilized in the program is based on Equations (17.3.10) and (17.3.11) and the
** discussion in the preceding paragraph. Following Example 17.3.1, GAUSS
** generates a random sample consistent with a version of the data sampling
** process in Equation (17.3.3) and then computes the LIML parameter estimates
** and the associated estimated standard errors. For comparison purposes,
** GAUSS also calculates and reports the 2SLS estimates for the same set of
** observations. Finally, the program conducts a hypothesis test and reports a
** confidence region for a subset of the model parameters. By increasing the
** sample size (n) used in the analysis, the user should find that the LIML and
** 2SLS estimates of model parameters converge to the true parameter values.
*/
/* Set up some graphics and and program parameters */
cls;
graphset;
/* Call the graphics and EF libraries */
library pgraph, ef;
begin:
print "LIML Estimation of the Parameters in the Linear Model";
print;
print " The purpose of this program is to demonstrate an application of the LIML";
print "estimator in the context of an over-determined simultaneous system with";
print "contemporaneously correlated errors. The observations are generated from a";
print "three-equation system defined as follows:";
print;
print " Y*Gama + X*B + E = 0";
print;
print "This system is similar, but not identical, to empirical models discussed in Tsurumi, H.,";
print "1990, ''Comparing Bayesian and Non-Bayesian Limited Information Estimators,''";
print "S. Geisser, J. S. Hodges, S. J. Press, and A. Zellner (ed.), Bayesian and Likelihood";
print "Methods in Statistics and Econometrics, North Holland Publishing, Amsterdam.";
print;
print " The objective will be to estimate parameters (and standard errors) for a single";
print "equation using LIML and 2SLS. The (3 x 1) noise vectors are generated iid from a";
print "bivariate normal distribution with a zero mean vector, and the user can choose the";
print "form of the contemporaneous covariance matrix of the noise elements. The matrix of";
print "instrumental variables for this exercise is defined as Z = X . The first column of X is";
print "a vector of ones and the remaining non-constant columns of X are generated iid from";
print "a uniform probability distribution on the interval [0,10]. The user will choose the sample";
print "size (n). Finally, the program conducts a hypothesis test and reports a confidence region";
print "for a subset of the model parameters.";
print;
print "Press any key to continue...";
waitc;
cls;
ncode:
print "What sample size (n) do you want to use for the demonstration?";
print " Please choose an integer in the set [25, 750]";
print;
n = con(1,1);
if n < 25 or n > 750;
cls;
print "The requested sample size is outside the stated range";
goto ncode;
elseif not n/int(n) == 1;
cls;
print "The requested sample size must be an integer in the set {25, ..., 750}";
goto ncode;
endif;
cls;
/* Specify the default noise covariance matrix */
let sigma[3,3] = 1.000 -1.000 -0.125
-1.000 4.000 0.0625
-0.125 0.0625 4.000;
SigmaInput:
print;
print " For your convenience, a default contemporaneous covariance matrix";
print "is supplied and is defined as:";
print;
print " let sigma[3,3] = 1.000 -1.0000 -0.1250";
print " -1.000 4.0000 0.0625";
print " -0.125 0.0625 4.0000;";
print;
print "Do you want to choose a DIFFERENT contemporaneous covariance matrix?";
print " Please answer Yes or No.";
print;
ans = cons;
if (ans $=="Yes") or (ans $=="yes") or (ans $=="Y") or (ans $=="y");
cls;
else;
cls;
goto endcor;
endif;
variances:
print "Please choose the values of the variances for the three noise elements:";
print " Choose values for sigma[1,1], sigma[2,2], and sigma[3,3]";
print " from the interval [0.1, 5]";
vars = con(3,1);
if vars[1] < 0.1 or vars[2] < 0.1 or vars[3] < 0.1 or vars[1] > 5 or vars[2] > 5 or vars[3] > 5;
cls;
print "At least one of the selected variance parameters is out of range";
goto variances;
endif;
correlation:
print;
print "Please choose the values of the CORRELATIONs between the three noise elements.";
print "For the sake of illustration, choose values for rho[1,2], rho[1,3], and rho[2,3]";
print "from the interval [-0.9, 0.9] to promote (but not ensure) numerical stability.";
print;
print "NOTE: The matrix you define must be positive definite to be a valid";
print "covariance matrix. GAUSS will verify this property before continuing.";
cor = con(3,1);
if cor[1] < -0.9 or cor[2] < -0.9 or cor[3] < -0.9 or cor[1] > 0.9 or cor[2] > 0.9 or cor[3] > 0.9;
cls;
print "At least one of the correlation coefficients is out of range";
goto correlation;
endif;
cls;
/* Form the specified covariance matrix */
sigma[1,1] = vars[1];
sigma[2,2] = vars[2];
sigma[3,3] = vars[3];
sigma[1,2] = cor[1]*sqrt(vars[1]*vars[2]);
sigma[1,3] = cor[2]*sqrt(vars[1]*vars[3]);
sigma[2,3] = cor[3]*sqrt(vars[2]*vars[3]);
sigma[2,1] = sigma[1,2];
sigma[3,1] = sigma[1,3];
sigma[3,2] = sigma[2,3];
/* Check to see that the user defined a valid, i.e., positive definite, Sigma matrix.
** If he/she did not, ask him/her to choose another one. */
Eigenvalues = eigrs(Sigma);
if Eigenvalues > 0;
goto passed;
else;
print "Your choice of Sigma matrix is NOT positive definite. Note the";
print "eigenvalues of the Sigma matrix you defined are given by:";
print;
print " " Eigenvalues';
print;
print "Please choose a POSITIVE DEFINITE Sigma matrix...";
print;
goto SigmaInput;
endif;
passed:
format /mat /rd 15,3;
print;
print " Here is the user-defined contemporaneous covariance matrix:";
print sigma;
print;
print "Do you want to choose a DIFFERENT contemporaneous covariance matrix?";
print " Please answer Yes or No";
ans = cons;
if (ans $=="Yes") or (ans $=="yes") or (ans $=="Y") or (ans $=="y");
cls;
goto variances;
endif;
endcor:
cls;
print;
print " We will now generate the sample data for the three-equation system using";
print "the sample size and contemporaneous covariance matrix specified by the user.";
print "The computer implementation of the data sampling process can be efficiently";
print "programmed in GAUSS code using the following statements:";
print;
print " noise = rndn(n,3)*chol(sigma);";
print ;
print " let Gama[3,3] = -1.000 0.267 0.087";
print " 0.222 -1.000 0.000";
print " 0.000 0.046 -1.000;";
print " invGama = inv(Gama);";
print;
print " let B[7,3] = 6.20 4.40 4.00";
print " 0.00 0.74 0.00";
print " 0.70 0.00 0.53";
print " 0.00 0.00 0.11";
print " 0.96 0.13 0.00";
print " 0.00 0.00 0.56";
print " 0.06 0.00 0.00";
print;
print " X = ones(n,1)~rndu(n,7);";
print " Y = X*B*invGama + noise*invGama;";
print " Z = X;";
print;
print "Press any key to continue...";
waitc;
cls;
/* Specify values for the structural coefficients */
let Gama[3,3] = -1.000 0.267 0.087
0.222 -1.000 0.000
0.000 0.046 -1.000;
let B[7,3] = 6.20 4.40 4.00
0.00 0.74 0.00
0.70 0.00 0.53
0.00 0.00 0.11
0.96 0.13 0.00
0.00 0.00 0.56
0.06 0.00 0.00;
btrue = Gama[2,1]|B[1,1]|B[3,1]|B[5,1]|B[7,1];
btrue = btrue|Gama[1,2]|Gama[3,2]|B[1,2]|B[2,2]|B[5,2];
btrue = btrue|Gama[1,3]|B[1,3]|B[3,3]|B[4,3]|B[6,3];
pii = (-B*inv(Gama));
pii1 = pii[.,2];
pii2 = pii[.,1]~pii[.,3];
pii3 = pii[.,1];
/* Generate the observations */
again:
x = ones(n,1)~rndn(n,6);
if rank(x'x) < cols(x'x);
goto again;
endif;
xmat = x;
x1 = x[.,1]~x[.,3]~x[.,5]~x[.,7];
x2 = x[.,1]~x[.,2]~x[.,5];
x3 = x[.,1]~x[.,3]~x[.,4]~x[.,6];
covx = sigma;
sqrtcovx = chol(covx);
estar = rndn(n,3)*sqrtcovx;
rfestar = -estar*inv(Gama);
ymat = -x*B*inv(Gama) - estar*inv(Gama);
emat = -ymat*Gama - x*B;
print;
print " We now attempt to generate LIML estimates of the parameters of the FIRST";
print "equation in the three-equation system, along with an estimate of the asymptotic";
print "covariance matrix of the estimator. For comparison purposes, we also calculate ";
print "estimates of the model parameters and the covariance matrix for 2SLS estimates.";
print;
print " Note the algorithm used to compute the LIML estimates is the non-gradient";
print "based Nelder-Meade algorithm (in the EF procs directory). This is a sufficiently";
print "complicated nonlinear optimization problem that, for extreme choices of Sigma";
print "and small sample sizes, solving the problem with gradient-based algorithms can";
print "be very challenging without very good starting values and methods for preventing";
print "matrix singularities during the iteration process. The Nelder-Meade approach,";
print "though not fool-proof, is reasonably stable in finding solutions.";
print;
print " If you are interested in further details, we suggest that the user view the";
print "GAUSS code for this example to better understand the data sampling process";
print "(DSP) and the programming of the estimation process.";
print;
print "Press any key to generate the sample data and continue...";
waitc;
cls;
/* First, compute the 2SLS estimates */
y = vec(ymat);
/* Compute the reduced-form estimates by LS */
b1 = olsqr(ymat[.,1], xmat);
b2 = olsqr(ymat[.,2], xmat);
b3 = olsqr(ymat[.,3], xmat);
/* Compute the predicted values */
yh1 = xmat*b1;
yh2 = xmat*b2;
yh3 = xmat*b3;
z1 = yh2~x1;
z2 = yh1~yh3~x2;
z3 = yh1~x3;
/* Compute the 2SLS estimates */
b2s1 = olsqr(ymat[.,1], z1);
b2s2 = olsqr(ymat[.,2], z2);
b2s3 = olsqr(ymat[.,3], z3);
b2sls = b2s1|b2s2|b2s3;
/* Compute the 2SLS predicted values */
zz1 = ymat[.,2]~x1;
zz2 = ymat[.,1]~ymat[.,3]~x2;
zz3 = ymat[.,1]~x3;
/* Estimate the contemporaneous covariance matrix from the 2SLS residuals */
sighat = zeros(3,3);
sighat[1,1] = (ymat[.,1] - zz1*b2s1)'(ymat[.,1] - zz1*b2s1)/n;
sighat[1,2] = (ymat[.,1] - zz1*b2s1)'(ymat[.,2] - zz2*b2s2)/n;
sighat[1,3] = (ymat[.,1] - zz1*b2s1)'(ymat[.,3] - zz3*b2s3)/n;
sighat[2,1] = (ymat[.,2] - zz2*b2s2)'(ymat[.,1] - zz1*b2s1)/n;
sighat[2,2] = (ymat[.,2] - zz2*b2s2)'(ymat[.,2] - zz2*b2s2)/n;
sighat[2,3] = (ymat[.,2] - zz2*b2s2)'(ymat[.,3] - zz3*b2s3)/n;
sighat[3,1] = (ymat[.,3] - zz3*b2s3)'(ymat[.,1] - zz1*b2s1)/n;
sighat[3,2] = (ymat[.,3] - zz3*b2s3)'(ymat[.,2] - zz2*b2s2)/n;
sighat[3,3] = (ymat[.,3] - zz3*b2s3)'(ymat[.,3] - zz3*b2s3)/n;
/* Compute the 2SLS covariance matrix in a way that saves storage space */
temp = inv(sighat);
block11 = temp[1,1]*z1'x*inv(x'x)*x'z1;
block12 = temp[1,2]*z1'x*inv(x'x)*x'z2;
block13 = temp[1,3]*z1'x*inv(x'x)*x'z3;
block21 = temp[2,1]*z2'x*inv(x'x)*x'z1;
block22 = temp[2,2]*z2'x*inv(x'x)*x'z2;
block23 = temp[2,3]*z2'x*inv(x'x)*x'z3;
block31 = temp[3,1]*z3'x*inv(x'x)*x'z1;
block32 = temp[3,2]*z3'x*inv(x'x)*x'z2;
block33 = temp[3,3]*z3'x*inv(x'x)*x'z3;
blockx = (block11~block12~block13)|(block21~block22~block23)|(block31~block32~block33);
cov2sls = inv(blockx);
s2sls = sqrt(diag(cov2sls));
/* Report the 2SLS estimates to the user */
format /mat /rd 16,5;
print "2SLS Parameter Estimates for Equation 1";
print;
print " 2SLS Standard Asymptotic";
print " True value Estimate Error Z-statistic";
btrue[1:5]~b2sls[1:5]~s2sls[1:5]~b2sls[1:5]./s2sls[1:5];
print;
print "Press any key to continue ...";
waitc;
cls;
print;
print "Computational NOTE: It is possible, as in all nonlinear estimation attempts,";
print "that the maximum likelihood algorithm will fail for a number of reasons,";
print "including failure to converge in the number of iterations allowed, or a";
print "matrix may be singular at non-optimal parameter values (i.e., during the";
print "search process). In these cases, rerun the example with a new data set.";
print "The starting values in this example are fixed. In practice, one might also";
print "try different starting values when a problem occurs.";
print;
print " If the program stops with an error message, simply close the GAUSS";
print "window and rerun the example.";
print;
print "Press any key to begin the LIML estimation attempt...";
waitc;
print;
print "GAUSS is now computing the LIML estimates for Equation 1 ... please be patient";
/* Initiate the starting values */
start = b2s1|b2;
zmat = ymat~xmat;
/* Run the Nelder-Meade minimization algorithm to try to find the FIML estimates.
** The ofn() procedure (at the end of the program file) links the objective function
** to the function that the Nelder-Meade algorithm attempts to minimize. See the
** Nelder-Meade code in the procs directory if more details are desired. */
neldtol = 1e-7;
bliml = nelder(start,&logliml,20000,0,neldtol,0);
cls;
covliml = asycov(bliml);
sliml = sqrt(diag(covliml));
/* Report the LIML estimates to the user */
format /mat /rd 16,5;
print "LIML Parameter Estimates for Equation 1";
print;
print " LIML Standard Asymptotic";
print " True value Estimate Error Z-statistic";
btrue[1:5]~bliml[1:5]~sliml[1:5]~bliml[1:5]./sliml[1:5];
print;
print "Press any key to continue ...";
waitc;
cls;
print;
print "Comparison of Coefficient Estimates for Equation 1";
print;
print " 2SLS LIML";
print b2s1~bliml[1:5];
format /mat /rd 5,0;
print;
print "Sample size, n = " n;
print;
print "Press any key to continue ...";
waitc;
cls;
format /mat /rd 16,5;
print;
print "Comparison of Standard Error Estimates for Equation 1";
print;
print " 2SLS LIML";
print s2sls[1:5]~sliml;
format /mat /rd 5,0;
print;
print "Sample size, n = " n;
print;
print "Press any key to continue ...";
waitc;
cls;
r = zeros(1,15);
r[1,13] = 1;
r[1,14] = 1;
r[1,15] = 1;
t11 = (b2sls[2] - 6)/s2sls[2];
t11pvalue = cdfn(t11);
t12 = b2sls[5]/s2sls[5];
t12pvalue = cdfn(t12);
t21 = (bliml[2] - 6)/sliml[2];
t21pvalue = cdfn(t21);
t22 = bliml[5]/sliml[5];
t22pvalue = cdfn(t22);
format /mat /rd 7,4;
print;
print " Next, we want to test the joint one-sided null hypothesis for the first equation,";
print "H0 : B[1,1] >= 6 and B[7,1] >= 0, which is true (the parameters are 6.2 and 0.06).";
print "Using the large sample results for the estimators, we compute the one-sided";
print "Z-test statistics and the associated p-values for this hypothesis:";
print;
print " B[1,1] B[7,1]";
print " Z-statistic P-value Z-statistic P-value";
print;
print " 2SLS " t11 " " t11pvalue " " t12 " " t12pvalue;
print " LIML " t21 " " t21pvalue " " t22 " " t22pvalue;
print;
print " Based on the reported p-values, would you reject the joint null hypothesis due";
print "to one or both components? Would you commit a Type II error in either case? How";
print "could you use the Bonferroni Inequality to control the overall size of the joint test?";
print;
print "Press any key to continue ...";
waitc;
cls;
print;
print " For demonstration purposes, we consider the problem of forming a 90%";
print "confidence region for the LIML and 2SLS model parameters, B[1,1] and B[3,1]";
print "of the first equation in the system. We use the CRegion() procedure in the EF";
print "library to derive the boundary of the joint confidence region for the parameter";
print "estimates, and the region is plotted in a two dimensional figure. ";
print;
print "Press any key to view the confidence regions ...";
waitc;
cls;
/* Form the required weight matrix for the call to CRegion() */
let crg[2,5] = 0 1 0 0 0 0 0 1 0 0;
/* Use the CRegion() procedure to compute the joint region */
ellipse1 = CRegion(crg, bliml[1:5], covliml[1:5,1:5], 0.9, 80);
ellipse2 = CRegion(crg, b2sls[1:5], cov2sls[1:5,1:5], 0.9, 80);
/* Place a large red dot at the true values of the parameters */
s = zeros(3,7);
s[1,1] = 6.2;
s[1,2] = 0.7;
s[1,3] = 8;
s[1,4] = 3.5;
s[1,5] = 4;
s[1,6] = 1;
s[1,7] = 10;
s[2,1] = bliml[2,1];
s[2,2] = bliml[3,1];
s[2,3] = 8;
s[2,4] = 3.5;
s[2,5] = 2;
s[2,6] = 1;
s[2,7] = 10;
s[3,1] = b2sls[2,1];
s[3,2] = b2sls[3,1];
s[3,3] = 8;
s[3,4] = 3.5;
s[3,5] = 14;
s[3,6] = 1;
s[3,7] = 10;
_psym = s;
/* Defining message control values here */
q = zeros(3,7);
q[1,1] = 2;
q[1,2] = 5.8;
q[1,3] = 0.2;
q[1,5] = 2;
q[1,6] = 2;
q[2,1] = 4;
q[2,2] = 5.8;
q[2,3] = 0.2;
q[2,5] = 2;
q[2,6] = 14;
q[3,1] = 6;
q[3,2] = 5.8;
q[3,3] = 0.2;
q[3,5] = 2;
q[3,6] = 4;
_pmsgctl = q;
_pmsgstr = "BhatLIML\000Bhat2SLS\000True Beta";
_plegctl = {2 5 6.9 1.1};
_plegstr = "LIML\0002SLS";
_plwidth = 5;
_ptitlht = 0.2;
_paxht = 0.25;
_pcolor = {2 14};
title("Wald Confidence Ellipse: 90% Coverage");
fonts("simplex simgrma");
xlabel("\202\098\201(1,1)");
ylabel("\202\098\201(1,3)");
xy(ellipse1[1,.]'~ellipse2[1,.]', ellipse1[2,.]'~ellipse2[2,.]');
locate(22,1);
print "The confidence regions for the 2SLS and LIML estimates are plotted in GREEN";
print "and YELLOW, and the estimates are indicated with large dots of the same color.";
print "The large RED dot (if visible) indicates the true parameter values.";
print;
print " Please close the plots and then press any key to continue...";
waitc;
cls;
print;
print "Do you want to repeat the exercise with a different sample";
print " size or contemporaneous covariance matrix?";
print " Please answer Yes or No";
ans = cons;
if (ans $=="Yes") or (ans $=="yes") or (ans $=="Y") or (ans $=="y");
cls;
goto begin;
endif;
cls;
ndpclex;
system;
end;
/* Procedure for the concentrated log-likelihood function */
proc logliml(bs);
local q,LL,sighat, resid1, resid2, resid, invsig, beta, gamm,gamm1,pii;
gamm1 = (-1|bs[1]);
gamm = gamm1~(0|(-1));
q = cols(gamm);
beta = bs[2]|0|bs[3]|0|bs[4]|0|bs[5];
pii = bs[6:12];
resid1 = -ymat[.,1:2]*gamm1 - xmat*beta;
resid2 = ymat[.,2] - xmat*pii;
resid = resid1~resid2;
sighat = resid'resid/n;
LL = -q*0.5*ln(2*pi) - 0.5*ln(det(sighat)) + ln(abs(det(gamm)));
retp(-LL);
endp;
/* Procedure for the estimated asymptotic covariance matrix */
proc asycov(b);
local sighat, resid1, beta, gamm1,z, cov,yhat,pii;
gamm1 = (-1|b[1]);
beta = b[2]|0|b[3]|0|b[4]|0|b[5];
pii = b[6:12];
resid1 = -ymat[.,1:2]*gamm1 - xmat*beta;
sighat = resid1'resid1/n;
yhat = xmat*pii;
z = yhat~x1;
cov = sighat*((inv(z'*x*inv(x'x)*x'z)));
retp(cov);
endp;
/* Procedure for the function required by the Nelder-Meade algorithm */
proc ofn(&logliml,start);
local logliml:proc;
retp(logliml(start));
endp;
具体你可以下载Econometric Foundations这本书的gauss code，里面有3sls估计
扫码加好友，拉您进群
全部版块
我的主页
收藏
