close all
clear all
clc

load('../Section4.1/results_Y_L_RECESSIONS');
mom=[Estimated_Values_1ststep(1) Estimated_Values_1ststep(3)  Estimated_Values_2ndstep(1) ...
     Estimated_Values_1ststep(2) Estimated_Values_1ststep(4)  Estimated_Values_2ndstep(2)];

UB=[0.99 0.99  3   3    3       3 ];
LB=[0.5  0.5  -3  -3  0.000001 0.000001 ];
X0=[0.6  0.6  0.1 0.1  0.1     0.2 ];
p = haltonset(6);
pt = net(p,100).*(UB-LB)+LB;
tpoints = CustomStartPointSet(pt);

ms = MultiStart('FunctionTolerance',1e-4,'XTolerance',1e-4)    
poolobj = gcp('nocreate');
delete(poolobj);
parpool
ms.UseParallel='always'

counter=0;
concC80=[];
concE80=[];

problem1=createOptimProblem('lsqnonlin','objective',@(xx) objective_function_MG3( xx, mom'),'x0',X0,'lb',LB,'ub',UB, 'options',optimset('Display','iter','MaxIter',3000,'TolX',1e-4,'TolFun',1e-4,'UseParallel',true));

[FIT,fval,exitflag,output,solutions]=run(ms,problem1,tpoints);

pp1=FIT(1);
pp2=FIT(2);
mu1=FIT(3);
mu2=pp1*mu1/(pp1-1);
mu3=FIT(4);
mu4=pp2*mu3/(pp2-1);
sigma1=FIT(5);
sigma2=FIT(6);

mu_C       =pp1*mu1+(1-pp1)*mu2;
variance_C =pp1*(mu1^2+sigma1)+(1-pp1)*(mu2^2+sigma2);
Third_mom_C=pp1*(mu1^3+3*mu1*sigma1)+(1-pp1)*(mu2^3+3*mu2*sigma2);
Fouth_mom_C=pp1*(mu1^4+6*(mu1^2)*sigma1+3*sigma1^2)+(1-pp1)*(mu2^4+6*(mu2^2)*sigma2+3*sigma2^2);

mu_E       =pp2*mu3+(1-pp2)*mu4;
variance_E =pp2*(mu3^2+sigma1)+(1-pp2)*(mu4^2+sigma2);
Third_mom_E=pp2*(mu3^3+3*mu3*sigma1)+(1-pp2)*(mu4^3+3*mu4*sigma2);
Fouth_mom_E=pp2*(mu3^4+6*(mu3^2)*sigma1+3*sigma1^2)+(1-pp2)*(mu4^4+6*(mu4^2)*sigma2+3*sigma2^2);

zmin=-5;
zmax=5;
N=10000;
%create the grid
dz=(zmax-zmin)/(N-1); %distance between two consecutive gridpoints
zgrid=(zmin:dz:zmax)';

y1=zeros(N,1);
y2=zeros(N,1);
for i=2:N-1
    y1(i)=pp2*(normcdf((zgrid(i)+dz/2),mu3,sqrt(sigma1))-normcdf((zgrid(i)-dz/2),mu3,sqrt(sigma1)))...
        +(1-pp2)*(normcdf((zgrid(i)+dz/2),mu4,sqrt(sigma2))-normcdf((zgrid(i)-dz/2),mu4,sqrt(sigma2)));
end

for i=2:N-1
    y2(i)=pp1*(normcdf((zgrid(i)+dz/2),mu1,sqrt(sigma1))-normcdf((zgrid(i)-dz/2),mu1,sqrt(sigma1)))...
        +(1-pp1)*(normcdf((zgrid(i)+dz/2),mu2,sqrt(sigma2))-normcdf((zgrid(i)-dz/2),mu2,sqrt(sigma2)));
end

idx1 = find(exp(zgrid) >=1.40);
probA_1=sum(y1(idx1));
idx1 = find(exp(zgrid) >=1.20);
probA_2=sum(y1(idx1));
idx1 = find(exp(zgrid) <=0.50);
probA_3=sum(y1(idx1));
idx1 = find(exp(zgrid) <=0.25);
probA_4=sum(y1(idx1));

expansions=round((100*[probA_1;probA_2;probA_3;probA_4]),2);

idx1 = find(exp(zgrid) >=1.40);
probB_1=sum(y2(idx1));
idx1 = find(exp(zgrid) >=1.20);
probB_2=sum(y2(idx1));
idx1 = find(exp(zgrid) <= 0.50);
probB_3=sum(y2(idx1));
idx1 = find(exp(zgrid) <= 0.25);
probB_4=sum(y2(idx1));

recessions=round((100*[probB_1;probB_2;probB_3;probB_4]),2);


zmin=-5;
zmax=5;
N=10000;
%create the grid
dz=(zmax-zmin)/(N-1); %distance between two consecutive gridpoints
zgrid=(zmin:dz:zmax)';

Prz1=zeros(N,1);
Prz2=zeros(N,1);
for i=2:N-1
    Prz1(i)=normcdf((zgrid(i)+dz/2),0,sqrt((zgrid.^2)'*y1))-normcdf((zgrid(i)-dz/2),0,sqrt((zgrid.^2)'*y1));
end

for i=2:N-1
    Prz2(i)=normcdf((zgrid(i)+dz/2),0,sqrt((zgrid.^2)'*y2))-normcdf((zgrid(i)-dz/2),0,sqrt((zgrid.^2)'*y2));
end

PP1=cumsum(Prz1);
PP2=cumsum(Prz2);

idx1 = find(exp(zgrid) > 1.40);
probA_1=sum(Prz1(idx1));
idx1 = find(exp(zgrid) > 1.20);
probA_2=sum(Prz1(idx1));
idx1 = find(exp(zgrid) < 0.50);
probA_3=sum(Prz1(idx1));
idx1 = find(exp(zgrid) < 0.25);
probA_4=sum(Prz1(idx1));

expansions_Normal=round((100*[probA_1;probA_2;probA_3;probA_4]),2);

idx1 = find(exp(zgrid) > 1.40);
probB_1=sum(Prz2(idx1));
idx1 = find(exp(zgrid) > 1.20);
probB_2=sum(Prz2(idx1));
idx1 = find(exp(zgrid) < 0.50);
probB_3=sum(Prz2(idx1));
idx1 = find(exp(zgrid) < 0.25);
probB_4=sum(Prz2(idx1));

recession_Normal=round((100*[probB_1;probB_2;probB_3;probB_4]),2);

name = {'income gain>40%';...
        'income gain>20%';...
        'income loss>50%';...
        'income loss>75%'};
    
TABLE_3 = table(recessions,expansions,recession_Normal,expansions_Normal,'RowNames',name);
disp('Table 3: Probabilities of big shocks')
disp(TABLE_3)

function [ out ] = objective_function_MG3( theta, mom)

pp1=theta(1);
pp2=theta(2);
mu1=theta(3);
mu2=pp1*mu1/(pp1-1);
mu3=theta(4);
mu4=pp2*mu3/(pp2-1);
sigma1=theta(5);
sigma2=theta(6);

mu_C       =pp1*mu1+(1-pp1)*mu2;
variance_C =pp1*(mu1^2+sigma1)+(1-pp1)*(mu2^2+sigma2);
Third_mom_C=pp1*(mu1^3+3*mu1*sigma1)+(1-pp1)*(mu2^3+3*mu2*sigma2);
Fouth_mom_C=pp1*(mu1^4+6*(mu1^2)*sigma1+3*sigma1^2)+(1-pp1)*(mu2^4+6*(mu2^2)*sigma2+3*sigma2^2);

mu_E       =pp2*mu3+(1-pp2)*mu4;
variance_E =pp2*(mu3^2+sigma1)+(1-pp2)*(mu4^2+sigma2);
Third_mom_E=pp2*(mu3^3+3*mu3*sigma1)+(1-pp2)*(mu4^3+3*mu4*sigma2);
Fouth_mom_E=pp2*(mu3^4+6*(mu3^2)*sigma1+3*sigma1^2)+(1-pp2)*(mu4^4+6*(mu4^2)*sigma2+3*sigma2^2);

theo_mom=[variance_C;Third_mom_C;Fouth_mom_C;variance_E;Third_mom_E;Fouth_mom_E];
out=(theo_mom-mom);

end
