/* TWO-STEP AND MAXIMUM LIKELIHOOD IN THE PRESENCE OF HETEROSCEDASTICITY */

/* Load Data and Define Variables */
load path = c:\gauss8.0\classes\econ5350\data\;
load x[364,30] = tescores.txt;
tescore = x[.,2];
modern = x[.,3];
experience = x[.,4];
gender = x[.,5];
age = x[.,6];
elemed = x[.,7];
seced = x[.,8];
somecollege = x[.,9];
college = x[.,10];
prof = x[.,11];
plots = x[.,12];
crops = x[.,13];
erosion = x[.,14];
fertility = x[.,15];
aptitude = x[.,16];
slope = x[.,17];
pests = x[.,18];
weeddens = x[.,19];
weedhght = x[.,20];
plantdis = x[.,21];
hydromorph = x[.,22];
lowland = x[.,23];
irrigated = x[.,24];
raindays = x[.,25];
rainqty = x[.,26];
region2 = x[.,27];
region3 = x[.,28];
yr94 = x[.,29];
yr95 = x[.,30];
nobs = rows(x);
constant = ones(nobs,1);

y = -ln(tescore);
xmat = constant~modern~modern^2~experience~experience^2~gender~age~age^2~elemed~
        seced~somecollege~college~prof~plots~plots^2~crops~crops^2~erosion~fertility~
        aptitude~slope~slope^2~pests~pests^2~weeddens~weeddens^2~weedhght~weedhght^2~
        plantdis~plantdis^2~hydromorph~lowland~irrigated~raindays~raindays^2~
        rainqty~rainqty^2~region2~region3~yr94~yr95;
k = cols(xmat);
vnum = seqa(1,1,k);

/* *** */
/* OLS */
/* *** */

b = inv(xmat'*xmat)*(xmat'*y);
resid = y - xmat*b;
s2 = resid'*resid/(nobs-k);

/* ******************* */
/* TWO-STEP ESTIMATION */
/* ******************* */

lne2 = ln(resid.*resid);
z = constant~region2~region3;
a = inv(z'*z)*(z'*lne2);
sigma2hat = exp(z*a);
omega = diagrv(eye(nobs),sigma2hat);
invomega = inv(omega);
betahat = inv(xmat'*invomega*xmat)*(xmat'*invomega*y);
residstar = y - xmat*betahat;
s2star = residstar'*invomega*residstar/(nobs-k);
varbetahat = s2star*inv(xmat'*invomega*xmat); 
sebetahat = sqrt(diag(varbetahat));
stdnstat = betahat./sebetahat;
print "Variable    Two-Step Estimate    Two-Step Standard Error     N(0,1) Statistic ";
print vnum~betahat~sebetahat~stdnstat;
print;

/* ***************************** */
/* MAXIMUM LIKELIHOOD ESTIMATION */
/* ***************************** */

library cml;
#include cml.ext;

/* Likelihood Procedure */
data = y~xmat;
proc lnlk(theta,data);
	local f, dev, logl, fmatinv, sigma2;
	sigma2 = theta[k+1];
	f = exp(z*theta[k+2:k+4]);
	dev = data[.,1] - data[.,2:k+1]*theta[1:k];
	logl = -0.5*(ln(2*pi) + ln(sigma2) + ln(f) + (1/sigma2)*(dev.*dev)./f);
	retp(logl);
endp;

/* Supplying an Analytical Gradient */
proc gradient(theta,data);
	local grad, dev, sigma2, f, g, g1, g2, g3;
	dev = data[.,1] - data[.,2:k+1]*theta[1:k];
	sigma2 = theta[k+1];
	f = exp(z*theta[k+2:k+4]);
	g = z.*f;
	g1 = (1/sigma2)*((data[.,2:k+1].*dev)./f);
	g2 = -0.5*((1/sigma2) - (1/(sigma2*sigma2))*(dev.*dev)./f); 
	g3 = -0.5*((g./f) - (1/sigma2)*(dev.*dev.*g./(f^2)));
	grad = g1~g2~g3;
	retp(grad);			
endp; 

/* Call the CML Module */
cmlset;
_cml_DirTol = 1e-3;
_cml_Algorithm = 3;
_cml_GradProc = &gradient; 
_cml_GradCheckTol = 0.01;
startv = betahat|s2star|a;
{mltheta,f0,g,cov,ret} = CMLPrt(CML(data,0,&lnlk,startv));