/* ESTIMATION WITH AUTOCORRELATED ERRORS -- APPLICATION TO THE CONSUMPTION FUNCTION */

/* LOAD DATA AND DEFINE VARIABLES */
load path = c:\gauss35\classes\econ5340\data\;
load x[37,3] = table6.10;
t = x[2:37,1];
y = x[2:37,2];
c = x[2:37,3];
nobs = rows(t);
constant = ones(nobs,1);
xmat = constant~t~y;
k = cols(xmat);

/* OLS ESTIMATION */
bols = inv(xmat'*xmat)*(xmat'*c);
b = bols;
print "     OLS Intercept      OLS Trend        OLS Slope";
print b';
print;

/* ******************** */
/* EFFICIENT ESTIMATION */
/* ******************** */

/* COCHRANE ORCUTT ESTIMATION */
ct1 = c[1:nobs-1,1];
ct = c[2:nobs,1];
xmatt1 = xmat[1:nobs-1,.];
xmatt = xmat[2:nobs,.];
cc = 1;
print "     CO Intercept        CO Trend        CO Slope         CO Rho";

do while cc>1e-3;
	resids = c - xmat*b;
	err1 = resids[1:nobs-1,1];
	err = resids[2:nobs,1];
	rhohat = inv(err1'*err1)*(err1'*err);
	cstar = ct - rhohat*ct1;
	xstar = xmatt - rhohat*xmatt1;
	bold = b;
	b = inv(xstar'*xstar)*(xstar'*cstar); 
	print b[1,1]~b[2,1]~b[3,1]~rhohat;
	cc = (bold - b)'*(bold - b);
endo;

/* PRAIS-WINSTEN ESTIMATION */
b = inv(xmat'*xmat)*(xmat'*c);
cc = 1;
print;
print "     PW Intercept        PW Trend        PW Slope         PW Rho";

do while cc>1e-3;
	resids = c - xmat*b;
	err1 = resids[1:nobs-1,1];
	err = resids[2:nobs,1];
	rhohat = (err'*err1)/(resids'*resids);
    cstar = (sqrt(1-rhohat^2)*c[1,1])|(ct - rhohat*ct1);
	xstar = (sqrt(1-rhohat^2)*xmat[1,.])|(xmatt - rhohat*xmatt1);
	bold = b;
	b = inv(xstar'*xstar)*(xstar'*cstar); 
	print b[1,1]~b[2,1]~b[3,1]~rhohat;
	cc = (bold - b)'*(bold - b);
endo;

/* MAXIMUM LIKELIHOOD ESTIMATION */
library cml;
#include cml.ext;

/* Likelihood Procedure */
data = c~xmat;
proc lnlk(theta,data);
	local ystar, ystar1, xstar, xstar1, epstar, beta, sigma2, rho, logl;
	beta = theta[1:3];
	sigma2 = theta[4];
	rho = theta[5];
	ystar1 = sqrt(1-rho^2)*data[1,1];
	ystar = ystar1|(data[2:nobs,1] - rho*data[1:nobs-1,1]);
	xstar1 = sqrt(1-rho^2)*data[1,2:4];
	xstar = xstar1|(data[2:nobs,2:4] - rho*data[1:nobs-1,2:4]);
	epstar = ystar - xstar*beta;
	logl = -0.5*(nobs*(ln(2*pi) + ln(sigma2)) + (1/sigma2)*(epstar'*epstar)) + 0.5*ln(1-rho^2);
	retp(logl);
endp;

/* Supplying an Analytical Gradient */
proc gradient(theta,data);
	local grad, g1, g2, g3, beta, sigma2, rho, ystar, ystar1, xstar, xstar1, eps, epstar;
	beta = theta[1:3];
	sigma2 = theta[4];
	rho = theta[5];
	ystar1 = sqrt(1-rho^2)*data[1,1];
	ystar = ystar1|(data[2:nobs,1] - rho*data[1:nobs-1,1]);
	xstar1 = sqrt(1-rho^2)*data[1,2:4];
	xstar = xstar1|(data[2:nobs,2:4] - rho*data[1:nobs-1,2:4]);
	eps = data[.,1] - data[.,2:4]*beta;
	epstar = ystar - xstar*beta;
	g1 = (1/sigma2)*((epstar'*xstar)); 
   	g2 = -(0.5*nobs/sigma2) + (0.5/sigma2^2)*((epstar'*epstar));
   	g3 = (1/sigma2)*(epstar[2:nobs]'*eps[1:nobs-1]) - (rho/(1-rho^2)) + 2*rho*eps[1,1]/sigma2;
   	grad = g1~g2~g3;
	retp(grad);
endp;

/* Call the CML Module */
cmlset;
_cml_GradProc = &gradient;
_cml_GradCheckTol = 0.01; 
resids = cstar - xstar*b;
s2 = (resids'*resids)/(nobs-k);
startval = bols|s2|rhohat;
{mltheta,f,g,cov,ret} = CMLPrt(CML(data,0,&lnlk,startval));