% Time Series. Econ232C. Spring 2006. Homework 1. Problem 2. Anton Cheremukhin
clear all;
close all;
clc;
% get data
[Year1, Month1, Day1, Y] = textread('INDPRO.txt','%d-%d-%d %f');
[Year2, Month2, Day2, H] = textread('AWHNONAG.txt','%d-%d-%d %f');
%construct productivity
a=length(Y)-length(H)+1;  b=length(Y);  Y=Y(a:b,1);     P=log(Y./H); 
%plot time series
time=Year2+(Month2-1/2)/12;
ON=ones(size(P));
X=[ON time];
trend=inv(X'*X)*(X'*P);
Pdtr=P-X*trend;
figure(1);
plot(time, P, time, X*trend);
title('Labor Productivity, Monthly, 1964:1-2006:3');
T=size(P,1);
dP=P(2:T,1)-P(1:(T-1),1);
%figure(4);
%plot(time(2:T,1), dP);
%title('Labor Productivity Growth, Monthly, 1964:1-2006:3');

%compute ACF
h=dP'-dP'*ones(size(dP))/size(dP,1);    T=size(h,2);    L=24;
gamma0=(h*h')/T;
for j=1:L
    gamma(j)=( h(1,(j+1):T) * (h(1,1:(T-j)))')/(T-j);
    ster(j)=2/sqrt(T);
end
ACF=gamma/gamma0;   
figure(2);  subplot(1,2,1); stem_hadles=stem(1:L,ACF,'fill');
hold on;    plot_hadles=plot(1:L,ster,'r',1:L,-ster,'r');   title('ACF');   hold off;

%compute PACF
for j=1:L
    Y=h(1,(j+1):T)';        X=zeros((T-j),j);
    for k=1:j
        X(:,k)=h(1,(j+1-k):(T-k))';
    end
    beta=inv(X'*X)*(X'*Y);      alpha(j)=beta(j);       ster(j)=2/sqrt(T);
end
PACF=alpha;     subplot(1,2,2);     stem_hadles=stem(1:L,PACF,'fill');
hold on;    plot_hadles=plot(1:L,ster,'r',1:L,-ster,'r');   title('PACF');  hold off;

%estimate AR(1)
T=size(dP,1);   Y=dP(2:T,1);        ON=ones(T-1,1); X=[ON dP(1:(T-1),1)];

beta=inv(X'*X)*(X'*Y);
V=(((Y-X*beta)'*(Y-X*beta))/(size(X,1)-size(X,2)))*inv(X'*X);
GE=eye(size(X,2));  ste=sqrt(spdiags((GE*V*GE').*eye(size(X,2))));
t=GE*beta./ste; pv=2-2*tcdf(abs(t),size(X,1)-size(X,2));
disp('     c        AR(1)');    disp([beta ste t pv]'); %disp('R-squared');  
R_2=zeros(1,2);
R_2(1,1)=1-((Y-X*beta)'*(Y-X*beta))/((Y-(Y'*ON)/(ON'*ON))'*(Y-(Y'*ON)/(ON'*ON)));
R_2(1,2)=1-(1-R_2(1,1))*(size(X,1)-1)/(size(X,1)-size(X,2));    %disp(R_2);

GE=[0 1];   t=(GE*beta)./sqrt(spdiags((GE*V*GE')));
disp(['P-value of hypothesis: AR(1)=0 ' num2str(2-2*tcdf(abs(t),size(X,1)-size(X,2)))]);

%ster=2/sqrt(T);
%disp(['Asy standard error = 2/sqrt(T) ' num2str(ster)]);
%disp(['P-value of hypothesis: AR(1)=0 ' num2str(2-2*normcdf(abs(beta(2)/ster)))]);

resid=Y-X*beta;
sigma=sqrt(resid'*resid/(T-2));
figure(3);
plot(time(2:T,1), resid,time(2:T,1), 2*sigma,time(2:T,1), -2*sigma);
title('Residuals');

%compute ACF
h=resid';    T=size(h,2);    L=24;
gamma0=(h*h')/T;
for j=1:L
    gamma(j)=( h(1,(j+1):T) * (h(1,1:(T-j)))')/(T-j);
    ster(j)=2/sqrt(T);
end
ACF=gamma/gamma0;   
figure(4);  subplot(1,2,1); stem_hadles=stem(1:L,ACF,'fill');
hold on;    plot_hadles=plot(1:L,ster,'r',1:L,-ster,'r');   title('ACF');   hold off;

%compute PACF
for j=1:L
    Y=h(1,(j+1):T)';        X=zeros((T-j),j);
    for k=1:j
        X(:,k)=h(1,(j+1-k):(T-k))';
    end
    beta=inv(X'*X)*(X'*Y);      alpha(j)=beta(j);       ster(j)=2/sqrt(T);
end
PACF=alpha;     subplot(1,2,2);     stem_hadles=stem(1:L,PACF,'fill');
hold on;    plot_hadles=plot(1:L,ster,'r',1:L,-ster,'r');   title('PACF');  hold off;

%Breusch-Godfrey test 
P=4; %according to ACF and PACF
Y=h(1,(P+1):T)';        X=zeros((T-P),P);   ON=ones(size(Y));
for k=1:P
    X(:,k)=h(1,(P+1-k):(T-k))';
end
X=[ON dP(1:(T-P),1) X];
beta=inv(X'*X)*(X'*Y);
V=(((Y-X*beta)'*(Y-X*beta))/(size(X,1)-size(X,2)))*inv(X'*X);
GE=eye(size(X,2));  ste=sqrt(spdiags((GE*V*GE').*eye(size(X,2))));
t=GE*beta./ste; pv=2-2*tcdf(abs(t),size(X,1)-size(X,2));
disp('     c        y(-1)     u(-1)       u(-2)       u(-3)       u(-4)   ');    
disp([beta ste t pv]'); disp('R-squared');  R_2=zeros(1,2);
R_2(1,1)=1-((Y-X*beta)'*(Y-X*beta))/((Y-(Y'*ON)/(ON'*ON))'*(Y-(Y'*ON)/(ON'*ON)));
R_2(1,2)=1-(1-R_2(1,1))*(size(X,1)-1)/(size(X,1)-size(X,2));    disp(R_2(1));
disp(T*R_2(1));
%GE=[0 0 1 1 1 1];   t=(GE*beta)./sqrt(spdiags((GE*V*GE')));
%disp(['P-value of hypothesis: No serial correlation ' num2str(2-2*tcdf(abs(t),size(X,1)-size(X,2)))]);
disp(['P-value of hypothesis: No serial correlation ' num2str(1-chi2cdf(T*R_2(1),P))]);

%Andrews 1993 based on Chow 1960
T=size(dP,1);   Y=dP(2:T,1);        ON=ones(T-1,1); X=[ON dP(1:(T-1),1)];

a1=25;
b1=T-a1;
for T0=a1:b1
    ON=ones(T-1,1); 
    X=[ON(1:T0,1)      zeros(T0,1)        dP(1:T0,1)
       zeros(T-T0-1,1) ON((T0+1):(T-1),1) dP((T0+1):(T-1),1)];
%OLS estimate
    beta=inv(X'*X)*(X'*Y);
%heteroscedasticity-robust errors
    err=(Y-X*beta);     err2=err.^2;
    V2=(inv(X'*X)*(X'*diag(err2)*X)*inv(X'*X));
%Wald asymptotic test for 1 restriction
    R=[1 -1 0];
    W(T0-a1+1)=(R*beta)'*inv(R*V2*R')*(R*beta);
    SSR(T0-a1+1)=err'*err;
%supF-test
    X2=[ON dP(1:(T-1),1)];
    beta2=inv(X2'*X2)*(X2'*Y);
    err2=(Y-X2*beta2);
    F(T0-a1+1)=(err2'*err2-err'*err)*(T-3)/(err'*err);
end

figure(5)
plot(time(a1:b1),W,time(a1:b1),F,'--',time(a1:b1),8.13,time(a1:b1),9.71,time(a1:b1),13.17,...
               'LineWidth',2)
legend('Heteroscedasticity-robust : Wald','Heteroscedasticity non-robust : F*k');
figure(6)
plot(time(a1:b1),SSR)
legend('SSR unrestricted');
[g1 f]=max(F);
[g f1]=max(-SSR);
if g1>8.13
    disp('break N1: ');
    disp(['Value of supF-statistic ' num2str(g1)]);
    d=(time(a1+f1)*12-1/2)/12-fix((time(a1+f1)*12-1/2)/12);
    disp(['break in ' num2str(fix((time(a1+f1)*12-1/2)/12)) ':' num2str(d*12)]);

%break N2
dP=dP(a1+f1+1:T);
T=size(dP,1);   Y=dP(2:T,1);        ON=ones(T-1,1); X=[ON dP(1:(T-1),1)];

a2=24;
b2=T-a2;
for T0=a2:b2
    ON=ones(T-1,1); 
    X=[ON(1:T0,1)      zeros(T0,1)        dP(1:T0,1)
       zeros(T-T0-1,1) ON((T0+1):(T-1),1) dP((T0+1):(T-1),1)];
%OLS estimate
    beta=inv(X'*X)*(X'*Y);
%heteroscedasticity-robust errors
    err=(Y-X*beta);     err2=err.^2;
    V2=(inv(X'*X)*(X'*diag(err2)*X)*inv(X'*X));
%Wald asymptotic test for 1 restriction
    R=[1 -1 0];
    W2(T0-a2+1)=(R*beta)'*inv(R*V2*R')*(R*beta);
    SSR2(T0-a2+1)=err'*err;
%supF-test
    X2=[ON dP(1:(T-1),1)];
    beta2=inv(X2'*X2)*(X2'*Y);
    err2=(Y-X2*beta2);
    F2(T0-a2+1)=(err2'*err2-err'*err)*(T-3)/(err'*err);
end

figure(7)
plot(time((a2+a1+f1):(b2+a1+f1)),W2,time((a2+a1+f1):(b2+a1+f1)),F2,'--',time((a2+a1+f1):(b2+a1+f1)),8.13,time((a2+a1+f1):(b2+a1+f1)),9.71,time((a2+a1+f1):(b2+a1+f1)),13.17,...
               'LineWidth',2)
legend('Heteroscedasticity-robust : Wald','Heteroscedasticity non-robust : F*k');
figure(8)
plot(time((a2+a1+f1):(b2+a1+f1)),SSR2)
legend('SSR unrestricted');
[g2 f]=max(F2);
[g f2]=max(-SSR2);
if g2>8.13
    disp('break N2: ')
    disp(['Value of supF-statistic ' num2str(g2)]);
    d=(time(a1+f1+a2+f2)*12-1/2)/12-fix((time(a1+f1+a2+f2)*12-1/2)/12);
    disp(['break in ' num2str(fix((time(a1+f1+a2+f2)*12-1/2)/12)) ':' num2str(d*12)]);
else
    disp('no break N2 detected ');
end
else
    disp('no break N1 detected ');    
end