new; /* erases everything in memory and closes files and windows 1 */

/*	This is an annotated Gauss batch program for estimating 
**      simultaneous equations by two systems methods: FIML  and 3SLS.
**      It is based in large part on Example 17.3.1 in Mittelhammer 
**      et al. (2000).
**  The data are generated from a three-equation system of the general form:
**
**  			       Y*Gamma + X*B + E = 0
**
**  Notice that the signs of B and E are the reverse of those in 
**  Mittelhammer et al. eq. (17.1.1).
**  The program computes 3SLS estimates and then uses these as starting
**  points for iterative FIML estimation. To obtain FIML estimates, the program
**	uses the concentrated log-likelihood function in Equations (17.3.7)
**	and (17.3.8) to solve for the parameters of a structural model and to
**	derive an estimate of the contemporaneous covariance matrix.  
**	DSP stated in Equation (17.3.3) and then computes the FIML parameter
**	estimates and the associated estimated standard errors.  The GAUSS
**	program also reports the 3SLS estimates for comparison purposes.
**	Finally, the program conducts a hypothesis test and reports a
**	confidence region for a subset of the model parameters. 
*/

/*	Set some graphics and program parameters  */
gausset; /* Resets global control variables declared in gauss.dec */ 
/* The next 3 lines set up an output file */
format /rd 11,5; /* right justified numbers of form xxxxxx.xxxxx */
outwidth 110; /* 110 columns in output */
output file = FIML-3SLS.out reset; /* Overwrite previous version of file */

print on;

begin:
/* Choose number of observations in simulated data set */
n = 100;

/*	Specify the a noise covariance matrix for the simulated data set */
let Sigma[3,3] =	 1.000    -1.000    -0.125
			-1.000     4.000     0.0625
			-0.125     0.0625    2.5000;

/*	Check to see that specified  sigma matrix is positive definite */
/*      The eigenvalues should all be positive. */

Eigenvalues = eigrs(Sigma);


/*  Generate the sample data for the three-equation system */
/*	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];

print off;
/*	Generate the observations  */
again:
x = ones(n,1)~sqrt(10)*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;

/* Compute 3SLS and then FIML estimates of the parameters of the 
** three-equation system along with an estimate of the asymptotic covariance
** matrix of the estimator. 
** The algorithm used to solve for the FIML estimates is the non-
** gradient based Nelder-Meade algorithm (a proc near the end of this program). 
*/

/*	First, compute the 3SLS estimates  */
y = vec(ymat);

/*	Compute the reduced-form estimates by LS  */
print on;
b1 = olsqr(ymat[.,1], xmat);
b2 = olsqr(ymat[.,2], xmat);
b3 = olsqr(ymat[.,3], xmat);
print off;

/*	Compute the predicted values  */
yh1 = xmat*b1;
yh2 = xmat*b2;
yh3 = xmat*b3;

z1 = yh2~x1;
z2 = yh1~yh3~x2;
z3 = yh1~x3;

print on;
/*	Compute the 2SLS estimates  */
b2s1 = olsqr(ymat[.,1], z1);
b2s2 = olsqr(ymat[.,2], z2);
b2s3 = olsqr(ymat[.,3], z3);
print off;

/*	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 3SLS estimates 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);

block1y = temp[1,1]*z1'x*inv(x'x)*x'ymat[.,1] + temp[1,2]*z1'x*inv(x'x)*x'ymat[.,2] + temp[1,3]*z1'x*inv(x'x)*x'ymat[.,3];
block2y = temp[2,1]*z2'x*inv(x'x)*x'ymat[.,1] + temp[2,2]*z2'x*inv(x'x)*x'ymat[.,2] + temp[2,3]*z2'x*inv(x'x)*x'ymat[.,3];
block3y = temp[3,1]*z3'x*inv(x'x)*x'ymat[.,1] + temp[3,2]*z3'x*inv(x'x)*x'ymat[.,2] + temp[3,3]*z3'x*inv(x'x)*x'ymat[.,3];
blocky = block1y|block2y|}block3y;

b3sls = inv(blockx)*blocky;
cov3sls = inv(blockx);
s3sls = sqrt(diag(cov3sls));

/*	Report the 3SLS estimates to the user  */
format /mat /rd 16,5;
print "3SLS Parameter Estimates";
print;
print "                          3SLS            Standard         Asymptotic";
print "      True value         Estimate           Error          Z-statistic";
btrue~b3sls~s3sls~b3sls./s3sls;
print;

/*  Computational Note:  It is possible, as in all nonlinear estimation attempts,
**  that the maximum likelihood algorithm will fail for a number of reasons.  For
**  example, the algorithm may not converge in the number of iterations allowed,
**  or a matrix may be singular at non-optimal parameter values (i.e., during
**  the search process).  In these cases, rerun the example with a new data set
**  or different starting values.
**  If the program stops with an error message, simply close the GAUSS
**  window and rerun the example.
*/

