(* ::Package:: *)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% This code decomposes MP shocks into exp rates and TP, and then into the factors outlines in Kaminska et al 2022
% Created for Hambur and Haque (2023) based on KMS
% Baseline estiamtes code with other events



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Options
clear all
close all
restoredefaultpath
addpath('Code\Nom ATSM') % Contains some of the sub-functions
addpath('\Code\Decomposing shocks\Subs') % Contains some of the sub-functions
addpath('Code\Decomposing shocks\additional');

clc
options = optimoptions('fsolve','Display','off');
month_data _start = "1994-01-31"; % exclude missing yearr
% setting up some infomration on crossover dates of the ATSM and the yield data
% paramS.NoF=4; % Number of fractors
paramSH.cross.start = 105;
paramSH.cross.end = 398;

yield_select = [1;12;24;36;60;120]; % selection of yields to include in th PC creation
% yield_select = [2:120]; can use this, but makes very little diffrence

% Select the model we want - number of parameters, survey or bootstrap etc
load('Input data\Nom ATSM\paramM93_S_ 4fac_cash _both _highsd _ 2022.mat');
load('Input data\Nom ATSM\M93_S_ 4fac_cash _both _highsd _ 2022.mat');
MS=M;

clear M;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Start

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Load and format data



% Load data
DataYields.presp= csvread('\Input data\HF zcr\output\2006-04-28.To .2020-02-05.Yield.pre_extra.csv'); % other events
DataYields.postsp= csvread('\Input data\HF zcr\output\2006-04-28.To .2020-02-05.Yield.post_extra.csv'); % other events
DataYields.premp= csvread('G\Input data\HF zcr\output\2001-04-04.To .2019-12-03.Yield.pre.csv');
DataYields.postmp= csvread('\Input data\HF zcr\output\2001-04-04.To .2019-12-03.Yield.post.csv');
DataYields.old= csvread('\Note replication\Input data\Nom ATSM\Yields_m - EXTRA.csv'); % full sample daily data for cross-ver sample comparison
% Make data matrices
data.Yields.premp=DataYields.premp(2:size(DataYields.premp,1),5:size(DataYields.premp,2));
data.Yields.premp=data.Yields.premp/100; % old data scaled down by 100
data.Yields.postmp=DataYields.postmp(2:size(DataYields.postmp,1),5:size(DataYields.postmp,2));
data.Yields.postmp=data.Yields.postmp/100;
data.Yields.presp=DataYields.presp(2:size(DataYields.presp,1),5:size(DataYields.presp,2));
data.Yields.presp=data.Yields.presp/100; % old data scaled down by 100
data.Yields.postsp=DataYields.postsp(2:size(DataYields.postsp,1),5:size(DataYields.postsp,2));
data.Yields.postsp=data.Yields.postsp/100;
data.Yields.pre = [data.Yields.premp; data.Yields.presp];
data.Yields.post = [data.Yields.postmp; data.Yields.postsp];

data.Yields.old=DataYields.old(2:size(DataYields.old,1),3:size(DataYields.old,2));
data.dates = [DataYields.premp(2:size(DataYields.premp,1),1); DataYields.presp(2:size(DataYields.presp,1),1)];
data.Yields.change = data.Yields.pre-data.Yields.post;
[t, q] = size(data.Yields.pre);

%% Now get the PCs fro mthe yeilds so can put them into the ATSM

% I take a differnt approach. They take diffrences then take the PC. But issue is loading this way may differ to levels
% I apply the loading from the ATSM to yeilds before and after. Fit ATSM
% before and after, and then look at changes
% Requires me to do some re-scaling of the yields/ factors given differign
% samples, so when standardise for PC creates issues.

%Standardise data matrix



Amean_old = mean(data.Yields.old); % mean from old data. Take this off so this data has negatvive mean
B=(data.Yields.pre-repmat(Amean_old,[t,1]));
M.X.pre.old = B*paramS.PC;

