%% PART 1: Yield Curve Modeling
%% Run bootstrapping code and NS regressions to find yields given bond quotes 1991-2014
clear;
clc;
%Initiate tables
coefficients_2014 = zeros(288,3);
rsquared_2014 = zeros(288,1);
data_2014 = zeros(17,288);

%The input is data on price quotes on noncallable bonds, noncallable notes 
%and T-bills filtered for securities with liquidity problems, i.e. bonds
%and notes with less than one year to maturity, and bills with less than
%one month to maturity. The prices are clean (i.e. not including accrued
%interest) and the output is continously compounded. We use an 
%actual/actual (ICMA) day-count convention in accordance with the US
%Treasury. Lastly, we use an end-of-month-rule that assumes that a bond's 
%coupon payment date is always the last actual day of the month. We will
%now run for loop which first estimates unsmoothed yields using 
%bootstrapping (zbtprice), then fits a smmoth NS function to these 
%unsmoothed yields.

for i=1:288%We have 300 settle dates
T = readtable('DataBonds5.xlsx','ReadVariableNames',true);
ind1 = table2array(T(:,1)) == i; %The settle dates have IDs which we loop on
B = T(ind1,2:3);%Create matrix B=(nxm), n=#maturities, m={maturity date, CR}
B.TTM = datetime(table2array(B(:,1)),'InputFormat','dd.MM.yyyy');
B(:,3) = array2table((zeros(size(B,1),1)+100));%Column 3: Face value = 100
B(:,4) = array2table((zeros(size(B,1),1)+2));%Column 4: Semiannual coupons
B(:,5) = array2table((zeros(size(B,1),1)+8));%Column 5: act/act(ICMA) day count
P = T(ind1,5);%Prices observed on settle date 
S = T(ind1,4);%The current settle (quote) date
S = S(1,1);
S = datetime(table2array(S),'InputFormat','dd.MM.yyyy');
P = table2array(P);
OutputCompounding = -1;%Continuous output compounding

[ZeroRates,CurveDates] = zbtprice(B,P,S,OutputCompounding);%Obtain zero yields
Mat = (daysact(S, CurveDates))./(365/12);%Obtain time to maturity. Use in NS
Mat2 = (daysact(S, CurveDates))./(365/12);%This will be used in interpolation
Yield = table(ZeroRates.*100, Mat);%Save both the rates and mat in a table. Use in NS
Yield2 = table(ZeroRates.*100, Mat2);%This will be used in interpolation
%Decide on the longest maturity to be included.
%ind5 = Yield.Mat >= 115 & Yield.Mat<=125;
%longestmat(i) = sum(ind5);
%Conclusion: We stop at 120 months and start at 3. Hence, we remove the
%observations where m<3 and m>120
ind2 = Yield.Mat < 3;%Remove bonds with <2.5m matutrity due to liq. probs.
Yield(ind2,:) = [];
Mat(ind2) = [];
ind3 = Yield.Mat > 120;%Remove bonds with >130m since we only want NS to fit in the range 3m-120m
Yield(ind3,:) = [];
Mat(ind3) = [];

%Nelson-Siegel
m = Yield.Mat;
Yield.X_1 = (1-exp(-m.*0.0609))./(m.*0.0609);
Yield.X_2 = ((1-exp(-m.*0.0609))./(m.*0.0609))-exp(-m.*0.0609);
model = fitlm(Yield, 'Var1 ~ X_1 + X_2');%Run NS regressions
coefficients_2014(i,:) = model.Coefficients.Estimate';%Save coeffs
rsquared_2014(i,:) = model.Rsquared.Ordinary%Save R2

%To make comparison between NS and actual data easier, we linearly
%interpolate the observed yields to get a fixed set of yields
x = Mat2;
xq1 = 3:1:125;
y = table2array(Yield2(:,1));
interpolated = pchip(x,y,xq1)';%Run PCHIP instead??
interpolated(:,2) = xq1;
m_2014 = [3 6 9 12 15 18 21 24 30 36 48 60 72 84 96 108 120]';%These are the fixed maturities
ind4 = ismember(interpolated(:,2), m_2014);%This identifies the fixed maturities
data_2014(:,i) = interpolated(ind4,1);%This extracts the fixed maturities from the interpolated data
end

b1_2014 = coefficients_2014(:,1);
b2_2014 = coefficients_2014(:,2);
b3_2014 = coefficients_2014(:,3);
mean(rsquared_2014)

%% Nelson-Siegel 2015-2020
clear data_2020
data_2020 = readtable('USTresData.xlsx','ReadVariableNames',false);
data_2020(1:300,:) = [];
%Set to 301 if you have Matlab R2019b. 
%It seems that Matlab R2020a imports the US Treasury data differently than
%R2019a. While R2019b imports the variable names as first row, R2020a does
%not. Hence, if you are using R2020a, please only delete the 300 first
%rows. If you are using R2019a, please delete the 301 first rows. 
data_2020(:,1:3) = [];
data_2020(:,9:10) =[];

% Run str2double if you have R2019b
data_2020.Var4 = str2double(data_2020.Var4);
data_2020.Var5 = str2double(data_2020.Var5);
data_2020.Var6 = str2double(data_2020.Var6);
data_2020.Var7 = str2double(data_2020.Var7);
data_2020.Var8 = str2double(data_2020.Var8);
data_2020.Var9 = str2double(data_2020.Var9);
data_2020.Var10 = str2double(data_2020.Var10);
data_2020.Var11 = str2double(data_2020.Var11);

data_2020 = rows2vars(data_2020);
data_2020(:,1) = [];

data_2020 = table2cell(data_2020);
data_2020 = cell2table(data_2020, 'VariableNames',{'JAN15','FEB15','MAR15','APR15','MAY15','JUN15','JUL15','AUG15'...
    ,'SEP15','OCT15','NOV15','DEC15','JAN16','FEB16','MAR16','APR16','MAY16','JUN16','JUL16','AUG16'...
    ,'SEP16','OCT16','NOV16','DEC16','JAN17','FEB17','MAR17','APR17','MAY17','JUN17','JUL17','AUG17'...
    ,'SEP17','OCT17','NOV17','DEC17','JAN18','FEB18','MAR18','APR18','MAY18','JUN18','JUL18','AUG18'...
    ,'SEP18','OCT18','NOV18','DEC18','JAN19','FEB19','MAR19','APR19','MAY19','JUN19','JUL19','AUG19'...
    ,'SEP19','OCT19','NOV19', 'DEC19', 'JAN20'});

m_2020 = [3 6 12 24 36 60 84 120]';
VarNames_2020 = data_2020.Properties.VariableNames(1,[1:(size(data_2020,2))]);

coefficients_2020 = zeros((size(data_2020,2)),3);
rsquared_2020 = zeros((size(data_2020,2)),1);

for v = 1:numel(VarNames_2020); % loop over number of elements in VarNames
    current_column = VarNames_2020{v}; % Get the varname for this loop

    data_2020.X_1_2020 = (1-exp(-m_2020.*0.0609))./(m_2020.*0.0609);
    data_2020.X_2_2020 = ((1-exp(-m_2020.*0.0609))./(m_2020.*0.0609))-exp(-m_2020.*0.0609);

    model_spec = [current_column ' ~ X_1_2020 + X_2_2020'];
    model = fitlm(data_2020, model_spec);
    
    coefficients_2020(v,:) = model.Coefficients.Estimate';
    rsquared_2020(v) = model.Rsquared.Ordinary'
end

b1_2020 = coefficients_2020(:,1);
b2_2020 = coefficients_2020(:,2);
b3_2020 = coefficients_2020(:,3);

mean(rsquared_2020(1:60))
%% Combine the model outputs

b1 = [b1_2014; b1_2020]
b2 = [b2_2014; b2_2020]
b3 = [b3_2014; b3_2020]

r_2 = [rsquared_2014; rsquared_2020]
mean(r_2)

%% NS fit statistics
mean(r_2)
mean(rsquared_2014)
mean(rsquared_2020)
min(r_2)
max(r_2)
worst = find(r_2 == min(r_2))
best = find(r_2 == max(r_2))
%The worst fit is in obs. 181
%The best fit is in obs. 292

ind0 = r_2 >= 0.9;
sum(ind0)
ind1 = r_2 < 0.9;
sum(ind1)
ind2 = r_2 < 0.8;
sum(ind2)
ind3 = r_2 < 0.7;
sum(ind3)
ind4 = r_2 < 0.6;
sum(ind4)
ind5 = r_2 < 0.5;
sum(ind5)
ind6 = r_2 < 0.4;
sum(ind6)
ind7 = r_2 < 0.3;
sum(ind7)
ind8 = r_2 < 0.2;
sum(ind8)
ind9 = r_2 < 0.1;
sum(ind9)
%Maybe show this in a table? Make in latex. 

%% Best and worst yield curve replication
badact = data_2014(:,197);
goodact = table2array(data_2020(:,(best-288)));
badb1 = b1(197);
badb2 = b2(197);
badb3 = b3(197);
badns = badb1 + badb2.*x1 + badb3.*x2;

goodb1 = b1(best);
goodb2 = b2(best);
goodb3 = b3(best);
goodns = goodb1 + goodb2.*x1 + goodb3.*x2;

%Worst fit 1: obs 197
badact = data_2014(:,197);
badb1 = b1(197);
badb2 = b2(197);
badb3 = b3(197);
badns = badb1 + badb2.*x1 + badb3.*x2
figure
plot(m, badns, 'k')
hold on
sz = 30
scatter(m_2014, badact, sz,'filled','d')
xlabel('Maturity in months')
ylabel('Yield in percent')
legend({'Nelson-Siegel', 'Actual'}, 'Location','southeast')
dim = [.72 .01 .3 .3];
str = 'R^2 = 0.110';
annotation('textbox',dim,'String',str,'FitBoxToText','on','EdgeColor', [1 1 1]);

%Worst fit 2: obs 181
badact = data_2014(:,181);
badb1 = b1(181);
badb2 = b2(181);
badb3 = b3(181);
badns = badb1 + badb2.*x1 + badb3.*x2;
figure
plot(m, badns, 'k')
hold on
sz = 30
scatter(m_2014, badact, sz,'filled','d')
xlabel('Maturity in months')
ylabel('Yield in percent')
legend({'Nelson-Siegel', 'Actual'}, 'Location','southeast')
dim = [.72 .01 .3 .3];
str = 'R^2 = 0.187';
annotation('textbox',dim,'String',str,'FitBoxToText','on','EdgeColor', [1 1 1]);

%Best fit
figure
plot(m, goodns, 'k')
hold on
sz = 30
scatter(m_2020, goodact,sz,'filled','d')
title('Best Nelson-Siegel fit: April 2015')
xlabel('Maturity in months')
ylabel('Yield To Maturity')
legend({'Actual', 'Nelson-Siegel'}, 'Location','southeast')
dim = [.72 .01 .3 .3];
str = 'R^2 = 0.9998';
annotation('textbox',dim,'String',str,'FitBoxToText','on','EdgeColor', [1 1 1]);

%% Find special shapes where NS do a good job
%First, create a matrix with each NS curve
m = [3:1:120]';
x1 = (1-exp(-m.*0.0609))./(m.*0.0609);
x2 = ((1-exp(-m.*0.0609))./(m.*0.0609))-exp(-m.*0.0609);
for i=1:349
    b1_ns = b1(i);
    b2_ns = b2(i);
    b3_ns = b3(i);
    
    NS(:,i) = b1_ns + b2_ns.*x1 + b3_ns.*x2;
    %plot(m, NS(:,i))
    %figure
end

% Find interesting yield curves by eye-balling the plots
interesting = [349 345 344 342 341 340 339 198 196 195 194 193 192 ...
    191 190 189 188 187 186 183 182 181 180 179 119 118 117 116 115 114 ...
    113 112 111 110 109 108 92 91 90 89]';

%Select only the interesting curves with R2>90
r_2_interesting = r_2(interesting);
superinteresting = interesting(find(r_2_interesting >= 0.85));
superinteresting_2014 = superinteresting(find(superinteresting <= 288));
superinteresting_2020 = superinteresting(find(superinteresting > 288));

%Find the actual yield curves where NS R2>2
superinterestingdata_2014 = data_2014(:,superinteresting_2014);
superinterestingdata_2020 = table2array(data_2020(:,(superinteresting_2020-288))); 

%Plot to look at the interesting curves with R2>0.85 2014
% for i=1:size(superinterestingdata_2014, 2)
% plot(m_2014, superinterestingdata_2014(:,i));
% figure
% end

%Plot to look at the interesting curves with R2>0.85 2020
% for i=1:size(superinterestingdata_2020, 2)
% plot(m_2020, superinterestingdata_2020(:,i));
% figure
% end

%Find the NS that replicate the interesting actual curves
NS_interesting_2014 = NS(:,superinteresting_2014);
NS_interesting_2020 = NS(:,superinteresting_2020);

%Create plots that compare interesting actual and NS 1991-2014
for i=1:size(superinterestingdata_2014, 2)
figure
scatter(m_2014, superinterestingdata_2014(:,i))
hold on
plot(m, NS_interesting_2014(:,i))
hold off
end

%Create plots that compare interesting actual and NS 2015-2020
for i=1:size(superinterestingdata_2020, 2)
scatter(m_2020, superinterestingdata_2020(:,i))
hold on
plot(m, NS_interesting_2020(:,i))
hold off
figure
end

%Let's identify the most interesting curves
%I find the following curves to be extra interesting
extract = [2 11 14 15];
%I must now find which columns in NS correspond to these observations
h = NS_interesting_2014(:,extract);
for i=1:size(extract,2);
[~, column] = find(NS == h(:,i));
colnum(:,i) = column(1,1);
end
%Let's identify which dates produce these intereresting most interesting curves 
data_2014(:,colnum);
r_2(colnum);

%Let's create plots of these 4 curves
data_2014_extract = data_2014(:,colnum);
%March 2007: U-shape
figure
plot(m, h(:,1),'k')
hold on
sz = 30
scatter(m_2014, data_2014_extract(:,1), sz,'filled','d')
xlabel('Maturity in months')
ylabel('Yield in percent')
legend({'Nelson-Siegel', 'Actual'}, 'Location','southeast')
dim = [.72 .2 .3 .3];
str = 'R^2 = 0.955';
annotation('textbox',dim,'String',str,'FitBoxToText','on','EdgeColor', [1 1 1]);

%November 2000: Inverted and convex
figure
plot(m, h(:,2),'k')
hold on
sz = 30
scatter(m_2014, data_2014_extract(:,2), sz,'filled','d')
xlabel('Maturity in months')
ylabel('Yield in percent')
legend({'Nelson-Siegel', 'Actual'}, 'Location','east')
dim = [.72 .35 .3 .3];
str = 'R^2 = 0.970';
annotation('textbox',dim,'String',str,'FitBoxToText','on','EdgeColor', [1 1 1]);

%August 2000: Inverted
figure
plot(m, h(:,3),'k')
hold on
sz = 30
scatter(m_2014, data_2014_extract(:,3),sz,'filled','d')
xlabel('Maturity in months')
ylabel('Yield in percent')
legend({'Nelson-Siegel', 'Actual'}, 'Location','southeast')
dim = [.72 .01 .3 .3];
str = 'R^2 = 0.941';
annotation('textbox',dim,'String',str,'FitBoxToText','on','EdgeColor', [1 1 1]);
axis([0 130 5 6.8])

%April 2000: Humped
figure
plot(m, h(:,4),'k')
hold on
sz = 30
scatter(m_2014, data_2014_extract(:,4),sz,'filled','d')
xlabel('Maturity in months')
ylabel('Yield in percent')
legend({'Nelson-Siegel', 'Actual'}, 'Location','southeast')
dim = [.72 .01 .3 .3];
str = 'R^2 = 0.912';
annotation('textbox',dim,'String',str,'FitBoxToText','on','EdgeColor', [1 1 1]);
axis([0 130 5.75 6.7])

%% Nelson Siegel Descriptive statistics, yield curves: 1991-2014
%First, we look at the NS modelled yields
%m_2014 = [3 6 9 12 15 18 21 24 30 36 48 60 72 83 96 108 120]';
NS_14 = NS(:,1:288);

for i=1:size(NS_14, 1)
    mean_NS_14(i,:) = mean(NS_14(i,:));
    std_NS_14(i,:) = std(NS_14(i,:));
    min_NS_14(i,:) = min(NS_14(i,:));
    max_NS_14(i,:) = max(NS_14(i,:));
    Autocorr = autocorr(NS_14(i,:), 'NumLags', 30);
    Autocorr_NS_14_1Lag(i,:) = Autocorr(:,2);
    Autocorr_NS_14_12Lag(i,:) = Autocorr(:,13);
    Autocorr_NS_14_30Lag(i,:) = Autocorr(:,31);
end

descript_NS_14(:,1) = mean_NS_14;
descript_NS_14(:,2) = std_NS_14;
descript_NS_14(:,3) = min_NS_14;
descript_NS_14(:,4) = max_NS_14;
descript_NS_14(:,5) = Autocorr_NS_14_1Lag;
descript_NS_14(:,6) = Autocorr_NS_14_12Lag;
descript_NS_14(:,7) = Autocorr_NS_14_30Lag;
descript_NS_14(:,8) = [3:1:120]';;
descript_NS_14 = array2table(descript_NS_14, 'VariableNames', ...
    {'Mean', 'StdDev', 'Min', 'Max', 'Rho1', 'Rho2', 'Rho3', 'Mat'} );

%% Actual Yield Descriptive Statistics 1991-2014
%We now turn to the actual observed yields 
%m_2014 = [3 6 9 12 15 18 21 24 30 36 48 60 72 83 96 108 120]';
clear descript_2014

for i=1:size(data_2014, 1)
    mean_2014(i,:) = mean(data_2014(i,:));
    std_2014(i,:) = std(data_2014(i,:));
    min_2014(i,:) = min(data_2014(i,:));
    max_2014(i,:) = max(data_2014(i,:));
    Autocorr = autocorr(data_2014(i,:), 'NumLags', 30);
    Autocorr_2014_1Lag(i,:) = Autocorr(:,2);
    Autocorr_2014_12Lag(i,:) = Autocorr(:,13);
    Autocorr_2014_30Lag(i,:) = Autocorr(:,31);
end

descript_2014(:,1) = mean_2014;
descript_2014(:,2) = std_2014;
descript_2014(:,3) = min_2014;
descript_2014(:,4) = max_2014;
descript_2014(:,5) = Autocorr_2014_1Lag;
descript_2014(:,6) = Autocorr_2014_12Lag;
descript_2014(:,7) = Autocorr_2014_30Lag;
descript_2014(:,8) = m_2014;
descript_2014 = array2table(descript_2014, 'VariableNames', ...
    {'Mean', 'StdDev', 'Min', 'Max', 'Rho1', 'Rho2', 'Rho3', 'Mat'} );

%Plot of average yield curve 1991-2014
plot(m_2014, mean_2014, '-o')
axis([0 130 2.8 5])
xlabel('Maturity in months')
ylabel('Yield in percent')

%% Nelson Siegel Descriptive statistics, yield curves: 2015-2019

%First, we look at the NS modelled yields
%m = [3 6 12 24 36 60 84 120]';
NS_19 = NS(:,289:348);

for i=1:size(NS_19, 1)
    mean_NS_19(i,:) = mean(NS_19(i,:));
    std_NS_19(i,:) = std(NS_19(i,:));
    min_NS_19(i,:) = min(NS_19(i,:));
    max_NS_19(i,:) = max(NS_19(i,:));
    Autocorr = autocorr(NS_19(i,:), 'NumLags', 30);
    Autocorr_NS_19_1Lag(i,:) = Autocorr(:,2);
    Autocorr_NS_19_12Lag(i,:) = Autocorr(:,13);
    Autocorr_NS_19_30Lag(i,:) = Autocorr(:,31);
end

descript_NS_19(:,1) = mean_NS_19;
descript_NS_19(:,2) = std_NS_19;
descript_NS_19(:,3) = min_NS_19;
descript_NS_19(:,4) = max_NS_19;
descript_NS_19(:,5) = Autocorr_NS_19_1Lag;
descript_NS_19(:,6) = Autocorr_NS_19_12Lag;
descript_NS_19(:,7) = Autocorr_NS_19_30Lag;
descript_NS_19(:,8) = [3:1:120]';;
descript_NS_19 = array2table(descript_NS_19, 'VariableNames', ...
    {'Mean', 'StdDev', 'Min', 'Max', 'Rho1', 'Rho2', 'Rho3', 'Mat'} );

%% Nelson Siegel Descriptive statistics, yield curves full sample 1991-2019
clear descript_NS
for i=1:size(NS, 1)
    mean_NS(i,:) = mean(NS(i,1:348));
    std_NS(i,:) = std(NS(i,1:348));
    min_NS(i,:) = min(NS(i,1:348));
    max_NS(i,:) = max(NS(i,1:348));
    Autocorr = autocorr(NS(i,1:348), 'NumLags', 30);
    Autocorr_NS_1Lag(i,:) = Autocorr(:,2);
    Autocorr_NS_12Lag(i,:) = Autocorr(:,13);
    Autocorr_NS_30Lag(i,:) = Autocorr(:,31);
