%Time Series. Econ 232C. homework 2 problem 2 Anton Cheremukhin, Christoph Bertsch and Martin Szydlowski
close all;
clear all;
clc;
g=15;
% get data
[GDP] = textread('realGDP.txt','%f');
[Year2, Month2, Day2, H] = textread('AWHNONAG.txt','%d-%d-%d %f');
H=H(1:size(H,1)-3);
H=log(H);
%hours convert monthly to quarterly
for i=1:(size(H,1)/3)
    H2(i)=mean(H(3*i-2:3*i));
end
H2=H2';
T=size(GDP,1);
%construct productivity
prod=log(GDP)-H2;
%growth rates
dprod=prod(2:T)-prod(1:T-1);
dhours=H2(2:T)-H2(1:T-1);
T=T-1;
%demean
%dprod=dprod-mean(dprod);
%dhours=dhours-mean(dhours);
%plot
figure(1)
plot(1:T,dprod,1:T,dhours);
legend('productivity','hours')
%correlation
%corrcoef(dprod,dhours)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%ACF-PACF for productivity
%compute ACF
h=dprod'-mean(dprod);    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 productivity');   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 productivity');  hold off;

%ACF-PACF for hours
%compute ACF
h=dhours'-mean(dhours);    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(3);  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 hours');   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 hours');  hold off;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%55
%VAR RF - OLS - regression-by-regression estimation
p=2;
%construct xs 
X=ones(T-p,1);
for i=1:p
    X=[X dprod(p+1-i:T-i) dhours(p+1-i:T-i)];
end
%OLS 1
Y=dprod(p+1:T);
beta1=inv(X'*X)*(X'*Y);
err1=Y-X*beta1;
%OLS 2
Y=dhours(p+1:T);
beta2=inv(X'*X)*(X'*Y);
err2=Y-X*beta2;
%construct As
cA=[beta1';beta2'];
c=cA(:,1);
A=zeros(2,2,g);
for i=1:p
    A(:,:,i)=cA(:,2*i:2*i+1);
end
%construct Bs
B=zeros(2,2,g);
B(:,:,1)=eye(2);
for i=1:g-1
    for j=1:i
        B(:,:,i+1)=B(:,:,i+1)+B(:,:,i+1-j)*A(:,:,j);
    end
end
%sum them up
B2=zeros(2);
for i=1:g
    B2=B2+B(:,:,i);
end
%construct Omega 
e=[err1, err2];     Om=zeros(2,2);
for i=1:size(e,1)
    Om=Om+e(i,:)'*e(i,:)/(size(e,1)-3*p+1);
end
%construct the lower-triangular choletsky matrix using a LR restriction
Om2=B2*Om*B2';      P=chol(Om2); %P=P';      %D=P.*eye(2);        A2=inv(D)*P;
A2=inv(B2)*P;
IRF=zeros(2,2,g);
for i=1:g
    IRF(:,:,i)=B(:,:,i)*A2;
end
IRF2(1,:)=reshape(IRF(1,1,:),[1 size(IRF,3)]);
IRF2(2,:)=reshape(IRF(1,2,:),[1 size(IRF,3)]);
IRF2(3,:)=reshape(IRF(2,1,:),[1 size(IRF,3)]);
IRF2(4,:)=reshape(IRF(2,2,:),[1 size(IRF,3)]);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%bootstrap - to bias-correct
N=1000;
k=1;
for k=1:N
%resample errors
%for i=1:size(e,1)
%    j=int8(rand*size(e,1)+1);
%    e_(i,:)=e(j,:);
%end
%reconstruct dprod and dhours
dprod2=zeros(size(dprod)+p);
dhours2=zeros(size(dhours)+p);
for t=(p+1):(size(dprod,1)+p)
    X=1;
    for i=1:p
        X=[X dprod2(t-i) dhours2(t-i)];
    end
    j=int8(rand*size(e,1)+1);
    e_=e(j,:);
    dprod2(t)=X*cA(1,:)'+e_(1);
    dhours2(t)=X*cA(2,:)'+e_(2);
end
dprod2=dprod2(p+1:size(dprod)+p)';
dhours2=dhours2(p+1:size(dhours)+p)';

%construct xs 
X=ones(T-p,1);
for i=1:p
    X=[X dprod2(p+1-i:T-i) dhours2(p+1-i:T-i)];
end
%OLS 1
Y=dprod2(p+1:T);
beta1=inv(X'*X)*(X'*Y);
err1=Y-X*beta1;
%OLS 2
Y=dhours2(p+1:T);
beta2=inv(X'*X)*(X'*Y);
err2=Y-X*beta2;
%construct cA
cA3(:,:,k)=[beta1';beta2'];
end
%the bias-correction procedure
for i=1:size(cA3,1)
    for j=1:size(cA3,2)
        cA(i,j)=2*cA(i,j)-mean(reshape(cA3(i,j,:),[1,size(cA3,3)]));
    end
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%recalculate the IRFs
X=ones(T-p,1);
for i=1:p
    X=[X dprod(p+1-i:T-i) dhours(p+1-i:T-i)];
end
%OLS 1
Y=dprod(p+1:T);
beta1=cA(1,:)';
err1=Y-X*beta1;
Y=dhours(p+1:T);
beta2=cA(2,:)';
err2=Y-X*beta2;
%construct As
A=zeros(2,2,g);
for i=1:p
    A(:,:,i)=cA(:,2*i:2*i+1);