B=(data.Yields.post-repmat(Amean_old,[t,1]));
% B=B.*repmat(std(yields.ext(1:6189,:)),[len,1]);
M.X.post.old = B*paramS.PC;


% might need change in bit above
%% Now run through ATSM
% Fitted yields
M.y.pre.old=fit_y(M.X.pre.old, paramS.An, paramS.Bn, t ,paramS.maxdate);
M.y.post.old=fit_y(M.X.post.old, paramS.An, paramS.Bn, t ,paramS.maxdate);
M.y.change.old = M.y.pre.old-M.y.post.old;

%% Risk free
M.yrf.pre.old=fit_y(M.X.pre.old, paramS.Anrf, paramS.Bnrf, t ,paramS.maxdate);
M.yrf.post.old=fit_y(M.X.post.old, paramS.Anrf, paramS.Bnrf, t ,paramS.maxdate);
M.yrf.change.old = M.yrf.pre.old-M.yrf.post.old;

%% Risk free  AVERAGE WITH NO CONVEXITY ADJ (means rate is average of short rates - preffered metric)
M.yrfexp.pre.old=fit_y(M.X.pre.old, paramS.Anrfexp, paramS.Bnrfexp, t ,paramS.maxdate);
M.yrfexp.post.old=fit_y(M.X.post.old, paramS.Anrfexp, paramS.Bnrfexp, t ,paramS.maxdate);
M.yrfexp.change.old = M.yrfexp.pre.old-M.yrfexp.post.old;
% Calckaute forward exp RF short rates

all_per=[1:paramS.maxdate];
M.yrffwd.pre.old=fit_short(M.X.pre.old,paramS.delta0,paramS.delta1,paramS.MuMLE,paramS.BMLE-eye(paramS.NoF),all_per,paramS.NoF,t);
M.yrffwd.post.old=fit_short(M.X.post.old,paramS.delta0,paramS.delta1,paramS.MuMLE,paramS.BMLE-eye(paramS.NoF),all_per,paramS.NoF,t);
M.yrffwd.change.old = M.yrffwd.pre.old-M.yrffwd.post.old;

% Calckaute forward exp  short rates

M.yfwd.pre.old=fit_short(M.X.pre.old,paramS.delta0,paramS.delta1,(paramS.MuMLE-paramS.L0),(paramS.BMLE-paramS.L1)-eye(paramS.NoF),all_per,paramS.NoF,t);
M.yfwd.post.old=fit_short(M.X.post.old,paramS.delta0,paramS.delta1,(paramS.MuMLE-paramS.L0),(paramS.BMLE-paramS.L1)-eye(paramS.NoF),all_per,paramS.NoF,t);
M.yfwd.change.old=M.yfwd.pre.old-M.yfwd.post.old;

% Average TP
M.tp.pre.old=M.y.pre.old-M.yrfexp.pre.old;
M.tp.post.old=M.y.post.old-M.yrfexp.post.old;
M.tp.change.old=M.tp.pre.old-M.tp.post.old;


% period TP 
M.tpfwd.pre.old=M.yfwd.pre.old-M.yrffwd.pre.old;
M.tpfwd.post.old=M.yfwd.post.old-M.yrffwd.post.old;
M.tpfwd.change.old=M.tpfwd.pre.old-M.tpfwd.post.old;


%% Make the PCs of the exp and TP hcanges
pc_input = standardise(M.yrfexp.change.old(:,yield_select));
%[PC, Mean,V,X]= pca1(test);
[pc2,W2]=extract(pc_input,2);
SK.X.exp.old=pc2;
 SK.PC.exp.old=W2';
SK.V.exp.old=mean(W2.^2,2);



