clear all
clc
filename = 'Data_NS_5YRF.xlsx';
Data_table = readtable(filename);
Data=table2array(Data_table);
%% Compute delta
rng(1)
Real_Forward_5Y=Data(:,1);
Policy=Data(:,2);
FOMC_dates=Data(:,3);
FOMC=find(FOMC_dates);
No_FOMC=find(FOMC_dates==0);

T_P=length(Real_Forward_5Y(FOMC));
T_C=length(Real_Forward_5Y(No_FOMC));
beta=0.47;
D_bar_C=mean(Policy(No_FOMC));
Sigma2_hat_D_C=var(Policy(No_FOMC));
rho_hat_D_C=0;
Sigma2_hat_D_P=var(Policy(FOMC));
Sigma2_hat_Y_C=var(Real_Forward_5Y(No_FOMC));
Sigma2_hat_Y_P=var(Real_Forward_5Y(FOMC));
Sigma2_hat_Y_C_prop=Sigma2_hat_Y_P*Sigma2_hat_D_C/Sigma2_hat_D_P;
rho_hat_Y_C=0;
rho_hat_y_P=0;
sims=5000;
rho=[0.1, 0.15, 0.2, 0.25, 0.30, 0.35, 0.4, 0.45, 0.5, 0.55, 0.6, 0.65, 0.7, 0.75, 0.80, 0.85, 0.90, 0.95, 1];
delta=[1.5, 4, 7.3, 12.1, 17.8, 24.8, 32.4, 41.2, 50.8, 61.5,...
    73.5, 86.5, 100.5, 115.3, 132.1,  148.2, 166.5, 186.7, 207.5]*Sigma2_hat_D_C;
%rho=0.40;
%delta=32.43*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
delta_star=diag(MSE_beta_ES_hat_delta);