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);
Policy=Policy(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;

Y=Real_Forward_2Y(FOMC);
D=Policy(FOMC);
Sigma2_hat_D_C=var(Policy(No_FOMC));
Sigma2_hat_D_P=var(Policy(FOMC));
ratio_sigmaD=Sigma2_hat_D_P/Sigma2_hat_D_C;
[estimator,SE,SE_robust,sigma2,resid] = OLS(Y,D,1);
beta_ES=estimator(2,1);
CI_beta_ES=[estimator(2,1)-1.96*SE_robust(2,1)/sqrt(ratio_sigmaD), estimator(2,1)+ 1.96*SE_robust(2,1)/sqrt(ratio_sigmaD)]; 
CI_beta_ES_Sigma2_D_P=[estimator(2,1)-1.96*SE_robust(2,1)/sqrt(Sigma2_hat_D_P), estimator(2,1)+ 1.96*SE_robust(2,1)/sqrt(Sigma2_hat_D_P)]; 


aVar=(SE(2,1)*sqrt(T_P))^2;
sqrt_J=sqrt(aVar);
Sigma2_D=var(D);
c=sqrt(T_P)/sqrt(Sigma2_D);
M=[0.0001 0.1 0.5 1 10 100 300 1000]';
%ratio=[0:0.00001:1]';
%ratio=[0:1:2]';
ratio=[0 0.00001 0.0001 0.001 0.01 0.1 1]';
ratio_actual=sqrt(Sigma2_hat_Y_C*Sigma2_hat_D_C)/sqrt_J;  %ratio_actual=0.00013
for i_M=1:length(M)
    for i_ratio=1:length(ratio)
        b_function(i_M, i_ratio) = normcdf(-1.645 - M(i_M)* ratio(i_ratio)) + normcdf(-1.645 + M(i_M)* ratio(i_ratio));
        aCoverage(i_M, i_ratio)=1- b_function(i_M, i_ratio);
    end
end

M_acutal=c;
aCoverage_actual=1- normcdf(-1.645 - M_acutal * ratio_actual) - normcdf(-1.645 + M_acutal * ratio_actual);
aCoverage
aCoverage_actual

r_actual=M_acutal*ratio_actual;
c_quatile=-[0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.645 1.72 1.8 1.85 1.96 2 2.20 2.3 2.54 2.6 2.8 3 4 4.3 5 6 7 8 9 10 11.3 30 31.3 100 101.3 301.3 1000 1001.3]';
for i_cquantile=1:length(c_quatile)
    for i_M=1:length(M)
    for i_ratio=1:length(ratio)
        b_function_find(i_M, i_ratio ,i_cquantile) = normcdf(c_quatile(i_cquantile) - M(i_M)* ratio(i_ratio)) + normcdf(c_quatile(i_cquantile) + M(i_M)* ratio(i_ratio));
    end
    end
end

cv_ba_val=NaN(8,7);
cv_ba_val(:, 1)=1.645;   %-c_quatile(17,1);
cv_ba_val(:, 2)=1.645;
cv_ba_val(:, 3)=1.645;
cv_ba_val(1:6, 4)=1.645;
cv_ba_val(7, 4)=1.72;
cv_ba_val(8, 4)=2.3;
cv_ba_val(1:5, 5)=1.645;
cv_ba_val(6, 5)=2.3;
cv_ba_val(7, 5)=4.3;
cv_ba_val(8, 5)=11.3;
cv_ba_val(1:4, 6)=1.645;
cv_ba_val(5, 6)=2.3;
cv_ba_val(6, 6)=11.3;
cv_ba_val(7, 6)=31.3;
cv_ba_val(8, 6)=101.3;
cv_ba_val(1:2, 7)=1.645;
cv_ba_val(3, 7)=1.85;
cv_ba_val(4, 7)=2.3;
cv_ba_val(5, 7)=11.3;
cv_ba_val(6, 7)=101.3;
cv_ba_val(7, 7)=301.3;
cv_ba_val(8, 7)=1001.3;



%CI_bias_aware=[beta_ES- val* SE_robust(2,1)/sqrt(ratio_sigmaD), beta_ES+ val*SE_robust(2,1)/sqrt(ratio_sigmaD)]; 
for i_cquantile=1:length(c_quatile)
    for i_M=1:length(M)
    for i_ratio=1:length(ratio)
        equations = @(x) [normcdf(x - M(i_M)* ratio(i_ratio)) + normcdf(x + M(i_M)* ratio(i_ratio))-0.1]
        b_function_find(i_M, i_ratio ,i_cquantile) = normcdf(c_quatile(i_cquantile) - M(i_M)* ratio(i_ratio)) + normcdf(c_quatile(i_cquantile) + M(i_M)* ratio(i_ratio));
        x0        = [-1.645];
        %x0        = [-11.5];
        solution(i_M, i_ratio)  = fsolve(equations, x0);
    end
    end
end

cv_ba=-solution;
cv_ba(8, 5)=cv_ba_val(8, 5);
cv_ba(6:8, 6)=cv_ba_val(6:8, 6);
cv_ba(5:8, 7)=cv_ba_val(5:8, 7);


for i_M=1:length(M)
    for i_ratio=1:length(ratio)
    CI_bias_aware=[beta_ES- cv_ba(i_M, i_ratio)* SE_robust(2,1)/sqrt(ratio_sigmaD), beta_ES+ cv_ba(i_M, i_ratio)*SE_robust(2,1)/sqrt(ratio_sigmaD)]; 
    left_ci(i_M, i_ratio)=CI_bias_aware(1,1);
    right_ci(i_M, i_ratio)=CI_bias_aware(1,2);
    CI_standard=[beta_ES- 1.645* SE_robust(2,1)/sqrt(ratio_sigmaD), beta_ES+ 1.645*SE_robust(2,1)/sqrt(ratio_sigmaD)]; 
    left_ci_std(i_M, i_ratio)=CI_standard(1,1);
    right_ci_std(i_M, i_ratio)=CI_standard(1,2);
    end 
end

lenght_CI_ba=right_ci-left_ci;
lenght_CI_standard=right_ci_std-left_ci_std;

ratio_lenght=lenght_CI_ba./lenght_CI_standard;
Table=ratio_lenght