% This is the main executable file replicating the third and fourth columns
% of Table 4 and Figures 4-5.
clear all; close all; clc

[A,B]  = xlsread('DataExtended.xlsx', 'Inflation');
transf = A(1,2:end);
data = A(2:end,2:end);
date = A(2:end,1);
h = (0:1:4);

% if plots are not needed then this indicators could be set to zero
wanttoplot1 = 1;
wanttoplot2 = 1;

T = size(data,1);
R = 60;

pdate = [1968.75:0.25:2008.75]';
ptruncdate = NaN(h(1,end)+1,2,T-R+1);
resultstoplot = NaN(h(1,end)+1,2,T-R+1);
resultstoplotcoef = NaN(h(1,end)+1,2,T-R+1);
critvalues = NaN(h(1,end)+1,2);

table = [];
for i = (h+1)
    seldata = data(:,ismember(transf,(i-1)));
    if (i-1) > 4
        N = 1;
    else
        N = 2;
    end
    
    for n = 1:N
        j = n + 2;
        truncindex = isfinite(seldata(:,[1 2 j]));
        truncdata = seldata(sum(truncindex(:,2:end),2)==2,[1 2 j]);
        truncdate = pdate(sum(truncindex(:,2:end),2)==2);
        
        gbn(i,n) = size(truncdata,1);    
        if (gbn(i,n) < R)
            disp('The rolling window size has been adjusted to the minimum number of observations')
            m = gbn(i,n);
        else
            m = R;
        end
        
        ptruncdate(i,n,1:size(truncdate(m:end,1))) = truncdate(m:end,1);
        y = truncdata(:,1) - truncdata(:,end);
        x = truncdata(:,2);
        
        % Traditional Test Results
        wald(i,n) = olsWaldtestNW(y,x,[ones(gbn(i,n),1) truncdata(:,end)]);
        waldp(i,n) = 1 - chi2cdf(wald(i,n),1);
             
        % Fluctuation Test Results
        if n == 1
            vname = strcat('BCEI at quarter h=',num2str(i-1));
            fname = 'BCEI';
        elseif n == 2;
            vname = strcat('SPF at quarter h=',num2str(i-1));
            fname = 'SPF';
        end
        
        result = FluctuationRationalityRecNW(m,y,x,[ones(gbn(i,n),1) truncdata(:,end)]);   
        
        table = [table; i-1   result.MZ   result.MZrej];
        resultstoplot(i,n,1:length(result.wseries)) = result.wseries;    
        resultstoplotcoef(i,n,1:length(result.wseries)) = result.coef(:,1);    
        critvalues(i,n) = result.cv; 
    end
    
    if wanttoplot1==1; 
        fn = figure; 
        plot(squeeze(ptruncdate(i,1,:)),squeeze(resultstoplot(i,1,:)),':k','LineWidth',2)
        hold on
        plot(squeeze(ptruncdate(i,2,:)),squeeze(resultstoplot(i,2,:)),'-.r','LineWidth',2)
        hold on
        plot(squeeze(ptruncdate(i,2,:)),max(critvalues(i,:))*ones(size(ptruncdate(i,2,:),3),1),'-b','LineWidth',2)
        legend('BCEI','SPF','5% crit. value','Location','best')
        title(['Forecast error at quarter h=',num2str(i-1)]);
        set(findall(gcf,'Type','text'),'FontSize',14)
        set(findall(gcf,'Type','axes'),'FontSize',14)
        axis tight
        saveas(fn,strcat('regtestfigure',num2str(i-1),'.eps'),'epsc')
        hold off
    end

    if wanttoplot2==1; 
        fn = figure; 
        plot(squeeze(ptruncdate(i,1,:)),squeeze(resultstoplotcoef(i,1,:)),':k','LineWidth',2)
        hold on
        plot(squeeze(ptruncdate(i,2,:)),squeeze(resultstoplotcoef(i,2,:)),'-.r','LineWidth',2)
        legend('BCEI','SPF','Location','best')
        title(['Forecast error at quarter h=',num2str(i-1)]);
        set(findall(gcf,'Type','text'),'FontSize',14)
        set(findall(gcf,'Type','axes'),'FontSize',14)
        axis tight
        saveas(fn,strcat('regcoeffigure',num2str(i-1),'.eps'),'epsc')
        hold off
    end
end


'BCEI - Full Sample'
[h' wald(:,1) waldp(:,1) gbn(:,1)]

'SPF - Full Sample'
[h' wald(:,2) waldp(:,2) gbn(:,2)]

table

'Program Complete'