% Out-of-Sample Forecast of 10-year Forward Rate
% Sample: 197110 - 202204
% maturities: 3m,6m,1-10y
% For BR model, we generate yields from GSW, and use end of quarter data
% and run the model in R code provided in BR's paper (1971Q4 - 2018Q1) 
clear all
close all
clc

% % Generate yields from GSW 
% mats = [3, 6, 12:12:120]';
% Yn = genYn('US_sv.csv',mats);
% 
% % Generate monthly Forward Rate 197110-202204
% load USdata.mat
% FR = USdata{:,3:14};

% Yields and Forward rate from LW: 197111-202204
% mats = [3, 6, 12:12:120]';
usyields = readtable("US_yields_LW.csv");
Yn = usyields{126:731,2:13};

% % Generate monthly Forward Rate 197111-202204
load USdata_LW2.mat
FR = USdata{:,3:14};

% % Generate end-of-quarter yields (1971Q4 - 2022Q1)
% mats2 = [3, 6, 12:12:180]';
% Yn2 = genYn('US_sv.csv',mats2);
% Yn_Quart = Yn2(3:3:end-1,:);   % quarter end
% % Yn_Quart = Yn2(1:3:end-3,:);   % quarter start
% % Yn_Quart = Yn2(3:3:end-49,:);
% Yn_Quart = nan((size(Yn2,1)-1)/3, size(Yn2,2));
% for i = 1:size(Yn2,1)-1
%     if mod(i,3)==0
%         Yn_Quart(i/3,:) = mean(Yn2(i-2:i,:));
%     end
% end
% writematrix(Yn_Quart, 'Y_GSW3.csv');
% writematrix(Yn_Quart, 'Y_GSW_QuarterStart.csv');


% Generate averaged quarter yields (1971Q4 - 2022Q1)
Yn2 = usyields{126:731,2:18};
Yn_Quart = nan(size(Yn2));
xdate = datetime(1971,11,1) + calmonths(0:605)';
xquarts = datetime(1971,12,1) + calquarters(0:201)';
Yn_Quart(1,:) = mean(Yn2(1:2,:),1);
for i = 5:3:(size(Yn2,1)-1)
    Yn_Quart(i,:) = mean(Yn2(i-2:i,:),1);
end
Yn_Quart(any(isnan(Yn_Quart),2),:) = []; 
writematrix(Yn_Quart, './results/Y_Quart_LW.csv');


% We compare our model (SDTSM) with 1.Random Walk (RW) 
%                                   2.Yield-only model (JSZ)
%                                   3. Bauer and Rudebusch ese model(BR)
% For SDTSM and BR, we accept pesudo out of sample forecast
% For JSZ and RW, the estimations are rolling forward

mcmc_ind = 50000;
dt = 1;
rlb = 0;
N =3;
mats = [3, 6, 12:12:120]';

% load('mcmc_US.mat', 'Omega','Phi','kinfQ','lamQ','sige2','gamma','Xbar','sigtau2');
load('temp100000.mat', 'Omega','Phi','kinfQ','lamQ','sige2','gamma','Xbar','sigtau2');

temp_omega = mean(Omega(mcmc_ind+1,:),1);
temp_omega = [temp_omega(1,1) temp_omega(1,2) temp_omega(1,3);...
    temp_omega(1,2) temp_omega(1,4) temp_omega(1,5);...
    temp_omega(1,3) temp_omega(1,5) temp_omega(1,6)];
temp_phi = reshape(mean(Phi(mcmc_ind+1:end,:),1),3,3);
temp_loads = SRTSM_loadings(mean(kinfQ(mcmc_ind+1:end)), diag(mean(lamQ(mcmc_ind+1:end,:),1)), temp_omega, mats, dt);
temp_sige2 =  mean(sige2(mcmc_ind+1:end));
temp_gamma = mean(gamma(mcmc_ind+1:end,:),1)';
temp_Xbar =  mean(Xbar(mcmc_ind+1:end,:),1)';
temp_sigtau2 = mean(sigtau2(mcmc_ind+1:end));