/*	Initiate the starting values  */
start = b3sls;

zmat = ymat~xmat;

/*	Run the Nelder-Meade minimization algorithm to try to find the FIML estimates.
**	The procedure ofn() links the objective function
**	to the function that the Nelder-Meade algorithm attempts to minimize.  
*/

neldtol = 1e-7;
bfiml = nelder(start,&logfiml,20000,0,neldtol,0); 

covfiml = asycov(bfiml);
sfiml = sqrt(diag(covfiml));

/*	Report the FIML estimates  */
format /mat /rd 16,5;
print "FIML Parameter Estimates";
print;
print "                           FIML            Standard          Asymptotic";
print "      True value         Estimate           Error            Z-statistic";
btrue~bfiml~sfiml~bfiml./sfiml;
print;

print "Comparison of Coefficient Estimates";
print;
print "      True value           3SLS            FIML";
print btrue~b3sls~bfiml; 
format /mat /rd 5,0;
print;
print "Sample size, n = " n;

format /mat /rd 16,5;
print;
print "Comparison of Standard Error Estimates";
print;
print "         3SLS              FIML";	
print s3sls~sfiml; 
format /mat /rd 5,0;
print;
print "Sample size, n = " n;
print;

r = zeros(1,15);
r[1,13] = 1;
r[1,14] = 1;
r[1,15] = 1;

