clearvars
% load 1.3.mat
% load initial_coef
% load SAcoef.mat
% save coef coef
% % %Open the file
% % 
load data

% [~,Y]=hpfilter(Y,1600);
% [~,X]=hpfilter(X,1600);

% Select expected variables from Y matrix
c = Y(3:size(Y,1),2);
pi_c = Y(3:size(Y,1),5);
k = Y(3:size(Y,1),3);
w_k = Y(3:size(Y,1),8);
s = Y(3:size(Y,1),7);
W = Y(3:size(Y,1),1);
pi = Y(3:size(Y,1),4);
pi_m = Y(3:size(Y,1),6);

% Select other variables from X matrix
intr = X(3:size(X,1),1);
bf = X(3:size(X,1),2);
w = X(3:size(X,1),3);
h = X(3:size(X,1),4);
q = X(3:size(X,1),5);
V_n = X(3:size(X,1),6);
p_o = X(3:size(X,1),7);
p_g = X(3:size(X,1),8);
p_m = X(3:size(X,1),9);
x_n = X(3:size(X,1),10);
mu = X(3:size(X,1),11);
e = X(3:size(X,1),12);
V = X(3:size(X,1),13);
z = X(3:size(X,1),14);
M_n = X(3:size(X,1),15);
c_n=X(3:size(X,1),16);
invt=X(3:size(X,1),17);

% calculate the lagged (t-1 subscript) variables
ctm1 = Y(2:size(Y,1)-1,2);
ktm1 = Y(2:size(Y,1)-1,3);
ktm2=Y(1:size(Y,1)-2,3);
Wtm1 = Y(2:size(Y,1)-1,1);
pitm1 = Y(2:size(Y,1)-1,4);
pi_mtm1 = Y(2:size(Y,1)-1,6);
intrtm1 = X(2:size(X,1)-1,1);
x_ntm1 = X(2:size(X,1)-1,10);

% Create VAR from the Y matrix
results=vare(Y,1);
for i=1:8
    beta(:,i)=results(i).beta;
end
Ynew=[Y ones(size(Y,1),1)];
YPred=Ynew*beta;
YPred=YPred(3:end,:);

ctp1 = YPred(:,2);
pi_ctp1 = YPred(:,5);
ktp1 = YPred(:,3);
w_ktp1 = YPred(:,8);
stp1 = YPred(:,7);
Wtp1 = YPred(:,1);
pitp1 = YPred(:,4);
pi_mtp1 = YPred(:,6);

% % %setting the parameters
load coef

theta_y=coef(1,1);
epsilon=coef(1,2);
chi_p=coef(1,3);
t_pm=coef(1,4);
xsi_pm=coef(1,5);
eta_x=coef(1,6);
psi_z=coef(1,7);
psi_hab=coef(1,8);
sigma_c=coef(1,9);
epsilon_k=coef(1,10);
xsi_w=coef(1,11);
psi_w=coef(1,12);
sigma_h=coef(1,13);
theta_rg=coef(1,14);
theta_pdot=coef(1,15);
chi_k=coef(1,16);
phiz=coef(1,22);
sigma_w=coef(1,23);
chi_bf=0.001;
chi_z= 0.0206;
alpha_q=0.1;%0.0528;
alpha_b=0.3154;
delta=0.013;
betap=0.9925;
Q1=0.5802;
Q2=4.7202;
Q4=0.0206*Q2;
Q5=0.2552;
Q6=0.1032;
% calculate indirectly known variables, mrs and B
mrs = (1/sigma_h)*h + (1/sigma_c)*(c+psi_hab*(sigma_c-1)*ctm1);
B = (1-alpha_b)*V_n + alpha_b*M_n;

% calculate residuals
epsb = (((psi_hab*(1-sigma_c))/(1+(psi_hab*(1-sigma_c))))*ctm1) + ((1/(1+psi_hab)*(1-sigma_c))*ctp1)...
       -c/(sigma_c/(1+psi_hab*(1-sigma_c))) -(intr-pi_ctp1-((1/betap)-1));

epsinv = (intr-pi_ctp1-((1/betap)-1)+epsb)...
    -((epsilon_k/(1-delta+chi_z))+(1+epsilon_k))*chi_k*ktm1...
    + (((1+epsilon_k)/(1-delta+chi_z))+1)*chi_k*k...
    -(chi_k/(1-delta+chi_z))*ktp1 - (chi_z/(1-delta+chi_z))*w_ktp1 + chi_k*epsilon_k*ktm2;

