clear all
clc
rng(1)
filename = 'Data_NS18_DYRF.xlsx';
Data_table = readtable(filename);
Data=table2array(Data_table);

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));



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));
Sigma2_hat_Y_P=var(Real_Forward_2Y(FOMC));

ratio_Sigma_D_fullsample=Sigma2_hat_D_P/Sigma2_hat_D_C;


%% Check Var(D) in sub-samples
bandwith=[4/5];
ratio_Sigma2_D=NaN(10,length(bandwith));
Real_Forward_2Y_P=Real_Forward_2Y(FOMC);
Policy=Policy(ind);
% Compute Var(D) for sample post 2004
Sigma2_hat_D_C_post04=var(Policy(No_FOMC));
Sigma2_hat_D_P_post04=var(Policy(FOMC));
ratio_Sigma_D_post04=Sigma2_hat_D_P_post04/Sigma2_hat_D_C_post04;

Policy_P=Policy(FOMC);
beta_hat = OLS(Real_Forward_2Y_P, Policy_P,1);
beta_ES_hat_NS=beta_hat(2);

for i_band=1:length(bandwith)
n_T_P=floor(T_P^bandwith(i_band)); %4/5
n_block_P=floor(T_P/n_T_P);
Sigma2_D_P_block(1,1)=var(Policy_P(1:n_T_P));
for i_block=2:n_block_P
    Sigma2_D_P_block(i_block,1)= var(Policy_P( (i_block-1)*n_T_P+1:i_block*n_T_P ));
end
Sigma2_D_P_block(n_block_P+1,1)= var(Policy_P( (n_block_P)*n_T_P+1:T_P));  

% Eliminate the last 10 FOMC days in order to balance the sample
Policy_P=Policy_P(1:end-10);
Real_Forward_2Y_P=Real_Forward_2Y_P(1:end-10);
T_P=T_P-10;
Policy_C=Policy(No_FOMC);
n_T_P=floor(T_P^bandwith(i_band)); %4/5
n_block_P=floor(T_P/n_T_P);
n_block_C=n_block_P;
n_T_C=floor(T_C/ n_block_C);

Sigma2_D_P_block(1,1)=var(Policy_P(1:n_T_P));
Sigma2_D_C_block(1,1)=var(Policy_C(1:n_T_C));
for i_block=2:n_block_P
    Sigma2_D_P_block(i_block,1)= var(Policy_P( (i_block-1)*n_T_P+1:i_block*n_T_P ));
end
Sigma2_D_P_block(n_block_P+1,1)= var(Policy_P( (n_block_P)*n_T_P+1:T_P));  
for i_block=2:n_block_C
    Sigma2_D_C_block(i_block,1)= var(Policy_C( (i_block-1)*n_T_C+1:i_block*n_T_C ));
end
Sigma2_D_C_block(n_block_C+1,1)= var(Policy_C( (n_block_C)*n_T_C+1:T_C));  

beta_hat = OLS(Real_Forward_2Y_P, Policy_P,1);
beta_ES_hat_sub=beta_hat(2);

ratio_Sigma2_D(1:length(Sigma2_D_P_block),i_band)=Sigma2_D_P_block./Sigma2_D_C_block;
end