end

descript_NS(:,1) = mean_NS;
descript_NS(:,2) = std_NS;
descript_NS(:,3) = min_NS;
descript_NS(:,4) = max_NS;
descript_NS(:,5) = Autocorr_NS_1Lag;
descript_NS(:,6) = Autocorr_NS_12Lag;
descript_NS(:,7) = Autocorr_NS_30Lag;
descript_NS(:,8) = [3:1:120]';;
descript_NS = array2table(descript_NS, 'VariableNames', ...
    {'Mean', 'StdDev', 'Min', 'Max', 'Rho1', 'Rho2', 'Rho3', 'Mat'} );

ind = ismember(table2array(descript_NS(:,8)), m_2020);
descript_NS = descript_NS(ind,:);

%% Actual Yield Descriptive Statistics 2015-2019

%We now turn to the actual observed yields 
clear descript_2019
data_2019 = table2array(data_2020(:,1:60));

for i=1:size(data_2019, 1)
    mean_2019(i,:) = mean(data_2019(i,:));
    std_2019(i,:) = std(data_2019(i,:));
    min_2019(i,:) = min(data_2019(i,:));
    max_2019(i,:) = max(data_2019(i,:));
    Autocorr = autocorr(data_2019(i,:), 'NumLags', 30);
    Autocorr_2019_1Lag(i,:) = Autocorr(:,2);
    Autocorr_2019_12Lag(i,:) = Autocorr(:,13);
    Autocorr_2019_30Lag(i,:) = Autocorr(:,31);
end

descript_2019(:,1) = mean_2019;
descript_2019(:,2) = std_2019;
descript_2019(:,3) = min_2019;
descript_2019(:,4) = max_2019;
descript_2019(:,5) = Autocorr_2019_1Lag;
descript_2019(:,6) = Autocorr_2019_12Lag;
descript_2019(:,7) = Autocorr_2019_30Lag;
descript_2019(:,8) = m_2020;
descript_2019 = array2table(descript_2019, 'VariableNames', ...
    {'Mean', 'StdDev', 'Min', 'Max', 'Rho1', 'Rho2', 'Rho3', 'Mat'} );

%Plot of average yield curve 2015-2019
plot(m_2020, mean_2020, '-o')
axis([0 130 1 2.5])
xlabel('Maturity in months')
ylabel('Yield in percent')

%% Actual Yield Descriptive Statistics 1991-2019
%We now turn to the actual observed yields 
clear descript 
data_2014_2 = data_2014;
data_2014_2(:,289) = m_2014;
ind = ismember(data_2014_2(:,289), m_2020);
data_full = [data_2014_2(ind,1:288) data_2019];

for i=1:size(data_full, 1)
    mean_full(i,:) = mean(data_full(i,:));
    std_full(i,:) = std(data_full(i,:));
    min_full(i,:) = min(data_full(i,:));
    max_full(i,:) = max(data_full(i,:));
    Autocorr = autocorr(data_full(i,:), 'NumLags', 30)
    Autocorr_full_1Lag(i,:) = Autocorr(:,2);
    Autocorr_full_12Lag(i,:) = Autocorr(:,13);
    Autocorr_full_30Lag(i,:) = Autocorr(:,31);
end

descript(:,1) = mean_full;
descript(:,2) = std_full;
descript(:,3) = min_full;
descript(:,4) = max_full;
descript(:,5) = Autocorr_full_1Lag;
descript(:,6) = Autocorr_full_12Lag;
descript(:,7) = Autocorr_full_30Lag;
descript(:,8) = m_2020;
descript = array2table(descript, 'VariableNames', ...
    {'Mean', 'StdDev', 'Min', 'Max', 'Rho1', 'Rho2', 'Rho3', 'Mat'} );

%Plot of average yield curve 1991-2019
plot(m_2020, mean_full, '-o')
axis([0 130 2.5 4.5])
xlabel('Maturity in months')
ylabel('Yield in percent')

function S
%% 3D plot of all historical yield curves 1991-2014: 
clear t1 t2 T Z X Xmesh Tmesh
t1 = datetime(1991,01,1,0,0,0);
t2 = datetime(2014,12,1,0,0,0);
T= t1:calmonths(1):t2;

Z = data_2014';
X = m_2014;
[Xmesh, Tmesh] = meshgrid(X,T);
Tmesh = datetime(Tmesh,'ConvertFrom','datenum');
figure
surf(Tmesh,Xmesh,Z);
zlabel('Yield in percent');
xlabel('Time');
ylabel('Maturity in months');

%% 3D plot of all historical yield curves 2015-2019: 
clear t1 t2 T Z X Xmesh Tmesh
t1 = datetime(2015,01,1,0,0,0);
t2 = datetime(2019,12,1,0,0,0);
T= t1:calmonths(1):t2;

Z = data_2019';
X = m_2020;
[Xmesh, Tmesh] = meshgrid(X,T);
Tmesh = datetime(Tmesh,'ConvertFrom','datenum');
figure
surf(Tmesh,Xmesh,Z);
zlabel('Yield in percent');
xlabel('Time');
ylabel('Maturity in months');

%% 3D plot of all historical yield curves 1991-2019: 
clear t1 t2 T Z X Xmesh Tmesh
t1 = datetime(1991,01,1,0,0,0);
t2 = datetime(2019,12,1,0,0,0);
T= t1:calmonths(1):t2;

Z = data_full';
X = m_2020;
[Xmesh, Tmesh] = meshgrid(X,T);
Tmesh = datetime(Tmesh,'ConvertFrom','datenum');
figure
surf(Tmesh,Xmesh,Z);
zlabel('Yield in percent');
xlabel('Time');
ylabel('Maturity in months');

%% Descriptive statistics betas
descript_betas = zeros(3,8);

%Temporary removal of the 349th obs (31.01.2020)
b1 = b1(1:348);
b2 = b2(1:348);
b3 = b3(1:348);

%First we find the average parameter values
avg_b1 = mean(b1);
avg_b2 = mean(b2);
avg_b3 = mean(b3);
descript_betas(:,1) = [avg_b1 avg_b2 avg_b3];

%std.dev of parameters
stdev_b1 = std(b1);
stdev_b2 = std(b2);
stdev_b3 = std(b3);
descript_betas(:,2) = [stdev_b1 stdev_b2 stdev_b3];

%Minimum of parameters
min_b1 = min(b1);
min_b2 = min(b2);
min_b3 = min(b3);
descript_betas(:,3) = [min_b1 min_b2 min_b3];

%Maximum of parameters
max_b1 = max(b1);
max_b2 = max(b2);
max_b3 = max(b3);
descript_betas(:,4) = [max_b1 max_b2 max_b3];

%Autocorrelations
Autocorr_b1 = autocorr(b1, 'NumLags', 30);
Autocorr_b1_1Lag = Autocorr_b1(2,:);
Autocorr_b1_12Lag = Autocorr_b1(13,:);
Autocorr_b1_30Lag = Autocorr_b1(31,:);

Autocorr_b2 = autocorr(b2, 'NumLags', 30);
Autocorr_b2_1Lag = Autocorr_b2(2,:);
Autocorr_b2_12Lag = Autocorr_b2(13,:);
Autocorr_b2_30Lag = Autocorr_b2(31,:);

Autocorr_b3 = autocorr(b3, 'NumLags', 30);
Autocorr_b3_1Lag = Autocorr_b3(2,:);
Autocorr_b3_12Lag = Autocorr_b3(13,:);
Autocorr_b3_30Lag = Autocorr_b3(31,:);

descript_betas(:,5) = [Autocorr_b1_1Lag Autocorr_b2_1Lag Autocorr_b3_1Lag];
descript_betas(:,6) = [Autocorr_b1_12Lag Autocorr_b2_12Lag Autocorr_b3_12Lag];
descript_betas(:,7) = [Autocorr_b1_30Lag Autocorr_b2_30Lag Autocorr_b3_30Lag];

%Augmented Dickey-Fuller (ADF) unit root test
[h,pValue] = adftest(b1);
Pval_b1= pValue;
[h,pValue] = adftest(b2);
Pval_b2= pValue;
[h,pValue] = adftest(b3);
Pval_b3= pValue;
descript_betas(:,8) = [Pval_b1 Pval_b2 Pval_b3];

descript_betas = array2table(descript_betas, 'RowNames', {'B1', 'B2', 'B3'},...
    'VariableNames', {'Mean', 'StdDev', 'Min', 'Max', 'Rho1', 'Rho12', 'Rho30', 'ADFpValue'} )


%% NS residual diagnostics
clear descript_res
%We obtain the residual by subtracting data_2014 from NS
NS_2 = NS(:,1:348);
NS_2(:,349) = m;
ind = ismember(NS_2(:,349), m_2020);
NS_2 = NS_2(ind,1:348);

residuals = data_full-NS_2;

for i=1:size(residuals, 1)
    mean_residual(i,:) = mean(residuals(i,:));
    std_residual(i,:) = std(residuals(i,:));
    min_residual(i,:) = min(residuals(i,:));
    max_residual(i,:) = max(residuals(i,:));
    MAE_residual(i,:) = mean(abs(residuals(i,:)));
    RMSE_residual(i,:) = sqrt(immse(data_full(i,:),...
        NS_2(i,:)));
    Autocorr = autocorr(residuals(i,:), 'NumLags', 30);
    Autocorr_res_1Lag(i,:) = Autocorr(:,2);
    Autocorr_res_12Lag(i,:) = Autocorr(:,13);
    Autocorr_res_30Lag(i,:) = Autocorr(:,31);
end

descript_res(:,1) = mean_residual;
descript_res(:,2) = std_residual;
descript_res(:,3) = min_residual;
descript_res(:,4) = max_residual;
descript_res(:,5) = MAE_residual;
descript_res(:,6) = RMSE_residual;
descript_res(:,7) = Autocorr_res_1Lag;
descript_res(:,8) = Autocorr_res_12Lag;
descript_res(:,9) = Autocorr_res_30Lag;
descript_res = array2table(descript_res, 'VariableNames', ...
    {'Mean', 'StdDev', 'Min', 'Max', 'MAE', 'RMSE', 'Rho1', 'Rho2', 'Rho3'} )

%Plot of NS vs Actual: 1991-2014
figure
plot(m, mean_NS_14, 'k')
hold on
sz = 30
scatter(m_2014, mean_2014, sz,'filled','d')
xlabel('Maturity in months')
ylabel('Yield in percent')
legend({'Nelson-Siegel', 'Actual'}, 'Location','southeast')
dim = [.72 .01 .3 .3];
str = 'R^2 = 0.933';
annotation('textbox',dim,'String',str,'FitBoxToText','on','EdgeColor', [1 1 1]);
axis([0 130 2.5 5])

%Plot of NS vs Actual: 2014-2019
figure
plot(m, mean_NS_19, 'k')
hold on
sz = 30
scatter(m_2020, mean_2019, sz,'filled','d')
xlabel('Maturity in months')
ylabel('Yield in percent')
legend({'Nelson-Siegel', 'Actual'}, 'Location','southeast')
dim = [.72 .01 .3 .3];
str = 'R^2 = 0.943';
annotation('textbox',dim,'String',str,'FitBoxToText','on','EdgeColor', [1 1 1]);
axis([0 130 1 2.4])

