% Check if our term premium keep the feature of Cochrane and Piazzesi
% (2005 AER). They use annual CRSP 1964:1-2001:12
% Our sample: monthly, 1971:10-2022:04,  607 periods

clear all
close all 
clc
addpath(genpath('./functions'));

mats = (1:1:120)';
mcmc_ind = 50000;
dt = 1;
rlb = 0;
mat_ind = (12:12:60)';
% mat_ind = (12:12:120)';


% load parameters and factors

load('./results/temp100000.mat','USdata','X','kinfQ',...
    'lamQ','Omega','sige2');
Date = USdata{:,'Date'};

cp_est2 = nan(50000,length(mat_ind)+1);
cp_est3 = nan(50000,length(mat_ind)+1);
cp_est4 = nan(50000,length(mat_ind)+1);
cp_est5 = nan(50000,length(mat_ind)+1);
cp_r2 = nan(50000,length(mat_ind)-1);


rng(616)
for i = mcmc_ind+1:100000
    if mod(i,1000)==0
        disp(i-50000);
    end
    omega = [Omega(i,1) Omega(i,2) Omega(i,3);...
        Omega(i,2) Omega(i,4) Omega(i,5);...
        Omega(i,3) Omega(i,5) Omega(i,6)];
    sim_data = y_sim(X(3*(i-1)+1:3*i,:), diag(lamQ(i,:)), kinfQ(i),...
        omega, sige2(i), rlb, mats, dt); 
    sim_yields = sim_data.yields;   % jjxT
    sim_fwdr = sim_data.fwdr;       % (jj+1)xT
    % excess return for n = 2,3,...,K; t= 13,14,...,T
    rx = excess_r(sim_yields, mat_ind); 
    % averaged excess return
    % regressors: 1 y_{1,t},f_{1,t}, f_{2,t}, ... , f_{k-1,t}
    Y1 = sim_yields(12,1:end-12)';
    myf = sim_fwdr([13 25 37 49],1:end-12)';
    myx = [Y1 myf];
%     % regress average excess return on myx
    cp_reg2 = fitlm(myx, rx(1,:)');
    cp_reg3 = fitlm(myx, rx(2,:)');
    cp_reg4 = fitlm(myx, rx(3,:)');
    cp_reg5 = fitlm(myx, rx(4,:)');
    cp_est2(i-mcmc_ind,:) = cp_reg2.Coefficients{:,1};
    cp_est3(i-mcmc_ind,:) = cp_reg3.Coefficients{:,1};
    cp_est4(i-mcmc_ind,:) = cp_reg4.Coefficients{:,1};
    cp_est5(i-mcmc_ind,:) = cp_reg5.Coefficients{:,1};
    cp_r2(i-mcmc_ind,1) = cp_reg2.Rsquared.Adjusted;
    cp_r2(i-mcmc_ind,2) = cp_reg3.Rsquared.Adjusted;
    cp_r2(i-mcmc_ind,3) = cp_reg4.Rsquared.Adjusted;
    cp_r2(i-mcmc_ind,4) = cp_reg5.Rsquared.Adjusted;
    
    
end

save('./results/CP_results_LW','cp_est2','cp_est3','cp_est4',...
    'cp_est5','cp_r2');


cp_est_mean = nan(4,length(mat_ind)+1);
cp_est_mean(1,:) = mean(cp_est2,1);
cp_est_mean(2,:) = mean(cp_est3,1);
cp_est_mean(3,:) = mean(cp_est4,1);
cp_est_mean(4,:) = mean(cp_est5,1);

cp_r2_mean = mean(cp_r2,1);




%% 

% Plot Coefficients of excess returns on forward rates
plot(cp_est_mean(1,2:end),'-o','LineWidth',1.5)
xticks([1 2 3 4 5])
hold on
plot(cp_est_mean(2,2:end),'-v','LineWidth',1.5)
hold on
plot(cp_est_mean(3,2:end),'-s','LineWidth',1.5)
hold on 
plot(cp_est_mean(4,2:end),'-d','LineWidth',1.5)
legend({'24','36','48','60'},'Location','southeast','FontSize',30)

    







 