% 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