% Out-of-Sample Forecast of 10-year Forward Rate
% Sample: 198301 - 202203
% maturities: 1-10y
clear all
close all
clc

% read yields 
mats = (12:12:120)';
Yn = genYn('UK_yields2.xlsx',mats/12);

% Generate monthly Forward Rate 
load UKdata.mat
FR = UKdata{122:end,3:12};


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

mcmc_ind = 50000;
dt = 1;
rlb = 0;
N =3;

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 %%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
rw6 = forecast_rw(FR(:,10),6);
rw10 = forecast_rw(FR(:,10),10);
rw20 = forecast_rw(FR(:,10),20);
rw30 = forecast_rw(FR(:,10),30);
rw40 = forecast_rw(FR(:,10),40);

writematrix(rw6.resid, 'rw6.csv');
writematrix(rw10.resid, 'rw10.csv');
writematrix(rw20.resid, 'rw20.csv');
writematrix(rw30.resid, 'rw30.csv');
writematrix(rw40.resid, 'rw40.csv');


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%% JSZ forecast %%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
jsz6 = forecast_jsz3(FR(:,10),Yn{:,2:11}, mats', 6);
jsz10 = forecast_jsz3(FR(:,10),Yn{:,2:11}, mats', 10);
jsz20 = forecast_jsz3(FR(:,10),Yn{:,2:11}, mats', 20);
jsz30 = forecast_jsz3(FR(:,10),Yn{:,2:11}, mats', 30);
jsz40 = forecast_jsz3(FR(:,10),Yn{:,2:11}, mats', 40);

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

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%% SDTSM forecast %%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
sdtsm6 = forecast_sdtsm(FR, temp_phi, temp_omega, temp_sige2, temp_gamma,...
    temp_Xbar, temp_sigtau2, 10, 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, 10, 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, 10, 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, 10, 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, 10, temp_loads.an, temp_loads.bn, temp_loads.SigmaQ, rlb, 40);

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



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

load UKdata.mat
FR = UKdata{122:end,3:12};

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

dt = 1;
rlb = 0;
N = 3;
mats = [12:12:120]';
mcmc_ind = 50000;
% plot(mean(tau(50000:end,:),1))

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_Sigma_tilde = mean(Sigma_tilde(mcmc_ind+1:end,:),1);

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;
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")
save("temp_theta2_rr2.mat","temp_theta")

load("temp_theta2.mat")

pars_mode = theta2pars(temp_theta,mats,dt);

% Sigma_tilde
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.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.mean_sige2 = log(pars_mode.sige2);
prior0.delta_e = 0.1^2 ;
% kinfQ
prior0.kmean = pars_mode.kinfQ;
prior0.kvar = 0.01^2;


sdtsm_rt = forecast_sdtsm_rt(FR, temp_theta, prior0, mats, 10, dt, rlb, [6;10;20;30;40]);
% save("./results/sdtsm_rt.mat","sdtsm_rt")
save("./results/sdtsm_rt_rr2.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');














