% Exam the relationship among r*
clc
clear
close all


%%
% ADF test of stationarity for r*. All fail to reject the null of unit root
[h, pvalue,~,~,reg] = adftest(usrstar,'model','ARD','lags',0:12);  
[h, pvalue,~,~,reg] = adftest(ukrstar,'model','ARD','lags',0:12);  
[h, pvalue,~,~,reg] = adftest(gmrstar,'model','ARD','lags',0:12); 

% test the 1st order difference. All differenced series are stationary
adftest(diff(usrstar), 'model','ARD','lags',0:12);
adftest(diff(ukrstar),'model','ARD','lags',0:12);
adftest(diff(gmrstar),'model','ARD','lags',0:12);
clear h pvalue reg

% Johansen cointegration test
% rstars since 199101 - 202203
rstars = [usrstar(231:end-1) ukrstar(97:end) gmrstar(220:end)];
[h, pValue] = jcitest(rstars, 'model','H*','lags',24, 'alpha',0.05);
% the rank of is 2, so there is 1 common trend

% Estimate VECM with UK GM US
% newcolors = [0 0.5 1; 0.5 0 1; 0.7 0.7 0.7];
newcolors = [0.2 0.5 0.8;   
          0.4 0.7 1.0;   
          0.7 0.9 1.0];
Mdl = vecm(3,2,24);
% Mdl.Trend = 0;
% Mdl.CointegrationTrend = 0;
[EstMdl, EstSE] = estimate(Mdl,[rstars(:,2) rstars(:,3) rstars(:,1)], 'model','H*');
[Decomp,Lr, Ur] = fevd(EstMdl,'NumObs',120);
% plot the variance decomposition
x = 0:120;
subplot(2,2,1);
a = area(x,[1 0 0;Decomp(:,:,1)], 'FaceColor', 'flat');
colormap(newcolors);
ylim([0,1]);
legend('UK','DE','US','Location','southwest','FontSize',15);
title('UK','FontSize',20);
xlabel('period');
subplot(2,2,2);
area(x,[0 1 0;Decomp(:,:,2)],'FaceColor', 'flat');
colormap(newcolors);
ylim([0,1]);
title('DE','FontSize',20);
xlabel('period');
subplot(2,2,3);
area(x,[0.07 0 0.93;Decomp(:,:,3)],'FaceColor', 'flat');
colormap(newcolors);
ylim([0,1]);
title('US','FontSize',20);
xlabel('period');
% sgtitle('Variance Decomposition of r*', 'FontSize',20);

% Estimate VECM with GM UK US
Mdl = vecm(3,2,24);
% Mdl.Trend = 0;
% Mdl.CointegrationTrend = 0;
[EstMdl, EstSE] = estimate(Mdl,[rstars(:,3) rstars(:,2) rstars(:,1)], 'model','H*');
[Decomp,Lr, Ur] = fevd(EstMdl,'NumObs',120);
% plot the variance decomposition
x = 0:120;
subplot(2,2,1);
area(x,[1 0 0;Decomp(:,:,1)],'FaceColor', 'flat');
colormap(newcolors);
ylim([0,1]);
legend('DE','UK','US','Location','southwest','FontSize',15);
title('DE','FontSize',20);
xlabel('period');
subplot(2,2,2);
area(x,[0 1 0;Decomp(:,:,2)],'FaceColor', 'flat');
colormap(newcolors);
ylim([0,1]);
title('UK','FontSize',20);
xlabel('period');
subplot(2,2,3);
area(x,[0.07 0 0.93;Decomp(:,:,3)],'FaceColor', 'flat');
colormap(newcolors);
ylim([0,1]);
title('US','FontSize',20);
xlabel('period');
% sgtitle('Variance Decomposition of r*', 'FontSize',20);



%% Plot term premiums
% US 197110 - 202204
% UK 198301 - 202203
% GM 197210 - 202203
load './data/USrp.mat'
load './data/UKrp.mat'
load './data/GMrp.mat'

ustp = mean(USrp,1)';
uktp = mean(UKrp,1)';
gmtp = mean(GMrp,1)';

xdate1 = datetime(1971,11,1) + calmonths(0:605);
xdate2 = datetime(1983,1,1) + calmonths(0:470);
xdate3 = datetime(1972,10,1) + calmonths(0:593);

