% 1.Forecast the conditioanl expectation of shadow rate and its variance at
% each point t. And plot the forecasted trajectory of shadow rate.
% 2.With the conditional mean and variance of future shadow rate, we
% calculate the probability of negative rate in the future.
% 3.We can also compose the confidence band of this probability using 50000
% draws from MCMC

clc
clear
close all
addpath(genpath('./functions'));
LB = [100, 207, 244]/255;  % light blue
% xdate = datetime(1971,10,1) + calmonths(0:606);
xdate = datetime(1971,11,1) + calmonths(0:605);

% load('./results/mcmc_US.mat','X','tau','Phi','Xbar','gamma','Omega','sigtau2');
load('./results/temp100000.mat','X','tau','Phi','Xbar','gamma','Omega','sigtau2');


[M, T] = size(tau);
mcmc_ind = M/2;
ltri = tril(ones(3)) ;
delta0 = 0;
delta1 = [1 1 1]';


%%
% 1-month ahead forecast at each point
% xdate1 = datetime(1971,11,1) + calmonths(0:606);
% xdate = datetime(1971,10,1) + calmonths(0:606);
xdate1 = datetime(1971,12,1) + calmonths(0:605);
ESR1 = nan(mcmc_ind,T);
varSR1 = nan(mcmc_ind,1);
% Ubound = nan(mcmc_ind,T);  % mean + 1.65sd
% Lbound = nan(mcmc_ind,T);  % mean - 1.65sd
% At each point, we can calculate the prob of SR<0
prob_SR1 = nan(mcmc_ind,T);
for i = 1:mcmc_ind
    if mod(i,1000)==0
        disp(i);
    end
    Xipt = X((mcmc_ind+i-1)*3+1:(mcmc_ind+i)*3,:);
    tauipt = tau(mcmc_ind+i,:);
    Phiipt = reshape(Phi(mcmc_ind+i,:),3,3);
    Xbaript = Xbar(mcmc_ind+i,:)';
    gammaipt = gamma(mcmc_ind+i,:)';
    Omegaipt = [Omega(mcmc_ind+i,1) Omega(mcmc_ind+i,2) Omega(mcmc_ind+i,3);...
        Omega(mcmc_ind+i,2) Omega(mcmc_ind+i,4) Omega(mcmc_ind+i,5);...
        Omega(mcmc_ind+i,3) Omega(mcmc_ind+i,5) Omega(mcmc_ind+i,6)];
    sig2ipt = sigtau2(mcmc_ind+i);
    [ESR1(i,:), varSR1(i)] = SRpredt(Xipt, tauipt, Phiipt, Xbaript,...
       gammaipt, Omegaipt, sig2ipt, delta0, delta1, 1);
%    Ubound(i,:) = ESR3(i,:) + 1.65*varSR3(i);
%    Lbound(i,:) = ESR3(i,:) - 1.65*varSR3(i);
   prob_SR1(i,:) = normcdf(0,ESR1(i,:), sqrt(varSR1(i)));
   
end

% plot 1-month ahead forecasted shadow rate with 90% CI
% y1 = max(Ubound,[],1);
% y2 = min(Lbound,[],1);
y1 = prctile(ESR1,5,1);
y2 = prctile(ESR1,95,1);
fill([xdate1 fliplr(xdate1)], [y1 fliplr(y2)],LB, 'LineStyle','none');
hold on 
plot(xdate1, mean(ESR1,1),'k');
legend({'90% CI','1-Month Forecast'},'AutoUpdate','off')
set(gca, 'FontSize',20, 'XTick',datetime([1975 1980 1985 1990 1995 2000 2005 2010 2015 2020],01,01),...
    'XTickLabel', {'1975', '1980', '1985', '1990', '1995', '2000', '2005', '2010', '2015', '2020'});
yline(0,'--');
% xline(datetime(2022,04,01),'-r')
title('1-Month Forecasted SR');

