clear all;
close all;

%%%%%%%%%%%%%%%%%%%%%%%
%% DGP
%%%%%%%%%%%%%%%%%%%%%%%
%s = [3332253790];
%randn('state',s)
year =1951:1990;
N=length(year);
X = 10+5*randn(N,1);
Y = zeros(N,1);
Y(1)=5;
for j=2:N
    Y(j)=2+2*X(j)-0.5*X(j-1)+0.7*Y(j-1)+5*randn(1);
end
XX=[ones(N-1,1) X(2:N) X(1:(N-1)) Y(1:(N-1))];
YY=Y(2:N);

fprintf(1,'\n%s\n','===================================')
fprintf(1,'%s\n',  'A correctly specified relation')
fprintf(1,'%s\n',  '===================================');
reg1 = ols(YY,XX);
prt(reg1)


fprintf(1,'\n%s\n','Test of first-order autocorrelation');
fprintf(1,'%s\n','************************************************');
%%%%%%%%%%%%%%%%%%%%%%%%%
%% Breush_Godfrey Test
%%%%%%%%%%%%%%%%%%%%%%%%%
resid = reg1.resid(2:(N-1));
XXu = [XX(2:(N-1),:) reg1.resid(1:(N-2))];
reg_res=ols(resid,XXu);
disp('LM-test regression results'); prt(reg_res)
lmstat = reg_res.nobs*reg_res.rsqr;
lmprob = 1 - chi2cdf(lmstat,1);
fprintf(1,'\n%s\n','Breush-Godfrey Serial Correlation LM Test');
fprintf(1,'%s\n','***************************************************');
fprintf(1,'Obs*R-squared = %4.6f    ',lmstat);
fprintf(1,'Probability = %4.6f \n',lmprob);


%%%%%%%%%%%%%%%%%%%%%%%%%
%% Test of the Cochrane-Orcutt specification
%%%%%%%%%%%%%%%%%%%%%%%%%
fprintf(1,'\n%s\n','===================================')
fprintf(1,'%s\n',  'A misspecified relation')
fprintf(1,'%s\n',  '===================================');

X = [ones(N,1) X];
reg2 = ols(Y,X);
prt(reg2)
reg_co = olsc(Y,X);
prt(reg_co)

es = reg_co.resid;
e = reg1.resid;
waldstat = N*(es'*es-e'*e)/(e'*e);
waldprob = 1 - chi2cdf(waldstat,1);
fprintf(1,'\n%s\n','A test of the Cochrane-Orcutt specification');
fprintf(1,'%s\n','************************************************');
fprintf(1,'%s\n','Wald Test');
fprintf(1,'%s\n','Null Hypothesis: C(3)+C(2)*C(4)=0');
fprintf(1,'Chi-square = %4.5f    ',waldstat);
fprintf(1,'Probability = %4.6f \n',waldprob);