plot(xdate1, ustp);
hold on
plot(xdate2, uktp,'LineStyle','--');
hold on
plot(xdate3, gmtp,'LineStyle','-.');
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,'--');
legend('U.S.','U.K.','DE');


%% 
% Explore the relationship between term premium and inflation uncertainty
clc
clear
close all

% US 197110 - 202204
% UK 198301 - 202203
% GM 197210 - 202203
load './data/USrp.mat'
load './data/usinfl.mat'
load './data/UKrp.mat'
load './data/ukinfl.mat'
load './data/GMrp.mat'
load './data/gminfl.mat'

% 198301 - 202203
% usrp = mean(USrp(:,136:end-1),1)';
usrp = mean(USrp(:,135:end-1),1)';
ukrp = mean(UKrp,1)';
gmrp = mean(GMrp(:,124:end),1)';
usinfl = usinfl(136:end-1);
gminfl = gminfl(124:end);

xdate = datetime(1983,1,1) + calmonths(0:470);


% ADF test for term premiums
adftest(usrp)
adftest(ukrp)
adftest(gmrp)

% Principal component decomposition
rp_pc = Pca(usrp', ukrp', gmrp');
% 3 PCs are: 66.45, 26.25, 7.30
% R^2 are: 0.91, 0.36, 0.22

%%%%%======================================================================
%%%%%============================ US inflation ============================
%%%%%======================================================================
[sd_eps_us, sd_eta_us, tau_us] = UCSV2(usinfl,1100,1);

% sd of /epsilon and /eta
plot(mean(sd_eps_us,1))
hold on
plot(mean(sd_eta_us,1))

% inflation trend
plot(usinfl)
hold on 
plot(mean(tau_us,1)'+ mean(usinfl(1:12)) )


save('./results/US_ucsv.mat', 'sd_eps_us', 'sd_eta_us', 'tau_us' );


%%%%%======================================================================
%%%%%============================ UK inflation ============================
%%%%%======================================================================
[sd_eps_uk, sd_eta_uk, tau_uk] = UCSV2(ukinfl,1100,2);

% sd of /epsilon and /eta
plot(xdate,mean(sd_eps_uk,1)')
hold on
plot(xdate,mean(sd_eta_uk,1)')

% inflation trend
plot(xdate,ukinfl)
hold on
plot(xdate,mean(tau_uk,1)'+ mean(ukinfl(1:12)) )

save('./results/UK_ucsv.mat', 'sd_eps_uk', 'sd_eta_uk', 'tau_uk' );


%%%%%======================================================================
%%%%%============================ GM inflation ============================
%%%%%======================================================================
[sd_eps_gm, sd_eta_gm, tau_gm] = UCSV2(gminfl,5100,3);

% xdate = datetime(1972,10,1) + calmonths(0:593);
% sd of /epsilon and /eta
plot(xdate,mean(sd_eps_gm,1)')
hold on
plot(xdate, mean(sd_eta_gm,1)')

% inflation trend
plot(xdate,gminfl)
hold on
plot(xdate,mean(tau_gm,1)'+ mean(gminfl(1:12)) )


save('./results/GM_ucsv.mat', 'sd_eps_gm', 'sd_eta_gm', 'tau_gm' );

[~,~,~,Response,Lower,Upper] = gct([myx gmrp],12,["inflation uncertainty";"term premium"]);
irfinf2tp = Response(:,2,1);    % IRF of inflation uncertainty to a shock to term premium
IRFinf2tp = [Lower(:,2,1) Upper(:,2,1)];
xdate4 = 0:59;
h2 = fill([xdate4 fliplr(xdate4)], [IRFinf2tp(:,1)' fliplr(IRFinf2tp(:,2)')],'blue');
h2.FaceColor = [0 0.8 1];
h2.EdgeColor = 'none';
hold on
h1 = plot(xdate4, irfinf2tp,'k');
xlabel("period");
ylabel("response");
title('GM IR of inflation uncertainty to term premium')



%%%%%======================================================================
%%%%%========================= Panel Regression ===========================
%%%%%======================================================================

% S.D. of permanent shock
infunct = [mean(sd_eta_us,1)', mean(sd_eta_uk,1)', mean(sd_eta_gm,1)'];
[n,~] = size(infunct);
rps = [usrp, ukrp, gmrp];
myx = cell(n,1);
for i = 1:n
    myx{i} = [eye(3), infunct(i,:)'];
end
[b,sig,E,CovB] = mvregress(myx,rps, 'algorithm','mvn');
b(4)/sqrt(CovB(4,4))  % b(4) = 1.20
tboots = bbs3(rps, infunct, 50000, 6);
prctile(tboots,96.6)  % 0.03

% S.D. of transitory shock
infunct2 = [mean(sd_eps_us,1)', mean(sd_eps_uk,1)', mean(sd_eps_gm,1)'];
[n,~] = size(infunct2);
rps = [usrp, ukrp, gmrp];
myx = cell(n,1);
for i = 1:n
    myx{i} = [eye(3), infunct2(i,:)'];
end
[b2,sig2,E2,CovB2] = mvregress(myx,rps, 'algorithm','mvn');
b2(4)/sqrt(CovB2(4,4))  % b(4) = -0.6716
tboots2 = bbs3(rps, infunct2, 50000, 6);
prctile(tboots2,0.8)  % 

% cycle s.d. and trend s.d.
myx = cell(n,1);
for i = 1:n
    myx{i} = [eye(3), infunct2(i,:)' infunct(i,:)'];
end
[b3,sig3,E3,CovB3] = mvregress(myx,rps, 'algorithm','mvn');
b3(4)/sqrt(CovB3(4,4))  % b(4) = -0.7341
b3(5)/sqrt(CovB3(5,5))  % b(5) = 1.3993
tb1 = bbs4(rps, infunct2, infunct, 50000, 6, 1);
tb2 = bbs4(rps, infunct2, infunct, 50000, 6, 2);
prctile(tb1,0.4)  % 0.007
prctile(tb2,98.2) % 0.01


%% Plot sd of permanent shock to inflation
load('./results/US_ucsv.mat')
load('./results/UK_ucsv.mat')
load('./results/GM_ucsv.mat')

xdate = datetime(1983,1,1) + calmonths(0:470);

subplot(2,2,1)
plot(xdate,mean(sd_eta_us,1));
title('US');
subplot(2,2,2)
plot(xdate,mean(sd_eta_uk,1));
title('UK');
subplot(2,2,3)
plot(xdate,mean(sd_eta_gm,1));
title('DE');


%% Wright's survey data 199001-200906
Surveydata = readtable("./data/Surveydata.xls");
usinfunc_suv = Surveydata{:,'USA'};
ukinfunc_suv = Surveydata{:,'UK'};
gminfunc_suv = Surveydata{:,'Germany'};

infunct3 = [usinfunc_suv, ukinfunc_suv, gminfunc_suv];  % 198910-200906 
[n,d] = size(infunct3);
% tps = [mean(USrp(:,217:453),1)', mean(UKrp(:,82:318),1)', mean(GMrp(:,205:441),1)'];
tps = [mean(USrp(:,216:452),1)', mean(UKrp(:,82:318),1)', mean(GMrp(:,205:441),1)'];
myx = cell(n,1);
for i = 1:n
    myx{i} = [eye(3), infunct3(i,:)'];
end
[b_w,sig_w,E_w,CovB_w] = mvregress(myx,tps, 'algorithm','mvn');
b_w(4)/sqrt(CovB_w(4,4))  % b =  -0.5440
tb_w = bbs3(tps, infunct3, 50000, 6);
prctile(tb_w, 3)  % 3%



%% Compare trends of shadow rates
clear all
load("./data/taus.mat");

xdate1 = datetime(1971,11,1) + calmonths(0:605);
xdate2 = datetime(1983,1,1) + calmonths(0:470);
xdate3 = datetime(1972,10,1) + calmonths(0:593);

plot(xdate1, ustau,'LineWidth',1);
hold on
plot(xdate2, uktau,'LineWidth',1,'LineStyle','--');
hold on
plot(xdate3, gmtau,'LineWidth',1,'LineStyle','-.');
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,'--');
legend('U.S.','U.K.','DE');


