clear all
clc
rng(1)
filename = 'Data_NS18_DYRF.xlsx';
Data_table = readtable(filename);
Data=table2array(Data_table);
%% Input from Nakamura and Steinsson
% 2-Year U.S. Treasury instantaneous real forward rates
% Sigma_2_D_C = 0.00003688782
% Mu_D_C = 0.006050142
% Sigma_2_Y_C = 0.003422089
% Sigma_2_D_T = 0.00119797
% Sigma_2_Y_T=0.00651286  Variance of Y in the policy sample
% rho_1_Y_C =0
% rho_1_D_C =0
%% Procedure
Policy=Data(:,2);
FOMC_dates=Data(:,3);
FOMC=find(FOMC_dates);
No_FOMC=find(FOMC_dates==0);
D_bar_C=mean(Policy(No_FOMC));
Sigma2_hat_D_C=var(Policy(No_FOMC));
Sigma2_hat_D_P=var(Policy(FOMC));
ratio=Sigma2_hat_D_P/Sigma2_hat_D_C;

Real_Forward_2Y=Data(:,4);
ind = ~isnan(Real_Forward_2Y);
Real_Forward_2Y=Real_Forward_2Y(ind);
FOMC_dates=FOMC_dates(ind);
FOMC=find(FOMC_dates);
No_FOMC=find(FOMC_dates==0);

T_P=length(Real_Forward_2Y(FOMC));
T_C=length(Real_Forward_2Y(No_FOMC));
beta=0.99;
rho_hat_D_C=0;
%Sigma2_hat_Y_C=var(Real_Forward_2Y(No_FOMC));
Sigma2_hat_Y_P=var(Real_Forward_2Y(FOMC));
Sigma2_hat_Y_C_prop=Sigma2_hat_Y_P*Sigma2_hat_D_C/Sigma2_hat_D_P;
Sigma2_hat_Y_C_prop_temp=Sigma2_hat_Y_C_prop;

rho_hat_Y_C=0;
sims=5000;
rho=[0:0.25:1];
delta=[0, 1, 8, 16, 32, 40, 48]*Sigma2_hat_D_C;
for i=1:sims
    D_tilde(1)=0;
    for t=2:T_P+1
        D_tilde(t,1)=D_bar_C + rho_hat_D_C*D_tilde(t-1,1) + sqrt((1-rho_hat_D_C^2)*Sigma2_hat_D_C)*randn;
    end
    D_tilde(1)=[];
    u_tilde(1)=0;
    for t=2:T_P+1
        u_tilde(t,1)= rho_hat_Y_C*u_tilde(t-1,1) + sqrt((1-rho_hat_Y_C^2)*Sigma2_hat_Y_C_prop)*randn;
    end
    u_tilde(1)=[];
    D_tilde_store(:,i)=D_tilde;
    Y_tilde = beta*D_tilde + u_tilde;
    beta_oracle_hat(i,1)=OLS(Y_tilde, D_tilde);
end
 Bias_beta_oracle=mean(beta_oracle_hat)-beta;
 MAE_beta_oracle=mean(abs(beta_oracle_hat-beta));
 MSE_beta_oracle=mean((beta_oracle_hat-beta).^2);

u_tilde_oracle=u_tilde;
Y_tilde_oracle=Y_tilde;
clear u_tilde Y_tilde

for i_rho=1:length(rho)
    for i_delta=1:length(delta)
for ii=1:sims
    for t=1:T_P
        u_tilde(t,1)= rho(i_rho)*sqrt(Sigma2_hat_Y_C_prop)* (sqrt(Sigma2_hat_D_C))^-1*(D_tilde_store(t,ii)-D_bar_C) + sqrt((1-rho(i_rho)^2)*Sigma2_hat_Y_C_prop)*randn;
    end

    Y_tilde=beta*D_tilde_store(:,ii) + u_tilde;
    beta_ES_hat(ii,i_rho)=OLS(Y_tilde, D_tilde_store(:,ii));
  
    D_tilde_delta = D_bar_C + (sqrt(Sigma2_hat_D_C))^-1*(D_tilde_store(:,ii) - D_bar_C)* sqrt(Sigma2_hat_D_C + delta(i_delta));
    Y_tilde=beta*D_tilde_delta + u_tilde;
    beta_ES_hat_delta(ii,i_rho,i_delta) = OLS(Y_tilde, D_tilde_delta);
    clear u_tilde u0_tilde
end
%% Bias, MAE, MSE
Bias_beta_ES_hat(i_rho,1)=mean(beta_ES_hat(:,i_rho))-beta;
MAE_beta_ES_hat(i_rho,1)=mean(abs(beta_ES_hat(:,i_rho)-beta));
MSE_beta_ES_hat(i_rho,1)=mean((beta_ES_hat(:,i_rho)-beta).^2);

Bias_beta_ES_hat_delta(i_rho,i_delta)=mean(beta_ES_hat_delta(:,i_rho,i_delta))-beta;
MAE_beta_ES_hat_delta(i_rho,i_delta)=mean(abs(beta_ES_hat_delta(:,i_rho,i_delta)-beta));
MSE_beta_ES_hat_delta(i_rho,i_delta)=mean((beta_ES_hat_delta(:,i_rho,i_delta)-beta).^2);
Variance_beta_ES_hat_delta(i_rho,i_delta)=var(beta_ES_hat_delta(:,i_rho,i_delta));
    end
end

%% Variance of Treasury Yields
%Sigma2_hat_Y_C=var(Real_Forward_5Y(No_FOMC))
Sigma2_hat_Y_P=var(Real_Forward_2Y(FOMC));
Real_Forward_2Y_P_0014=Real_Forward_2Y(FOMC);
Sigma2_hat_Y_P=var(Real_Forward_2Y_P_0014)

Real_Forward_2Y_P_0007=Real_Forward_2Y(1:416);
Real_Forward_2Y_P_0007=Real_Forward_2Y_P_0007(FOMC(1:32));
Sigma2_hat_Y_P_0007=var(Real_Forward_2Y_P_0007)

