close all
clear all
clc

Boot=3000;
NumObs=0;load('CLEAN')
load('CDM_RECESSIONS')

iF_YEAR = 1991;   %  The first year of differences
iL_YEAR = 2008;   %  The last year of differences

iNUM_YEARS = iL_YEAR - iF_YEAR + 1;
iOBS = 3;  %The columns including the data
iAGE =  5;  %The columns indicating whether the data exist
mDATA_BHPS_IN1 = load('-ascii','../data/residuals_Y_L_U.out');
mDATA_BHPS_IN2 = load('-ascii','../data/residuals_Y_L_NU.out');
mDATA_BHPS_IN = [mDATA_BHPS_IN1; mDATA_BHPS_IN2];

log_DATA_TO_KEEP = mDATA_BHPS_IN(:,2) >= iF_YEAR & mDATA_BHPS_IN(:,2) <= iL_YEAR ; 
mDATA_BHPS = mDATA_BHPS_IN(log_DATA_TO_KEEP,:);
num_entries = size(mDATA_BHPS,1)/iNUM_YEARS;

log_DATA_TO_KEEP1 = mDATA_BHPS_IN1(:,2) >= iF_YEAR & mDATA_BHPS_IN1(:,2) <= iL_YEAR ; 
mDATA_BHPS1 = mDATA_BHPS_IN1(log_DATA_TO_KEEP1,:);
num_entries1 = size(mDATA_BHPS1,1)/iNUM_YEARS;

log_DATA_TO_KEEP2 = mDATA_BHPS_IN2(:,2) >= iF_YEAR & mDATA_BHPS_IN2(:,2) <= iL_YEAR ; 
mDATA_BHPS2 = mDATA_BHPS_IN2(log_DATA_TO_KEEP2,:);
num_entries2 = size(mDATA_BHPS2,1)/iNUM_YEARS;

initial_points=[];

 for k=1:3
     for i=1:3
         for j=1:3
             X0(1)=log(0.01+(i-1)*0.03); %sig2etaR
             X0(2)=log(0.01+(j-1)*0.03);%sig2etaE
             X0(3)=-0.05+(i-1)*0.05; %sk_etaR
             X0(4)=-0.05+(k-1)*0.05; %sk_etaE
             X0(5)=log(1+0.60+(j-1)*0.30)-log(1-0.60+(j-1)*0.30); %log((1+rho)/(1-rho))
             X0(6)=log(0.01+(i-1)*0.03); %sig2eps
             X0(7)=-0.05+(i-1)*0.05; %sk_eps
             X0(8)=log(0.05+(j-1)*0.05);%sd_fx
             X0(9)=-0.05+(k-1)*0.05; %sk_fx
             X0(10)=-0.01+(j-1)*0.01; %sk_fx
             X0(11)=-0.01+(i-1)*0.01; %sk_fx
             initial_points=[initial_points;X0];
         end
     end
 end
 
 initial_points2=[];

 for k=1:3
     for i=1:3
         for j=1:3
             X1(1)=log(0.02+(k-1)*0.04); %kurtetaR
             X1(2)=log(0.02+(j-1)*0.04); %kurtetaE
             X1(3)=log(0.02+(i-1)*0.06); %kurteps
             X1(4)=log(0.02+(i-1)*0.04); %kurtchi
             X1(5)=-0.01+(j-1)*0.01; %sk_fx
             initial_points2=[initial_points2;X1];
         end
     end
 end

%%%%%%%%%%%%%%%%%%%%%%%% 
%%%%%%       Original Data         %%%%%% 
%%%%%%%%%%%%%%%%%%%%%%%% 
DD = mDATA_BHPS;
%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%  Empirical Moments   %%%%%%
%%%%%%%%%%%%%%%%%%%%%%%
AGEi = reshape(DD(:,iAGE),iNUM_YEARS,num_entries);
Ui=reshape(DD(:,iOBS),iNUM_YEARS,num_entries);
[m2htk,m3htk,m4htk,Ihtk]=empirical_moments_double_sq(AGEi,Ui,CLEAN,iNUM_YEARS,num_entries);
%%%%%%%%%%%%%%%%%%%%%%
%%%%%%      Contractions      %%%%%%
%%%%%%%%%%%%%%%%%%%%%%

Contraction_index=cdm;
%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%    Estimation    %%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%
tic
[P,f]        = MDE_Timevarying_main_boot6(iNUM_YEARS,m2htk,m3htk,m4htk,Ihtk,initial_points,Contraction_index,NumObs);
toc
Estimated_Values_1ststep=P;
tic
[P,f]        = MDE_Timevarying_main_boot3(iNUM_YEARS,m2htk,m3htk,m4htk,Ihtk,initial_points2,Contraction_index,NumObs,Estimated_Values_1ststep);
toc
Estimated_Values_2ndstep=P;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%BOOTSTRAP%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Bootstrapped_Values1=zeros(Boot,11);
Bootstrapped_Values2=zeros(Boot,5);
Bootstrapped_f1=zeros(Boot,1);
Bootstrapped_f2=zeros(Boot,1);
rng('default')
 
ADATA1=mDATA_BHPS_IN1(:,iAGE);
MDATA1=mDATA_BHPS_IN1(:,iOBS);
DD1=[ADATA1 MDATA1];
ADATA2=mDATA_BHPS_IN2(:,iAGE);
MDATA2=mDATA_BHPS_IN2(:,iOBS);
DD2=[ADATA2 MDATA2];
    
parfor b=1:1:Boot
    b
 %%%%%%%%%%%%%%%%%%%%%%%% 
%%%%%%%%%%%%%%%%%%%%%%%% 
BSample1 = reshape(DD1,[iNUM_YEARS,num_entries1,2]);
BSample2 = reshape(DD2,[iNUM_YEARS,num_entries2,2]);
y1 = datasample(BSample1,size(BSample1,2),2);
y2 = datasample(BSample2,size(BSample2,2),2);
A1 = squeeze(y1(:,:,1));
G1 = squeeze(y1(:,:,2));
A2 = squeeze(y2(:,:,1));
G2 = squeeze(y2(:,:,2));
AGEi=[A1 A2];
Ui=[G1 G2];

[m2htk,m3htk,m4htk,Ihtk]=empirical_moments_double_sq(AGEi,Ui,CLEAN,iNUM_YEARS,num_entries);

[P,f]        = MDE_Timevarying_main_boot6(iNUM_YEARS,m2htk,m3htk,m4htk,Ihtk,initial_points,Contraction_index,NumObs);

V=real(P);

Bootstrapped_Values1(b,:)=V;
Bootstrapped_f1(b)=f;

[P,f]        = MDE_Timevarying_main_boot3(iNUM_YEARS,m2htk,m3htk,m4htk,Ihtk,initial_points2,Contraction_index,NumObs,V);

V=real(P);
    
Bootstrapped_Values2(b,:)=V;
Bootstrapped_f2(b)=f;

end

save('results_Y_L_RECESSIONS','Estimated_Values_1ststep','Estimated_Values_2ndstep','Bootstrapped_Values1','Bootstrapped_Values2','Bootstrapped_f1','Bootstrapped_f2')
