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