clear all; close all; clc; % get data [Year1, Month1, Day1, Um] = textread('UNEMP.txt','%d-%d-%d %f'); [Year2, Quarter2, Uf_1, Uf0, Uf1, Uf2] = textread('SFPMedian.txt','%d %d %f %f %f %f'); Um=Um(253:(size(Year1,1)-2)); for k=1:(size(Um,1)/3) Uq(k)=(Um(3*k-2)+Um(3*k-1)+Um(3*k))/3; end Uq=Uq'; % real data - 1969:2 till 2006:1 - we had 1969:1 till 2006:1 - 149 points Uq=Uq(2:149); % 1 period ahead forecast - 1969:1 till 2005:4 - we have 1968:4 till 2006:2 - 151 points Uf1=Uf1(2:149); %2:149 % 2 period ahead forecast - 1968:4 till 2005:3 Uf2=Uf2(1:148); %1:148 % now we have 148 points of quarterly data - real data versus two forecasts time=0.125:0.25:(148/4); time=1969.25+time; % plot time series figure(1) plot(time,Uq,time,Uf1,time,Uf2) legend('true','1 period ahead','2 periods ahead') e1=Uq-Uf1; e2=Uq-Uf2; % plot errors figure(2) plot(time,e1,time,e2) legend('1 period ahead','2 periods ahead') % unbiasedness test: regress error on information, known at the day of % forecast: e(i)=a+b*y(t-i), i=1,2 H0:a=b=0 % case i=1 disp('1 step ahead'); T=size(e1,1); Y=e1(2:T); X=[ones(size(Y))];% Uf1(1:T-1)]; Beta_hat=inv(X'*X)*(X'*Y); u_hat=(Y-X*Beta_hat); SSR=u_hat'*u_hat; %using White's heteroskedasticity consistent standard errors %Var_Beta=inv(X'*X)*((X.*[u_hat, u_hat])'*(X.*[u_hat, u_hat]))*inv(X'*X); Var_Beta=inv(X'*X)*((X.*u_hat)'*(X.*u_hat))*inv(X'*X); %not using HAC: i-1=0 %Var_Beta=inv(X'*X)*((X.*[u_hat, u_hat])'*(X.*[u_hat, u_hat]))*inv(X'*X); GE=eye(size(X,2)); ste=sqrt(spdiags((GE*Var_Beta*GE').*eye(size(X,2)))); t=GE*Beta_hat./ste; pv=2-2*tcdf(abs(t),size(X,1)-size(X,2)); disp(' a b '); disp([Beta_hat ste t pv]'); % Wald test with White R=1;%[1 0; 0 1]; %Wald test W=(R*Beta_hat)'*inv(R*Var_Beta*(R'))*(R*Beta_hat) disp(['P-value of hypothesis: unbiased ' num2str(1-chi2cdf(W,1))]); %2 %Breusch-Godfrey test P=5; %according to ACF and PACF h=e1'; T=size(h,2); 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 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 u(-1) u(-2) u(-3) u(-4) u(-5)'); 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))]); %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % case i=2 disp('2 step ahead'); T=size(e2,1); Y=e2(3:T); X=[ones(size(Y))];% Uf2(1:T-2)]; T=T-2; Beta_hat=inv(X'*X)*(X'*Y); u_hat=(Y-X*Beta_hat); SSR=u_hat'*u_hat; %not using White's heteroskedasticity consistent standard errors %Var_Beta=inv(X'*X)*((X.*[u_hat, u_hat])'*(X.*[u_hat, u_hat]))*inv(X'*X); %using HAC: i-1=1 %Var_Beta=inv(X'*X)*((X.*[u_hat, u_hat])'*(X.*[u_hat, u_hat])+0.5*((X(2:T,:).*[u_hat(2:T), u_hat(2:T)])'*(X(1:(T-1),:).*[u_hat(1:(T-1)), u_hat(1:(T-1))])+(X(1:T-1,:).*[u_hat(1:T-1), u_hat(1:T-1)])'*(X(2:T,:).*[u_hat(2:T), u_hat(2:T)])))*inv(X'*X); Var_Beta=inv(X'*X)*((X.*u_hat)'*(X.*u_hat)+0.5*((X(2:T,:).*u_hat(2:T))'*(X(1:(T-1),:).*u_hat(1:(T-1)))+(X(1:T-1,:).*u_hat(1:T-1))'*(X(2:T,:).*u_hat(2:T))))*inv(X'*X); GE=eye(size(X,2)); ste=sqrt(spdiags((GE*Var_Beta*GE').*eye(size(X,2)))); t=GE*Beta_hat./ste; pv=2-2*tcdf(abs(t),size(X,1)-size(X,2)); disp(' a b '); disp([Beta_hat ste t pv]'); % Wald test with HAC R=1;%[1 0; 0 1]; %Wald test W=(R*Beta_hat)'*inv(R*Var_Beta*(R'))*(R*Beta_hat) disp(['P-value of hypothesis: unbiased ' num2str(1-chi2cdf(W,1))]); %2 %compute ACF dP=e2; 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(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'); 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; %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %forecasts of changes T=size(Uq,1); dfor=Uf1(2:T)-Uq(1:T-1); dtru=Uq(2:T)-Uq(1:T-1); figure(4) plot(time(2:T),dfor,time(2:T),dtru) legend('forcasted change','true change') % sign predictability test L=1(sign(dfor)==sign(dtru)) %construct dummies sign(dfor)= a + b*sign(dtru) disp('test for sign predictability'); Y=(dfor>0); X=(dtru>0); X=[ones(size(X)) X]; Beta_hat=inv(X'*X)*(X'*Y); u_hat=(Y-X*Beta_hat); SSR=u_hat'*u_hat; %using White's heteroskedasticity consistent standard errors Var_Beta=inv(X'*X)*((X.*[u_hat, u_hat])'*(X.*[u_hat, u_hat]))*inv(X'*X); %not using HAC: i-1=0 %Var_Beta=inv(X'*X)*((X.*[u_hat, u_hat])'*(X.*[u_hat, u_hat]))*inv(X'*X); GE=eye(size(X,2)); ste=sqrt(spdiags((GE*Var_Beta*GE').*eye(size(X,2)))); t=GE*Beta_hat./ste; pv=2-2*tcdf(abs(t),size(X,1)-size(X,2)); disp(' a b '); disp([Beta_hat ste t pv]'); % Wald test with White H0: b=0 R=[0 1]; %Wald test W=(R*Beta_hat)'*inv(R*Var_Beta*(R'))*(R*Beta_hat) disp(['P-value of hypothesis: does not predict the sign of change ' num2str(1-chi2cdf(W,1))]); % set window to m=50 m=50; a=1; b=a+m-1; % estimate AR(p) using first m observations for p=1:20 Y=Uq(a+p:b); X=ones(size(Y)); for i=1:p X=[X Uq(a+p-i:b-i)]; end beta=inv(X'*X)*(X'*Y); u_hat=(Y-X*beta); SSR=u_hat'*u_hat; % calculate BIC BIC(p)=(p+1)*log(m-p)+(m-p)*log(SSR/(m-p)); end % choose p by minimum BIC figure(5) plot(1:20,BIC) [hjg,p]=min(BIC); %rolling window estimation and forecasting % Ufm1 - new forecast m=50; T=size(Uq,1); for a=1:T-m b=a+m-1; % estimate AR(p) using m observations Y=Uq(a+p:b); X=ones(size(Y)); for i=1:p X=[X Uq(a+p-i:b-i)]; end beta=inv(X'*X)*(X'*Y); %forecast X1=1; for i=1:p X1=[X1 Uq(a+m-i)]; end Ufm1(a)=X1*beta; end Ufm1=Ufm1'; Ufn1=Uf1(m+1:T); Uq1=Uq(m+1:T); %plot my forecast and SPF forecast figure(6) plot(time(m+1:T),Ufm1,time(m+1:T),Ufn1,time(m+1:T),Uq1) legend('AR(2) forecast','SPF forecast','true values'); em=Ufm1-Uq1; Lm=em.^2; en=Ufn1-Uq1; Ln=en.^2; dL=Ln-Lm; err=dL-mean(dL); % HAC variance of dL with 6 lags inside q=8; T2=size(err,1); s=err'*err; for i=1:q s=s+(1-i/(q+1))*2*err(1:T2-i)'*err(1+i:T2); end disp('Diebold Mariano test statistic'); DM=T2*mean(dL)/sqrt(s) disp(['P-value of hypothesis: AR(p) forecast is no better than SPF forecast ' num2str(1-1*normcdf(DM))]); figure(7) plot(time(m+1:T),dL,time(m+1:T),mean(dL)*ones(size(dL))) legend('dL(t)','mean')