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);
[estimator,SE,SE_robust,sigma2,resid] = OLS(Y,D,1);
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