% random walk forecast
rw3 = forecast_rw(FR(:,12),3);
rw6 = forecast_rw(FR(:,12),6);
rw10 = forecast_rw(FR(:,12),10);
rw20 = forecast_rw(FR(:,12),20);
rw30 = forecast_rw(FR(:,12),30);
rw40 = forecast_rw(FR(:,12),40);

% writematrix(rw3.resid, './results/rw3.csv');
writematrix(rw6.resid, './results/rw6.csv');
writematrix(rw10.resid, './results/rw10.csv');
writematrix(rw20.resid, './results/rw20.csv');
writematrix(rw30.resid, './results/rw30.csv');
writematrix(rw40.resid, './results/rw40.csv');

% jsz forecast
% jsz3 = forecast_jsz3(FR(:,12),Yn, mats', 3);
jsz6 = forecast_jsz3(FR(:,12),Yn, mats', 6);
jsz10 = forecast_jsz3(FR(:,12),Yn, mats', 10);
jsz20 = forecast_jsz3(FR(:,12),Yn, mats', 20);
jsz30 = forecast_jsz3(FR(:,12),Yn, mats', 30);
jsz40 = forecast_jsz3(FR(:,12),Yn, mats', 40);

writematrix(jsz6.resid, './results/jsz6.csv');
writematrix(jsz10.resid, './results/jsz10.csv');
writematrix(jsz20.resid, './results/jsz20.csv');
writematrix(jsz30.resid, './results/jsz30.csv');
writematrix(jsz40.resid, './results/jsz40.csv');

% SDTSM forecast
sdtsm3 = forecast_sdtsm(FR, temp_phi, temp_omega, temp_sige2, temp_gamma,...
    temp_Xbar, temp_sigtau2, 12, temp_loads.an, temp_loads.bn, temp_loads.SigmaQ, rlb, 3);

sdtsm6 = forecast_sdtsm(FR, temp_phi, temp_omega, temp_sige2, temp_gamma,...
    temp_Xbar, temp_sigtau2, 12, temp_loads.an, temp_loads.bn, temp_loads.SigmaQ, rlb, 6);

sdtsm10 = forecast_sdtsm(FR, temp_phi, temp_omega, temp_sige2, temp_gamma,...
    temp_Xbar, temp_sigtau2, 12, temp_loads.an, temp_loads.bn, temp_loads.SigmaQ, rlb, 10);

sdtsm20 = forecast_sdtsm(FR, temp_phi, temp_omega, temp_sige2, temp_gamma,...
    temp_Xbar, temp_sigtau2, 12, temp_loads.an, temp_loads.bn, temp_loads.SigmaQ, rlb, 20);

sdtsm30 = forecast_sdtsm(FR, temp_phi, temp_omega, temp_sige2, temp_gamma,...
    temp_Xbar, temp_sigtau2, 12, temp_loads.an, temp_loads.bn, temp_loads.SigmaQ, rlb, 30);

sdtsm40 = forecast_sdtsm(FR, temp_phi, temp_omega, temp_sige2, temp_gamma,...
    temp_Xbar, temp_sigtau2, 12, temp_loads.an, temp_loads.bn, temp_loads.SigmaQ, rlb, 40);

temp = forecast_sdtsm(FR, temp_phi, temp_omega, temp_sige2, temp_gamma,...
    temp_Xbar, temp_sigtau2, 12, temp_loads.an, temp_loads.bn, temp_loads.SigmaQ, rlb, 1);



writematrix(sdtsm6.resid, './results/eseM6.csv');
writematrix(sdtsm10.resid, './results/eseM10.csv');
writematrix(sdtsm20.resid, './results/eseM20.csv');
writematrix(sdtsm30.resid, './results/eseM30.csv');
writematrix(sdtsm40.resid, './results/eseM40.csv');


sqrt(mean(sdtsm6.resid.^2))
sqrt(mean(rw40.resid.^2))
%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%% Quarterly Forecast of 10-y Foreward Rate Using SDTSM %%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% The BR model uses quarterly data while our SDTSM uses monthly data
% As a result, we adopt the averaged value of three yields in each 
% quarter as a quarterly data

% % Overlapping periods: 197110 - 201803
% % Forecast starts at 197609 and ends at 201603
% % Yn2 = Yn(1:558,:);
% FR2 = FR(1:558,:);
% dates = USdata{1:558,1};

% Overlapping periods: 197201 - 201803
% Forecast starts at 197612 and ends at 201603
FR2 = FR(3:557,:);
dates = USdata{3:557,1};


% %%%%%%% 1-Quarter Forecast
% % sdtsm_Q1 = forecast_sdtsm_Quarter(FR2, temp_phi, temp_omega, temp_sige2, temp_gamma,...
% %     temp_Xbar, temp_sigtau2, 12, temp_loads.an, temp_loads.bn, temp_loads.SigmaQ, rlb, 3);
% 
% ese1 = forecast_sdtsm2(FR2, temp_phi, temp_omega, temp_sige2, temp_gamma,...
%     temp_Xbar, temp_sigtau2, 12, temp_loads.an, temp_loads.bn, temp_loads.SigmaQ, rlb, 1);
% ese2 = forecast_sdtsm2(FR2, temp_phi, temp_omega, temp_sige2, temp_gamma,...
%     temp_Xbar, temp_sigtau2, 12, temp_loads.an, temp_loads.bn, temp_loads.SigmaQ, rlb, 2);
% ese3 = forecast_sdtsm2(FR2, temp_phi, temp_omega, temp_sige2, temp_gamma,...
%     temp_Xbar, temp_sigtau2, 12, temp_loads.an, temp_loads.bn, temp_loads.SigmaQ, rlb, 3);
% ese = (ese1.resid + ese2.resid + ese3.resid)/3;
% ese_Q1 = nan(length(ese)/3,1);
% for i = 1:length(ese)
%     if mod(i,3)==0
%         ese_Q1(i/3,1) = ese(i,1);
%     end
% end
% sqrt(sum(ese_Q1.^2)/length(ese_Q1))
% writematrix(ese_Q1, './results/sdtsm_Q1_full.csv');
% 
% %%%%%%%%%% 2-Quarter Forecast
% ese1 = forecast_sdtsm2(FR2, temp_phi, temp_omega, temp_sige2, temp_gamma,...
%     temp_Xbar, temp_sigtau2, 12, temp_loads.an, temp_loads.bn, temp_loads.SigmaQ, rlb, 4);
% ese2 = forecast_sdtsm2(FR2, temp_phi, temp_omega, temp_sige2, temp_gamma,...
%     temp_Xbar, temp_sigtau2, 12, temp_loads.an, temp_loads.bn, temp_loads.SigmaQ, rlb, 5);
% ese3 = forecast_sdtsm2(FR2, temp_phi, temp_omega, temp_sige2, temp_gamma,...
%     temp_Xbar, temp_sigtau2, 12, temp_loads.an, temp_loads.bn, temp_loads.SigmaQ, rlb, 6);
% ese = (ese1.resid + ese2.resid + ese3.resid)/3;
% ese_Q2 = nan(length(ese)/3,1);
% for i = 1:length(ese)
%     if mod(i,3)==0
%         ese_Q2(i/3,1) = ese(i,1);
%     end
% end
% sqrt(sum(ese_Q2.^2)/length(ese_Q2))
% writematrix(ese_Q2, './results/sdtsm_Q2_full.csv');
% 
% %%%%%% 4-Quarter Forecast
% ese1 = forecast_sdtsm2(FR2, temp_phi, temp_omega, temp_sige2, temp_gamma,...
%     temp_Xbar, temp_sigtau2, 12, temp_loads.an, temp_loads.bn, temp_loads.SigmaQ, rlb, 10);
% ese2 = forecast_sdtsm2(FR2, temp_phi, temp_omega, temp_sige2, temp_gamma,...
%     temp_Xbar, temp_sigtau2, 12, temp_loads.an, temp_loads.bn, temp_loads.SigmaQ, rlb, 11);
% ese3 = forecast_sdtsm2(FR2, temp_phi, temp_omega, temp_sige2, temp_gamma,...
%     temp_Xbar, temp_sigtau2, 12, temp_loads.an, temp_loads.bn, temp_loads.SigmaQ, rlb, 12);
% ese = (ese1.resid + ese2.resid + ese3.resid)/3;
% ese_Q4 = nan(length(ese)/3,1);
% for i = 1:length(ese)
%     if mod(i,3)==0
%         ese_Q4(i/3,1) = ese(i,1);
%     end
% end
% sqrt(sum(ese_Q4.^2)/length(ese_Q4))
% writematrix(ese_Q4, './results/sdtsm_Q4_full.csv');
% 
% %%%%%% 8-Quarter Forecast
% ese1 = forecast_sdtsm2(FR2, temp_phi, temp_omega, temp_sige2, temp_gamma,...
%     temp_Xbar, temp_sigtau2, 12, temp_loads.an, temp_loads.bn, temp_loads.SigmaQ, rlb, 22);
% ese2 = forecast_sdtsm2(FR2, temp_phi, temp_omega, temp_sige2, temp_gamma,...
%     temp_Xbar, temp_sigtau2, 12, temp_loads.an, temp_loads.bn, temp_loads.SigmaQ, rlb, 23);
% ese3 = forecast_sdtsm2(FR2, temp_phi, temp_omega, temp_sige2, temp_gamma,...
%     temp_Xbar, temp_sigtau2, 12, temp_loads.an, temp_loads.bn, temp_loads.SigmaQ, rlb, 24);
% ese = (ese1.resid + ese2.resid + ese3.resid)/3;
% ese_Q8 = nan(length(ese)/3,1);
% for i = 1:length(ese)
%     if mod(i,3)==0
%         ese_Q8(i/3,1) = ese(i,1);
%     end
% end
% sqrt(sum(ese_Q8.^2)/length(ese_Q8))
% writematrix(ese_Q8, './results/sdtsm_Q8_full.csv');
% 
% 
% %%%%%% 10-Quarter Forecast
% % ese1 = forecast_sdtsm2(FR2, temp_phi, temp_omega, temp_sige2, temp_gamma,...
% %     temp_Xbar, temp_sigtau2, 12, temp_loads.an, temp_loads.bn, temp_loads.SigmaQ, rlb, 28);
% % ese2 = forecast_sdtsm2(FR2, temp_phi, temp_omega, temp_sige2, temp_gamma,...
% %     temp_Xbar, temp_sigtau2, 12, temp_loads.an, temp_loads.bn, temp_loads.SigmaQ, rlb, 29);
% % ese3 = forecast_sdtsm2(FR2, temp_phi, temp_omega, temp_sige2, temp_gamma,...
% %     temp_Xbar, temp_sigtau2, 12, temp_loads.an, temp_loads.bn, temp_loads.SigmaQ, rlb, 30);
% % ese = (ese1.resid + ese2.resid + ese3.resid)/3;
% % ese_Q10 = nan(length(ese)/3,1);
% % for i = 1:length(ese)
% %     if mod(i,3)==0
% %         ese_Q10(i/3,1) = ese(i,1);
% %     end
% % end
% % sqrt(sum(ese_Q10.^2)/length(ese_Q10))
% % writematrix(ese_Q10, 'eseQ10.csv');


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Longer forecasting period
% Overlapping periods: 197110 - 201803
% Forecast starts at 197609 and ends at 200803


%%%%%% 4-Quarter Forecast
ese1 = forecast_sdtsm3(FR2, temp_phi, temp_omega, temp_sige2, temp_gamma,...
    temp_Xbar, temp_sigtau2, 12, temp_loads.an, temp_loads.bn, temp_loads.SigmaQ, rlb, 10);
ese2 = forecast_sdtsm3(FR2, temp_phi, temp_omega, temp_sige2, temp_gamma,...
    temp_Xbar, temp_sigtau2, 12, temp_loads.an, temp_loads.bn, temp_loads.SigmaQ, rlb, 11);
ese3 = forecast_sdtsm3(FR2, temp_phi, temp_omega, temp_sige2, temp_gamma,...
    temp_Xbar, temp_sigtau2, 12, temp_loads.an, temp_loads.bn, temp_loads.SigmaQ, rlb, 12);
ese = (ese1.resid + ese2.resid + ese3.resid)/3;
temp4 = nan(length(ese)/3,1);
ese_Q4 = nan(length(ese)/3,1);
for i = 1:length(ese)
    if mod(i,3)==0
        ese_Q4(i/3,1) = ese(i,1);
        temp4(i/3,1) = ese1.resid(i,1);
    end
end
sqrt(sum(ese_Q4.^2)/length(ese_Q4))
writematrix(ese_Q4, './results/eseQ4.csv');
% sqrt(sum(temp4.^2)/length(temp4))
% writematrix(ese_Q4, 'eseQ4_end.csv') 
% writematrix(temp4, 'eseQ4_start.csv') 

%%%%%% 10-Quarter Forecast
ese1 = forecast_sdtsm3(FR2, temp_phi, temp_omega, temp_sige2, temp_gamma,...
    temp_Xbar, temp_sigtau2, 12, temp_loads.an, temp_loads.bn, temp_loads.SigmaQ, rlb, 28);
ese2 = forecast_sdtsm3(FR2, temp_phi, temp_omega, temp_sige2, temp_gamma,...
    temp_Xbar, temp_sigtau2, 12, temp_loads.an, temp_loads.bn, temp_loads.SigmaQ, rlb, 29);
ese3 = forecast_sdtsm3(FR2, temp_phi, temp_omega, temp_sige2, temp_gamma,...
    temp_Xbar, temp_sigtau2, 12, temp_loads.an, temp_loads.bn, temp_loads.SigmaQ, rlb, 30);
ese = (ese1.resid + ese2.resid + ese3.resid)/3;
ese_Q10 = nan(length(ese)/3,1);
temp10 = nan(length(ese)/3,1);
for i = 1:length(ese)
    if mod(i,3)==0
        ese_Q10(i/3,1) = ese(i,1);
        temp10(i/3,1) = ese1.resid(i,1);
    end
end
sqrt(sum(ese_Q10.^2)/length(ese_Q10))
writematrix(ese_Q10, './results/eseQ10.csv');
% sqrt(sum(temp10.^2)/length(temp10))
% writematrix(ese_Q10, 'eseQ10_end.csv')
% writematrix(temp10, 'eseQ10_start.csv') 

%%%%%% 20-Quarter Forecast
ese1 = forecast_sdtsm3(FR2, temp_phi, temp_omega, temp_sige2, temp_gamma,...
    temp_Xbar, temp_sigtau2, 12, temp_loads.an, temp_loads.bn, temp_loads.SigmaQ, rlb, 58);
ese2 = forecast_sdtsm3(FR2, temp_phi, temp_omega, temp_sige2, temp_gamma,...
    temp_Xbar, temp_sigtau2, 12, temp_loads.an, temp_loads.bn, temp_loads.SigmaQ, rlb, 59);
ese3 = forecast_sdtsm3(FR2, temp_phi, temp_omega, temp_sige2, temp_gamma,...
    temp_Xbar, temp_sigtau2, 12, temp_loads.an, temp_loads.bn, temp_loads.SigmaQ, rlb, 60);
ese = (ese1.resid + ese2.resid + ese3.resid)/3;
ese_Q20 = nan(length(ese)/3,1);
temp20 = nan(length(ese)/3,1);
for i = 1:length(ese)
    if mod(i,3)==0
        ese_Q20(i/3,1) = ese(i,1);
        temp20(i/3,1) = ese1.resid(i,1);
    end
end
sqrt(sum(ese_Q20.^2)/length(ese_Q20))
writematrix(ese_Q20, './results/eseQ20.csv');
% sqrt(sum(temp20.^2)/length(temp20))
% writematrix(ese_Q20, 'eseQ20_end.csv')
% writematrix(temp20, 'eseQ20_start.csv') 


%%%%%% 30-Quarter Forecast
ese1 = forecast_sdtsm3(FR2, temp_phi, temp_omega, temp_sige2, temp_gamma,...
    temp_Xbar, temp_sigtau2, 12, temp_loads.an, temp_loads.bn, temp_loads.SigmaQ, rlb, 88);
ese2 = forecast_sdtsm3(FR2, temp_phi, temp_omega, temp_sige2, temp_gamma,...
    temp_Xbar, temp_sigtau2, 12, temp_loads.an, temp_loads.bn, temp_loads.SigmaQ, rlb, 89);
ese3 = forecast_sdtsm3(FR2, temp_phi, temp_omega, temp_sige2, temp_gamma,...
    temp_Xbar, temp_sigtau2, 12, temp_loads.an, temp_loads.bn, temp_loads.SigmaQ, rlb, 90);
ese = (ese1.resid + ese2.resid + ese3.resid)/3;
ese_Q30 = nan(length(ese)/3,1);
temp30 = nan(length(ese)/3,1);
for i = 1:length(ese)
    if mod(i,3)==0
        ese_Q30(i/3,1) = ese(i,1);
        temp30(i/3,1) = ese1.resid(i,1);
    end
end
sqrt(sum(ese_Q30.^2)/length(ese_Q30))
writematrix(ese_Q30, './results/eseQ30.csv');
% sqrt(sum(temp30.^2)/length(temp30))
% writematrix(ese_Q30, 'eseQ30_end.csv') 
% writematrix(temp30, 'eseQ30_start.csv') 

%%%%%% 40-Quarter Forecast
ese1 = forecast_sdtsm3(FR2, temp_phi, temp_omega, temp_sige2, temp_gamma,...
    temp_Xbar, temp_sigtau2, 12, temp_loads.an, temp_loads.bn, temp_loads.SigmaQ, rlb, 118);
ese2 = forecast_sdtsm3(FR2, temp_phi, temp_omega, temp_sige2, temp_gamma,...
    temp_Xbar, temp_sigtau2, 12, temp_loads.an, temp_loads.bn, temp_loads.SigmaQ, rlb, 119);
ese3 = forecast_sdtsm3(FR2, temp_phi, temp_omega, temp_sige2, temp_gamma,...
    temp_Xbar, temp_sigtau2, 12, temp_loads.an, temp_loads.bn, temp_loads.SigmaQ, rlb, 120);
ese = (ese1.resid + ese2.resid + ese3.resid)/3;
ese_Q40 = nan(length(ese)/3,1);
temp40 = nan(length(ese)/3,1);
for i = 1:length(ese)
    if mod(i,3)==0
        ese_Q40(i/3,1) = ese(i,1);
        temp40(i/3,1) = ese1.resid(i,1);
    end
end
sqrt(sum(ese_Q40.^2)/length(ese_Q40))
writematrix(ese_Q40, './results/eseQ40.csv');
% sqrt(sum(temp40.^2)/length(temp40))
% writematrix(ese_Q40, 'eseQ40_end.csv')  
% writematrix(temp40, 'eseQ40_start.csv') 

sqrt(mean(ese_Q40.^2))
%%
% Forecast 1-year forward rate in 10 years
% Compare SDTSM and BR

%%%%%% 40-Quarter Forecast
ese1 = forecast_sdtsm4(FR2, temp_phi, temp_omega, temp_sige2, temp_gamma,...
    temp_Xbar, temp_sigtau2, 12, temp_loads.an, temp_loads.bn, temp_loads.SigmaQ, rlb, 118);
ese2 = forecast_sdtsm4(FR2, temp_phi, temp_omega, temp_sige2, temp_gamma,...
    temp_Xbar, temp_sigtau2, 12, temp_loads.an, temp_loads.bn, temp_loads.SigmaQ, rlb, 119);
ese3 = forecast_sdtsm4(FR2, temp_phi, temp_omega, temp_sige2, temp_gamma,...
    temp_Xbar, temp_sigtau2, 12, temp_loads.an, temp_loads.bn, temp_loads.SigmaQ, rlb, 120);
ese = (ese1.resid + ese2.resid + ese3.resid)/3;
ese_Q40 = nan(length(ese)/3,1);
temp40 = nan(length(ese)/3,1);
for i = 1:length(ese)
    if mod(i,3)==0
        ese_Q40(i/3,1) = ese(i,1);
        temp40(i/3,1) = ese1.resid(i,1);
    end
end
sqrt(sum(ese_Q40.^2)/length(ese_Q40))
writematrix(ese_Q40, './results/SDTSM_eseQ40_F1.csv');






%% generate quartely smoothed forward rate from LW

% end-of-quarter
temp = USdata(2:end,2:end);
temp2 = temp{1:3:604,1:18};
writematrix(temp2, 'LW_forward_quarterend.csv');

% quarter mean
temp3 = nan(size(USdata,1),18);
temp3(:,1) = USdata.yyyymm;
temp3(2,2:end) = mean(USdata{1:2,3:19},1);
for i = 5:3:605
    temp3(i,2:end) = mean(USdata{i-2:i,3:19},1);
end
temp3 = temp3(2:3:605,:);
writematrix(temp3, 'LW_forward_quartermean.csv');


%% Real OOS forecast by SDTSM
clear all
close all
clc

% Yields and Forward rate from LW: 197111-202204
% mats = [3, 6, 12:12:120]';
usyields = readtable("US_yields_LW.csv");
Yn = usyields{126:731,2:13};

% % Generate monthly Forward Rate 197111-202204
load USdata_LW2.mat
FR = USdata{:,3:14};

load("./data/temp_theta.mat")
dt = 1;
rlb = 0;
N = 3;
mats = [3, 6, 12:12:120]';
mcmc_ind = 50000;


% ind = 61829
% temp_theta = nan(25,1);
% temp_theta(1:9) = mean(Phi(ind,:),1);
% temp_theta(10:12) = mean(lamQ(ind,:),1);
% temp_omega = mean(Omega(ind,:),1);
% temp_omega = [temp_omega(1) temp_omega(2) temp_omega(3);...
%     temp_omega(2) temp_omega(4) temp_omega(5);...
%     temp_omega(3) temp_omega(5) temp_omega(6)];
temp_omega2 = mean(Sigma_tilde(ind,:),1);
temp_omega2 = [temp_omega2(1) temp_omega2(2) temp_omega2(3);...
    temp_omega2(2) temp_omega2(4) temp_omega2(5);...
    temp_omega2(3) temp_omega2(5) temp_omega2(6)];
temp_Sigma_tilde = temp_omega2;
% % temp_omega - mean(sigtau2(ind,:),1)*mean(gamma(ind,:),1)' * mean(gamma(ind,:),1)
% ltri = tril(ones(3));
% % temp_sigma = chol(temp_omega,'lower') ;
% temp_theta(13:18) = temp_omega2(ltri==1) ;
% temp_theta(19) =  mean(kinfQ(ind,:),1);
% temp_theta(20) = mean(sige2(ind,:),1);
% temp_theta(21:22) = mean(gamma(ind,2:3),1);
% temp_theta(23:24) = mean(Xbar(ind,2:3),1);
% temp_theta(25) =  mean(sigtau2(ind,:),1);
% save("temp_theta.mat","temp_theta")

load('temp100000.mat', 'Omega','Phi','kinfQ','lamQ','sige2','gamma','Xbar','sigtau2');

temp_omega = mean(Omega(mcmc_ind+1,:),1);
temp_omega = [temp_omega(1,1) temp_omega(1,2) temp_omega(1,3);...
    temp_omega(1,2) temp_omega(1,4) temp_omega(1,5);...
    temp_omega(1,3) temp_omega(1,5) temp_omega(1,6)];
temp_phi = reshape(mean(Phi(mcmc_ind+1:end,:),1),3,3);
temp_loads = SRTSM_loadings(mean(kinfQ(mcmc_ind+1:end)), diag(mean(lamQ(mcmc_ind+1:end,:),1)), temp_omega, mats, dt);
temp_sige2 =  mean(sige2(mcmc_ind+1:end));
temp_gamma = mean(gamma(mcmc_ind+1:end,:),1)';
temp_Xbar =  mean(Xbar(mcmc_ind+1:end,:),1)';
temp_sigtau2 = mean(sigtau2(mcmc_ind+1:end));

temp_theta = nan(25,1);
temp_theta(1:9) = temp_phi;
temp_theta(10:12) = mean(lamQ(mcmc_ind+1:end,:),1);
ltri = tril(ones(N));
temp_theta(13:18) = temp_Sigma_tilde(ltri==1);
temp_theta(19) = mean(kinfQ(mcmc_ind+1:end));
temp_theta(20) = mean(temp_sige2,1);
temp_theta(21:22) = mean(temp_gamma(2:3,1),1);
temp_theta(23:24) = mean(temp_Xbar(2:3,1),1);
temp_theta(25) =  mean(temp_sigtau2,1);
save("temp_theta2.mat","temp_theta")


% temp.Phi = reshape(mean(Phi(mcmc_ind+1:end,:),1),3,3);
% temp.lamQ = mean(lamQ(mcmc_ind+1:end,:),1)';
% temp.dlamQ = [1-temp.lamQ(1); temp.lamQ(1)-temp.lamQ(2); temp.lamQ(2)-temp.lamQ(3)];
% temp_omega = mean(Omega(mcmc_ind+1:end,:),1);
% temp.omega = [temp_omega(1) temp_omega(2) temp_omega(3);...
%     temp_omega(2) temp_omega(4) temp_omega(5);...
%     temp_omega(3) temp_omega(5) temp_omega(6)];
% temp.kinfQ = mean(kinfQ(mcmc_ind+1:end,:),1);
% temp.sige2 = mean(sige2(mcmc_ind+1:end,:),1);
% temp.gamma = mean(gamma(mcmc_ind+1:end,:),1);
% temp.Xbar = mean(Xbar(mcmc_ind+1:end,:),1);
% temp.sigtau2 = mean(sigtau2(mcmc_ind+1:end,:),1);
% theta_temp = pars2theta(temp);

pars_mode = theta2pars(temp_theta,mats,dt);



% Sigma_tilde
% prior0.kappa = 1000 ;    % low precision
prior0.Psi = pars_mode.Sigma_tilde ;
prior0.PsiCov = 0.1^2*eye(N^2);
% ap
prior0.ap_mean = [pars_mode.a; pars_mode.p];
prior0.ap_cov = diag([0.01; 0.01; 0.1; 0.1].^2) ;
prior0.ap_cov_inv = prior0.ap_cov\eye(4);
% sigtau2
% prior0.alpha_eta = 1000;  % tightness. Del Negro et al use 100
% prior0.mean_sigtau2  = pars_mode.sigtau2;   
% prior0.delta_eta = prior0.mean_sigtau2 * (prior0.alpha_eta - 2) ;
prior0.mean_sigtau2  = log(pars_mode.sigtau2);   
prior0.delta_eta = 0.1^2;
% Phi
prior0.Phi = pars_mode.Phi;
prior0.omega = 0.1^2*ones(N^2,1);
prior0.Omega = diag(prior0.omega) ;
% sige2
% prior0.alpha_e = 4;
% prior0.alpha_e = 1000;
% prior0.mean_sige2 = pars_mode.sige2;
% prior0.delta_e = prior0.mean_sige2 * (prior0.alpha_e - 2) ;
prior0.mean_sige2 = log(pars_mode.sige2);
prior0.delta_e = 0.1^2 ;
% kinfQ
prior0.kmean = pars_mode.kinfQ;
prior0.kvar = 0.01^2;


temp = llk_sdtsm(temp_theta, FR', 0, rlb, mats, dt);
sdtsm_rt = forecast_sdtsm_rt(FR, temp_theta, prior0, mats, 12, dt, rlb, [6;10;20;30;40]);
save("./results/sdtsm_rt.mat","sdtsm_rt")


writematrix(sdtsm_rt.resid(:,1), './results/eseM6_rt.csv');
writematrix(sdtsm_rt.resid(:,2), './results/eseM10_rt.csv');
writematrix(sdtsm_rt.resid(:,3), './results/eseM20_rt.csv');
writematrix(sdtsm_rt.resid(:,4), './results/eseM30_rt.csv');
writematrix(sdtsm_rt.resid(:,5), './results/eseM40_rt.csv');





