/* 

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

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

*/

new; 

format /m1 /rd 7,2;     @sets format for output printed to screen or file@

library pgraph;

@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); 
dy=y[2:rows(y)]-y[1:rows(y)-1];  
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@    




@Parametric Bootstraps to get critical value and p-value for test statistic
 and confidence interval based on percentile bootstrap@

/*
Steps:
1. Find estimate and test statistic for actual data and store

2. Find parameters for boostrap distribution #1 (i.e., estimates imposing null)

3. Repeat the following for j=1,...,B
    a. Simulate samples from BDGP #1
    b. Find estimate for bootstrap sample under H1
    c. Find test statistic for bootstrap sample under H0
    d. Store test stat, set j=j+1, return to a until j>B

4. Based on bootstrap distribution of test statistic, calculate critical value(s) and p-value
for actual test stat

5. Repeat the following for j=1,...,B
    a. Simulate samples from BDGP #2 (under alternative)
    b. Find estimate for bootstrap sample under H1
    c. Store estimate, set j=j+1, return to a until j>B

6. Based on bootstrap distribution of estimate, construct percentile bootstrap confidence interval
*/

@Globals@
alpha=0.05; @level at which to set Type I error@
prm_num=1;  @3=ar(1) coefficient@
b=99; @number of bootstrap samples@

rndseed 444444;  @sets a seed for all generated random numbers. Thus, you will get the same results each time you run th
                    program.@




@Step 1@
@MLE for H1 Model on Actual Data@
@if constructing LR, you will need to do MLE for both H0 and H1@

@initial values for optimization@
prmtr_in=zeros(3,1); 
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@

{xout,fout,gout,cout}=qnewton(&lik_fcn,prmtr_in); 

prm_fnl=trans(xout); @final constrained point estimates@
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@

b_hat=prm_fnl[3]; @change which parameter for different tests. Here we will consider AR(1) parameter@
t_hat=abs(prm_fnl[3]/se_fnl[3]); @change for tests of different parameters. Also, this makes the t-stat of interest a t-test. You will also consider an LR test@

@Store data, estimate, and test stat for actual data@
dy_=dy;
b_hat_=b_hat;   
t_hat_=t_hat;
prm_fnl_=prm_fnl;
se_fnl_=se_fnl;

mu_h1=prm_fnl[1];
var_h1=prm_fnl[2];
phi_h1=prm_fnl[3];




@Step 2@
@Estimate model under H0 to get bootstrap DGP#1@
@You do not have to repeat this for LR, given that you have already estimated DGP under H0@

prmtr_in=zeros(2,1); 
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@

{xout,fout,gout,cout}=qnewton(&lik_fcn_h0,prmtr_in); 

prm_fnl=trans_h0(xout); @final constrained point estimates@

mu_h0=prm_fnl[1];
var_h0=prm_fnl[2];



screen off; @stops GAUSS from writing all of the qnewton steps to the screen; reversed below@


@Step 3 Loop over Bootstrap experiments for test stat@
out_mat1={};

j=1;
do until j>b;

@3a. Simulate from Bootstrap DGP H0@
@BDGP#1: H0@
dy=mu_h0+sqrt(var_h0)*rndn(t,1);


@3b. and 3c. Find estimate and test stat for bootstrap DGP@ 
@MLE for H1 Model on Bootstrap Data@
@if constructing LR, you will need to do MLE for both H0 and H1@

@initial values for optimization@
prmtr_in=zeros(3,1); 
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@

{xout,fout,gout,cout}=qnewton(&lik_fcn,prmtr_in); 

prm_fnl=trans(xout); @final constrained point estimates@
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@

b_hat=prm_fnl[prm_num]; 
t_hat=abs(prm_fnl[prm_num]/se_fnl[prm_num]); 

@3d Store and update loop@
out_mat1=out_mat1|t_hat;

j=j+1;
endo;




@Step 4 Loop over Bootstrap experiments for estimate to construct confidence interval@

out_mat2={};

j=1;
do until j>b;

@4a. Simulate from Bootstrap DGP H1@

@BDGP#2: H1@
dy=recserar(sqrt(var_h1)*rndn(t+100,1),0,phi_h1);   @built in GAUSS function for constructing AR series@
dy=dy[101:rows(dy)]; @toss out first 100 values to eliminate effect of initial value, which is set to 0@
dy=dy+mu_h0; @add unconditional mean to series@



@4b. and 4c. Find estimate for bootstrap DGP@ 
@MLE for H1 Model on Bootstrap Data@

@initial values for optimization@
prmtr_in=zeros(3,1); 
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@

{xout,fout,gout,cout}=qnewton(&lik_fcn,prmtr_in); 

prm_fnl=trans(xout); @final constrained point estimates@
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@

b_hat=prm_fnl[prm_num]; 
t_hat=abs(prm_fnl[prm_num]/se_fnl[prm_num]); 

@3d Store and update loop@
out_mat2=out_mat2|b_hat;

j=j+1;
endo;

screen on;





@Step 4: Find critical value and p-value based on bootstrap distribution under H0@
t_boot=abs(out_mat1); @by taking abs, we assume distribution is symmetric for purposes of calculating critical value and p-value@
t_boot=sortc(t_boot,1); @sorts the test statistics from smallest to largest@

crit=t_boot[ceil((1-alpha)*rows(t_boot)),1];
pv_mat=(t_hat_~1)|(t_boot~zeros(rows(t_boot),1));
pv_mat=sortc(pv_mat,1);
pv=(rows(pv_mat)-indexcat(pv_mat[.,2],1))/b;




@Step 5: Find confidence interval based on percentile approach based on bootstrap under H1@
@Find boostrap distribution of b_hat@
b_boot=out_mat2;
b_boot=sortc(b_boot,1);

q_l=b_boot[ceil((alpha/2)*rows(b_boot))];
q_u=b_boot[ceil((1-alpha/2)*rows(b_boot))];

ci_l=2*b_hat_-q_u;    @remember that the quantiles should be flipped in case the bootstrap distribution is asymmetric@
ci_u=2*b_hat_-q_l;    @again quantile is flipped@


cls;

@Output@
output file=gauss_example4.out reset;
output on;
"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(b=0)|=";;abs(prm_fnl_[prm_num]/se_fnl_[prm_num]);;"  (given H0, which may be on another parameter)";
100*(1-alpha);;"% bootstrap critical value (two-tailed)=";;crit;
"bootstrap p-value=";;pv;
"";
100*(1-alpha);;"% confidence interval based on percentile bootstrap";
"[";;ci_l;;",";;ci_u;;"]";
"";
output off;

end;  

@========================================================================@
@========================================================================@
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]@0; @ar coefficient@
	theta1=@prmtr[5]@0; @ma coefficient@
	beta1=@prmtr[6]@0; @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;

@========================================================================@
proc LIK_FCN_h0(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_h0(prmtr1); @constrain parameters@

@PARAMETERS@
	mu=prmtr[1]; @unconditional mean@
	var=prmtr[2]; @standard deviation of error term@
	phi1=@prmtr[3]@0; @ar coefficient@
	phi2=@prmtr[4]@0; @ar coefficient@
	theta1=@prmtr[5]@0; @ma coefficient@
	beta1=@prmtr[6]@0; @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_h0(c0); 

local c1,aaa,ccc;

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


retp(c1);
endp;