% Plot 90% confidence band of probability that SR<0
y1 = prctile(prob_SR1,5,1);
y2 = prctile(prob_SR1,95,1);
fill([xdate fliplr(xdate)], [y1 fliplr(y2)],LB, 'LineStyle','none');
hold on 
plot(xdate, mean(prob_SR1,1),'k');
legend({'90% CI','Mean'},'AutoUpdate','off', 'location','northwest')
set(gca, 'FontSize',20, 'XTick',datetime([1975 1980 1985 1990 1995 2000 2005 2010 2015 2020],01,01),...
    'XTickLabel', {'1975', '1980', '1985', '1990', '1995', '2000', '2005', '2010', '2015', '2020'});
yline(0,'--');
% xline(datetime(2022,04,01),'-r')
title('1-Month Forecast Probability');

%%
% 3-month ahead forecast
% xdate3 = datetime(1972,1,1) + calmonths(0:606);
xdate3 = datetime(1972,2,1) + calmonths(0:605);
ESR3 = nan(mcmc_ind,T);
varSR3 = nan(mcmc_ind,1);
Ubound = nan(mcmc_ind,T);  % mean + 1.65sd
Lbound = nan(mcmc_ind,T);  % mean - 1.65sd
% At each point, we can calculate the prob of SR<0
prob_SR3 = nan(mcmc_ind,T);
for i = 1:mcmc_ind
    if mod(i,1000)==0
        disp(i);
    end
    Xipt = X((mcmc_ind+i-1)*3+1:(mcmc_ind+i)*3,:);
    tauipt = tau(mcmc_ind+i,:);
    Phiipt = reshape(Phi(mcmc_ind+i,:),3,3);
    Xbaript = Xbar(mcmc_ind+i,:)';
    gammaipt = gamma(mcmc_ind+i,:)';
    Omegaipt = [Omega(mcmc_ind+i,1) Omega(mcmc_ind+i,2) Omega(mcmc_ind+i,3);...
        Omega(mcmc_ind+i,2) Omega(mcmc_ind+i,4) Omega(mcmc_ind+i,5);...
        Omega(mcmc_ind+i,3) Omega(mcmc_ind+i,5) Omega(mcmc_ind+i,6)];
    sig2ipt = sigtau2(mcmc_ind+i);
    [ESR3(i,:), varSR3(i)] = SRpredt(Xipt, tauipt, Phiipt, Xbaript,...
       gammaipt, Omegaipt, sig2ipt, delta0, delta1, 3);
   Ubound(i,:) = ESR3(i,:) + 1.65*varSR3(i);
   Lbound(i,:) = ESR3(i,:) - 1.65*varSR3(i);
   prob_SR3(i,:) = normcdf(0,ESR3(i,:), sqrt(varSR3(i)));
   
end

% plot 3-month ahead forecasted shadow rate with 90% CI
% y1 = max(Ubound,[],1);
% y2 = min(Lbound,[],1);
y1 = prctile(ESR3,5,1);
y2 = prctile(ESR3,95,1);
fill([xdate3 fliplr(xdate3)], [y1 fliplr(y2)],LB, 'LineStyle','none');
hold on 
plot(xdate3, mean(ESR3,1),'k');
legend({'90% CI','3-Month Forecast'},'AutoUpdate','off')
set(gca, 'FontSize',20, 'XTick',datetime([1975 1980 1985 1990 1995 2000 2005 2010 2015 2020],01,01),...
    'XTickLabel', {'1975', '1980', '1985', '1990', '1995', '2000', '2005', '2010', '2015', '2020'});
yline(0,'--');
% xline(datetime(2022,04,01),'-r')
title('3-Month Forecasted SR');

% Plot 90% confidence band of probability that SR<0
y1 = prctile(prob_SR3,5,1);
y2 = prctile(prob_SR3,95,1);
fill([xdate fliplr(xdate)], [y1 fliplr(y2)],LB, 'LineStyle','none');
hold on 
plot(xdate, mean(prob_SR3,1),'k');
legend({'90% CI','Mean'},'AutoUpdate','off', 'location','northwest')
set(gca, 'FontSize',20, 'XTick',datetime([1975 1980 1985 1990 1995 2000 2005 2010 2015 2020],01,01),...
    'XTickLabel', {'1975', '1980', '1985', '1990', '1995', '2000', '2005', '2010', '2015', '2020'});
