(* ::Package:: *)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% This code does the KMS decomposition over all days in the sample
% Created for Hambur and Haque (2023) based on KMS

%% Very similair to  Code_for _HF _shock _through _ACM1



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 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');
day_data _start = "2001-04-1"; % exclude missing yearr
day_data _end = "2019-12-31"; % End at end of shockc

clc
options = optimoptions('fsolve','Display','off');

% paramS.NoF=4; % Number of fractors
irf_pc = 2;
yield_select = [1;12;24;36;60;120];

paramSH.cross.start = 105;
paramSH.cross.end = 166;

% yield_select = [2:120];

% 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

opts = spreadsheetImportOptions;
opts = detectImportOptions('\\ Input data\HF zcr\input\full_hist_day.xls');
opts = setvaropts(opts,{'Date'},'InputFormat','dd/MM/uuuu'); 
opts.Sheet='full_hist _day';
allday = readtable('\\ Input data\HF zcr\input\full_hist_day.xls',opts);
dates = allday(:,1);
allday = table2timetable(allday);
allday=allday{timerange(day_data _start,day_data _end),:};

DataYields.old= csvread('\\ Input data\Nom ATSM\Yields_m - EXTRA.csv');
data.Yields.old=DataYields.old(2:size(DataYields.old,1),3:size(DataYields.old,2));


%% Set up storage to turn to monthly data as only hae quarterly maturities
len = size( allday,1);

yields.ext = zeros(len,121);


%% Interpolate for missing months as only q freq 
yields.ext(:,1)=(allday(:,1))./100;


for i = 1:40
    yields.ext(:,1+i*3)=(allday(:,1+i))./100;
end

for i = 1:40
	gap = (yields.ext(:,i*3+1)-yields.ext(:,(i-1)*3+1))/3;
    yields.ext(:,(i-1)*3+1+1)=yields.ext(:,(i-1)*3+1)+gap;
    yields.ext(:,(i-1)*3+1+2)=yields.ext(:,(i-1)*3+1)+2*gap;

end
yields.ext = yields.ext(:,3:end);




% Standardise data matrix - need to make sure these data have shifted down
% similairy to old data
% old data over crossover period
Amean_old = mean(data.Yields.old); % mean from old data. Take this off so this data has negatvive mean
B=(yields.ext-repmat(Amean_old,[len,1]));
% B=B.*repmat(std(yields.ext(1:6189,:)),[len,1]);
M.pcs=(B*paramS.PC);


%% Run through ATSM

M.y.lvl=fit_y(M.pcs, paramS.An, paramS.Bn, len ,paramS.maxdate);
M.yrfexp.lvl=fit_y(M.pcs, paramS.Anrfexp, paramS.Bnrfexp, len ,paramS.maxdate);

M.y.pre=M.y.lvl(1:end-1,:);
M.yrfexp.pre=M.yrfexp.lvl(1:end-1,:);

M.y.post=M.y.lvl(2:end,:);
M.yrfexp.post=M.yrfexp.lvl(2:end,:);

M.y.change=M.y.post-M.y.pre;
M.yrfexp.change=M.yrfexp.post-M.yrfexp.pre;

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

%% Take PC and decompose
t=len-1; % lose 1 because in changes over days
%% Approach 3
pc_input = standardise(M.yrfexp.change(:,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);
% 

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].';
temp.a = M.tp.change(:,2:end)'*temp.z.'*inv(temp.z*temp.z.');
temp.residuals= M.tp.change(:,2:end)-(temp.a*temp.z).'; 

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


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;



%% reg of curve changes on PC to find shares expalined
clear temp
temp.z = [ones(t,1), SK.XR.exp.old(:,1)].';
temp.a = M.y.change(:,yield_select)'*temp.z.'*inv(temp.z*temp.z.');
temp.fit= (temp.a*temp.z).'; 
temp.residuals= M.y.change(:,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;

clear temp
temp.z = [ones(t,1), SK.XR.exp.old(:,2)].';
temp.a = M.y.change(:,yield_select)'*temp.z.'*inv(temp.z*temp.z.');
temp.fit= (temp.a*temp.z).'; 
temp.residuals= M.y.change(:,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;

clear temp
temp.z = [ones(t,1), SK.XR.tp.old(:,2)].';
temp.a = M.y.change(:,yield_select)'*temp.z.'*inv(temp.z*temp.z.');
temp.fit= (temp.a*temp.z).'; 
temp.residuals= M.y.change(:,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;

clear temp
temp.z = [ones(t,1), SK.XR.tp.old(:,1)].';
temp.a = M.y.change(:,yield_select)'*temp.z.'*inv(temp.z*temp.z.');
temp.fit= (temp.a*temp.z).'; 
temp.residuals= M.y.change(:,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;


clear temp
temp.z = [ones(t,1), SK.X.tp.old(:,1)].';
temp.a = M.y.change(:,yield_select)'*temp.z.'*inv(temp.z*temp.z.');
temp.fit= (temp.a*temp.z).'; 
temp.residuals= M.y.change(:,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;


contrs=zeros(t,5,6);
for i = 1:6
    contrs(:,:,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_all.csv')
writematrix(contrs(:,:,2), '\\ Outputs\HF decomp\contrs_1year_all.csv')
writematrix(contrs(:,:,3), '\\ Outputs\HF decomp\contrs_2year_all.csv')
writematrix(contrs(:,:,4), '\\ Outputs\HF decomp\contrs_3year_all.csv')
writematrix(contrs(:,:,5), '\\ Outputs\HF decomp\contrs_5year_all.csv')
writematrix(contrs(:,:,6), '\\ Outputs\HF decomp\contrs_10year_all.csv')



% write out the tp versus exp decomp

writematrix(M.yrfexp.change(:, yield_select),'\\ Outputs\HF decomp\Curve_all_day.xlsx', 'Sheet', 'yrf Change');

writematrix(M.tp.change(:, yield_select),'\\ Outputs\HF decomp\Curve_all_day.xlsx', 'Sheet', 'tp Change');