end
%construct Bs
B=zeros(2,2,g);
B(:,:,1)=eye(2);
for i=1:g-1
    for j=1:i
        B(:,:,i+1)=B(:,:,i+1)+B(:,:,i+1-j)*A(:,:,j);
    end
end
%sum them up
B2=zeros(2);
for i=1:g
    B2=B2+B(:,:,i);
end
%construct Omega 
e=[err1, err2];     Om=zeros(2,2);
for i=1:size(e,1)
    Om=Om+e(i,:)'*e(i,:)/(size(e,1)-3*p+1);
end
%construct the lower-triangular choletsky matrix using a LR restriction
Om2=B2*Om*B2';      P=chol(Om2); %P=P';      %D=P.*eye(2);        A2=inv(D)*P;
A2=inv(B2)*P;   %=Co
IRF=zeros(2,2,g);
for i=1:g
    IRF(:,:,i)=B(:,:,i)*A2;
end
IRF2(1,:)=reshape(IRF(1,1,:),[1 size(IRF,3)]);
IRF2(2,:)=reshape(IRF(1,2,:),[1 size(IRF,3)]);
IRF2(3,:)=reshape(IRF(2,1,:),[1 size(IRF,3)]);
IRF2(4,:)=reshape(IRF(2,2,:),[1 size(IRF,3)]);
% orthogonalized errors
u=(inv(A2)*e')';
inv(A2)
inv(A2)*cA

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%55
%TEST ERRORS
%compute ACF
h=u(:,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(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 shocks to productivity');   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 dprod(1:(T-P),1) dhours(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))]);

%compute ACF
h=u(:,2)'; 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(5);  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  shocks to hours');   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 dprod(1:(T-P),1) dhours(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))]);



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%5
%bootstrap to estimate confidence intervals
N=1000;
k=1;
for k=1:N
%resample errors
%for i=1:size(e,1)
%    j=int8(rand*size(e,1)+1);
%    e_(i,:)=e(j,:);
%end
%reconstruct dprod and dhours
dprod2=zeros(size(dprod)+p);
dhours2=zeros(size(dhours)+p);
for t=(p+1):(size(dprod,1)+p)
    X=1;
    for i=1:p
        X=[X dprod2(t-i) dhours2(t-i)];
    end
    j=int8(rand*size(e,1)+1);
    e_=e(j,:);
    dprod2(t)=X*cA(1,:)'+e_(1);
    dhours2(t)=X*cA(2,:)'+e_(2);
end
dprod2=dprod2(p+1:size(dprod)+p)';
dhours2=dhours2(p+1:size(dhours)+p)';

%construct xs 
X=ones(T-p,1);
for i=1:p
    X=[X dprod2(p+1-i:T-i) dhours2(p+1-i:T-i)];
end
%OLS 1
Y=dprod2(p+1:T);
beta1=inv(X'*X)*(X'*Y);
err1=Y-X*beta1;
%OLS 2
Y=dhours2(p+1:T);
beta2=inv(X'*X)*(X'*Y);
err2=Y-X*beta2;
%construct As
cA2=[beta1';beta2'];
c=cA2(:,1);
A=zeros(2,2,g);
for i=1:p
    A(:,:,i)=cA2(:,2*i:2*i+1);
end
%construct Bs
B=zeros(2,2,g);
B(:,:,1)=eye(2);
for i=1:g-1
    for j=1:i
        B(:,:,i+1)=B(:,:,i+1)+B(:,:,i+1-j)*A(:,:,j);
    end
end
%sum them up
B2=zeros(2);
for i=1:g
    B2=B2+B(:,:,i);
end
%construct Omega 
e2=[err1, err2];     Om=zeros(2,2);
for i=1:size(e2,1)
    Om=Om+e2(i,:)'*e2(i,:)/(size(e2,1)-3*p+1);
end
%construct the lower-triangular choletsky matrix using a LR restriction
Om2=B2*Om*B2';      P=chol(Om2); %P=P';      %D=P.*eye(2);        A2=inv(D)*P;
A2=inv(B2)*P;
IRF=zeros(2,2,g);
for i=1:g
    IRF(:,:,i)=B(:,:,i)*A2;
end
IRF3(1,:,k)=reshape(IRF(1,1,:),[1 size(IRF,3)]);
IRF3(2,:,k)=reshape(IRF(1,2,:),[1 size(IRF,3)]);
IRF3(3,:,k)=reshape(IRF(2,1,:),[1 size(IRF,3)]);
IRF3(4,:,k)=reshape(IRF(2,2,:),[1 size(IRF,3)]);
end

for i=1:size(IRF3,1)
    for j=1:size(IRF3,2)
        temp=sort(reshape(IRF3(i,j,:),[1,size(IRF3,3)]));
        IRF2l(i,j)=temp(N*0.05);
        IRF2u(i,j)=temp(N*0.95);
    end
end

%plot results
m=size(IRF2,2);
figure(6)
subplot(2,2,1)
plot(1:m,IRF2(1,:),1:m,IRF2u(1,:),1:m,IRF2l(1,:))
title('prod response to innov in prod');
subplot(2,2,2)
plot(1:m,IRF2(2,:),1:m,IRF2u(2,:),1:m,IRF2l(2,:))
title('prod response to innov in hours');
subplot(2,2,3)
plot(1:m,IRF2(3,:),1:m,IRF2u(3,:),1:m,IRF2l(3,:))
title('hours response to innov in prod');
subplot(2,2,4)
plot(1:m,IRF2(4,:),1:m,IRF2u(4,:),1:m,IRF2l(4,:))
title('hours response to innov in hours');