t1 = (r*b3sls - 1)/sqrt(r*cov3sls*r');
t1pvalue = cdfnc(t1);
t2 = (r*bfiml - 1)/sqrt(r*covfiml*r');
t2pvalue = cdfnc(t2);

format /mat /rd 7,4;

/* Test the one-sided null hypothesis for the third equation,
** H0 : B[3,3] + B[4,3] + B[6,3] <= 1, which is not true (the parameters sum to 1.2).
** Using the large sample results for the estimators, we compute the one-sided
** Z-test statistics and the associated p-values for this hypothesis:
*/
print;
print "                       Test statistic      P-value";
print;
print "           3SLS          " t1 "          " t1pvalue;
print "           FIML          " t2 "          " t2pvalue;
print;

end; 

/*	Procedure for the concentrated log-likelihood function  */
proc logfiml(b);
	local q,LL,sighat, resid, invsig, beta, gamm, cnst;
     
	gamm = -1|b[1]|0;
	gamm = gamm~(b[6]|-1|b[7]);
	gamm=gamm~(b[11]|0|-1);
	q = cols(gamm);

	beta = b[2]|0|b[3]|0|b[4]|0|b[5];
	beta = beta~(b[8:9]|0|0|b[10]|0|0);
	beta = beta~(b[12]|0|b[13:14]|0|b[15]|0);

	resid = -ymat*gamm - xmat*beta;

	sighat = resid'resid/n;

	resid = resid * chol(inv(sighat))';
  
	LL = - 0.5*ln(det(sighat)) + ln(abs(det(gamm)));

	retp(-LL);
endp;



/*	Procedure for the asymptotic covariance matrix  */
proc asycov(b);
	local z1,z2,z3,q,cov,sighat, resid, invsig, beta, gamm, yhat;
	local block11,block12,block13,block21,block22,block23;
	local block31,block32,block33,blockx;

	gamm = -1|b[1]|0;
	gamm = gamm~(b[6]|-1|b[7]);
	gamm = gamm~(b[11]|0|-1);
	q = cols(gamm);

	beta = b[2]|0|b[3]|0|b[4]|0|b[5];
	beta = beta~(b[8:9]|0|0|b[10]|0|0);
	beta = beta~(b[12]|0|b[13:14]|0|b[15]|0);

	resid = -ymat*gamm - xmat*beta;
	sighat = resid'resid/n;

	yhat = ymat - resid;
	z1 = yhat[.,2]~x1;
	z2 = yhat[.,1]~yhat[.,3]~x2;
	z3 = yhat[.,1]~x3;        

	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);

	cov = pinv(blockx);

	retp(cov);
endp;




/*	Procedure required by the Nelder-Meade algorithm  */
proc ofn(&logfiml,start);
	local logfiml:proc;

	retp(logfiml(start));
endp;


proc nelder(theta,&fn,itlimit,neldout,ftol,scaleit);
	local changes,alpha,beta,gam,npar,hi,par,ii,y,yp,iter,rtol,
		fn:proc,x0,indxmax,indxmin,indxmax2,alphap1,
		onemgam,onembeta,mxr,mxe,mxc,xr,xe,xc,maxval,minval,maxval2,e;

	if neldout == 0;
		screen off;
	endif;
	if scaleit == 0;
		SCALEIT = 0.1;
	endif;
	if ftol == 0;
		ftol=1e-14;
	endif;

	changes = theta.*SCALEIT + (theta.==0).*SCALEIT;
	alpha = 1.0;
	beta = 0.5;
	gam = 2.0;

	/*	Calculate the objective function at (npar + 1) initial values of theta.
	**	Store the ofn() results in y, which is a ((npar+1) X 1) column vector.
	**	Place the values of theta as columns of par, a (npar X (npar+1)) matrix.  */
	npar = rows(theta);
	hi = npar + 1;
	par = theta~(theta+eye(npar).*changes);
	ii = 1;
	y = {};
	do while ii <= hi;
		y = y|ofn(&fn,par[.,ii]);
		ii = ii + 1;
	endo;

	/*	Sort the values of y in ascending order, and sort the corresponding
	**	columns of par  */
	yp = sortc(y~par',1);
	y = yp[.,1];
	par = yp[.,2:hi];

	x0 = meanc(par[1:npar,.]);
	indxmax = hi;
	indxmin = 1;
	indxmax2 = hi-1;
	alphap1 = alpha+1;
	onemgam = 1 - gam;
	onembeta = 1 - beta;

	iter = 1;
	rtol = 2*(abs(y[hi]-y[1]))/(abs(y[hi])+abs(y[1]));
	do while iter <= itlimit;
		xr = alphap1*x0 - (alpha*par[indxmax,.])';
		mxr = ofn(&fn,xr);
		if mxr < y[indxmin];
			xe = gam*xr + onemgam*x0;
			mxe = ofn(&fn,xe);
			if mxe < y[indxmin];
                		y[indxmax] = mxe;
				par[indxmax,.] = xe';
			else;
				y[indxmax] = mxr;
				par[indxmax,.] = xr';
			endif;
			goto testconv;
		elseif y[indxmax2] >= mxr;
			y[indxmax] = mxr;
			par[indxmax,.] = xr';
 			goto testconv;
		elseif mxr >= y[indxmax];
			xc = (beta*par[indxmax,.]')+ onembeta*x0;
			mxc = ofn(&fn,xc);
			if mxc >= y[indxmax];
				par = (.5)*(par+par[indxmin,.]);
				ii = 1;
				do while ii <= hi;
					y[ii] = ofn(&fn,par[ii,.]');    
					ii = ii + 1;
				endo;
				goto testconv;
			else;
				par[indxmax,.] = xc';
				y[indxmax] = mxc;
				goto testconv;
			endif;
		else;
			par[indxmax,.] = xr';
			y[indxmax] = mxr';
			xc = (beta*par[indxmax,.])'+ onembeta*x0;
			mxc = ofn(&fn,xc);
			if mxc >= y[indxmax];
				par = (.5)*(par+par[indxmin,.]);
				ii = 1;
				do while ii <= hi;
					y[ii] = ofn(&fn,par[ii,.]');     
					ii = ii + 1;
				endo;
				goto testconv;
			else;
				par[indxmax,.] = xc';
				y[indxmax] = mxc;
				goto testconv;
			endif;
		endif;
		testconv:
		rtol = 2*abs(y[indxmax]-y[indxmin])/abs(y[indxmax]+y[indxmin]);

		if rtol < ftol;
			goto finishup;
		else;
			maxval = maxc(y);
			minval = minc(y);
			indxmax = indexcat(y,maxval);
			indxmax = indxmax[1];             
			indxmin = indexcat(y,minval);
			indxmin = indxmin[1];             
			e = zeros(hi,1);
			e[indxmax[1]] = 1;
			maxval2 = maxc(delif(y,e));
			indxmax2 = indexcat(y,maxval2);
			indxmax2 = indxmax2[1];          
			x0 = meanc(delif(par,e));

			   " ";
			   "OBJECTIVE FUNCTION =";; y[indxmin];
			   locate 5,5;
			   "Nelder ITERATION =";; iter;        
		endif;  

	iter = iter+1; 
endo;

finishup:
if iter > itlimit;
	" ";
	" ";
	"FAILED TO CONVERGE";
	"PARAMETER VECTOR:";
	par[indxmin,.];
	"OBJECTIVE FUNCTION =";; y[indxmin];
	"RTOL=";; rtol;
else;
	" ";
	" ";
	"PROBABLE CONVERGENCE (TRY RESTARTING)";
	"PARAMETER VECTOR:";
	par[indxmin,.];
	"OBJECTIVE FUNCTION =";; y[indxmin];
	"RTOL =";; rtol;
	"ITERATIONS =";; iter-1;
endif;
screen on;

retp(par[indxmin,.]');

endp;