epsrf =stp1 - s + (intr-((1/betap)-1)+epsb)+chi_bf*bf;

epsw = W - (xsi_w/(1+betap*xsi_w))*Wtm1 - (betap/(1+betap*xsi_w))*Wtp1 + (psi_w ...
    *(1-betap*(1-psi_w))/((1+(sigma_w/sigma_h))*(1-psi_w)*(1+betap*xsi_w)))*(w-mrs);

epsa = q - (1-alpha_q)*B - alpha_q*e;

epsmu = pi - (betap/(1+(betap*epsilon)))*pitp1 - (epsilon/(1+(betap*epsilon)))*pitm1 ...
    - ((1-chi_p)*(1-betap*chi_p)/((1+betap*epsilon)*chi_p))*mu;

epsi = intr - ((1/betap)-1) - theta_rg*(intrtm1-((1/betap)-1)) - (1-theta_rg)...
   *(theta_pdot*pi_c+theta_y*V);

epspo = p_o + s;

epspg = p_g + s;

epspmf = (pi_m - (t_pm/(1+betap*t_pm))*pi_mtm1 - (betap/(1+betap*t_pm))*pi_mtp1) / ...
    ((1-xsi_pm)*(1-(betap*xsi_pm))/((1+betap*t_pm)*xsi_pm)) +s +p_m;

epsyf = x_n - (psi_z*x_ntm1) + (1-psi_z)*eta_x*s;

epsg = (q - (Q1*c_n) - (Q2*k)+((1-delta)*Q2*ktm1) - (Q4*z) - (Q5*x_n));

Residuals = [epsa epsb epsg epsi epsinv epsmu epspg epspmf epspo epsrf epsw epsyf];

rho = zeros(1,size(Residuals,2));
innov = zeros(size(Y,1)-3,size(Residuals,2));

for i=1:size(Residuals,2)
    y=Residuals(2:end,i);
    x=[ones(size(Residuals,1)-1,1) Residuals(1:end-1,i)];
    results=ols(y,x);
    rho(i)=results.beta(2);
    innov(:,i)=results.resid;
end
    
save rhocoef rho
save innovations innov
stds=std(innov);

save stds stds
% %run dynare to get solutions of the model
dynare Model noclearall
%   Mill;
load EDmatrix
psi = D;
omega = E;
m=27; %Number of State Variables
n=33; %Number of Other Endogenous Variables (controls)

%model simulation after generating innovations from the shocks and 
% using dynare output of the model
T = size(innov,1);
Xsim=zeros(m+n,T,1000);
Xsim_tmp=zeros(m+n,T);
Xsim_tmp_init(:,1) = [B(1); 0; c_n(1); 0; 0; e(1); h(1); 0; 0; 0; 0; M_n(1); mrs(1); mu(1); 0; 0; 0; p_g(1); p_o(1); 0; 0; q(1); 0; 0; V(1); V_n(1); 0; 0; 0; 0; z(1); bf(1); epsa(1); epsb(1); epsg(1); epsi(1); epsinv(1); epsmu(1); epspg(1); epspmf(1); epspo(1); epsrf(1); epsw(1); epsyf(1); intr(1); p_m(1); 0; 0; w(1); x_n(1); c(1); pi_m(1); pi(1); 0; 0; W(1); k(1); pi_c(1); s(1); w_k(1)];           
Xsim_tmp_init(:,2) = [B(2); 0; c_n(2); 0; 0; e(2); h(2); 0; 0; 0; 0; M_n(2); mrs(2); mu(2); 0; 0; 0; p_g(2); p_o(2); 0; 0; q(2); 0; 0; V(2); V_n(2); 0; 0; 0; 0; z(2); bf(2); epsa(2); epsb(2); epsg(2); epsi(2); epsinv(2); epsmu(2); epspg(2); epspmf(2); epspo(2); epsrf(2); epsw(2); epsyf(2); intr(2); p_m(2); 0; 0; w(2); x_n(2); c(2); pi_m(2); pi(2); 0; 0; W(2); k(2); pi_c(2); s(2); w_k(2)];           
Xsim_tmp(:,1:2) = Xsim_tmp_init;
Xsim_tmp(:,1:2) = 0;

