/* 

The following program represents the third in a series of primers on basic commands in GAUSS.
Written by James Morley, Feb. 10th, 2008

In this case, we will use the Kalman filter to construct the prediction error decomposition
of the likelihood funciton for ARMA models. We will use a numerical optimization routine
to maximize the likelihood with respect to the ARMA parameters.

This code also calculates test statistics, critical values, p-values, and confidence intervals
based on asymptotic distributions 

*/

new; 


@LOAD DATA@
load data[243,1]=gdp.txt; @put your data file of raw US real GDP here; make sure dimension is correct@

y=100*ln(data); @Take natural logs and multiply by 100@

dy=y[2:rows(y)]-y[1:rows(y)-1];  @Take differences of 100 times logs to get continuously compounded growth rates@
t=rows(dy); @sample size@

load data[243,1]=dhstar.txt; @put your data file of inventories here. In terms of the file, I've already transformed it as suggested in question 1 part i@

z=data[1:t];    @sets z to lag of dhstar@    


@MAXIMUM LIKELIHOOD ESTIMATION@
prmtr_in=zeros(6,1); @starting values for parameters; make sure # of parameters is correct@
prmtr_in[1]=meanc(dy);  @Use the sample mean as starting value for mu@
prmtr_in[2]=ln((stdc(dy))^2);   @Use sample variance of dy as starting value for var_e@
@!!!make sure you have the same number of starting values as parameters when you alter the model you want to estimate!!!@

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

/*
qnewton is a built in numerical optimization routine. It is a minimization procedure.
"xout"=parameter values at minimum
"fout"=value of function at minimum
"lik_fcn"=user defined procedure that evaluates function given parameter values
"&lik_fcn"=qnewton calls lik_fcn procedure multiple times for different parameter values
"prmtr_in"=initial guess for the parameters

Note: prmtr_in is in terms of unconstrained parameter values that can range between -infinity and +infinity
We constrain these parameters in the lik_fcn to correspond to the model parameters we are interested in.
For example, consider the variance parameter. The associated unconstrained parameter varies from
-infinity to +infinity. However, this unconstrained parameter gets constrained by the trans function
as follow: constrained = exp(unconstrained). Thus, the constrainted parameter will always be 
non-negative. Then, when the likelihood is constructed the variance is set to the constrained value
whenever it shows up in the calculation of the likelihood. 

*/

prm_fnl=trans(xout); @final constrained point estimates@
lik=-fout; @final likelihood value@
cov0=inv(hessp(&lik_fcn,xout)); @calculate inverse negative hessian; note that qnewton returns negative hessian since it conducts minimization@
grdn_fnl=gradp(&TRANS,xout);
cov=grdn_fnl*cov0*grdn_fnl'; @delta method to get var-cov for constrained parameters@ 
se_fnl =sqrt(diag(cov)); @Standard Errors of Parameter Estimates@


alpha=0.05; @level at which to set Type I error@
prm_num=3; @parmater to be tested@

output file=gauss_example3.out reset;
output on;
"==FINAL OUTPUT========================================================";
"";
"likelihood value is ";; lik;
"Estimated parameters are:";
prm_fnl';
"(";;se_fnl';;")";
"standard errors in parentheses";
"";
"";
"Estimate, SE, and T-stat for parameter #";;prm_num;;"(e.g., 3=AR(1) coefficient):";
"";
"b_hat=";;prm_fnl[prm_num];
"se(b_hat)=";;se_fnl[prm_num];
"";
"|t(phi1=0)|=";;abs(prm_fnl[prm_num]/se_fnl[prm_num]);
100*(1-alpha);;"% critical value (two-tailed)=";;cdfni(1-(alpha/2));
"p-value=";;2*(1-cdfn(abs(prm_fnl[prm_num]/se_fnl[prm_num])));
"";
100*(1-alpha);;"% confidence interval based on inverting t-stat";
"[";;prm_fnl[prm_num]-cdfni(1-(alpha/2))*se_fnl[prm_num];;",";;prm_fnl[prm_num]+cdfni(1-(alpha/2))*se_fnl[prm_num];;"]";
"";
"===============================================================";
output off;

end;  @remove this end in order to @

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

local prmtr,mu,var,phi1,phi2,theta1,f,q,h,x_ll,p_ll,
	llikv,j,x_tl,p_tl,eta_tl,f_tl,x_tt,p_tt,vP_ll,A,beta1;

     prmtr=trans(prmtr1); @constrain parameters@

@PARAMETERS@
	mu=prmtr[1]; @unconditional mean@
	var=prmtr[2]; @standard deviation of error term@
	phi1=prmtr[3]; @ar coefficient@
	phi2=prmtr[4]; @ar coefficient@
	theta1=prmtr[5]; @ma coefficient@
	beta1=prmtr[6]; @coefficient on Z variable@


@SET UP STATE SPACE FORM@

F=	phi1~phi2|
	1~0;

Q= var~0|
	0~0;

H= 1~theta1;

A= beta1;


@INITIAL STATE ESTIMATES@
X_LL=	zeros(rows(F),1); @unconditional expectation of state vector@
vP_LL = inv( (eye(rows(F)^2)-F.*.F) )*vec(Q); 
P_LL=  reshape(vP_LL,rows(F),rows(F));

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

@KALMAN FILTER@
X_TL = F * X_LL;
P_TL = F* P_LL * F' + Q;

ETA_TL = dy[j] - mu - H * X_TL - A*z[j];
F_TL = H * P_TL * H';

X_TT = X_TL + P_TL*H'*INV(F_TL) * ETA_TL;
P_TT = P_TL - P_TL * H'*INV(F_TL) * H * P_TL;


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

@UPDATE FOR NEXT ITERATION@
X_LL=X_TT;  
P_LL=P_TT;

j=j+1;
endo;

retp(-llikv); @return negative of likelihood because qnewton is a minimization routine@
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(p,1)@

c1[5]=c0[5]/(1+abs(c0[5]));  @comment out if estimating AR model only@



@USE THE FOLLOWING CONSTRAINTS FOR ARMA(1,q)@

	@c1[3]=c0[3]/(1+abs(c0[3]));@  @comment out if estimating AR(2) or ARMA(2,q) models; remove comments from below@


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


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


retp(c1);
endp;