clear all
clc
%% 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
%% Compute delta
rng(1)
T_P=74;
T_C=762;
beta=0.99;
D_bar_C=0.006050142;
Sigma2_hat_D_C=0.00003688782;
rho_hat_D_C=0;
rho_hat_Y_C=0;
Sigma2_hat_D_P=0.00119797;
Sigma2_hat_Y_P=0.00651286;
Sigma2_hat_Y_C_prop=Sigma2_hat_Y_P*Sigma2_hat_D_C/Sigma2_hat_D_P;
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=[0.78, 2.5, 4.60, 7.60, 11.40, 15.80, 21.00, 26.80, 33.50, 40.80, 49.3, 56.8, 66, 75.9, 86.5, 97.8, 109.7, 122.3... 
    135.5]*Sigma2_hat_D_C;
%rho=0.49;
%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);
