% This is an executable file which reads matlab data structures PITSgdp1, PITSpgdp1, PITSgdp4 and PITSpgdp4
% which include the stored versions of the PITS for one-step-ahead and
% four-step ahead horizons.
% It further constructs figures consistent with Diebold, Gunther and Tay (1998) testing of uniformity
% via histograms. The 2.5th and 97.5th percentiles are calculated based on a normal approximation of
% binomial distributions.

clear all; close all
load PITSgdp1.mat
load PITSpgdp1.mat

vnamegdp = [vnamegdp(1:32), 'Simple Average', 'BMA-OLS', 'BMA', 'BMA-MC3', 'Factor', 'BVAR'];
vnamepgdp = [vnamepgdp(1:32), 'Simple Average', 'BMA-OLS', 'BMA', 'BMA-MC3', 'Factor', 'BVAR'];

% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% % Diebold, Gunther and Tay type plots
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
bin = 5;
m = size(densforecastgdp,2);
for i = [1 7 20 26];
    fn = figure;
    hs = histc(densforecastgdp(i,:),0:1/bin:1);
    bar(0:1/bin:1,hs./m*100,'histc')
    hold on
    plot(0:1/bin:1,(1/5*100)*ones(bin+1,1), '-r','LineWidth',2)
    hold on
    plot(0:1/bin:1,(1/5+1.96*sqrt((1/5*(1-1/5))/m))*100*ones(bin+1,1), '--r','LineWidth',2)
    hold on
    plot(0:1/bin:1,(1/5-1.96*sqrt((1/5*(1-1/5))/m))*100*ones(bin+1,1), '--r','LineWidth',2)
    xlim([0 1])
    ylim([0 100])
    title('PIT for output growth - horizon(1)')
    legend(vnamegdp(:,i),'Location','Best')
    box on
    saveas(fn,strcat('pitgdpgrowth',char(vnamegdp(:,i))),'eps')
    close
end

for i = 1:size(pooldensforecastgdp,1)
    fn = figure;
    hs = histc(pooldensforecastgdp(i,:),0:1/bin:1);
    bar(0:1/bin:1,hs./m*100,'histc')
    hold on
    plot(0:1/bin:1,(1/5*100)*ones(bin+1,1), '-r','LineWidth',2)
    hold on
    plot(0:1/bin:1,(1/5+1.96*sqrt((1/5*(1-1/5))/m))*100*ones(bin+1,1), '--r','LineWidth',2)
    hold on
    plot(0:1/bin:1,(1/5-1.96*sqrt((1/5*(1-1/5))/m))*100*ones(bin+1,1), '--r','LineWidth',2)
    xlim([0 1])
    ylim([0 100])
    title('PIT for output growth - horizon(1)')
    legend(vnamegdp(:,32+i),'Location','Best')
    box on
    saveas(fn,strcat('pitgdpgrowth',char(vnamegdp(:,i+32))),'eps')
    close
end


fn = figure;
hs = histc(factordensforecastgdp,0:1/bin:1);
bar(0:1/bin:1,hs./m*100,'histc')
hold on
plot(0:1/bin:1,(1/5*100)*ones(bin+1,1), '-r','LineWidth',2)
hold on
plot(0:1/bin:1,(1/5+1.96*sqrt((1/5*(1-1/5))/m))*100*ones(bin+1,1), '--r','LineWidth',2)
hold on
plot(0:1/bin:1,(1/5-1.96*sqrt((1/5*(1-1/5))/m))*100*ones(bin+1,1), '--r','LineWidth',2)
xlim([0 1])
ylim([0 100])
title('PIT for output growth - horizon(1)')
legend(vnamegdp(:,end-1),'Location','Best')
box on
saveas(fn,strcat('pitgdpgrowth',char(vnamegdp(:,end-1))),'eps')
close