yline(0,'--');
%xline(datetime(2022,04,01),'-r')
title('3-Month Forecast Probability');


%%
% 6-month ahead forecast
% xdate6 = datetime(1972,4,1) + calmonths(0:606);
xdate6 = datetime(1972,5,1) + calmonths(0:605);
ESR6 = nan(mcmc_ind,T);
varSR6 = nan(mcmc_ind,1);
% Ubound = nan(mcmc_ind,T);  % mean + 1.65sd
% Lbound = nan(mcmc_ind,T);  % mean - 1.65sd
% At each point, we can calculate the prob of SR<0
prob_SR6 = nan(mcmc_ind,T);
for i = 1:mcmc_ind
    if mod(i,1000)==0
        disp(i);
    end
    Xipt = X((mcmc_ind+i-1)*3+1:(mcmc_ind+i)*3,:);
    tauipt = tau(mcmc_ind+i,:);
    Phiipt = reshape(Phi(mcmc_ind+i,:),3,3);
    Xbaript = Xbar(mcmc_ind+i,:)';
    gammaipt = gamma(mcmc_ind+i,:)';
    Omegaipt = [Omega(mcmc_ind+i,1) Omega(mcmc_ind+i,2) Omega(mcmc_ind+i,3);...
        Omega(mcmc_ind+i,2) Omega(mcmc_ind+i,4) Omega(mcmc_ind+i,5);...
        Omega(mcmc_ind+i,3) Omega(mcmc_ind+i,5) Omega(mcmc_ind+i,6)];
    sig2ipt = sigtau2(mcmc_ind+i);
    [ESR6(i,:), varSR6(i)] = SRpredt(Xipt, tauipt, Phiipt, Xbaript,...
       gammaipt, Omegaipt, sig2ipt, delta0, delta1, 6);
%    Ubound(i,:) = ESR6(i,:) + 1.65*varSR6(i);
%    Lbound(i,:) = ESR6(i,:) - 1.65*varSR6(i);
   prob_SR6(i,:) = normcdf(0,ESR6(i,:), sqrt(varSR6(i)));
   
end


% plot 6-month ahead forecasted shadow rate with 90% CI
% y1 = max(Ubound,[],1);
% y2 = min(Lbound,[],1);
y1 = prctile(ESR6,5,1);
y2 = prctile(ESR6,95,1);
fill([xdate6 fliplr(xdate6)], [y1 fliplr(y2)],LB, 'LineStyle','none');
hold on 
plot(xdate6, mean(ESR6,1),'k');
legend({'90% CI','6-Month Forecast'},'AutoUpdate','off')
set(gca, 'FontSize',20, 'XTick',datetime([1975 1980 1985 1990 1995 2000 2005 2010 2015 2020],01,01),...
    'XTickLabel', {'1975', '1980', '1985', '1990', '1995', '2000', '2005', '2010', '2015', '2020'});
yline(0,'--');
% xline(datetime(2022,04,01),'-r')
title('6-Month Forecasted SR');

% Plot 90% confidence band of probability that SR<0
y1 = prctile(prob_SR6,5,1);
y2 = prctile(prob_SR6,95,1);
fill([xdate fliplr(xdate)], [y1 fliplr(y2)],LB, 'LineStyle','none');
hold on 
plot(xdate, mean(prob_SR6,1),'k');
legend({'90% CI','Mean'},'AutoUpdate','off','location','northwest')
set(gca, 'FontSize',20, 'XTick',datetime([1975 1980 1985 1990 1995 2000 2005 2010 2015 2020],01,01),...
    'XTickLabel', {'1975', '1980', '1985', '1990', '1995', '2000', '2005', '2010', '2015', '2020'});
