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