% This is the main executable file replicating the third and fourth columns
% of Table 4 and Figures 1-3.
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:(transf(:,end)));
 
% if plots are not needed then this indicators could be set to zero
wanttoplot1 = 1;
wanttoplot2 = 1;

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

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

table = [];
for i = (h+1)
    seldata = data(:,ismember(transf,(i-1)));
    
    if wanttoplot1 == 1
        fn = figure;
        plot(pdate,seldata(:,1), 'LineWidth',2)
        hold on
        plot(pdate,seldata(:,2), '--g', 'LineWidth',2)
        hold on
        plot(pdate,seldata(:,3),':k','LineWidth',2)
        if (i-1 < 5)
            hold on
            plot(pdate,seldata(:,4), '-.r', 'LineWidth',2)
            legend('actual','Greenbook','BCEI','SPF','Location','best')
        else
            legend('actual','Greenbook','BCEI','Location','best')
        end
        title(strcat('Inflation at quarter h=',num2str(i-1)),'Fontsize',14);
        set(findall(gcf,'Type','text'),'FontSize',14)
        set(findall(gcf,'Type','axes'),'FontSize',14)
        axis tight
        saveas(fn,strcat('datafigure',num2str(i-1),'.eps'),'epsc')
        hold off
    end
    
    if (i-1) > 4
        N = 2;
    else
        N = 3;
    end

    for n = 1:N
        j = n + 1;
        truncindex = isfinite(seldata(:,[1 j]));
        truncdata = seldata(truncindex(:,2),[1 j]);
        truncdate = pdate(truncindex(:,2));
        
        gbn(i,n) = size(truncdata,1);
        if (gbn(i,n) < R)
            disp(strcat('The rolling window size has been adjusted to the minimum number of observations for n=', num2str(n), ' and h=',num2str(i-1)))
            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(:,2);
        x = [ones(gbn(i,n),1) truncdata(:,2)];

        % Traditional Test Results
        wald(i,n) = olsWaldtestNW(y,x);
        waldp(i,n) = 1 - chi2cdf(wald(i,n),2);
        
        nlag = round(size(y,1)^(1/4));
        nwresult = nwest(y,x,nlag);
        varbetahat = nwresult.vcv;
        sehat = sqrt(diag(varbetahat));

        % Fluctuation Test Results
        if n == 1
            vname = strcat('at quarter h=',num2str(i-1));
        elseif n == 2;
            vname = strcat('at quarter h=',num2str(i-1));
        elseif n == 3;
            vname = strcat('at quarter h=',num2str(i-1));
        end
        
        result = FluctuationRationalityRecNW(m,y,x,[]);
        
        table = [table; (i-1)   result.MZ   result.MZrej];
        resultstoplot(i,n,1:length(result.wseries)) = result.wseries;    
        critvalues(i,n) = result.cv;
    end
    
    % the critical values would depend on a sample size so one should be
    % careful in choosing which one to report. We report the max.
    if wanttoplot2==1; 
        fn = figure; 
        plot(squeeze(ptruncdate(i,1,:)),squeeze(resultstoplot(i,1,:)),'--g','LineWidth',2)
        hold on
        plot(squeeze(ptruncdate(i,2,:)),squeeze(resultstoplot(i,2,:)),':k','LineWidth',2)
        if (i-1 < 5)
            hold on
            plot(squeeze(ptruncdate(i,3,:)),squeeze(resultstoplot(i,3,:)), '-.r', 'LineWidth',2)
            hold on
            plot(squeeze(ptruncdate(i,3,:)),max(critvalues(i,:))*ones(size(ptruncdate(i,3,:),3),1),'-b','LineWidth',2)
            legend('Greenbook','BCEI','SPF','5% crit. value','Location','best')
        else
            hold on
            plot(squeeze(ptruncdate(i,1,:)),max(critvalues(i,:))*ones(size(ptruncdate(i,1,:),3),1),'-b','LineWidth',2)
            legend('Greenbook','BCEI','5% crit. value','Location','best')
        end
        title(['Forecast 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('testfigure',num2str(i-1),'.eps'),'epsc')
        hold off
    end
end

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

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

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

'Fluctuation Test Results'
table

'Program Complete'