yline(0,'--');
% xline(datetime(2022,04,01),'-r')
title('6-Month Forecast Probability');

%%
% 12-month ahead forecast
% xdate12 = datetime(1972,10,1) + calmonths(0:606);
xdate12 = datetime(1972,11,1) + calmonths(0:605);
ESR12 = nan(mcmc_ind,T);
varSR12 = nan(mcmc_ind,1);
% Ubound = nan(mcmc_ind,T);  % mean + 1.65sd
% Lbound = nan(mcmc_ind,T);  % mean - 1.65sd
% At each point, we can calculate the prob of SR<0
prob_SR12 = nan(mcmc_ind,T);
for i = 1:mcmc_ind
    if mod(i,1000)==0
        disp(i);
    end
    Xipt = X((mcmc_ind+i-1)*3+1:(mcmc_ind+i)*3,:);
    tauipt = tau(mcmc_ind+i,:);
    Phiipt = reshape(Phi(mcmc_ind+i,:),3,3);
    Xbaript = Xbar(mcmc_ind+i,:)';
    gammaipt = gamma(mcmc_ind+i,:)';
    Omegaipt = [Omega(mcmc_ind+i,1) Omega(mcmc_ind+i,2) Omega(mcmc_ind+i,3);...
        Omega(mcmc_ind+i,2) Omega(mcmc_ind+i,4) Omega(mcmc_ind+i,5);...
        Omega(mcmc_ind+i,3) Omega(mcmc_ind+i,5) Omega(mcmc_ind+i,6)];
    sig2ipt = sigtau2(mcmc_ind+i);
    [ESR12(i,:), varSR12(i)] = SRpredt(Xipt, tauipt, Phiipt, Xbaript,...
       gammaipt, Omegaipt, sig2ipt, delta0, delta1, 12);
%    Ubound(i,:) = ESR6(i,:) + 1.65*varSR6(i);
%    Lbound(i,:) = ESR6(i,:) - 1.65*varSR6(i);
   prob_SR12(i,:) = normcdf(0,ESR12(i,:), sqrt(varSR12(i)));
   
end

% plot 12-month ahead forecasted shadow rate with 90% CI
y1 = prctile(ESR12,5,1);
y2 = prctile(ESR12,95,1);
fill([xdate12 fliplr(xdate12)], [y1 fliplr(y2)],LB, 'LineStyle','none');
hold on 
plot(xdate12, mean(ESR12,1),'k');
legend({'90% CI','12-Month Forecast'},'AutoUpdate','off')
set(gca, 'FontSize',20, 'XTick',datetime([1975 1980 1985 1990 1995 2000 2005 2010 2015 2020],01,01),...
    'XTickLabel', {'1975', '1980', '1985', '1990', '1995', '2000', '2005', '2010', '2015', '2020'});
yline(0,'--');
% xline(datetime(2022,04,01),'-r')
title('12-Month Forecasted SR');

% Plot 90% confidence band of probability that SR<0
y1 = prctile(prob_SR12,5,1);
y2 = prctile(prob_SR12,95,1);
fill([xdate fliplr(xdate)], [y1 fliplr(y2)],LB, 'LineStyle','none');
hold on 
plot(xdate, mean(prob_SR12,1),'k');
legend({'90% CI','Mean'},'AutoUpdate','off','location','northwest')
set(gca, 'FontSize',20, 'XTick',datetime([1975 1980 1985 1990 1995 2000 2005 2010 2015 2020],01,01),...
    'XTickLabel', {'1975', '1980', '1985', '1990', '1995', '2000', '2005', '2010', '2015', '2020'});
yline(0,'--');
% xline(datetime(2022,04,01),'-r')
title('12-Month Forecast Probability');