rints=randi([1,T],T,1000);
for i=1:1000
    boot_es=innov(rints(:,i),:);
    for j=3:T
        Xsim_tmp(:,j)=psi*Xsim_tmp(:,j-1)+omega*boot_es(j,:)'+D_k.*Xsim_tmp(57,j-2);
    end
    Xsim(:,:,i)=Xsim_tmp;
end

% True VAR for only mu, pi, intr
true_Y = [V pi intr];
vectruevar=[];
truevar=vare(true_Y,1);
for j=1:3
    vectruevar=[vectruevar;truevar(j).beta(1:3)];
end
vectruevar=vectruevar';
% Create 1000 series of mu, pi, and intr with bootstrapped error associated
% with each
boot_V = zeros(1000, size(Xsim,2));
boot_pi = zeros(1000, size(Xsim,2));
boot_intr = zeros(1000, size(Xsim,2));

for i = 1:1000
    boot_V(i,:) = Xsim(25,:,i);
    boot_pi(i,:) = Xsim(53,:,i);
    boot_intr(i,:) = Xsim(45,:,i);
end
%  % %calculate the variances
varV = var(truevar(1,1).resid);
varpi = var(truevar(1,2).resid);
varintr = var(truevar(1,3).resid);
true_variance=[varV varpi varintr];
vectruevar = [vectruevar true_variance];
% vectruevar = [true_variance];

% Create VARs from bootstrapped data
vecbootVARs=[];
for i=1:1000;
    boot_Y = [ boot_V(i,:)' boot_pi(i,:)' boot_intr(i,:)'];
    boots=vare(boot_Y,1);
    tmp=[];
    for j=1:3
        tmp=[tmp;boots(j).beta(1:3)];
    end
%     vecbootVARs(i,:)=tmp';
    varbootV=(var(boots(1,1).resid));
    varbootpi=(var(boots(1,2).resid));
    varbootintr=(var(boots(1,3).resid));
    vecbootVARs(i,:)=[tmp;varbootV;varbootpi;varbootintr]';
%     vecbootVARs(i,:)=[varbootV;varbootpi;varbootintr]';
end

%calculating the Wald and t-stat
betabar=mean(vecbootVARs);
betacov=cov(vecbootVARs);
for i=1:1000
    metric(i,:)=(vecbootVARs(i,:)-betabar)*inv(betacov)*(vecbootVARs(i,:)-betabar)';
end

truemetric = (vectruevar-betabar)*inv(betacov)*(vectruevar-betabar)'
number=zeros(size(truemetric,1),1);
smmetric=sort(metric);
smmetric_950=smmetric(950,1);
smmetric_1000=smmetric(1000,1);
n95=round(i*0.95);
Trans_Mdis_norm=((sqrt(2*truemetric)-sqrt(2*size(vectruevar,2)-1))/(sqrt(2*smmetric(n95,1))-sqrt(2*size(vectruevar,2)-1)))*1.645
% disp(coef');
delete('Mill1.log','Mill1.m','Mill1_fff.m','Mill1_ff.m','Mill1_results.mat');

%%
% figure(4)
% subplot(2,2,1);
% plot(Xsim_tmp(21,:));hold on;plot((q),'r');axis tight
% title('q', 'fontweight', 'bold');
% subplot(2,2,2);
% plot(Xsim_tmp(53,:));hold on;plot((pi),'r');axis tight
% title('infl', 'fontweight', 'bold');
% subplot(2,2,3);
% plot(Xsim_tmp(45,:));hold on;plot((intr),'r');axis tight
% title('intr', 'fontweight', 'bold');
% subplot(2,2,4);
% plot(Xsim_tmp(25,:));hold on;plot((V),'r');axis tight
% title('Va', 'fontweight', 'bold');

% for i=1:20 
%    figure(1); subplot(4,5,i);plot(boot_V(i,:),'r');hold on;plot(V);axis tight;
%    figure(2); subplot(4,5,i);plot(boot_pi(i,:),'r');hold on;plot(pi);axis tight;
%    figure(3); subplot(4,5,i);plot(boot_intr(i,:),'r');hold on;plot(intr);axis tight;
% end