%  clear temp
% Rotate the PC so one loads on shortrate, and one doesn't
chck=-1;
while chck<0
x0=vec(getqr(randn(2,2)));
f = @(x)rotate_pc(x,W2(:,1));
[out,fval,exitflag]=fsolve(f,x0,options);
[F,U]=rotate_pc(out,W2(:,1));
pc0=pc2*U;
w20=pc0\pc_input;
if w20(1,1)>0 && w20(2,end)>0
    chck=1;
end
end

SK.PCR.exp.old=w20';
SK.XR.exp.old=pc0;


% Now do regression of TP on expected
temp.z = [ones(t,1), SK.XR.exp.old].';
% Nopte droppping SR TP as always zero by construction in ATSM (no term)
temp.a = M.tp.change.old(:,2:end)'*temp.z.'*inv(temp.z*temp.z.');
temp.residuals= M.tp.change.old(:,2:end)-(temp.a*temp.z).'; 
%% and take the STDEV


pc_input = standardise(temp.residuals);
%[PC, Mean,V,X]= pca1(test);
[pc3,W3]=extract(pc_input,2);
SK.X.tp.old=pc3;
 SK.PC.tp.old=W3';
SK.V.tp.old=mean(W3.^2,2);
% TP

%% and rotate
chck=-1;
while chck<0
x0=vec(getqr(randn(2,2)));
f = @(x)rotate_pc(x,W3(:,1));
[out,fval,exitflag]=fsolve(f,x0,options);
[F,U]=rotate_pc(out,W3(:,1));
tp0=pc3*U;
w30=tp0\pc_input;
if w30(1,1)>0 && w30(2,end)>0
    chck=1;
end
end

SK.PCR.tp.old=w30';
SK.XR.tp.old=tp0;


%% Now outputting
% note my TP are switched, so both 1 tp are the ones loaded on current
% XR are the rotted factors, X are non rotated
term_str _pc.old = exp_f(SK.XR.exp.old, SK.XR.tp.old, SK.X.exp.old, SK.X.tp.old, data.dates);

writetable(term_str _pc.old, 'Input data\VAR data\pc_yield__allevent.csv')

%% Now estimate shars of yields accounted for by each of the factors via regressions
% choice of yields to look at (all shifted by one)
data.Yields.select.change = data.Yields.change(:, [1;11;23;35;59;119]);

% No reg, so just gett SS variables
clear temp
temp.z = [ones(t,1)].';
temp.a = M.y.change.old(:,yield_select)'*temp.z.'*inv(temp.z*temp.z.');
temp.fit= (temp.a*temp.z).'; 
temp.residuals= M.y.change.old(:,yield_select)-(temp.a*temp.z).'; 
temp.residuals2 = temp.residuals.^2;
temp.SSE=sum(sum(temp.residuals2,1),2);

decomp.int.old.beta=temp.a;
decomp.int.old.fit=temp.fit;
decomp.int.old.sse=temp.SSE;

% Now action factor
clear temp
temp.z = [ones(t,1), SK.XR.exp.old(:,1)].';
temp.a = M.y.change.old(:,yield_select)'*temp.z.'*inv(temp.z*temp.z.');
temp.fit= (temp.a*temp.z).'; 
temp.residuals= M.y.change.old(:,yield_select)-(temp.a*temp.z).'; 
temp.residuals2 = temp.residuals.^2;
temp.SSE=sum(sum(temp.residuals2,1),2);

decomp.exppc1.old.beta=temp.a;
decomp.exppc1.old.fit=temp.fit;
decomp.exppc1.old.sse=temp.SSE;

% Now path factor
clear temp
temp.z = [ones(t,1), SK.XR.exp.old(:,2)].';
temp.a = M.y.change.old(:,yield_select)'*temp.z.'*inv(temp.z*temp.z.');
temp.fit= (temp.a*temp.z).'; 
temp.residuals= M.y.change.old(:,yield_select)-(temp.a*temp.z).'; 
temp.residuals2 = temp.residuals.^2;
temp.SSE=sum(sum(temp.residuals2,1),2);

decomp.exppc2.old.beta=temp.a;
decomp.exppc2.old.fit=temp.fit;
decomp.exppc2.old.sse=temp.SSE;

% Now action factor first rotated TP
clear temp
temp.z = [ones(t,1), SK.XR.tp.old(:,2)].';
temp.a = M.y.change.old(:,yield_select)'*temp.z.'*inv(temp.z*temp.z.');
temp.fit= (temp.a*temp.z).'; 
temp.residuals= M.y.change.old(:,yield_select)-(temp.a*temp.z).'; 
temp.residuals2 = temp.residuals.^2;
temp.SSE=sum(sum(temp.residuals2,1),2);

decomp.tppc1.old.beta=temp.a;
decomp.tppc1.old.fit=temp.fit;
decomp.tppc1.old.sse=temp.SSE;

% Second rotated TP
clear temp
temp.z = [ones(t,1), SK.XR.tp.old(:,1)].';
temp.a = M.y.change.old(:,yield_select)'*temp.z.'*inv(temp.z*temp.z.');
temp.fit= (temp.a*temp.z).'; 
temp.residuals= M.y.change.old(:,yield_select)-(temp.a*temp.z).'; 
temp.residuals2 = temp.residuals.^2;
temp.SSE=sum(sum(temp.residuals2,1),2);

decomp.tppc2.old.beta=temp.a;
decomp.tppc2.old.fit=temp.fit;
decomp.tppc2.old.sse=temp.SSE;

% Un roated first TP factor
clear temp
temp.z = [ones(t,1), SK.X.tp.old(:,1)].';
temp.a = M.y.change.old(:,yield_select)'*temp.z.'*inv(temp.z*temp.z.');
temp.fit= (temp.a*temp.z).'; 
temp.residuals= M.y.change.old(:,yield_select)-(temp.a*temp.z).'; 
temp.residuals2 = temp.residuals.^2;
temp.SSE=sum(sum(temp.residuals2,1),2);

decomp.tpdir.old.beta=temp.a;
decomp.tpdir.old.fit=temp.fit;
decomp.tpdir.old.sse=temp.SSE;

% Everything for SSR and SSE
clear temp
temp.z = [ones(t,1), SK.XR.exp.old, SK.XR.tp.old].';
temp.a = M.y.change.old(:,yield_select)'*temp.z.'*inv(temp.z*temp.z.');
temp.fit= (temp.a*temp.z).'; 
temp.residuals= M.y.change.old(:,yield_select)-(temp.a*temp.z).'; 
temp.residuals2 = temp.residuals.^2;
temp.SSE=sum(sum(temp.residuals2,1),2);

decomp.all.old.beta=temp.a;
decomp.all.old.fit=temp.fit;
decomp.all.old.sse=temp.SSE;

var = sum(sum(M.y.change.old(:,yield_select).^2,1),2);
% R squared
r2=zeros(6,1);
r2(1,1) = 1-decomp.int.old.sse/var;
r2(2,1) = 1-decomp.exppc1.old.sse/var;
r2(3,1)= 1-decomp.exppc2.old.sse/var;
r2(4,1)= 1-decomp.tppc2.old.sse/var;
r2(5,1)= 1-decomp.tpdir.old.sse/var;
r2(6,1)= 1-decomp.all.old.sse/var;


% output
contrs=zeros(t,6,6);
for i = 1:6
    contrs(:,:,i)=[data.Yields.select.change(:,i), decomp.exppc1.old.fit(:,i), decomp.exppc2.old.fit(:,i), decomp.tppc1.old.fit(:,i), decomp.tppc2.old.fit(:,i) decomp.tpdir.old.fit(:,i) ];
end

writematrix(contrs(:,:,1), '\Outputs\HF decomp\contrs_1m_adj _extra _stack.csv')
writematrix(contrs(:,:,2), '\Outputs\HF decomp\contrs_1year_adj _extra _stack.csv')
writematrix(contrs(:,:,3), '\Outputs\HF decomp\contrs_2year_adj _extra _stack.csv')
writematrix(contrs(:,:,4), '\Outputs\HF decomp\contrs_3year_adj _extra _stack.csv')
writematrix(contrs(:,:,5), '\Outputs\HF decomp\contrs_5year_adj _extra _stack.csv')
writematrix(contrs(:,:,6), '\Outputs\HF decomp\contrs_10year_adj _extra _stack.csv')


%% Now just seperate changes into the TP and the exp rates for cahrtin

% Yield
data.Yields.change(:, [1;11;23;35;59;119]);


writematrix(M.y.pre.old(:, yield_select),'\Outputs\HF decomp\Curve_extra_stack.xlsx', 'Sheet', 'Pre');
writematrix(M.y.post.old(:, yield_select),'\Outputs\HF decomp\Curve_extra_stack.xlsx', 'Sheet', 'Post');
writematrix(M.y.change.old(:, yield_select),'\Outputs\HF decomp\Curve_extra_stack.xlsx', 'Sheet', 'Change');

writematrix(M.yrfexp.pre.old(:, yield_select),'\Outputs\HF decomp\Curve_extra_stack.xlsx', 'Sheet', 'yrf Pre');
writematrix(M.yrfexp.post.old(:, yield_select),'\Outputs\HF decomp\Curve_extra_stack.xlsx', 'Sheet', 'yrf Post');
writematrix(M.yrfexp.change.old(:, yield_select),'\Outputs\HF decomp\Curve_extra_stack.xlsx', 'Sheet', 'yrf Change');

writematrix(M.tp.pre.old(:, yield_select),'\Outputs\HF decomp\Curve_extra_stack.xlsx', 'Sheet', 'tp Pre');
writematrix(M.tp.post.old(:, yield_select),'\Outputs\HF decomp\Curve_extra_stack.xlsx', 'Sheet', 'tp Post');
writematrix(M.tp.change.old(:, yield_select),'\Outputs\HF decomp\Curve_extra_stack.xlsx', 'Sheet', 'tp Change');


%% Now estimate shars of yields accounted for by each of the factors via regressions - just for non-MP dates
% choice of yields to look at (all shifted by one)
data.Yields.select.change = data.Yields.change(:, [1;11;23;35;59;119]);
t=size(M.y.change.old(200:end,1),1);
% No reg, so just gett SS variables
clear temp
temp.z = [ones(t,1)].';
temp.a = M.y.change.old(200:end,yield_select)'*temp.z.'*inv(temp.z*temp.z.');
temp.fit= (temp.a*temp.z).'; 
temp.residuals= M.y.change.old(200:end,yield_select)-(temp.a*temp.z).'; 
temp.residuals2 = temp.residuals.^2;
temp.SSE=sum(sum(temp.residuals2,1),2);

decomp.int.oldev.beta=temp.a;
decomp.int.oldev.fit=temp.fit;
decomp.int.oldev.sse=temp.SSE;

% Now action factor
clear temp
temp.z = [ones(t,1), SK.XR.exp.old(200:end,1)].';
temp.a = M.y.change.old(200:end,yield_select)'*temp.z.'*inv(temp.z*temp.z.');
temp.fit= (temp.a*temp.z).'; 
temp.residuals= M.y.change.old(200:end,yield_select)-(temp.a*temp.z).'; 
temp.residuals2 = temp.residuals.^2;
temp.SSE=sum(sum(temp.residuals2,1),2);

decomp.exppc1.oldev.beta=temp.a;
decomp.exppc1.oldev.fit=temp.fit;
decomp.exppc1.oldev.sse=temp.SSE;

% Now path factor
clear temp
temp.z = [ones(t,1), SK.XR.exp.old(200:end,2)].';
temp.a = M.y.change.old(200:end,yield_select)'*temp.z.'*inv(temp.z*temp.z.');
temp.fit= (temp.a*temp.z).'; 
temp.residuals= M.y.change.old(200:end,yield_select)-(temp.a*temp.z).'; 
temp.residuals2 = temp.residuals.^2;
temp.SSE=sum(sum(temp.residuals2,1),2);

decomp.exppc2.oldev.beta=temp.a;
decomp.exppc2.oldev.fit=temp.fit;
decomp.exppc2.oldev.sse=temp.SSE;

% Now action factor first rotated TP
clear temp
temp.z = [ones(t,1), SK.XR.tp.old(200:end,2)].';
temp.a = M.y.change.old(200:end,yield_select)'*temp.z.'*inv(temp.z*temp.z.');
temp.fit= (temp.a*temp.z).'; 
temp.residuals= M.y.change.old(200:end,yield_select)-(temp.a*temp.z).'; 
temp.residuals2 = temp.residuals.^2;
temp.SSE=sum(sum(temp.residuals2,1),2);

decomp.tppc1.oldev.beta=temp.a;
decomp.tppc1.oldev.fit=temp.fit;
decomp.tppc1.oldev.sse=temp.SSE;

% Second rotated TP
clear temp
temp.z = [ones(t,1), SK.XR.tp.old(200:end,1)].';
temp.a = M.y.change.old(200:end,yield_select)'*temp.z.'*inv(temp.z*temp.z.');
temp.fit= (temp.a*temp.z).'; 
temp.residuals= M.y.change.old(200:end,yield_select)-(temp.a*temp.z).'; 
temp.residuals2 = temp.residuals.^2;
temp.SSE=sum(sum(temp.residuals2,1),2);

decomp.tppc2.oldev.beta=temp.a;
decomp.tppc2.oldev.fit=temp.fit;
decomp.tppc2.oldev.sse=temp.SSE;

% Un roated first TP factor
clear temp
temp.z = [ones(t,1), SK.X.tp.old(200:end,1)].';
temp.a = M.y.change.old(200:end,yield_select)'*temp.z.'*inv(temp.z*temp.z.');
temp.fit= (temp.a*temp.z).'; 
temp.residuals= M.y.change.old(200:end,yield_select)-(temp.a*temp.z).'; 
temp.residuals2 = temp.residuals.^2;
temp.SSE=sum(sum(temp.residuals2,1),2);

decomp.tpdir.oldev.beta=temp.a;
decomp.tpdir.oldev.fit=temp.fit;
decomp.tpdir.oldev.sse=temp.SSE;

% Everything for SSR and SSE
clear temp
temp.z = [ones(t,1), SK.XR.exp.old(200:end,:), SK.XR.tp.old(200:end,:)].';
temp.a = M.y.change.old(200:end,yield_select)'*temp.z.'*inv(temp.z*temp.z.');
temp.fit= (temp.a*temp.z).'; 
temp.residuals= M.y.change.old(200:end,yield_select)-(temp.a*temp.z).'; 
temp.residuals2 = temp.residuals.^2;
temp.SSE=sum(sum(temp.residuals2,1),2);

decomp.all.oldev.beta=temp.a;
decomp.all.oldev.fit=temp.fit;
decomp.all.oldev.sse=temp.SSE;

var = sum(sum(M.y.change.old(200:end,yield_select).^2,1),2);
% R squared
r2ev=zeros(6,1);
r2ev(1,1) = 1-decomp.int.oldev.sse/var;
r2ev(2,1) = 1-decomp.exppc1.oldev.sse/var;
r2ev(3,1)= 1-decomp.exppc2.oldev.sse/var;
r2ev(4,1)= 1-decomp.tppc2.oldev.sse/var;
r2ev(5,1)= 1-decomp.tpdir.oldev.sse/var;
r2ev(6,1)= 1-decomp.all.oldev.sse/var;