%%
% 24-month ahead forecast
% xdate24 = datetime(1973,10,1) + calmonths(0:606);
xdate24 = datetime(1973,11,1) + calmonths(0:605);
ESR24 = nan(mcmc_ind,T);
varSR24 = nan(mcmc_ind,1);
% Ubound = nan(mcmc_ind,T);  % mean + 1.65sd
% Lbound = nan(mcmc_ind,T);  % mean - 1.65sd
% At each point, we can calculate the prob of SR<0
prob_SR24 = nan(mcmc_ind,T);
for i = 1:mcmc_ind
    if mod(i,1000)==0
        disp(i);
    end
    Xipt = X((mcmc_ind+i-1)*3+1:(mcmc_ind+i)*3,:);
    tauipt = tau(mcmc_ind+i,:);
    Phiipt = reshape(Phi(mcmc_ind+i,:),3,3);
    Xbaript = Xbar(mcmc_ind+i,:)';
    gammaipt = gamma(mcmc_ind+i,:)';
    Omegaipt = [Omega(mcmc_ind+i,1) Omega(mcmc_ind+i,2) Omega(mcmc_ind+i,3);...
        Omega(mcmc_ind+i,2) Omega(mcmc_ind+i,4) Omega(mcmc_ind+i,5);...
        Omega(mcmc_ind+i,3) Omega(mcmc_ind+i,5) Omega(mcmc_ind+i,6)];
    sig2ipt = sigtau2(mcmc_ind+i);
    [ESR24(i,:), varSR24(i)] = SRpredt(Xipt, tauipt, Phiipt, Xbaript,...
       gammaipt, Omegaipt, sig2ipt, delta0, delta1, 24);
%    Ubound(i,:) = ESR6(i,:) + 1.65*varSR6(i);
%    Lbound(i,:) = ESR6(i,:) - 1.65*varSR6(i);
   prob_SR24(i,:) = normcdf(0,ESR24(i,:), sqrt(varSR24(i)));   
end

% plot 24-month ahead forecasted shadow rate with 90% CI
y1 = prctile(ESR24,5,1);
y2 = prctile(ESR24,95,1);
fill([xdate24 fliplr(xdate24)], [y1 fliplr(y2)],LB, 'LineStyle','none');
hold on 
plot(xdate24, mean(ESR24,1),'k');
legend({'90% CI','24-Month Forecast'},'AutoUpdate','off')
set(gca, 'FontSize',20, 'XTick',datetime([1975 1980 1985 1990 1995 2000 2005 2010 2015 2020],01,01),...
    'XTickLabel', {'1975', '1980', '1985', '1990', '1995', '2000', '2005', '2010', '2015', '2020'});
yline(0,'--');
% xline(datetime(2022,04,01),'-r')
title('24-Month Forecasted SR');

% Plot 90% confidence band of probability that SR<0
y1 = prctile(prob_SR24,5,1);
y2 = prctile(prob_SR24,95,1);
fill([xdate fliplr(xdate)], [y1 fliplr(y2)],LB, 'LineStyle','none');
hold on 
plot(xdate, mean(prob_SR24,1),'k');
legend({'90% CI','Mean'},'AutoUpdate','off','location','northwest')
set(gca, 'FontSize',20, 'XTick',datetime([1975 1980 1985 1990 1995 2000 2005 2010 2015 2020],01,01),...
    'XTickLabel', {'1975', '1980', '1985', '1990', '1995', '2000', '2005', '2010', '2015', '2020'});
yline(0,'--');
% xline(datetime(2022,04,01),'-r')
title('24-Month Forecast Probability');

save("./results/Our_SR_Prob.mat","prob_SR1","prob_SR3","prob_SR6","prob_SR12","prob_SR24");

%% Compare our model with Wu-Xia shadow rate model 
%  in forecasting ZLB porbability
load("./results/Our_SR_Prob.mat")
load("./results/WX_SR_Prob.mat");

% 1-month
plot(xdate, WX_prob_SR1,'LineStyle','--')
hold on
plot(xdate,mean(prob_SR1,1))
legend("WX model",'Our model',"Location","northwest")
set(gca, 'FontSize',20, 'XTick',datetime([1975 1980 1985 1990 1995 2000 2005 2010 2015 2020],03,01),...
    'XTickLabel', {'1975', '1980', '1985', '1990', '1995', '2000', '2005', '2010', '2015', '2020'});
