/* PROGRAM FINDS EXACT MLE FOR ARMA MODELS OF
GDP GROWTH 

PROGRAM DESIGNED BY JAMES MORLEY 2/9/00
Modified 9/5/08 */

@INITIAL GAUSS STUFF@
new; @clear memory@
library  pgraph; @call up GAUSS modules for optimization and graphing@
format /m1 /rd 9,6; @set details of output format@

@LOAD DATA@
load data[215,1]=gdpq.prn; @1954:1 -2007:3@
y=ln(data);
ty=rows(y);

dy=(y[2:ty]-y[1:(ty-1)]);  
T=rows(dy); @sample size@

@MAXIMUM LIKELIHOOD ESTIMATION@
prmtr_in={0 0 1 -1 -1.5 0.5}; @starting values for parameters; try different ones to get something close to CMLE results in EViews@
prmtr_in[1]=meanc(dy); @use sample average as starting value for unconditional mean@

@!!!make sure you have the same number of starting values as parameters!!!@

prmtr_in=prmtr_in';
{xout,fout,gout,cout}=qnewton(&lik_fcn,prmtr_in); @calls optmum to calculate MLEs@

prm_fnl=trans(xout); @final constrained point estimates@
lik=-fout; @final likelihood value@
cov0=inv(hessp(&lik_fcn,xout)); @calculate hessian@
grdn_fnl=gradp(&TRANS,xout);
cov=grdn_fnl*cov0*grdn_fnl'; @delta method @ 
sd_fnl =sqrt(diag(cov)); @Standard Errors of Parameter Estimates@

output file=hw3.out reset;
output on;
"==FINAL OUTPUT========================================================";
"likelihood value is ";; lik;
"Estimated parameters are:";
prm_fnl';
"(";;sd_fnl';;")";
"===============================================================";
output off;

END;

@========================================================================@
@========================================================================@
PROC LIK_FCN(PRMTR1);

local prmtr,mu,sig,phi1,phi2,theta1,theta2,mu_tilda,f,q,h,beta_ll,p_ll,
	llikv,j,beta_tl,p_tl,eta_tl,f_tl,beta_tt,p_tt,A,R,z;

@PARAMETERS@
     prmtr=trans(prmtr1); @constrain parameters@
	mu=prmtr[1]; @unconditional mean@
	sig=prmtr[2]; @standard deviation of error term@
	phi1=prmtr[3]; @ar coefficient@
	phi2=prmtr[4]; @ar coefficient@
	theta1=prmtr[5]; @ma coefficient@
	theta2=prmtr[6]; @ma coefficient@


@SET UP STATE SPACE FORM@

@Observation Equation@

H= 1~theta1~theta2;

A= mu;
z=ones(t,1);

R=0;

@State Equation@

MU_TILDA= 0|0|0;

F=	phi1~phi2~0|
	1~0~0|
	0~1~0;

Q= sig^2~0~0|
	0~0~0|
	0~0~0;



@INITIAL STATE ESTIMATES@
BETA_LL=	0|
		 	0|
		 	0; @unconditional expectation of state vector@
P_LL = inv( (eye(rows(F)^2)-F.*.F) )*vec(Q); @see Kim & Nelson, ch. 3@
P_LL=  reshape(P_LL,rows(F),rows(F));

llikv = 0;
j=1;
do until j>T;

@KALMAN FILTER@
BETA_TL = MU_TILDA + F * BETA_LL;
P_TL = F* P_LL * F' + Q;

ETA_TL = dy[j,.]' - H * BETA_TL - A*z[j,.]';
F_TL = H * P_TL * H' + R;

BETA_TT = BETA_TL + P_TL*H'*INV(F_TL) * ETA_TL;
P_TT = P_TL - P_TL * H'*INV(F_TL) * H * P_TL;


@LOG LIKELIHOOD FUNCTION@
llikv=llikv - 0.5*ln(2*PI*F_TL) - 0.5*(ETA_TL^2)/F_TL;

@UPDATE FOR NEXT ITERATION@
BETA_LL=BETA_TT;  
P_LL=P_TT;

j=j+1;
endo;

retp(-llikv);
endp;
@========================================================================@

/*-------------------*/
@PROCEDURE TO CONSTRAIN PARAMETERS@
proc trans(c0); 

local c1,aaa,ccc;

      c1 = c0;
	c1[2]=exp(-c0[2]); @positive variance@


@USE THE FOLLOWING CONSTRAINTS for ARMA(2,2)@

        aaa=c0[3]./(1+abs(c0[3]));
     	ccc=(1-abs(aaa))*c0[4]./(1+abs(c0[4]))+abs(aaa)-aaa^2;

     c1[3]=2*aaa;
     c1[4]= -1* (aaa^2+ccc) ; 

c1[6]=c0[6]/(1+abs(c0[6]));  


@USE THE FOLLOWING CONSTRAINTS for ARMA(2,1)@
/*
        aaa=c0[3]./(1+abs(c0[3]));
     	ccc=(1-abs(aaa))*c0[4]./(1+abs(c0[4]))+abs(aaa)-aaa^2;

     c1[3]=2*aaa;
     c1[4]= -1* (aaa^2+ccc) ; 

c1[5]=c0[5]/(1+abs(c0[5]));  
*/


@USE THE FOLLOWING CONSTRAINTS for ARMA(1,2)@
/*
	c1[3]=c0[3]/(1+abs(c0[3])); 
	c1[5]=c0[5]/(1+abs(c0[5])); 
*/

@USE THE FOLLOWING CONSTRAINTS for ARMA(1,1)@
/*
	c1[3]=c0[3]/(1+abs(c0[3])); 
	c1[4]=c0[4]/(1+abs(c0[4])); 
*/

@USE THE FOLLOWING CONSTRAINTS for AR(1) or MA(1)@
/*
	c1[3]=c0[3]/(1+abs(c0[3])); 
*/


retp(c1);
endp;
/*-------------------*/