%% PCA for full sample
[coeff,score,latent] = pca(data_full');
PC1 = score(:,1);
PC2 = score(:,2);
PC3 = score(:,3);

l = data_full(8,:)';
s = (data_full(8,:) - data_full(1,:))';
c = (2.*data_full(4,:)-data_full(1,:)-data_full(8,:))';

corr(b1,b2)
corr(b1,b3)
corr(b2,b3)
%Beta1, PC1, and level
corr(PC1, b1)
corr(b1, l)
corr(PC1,l)

figure
plot(T, PC1+14, 'k')
hold on
plot(T, 3.*b1, 'b')
hold on
plot(T, 3.*l, 'r')
hold off
xlabel('Time')
legend({'PC_{1t}+14', '3\beta_{lt}', '3l_t'}, 'Location','northeast')
pbaspect([2.5 1 1])

%Beta2, PC2, slope
corr(PC2, b2)
corr(b2, s)
corr(PC2,s)

figure
plot(T, PC2, 'k')
hold on
plot(T, -b2, 'b')
hold on
plot(T, s, 'r')
hold off
xlabel('Time')
legend({'PC_{2t}', '-\beta_{2t}', 's_t'}, 'Location','southwest')
pbaspect([2.5 1 1])
ylim([-6.5 6])

%Beta3, PC3, curvature
corr(PC3, b3)
corr(b3, c)
corr(PC3,c)

figure
plot(T, -2.*PC3, 'k')
hold on
plot(T, 0.3.*b3, 'b')
hold on
plot(T, c, 'r')
hold off
xlabel('Time')
legend({'-2PC_{3t}', '0.3\beta_{3t}', 'c_t'}, 'Location','southwest')
pbaspect([2.5 1 1])
ylim([-5 2])

%% Part II: Out-of-sample forecasting 

%% Find factor changes
%Find monthly change in factors
b1_change = diff(b1)./b1(1:end-1,:);
b2_change = diff(b2)./b2(1:end-1,:);
b3_change = diff(b3)./b3(1:end-1,:);
%Note: First element of b_change is the change between 31.01.91 and 
%28.02.91, meaning the change during February

%Find 6-month change in factors
b1_change_6m = b1./lagmatrix(b1,6)-1;
b2_change_6m = b2./lagmatrix(b2,6)-1;
b3_change_6m = b3./lagmatrix(b3,6)-1;
b1_change_6m(isnan(b1_change_6m)) = [];
b2_change_6m(isnan(b2_change_6m)) = [];
b3_change_6m(isnan(b3_change_6m)) = [];
%Note: We lose the 6 first observations, such that the first obs. is  the
%change between 31.01.91 and 31.07.91.

%Find 12-month change in factors 
b1_change_12m = b1./lagmatrix(b1,12)-1;
b2_change_12m = b2./lagmatrix(b2,12)-1;
b3_change_12m = b3./lagmatrix(b3,12)-1;
b1_change_12m(isnan(b1_change_12m)) = [];
b2_change_12m(isnan(b2_change_12m)) = [];
b3_change_12m(isnan(b3_change_12m)) = [];
%Note: We lose the 12 first observations, such that the first obs. is the
%change between 31.01.91 and 31.01.92.

%% Importing regressors in the huge matrix X
clear X X_stat X_stat_1m X_stat_6m X_stat_12m
%Import and transform variables
X_main = readtable('VariableMasterfileRecursive3.xlsx','ReadVariableNames',true);
X = table2array(X_main(24:size(X_main,1),2:size(X_main,2)));

%We then have to transform the data to make it stationary
clear X_2 trans
for i=1:size(X,2)
    if adftest(X(:,i))==0 & sum((X(:,i))<=0)==0
            X_2(:,i) = log(X(:,i));%Then take the log of x_it
            trans(1,i) = "Log";
            trans(2,i) = 1;
            if adftest(X_2(:,i)) == 0%If still not stationary
                X_2(2:size(X_2,1),i) = diff(X_2(:,i));%Then take log first diff
                X_2(1,i) = NaN;%The first obs becomes NaN because of first diff
                trans(1,i) = "LogFirst";
                trans(2,i) = 2;
                if adftest(X_2(:,i))==0%If first diff is not stationary
                    X_2(2:size(X_2,1),i) = diff(X_2(:,i));%Then take second diff
                    trans(1,i) = "LogSecond";
                    trans(2,i) = 3;
                else %If first diff is stat, do nothing
                end
            else %If log stationary, do nothing
                X_2(:,i) = X_2(:,i);
            end
    elseif adftest(X(:,i))==0 & sum((X(:,i))<=0)~=0
                if sum(isinf(diff(X(:,i))./X(1:end-1,i)))==0
                    X_2(1,i) = NaN;    
                    X_2(2:size(X,1),i) = diff(X(:,i))./X(1:end-1,i);
                    trans(1,i) = "FirstInPercentage";
                    trans(2,i) = 4;
                else
                    X_2(1,i) = NaN;    
                    X_2(2:size(X,1),i) = diff(X(:,i));
                    trans(1,i) = "FirstInAbsolute";
                    trans(2,i) = 5;
                end
            if adftest(X_2(:,i)) == 0%If first diff is not stationary
                X_2(2:size(X_2,1),i) = diff(X_2(:,i));%Take second diff
                trans(1,i) = "SecondInAbsolute";
                trans(2,i) = 6;
            else %If first diff is stationary, do nothing
            end
    else%If original variable is stationary, do nothing
        X_2(:,i) = X(:,i);
        trans(1,i) = "Level(NoTrans)";
        trans(2,i) = 7;
    end
end

X_stat = X_2(3:size(X_2,1),:);

trans_nr = str2double(trans(2,:));
sum(trans_nr==1)%Number of log 
sum(trans_nr==2)%Number of log first
sum(trans_nr==3)%Number of log second
sum(trans_nr==4)%Number of first diff
sum(trans_nr==5)%Number of first diff abs
sum(trans_nr==6)%Number of second diff
sum(trans_nr==7)%Number of non-trans

%Control that transformation does not yield inf or nan
sum(sum(isnan(X_stat)))
sum(sum(isinf(X_stat)))

%Control that all transformed variables are stationary
clear h
for i=1:size(X_2,2)
    h(i) = adftest(X_2(:,i));
end
sum(h==0)
%All transformed variables are stationary! Great success. 

%Include squared terms. 
clear X_sq
X_sq = X_stat.^2;
%Concatenate
X_stat = [X_stat X_sq];

%The first observation in X is 01.01.1989 and the last is 01.12.2019. 
%By taking first differences, we find the change during January 1989 to the 
%change during November 2019. The change during January 89 is recorded as
%01.02.1989 and the change during November 19 is recorded as 01.12.2019.
%By taking second differences, we find the change in the changes; i.e. the
%change in the change during January 89 and the change during Febuary 89.
%This will be recorded as 01.03.1989. The change in the change during
%October and change during November 19 will be recorded as 01.12.2019. 
%The first observation in b1_change is the change between 31.01.91 and 
%28.02.91, i.e. the change during February! The last observation is the
%change during January 2020. Since we do not have the change during
%December 2019 for the regressors, we cannot forcast the change in beta1
%during January. The last observation in b1_change is hence deleted. 
b1_change(348)=[];
b2_change(348)=[];
b3_change(348)=[];

b1_change_6m(343)=[];
b2_change_6m(343)=[];
b3_change_6m(343)=[];

b1_change_12m(337)=[];
b2_change_12m(337)=[];
b3_change_12m(337)=[];

%For 1 month ahead:
%The first observation in X_stat is 01.03.1989. The first forecast 
%(in-sample) of b1 is going to be 01.02.1991. We need to extract the
%following sample: 01.02.1991 (1m prior to the 01.02.1991 obs in b1!) 
%to 01.02.2001. Row 24 of X_stat corresponds to 01.02.1991
X_stat_1m = X_stat(1:size(X_stat,1),:);

%For 6 months ahead:
%The first observation in X_stat is 01.03.1989. The first forecast 
%(in-sample) of b1 is going to be 01.07.1991 (Which means the change
%between 31.01.91 and 31.07.91). The 31.01.91 (01.02.91) value of X is
%going to forecast the first observation in b1_change_6m. 
%The last obs in b1_change_6m is the change between 31.06.2019 and
%31.12.2019. The 31.06.2019 (01.07.2019) value of X is going to forecast
%this. Row 24 of X corresponds to 01.02.91, and row 342 corresponds to
%01.07.2019. 
X_stat_6m = X_stat(1:342,:);

%For 12 months ahead:
%The first observation in X_stat is 01.03.1989. The first forecast 
%(in-sample) of b1 is going to be 01.01.1992 (Which means the change
%between 31.01.91 and 31.01.92). The 31.01.91 (01.02.91) value of X is
%going to forecast the first observation in b1_change_12m. 
%The last obs in b1_change_12m is the change between 31.12.2018 and
%31.12.2019. The 31.12.2018 (01.01.2019) value of X is going to forecast
%this. Row 24 of X corresponds to 01.02.91, and row 336 corresponds to
%01.01.2019. 
X_stat_12m = X_stat(1:336,:);

%Create vectors of dates to keep track of which row represents which dates
X_stat_dates = table2array(X_main(26:size(X_main,1),1));
X_stat_1m_dates = X_stat_dates(1:size(X_stat,1),:);
X_stat_6m_dates = X_stat_dates(1:342,:);
X_stat_12m_dates = X_stat_dates(1:336,:);

%We need lagged values of b1 (up to four lags) which we use as control
%variables.
clear b1_change_lagged b1_change_6m_lagged b1_change_12m_lagged ...
    b2_change_lagged b2_change_6m_lagged b2_change_12m_lagged ...
    b3_change_lagged b3_change_6m_lagged b3_change_12m_lagged
for i=1:4
b1_change_lagged(:,i) = lagmatrix(b1_change,i);
b2_change_lagged(:,i) = lagmatrix(b2_change,i);
b3_change_lagged(:,i) = lagmatrix(b3_change,i);
end
for i=1:4
b1_change_6m_lagged(:,i) = lagmatrix(b1_change_6m,i);
b2_change_6m_lagged(:,i) = lagmatrix(b2_change_6m,i);
b3_change_6m_lagged(:,i) = lagmatrix(b3_change_6m,i);
end
for i=1:4
b1_change_12m_lagged(:,i) = lagmatrix(b1_change_12m,i);
b2_change_12m_lagged(:,i) = lagmatrix(b2_change_12m,i);
b3_change_12m_lagged(:,i) = lagmatrix(b3_change_12m,i);
end

%% Out-of-sample, recursive forecasting using t-stat threshold
%% 1-month ahead forecasts: beta1 
%% Holdout-period: 2001:3 to 2019:12

clear b1_estimate T1 variable_count_b1_1m ind4
T1 = 122 %T1 is the first out-of-sample prediction (01.03.2001)
for i=(T1-1):(size(b1_change,1)-1)
X_stat_insample = X_stat_1m(1:i,:);%This will be the first predictor-matrix
%We now use the above sample to construct the forecasting model   
%All variables with a t-stat below alpha is discarded

%The first observation in b1_changes is the changes during February,
%reported as 01.02.1991. We need data on b1_changes until 01.02.2001.
b1_change_insample = b1_change(1:i);
%We need lagged values of b1 (up to four lags) which we use as control
%variables.
clear b1_change_lagged_insample
b1_change_lagged_insample = b1_change_lagged(1:i,:);

%We then run individual regressions of X on b1 changes controlling for lags
%of b1. Possible to initiate vectors to increase speed?
clear X_k tbl mdl t ind1 ind2 X_reduced_insample
t = zeros(size(X_stat_insample, 2),1);%initiate t vector
for k=1:size(X_stat_insample, 2);
    X_k = X_stat_insample(:,k);
    tbl = [b1_change_lagged_insample X_k];
    mdl = fitlm(tbl, b1_change_insample);
    t(k)= abs(mdl.Coefficients.tStat(6))';
end
%Hard thresholding based on ranking of t-stat
mean(t);
max(t);
min(t);
ind1 = find(t>=1.65);%aplha=5%
ind2 = t>=1.65;
X_reduced_insample = X_stat_insample(:,ind1);
variable_count_b1_1m(i-120,:) = ind2;

%PCA 
[coeff,score,latent] = pca(X_reduced_insample);
%We are only going to use the first 100  PC's since these will capture almost
%all variation in X_reduced
score = score(:,1:10);

%Use stepwiselm
clear Y
Y = [b1_change_lagged_insample score];
mdl = stepwiselm(Y, b1_change_insample, 'Upper', 'linear', 'Criterion', 'BIC');
%This model wil pick the best variables to minimize BIC
%Save R2 and BIC!

%Now we have to use the model using X_stat date for 01.03.2001 to forecast
%the change in b1 during March 2001.
clear ind3 X_reduced Y_2 Y_2_act coefficients 
ind3 = mdl.Formula.InModel;%This is an indicator of the variables used in the
%model.
%First, we need to find the targeted predictor values until 01.03.2001:
X_reduced = X_stat_1m(1:i+1,ind1);
[~,score,~] = pca(X_reduced);%We then calculate the updated PCs
score = score(:,1:10);%Select the first 10 PCs
%We then reconstruct Y_2 using data including 01.03.2001
Y_2 = [b1_change_lagged(1:i+1,:) score];
%We now wish to extract the variables choosen by stepwiselm for 01.03.2001:
Y_2_act = Y_2(i+1,ind3)';
%We need the coefficent from stepwiselm:
coefficients = mdl.Coefficients.Estimate;
b1_estimate(i-120) = (coefficients(1)+coefficients(2:size(coefficients,1))' ...
    *Y_2_act)'
end
b1_estimate = b1_estimate';
%FERDIG
%% 1-month ahead forecasts: beta2 
%% Holdout-period: 2001:3 to 2019:12

clear BIC R2 b2_estimate T1 variable_count_b2_1m ind4
T1 = 122 %T1 is the first out-of-sample prediction
for i=(T1-1):(size(b2_change,1)-1)
X_stat_insample = X_stat_1m(1:i,:);%This will be the first predictor-matrix
%We now use the above sample to construct the forecasting model   
%All variables with a t-stat below alpha is discarded

%The first observation in b1_changes is the changes during February,
%reported as 01.02.1991. We need data on b1_changes until 01.02.2001.
b2_change_insample = b2_change(1:i);
%We need lagged values of b1 (up to four lags) which we use as control
%variables.
clear b2_change_lagged_insample
b2_change_lagged_insample = b2_change_lagged(1:i,:);
%We then run individual regressions of X on b1 changes controlling for lags
%of b1. 
clear X_k tbl mdl t ind1 ind2 X_reduced_insample
t = zeros(size(X_stat_insample, 2),1);%initiate t vector
for k=1:size(X_stat_insample, 2);
    X_k = X_stat_insample(:,k);
    tbl = [b2_change_lagged_insample X_k];
    mdl = fitlm(tbl, b2_change_insample);
    t(k)= abs(mdl.Coefficients.tStat(6))';
end
%Hard thresholding based on ranking of t-stat
mean(t);
max(t);
min(t);
ind1 = find(t>=1.65);%aplha=5%
ind2 = t>=1.65;
X_reduced_insample = X_stat_insample(:,ind1);
variable_count_b2_1m(i-120,:) = ind2;

%PCA 
[coeff,score,latent] = pca(X_reduced_insample);
%We are only going to use the first 100  PC's since these will capture almost
%all variation in X_reduced
score = score(:,1:10);

%Model Alternative 1: Use stepwiselm
clear Y
Y = [b2_change_lagged_insample score];
mdl = stepwiselm(Y, b2_change_insample, 'Upper', 'linear', 'Criterion', 'BIC');
%This model wil pick the best variables to minimize BIC

%Now we have to use the model using X_stat date for 01.03.2001 to forecast
%the change in b1 during March 2001.
clear ind3 X_reduced Y_2 Y_2_act coefficients 
ind3 = mdl.Formula.InModel;%This is an indicator of the variables used in the
ind4(i-120,:) = ind3;%Save ind3 to see how often scores are included
%model.
%First, we need to find the targeted predictor values until 01.03.2001:
X_reduced = X_stat_1m(1:i+1,ind1);
[coeff,score,latent] = pca(X_reduced);%We then calculate the updated PCs
score = score(:,1:10);%Select the first 50 PCs
%We then reconstruct Y_2 using data including 01.03.2001
Y_2 = [b2_change_lagged(1:i+1,:) score];
%We now wish to extract the variables choosen by stepwiselm for 01.03.2001:
Y_2_act = Y_2(i+1,ind3)';
%We need the coefficent from stepwiselm:
coefficients = mdl.Coefficients.Estimate;
if sum(ind3) == 0 
    b2_estimate(i-120) = coefficients(1);
else
    b2_estimate(i-120) = (coefficients(1)+coefficients(2:size(coefficients,1))' ...
    *Y_2_act)'
end
end 
b2_estimate = b2_estimate';
% FERDIG
%% 1-month ahead forecasts: beta3
%% Holdout-period: 2001:3 to 2019:12

%The firs t observation in X_stat is 01.02.1991. We need to extract the
%following sample: 01.02.1991 to 01.02.2001. Row 1 of X_stat corresponds to
%01.02.1991 and row 121 correspond to 01.02.2001:

clear BIC R2 b3_estimate T1 variable_count_b3_1m
T1 = 122; %T1 is the first out-of-sample prediction
for i=(T1-1):(size(b3_change,1)-1)
X_stat_insample = X_stat_1m(1:i,:);%This will be the first predictor-matrix
%We now use the above sample to construct the forecasting model   
%All variables with a t-stat below alpha is discarded

%The first observation in b1_changes is the changes during February,
%reported as 01.02.1991. We need data on b1_changes until 01.02.2001.
b3_change_insample = b3_change(1:i);
%We need lagged values of b1 (up to four lags) which we use as control
%variables.
clear b3_change_lagged_insample
b3_change_lagged_insample = b3_change_lagged(1:i,:);
%We then run individual regressions of X on b1 changes controlling for lags
%of b1. 
clear X_k tbl mdl t ind1 ind2 X_reduced_insample
t = zeros(size(X_stat_insample, 2),1);%initiate t vector
for k=1:size(X_stat_insample, 2);
    X_k = X_stat_insample(:,k);
    tbl = [b3_change_lagged_insample X_k];
    mdl = fitlm(tbl, b3_change_insample);
    t(k)= abs(mdl.Coefficients.tStat(6))';
end
%Hard thresholding based on ranking of t-stat
mean(t);
max(t);
min(t);
ind1 = find(t>=1.65);%aplha=5%
ind_varcount = t>=1.65;
X_reduced_insample = X_stat_insample(:,ind1);
variable_count_b3_1m(i-120,:) = ind_varcount;

%PCA 
[coeff,score,latent] = pca(X_reduced_insample);
%We are only going to use the first 100  PC's since these will capture almost
%all variation in X_reduced
score = score(:,1:10);

%Model Alternative 1: Use stepwiselm
clear Y
Y = [b3_change_lagged_insample score];
mdl = stepwiselm(Y, b3_change_insample, 'Upper', 'linear', 'Criterion', 'BIC');
%This model wil pick the best variables to minimize BIC
%Save R2 and BIC!

%Now we have to use the model using X_stat date for 01.03.2001 to forecast
%the change in b1 during March 2001.
clear ind2 X_reduced Y_2 Y_2_act coefficients 
ind2 = mdl.Formula.InModel;%This is an indicator of the variables used in the
%model.
%First, we need to find the targeted predictor values until 01.03.2001:
X_reduced = X_stat_1m(1:i+1,ind1);
[coeff,score,latent] = pca(X_reduced);%We then calculate the updated PCs
score = score(:,1:10);%Select the first 50 PCs
%We then reconstruct Y_2 using data including 01.03.2001
Y_2 = [b3_change_lagged(1:i+1,:) score];
%We now wish to extract the variables choosen by stepwiselm for 01.03.2001:
Y_2_act = Y_2(i+1,ind2)';
%We need the coefficent from stepwiselm:
coefficients = mdl.Coefficients.Estimate;
if sum(ind2) == 0
    b3_estimate(i-120) = coefficients(1);
else
    b3_estimate(i-120) = (coefficients(1)+coefficients(2:size(coefficients,1))' ...
    *Y_2_act)
end
end

b3_estimate = b3_estimate';
% FERDIG
%% 6-months ahead forecasts: beta1 
%% Holdout-period: 2001:3 to 2019:12
function sixMbeta1
%The first in-sample forecast for b1_6m is 01.07.1991 (change between 01.01.91 and 01.07.91).
clear BIC R2 b1_6m_estimate T1 variable_count_b1_6m
T1 = 122; %Corresponds to 01.08.2001, which is the first out-of-sample 
%forecast of b1_changes_6m.
for i=(T1-1):(size(b1_change_6m,1)-1)%In sample period begins with row 1 to 121,
%and ends with row 341. At the last iteration (i=341), and an out-of-sample
%forecast for 342 is made. 
X_stat_insample = X_stat_6m(1:i,:);%This will be the first predictor-matrix
%We now use the above sample to construct the forecasting model   
%All variables with a t-stat below alpha is discarded.

%The first observation in b1_change_6m is the change from 31.01.91 to
%31.07.91, reported as 01.07.1991. We need data on b1_changes_6m until 01.07.2001.
b1_change_insample = b1_change_6m(1:i);
%We need lagged values of b1 (up to four lags) which we use as control
%variables.
clear b1_change_lagged_insample
b1_change_lagged_insample = b1_change_6m_lagged(1:i,:);
%We then run individual regressions of X on b1 changes controlling for lags
%of b1. 
clear X_k tbl mdl t ind1 ind2 X_reduced_insample
t = zeros(size(X_stat_insample, 2),1);%initiate t vector
for k=1:size(X_stat_insample, 2);
    X_k = X_stat_insample(:,k);
    tbl = [b1_change_lagged_insample X_k];
    mdl = fitlm(tbl, b1_change_insample);
    t(k)= abs(mdl.Coefficients.tStat(6))';
end
%Hard thresholding based on ranking of t-stat
mean(t);
max(t);
min(t);
ind1 = find(t>=1.65);%aplha=5%
ind2 = t>=1.65;
X_reduced_insample = X_stat_insample(:,ind1);
variable_count_b1_6m(i-120,:) = ind2;

%PCA 
[coeff,score,latent] = pca(X_reduced_insample);
%We are only going to use the first 50 PC's since these will capture almost
%all variation in X_reduced
score = score(:,1:10);

%Use stepwiselm
clear Y
Y = [b1_change_lagged_insample score];
mdl = stepwiselm(Y, b1_change_insample, 'Upper', 'linear', 'Criterion', 'BIC');
BIC(i-120) = mdl.ModelCriterion.BIC;
R2(i-120) = mdl.Rsquared.Ordinary;
%This model wil pick the best variables to minimize BIC
%Save R2 and BIC!

%Now we have to use the model using X_stat_6m until 01.03.2001 to forecast
%b1_change_6m on 01.08.2001 (which means the change between 31.02.2001 and
%31.08.2001)
clear ind2 X_reduced Y_2 Y_2_act coefficients 
ind2 = mdl.Formula.InModel;%This is an indicator of the variables used in the
%model.
%First, we need to find the targeted predictor values until 01.03.2001:
X_reduced = X_stat_6m(1:i+1,ind1);
[coeff,score,latent] = pca(X_reduced);%We then calculate the updated PCs
score = score(:,1:10);%Select the first 50 PCs
%We then reconstruct Y_2 using data including 01.03.2001
Y_2 = [b1_change_6m_lagged(1:i+1,:) score];
%We now wish to extract the variables choosen by stepwiselm for 01.03.2001:
Y_2_act = Y_2(i+1,ind2)';
%We need the coefficent from stepwiselm:
coefficients = mdl.Coefficients.Estimate;
if sum(ind2) == 0
    b1_6m_estimate(i-120) = coefficients(1);
else
    b1_6m_estimate(i-120) = (coefficients(1)+coefficients(2:size(coefficients,1))' ...
    *Y_2_act)
end
end

b1_6m_estimate = b1_6m_estimate';
% FERDIG!
%% 6-month ahead forecasts: beta2 
%% Holdout-period: 2001:3 to 2019:12

function sixMbeta2
%The first in-sample forecast for b2_6m is 01.07.1991 (change between 01.01.91 and 01.07.91).
clear BIC R2 b2_6m_estimate T1
T1 = 122 %Corresponds to 01.08.2001, which is the first out-of-sample 
%forecast of b2_changes_6m.
for i=(T1-1):(size(b2_change_6m,1)-1)%In sample period begins with row 1 to 121,
%and ends with row 341. At the last iteration (i=341), and an out-of-sample
%forecast for 342 is made. 
X_stat_insample = X_stat_6m(1:i,:);%This will be the first predictor-matrix
%We now use the above sample to construct the forecasting model   
%All variables with a t-stat below alpha is discarded.

%The first observation in b2_changes is the change from 31.01.91 to
%31.07.91, reported as 01.07.1991. We need data on b2_changes_6m until 01.07.2001.
b2_change_insample = b2_change_6m(1:i);
%We need lagged values of b2 (up to four lags) which we use as control
%variables.
clear b2_change_lagged_insample
b2_change_lagged_insample = b2_change_6m_lagged(1:i,:);
%We then run individual regressions of X on b1 changes controlling for lags
%of b2. 
clear X_k tbl mdl t ind1 ind2 X_reduced_insample
t = zeros(size(X_stat_insample, 2),1);%initiate t vector
for k=1:size(X_stat_insample, 2);
    X_k = X_stat_insample(:,k);
    tbl = [b2_change_lagged_insample X_k];
    mdl = fitlm(tbl, b2_change_insample);
    t(k)= abs(mdl.Coefficients.tStat(6))';
end
%Hard thresholding based on ranking of t-stat
mean(t);
max(t);
min(t);
ind1 = find(t>=1.65);%aplha=5%
ind2 = t>=1.65;
X_reduced_insample = X_stat_insample(:,ind1);
variable_count_b2_6m(i,:) = ind2;

%PCA 
[coeff,score,latent] = pca(X_reduced_insample);
%We are only going to use the first 50 PC's since these will capture almost
%all variation in X_reduced
score = score(:,1:10);

%Use stepwiselm
clear Y
Y = [b2_change_lagged_insample score];
mdl = stepwiselm(Y, b2_change_insample, 'Upper', 'linear', 'Criterion', 'BIC');
%This model wil pick the best variables to minimize BIC

%Now we have to use the model using X_stat_6m until 01.03.2001 to forecast
%b2_change_6m on 01.08.2001 (which means the change between 31.02.2001 and
%31.08.2001)
clear ind2 X_reduced Y_2 Y_2_act coefficients 
ind2 = mdl.Formula.InModel;%This is an indicator of the variables used in the
%model.
%First, we need to find the targeted predictor values until 01.03.2001:
X_reduced = X_stat_6m(1:i+1,ind1);
[coeff,score,latent] = pca(X_reduced);%We then calculate the updated PCs
score = score(:,1:10);%Select the first 50 PCs
%We then reconstruct Y_2 using data including 01.03.2001
Y_2 = [b2_change_6m_lagged(1:i+1,:) score];
%We now wish to extract the variables choosen by stepwiselm for 01.03.2001:
Y_2_act = Y_2(i+1,ind2)';
%We need the coefficent from stepwiselm:
coefficients = mdl.Coefficients.Estimate;
if sum(ind2) == 0
    b2_6m_estimate(i-120) = coefficients(1);
else
    b2_6m_estimate(i-120) = (coefficients(1)+coefficients(2:size(coefficients,1))' ...
    *Y_2_act)
end
end
b2_6m_estimate = b2_6m_estimate';


function sixMbeta3
%The first in-sample forecast for b3_6m is 01.07.1991 (change between 01.01.91 and 01.07.91).
clear BIC R2 b3_6m_estimate T1
T1 = 122 %Corresponds to 01.08.2001, which is the first out-of-sample 
%forecast of b3_changes_6m.
for i=(T1-1):(size(b3_change_6m,1)-1)%In sample period begins with row 1 to 121,
%and ends with row 341. At the last iteration (i=341), and an out-of-sample
%forecast for 342 is made. 
X_stat_insample = X_stat_6m(1:i,:);%This will be the first predictor-matrix
%We now use the above sample to construct the forecasting model   
%All variables with a t-stat below alpha is discarded.

%The first observation in b2_changes is the change from 31.01.91 to
%31.07.91, reported as 01.07.1991. We need data on b2_changes_6m until 01.07.2001.
b3_change_insample = b3_change_6m(1:i);
%We need lagged values of b3 (up to four lags) which we use as control
%variables.
clear b3_change_lagged_insample
b3_change_lagged_insample = b3_change_6m_lagged(1:i,:);
%We then run individual regressions of X on b1 changes controlling for lags
%of b3. 
clear X_k tbl mdl t ind1 ind2 X_reduced_insample
for k=1:size(X_stat_insample, 2);
    X_k = X_stat_insample(:,k);
    tbl = [b3_change_lagged_insample X_k];
    mdl = fitlm(tbl, b3_change_insample);
    t(k)= abs(mdl.Coefficients.tStat(6))';
end
%Hard thresholding based on ranking of t-stat
mean(t);
max(t);
min(t);
ind1 = find(t>=1.65);%aplha=5%
ind2 = t>=1.65;
X_reduced_insample = X_stat_insample(:,ind1);
variable_count_b3_6m(i,:) = ind2;

%PCA 
[coeff,score,latent] = pca(X_reduced_insample);
%We are only going to use the first 50 PC's since these will capture almost
%all variation in X_reduced
score = score(:,1:10);

%Use stepwiselm
clear Y
Y = [b3_change_lagged_insample score];
mdl = stepwiselm(Y, b3_change_insample, 'Upper', 'linear', 'Criterion', 'BIC');
BIC(i-120) = mdl.ModelCriterion.BIC;
R2(i-120) = mdl.Rsquared.Ordinary;
%This model wil pick the best variables to minimize BIC
%Save R2 and BIC!

%Now we have to use the model using X_stat_6m until 01.03.2001 to forecast
%b3_change_6m on 01.08.2001 (which means the change between 31.02.2001 and
%31.08.2001)
clear ind2 X_reduced Y_2 Y_2_act coefficients 
ind2 = mdl.Formula.InModel;%This is an indicator of the variables used in the
%model.
%First, we need to find the targeted predictor values until 01.03.2001:
X_reduced = X_stat_6m(1:i+1,ind1);
[coeff,score,latent] = pca(X_reduced);%We then calculate the updated PCs
score = score(:,1:10);%Select the first 50 PCs
%We then reconstruct Y_2 using data including 01.03.2001
Y_2 = [b3_change_6m_lagged(1:i+1,:) score];
%We now wish to extract the variables choosen by stepwiselm for 01.03.2001:
Y_2_act = Y_2(i+1,ind2)';
%We need the coefficent from stepwiselm:
coefficients = mdl.Coefficients.Estimate;
if sum(ind2) == 0
    b3_6m_estimate(i-120) = coefficients(1);
else
    b3_6m_estimate(i-120) = (coefficients(1)+coefficients(2:size(coefficients,1))' ...
    *Y_2_act)
end
end
b3_6m_estimate = b3_6m_estimate';
% KJØR!
%% 12-month ahead forecasts: beta1 
%% Holdout-period: 2001:3 to 2019:12
function twelveMbeta1
%The first in-sample forecast for b1_12m is 01.01.1992 (change between 31.01.91 and 31.01.92).
clear BIC R2 b1_12m_estimate T1
T1 = 122 %Corresponds to 01.02.2002, which is the first out-of-sample 
%forecast of b1_changes_12m.
for i=(T1-1):(size(b1_change_12m,1)-1)%In sample period begins with row 1 to 121,
%and ends with row 335. At the last iteration (i=335), and an out-of-sample
%forecast for 336 is made. 
X_stat_insample = X_stat_12m(1:i,:);%This will be the first predictor-matrix
%We now use the above sample to construct the forecasting model   
%All variables with a t-stat below alpha is discarded.

%The first observation in b1_changes is the change from 31.01.91 to
%31.01.92, reported as 01.01.1992. We need data on b1_changes_6m until 01.01.2002.
b1_change_insample = b1_change_12m(1:i);
%We need lagged values of b1 (up to four lags) which we use as control
%variables.
clear b1_change_lagged_insample
b1_change_lagged_insample = b1_change_12m_lagged(1:i,:);
%We then run individual regressions of X on b1 changes controlling for lags
%of b1. 
clear X_k tbl mdl t ind1 ind2 X_reduced_insample
t = zeros(size(X_stat_insample, 2),1);%initiate t vector
for k=1:size(X_stat_insample, 2);
    X_k = X_stat_insample(:,k);
    tbl = [b1_change_lagged_insample X_k];
    mdl = fitlm(tbl, b1_change_insample);
    t(k)= abs(mdl.Coefficients.tStat(6))';
end
%Hard thresholding based on ranking of t-stat
mean(t);
max(t);
min(t);
ind1 = find(t>=1.65);%aplha=5%
ind2 = t>=1.65;
X_reduced_insample = X_stat_insample(:,ind1);
variable_count_b1_12m(i,:) = ind2;

%PCA 
[coeff,score,latent] = pca(X_reduced_insample);
%We are only going to use the first 50 PC's since these will capture almost
%all variation in X_reduced
score = score(:,1:10);

%Use stepwiselm
clear Y
Y = [b1_change_lagged_insample score];
mdl = stepwiselm(Y, b1_change_insample, 'Upper', 'linear', 'Criterion', 'BIC');
BIC(i-120) = mdl.ModelCriterion.BIC;
R2(i-120) = mdl.Rsquared.Ordinary;
%This model wil pick the best variables to minimize BIC
%Save R2 and BIC!

%Now we have to use the model using X_stat_6m until 01.03.2001 to forecast
%b1_change_12m on 01.02.2002 (which means the change between 31.02.2001 and
%31.02.2002)
clear ind2 X_reduced Y_2 Y_2_act coefficients 
ind2 = mdl.Formula.InModel;%This is an indicator of the variables used in the
%model.
%First, we need to find the targeted predictor values until 01.03.2001:
X_reduced = X_stat_12m(1:i+1,ind1);
[coeff,score,latent] = pca(X_reduced);%We then calculate the updated PCs
score = score(:,1:10);%Select the first 50 PCs
%We then reconstruct Y_2 using data including 01.03.2001
Y_2 = [b1_change_12m_lagged(1:i+1,:) score];
%We now wish to extract the variables choosen by stepwiselm for 01.03.2001:
Y_2_act = Y_2(i+1,ind2)';
%We need the coefficent from stepwiselm:
coefficients = mdl.Coefficients.Estimate;
if sum(ind2) == 0
    b1_12m_estimate(i-120) = coefficients(1);
else
    b1_12m_estimate(i-120) = (coefficients(1)+coefficients(2:size(coefficients,1))' ...
    *Y_2_act)
end
end
b1_12m_estimate = b1_12m_estimate';
% KJØR!
%% 12-month ahead forecasts: beta2 
%% Holdout-period: 2001:3 to 2019:12
function twelveMbeta2
%The first in-sample forecast for b2_12m is 01.01.1992 (change between 31.01.91 and 31.01.92).
clear BIC R2 b2_12m_estimate T1
T1 = 122 %Corresponds to 01.02.2002, which is the first out-of-sample 
%forecast of b2_changes_12m.
for i=(T1-1):(size(b2_change_12m,1)-1)%In sample period begins with row 1 to 121,
%and ends with row 335. At the last iteration (i=335), and an out-of-sample
%forecast for 336 is made. 
X_stat_insample = X_stat_12m(1:i,:);%This will be the first predictor-matrix
%We now use the above sample to construct the forecasting model   
%All variables with a t-stat below alpha is discarded.

%The first observation in b2_changes is the change from 31.01.91 to
%31.01.92, reported as 01.01.1992. We need data on b2_changes_6m until 01.01.2002.
b2_change_insample = b2_change_12m(1:i);
%We need lagged values of b2 (up to four lags) which we use as control
%variables.
clear b2_change_lagged_insample
b2_change_lagged_insample = b2_change_12m_lagged(1:i,:);
%We then run individual regressions of X on b1 changes controlling for lags
%of b2. 
clear X_k tbl mdl t ind1 ind2 X_reduced_insample
t = zeros(size(X_stat_insample, 2),1);%initiate t vector
for k=1:size(X_stat_insample, 2);
    X_k = X_stat_insample(:,k);
    tbl = [b2_change_lagged_insample X_k];
    mdl = fitlm(tbl, b2_change_insample);
    t(k)= abs(mdl.Coefficients.tStat(6))';
end
%Hard thresholding based on ranking of t-stat
mean(t);
max(t);
min(t);
ind1 = find(t>=1.65);%aplha=5%
ind2 = t>=1.65;
X_reduced_insample = X_stat_insample(:,ind1);
variable_count_b2_12m(i,:) = ind2;

%PCA 
[coeff,score,latent] = pca(X_reduced_insample);
%We are only going to use the first 50 PC's since these will capture almost
%all variation in X_reduced
score = score(:,1:10);

%Use stepwiselm
clear Y
Y = [b2_change_lagged_insample score];
mdl = stepwiselm(Y, b2_change_insample, 'Upper', 'linear', 'Criterion', 'BIC');
BIC(i-120) = mdl.ModelCriterion.BIC;
R2(i-120) = mdl.Rsquared.Ordinary;
%This model wil pick the best variables to minimize BIC
%Save R2 and BIC!

%Now we have to use the model using X_stat_6m until 01.03.2001 to forecast
%b2_change_12m on 01.02.2002 (which means the change between 31.02.2001 and
%31.02.2002)
clear ind2 X_reduced Y_2 Y_2_act coefficients 
ind2 = mdl.Formula.InModel;%This is an indicator of the variables used in the
%model.
%First, we need to find the targeted predictor values until 01.03.2001:
X_reduced = X_stat_12m(1:i+1,ind1);
[coeff,score,latent] = pca(X_reduced);%We then calculate the updated PCs
score = score(:,1:50);%Select the first 50 PCs
%We then reconstruct Y_2 using data including 01.03.2001
Y_2 = [b2_change_12m_lagged(1:i+1,:) score];
%We now wish to extract the variables choosen by stepwiselm for 01.03.2001:
Y_2_act = Y_2(i+1,ind2)';
%We need the coefficent from stepwiselm:
coefficients = mdl.Coefficients.Estimate;
if sum(ind2) == 0
    b2_12m_estimate(i-120) = coefficients(1);
else
    b2_12m_estimate(i-120) = (coefficients(1)+coefficients(2:size(coefficients,1))' ...
    *Y_2_act)
end
end
b2_12m_estimate = b2_12m_estimate';
% KJØR!
%% 12-month ahead forecasts: beta3
%% Holdout-period: 2001:3 to 2019:12
function twelveMbeta3
%The first in-sample forecast for b3_12m is 01.01.1992 (change between 31.01.91 and 31.01.92).
clear BIC R2 b3_12m_estimate T1
T1 = 122 %Corresponds to 01.02.2002, which is the first out-of-sample 
%forecast of b3_changes_12m.
for i=(T1-1):(size(b3_change_12m,1)-1)%In sample period begins with row 1 to 121,
%and ends with row 335. At the last iteration (i=335), and an out-of-sample
%forecast for 336 is made. 
X_stat_insample = X_stat_12m(1:i,:);%This will be the first predictor-matrix
%We now use the above sample to construct the forecasting model   
%All variables with a t-stat below alpha is discarded.

%The first observation in b3_changes is the change from 31.01.91 to
%31.01.92, reported as 01.01.1992. We need data on b3_changes_6m until 01.01.2002.
b3_change_insample = b3_change_12m(1:i);
%We need lagged values of b3 (up to four lags) which we use as control
%variables.
clear b3_change_lagged_insample
b3_change_lagged_insample = b3_change_12m_lagged(1:i,:);
%We then run individual regressions of X on b1 changes controlling for lags
%of b3. 
clear X_k tbl mdl t ind1 ind2 X_reduced_insample
t = zeros(size(X_stat_insample, 2),1);%initiate t vector
for k=1:size(X_stat_insample, 2);
    X_k = X_stat_insample(:,k);
    tbl = [b3_change_lagged_insample X_k];
    mdl = fitlm(tbl, b3_change_insample);
    t(k)= abs(mdl.Coefficients.tStat(6))';
end
%Hard thresholding based on ranking of t-stat
mean(t);
max(t);
min(t);
ind1 = find(t>=1.65);%aplha=5%
ind2 = t>=1.65;
X_reduced_insample = X_stat_insample(:,ind1);
variable_count_b3_12m(i,:) = ind2;

%PCA 
[coeff,score,latent] = pca(X_reduced_insample);
%We are only going to use the first 50 PC's since these will capture almost
%all variation in X_reduced
score = score(:,1:10);

%Use stepwiselm
clear Y
Y = [b3_change_lagged_insample score];
mdl = stepwiselm(Y, b3_change_insample, 'Upper', 'linear', 'Criterion', 'BIC');
BIC(i-120) = mdl.ModelCriterion.BIC;
R2(i-120) = mdl.Rsquared.Ordinary;
%This model wil pick the best variables to minimize BIC
%Save R2 and BIC!

%Now we have to use the model using X_stat_6m until 01.03.2001 to forecast
%b3_change_12m on 01.02.2002 (which means the change between 31.02.2001 and
%31.02.2002)
clear ind2 X_reduced Y_2 Y_2_act coefficients 
ind2 = mdl.Formula.InModel;%This is an indicator of the variables used in the
%model.
%First, we need to find the targeted predictor values until 01.03.2001:
X_reduced = X_stat_12m(1:i+1,ind1);
[coeff,score,latent] = pca(X_reduced);%We then calculate the updated PCs
score = score(:,1:10);%Select the first 50 PCs
%We then reconstruct Y_2 using data including 01.03.2001
Y_2 = [b3_change_12m_lagged(1:i+1,:) score];
%We now wish to extract the variables choosen by stepwiselm for 01.03.2001:
Y_2_act = Y_2(i+1,ind2)';
%We need the coefficent from stepwiselm:
coefficients = mdl.Coefficients.Estimate;
if sum(ind2) == 0
    b3_12m_estimate(i-120) = coefficients(1);
else
    b3_12m_estimate(i-120) = (coefficients(1)+coefficients(2:size(coefficients,1))' ...
    *Y_2_act)
end
end
b3_12m_estimate = b3_12m_estimate';

%% 1-month ahead forecasts: Constructing forecasted yield curves
%Use the 1-month ahead forecasted values of b1, b2, and b3 to form model
%yield curves, and compare the forecasted and actual yields. 
clear RMSE MSE Yields_act forecast_dates ...
    b3_OoS b2_OoS b1_OoS b3_forecasted b2_forecasted b1_forecasted ...
    b3_tmin1 b2_tmin1 b1_tmin1 b_dates t2 t1
%Find the b1, b2 and b3 dates:
t1 = datetime(1991,1,31);
t2 = datetime(2019,12,31);
b_dates = (t1:calmonths:t2)';

%The first b1_changes forecast is the change between 31.02.2001 and
%31.03.2001 (reported as 01.03.2001). We need to multiply the 31.02.2001
%value of b1 with this growth factor in order to obtain the value of b1 on
%31.03.2001. We thus need to extract the b1 values from 31.02.2001
%(01.02.2001) to 31.11.2019 (01.11.2019). Row 347 of b1 corresponds to the
%31.11.2019 value of b1, while row 122 corresponds to 31.02.2001:
b1_tmin1 = b1(122:347);
b2_tmin1 = b2(122:347);
b3_tmin1 = b3(122:347);

clear b1_estimate b2_estimate b3_estimate b_estimate
b_estimate = readtable('betaEstimates1mNEW.xlsx','ReadVariableNames',false);
b1_estimate = table2array(b_estimate(:,1));
b2_estimate = table2array(b_estimate(:,2));
b3_estimate = table2array(b_estimate(:,3));

%We now find the forecasted factors:
b1_forecasted = b1_tmin1.*(1 + b1_estimate);
b2_forecasted = b2_tmin1.*(1 + b2_estimate);
b3_forecasted = b3_tmin1.*(1 + b3_estimate);

%We create a matrix of forecasted NS curves 
%whose rows correpsond to maturities and columns correspond to dates. 
%Column 1 thus represents the forecasted Nelson-Siegel curve on 01.03.2001.
clear m x1 x2 b1_ns b2_ns b3_ns NS_forecasted
m = [3 6 12 24 36 60 84 120]';
x1 = (1-exp(-m.*0.0609))./(m.*0.0609);
x2 = ((1-exp(-m.*0.0609))./(m.*0.0609))-exp(-m.*0.0609);
for i=1:size(b1_forecasted,1)
    b1_ns = b1_forecasted(i);
    b2_ns = b2_forecasted(i);
    b3_ns = b3_forecasted(i);
    NS_forecasted(:,i) = b1_ns + b2_ns.*x1 + b3_ns.*x2;
end

%We then need the actual yields for these samples: 
Yields_act = data_full(:,123:size(data_full,2));

%We can now caluclate performance statistics
for i=1:size(NS_forecasted,1)
    MSE(i) = (sum((Yields_act(i,:) - NS_forecasted(i,:))...
        .^2))./size(NS_forecasted,2);
    RMSE(i,:) = sqrt(MSE(i));
end

%% 6-month ahead forecasts: Constructing forecasted yield curves
%Use the 6-month ahead forecasted values of b1, b2, and b3 to form model
%yield curves, and compare the forecasted and actual yields. 
clear RMSE MSE Yields_act forecast_dates ...
    b3_OoS b2_OoS b1_OoS b3_forecasted b2_forecasted b1_forecasted ...
    b3_tmin1 b2_tmin1 b1_tmin1 b_dates t2 t1
%Find the b1, b2 and b3 dates:
t1 = datetime(1991,1,31);
t2 = datetime(2019,1,31);
b_dates = (t1:calmonths:t2)';

%The first b1_changes_6m forecast is the change between 31.02.2001 and
%31.08.2001 (reported as 01.08.2001). We need to multiply the 31.02.2001
%value of b1 with this growth factor in order to obtain the value of b1 on
%31.08.2001. We thus need to extract the b1 values from 31.02.2001
%(01.02.2001) to 31.06.2019 (01.06.2019). Row 342 of b1 corresponds to the
%31.06.2019 value of b1, while row 122 corresponds to 31.02.2001:

b1_tmin1 = b1(122:342);
b2_tmin1 = b2(122:342);
b3_tmin1 = b3(122:342);

clear b1_6m_estimate b2_6m_estimate b3_6m_estimate b_estimate
b_estimate = readtable('betaEstimates6m.xlsx','ReadVariableNames',false);
b1_6m_estimate = table2array(b_estimate(:,1));
b2_6m_estimate = table2array(b_estimate(:,2));
b3_6m_estimate = table2array(b_estimate(:,3));

%We now find the forecasted factors:
b1_forecasted = b1_tmin1.*(1 + b1_6m_estimate);
b2_forecasted = b2_tmin1.*(1 + b2_6m_estimate);
b3_forecasted = b3_tmin1.*(1 + b3_6m_estimate);

%We create a matrix of forecasted NS curves 
%whose rows correpsond to maturities and columns correspond to dates. 
%Column 1 thus represents the forecasted Nelson-Siegel curve on 01.03.2001.
clear m x1 x2 b1_ns b2_ns b3_ns NS_forecasted
m = [3 6 12 24 36 60 84 120]';
x1 = (1-exp(-m.*0.0609))./(m.*0.0609);
x2 = ((1-exp(-m.*0.0609))./(m.*0.0609))-exp(-m.*0.0609);
for i=1:size(b1_forecasted,1)
    b1_ns = b1_forecasted(i);
    b2_ns = b2_forecasted(i);
    b3_ns = b3_forecasted(i);
    NS_forecasted(:,i) = b1_ns + b2_ns.*x1 + b3_ns.*x2;
end

%We then need the actual yields for the sample: 
Yields_act = data_full(:,128:size(data_full,2));

%We can now caluclate performance statistics
for i=1:size(NS_forecasted,1)
    MSE(i) = (sum((Yields_act(i,:) - NS_forecasted(i,:))...
        .^2))./size(NS_forecasted,2);
    RMSE(i,:) = sqrt(MSE(i));
end

%% 12-month ahead forecasts: Constructing forecasted yield curves
%Use the 12-month ahead forecasted values of b1, b2, and b3 to form model
%yield curves, and compare the forecasted and actual yields. 
clear RMSE MSE Yields_act forecast_dates ...
    b3_OoS b2_OoS b1_OoS b3_forecasted b2_forecasted b1_forecasted ...
    b3_tmin1 b2_tmin1 b1_tmin1 b_dates t2 t1
%Find the b1, b2 and b3 dates:
t1 = datetime(1991,1,1);
t2 = datetime(2020,1,1);
b_dates = (t1:calmonths:t2)';

%The first b1_changes_6m forecast is the change between 31.02.2001 and
%31.02.2002 (reported as 01.02.2002). We need to multiply the 31.02.2001
%value of b1 with this growth factor in order to obtain the value of b1 on
%31.02.2001. We thus need to extract the b1 values from 31.02.2001
%(01.02.2001) to 31.12.2018 (01.12.2018). Row 336 of b1 corresponds to the
%31.06.2019 value of b1, while row 122 corresponds to 31.02.2001:

b1_tmin1 = b1(122:336);
b2_tmin1 = b2(122:336);
b3_tmin1 = b3(122:336);

clear b1_12m_estimate b2_12m_estimate b3_12m_estimate b_estimate
b_estimate = readtable('betaEstimates12m.xlsx','ReadVariableNames',false);
b1_12m_estimate = table2array(b_estimate(:,1));
b2_12m_estimate = table2array(b_estimate(:,2));
b3_12m_estimate = table2array(b_estimate(:,3));

%We now find the forecasted factors:
b1_forecasted = b1_tmin1.*(1 + b1_12m_estimate);
b2_forecasted = b2_tmin1.*(1 + b2_12m_estimate);
b3_forecasted = b3_tmin1.*(1 + b3_12m_estimate);

%We create a matrix of forecasted NS curves 
%whose rows correpsond to maturities and columns correspond to dates. 
%Column 1 thus represents the forecasted Nelson-Siegel curve on 01.03.2001.
clear m x1 x2 b1_ns b2_ns b3_ns NS_forecasted
m = [3 6 12 24 36 60 84 120]';
x1 = (1-exp(-m.*0.0609))./(m.*0.0609);
x2 = ((1-exp(-m.*0.0609))./(m.*0.0609))-exp(-m.*0.0609);
for i=1:size(b1_forecasted,1)
    b1_ns = b1_forecasted(i);
    b2_ns = b2_forecasted(i);
    b3_ns = b3_forecasted(i);
    NS_forecasted(:,i) = b1_ns + b2_ns.*x1 + b3_ns.*x2;
end

clear Yields_act
%We then need the actual yields for the sample: 
Yields_act = data_full(:,134:size(data_full,2));

%We can now caluclate performance statistics
for i=1:size(NS_forecasted,1)
    MSE(i) = (sum((Yields_act(i,:) - NS_forecasted(i,:))...
        .^2))./size(NS_forecasted,2);
    RMSE(i,:) = sqrt(MSE(i));
end

%% PART III: In-Sample Analysis

%% 1-month ahead forecasts: beta1 

%We then run individual regressions of X on b1 changes controlling for lags
%of b1.
clear X_k tbl mdl1 mdl2 t ind1 ind2 X_reduced
t = zeros(size(X_stat_1m, 2),1);%initiate t vector
for k=1:size(X_stat_1m, 2);
    X_k = X_stat_1m(:,k);
    tbl = [b1_change_lagged X_k];
    mdl = fitlm(tbl, b1_change);
    t(k)= abs(mdl.Coefficients.tStat(6))';
end
%Hard thresholding based on ranking of t-stat
mean(t);
median(t);
max(t);
min(t);
ind1 = find(t>=1.65);%aplha=5%
ind2 = t>=1.65;
X_reduced = X_stat_1m(:,ind1);
variable_count_b1_1m = ind2;

%Targeted predictors
clear score table_insample
[coeff,score,latent] = pca(X_reduced);
%We are only going to use the first 100  PC's since these will capture almost
%all variation in X_reduced
score = score(:,1:10);
clear Y
Y = [b1_change_lagged score];
%First we include PCs
mdl1 = stepwiselm(Y, b1_change, 'Upper', 'linear', 'Criterion', 'BIC');
table_insample(:,1) = cell2table(mdl1.Coefficients.Row);
table_insample(:,2) = array2table(mdl1.Coefficients.Estimate);
table_insample(:,3) = array2table(mdl1.Coefficients.tStat)
mdl1.ModelCriterion.BIC
mdl1.Rsquared.Ordinary
mdl1.Rsquared.Adjusted
%Then only lags 
clear table_insample mdl2
mdl2 = stepwiselm(b1_change_lagged, b1_change, 'Upper', 'linear', 'Criterion', 'BIC');
mdl2.ModelCriterion.BIC
table_insample(:,1) = cell2table(mdl2.Coefficients.Row);
table_insample(:,2) = array2table(mdl2.Coefficients.Estimate);
table_insample(:,3) = array2table(mdl2.Coefficients.tStat)
mdl2.ModelCriterion.BIC
mdl2.Rsquared.Ordinary
mdl2.Rsquared.Adjusted

%We regress each of the 10 PCs on b1_change
joint = fitlm(score, b1_change);
clear table
table(1,:) = joint.Coefficients.Estimate(2:11)';
table(2,:) = joint.Coefficients.tStat(2:11)';
joint.Rsquared.Ordinary
joint.Rsquared.Adjusted

%We regress each of the 10 non-targeted PCs on b1_change
clear score
[~,score,~] = pca(X_stat_1m);
score = score(:,1:10);
joint = fitlm(score, b1_change);
clear table
table(1,:) = joint.Coefficients.Estimate(2:11)';
table(2,:) = joint.Coefficients.tStat(2:11)';
joint.Rsquared.Ordinary
joint.Rsquared.Adjusted

%Not targeted predictors 
clear score table_insample
[coeff,score,latent] = pca(X_stat_1m);
score = score(:,1:10);
clear Y
Y = [b1_change_lagged score];
mdl = stepwiselm(Y, b1_change, 'Upper', 'linear', 'Criterion', 'BIC');
mdl.ModelCriterion.BIC
table_insample(:,1) = cell2table(mdl.Coefficients.Row);
table_insample(:,2) = array2table(mdl.Coefficients.Estimate);
table_insample(:,3) = array2table(mdl.Coefficients.tStat)

%Regress each of the 10 PCs on all explanatory variables
clear x_on_p x_on_p_R2
x_on_p_R2 = zeros(1196,10);
for i = 1:(size(X_stat_1m,2)./2)
    for j = 1:10
    x_on_p = fitlm(X_stat_1m(:,i),score(:,j));
    x_on_p_R2(i,j) = x_on_p.Rsquared.Ordinary;
    end
end

%% 1-month ahead forecasts: beta2 

%We then run individual regressions of X on b2 changes controlling for lags
%of b1.
clear X_k tbl mdl t ind1 ind2 X_reduced
t = zeros(size(X_stat_1m, 2),1);%initiate t vector
for k=1:size(X_stat_1m, 2);
    X_k = X_stat_1m(:,k);
    tbl = [b2_change_lagged X_k];
    mdl = fitlm(tbl, b2_change);
    t(k)= abs(mdl.Coefficients.tStat(6))';
end
%Hard thresholding based on ranking of t-stat
mean(t);
median(t);
max(t);
min(t);
ind1 = find(t>=1.65);%aplha=5%
ind2 = t>=1.65;
X_reduced = X_stat_1m(:,ind1);
variable_count_b2_1m = ind2;

%Targeted predictors
clear score table_insample
[coeff,score,latent] = pca(X_reduced);
%We are only going to use the first 100  PC's since these will capture almost
%all variation in X_reduced
score = score(:,1:10);
clear Y
Y = [b2_change_lagged score];
%First we include PCs
mdl1 = stepwiselm(Y, b2_change, 'Upper', 'linear', 'Criterion', 'BIC');
table_insample(:,1) = cell2table(mdl1.Coefficients.Row);
table_insample(:,2) = array2table(mdl1.Coefficients.Estimate);
table_insample(:,3) = array2table(mdl1.Coefficients.tStat)
mdl1.ModelCriterion.BIC
mdl1.Rsquared.Ordinary
mdl1.Rsquared.Adjusted
%Then only lags 
clear table_insample mdl2
mdl2 = stepwiselm(b2_change_lagged, b2_change, 'Upper', 'linear', 'Criterion', 'BIC');
mdl2.ModelCriterion.BIC
table_insample(:,1) = cell2table(mdl2.Coefficients.Row);
table_insample(:,2) = array2table(mdl2.Coefficients.Estimate);
table_insample(:,3) = array2table(mdl2.Coefficients.tStat)
mdl2.ModelCriterion.BIC
mdl2.Rsquared.Ordinary
mdl2.Rsquared.Adjusted

%We regress each of the 10 PCs on b2_change
joint = fitlm(score, b2_change);%All
clear table
table(1,:) = joint.Coefficients.Estimate(2:11)';
table(2,:) = joint.Coefficients.tStat(2:11)';
joint.Rsquared.Ordinary
joint.Rsquared.Adjusted

%We regress each of the 10 non-targeted PCs on b2_change
clear score
[~,score,~] = pca(X_stat_1m);
score = score(:,1:10);
joint = fitlm(score, b2_change);
clear table
table(1,:) = joint.Coefficients.Estimate(2:11)';
table(2,:) = joint.Coefficients.tStat(2:11)';
joint.Rsquared.Ordinary
joint.Rsquared.Adjusted


%Not targeted predictors
clear score table_insample
[coeff,score,latent] = pca(X_stat_1m);
score = score(:,1:10);
clear Y
Y = [b2_change_lagged score];
mdl = stepwiselm(Y, b2_change, 'Upper', 'linear', 'Criterion', 'BIC');
mdl.ModelCriterion.BIC
table_insample(:,1) = cell2table(mdl.Coefficients.Row);
table_insample(:,2) = array2table(mdl.Coefficients.Estimate);
table_insample(:,3) = array2table(mdl.Coefficients.tStat)

%Regress each of the 10 PCs on all explanatory variables
clear x_on_p x_on_p_R2
x_on_p_R2 = zeros(1196,10);
for i = 1:(size(X_stat_1m,2)./2)
    for j = 1:10
    x_on_p = fitlm(X_stat_1m(:,i),score(:,j));
    x_on_p_R2(i,j) = x_on_p.Rsquared.Ordinary;
    end
end

%% 1-month ahead forecasts: beta3 

%We then run individual regressions of X on b1 changes controlling for lags
%of b1.
clear X_k tbl mdl t ind1 ind2 X_reduced
t = zeros(size(X_stat_1m, 2),1);%initiate t vector
for k=1:size(X_stat_1m, 2);
    X_k = X_stat_1m(:,k);
    tbl = [b3_change_lagged X_k];
    mdl = fitlm(tbl, b3_change);
    t(k)= abs(mdl.Coefficients.tStat(6))';
end
%Hard thresholding based on ranking of t-stat
mean(t);
median(t);
max(t);
min(t);
ind1 = find(t>=1.65);%aplha=5%
ind2 = t>=1.65;
X_reduced = X_stat_1m(:,ind1);
variable_count_b3_1m = ind2;

%Targeted predictors
clear score table_insample
[coeff,score,latent] = pca(X_reduced);
%We are only going to use the first 100  PC's since these will capture almost
%all variation in X_reduced
score = score(:,1:10);
clear Y
Y = [b3_change_lagged score];
%First with PCs
mdl1 = stepwiselm(Y, b3_change, 'Upper', 'linear', 'Criterion', 'BIC');
table_insample(:,1) = cell2table(mdl1.Coefficients.Row);
table_insample(:,2) = array2table(mdl1.Coefficients.Estimate);
table_insample(:,3) = array2table(mdl1.Coefficients.tStat)
mdl1.ModelCriterion.BIC
mdl1.Rsquared.Ordinary
mdl1.Rsquared.Adjusted
%Then only lags 
clear table_insample mdl2
mdl2 = stepwiselm(b3_change_lagged, b3_change, 'Upper', 'linear', 'Criterion', 'BIC');
mdl2.ModelCriterion.BIC
table_insample(:,1) = cell2table(mdl2.Coefficients.Row);
table_insample(:,2) = array2table(mdl2.Coefficients.Estimate);
table_insample(:,3) = array2table(mdl2.Coefficients.tStat)
mdl2.ModelCriterion.BIC
mdl2.Rsquared.Ordinary
mdl2.Rsquared.Adjusted

%We regress each of the 10 PCs on b3_change
joint = fitlm(score, b3_change);
clear table
table(1,:) = joint.Coefficients.Estimate(2:11)';
table(2,:) = joint.Coefficients.tStat(2:11)';
joint.Rsquared.Ordinary
joint.Rsquared.Adjusted

%We regress each of the 10 non-targeted PCs on b3_change
clear score
[~,score,~] = pca(X_stat_1m);
score = score(:,1:10);
joint = fitlm(score, b3_change);
clear table
table(1,:) = joint.Coefficients.Estimate(2:11)';
table(2,:) = joint.Coefficients.tStat(2:11)';
joint.Rsquared.Ordinary
joint.Rsquared.Adjusted

%Not targeted predictors
clear score table_insample
[coeff,score,latent] = pca(X_stat_1m);
score = score(:,1:10);
clear Y
Y = [b1_change_lagged score];
mdl = stepwiselm(Y, b3_change, 'Upper', 'linear', 'Criterion', 'BIC');
mdl.ModelCriterion.BIC
table_insample(:,1) = cell2table(mdl.Coefficients.Row);
table_insample(:,2) = array2table(mdl.Coefficients.Estimate);
table_insample(:,3) = array2table(mdl.Coefficients.tStat)

%Regress each of the 10 PCs on all explanatory variables
clear x_on_p x_on_p_R2
x_on_p_R2 = zeros(1196,10);
for i = 1:(size(X_stat_1m,2)./2)
    for j = 1:10
    x_on_p = fitlm(X_stat_1m(:,i),score(:,j));
    x_on_p_R2(i,j) = x_on_p.Rsquared.Ordinary;
    end
end

%% 6-months ahead forecasts: beta1 

%We then run individual regressions of X on b1 changes controlling for lags
%of b1.
clear X_k tbl mdl t ind1 ind2 X_reduced
t = zeros(size(X_stat_6m, 2),1);%initiate t vector
for k=1:size(X_stat_6m, 2);
    X_k = X_stat_6m(:,k);
    tbl = [b1_change_6m_lagged X_k];
    mdl = fitlm(tbl, b1_change_6m);
    t(k)= abs(mdl.Coefficients.tStat(6))';
end
%Hard thresholding based on ranking of t-stat
mean(t);
median(t);
max(t);
min(t);
ind1 = find(t>=1.65);%aplha=5%
ind2 = t>=1.65;
X_reduced = X_stat_6m(:,ind1);
variable_count_b1_6m = ind2;

%Targeted predictors
clear score
[coeff,score,latent] = pca(X_reduced);
%We are only going to use the first 100  PC's since these will capture almost
%all variation in X_reduced
score = score(:,1:10);
clear Y table_insample
Y = [b1_change_6m_lagged score];
%First with PCs
mdl1 = stepwiselm(Y, b1_change_6m, 'Upper', 'linear', 'Criterion', 'BIC');
table_insample(:,1) = cell2table(mdl1.Coefficients.Row);
table_insample(:,2) = array2table(mdl1.Coefficients.Estimate);
table_insample(:,3) = array2table(mdl1.Coefficients.tStat)
mdl1.ModelCriterion.BIC
mdl1.Rsquared.Ordinary
mdl1.Rsquared.Adjusted
%Then only lags 
clear table_insample mdl2
mdl2 = stepwiselm(b1_change_6m_lagged, b1_change_6m, 'Upper', 'linear', 'Criterion', 'BIC');
mdl2.ModelCriterion.BIC
table_insample(:,1) = cell2table(mdl2.Coefficients.Row);
table_insample(:,2) = array2table(mdl2.Coefficients.Estimate);
table_insample(:,3) = array2table(mdl2.Coefficients.tStat)
mdl2.ModelCriterion.BIC
mdl2.Rsquared.Ordinary
mdl2.Rsquared.Adjusted

%We regress each of the 10 PCs on b1_change_6m
joint = fitlm(score, b1_change_6m);
clear table
table(1,:) = joint.Coefficients.Estimate(2:11)';
table(2,:) = joint.Coefficients.tStat(2:11)';
joint.Rsquared.Ordinary
joint.Rsquared.Adjusted

%We regress each of the 10 non-targeted PCs on b1_change_6m
clear score
[~,score,~] = pca(X_stat_6m);
score = score(:,1:10);
joint = fitlm(score, b1_change_6m);
clear table
table(1,:) = joint.Coefficients.Estimate(2:11)';
table(2,:) = joint.Coefficients.tStat(2:11)';
joint.Rsquared.Ordinary
joint.Rsquared.Adjusted

%Not targeted predictors 
clear score table_insample
[coeff,score,latent] = pca(X_stat_6m);
score = score(:,1:10);
clear Y
Y = [b1_change_6m_lagged score];
mdl = stepwiselm(Y, b1_change_6m, 'Upper', 'linear', 'Criterion', 'BIC');
mdl.ModelCriterion.BIC
table_insample(:,1) = cell2table(mdl.Coefficients.Row);
table_insample(:,2) = array2table(mdl.Coefficients.Estimate);
table_insample(:,3) = array2table(mdl.Coefficients.tStat)

%% 6-months ahead forecasts: beta2

%We then run individual regressions of X on b1 changes controlling for lags
%of b1.
clear X_k tbl mdl t ind1 ind2 X_reduced
t = zeros(size(X_stat_6m, 2),1);%initiate t vector
for k=1:size(X_stat_6m, 2);
    X_k = X_stat_6m(:,k);
    tbl = [b2_change_6m_lagged X_k];
    mdl = fitlm(tbl, b2_change_6m);
    t(k)= abs(mdl.Coefficients.tStat(6))';
end
%Hard thresholding based on ranking of t-stat
mean(t);
median(t);
max(t);
min(t);
ind1 = find(t>=1.65);%aplha=5%
ind2 = t>=1.65;
X_reduced = X_stat_6m(:,ind1);
variable_count_b2_6m = ind2;

%Targeted predictors
clear score table_insample
[coeff,score,latent] = pca(X_reduced);
%We are only going to use the first 100  PC's since these will capture almost
%all variation in X_reduced
score = score(:,1:10);
clear Y
Y = [b2_change_6m_lagged score];
%First with PCs
mdl1 = stepwiselm(Y, b2_change_6m, 'Upper', 'linear', 'Criterion', 'BIC');
table_insample(:,1) = cell2table(mdl1.Coefficients.Row);
table_insample(:,2) = array2table(mdl1.Coefficients.Estimate);
table_insample(:,3) = array2table(mdl1.Coefficients.tStat)
mdl1.ModelCriterion.BIC
mdl1.Rsquared.Ordinary
mdl1.Rsquared.Adjusted
%Then only lags 
clear table_insample mdl2
mdl2 = stepwiselm(b2_change_6m_lagged, b2_change_6m, 'Upper', 'linear', 'Criterion', 'BIC');
mdl2.ModelCriterion.BIC
table_insample(:,1) = cell2table(mdl2.Coefficients.Row);
table_insample(:,2) = array2table(mdl2.Coefficients.Estimate);
table_insample(:,3) = array2table(mdl2.Coefficients.tStat)
mdl2.ModelCriterion.BIC
mdl2.Rsquared.Ordinary
mdl2.Rsquared.Adjusted

%We regress each of the 10 PCs on b2_change_6m
joint = fitlm(score, b2_change_6m);
clear table
table(1,:) = joint.Coefficients.Estimate(2:11)';
table(2,:) = joint.Coefficients.tStat(2:11)';
joint.Rsquared.Ordinary
joint.Rsquared.Adjusted

%We regress each of the 10 non-targeted PCs on b2_change_6m
clear score
[~,score,~] = pca(X_stat_6m);
score = score(:,1:10);
joint = fitlm(score, b2_change_6m);
clear table
table(1,:) = joint.Coefficients.Estimate(2:11)';
table(2,:) = joint.Coefficients.tStat(2:11)';
joint.Rsquared.Ordinary
joint.Rsquared.Adjusted

%Not targeted predictors 
clear score table_insample
[coeff,score,latent] = pca(X_stat_6m);
score = score(:,1:10);
clear Y
Y = [b2_change_6m_lagged score];
mdl = stepwiselm(Y, b2_change_6m, 'Upper', 'linear', 'Criterion', 'BIC');
mdl.ModelCriterion.BIC
table_insample(:,1) = cell2table(mdl.Coefficients.Row);
table_insample(:,2) = array2table(mdl.Coefficients.Estimate);
table_insample(:,3) = array2table(mdl.Coefficients.tStat)

%% 6-months ahead forecasts: beta3

%We then run individual regressions of X on b1 changes controlling for lags
%of b1.
clear X_k tbl mdl t ind1 ind2 X_reduced
t = zeros(size(X_stat_6m, 2),1);%initiate t vector
for k=1:size(X_stat_6m, 2);
    X_k = X_stat_6m(:,k);
    tbl = [b3_change_6m_lagged X_k];
    mdl = fitlm(tbl, b3_change_6m);
    t(k)= abs(mdl.Coefficients.tStat(6))';
end
%Hard thresholding based on ranking of t-stat
mean(t);
median(t);
max(t);
min(t);
ind1 = find(t>=1.65);%aplha=5%
ind2 = t>=1.65;
X_reduced = X_stat_6m(:,ind1);
variable_count_b3_6m = ind2;

%Targeted predictors
clear score table_insample
[coeff,score,latent] = pca(X_reduced);
%We are only going to use the first 100  PC's since these will capture almost
%all variation in X_reduced
score = score(:,1:10);
clear Y
Y = [b3_change_6m_lagged score];
%First with PCs
mdl1 = stepwiselm(Y, b3_change_6m, 'Upper', 'linear', 'Criterion', 'BIC');
table_insample(:,1) = cell2table(mdl1.Coefficients.Row);
table_insample(:,2) = array2table(mdl1.Coefficients.Estimate);
table_insample(:,3) = array2table(mdl1.Coefficients.tStat)
mdl1.ModelCriterion.BIC
mdl1.Rsquared.Ordinary
mdl1.Rsquared.Adjusted
%Then only lags 
clear table_insample mdl2
mdl2 = stepwiselm(b3_change_6m_lagged, b3_change_6m, 'Upper', 'linear', 'Criterion', 'BIC');
mdl2.ModelCriterion.BIC
table_insample(:,1) = cell2table(mdl2.Coefficients.Row);
table_insample(:,2) = array2table(mdl2.Coefficients.Estimate);
table_insample(:,3) = array2table(mdl2.Coefficients.tStat)
mdl2.ModelCriterion.BIC
mdl2.Rsquared.Ordinary
mdl2.Rsquared.Adjusted

%We regress each of the 10 PCs on b3_change_6m
joint = fitlm(score, b3_change_6m);
clear table
table(1,:) = joint.Coefficients.Estimate(2:11)';
table(2,:) = joint.Coefficients.tStat(2:11)';
joint.Rsquared.Ordinary
joint.Rsquared.Adjusted

%We regress each of the 10 non-targeted PCs on b3_change_6m
clear score
[~,score,~] = pca(X_stat_6m);
score = score(:,1:10);
joint = fitlm(score, b3_change_6m);
clear table
table(1,:) = joint.Coefficients.Estimate(2:11)';
table(2,:) = joint.Coefficients.tStat(2:11)';
joint.Rsquared.Ordinary
joint.Rsquared.Adjusted

%Not targeted predictors 
clear score table_insample
[coeff,score,latent] = pca(X_stat_6m);
score = score(:,1:10);
clear Y
Y = [b3_change_6m_lagged score];
mdl = stepwiselm(Y, b3_change_6m, 'Upper', 'linear', 'Criterion', 'BIC');
mdl.ModelCriterion.BIC
table_insample(:,1) = cell2table(mdl.Coefficients.Row);
table_insample(:,2) = array2table(mdl.Coefficients.Estimate);
table_insample(:,3) = array2table(mdl.Coefficients.tStat)

%% 12-months ahead forecasts: beta1 

%We then run individual regressions of X on b1 changes controlling for lags
%of b1.
clear X_k tbl mdl t ind1 ind2 X_reduced
t = zeros(size(X_stat_12m, 2),1);%initiate t vector
for k=1:size(X_stat_12m, 2);
    X_k = X_stat_12m(:,k);
    tbl = [b1_change_12m_lagged X_k];
    mdl = fitlm(tbl, b1_change_12m);
    t(k)= abs(mdl.Coefficients.tStat(6))';
end
%Hard thresholding based on ranking of t-stat
mean(t);
median(t);
max(t);
min(t);
ind1 = find(t>=1.65);%aplha=5%
ind2 = t>=1.65;
X_reduced = X_stat_12m(:,ind1);
variable_count_b1_12m = ind2;

%Targeted predictors
clear score table_insample
[coeff,score,latent] = pca(X_reduced);
%We are only going to use the first 100  PC's since these will capture almost
%all variation in X_reduced
score = score(:,1:10);
clear Y
Y = [b1_change_12m_lagged score];
%First with PCs
mdl1 = stepwiselm(Y, b1_change_12m, 'Upper', 'linear', 'Criterion', 'BIC');
table_insample(:,1) = cell2table(mdl1.Coefficients.Row);
table_insample(:,2) = array2table(mdl1.Coefficients.Estimate);
table_insample(:,3) = array2table(mdl1.Coefficients.tStat)
mdl1.ModelCriterion.BIC
mdl1.Rsquared.Ordinary
mdl1.Rsquared.Adjusted
%Then only lags 
clear table_insample mdl2
mdl2 = stepwiselm(b1_change_12m_lagged, b1_change_12m, 'Upper', 'linear', 'Criterion', 'BIC');
mdl2.ModelCriterion.BIC
table_insample(:,1) = cell2table(mdl2.Coefficients.Row);
table_insample(:,2) = array2table(mdl2.Coefficients.Estimate);
table_insample(:,3) = array2table(mdl2.Coefficients.tStat)
mdl2.ModelCriterion.BIC
mdl2.Rsquared.Ordinary
mdl2.Rsquared.Adjusted

%We regress each of the 10 PCs on b1_change_12m
joint = fitlm(score, b1_change_12m);
clear table
table(1,:) = joint.Coefficients.Estimate(2:11)';
table(2,:) = joint.Coefficients.tStat(2:11)';
joint.Rsquared.Ordinary
joint.Rsquared.Adjusted

%We regress each of the 10 non-targeted PCs on b1_change_12m
clear score
[~,score,~] = pca(X_stat_12m);
score = score(:,1:10);
joint = fitlm(score, b1_change_12m);
clear table
table(1,:) = joint.Coefficients.Estimate(2:11)';
table(2,:) = joint.Coefficients.tStat(2:11)';
joint.Rsquared.Ordinary
joint.Rsquared.Adjusted

%Not targeted predictors 
clear score table_insample
[coeff,score,latent] = pca(X_stat_12m);
score = score(:,1:10);
clear Y
Y = [b1_change_12m_lagged score];
mdl = stepwiselm(Y, b1_change_12m, 'Upper', 'linear', 'Criterion', 'BIC');
mdl.ModelCriterion.BIC
table_insample(:,1) = cell2table(mdl.Coefficients.Row);
table_insample(:,2) = array2table(mdl.Coefficients.Estimate);
table_insample(:,3) = array2table(mdl.Coefficients.tStat)


%% 12-months ahead forecasts: beta2

%We then run individual regressions of X on b1 changes controlling for lags
%of b1.
clear X_k tbl mdl t ind1 ind2 X_reduced
t = zeros(size(X_stat_12m, 2),1);%initiate t vector
for k=1:size(X_stat_12m, 2);
    X_k = X_stat_12m(:,k);
    tbl = [b2_change_12m_lagged X_k];
    mdl = fitlm(tbl, b2_change_12m);
    t(k)= abs(mdl.Coefficients.tStat(6))';
end
%Hard thresholding based on ranking of t-stat
mean(t);
median(t);
max(t);
min(t);
ind1 = find(t>=1.65);%aplha=5%
ind2 = t>=1.65;
X_reduced = X_stat_12m(:,ind1);
variable_count_b2_12m = ind2;

%Targeted predictors
clear score table_insample
[coeff,score,latent] = pca(X_reduced);
%We are only going to use the first 100  PC's since these will capture almost
%all variation in X_reduced
score = score(:,1:10);
clear Y
Y = [b2_change_12m_lagged score];
mdl1 = stepwiselm(Y, b2_change_12m, 'Upper', 'linear', 'Criterion', 'BIC');
table_insample(:,1) = cell2table(mdl1.Coefficients.Row);
table_insample(:,2) = array2table(mdl1.Coefficients.Estimate);
table_insample(:,3) = array2table(mdl1.Coefficients.tStat)
mdl1.ModelCriterion.BIC
mdl1.Rsquared.Ordinary
mdl1.Rsquared.Adjusted
%Then only lags 
clear table_insample mdl2
mdl2 = stepwiselm(b2_change_12m_lagged, b2_change_12m, 'Upper', 'linear', 'Criterion', 'BIC');
mdl2.ModelCriterion.BIC
table_insample(:,1) = cell2table(mdl2.Coefficients.Row);
table_insample(:,2) = array2table(mdl2.Coefficients.Estimate);
table_insample(:,3) = array2table(mdl2.Coefficients.tStat)
mdl2.ModelCriterion.BIC
mdl2.Rsquared.Ordinary
mdl2.Rsquared.Adjusted

%We regress each of the 10 PCs on b2_change_12m
joint = fitlm(score, b2_change_12m);
clear table
table(1,:) = joint.Coefficients.Estimate(2:11)';
table(2,:) = joint.Coefficients.tStat(2:11)';
joint.Rsquared.Ordinary
joint.Rsquared.Adjusted

%We regress each of the 10 non-targeted PCs on b2_change_12m
clear score
[~,score,~] = pca(X_stat_12m);
score = score(:,1:10);
joint = fitlm(score, b2_change_12m);
clear table
table(1,:) = joint.Coefficients.Estimate(2:11)';
table(2,:) = joint.Coefficients.tStat(2:11)';
joint.Rsquared.Ordinary
joint.Rsquared.Adjusted

%Not targeted predictors 
clear score table_insample
[coeff,score,latent] = pca(X_stat_12m);
score = score(:,1:10);
clear Y
Y = [b2_change_12m_lagged score];
mdl = stepwiselm(Y, b2_change_12m, 'Upper', 'linear', 'Criterion', 'BIC');
mdl.ModelCriterion.BIC
table_insample(:,1) = cell2table(mdl.Coefficients.Row);
table_insample(:,2) = array2table(mdl.Coefficients.Estimate);
table_insample(:,3) = array2table(mdl.Coefficients.tStat)

%% 12-months ahead forecasts: beta3

%We then run individual regressions of X on b1 changes controlling for lags
%of b1.
clear X_k tbl mdl t ind1 ind2 X_reduced
t = zeros(size(X_stat_12m, 2),1);%initiate t vector
for k=1:size(X_stat_12m, 2);
    X_k = X_stat_12m(:,k);
    tbl = [b3_change_12m_lagged X_k];
    mdl = fitlm(tbl, b3_change_12m);
    t(k)= abs(mdl.Coefficients.tStat(6))';
end
%Hard thresholding based on ranking of t-stat
mean(t);
median(t);
max(t);
min(t);
ind1 = find(t>=1.65);%aplha=5%
ind2 = t>=1.65;
X_reduced = X_stat_12m(:,ind1);
variable_count_b3_12m = ind2;

%Targeted predictors
clear score table_insample
[coeff,score,latent] = pca(X_reduced);
%We are only going to use the first 100  PC's since these will capture almost
%all variation in X_reduced
score = score(:,1:10);
clear Y
Y = [b3_change_12m_lagged score];
mdl1 = stepwiselm(Y, b3_change_12m, 'Upper', 'linear', 'Criterion', 'BIC');
table_insample(:,1) = cell2table(mdl1.Coefficients.Row);
table_insample(:,2) = array2table(mdl1.Coefficients.Estimate);
table_insample(:,3) = array2table(mdl1.Coefficients.tStat)
mdl1.ModelCriterion.BIC
mdl1.Rsquared.Ordinary
mdl1.Rsquared.Adjusted
%Then only lags 
clear table_insample mdl2
mdl2 = stepwiselm(b3_change_12m_lagged, b3_change_12m, 'Upper', 'linear', 'Criterion', 'BIC');
mdl2.ModelCriterion.BIC
table_insample(:,1) = cell2table(mdl2.Coefficients.Row);
table_insample(:,2) = array2table(mdl2.Coefficients.Estimate);
table_insample(:,3) = array2table(mdl2.Coefficients.tStat)
mdl2.ModelCriterion.BIC
mdl2.Rsquared.Ordinary
mdl2.Rsquared.Adjusted

%We regress each of the 10 PCs on b3_change_12m
joint = fitlm(score, b3_change_12m);
clear table
table(1,:) = joint.Coefficients.Estimate(2:11)';
table(2,:) = joint.Coefficients.tStat(2:11)';
joint.Rsquared.Ordinary
joint.Rsquared.Adjusted

%We regress each of the 10 non-targeted PCs on b3_change_12m
clear score
[~,score,~] = pca(X_stat_12m);
score = score(:,1:10);
joint = fitlm(score, b3_change_12m);
clear table
table(1,:) = joint.Coefficients.Estimate(2:11)';
table(2,:) = joint.Coefficients.tStat(2:11)';
joint.Rsquared.Ordinary
joint.Rsquared.Adjusted

%Not targeted predictors 
clear score table_insample
[coeff,score,latent] = pca(X_stat_12m);
score = score(:,1:10);
clear Y
Y = [b3_change_12m_lagged score];
mdl = stepwiselm(Y, b3_change_12m, 'Upper', 'linear', 'Criterion', 'BIC');
mdl.ModelCriterion.BIC
table_insample(:,1) = cell2table(mdl.Coefficients.Row);
table_insample(:,2) = array2table(mdl.Coefficients.Estimate);
table_insample(:,3) = array2table(mdl.Coefficients.tStat)


%% Bars of Xit regressed on PCs for b1
clear X_on_p_b1
X_on_p_b1 = table2array(readtable('XonPCb1.xlsx','ReadVariableNames',false));
%The most important PCs for b1 1 month ahead are P1 and P6
%P1
bar(X_on_p_b1(:,2))
ylabel('R^2')
xline(9)
xline(64)
xline(142)
xline(239)
xline(417)
xline(456)
xline(471)
xline(491)
xline(495)
xline(505)
xline(563)
xline(624)
xline(715)
xline(816)
xline(874)
xline(937)
xline(1003)
xline(1024)
xline(1074)
xline(1101)
xline(1125)
xline(1139)
xline(1151)
xline(1158)
xline(1172)
pbaspect([2.5 1 1])

%P6 
bar(X_on_p_b1(:,7))
ylabel('R^2')
xline(9)
xline(64)
xline(142)
xline(239)
xline(417)
xline(456)
xline(471)
xline(491)
xline(495)
xline(505)
xline(563)
xline(624)
xline(715)
xline(816)
xline(874)
xline(937)
xline(1003)
xline(1024)
xline(1074)
xline(1101)
xline(1125)
xline(1139)
xline(1151)
xline(1158)
xline(1172)
pbaspect([2.5 1 1])

%% Bars of Xit regressed on PCs for b3
clear X_on_p_b3
X_on_p_b3 = table2array(readtable('XonPCb3.xlsx','ReadVariableNames',false));
%The most important PCs for b3 1 month ahead are P1 and P4
%P1
bar(X_on_p_b3(:,2))
ylabel('R^2')
xline(9)
xline(64)
xline(142)
xline(239)
xline(417)
xline(456)
xline(471)
xline(491)
xline(495)
xline(505)
xline(563)
xline(624)
xline(715)
xline(816)
xline(874)
xline(937)
xline(1003)
xline(1024)
xline(1074)
xline(1101)
xline(1125)
xline(1139)
xline(1151)
xline(1158)
xline(1172)
pbaspect([2.5 1 1])

%P6 
bar(X_on_p_b3(:,5))
ylabel('R^2')
xline(9)
xline(64)
xline(142)
xline(239)
xline(417)
xline(456)
xline(471)
xline(491)
xline(495)
xline(505)
xline(563)
xline(624)
xline(715)
xline(816)
xline(874)
xline(937)
xline(1003)
xline(1024)
xline(1074)
xline(1101)
xline(1125)
xline(1139)
xline(1151)
xline(1158)
xline(1172)
pbaspect([2.5 1 1])

%% Benchmarking Diebold-Li model

%% AR(1) 
%% 1-month ahead forecasts: beta1 
%% Holdout-period: 2001:3 to 2019:12
clear b1_estimate T1
T1 = 122 %T1 is the first out-of-sample prediction (01.03.2001)
for i=(T1-1):(size(b1_change,1)-1)

%The first observation in b1_changes is the changes during February,
%reported as 01.02.1991. We need data on b1_changes until 01.02.2001.
b1_change_insample = b1_change(1:i);
b1_change_lagged_insample = b1_change_lagged(1:i,1);
mdl = fitlm(b1_change_lagged_insample, b1_change_insample);

%Now we have to use the model using b1_change for 01.03.2001 to forecast
%the change in b1 during March 2001.
b1_estimate(i-120) = mdl.Coefficients.Estimate(1) + ...
                mdl.Coefficients.Estimate(2).*b1_change(i);
end
b1_estimate = b1_estimate';

%% 1-month ahead forecasts: beta2
%% Holdout-period: 2001:3 to 2019:12
clear b2_estimate T1
T1 = 122 %T1 is the first out-of-sample prediction (01.03.2001)
for i=(T1-1):(size(b2_change,1)-1)

%The first observation in b1_changes is the changes during February,
%reported as 01.02.1991. We need data on b1_changes until 01.02.2001.
b2_change_insample = b2_change(1:i);
b2_change_lagged_insample = b2_change_lagged(1:i,1);
mdl = fitlm(b2_change_lagged_insample, b2_change_insample);

%Now we have to use the model using b1_change for 01.03.2001 to forecast
%the change in b1 during March 2001.
b2_estimate(i-120) = mdl.Coefficients.Estimate(1) + ...
                mdl.Coefficients.Estimate(2).*b2_change(i);
end
b2_estimate = b2_estimate';

%% 1-month ahead forecasts: beta3
%% Holdout-period: 2001:3 to 2019:12
clear b3_estimate T1
T1 = 122; %T1 is the first out-of-sample prediction (01.03.2001)
for i=(T1-1):(size(b3_change,1)-1)

%The first observation in b1_changes is the changes during February,
%reported as 01.02.1991. We need data on b1_changes until 01.02.2001.
b3_change_insample = b3_change(1:i);
b3_change_lagged_insample = b3_change_lagged(1:i,1);
mdl = fitlm(b3_change_lagged_insample, b3_change_insample);

%Now we have to use the model using b1_change for 01.03.2001 to forecast
%the change in b1 during March 2001.
b3_estimate(i-120) = mdl.Coefficients.Estimate(1) + ...
                mdl.Coefficients.Estimate(2).*b3_change(i);
end
b3_estimate = b3_estimate';

%% 6-months ahead forecasts: beta1 
%% Holdout-period: 2001:3 to 2019:12
clear b1_6m_estimate T1
T1 = 122 %Corresponds to 01.08.2001, which is the first out-of-sample 
%forecast of b1_changes_6m.
for i=(T1-1):(size(b1_change_6m,1)-1)

%The first observation in b1_change_6m is the change from 31.01.91 to
%31.07.91, reported as 01.07.1991. We need data on b1_changes_6m until 01.07.2001.
b1_change_insample = b1_change_6m(1:i);
b1_change_6m_lagged_insample = b1_change_6m_lagged(1:i,1);
mdl = fitlm(b1_change_6m_lagged_insample, b1_change_insample);

%Now we have to use the model using b1_change for 01.03.2001 to forecast
%the change in b1 during March 2001.
b1_6m_estimate(i-120) = mdl.Coefficients.Estimate(1) + ...
                mdl.Coefficients.Estimate(2).*b1_change_6m(i);
end
b1_6m_estimate = b1_6m_estimate';

%% 6-months ahead forecasts: beta2
%% Holdout-period: 2001:3 to 2019:12
clear b2_estimate_6m T1
T1 = 122 %Corresponds to 01.08.2001, which is the first out-of-sample 
%forecast of b1_changes_6m.
for i=(T1-1):(size(b2_change_6m,1)-1)

%The first observation in b1_change_6m is the change from 31.01.91 to
%31.07.91, reported as 01.07.1991. We need data on b1_changes_6m until 01.07.2001.
b2_change_insample = b2_change_6m(1:i);
b2_change_6m_lagged_insample = b2_change_6m_lagged(1:i,1);
mdl = fitlm(b2_change_6m_lagged_insample, b2_change_insample);

%Now we have to use the model using b1_change for 01.03.2001 to forecast
%the change in b1 during March 2001.
b2_6m_estimate(i-120) = mdl.Coefficients.Estimate(1) + ...
                mdl.Coefficients.Estimate(2).*b2_change_6m(i);
end
b2_6m_estimate = b2_6m_estimate';

%% 6-month ahead forecasts: beta3
%% Holdout-period: 2001:3 to 2019:12
clear b3_estimate_6m T1
T1 = 122 %Corresponds to 01.08.2001, which is the first out-of-sample 
%forecast of b3_changes_6m.
for i=(T1-1):(size(b3_change_6m,1)-1)

%The first observation in b1_change_6m is the change from 31.01.91 to
%31.07.91, reported as 01.07.1991. We need data on b1_changes_6m until 01.07.2001.
b3_change_insample = b3_change_6m(1:i);
b3_change_6m_lagged_insample = b3_change_6m_lagged(1:i,1);
mdl = fitlm(b3_change_6m_lagged_insample, b3_change_insample);

%Now we have to use the model using b1_change for 01.03.2001 to forecast
%the change in b1 during March 2001.
b3_6m_estimate(i-120) = mdl.Coefficients.Estimate(1) + ...
                mdl.Coefficients.Estimate(2).*b3_change_6m(i);
end
b3_6m_estimate = b3_6m_estimate';

%% 12-month ahead forecasts: beta1 
%% Holdout-period: 2001:3 to 2019:12

clear b1_12m_estimate T1
T1 = 122 %Corresponds to 01.02.2002, which is the first out-of-sample 
%forecast of b1_changes_12m.
for i=(T1-1):(size(b1_change_12m,1)-1)

%The first observation in b1_change_6m is the change from 31.01.91 to
%31.07.91, reported as 01.07.1991. We need data on b1_changes_6m until 01.07.2001.
b1_change_insample = b1_change_12m(1:i);
b1_change_12m_lagged_insample = b1_change_12m_lagged(1:i,1);
mdl = fitlm(b1_change_12m_lagged_insample, b1_change_insample);
%Now we have to use the model using b1_change for 01.03.2001 to forecast
%the change in b1 during March 2001.
b1_12m_estimate(i-120) = mdl.Coefficients.Estimate(1) + ...
                mdl.Coefficients.Estimate(2).*b1_change_12m(i);
end
b1_12m_estimate = b1_12m_estimate';

%% 12-month ahead forecasts: beta2
%% Holdout-period: 2001:3 to 2019:12
clear b2_12m_estimate T1
T1 = 122 %Corresponds to 01.02.2002, which is the first out-of-sample 
%forecast of b1_changes_12m.
for i=(T1-1):(size(b2_change_12m,1)-1)

%The first observation in b1_change_6m is the change from 31.01.91 to
%31.07.91, reported as 01.07.1991. We need data on b2_changes_12m until 01.07.2001.
b2_change_insample = b2_change_12m(1:i);
b2_change_12m_lagged_insample = b2_change_12m_lagged(1:i,1);
mdl = fitlm(b2_change_12m_lagged_insample, b2_change_insample);
%Now we have to use the model using b1_change for 01.03.2001 to forecast
%the change in b1 during March 2001.
b2_12m_estimate(i-120) = mdl.Coefficients.Estimate(1) + ...
                mdl.Coefficients.Estimate(2).*b2_change_12m(i);
end
b2_12m_estimate = b2_12m_estimate';

%% 12-month ahead forecasts: beta3
%% Holdout-period: 2001:3 to 2019:12
clear b3_12m_estimate T1
T1 = 122 %Corresponds to 01.08.2001, which is the first out-of-sample 
%forecast of b3_changes_6m.
for i=(T1-1):(size(b3_change_12m,1)-1)

%The first observation in b1_change_6m is the change from 31.01.91 to
%31.07.91, reported as 01.07.1991. We need data on b1_changes_6m until 01.07.2001.
b3_change_insample = b3_change_12m(1:i);
b3_change_12m_lagged_insample = b3_change_12m_lagged(1:i,1);
mdl = fitlm(b3_change_12m_lagged_insample, b3_change_insample);
%Now we have to use the model using b1_change for 01.03.2001 to forecast
%the change in b1 during March 2001.
b3_12m_estimate(i-120) = mdl.Coefficients.Estimate(1) + ...
                mdl.Coefficients.Estimate(2).*b3_change_12m(i);
end
b3_12m_estimate = b3_12m_estimate';

%% 1-month ahead forecasts: Constructing forecasted yield curves
%Use the 1-month ahead forecasted values of b1, b2, and b3 to form model
%yield curves, and compare the forecasted and actual yields. 
clear RMSE MSE Yields_act forecast_dates ...
    b3_forecasted b2_forecasted b1_forecasted ...
    b3_tmin1 b2_tmin1 b1_tmin1 b_dates t2 t1
%Find the b1, b2 and b3 dates:
t1 = datetime(1991,1,31);
t2 = datetime(2019,12,31);
b_dates = (t1:calmonths:t2)';

%The first b1_changes forecast is the change between 31.02.2001 and
%31.03.2001 (reported as 01.03.2001). We need to multiply the 31.02.2001
%value of b1 with this growth factor in order to obtain the value of b1 on
%31.03.2001. We thus need to extract the b1 values from 31.02.2001
%(01.02.2001) to 31.11.2019 (01.11.2019). Row 347 of b1 corresponds to the
%31.11.2019 value of b1, while row 122 corresponds to 31.02.2001:
b1_tmin1 = b1(122:347);
b2_tmin1 = b2(122:347);
b3_tmin1 = b3(122:347);

%We now find the forecasted factors:
b1_forecasted = b1_tmin1.*(1 + b1_estimate);
b2_forecasted = b2_tmin1.*(1 + b2_estimate);
b3_forecasted = b3_tmin1.*(1 + b3_estimate);

%We create a matrix of forecasted NS curves 
%whose rows correpsond to maturities and columns correspond to dates. 
%Column 1 thus represents the forecasted Nelson-Siegel curve on 01.03.2001.
clear m x1 x2 b1_ns b2_ns b3_ns NS_forecasted
m = [3 6 12 24 36 60 84 120]';
x1 = (1-exp(-m.*0.0609))./(m.*0.0609);
x2 = ((1-exp(-m.*0.0609))./(m.*0.0609))-exp(-m.*0.0609);
for i=1:size(b1_forecasted,1)
    b1_ns = b1_forecasted(i);
    b2_ns = b2_forecasted(i);
    b3_ns = b3_forecasted(i);
    NS_forecasted(:,i) = b1_ns + b2_ns.*x1 + b3_ns.*x2;
end

%We then need the actual yields for these samples: 
Yields_act = data_full(:,123:size(data_full,2));

%We can now caluclate performance statistics
for i=1:size(NS_forecasted,1)
    MSE(i) = (sum((Yields_act(i,:) - NS_forecasted(i,:))...
        .^2))./size(NS_forecasted,2);
    RMSE(i,:) = sqrt(MSE(i));
end

%% 6-month ahead forecasts: Constructing forecasted yield curves
%Use the 6-month ahead forecasted values of b1, b2, and b3 to form model
%yield curves, and compare the forecasted and actual yields. 
clear RMSE MSE Yields_act forecast_dates ...
    b3_OoS b2_OoS b1_OoS b3_forecasted b2_forecasted b1_forecasted ...
    b3_tmin1 b2_tmin1 b1_tmin1 b_dates t2 t1
%Find the b1, b2 and b3 dates:
t1 = datetime(1991,1,31);
t2 = datetime(2019,1,31);
b_dates = (t1:calmonths:t2)';

%The first b1_changes_6m forecast is the change between 31.02.2001 and
%31.08.2001 (reported as 01.08.2001). We need to multiply the 31.02.2001
%value of b1 with this growth factor in order to obtain the value of b1 on
%31.08.2001. We thus need to extract the b1 values from 31.02.2001
%(01.02.2001) to 31.06.2019 (01.06.2019). Row 342 of b1 corresponds to the
%31.06.2019 value of b1, while row 122 corresponds to 31.02.2001:

b1_tmin1 = b1(122:342);
b2_tmin1 = b2(122:342);
b3_tmin1 = b3(122:342);

%We now find the forecasted factors:
b1_forecasted = b1_tmin1.*(1 + b1_6m_estimate);
b2_forecasted = b2_tmin1.*(1 + b2_6m_estimate);
b3_forecasted = b3_tmin1.*(1 + b3_6m_estimate);

%We create a matrix of forecasted NS curves 
%whose rows correpsond to maturities and columns correspond to dates. 
%Column 1 thus represents the forecasted Nelson-Siegel curve on 01.03.2001.
clear m x1 x2 b1_ns b2_ns b3_ns NS_forecasted
m = [3 6 12 24 36 60 84 120]';
x1 = (1-exp(-m.*0.0609))./(m.*0.0609);
x2 = ((1-exp(-m.*0.0609))./(m.*0.0609))-exp(-m.*0.0609);
for i=1:size(b1_forecasted,1)
    b1_ns = b1_forecasted(i);
    b2_ns = b2_forecasted(i);
    b3_ns = b3_forecasted(i);
    NS_forecasted(:,i) = b1_ns + b2_ns.*x1 + b3_ns.*x2;
end

%We then need the actual yields for the sample: 
Yields_act = data_full(:,128:size(data_full,2));

%We can now caluclate performance statistics
for i=1:size(NS_forecasted,1)
    MSE(i) = (sum((Yields_act(i,:) - NS_forecasted(i,:))...
        .^2))./size(NS_forecasted,2);
    RMSE(i,:) = sqrt(MSE(i));
end

%% 12-month ahead forecasts: Constructing forecasted yield curves
%Use the 12-month ahead forecasted values of b1, b2, and b3 to form model
%yield curves, and compare the forecasted and actual yields. 
clear RMSE MSE Yields_act forecast_dates ...
    b3_OoS b2_OoS b1_OoS b3_forecasted b2_forecasted b1_forecasted ...
    b3_tmin1 b2_tmin1 b1_tmin1 b_dates t2 t1

%The first b1_changes_6m forecast is the change between 31.02.2001 and
%31.02.2002 (reported as 01.02.2002). We need to multiply the 31.02.2001
%value of b1 with this growth factor in order to obtain the value of b1 on
%31.02.2001. We thus need to extract the b1 values from 31.02.2001
%(01.02.2001) to 31.12.2018 (01.12.2018). Row 336 of b1 corresponds to the
%31.06.2019 value of b1, while row 122 corresponds to 31.02.2001:

b1_tmin1 = b1(122:336);
b2_tmin1 = b2(122:336);
b3_tmin1 = b3(122:336);

%We now find the forecasted factors:
b1_forecasted = b1_tmin1.*(1 + b1_12m_estimate);
b2_forecasted = b2_tmin1.*(1 + b2_12m_estimate);
b3_forecasted = b3_tmin1.*(1 + b3_12m_estimate);

%We create a matrix of forecasted NS curves 
%whose rows correpsond to maturities and columns correspond to dates. 
%Column 1 thus represents the forecasted Nelson-Siegel curve on 01.03.2001.
clear m x1 x2 b1_ns b2_ns b3_ns NS_forecasted
m = [3 6 12 24 36 60 84 120]';
x1 = (1-exp(-m.*0.0609))./(m.*0.0609);
x2 = ((1-exp(-m.*0.0609))./(m.*0.0609))-exp(-m.*0.0609);
for i=1:size(b1_forecasted,1)
    b1_ns = b1_forecasted(i);
    b2_ns = b2_forecasted(i);
    b3_ns = b3_forecasted(i);
    NS_forecasted(:,i) = b1_ns + b2_ns.*x1 + b3_ns.*x2;
end

%We then need the actual yields for the sample: 
Yields_act = data_full(:,134:size(data_full,2));

%We can now caluclate performance statistics
for i=1:size(NS_forecasted,1)
    MSE(i) = (sum((Yields_act(i,:) - NS_forecasted(i,:))...
        .^2))./size(NS_forecasted,2);
    RMSE(i,:) = sqrt(MSE(i));
end

%% Benchmarking AR(p) with DNS
%% 1-month ahead forecasts: beta1 
%% Holdout-period: 2001:3 to 2019:12
clear b1_estimate T1
T1 = 122 %T1 is the first out-of-sample prediction (01.03.2001)
for i=(T1-1):(size(b1_change,1)-1)
%The first observation in b1_changes is the changes during February,
%reported as 01.02.1991. We need data on b1_changes until 01.02.2001.
b1_change_insample = b1_change(1:i);
b1_change_lagged_insample = b1_change_lagged(1:i,:);
mdl = stepwiselm(b1_change_lagged_insample, b1_change_insample, ...
    'Upper', 'linear', 'Criterion', 'BIC');

clear ind X
%Now we have to use the model using b1_change for 01.03.2001 to forecast
%the change in b1 during March 2001.
ind = mdl.Formula.InModel;%This is an indicator of the variables used in the
%model.
X = b1_change_lagged(i+1,ind)';%We get the updated regressor values for the lags
%picked by stepwiselm.
coefficients = mdl.Coefficients.Estimate;
if sum(ind) == 0 
    b1_estimate(i-120) = coefficients(1);
else
    b1_estimate(i-120) = (coefficients(1)+coefficients(2:size(coefficients,1))' ...
    *X)';
end
end
b1_estimate = b1_estimate';

%% 1-month ahead forecasts: beta2
%% Holdout-period: 2001:3 to 2019:12
clear b2_estimate T1
T1 = 122 %T1 is the first out-of-sample prediction (01.03.2001)
for i=(T1-1):(size(b2_change,1)-1)
%The first observation in b2_changes is the changes during February,
%reported as 01.02.1991. We need data on b1_changes until 01.02.2001.
b2_change_insample = b2_change(1:i);
b2_change_lagged_insample = b2_change_lagged(1:i,:);
mdl = stepwiselm(b2_change_lagged_insample, b2_change_insample, ...
    'Upper', 'linear', 'Criterion', 'BIC');

clear ind X
%Now we have to use the model using b1_change for 01.03.2001 to forecast
%the change in b1 during March 2001.
ind = mdl.Formula.InModel;%This is an indicator of the variables used in the
%model.
X = b2_change_lagged(i+1,ind)';%We get the updated regressor values for the lags
%picked by stepwiselm.
coefficients = mdl.Coefficients.Estimate;
if sum(ind) == 0 
    b2_estimate(i-120) = coefficients(1);
else
    b2_estimate(i-120) = (coefficients(1)+coefficients(2:size(coefficients,1))' ...
    *X)';
end
end
b2_estimate = b2_estimate';
%% 1-month ahead forecasts: beta3
%% Holdout-period: 2001:3 to 2019:12
clear b3_estimate T1
T1 = 122 %T1 is the first out-of-sample prediction (01.03.2001)
for i=(T1-1):(size(b3_change,1)-1)
%The first observation in b3_changes is the changes during February,
%reported as 01.02.1991. We need data on b1_changes until 01.02.2001.
b3_change_insample = b3_change(1:i);
b3_change_lagged_insample = b3_change_lagged(1:i,:);
mdl = stepwiselm(b3_change_lagged_insample, b3_change_insample, ...
    'Upper', 'linear', 'Criterion', 'BIC');

clear ind X
%Now we have to use the model using b1_change for 01.03.2001 to forecast
%the change in b1 during March 2001.
ind = mdl.Formula.InModel;%This is an indicator of the variables used in the
%model.
X = b3_change_lagged(i+1,ind)';%We get the updated regressor values for the lags
%picked by stepwiselm.
coefficients = mdl.Coefficients.Estimate;
if sum(ind) == 0 
    b3_estimate(i-120) = coefficients(1);
else
    b3_estimate(i-120) = (coefficients(1)+coefficients(2:size(coefficients,1))' ...
    *X)';
end
end
b3_estimate = b3_estimate';
%% 6-months ahead forecasts: beta1 
%% Holdout-period: 2001:3 to 2019:12

clear b1_6m_estimate T1
T1 = 122 %T1 is the first out-of-sample prediction (01.03.2001)
for i=(T1-1):(size(b1_change_6m,1)-1)
%The first observation in b1_changes is the changes during February,
%reported as 01.02.1991. We need data on b1_changes until 01.02.2001.
b1_change_insample = b1_change_6m(1:i);
b1_change_lagged_insample = b1_change_6m_lagged(1:i,:);
mdl = stepwiselm(b1_change_lagged_insample, b1_change_insample, ...
    'Upper', 'linear', 'Criterion', 'BIC');

clear ind X
%Now we have to use the model using b1_change for 01.03.2001 to forecast
%the change in b1 during March 2001.
ind = mdl.Formula.InModel;%This is an indicator of the variables used in the
%model.
X = b1_change_6m_lagged(i+1,ind)';%We get the updated regressor values for the lags
%picked by stepwiselm.
coefficients = mdl.Coefficients.Estimate;
if sum(ind) == 0 
    b1_6m_estimate(i-120) = coefficients(1);
else
    b1_6m_estimate(i-120) = (coefficients(1)+coefficients(2:size(coefficients,1))' ...
    *X)';
end
end
b1_6m_estimate = b1_6m_estimate';

%% 6-month ahead forecasts: beta2
%% Holdout-period: 2001:3 to 2019:12
clear b2_6m_estimate T1
T1 = 122 %T1 is the first out-of-sample prediction (01.03.2001)
for i=(T1-1):(size(b2_change_6m,1)-1)
%The first observation in b1_changes is the changes during February,
%reported as 01.02.1991. We need data on b1_changes until 01.02.2001.
b2_change_insample = b2_change_6m(1:i);
b2_change_lagged_insample = b2_change_6m_lagged(1:i,:);
mdl = stepwiselm(b2_change_lagged_insample, b2_change_insample, ...
    'Upper', 'linear', 'Criterion', 'BIC');

clear ind X
%Now we have to use the model using b1_change for 01.03.2001 to forecast
%the change in b1 during March 2001.
ind = mdl.Formula.InModel;%This is an indicator of the variables used in the
%model.
X = b2_change_6m_lagged(i+1,ind)';%We get the updated regressor values for the lags
%picked by stepwiselm.
coefficients = mdl.Coefficients.Estimate;
if sum(ind) == 0 
    b2_6m_estimate(i-120) = coefficients(1);
else
    b2_6m_estimate(i-120) = (coefficients(1)+coefficients(2:size(coefficients,1))' ...
    *X)';
end
end
b2_6m_estimate = b2_6m_estimate';

%% 6-month ahead forecasts: beta3
%% Holdout-period: 2001:3 to 2019:12
clear b3_6m_estimate T1
T1 = 122 %T1 is the first out-of-sample prediction (01.03.2001)
for i=(T1-1):(size(b3_change_6m,1)-1)
%The first observation in b1_changes is the changes during February,
%reported as 01.02.1991. We need data on b1_changes until 01.02.2001.
b3_change_insample = b3_change_6m(1:i);
b3_change_lagged_insample = b3_change_6m_lagged(1:i,:);
mdl = stepwiselm(b3_change_lagged_insample, b3_change_insample, ...
    'Upper', 'linear', 'Criterion', 'BIC');

clear ind X
%Now we have to use the model using b1_change for 01.03.2001 to forecast
%the change in b1 during March 2001.
ind = mdl.Formula.InModel;%This is an indicator of the variables used in the
%model.
X = b3_change_6m_lagged(i+1,ind)';%We get the updated regressor values for the lags
%picked by stepwiselm.
coefficients = mdl.Coefficients.Estimate;
if sum(ind) == 0 
    b3_6m_estimate(i-120) = coefficients(1);
else
    b3_6m_estimate(i-120) = (coefficients(1)+coefficients(2:size(coefficients,1))' ...
    *X)';
end
end
b3_6m_estimate = b3_6m_estimate';

%% 12-month ahead forecasts: beta1 
%% Holdout-period: 2001:3 to 2019:12

clear b1_12m_estimate T1
T1 = 122 %T1 is the first out-of-sample prediction (01.03.2001)
for i=(T1-1):(size(b1_change_12m,1)-1)
%The first observation in b1_changes is the changes during February,
%reported as 01.02.1991. We need data on b1_changes until 01.02.2001.
b1_change_insample = b1_change_12m(1:i);
b1_change_lagged_insample = b1_change_12m_lagged(1:i,:);
mdl = stepwiselm(b1_change_lagged_insample, b1_change_insample, ...
    'Upper', 'linear', 'Criterion', 'BIC');

clear ind X
%Now we have to use the model using b1_change for 01.03.2001 to forecast
%the change in b1 during March 2001.
ind = mdl.Formula.InModel;%This is an indicator of the variables used in the
%model.
X = b1_change_12m_lagged(i+1,ind)';%We get the updated regressor values for the lags
%picked by stepwiselm.
coefficients = mdl.Coefficients.Estimate;
if sum(ind) == 0 
    b1_12m_estimate(i-120) = coefficients(1);
else
    b1_12m_estimate(i-120) = (coefficients(1)+coefficients(2:size(coefficients,1))' ...
    *X)';
end
end
b1_12m_estimate = b1_12m_estimate';

%% 12-month ahead forecasts: beta2
%% Holdout-period: 2001:3 to 2019:12
clear b2_12m_estimate T1
T1 = 122 %T1 is the first out-of-sample prediction (01.03.2001)
for i=(T1-1):(size(b2_change_12m,1)-1)
%The first observation in b1_changes is the changes during February,
%reported as 01.02.1991. We need data on b1_changes until 01.02.2001.
b2_change_insample = b2_change_12m(1:i);
b2_change_lagged_insample = b2_change_12m_lagged(1:i,:);
mdl = stepwiselm(b2_change_lagged_insample, b2_change_insample, ...
    'Upper', 'linear', 'Criterion', 'BIC');

clear ind X
%Now we have to use the model using b1_change for 01.03.2001 to forecast
%the change in b1 during March 2001.
ind = mdl.Formula.InModel;%This is an indicator of the variables used in the
%model.
X = b2_change_12m_lagged(i+1,ind)';%We get the updated regressor values for the lags
%picked by stepwiselm.
coefficients = mdl.Coefficients.Estimate;
if sum(ind) == 0 
    b2_12m_estimate(i-120) = coefficients(1);
else
    b2_12m_estimate(i-120) = (coefficients(1)+coefficients(2:size(coefficients,1))' ...
    *X)';
end
end
b2_12m_estimate = b2_12m_estimate';

%% 12-month ahead forecasts: beta3
%% Holdout-period: 2001:3 to 2019:12
clear b3_12m_estimate T1
T1 = 122 %T1 is the first out-of-sample prediction (01.03.2001)
for i=(T1-1):(size(b3_change_12m,1)-1)
%The first observation in b1_changes is the changes during February,
%reported as 01.02.1991. We need data on b1_changes until 01.02.2001.
b3_change_insample = b3_change_12m(1:i);
b3_change_lagged_insample = b3_change_12m_lagged(1:i,:);
mdl = stepwiselm(b3_change_lagged_insample, b3_change_insample, ...
    'Upper', 'linear', 'Criterion', 'BIC');

clear ind X
%Now we have to use the model using b1_change for 01.03.2001 to forecast
%the change in b1 during March 2001.
ind = mdl.Formula.InModel;%This is an indicator of the variables used in the
%model.
X = b3_change_12m_lagged(i+1,ind)';%We get the updated regressor values for the lags
%picked by stepwiselm.
coefficients = mdl.Coefficients.Estimate;
if sum(ind) == 0 
    b3_12m_estimate(i-120) = coefficients(1);
else
    b3_12m_estimate(i-120) = (coefficients(1)+coefficients(2:size(coefficients,1))' ...
    *X)';
end
end
b3_12m_estimate = b3_12m_estimate';

%% 1-month ahead forecasts: Constructing forecasted yield curves
%Use the 1-month ahead forecasted values of b1, b2, and b3 to form model
%yield curves, and compare the forecasted and actual yields. 
clear RMSE MSE Yields_act forecast_dates ...
    b3_OoS b2_OoS b1_OoS b3_forecasted b2_forecasted b1_forecasted ...
    b3_tmin1 b2_tmin1 b1_tmin1 b_dates t2 t1
%Find the b1, b2 and b3 dates:
t1 = datetime(1991,1,31);
t2 = datetime(2019,12,31);
b_dates = (t1:calmonths:t2)';

%The first b1_changes forecast is the change between 31.02.2001 and
%31.03.2001 (reported as 01.03.2001). We need to multiply the 31.02.2001
%value of b1 with this growth factor in order to obtain the value of b1 on
%31.03.2001. We thus need to extract the b1 values from 31.02.2001
%(01.02.2001) to 31.11.2019 (01.11.2019). Row 347 of b1 corresponds to the
%31.11.2019 value of b1, while row 122 corresponds to 31.02.2001:
b1_tmin1 = b1(122:347);
b2_tmin1 = b2(122:347);
b3_tmin1 = b3(122:347);

%We now find the forecasted factors:
b1_forecasted = b1_tmin1.*(1 + b1_estimate);
b2_forecasted = b2_tmin1.*(1 + b2_estimate);
b3_forecasted = b3_tmin1.*(1 + b3_estimate);

%We create a matrix of forecasted NS curves 
%whose rows correpsond to maturities and columns correspond to dates. 
%Column 1 thus represents the forecasted Nelson-Siegel curve on 01.03.2001.
clear m x1 x2 b1_ns b2_ns b3_ns NS_forecasted
m = [3 6 12 24 36 60 84 120]';
x1 = (1-exp(-m.*0.0609))./(m.*0.0609);
x2 = ((1-exp(-m.*0.0609))./(m.*0.0609))-exp(-m.*0.0609);
for i=1:size(b1_forecasted,1)
    b1_ns = b1_forecasted(i);
    b2_ns = b2_forecasted(i);
    b3_ns = b3_forecasted(i);
    NS_forecasted(:,i) = b1_ns + b2_ns.*x1 + b3_ns.*x2;
end

%Find the forecasted dates:
t1 = datetime(2001,3,1);
t2 = datetime(2019,12,1);
forecast_dates = (t1:calmonths:t2)';

%We then need the actual yields for these samples: 
Yields_act = data_full(:,123:size(data_full,2));

%We can now caluclate performance statistics
for i=1:size(NS_forecasted,1)
    MSE(i) = (sum((Yields_act(i,:) - NS_forecasted(i,:))...
        .^2))./size(NS_forecasted,2);
    RMSE(i,:) = sqrt(MSE(i));
end

%% 6-month ahead forecasts: Constructing forecasted yield curves
%Use the 6-month ahead forecasted values of b1, b2, and b3 to form model
%yield curves, and compare the forecasted and actual yields. 
clear RMSE MSE Yields_act forecast_dates ...
    b3_OoS b2_OoS b1_OoS b3_forecasted b2_forecasted b1_forecasted ...
    b3_tmin1 b2_tmin1 b1_tmin1 b_dates t2 t1
%Find the b1, b2 and b3 dates:
t1 = datetime(1991,1,31);
t2 = datetime(2019,1,31);
b_dates = (t1:calmonths:t2)';

%The first b1_changes_6m forecast is the change between 31.02.2001 and
%31.08.2001 (reported as 01.08.2001). We need to multiply the 31.02.2001
%value of b1 with this growth factor in order to obtain the value of b1 on
%31.08.2001. We thus need to extract the b1 values from 31.02.2001
%(01.02.2001) to 31.06.2019 (01.06.2019). Row 342 of b1 corresponds to the
%31.06.2019 value of b1, while row 122 corresponds to 31.02.2001:

b1_tmin1 = b1(122:342);
b2_tmin1 = b2(122:342);
b3_tmin1 = b3(122:342);

%We now find the forecasted factors:
b1_forecasted = b1_tmin1.*(1 + b1_6m_estimate);
b2_forecasted = b2_tmin1.*(1 + b2_6m_estimate);
b3_forecasted = b3_tmin1.*(1 + b3_6m_estimate);

%We create a matrix of forecasted NS curves 
%whose rows correpsond to maturities and columns correspond to dates. 
%Column 1 thus represents the forecasted Nelson-Siegel curve on 01.03.2001.
clear m x1 x2 b1_ns b2_ns b3_ns NS_forecasted
m = [3 6 12 24 36 60 84 120]';
x1 = (1-exp(-m.*0.0609))./(m.*0.0609);
x2 = ((1-exp(-m.*0.0609))./(m.*0.0609))-exp(-m.*0.0609);
for i=1:size(b1_forecasted,1)
    b1_ns = b1_forecasted(i);
    b2_ns = b2_forecasted(i);
    b3_ns = b3_forecasted(i);
    NS_forecasted(:,i) = b1_ns + b2_ns.*x1 + b3_ns.*x2;
end

%We then need the actual yields for the sample: 
Yields_act = data_full(:,128:size(data_full,2));

%We can now caluclate performance statistics
for i=1:size(NS_forecasted,1)
    MSE(i) = (sum((Yields_act(i,:) - NS_forecasted(i,:))...
        .^2))./size(NS_forecasted,2);
    RMSE(i,:) = sqrt(MSE(i));
end

%% 12-month ahead forecasts: Constructing forecasted yield curves
%Use the 12-month ahead forecasted values of b1, b2, and b3 to form model
%yield curves, and compare the forecasted and actual yields. 
clear RMSE MSE Yields_act forecast_dates ...
    b3_OoS b2_OoS b1_OoS b3_forecasted b2_forecasted b1_forecasted ...
    b3_tmin1 b2_tmin1 b1_tmin1 b_dates t2 t1
%Find the b1, b2 and b3 dates:
t1 = datetime(1991,1,1);
t2 = datetime(2020,1,1);
b_dates = (t1:calmonths:t2)';

%The first b1_changes_6m forecast is the change between 31.02.2001 and
%31.02.2002 (reported as 01.02.2002). We need to multiply the 31.02.2001
%value of b1 with this growth factor in order to obtain the value of b1 on
%31.02.2001. We thus need to extract the b1 values from 31.02.2001
%(01.02.2001) to 31.12.2018 (01.12.2018). Row 336 of b1 corresponds to the
%31.06.2019 value of b1, while row 122 corresponds to 31.02.2001:

b1_tmin1 = b1(122:336);
b2_tmin1 = b2(122:336);
b3_tmin1 = b3(122:336);

%We now find the forecasted factors:
b1_forecasted = b1_tmin1.*(1 + b1_12m_estimate);
b2_forecasted = b2_tmin1.*(1 + b2_12m_estimate);
b3_forecasted = b3_tmin1.*(1 + b3_12m_estimate);

%We create a matrix of forecasted NS curves 
%whose rows correpsond to maturities and columns correspond to dates. 
%Column 1 thus represents the forecasted Nelson-Siegel curve on 01.03.2001.
clear m x1 x2 b1_ns b2_ns b3_ns NS_forecasted
m = [3 6 12 24 36 60 84 120]';
x1 = (1-exp(-m.*0.0609))./(m.*0.0609);
x2 = ((1-exp(-m.*0.0609))./(m.*0.0609))-exp(-m.*0.0609);
for i=1:size(b1_forecasted,1)
    b1_ns = b1_forecasted(i);
    b2_ns = b2_forecasted(i);
    b3_ns = b3_forecasted(i);
    NS_forecasted(:,i) = b1_ns + b2_ns.*x1 + b3_ns.*x2;
end

%We then need the actual yields for the sample: 
Yields_act = data_full(:,134:size(data_full,2));

%We can now caluclate performance statistics
for i=1:size(NS_forecasted,1)
    MSE(i) = (sum((Yields_act(i,:) - NS_forecasted(i,:))...
        .^2))./size(NS_forecasted,2);
    RMSE(i,:) = sqrt(MSE(i));
end

%% Benchmarking AR(1) on yields directly

y_3m = data_full(1,:)';
y_6m = data_full(2,:)';
y_12m = data_full(3,:)';
y_24m = data_full(4,:)';
y_36m = data_full(5,:)';
y_60m = data_full(6,:)';
y_84m = data_full(7,:)';
y_120m = data_full(8,:)';

y_3m_1m = lagmatrix(y_3m, 1);
y_6m_1m = lagmatrix(y_6m, 1);
y_12m_1m = lagmatrix(y_12m, 1);
y_24m_1m = lagmatrix(y_24m, 1);
y_36m_1m = lagmatrix(y_36m, 1);
y_60m_1m = lagmatrix(y_60m, 1);
y_84m_1m = lagmatrix(y_84m, 1);
y_120m_1m = lagmatrix(y_120m, 1);

y_3m_6m = lagmatrix(y_3m, 6);
y_6m_6m = lagmatrix(y_6m, 6);
y_12m_6m = lagmatrix(y_12m, 6);
y_24m_6m = lagmatrix(y_24m, 6);
y_36m_6m = lagmatrix(y_36m, 6);
y_60m_6m = lagmatrix(y_60m, 6);
y_84m_6m = lagmatrix(y_84m, 6);
y_120m_6m = lagmatrix(y_120m, 6);

y_3m_12m = lagmatrix(y_3m, 12);
y_6m_12m = lagmatrix(y_6m, 12);
y_12m_12m = lagmatrix(y_12m, 12);
y_24m_12m = lagmatrix(y_24m, 12);
y_36m_12m = lagmatrix(y_36m, 12);
y_60m_12m = lagmatrix(y_60m, 12);
y_84m_12m = lagmatrix(y_84m, 12);
y_120m_12m = lagmatrix(y_120m, 12);


%% 1-month ahead forecasts: Yields
%% Holdout-period: 2001:3 to 2019:12
clear y_AR1 T1
T1 = 122 %T1 is the first out-of-sample prediction (01.03.2001)
for i= T1:(size(y_3m,1)-1)
%The first observation in b1_changes is the changes during February,
%reported as 01.02.1991. We need data on b1_changes until 01.02.2001.

y_3m_insample = y_3m(1:i);
y_3m_1m_insample = y_3m_1m(1:i);
mdl3m = fitlm(y_3m_1m_insample, y_3m_insample);
y_AR1(1,i-121) = mdl3m.Coefficients.Estimate(1)...
                          + mdl3m.Coefficients.Estimate(2).*y_3m_1m(i+1);

y_6m_insample = y_6m(1:i);
y_6m_1m_insample = y_6m_1m(1:i);
mdl6m = fitlm(y_6m_1m_insample, y_6m_insample);
y_AR1(2,i-121) = mdl6m.Coefficients.Estimate(1)...
                          + mdl6m.Coefficients.Estimate(2).*y_6m_1m(i+1);

y_12m_insample = y_12m(1:i);
y_12m_1m_insample = y_12m_1m(1:i);
mdl12m = fitlm(y_12m_1m_insample, y_12m_insample);
y_AR1(3,i-121) = mdl12m.Coefficients.Estimate(1)...
                          + mdl12m.Coefficients.Estimate(2).*y_12m_1m(i+1);

y_24m_insample = y_24m(1:i);
y_24m_1m_insample = y_24m_1m(1:i);
mdl24m = fitlm(y_24m_1m_insample, y_24m_insample);
y_AR1(4,i-121) = mdl24m.Coefficients.Estimate(1)...
                          + mdl24m.Coefficients.Estimate(2).*y_24m_1m(i+1);

y_36m_insample = y_36m(1:i);
y_36m_1m_insample = y_36m_1m(1:i);
mdl36m = fitlm(y_36m_1m_insample, y_36m_insample);
y_AR1(5,i-121) = mdl36m.Coefficients.Estimate(1)...
                          + mdl36m.Coefficients.Estimate(2).*y_36m_1m(i+1);

y_60m_insample = y_60m(1:i);
y_60m_1m_insample = y_60m_1m(1:i);
mdl60m = fitlm(y_60m_1m_insample, y_60m_insample);
y_AR1(6,i-121) = mdl60m.Coefficients.Estimate(1)...
                          + mdl60m.Coefficients.Estimate(2).*y_60m_1m(i+1);

y_84m_insample = y_84m(1:i);
y_84m_1m_insample = y_84m_1m(1:i);
mdl84m = fitlm(y_84m_1m_insample, y_84m_insample);
y_AR1(7,i-121) = mdl84m.Coefficients.Estimate(1)...
                          + mdl84m.Coefficients.Estimate(2).*y_84m_1m(i+1);

y_120m_insample = y_120m(1:i);
y_120m_1m_insample = y_120m_1m(1:i);
mdl120m = fitlm(y_120m_1m_insample, y_120m_insample);
y_AR1(8,i-121) = mdl120m.Coefficients.Estimate(1)...
                          + mdl120m.Coefficients.Estimate(2).*y_120m_1m(i+1);
end

%% 6-months ahead forecasts: Yields
%% Holdout-period: 2001:3 to 2019:12
clear y_AR1_6m T1
T1 = 127 %T1 is the first out-of-sample prediction (01.03.2001)
for i= T1:(size(y_3m,1)-1)
%The first observation in b1_changes is the changes during February,
%reported as 01.02.1991. We need data on b1_changes until 01.02.2001.

y_3m_insample = y_3m(1:i);
y_3m_6m_insample = y_3m_6m(1:i);
mdl3m = fitlm(y_3m_6m_insample, y_3m_insample);
y_AR1_6m(1,i-126) = mdl3m.Coefficients.Estimate(1)...
                          + mdl3m.Coefficients.Estimate(2).*y_3m_6m(i+1);

y_6m_insample = y_6m(1:i);
y_6m_6m_insample = y_6m_6m(1:i);
mdl6m = fitlm(y_6m_6m_insample, y_6m_insample);
y_AR1_6m(2,i-126) = mdl6m.Coefficients.Estimate(1)...
                          + mdl6m.Coefficients.Estimate(2).*y_6m_6m(i+1);

y_12m_insample = y_12m(1:i);
y_12m_6m_insample = y_12m_6m(1:i);
mdl12m = fitlm(y_12m_6m_insample, y_12m_insample);
y_AR1_6m(3,i-126) = mdl12m.Coefficients.Estimate(1)...
                          + mdl12m.Coefficients.Estimate(2).*y_12m_6m(i+1);

y_24m_insample = y_24m(1:i);
y_24m_6m_insample = y_24m_6m(1:i);
mdl24m = fitlm(y_24m_6m_insample, y_24m_insample);
y_AR1_6m(4,i-126) = mdl24m.Coefficients.Estimate(1)...
                          + mdl24m.Coefficients.Estimate(2).*y_24m_6m(i+1);

y_36m_insample = y_36m(1:i);
y_36m_6m_insample = y_36m_6m(1:i);
mdl36m = fitlm(y_36m_6m_insample, y_36m_insample);
y_AR1_6m(5,i-126) = mdl36m.Coefficients.Estimate(1)...
                          + mdl36m.Coefficients.Estimate(2).*y_36m_6m(i+1);

y_60m_insample = y_60m(1:i);
y_60m_6m_insample = y_60m_6m(1:i);
mdl60m = fitlm(y_60m_6m_insample, y_60m_insample);
y_AR1_6m(6,i-126) = mdl60m.Coefficients.Estimate(1)...
                          + mdl60m.Coefficients.Estimate(2).*y_60m_6m(i+1);

y_84m_insample = y_84m(1:i);
y_84m_6m_insample = y_84m_6m(1:i);
mdl84m = fitlm(y_84m_6m_insample, y_84m_insample);
y_AR1_6m(7,i-126) = mdl84m.Coefficients.Estimate(1)...
                          + mdl84m.Coefficients.Estimate(2).*y_84m_6m(i+1);

y_120m_insample = y_120m(1:i);
y_120m_6m_insample = y_120m_6m(1:i);
mdl120m = fitlm(y_120m_6m_insample, y_120m_insample);
y_AR1_6m(8,i-126) = mdl120m.Coefficients.Estimate(1)...
                          + mdl120m.Coefficients.Estimate(2).*y_120m_6m(i+1);
end

%% 12-months ahead forecasts: Yields
%% Holdout-period: 2001:3 to 2019:12
clear y_AR1_12m T1
T1 = 122 %T1 is the first out-of-sample prediction (01.03.2001)
for i=133:(size(y_3m,1)-1)
%The first observation in b1_changes is the changes during February,
%reported as 01.02.1991. We need data on b1_changes until 01.02.2001.

y_3m_insample = y_3m(1:i);
y_3m_12m_insample = y_3m_12m(1:i);
mdl3m = fitlm(y_3m_12m_insample, y_3m_insample);
y_AR1_12m(1,i-132) = mdl3m.Coefficients.Estimate(1)...
                          + mdl3m.Coefficients.Estimate(2).*y_3m_12m(i+1);

y_6m_insample = y_6m(1:i);
y_6m_12m_insample = y_6m_12m(1:i);
mdl6m = fitlm(y_6m_12m_insample, y_6m_insample);
y_AR1_12m(2,i-132) = mdl6m.Coefficients.Estimate(1)...
                          + mdl6m.Coefficients.Estimate(2).*y_6m_12m(i+1);

y_12m_insample = y_12m(1:i);
y_12m_12m_insample = y_12m_12m(1:i);
mdl12m = fitlm(y_12m_12m_insample, y_12m_insample);
y_AR1_12m(3,i-132) = mdl12m.Coefficients.Estimate(1)...
                          + mdl12m.Coefficients.Estimate(2).*y_12m_12m(i+1);

y_24m_insample = y_24m(1:i);
y_24m_12m_insample = y_24m_12m(1:i);
mdl24m = fitlm(y_24m_12m_insample, y_24m_insample);
y_AR1_12m(4,i-132) = mdl24m.Coefficients.Estimate(1)...
                          + mdl24m.Coefficients.Estimate(2).*y_24m_12m(i+1);

y_36m_insample = y_36m(1:i);
y_36m_12m_insample = y_36m_12m(1:i)
mdl36m = fitlm(y_36m_12m_insample, y_36m_insample)
y_AR1_12m(5,i-132) = mdl36m.Coefficients.Estimate(1)...
                          + mdl36m.Coefficients.Estimate(2).*y_36m_12m(i+1);

y_60m_insample = y_60m(1:i);
y_60m_12m_insample = y_60m_12m(1:i);
mdl60m = fitlm(y_60m_12m_insample, y_60m_insample);
y_AR1_12m(6,i-132) = mdl60m.Coefficients.Estimate(1)...
                          + mdl60m.Coefficients.Estimate(2).*y_60m_12m(i+1);

y_84m_insample = y_84m(1:i);
y_84m_12m_insample = y_84m_12m(1:i);
mdl84m = fitlm(y_84m_12m_insample, y_84m_insample);
y_AR1_12m(7,i-132) = mdl84m.Coefficients.Estimate(1)...
                          + mdl84m.Coefficients.Estimate(2).*y_84m_12m(i+1);

y_120m_insample = y_120m(1:i);
y_120m_12m_insample = y_120m_12m(1:i);
mdl120m = fitlm(y_120m_12m_insample, y_120m_insample);
y_AR1_12m(8,i-132) = mdl120m.Coefficients.Estimate(1)...
                          + mdl120m.Coefficients.Estimate(2).*y_120m_12m(i+1);
end

%% 1-month ahead forecasts: Constructing forecasted yield curves
clear Yields_act MSE RMSE
%We then need the actual yields for these samples: 
Yields_act = data_full(:,123:size(data_full,2));

%We can now caluclate performance statistics
for i=1:size(y_AR1,1)
    MSE(i) = (sum((Yields_act(i,:) - y_AR1(i,:))...
        .^2))./size(y_AR1,2);
    RMSE(i,:) = sqrt(MSE(i));
end

%% 6-month ahead forecasts: Constructing forecasted yield curves
clear Yields_act MSE RMSE
%We then need the actual yields for the sample: 
Yields_act = data_full(:,128:size(data_full,2));

%We can now caluclate performance statistics
for i=1:size(y_AR1_6m,1)
    MSE(i) = (sum((Yields_act(i,:) - y_AR1_6m(i,:))...
        .^2))./size(y_AR1_6m,2);
    RMSE(i,:) = sqrt(MSE(i));
end

%% 12-month ahead forecasts: Constructing forecasted yield curves
clear Yields_act MSE RMSE
%We then need the actual yields for the sample: 
Yields_act = data_full(:,134:size(data_full,2));

%We can now caluclate performance statistics
for i=1:size(y_AR1_12m,1)
    MSE(i) = (sum((Yields_act(i,:) - y_AR1_12m(i,:))...
        .^2))./size(y_AR1_12m,2);
    RMSE(i,:) = sqrt(MSE(i));
end

%% Benchmarking Random Walk
%% 1 month ahead

%We are going to forecast yields from 31.03.2001 to 31.12.2019 by yields
%from 31.02.2001 to 31.11.2019
clear Yields_act Yields_RW RMSE MSE
Yields_RW = data_full(:,122:size(data_full,2)-1);
Yields_act = data_full(:,123:size(data_full,2));

%We can now caluclate performance statistics
for i=1:size(Yields_RW,1)
    MSE(i) = (sum((Yields_act(i,:) - Yields_RW(i,:))...
        .^2))./size(Yields_RW,2);
    RMSE(i,:) = sqrt(MSE(i));
end

%% 6 months ahead

%We are going to forecast yields from 31.03.2001 to 31.12.2019 by yields
%from 31.02.2001 to 31.11.2019
clear Yields_act Yields_RW RMSE MSE
Yields_RW = data_full(:,122:size(data_full,2)-6);
Yields_act = data_full(:,128:size(data_full,2));

%We can now caluclate performance statistics
for i=1:size(Yields_RW,1)
    MSE(i) = (sum((Yields_act(i,:) - Yields_RW(i,:))...
        .^2))./size(Yields_RW,2);
    RMSE(i,:) = sqrt(MSE(i));
end

%% 12 months ahead

%We are going to forecast yields from 31.03.2001 to 31.12.2019 by yields
%from 31.02.2001 to 31.11.2019
clear Yields_act Yields_RW RMSE MSE
Yields_RW = data_full(:,122:size(data_full,2)-12);
Yields_act = data_full(:,134:size(data_full,2));

%We can now caluclate performance statistics
for i=1:size(Yields_RW,1)
    MSE(i) = (sum((Yields_act(i,:) - Yields_RW(i,:))...
        .^2))./size(Yields_RW,2);
    RMSE(i,:) = sqrt(MSE(i));
end