title('1-Month Forecasted SR');

% 3-month
plot(xdate, WX_prob_SR3,'LineStyle','--')
hold on
plot(xdate,mean(prob_SR3,1))
legend("WX model",'Our model',"Location","northwest")
set(gca, 'FontSize',20, 'XTick',datetime([1975 1980 1985 1990 1995 2000 2005 2010 2015 2020],03,01),...
    'XTickLabel', {'1975', '1980', '1985', '1990', '1995', '2000', '2005', '2010', '2015', '2020'});
title('3-Month Forecasted SR');

% 6-month
plot(xdate, WX_prob_SR6,'LineStyle','--')
hold on
plot(xdate,mean(prob_SR6,1))
legend("WX model",'Our model',"Location","northwest")
set(gca, 'FontSize',20, 'XTick',datetime([1975 1980 1985 1990 1995 2000 2005 2010 2015 2020],03,01),...
    'XTickLabel', {'1975', '1980', '1985', '1990', '1995', '2000', '2005', '2010', '2015', '2020'});
title('6-Month Forecasted SR');

% 12-month
plot(xdate, WX_prob_SR12,'LineStyle','--')
hold on
plot(xdate,mean(prob_SR12,1))
legend("WX model",'Our model',"Location","northwest")
set(gca, 'FontSize',20, 'XTick',datetime([1975 1980 1985 1990 1995 2000 2005 2010 2015 2020],03,01),...
    'XTickLabel', {'1975', '1980', '1985', '1990', '1995', '2000', '2005', '2010', '2015', '2020'});
title('12-Month Forecasted SR');

% 24-month
plot(xdate, WX_prob_SR24,'LineStyle','--')
hold on
plot(xdate,mean(prob_SR24,1))
legend("WX model",'Our model',"Location","northwest")
set(gca, 'FontSize',20, 'XTick',datetime([1975 1980 1985 1990 1995 2000 2005 2010 2015 2020],03,01),...
    'XTickLabel', {'1975', '1980', '1985', '1990', '1995', '2000', '2005', '2010', '2015', '2020'});
title('24-Month Forecasted SR');


%% compare our ZLB duration with NY Fed Survey
clc
clear all

% probability of leaving ZLB within 1-year
Q0 = readtable("./data/fedsurvey/Q0.xlsx");
% plot(Q0.Date,1-Q0.within1Year)

% probability of returning to ZLB within 3-year
Q2 = readtable("./data/fedsurvey/Q2.xlsx");
% plot(Q2.Date, Q2.median,'k')
% hold on
% plot(Q2.Date, Q2.x25th,'r--')
% hold on
% plot(Q2.Date, Q2.x75th,'r--')

load("./results/Our_SR_Prob.mat")

xdate = datetime(1971,11,1) + calmonths(0:605);

plot(xdate,mean(prob_SR12,1))
hold on
plot(Q0.Date,1-Q0.within1Year,'--')
legend('Our model','NY Fed Survey',"Location","northwest")
set(gca, 'FontSize',20, 'XTick',datetime([1975 1980 1985 1990 1995 2000 2005 2010 2015 2020],03,01),...
    'XTickLabel', {'1975', '1980', '1985', '1990', '1995', '2000', '2005', '2010', '2015', '2020'});
title('12-Month Forecasted SR');

plot(xdate,mean(prob_SR24,1))
hold on
plot(Q2.Date, Q2.median,'--')
legend('Our model','NY Fed Survey',"Location","northwest")
set(gca, 'FontSize',20, 'XTick',datetime([1975 1980 1985 1990 1995 2000 2005 2010 2015 2020],03,01),...
    'XTickLabel', {'1975', '1980', '1985', '1990', '1995', '2000', '2005', '2010', '2015', '2020'});
title('24-Month Forecasted SR');