fn = figure;
hs = histc(vardensforecastgdp,0:1/bin:1);
bar(0:1/bin:1,hs./m*100,'histc')
hold on
plot(0:1/bin:1,(1/5*100)*ones(bin+1,1), '-r','LineWidth',2)
hold on
plot(0:1/bin:1,(1/5+1.96*sqrt((1/5*(1-1/5))/m))*100*ones(bin+1,1), '--r','LineWidth',2)
hold on
plot(0:1/bin:1,(1/5-1.96*sqrt((1/5*(1-1/5))/m))*100*ones(bin+1,1), '--r','LineWidth',2)
xlim([0 1])
ylim([0 100])
title('PIT for output growth - horizon(1)')
legend(vnamegdp(:,end),'Location','Best')
box on
saveas(fn,strcat('pitgdpgrowth',char(vnamegdp(:,end))),'eps')
close

m = size(densforecastpgdp,2);
for i = [1 7 21 26];
    fn = figure;
    hs = histc(densforecastpgdp(i,:),0:1/bin:1);
    bar(0:1/bin:1,hs./m*100,'histc')
    hold on
    plot(0:1/bin:1,(1/5*100)*ones(bin+1,1), '-r','LineWidth',2)
    hold on
    plot(0:1/bin:1,(1/5+1.96*sqrt((1/5*(1-1/5))/m))*100*ones(bin+1,1), '--r','LineWidth',2)
    hold on
    plot(0:1/bin:1,(1/5-1.96*sqrt((1/5*(1-1/5))/m))*100*ones(bin+1,1), '--r','LineWidth',2)
    xlim([0 1])
    ylim([0 100])
    title('PIT for inflation - horizon(1)')
    legend(vnamepgdp(:,i),'Location','Best')
    box on
    saveas(fn,strcat('pitpgdpgrowth',char(vnamepgdp(:,i))),'eps')
    close
end

for i = 1:size(pooldensforecastpgdp,1)
    fn = figure;
    hs = histc(pooldensforecastpgdp(i,:),0:1/bin:1);
    bar(0:1/bin:1,hs./m*100,'histc')
    hold on
    plot(0:1/bin:1,(1/5*100)*ones(bin+1,1), '-r','LineWidth',2)
    hold on
    plot(0:1/bin:1,(1/5+1.96*sqrt((1/5*(1-1/5))/m))*100*ones(bin+1,1), '--r','LineWidth',2)
    hold on
    plot(0:1/bin:1,(1/5-1.96*sqrt((1/5*(1-1/5))/m))*100*ones(bin+1,1), '--r','LineWidth',2)
    xlim([0 1])
    ylim([0 100])
    title('PIT for inflation - horizon(1)')
    legend(vnamepgdp(:,32+i),'Location','Best')
    box on
    saveas(fn,strcat('pitpgdpgrowth',char(vnamepgdp(:,i+32))),'eps')
    close
end

fn = figure;
hs = histc(factordensforecastpgdp,0:1/bin:1);
bar(0:1/bin:1,hs./m*100,'histc')
hold on
plot(0:1/bin:1,(1/5*100)*ones(bin+1,1), '-r','LineWidth',2)
hold on
plot(0:1/bin:1,(1/5+1.96*sqrt((1/5*(1-1/5))/m))*100*ones(bin+1,1), '--r','LineWidth',2)
hold on
plot(0:1/bin:1,(1/5-1.96*sqrt((1/5*(1-1/5))/m))*100*ones(bin+1,1), '--r','LineWidth',2)
xlim([0 1])
ylim([0 100])
title('PIT for inflation - horizon(1)')
legend(vnamepgdp(:,end-1),'Location','Best')
box on
saveas(fn,strcat('pitpgdpgrowth',char(vnamepgdp(:,end-1))),'eps')
close

fn = figure;
hs = histc(vardensforecastpgdp,0:1/bin:1);
bar(0:1/bin:1,hs./m*100,'histc')
hold on
plot(0:1/bin:1,(1/5*100)*ones(bin+1,1), '-r','LineWidth',2)
hold on
plot(0:1/bin:1,(1/5+1.96*sqrt((1/5*(1-1/5))/m))*100*ones(bin+1,1), '--r','LineWidth',2)
hold on
plot(0:1/bin:1,(1/5-1.96*sqrt((1/5*(1-1/5))/m))*100*ones(bin+1,1), '--r','LineWidth',2)
xlim([0 1])
ylim([0 100])
title('PIT for inflation - horizon(1)')
legend(vnamepgdp(:,end),'Location','Best')
box on
saveas(fn,strcat('pitpgdpgrowth',char(vnamepgdp(:,end))),'eps')
close