
run_simulation=0;

    %read in estimated parameters
    opt = {'CommentStyle','%', 'CollectOutput',true};
    [fid,msg] = fopen('param_est.txt','rt');
    assert(fid>=3,msg)
    C = textscan(fid, '%f', opt{:});
    fclose(fid);
    params = C{1};

    PARAMHAT       = params;
    PARAMSE        = [0.071 0.126 0.327 0.191 0.078 0.042];
    PARAMHAT_LB    = PARAMHAT - 1.96*PARAMSE;
    PARAMHAT_UB    = PARAMHAT + 1.96*PARAMSE;
    PARAMHAT_UB(5) = 1.0; %q cannot be above 1

    nbeta=2;

    stepr = [0.00 0.25 0.50 0.75 0.99];
    stepw = [0.01 0.50 1.00 1.50 2.00];
    stepe = [0.01 0.50 1.00 1.50 2.00];
    stepv = [0.01 0.50 1.00 1.50 2.00];
    stepp = [1.0 0.75 0.50 0.25 0.0];

if run_simulation==1

    %% RHO
    rhomatmom = zeros(nbeta,11,5);

    %Solve at estimated values
    thetavec=PARAMHAT;
    for stepct=1:length(stepr)
        thetavec(1) = stepr(stepct);
        rhomatmom(:,1,stepct) = SMM_bgms(thetavec);
    end

    ct=1;
    paramct=1;  
    while paramct<=5

        paramct=paramct+1; 

        %Solve at LB
        ct=ct+1;
        thetavec=PARAMHAT;
        thetavec(paramct)          = PARAMHAT_LB(paramct);
        for stepct=1:length(stepr)
            thetavec(1)            = stepr(stepct);
            rhomatmom(:,ct,stepct) = SMM_bgms(thetavec);
        end 

        %Solve at UB
        ct=ct+1;
        thetavec=PARAMHAT;
        thetavec(paramct)          = PARAMHAT_UB(paramct);
        for stepct=1:length(stepr)
            thetavec(1)            = stepr(stepct);
            rhomatmom(:,ct,stepct) = SMM_bgms(thetavec);
        end 
    end

    %% SIGMA_W
    sigwmatmom = zeros(nbeta,11,5);

    %Solve at estimated values
    thetavec                   = PARAMHAT;
    for stepct=1:length(stepw)
        thetavec(2)            = stepw(stepct);
        sigwmatmom(:,1,stepct) = SMM_bgms(thetavec);
    end

    ct=1;
    paramct=0;  
    while paramct<=5

        if paramct==2
            paramct=paramct+1;
            continue
        end

        paramct=paramct+1; 

        %Solve at LB
        ct=ct+1;
        thetavec                    = PARAMHAT;
        thetavec(paramct)           = PARAMHAT_LB(paramct);
        for stepct=1:length(stepw)
            thetavec(2)             = stepw(stepct);
            sigwmatmom(:,ct,stepct) = SMM_bgms(thetavec);
        end

        %Solve at UB
        ct=ct+1;
        thetavec                    = PARAMHAT;
        thetavec(paramct)           = PARAMHAT_UB(paramct);
        for stepct=1:length(stepw)
            thetavec(2)             = stepw(stepct);
            sigwmatmom(:,ct,stepct) = SMM_bgms(thetavec);
        end 
    end


    %% SIGMA_E
    sigematmom = zeros(nbeta,11,5);

    %Solve at estimated values
    thetavec                   = PARAMHAT;
    for stepct=1:length(stepe)
        thetavec(3)            = stepe(stepct);
        sigematmom(:,1,stepct) = SMM_bgms(thetavec);
    end

    ct=1;
    paramct=0;  
    while paramct<=5

        if paramct==3
            paramct=paramct+1;
            continue
        end

        paramct=paramct+1; 

        %Solve at LB
        ct=ct+1;
        thetavec                    = PARAMHAT;
        thetavec(paramct)           = PARAMHAT_LB(paramct);
        for stepct=1:length(stepe)
            thetavec(3)             = stepe(stepct);
            sigematmom(:,ct,stepct) = SMM_bgms(thetavec);
        end

        %Solve at UB
        ct=ct+1;
        thetavec=PARAMHAT;
        thetavec(paramct)           = PARAMHAT_UB(paramct);
        for stepct=1:length(stepe)
            thetavec(3)             = stepe(stepct);
            sigematmom(:,ct,stepct) = SMM_bgms(thetavec);
        end 
    end

    %% SIGMA_V
    sigvmatmom = zeros(nbeta,11,5);

    %Solve at estimated values
    thetavec = PARAMHAT;
    for stepct=1:length(stepv)
        thetavec(4)            = stepv(stepct);
        sigvmatmom(:,1,stepct) = SMM_bgms(thetavec);
    end

    ct=1;
    paramct=0;  
    while paramct<=5

        if paramct==4
            paramct=paramct+1;
            continue
        end

        paramct=paramct+1; 

        %Solve at LB
        ct=ct+1;
        thetavec=PARAMHAT;
        thetavec(paramct)           = PARAMHAT_LB(paramct);
        for stepct=1:length(stepv)
            thetavec(4)             = stepv(stepct);
            sigvmatmom(:,ct,stepct) = SMM_bgms(thetavec);
        end

        %Solve at UB
        ct=ct+1;
        thetavec = PARAMHAT;
        thetavec(paramct)           = PARAMHAT_UB(paramct);
        for stepct=1:length(stepv)
            thetavec(4)             = stepv(stepct);
            sigvmatmom(:,ct,stepct) = SMM_bgms(thetavec);
        end 
    end

    save('figure2.mat','rhomatmom','sigwmatmom','sigematmom','sigvmatmom');

    %% q(1-p)
    pqmatmom = zeros(nbeta,9,5);

    %Solve at estimated values
    thetavec=PARAMHAT;
    thetavec(5) = 1.0;
    for stepct=1:length(stepp)
        thetavec(6) = stepp(stepct);
        pqmatmom(:,1,stepct) = SMM_bgms(thetavec);
    end

    ct=1;
    paramct=0;  
    while paramct<5

        paramct=paramct+1; 

        %Solve at LB
        ct=ct+1;
        thetavec=PARAMHAT;
        thetavec(5)=1.0;
        thetavec(paramct) = PARAMHAT_LB(paramct);
        for stepct=1:length(stepp)
            thetavec(6) = stepp(stepct);
            pqmatmom(:,ct,stepct) = SMM_bgms(thetavec);
        end

        %Solve at UB
        ct=ct+1;
        thetavec=PARAMHAT;
        thetavec(5)=1.0;
        thetavec(paramct) = PARAMHAT_UB(paramct);
        for stepct=1:length(stepp)
            thetavec(6) = stepp(stepct);
            pqmatmom(:,ct,stepct) = SMM_bgms(thetavec);
        end 

    end
    pqmatmom2=pqmatmom;
    pqmatmom2(1,:,1)=0;%adjusting for numerical error
    save('figure3.mat','pqmatmom2');

else
    
    load('figure2.mat');
    load('figure3.mat');
    
end
    
%% PLOT FIGURES 

% - Colors -
green = [0.0 0.6 0.0];
blue  = [0.0 0.0 0.6];
red   = [0.6 0.0 0.0];
myred= [0.6350 0.0780 0.1840];
myblue = [0 0.4470 0.7410];
grayColor = [.7 .7 .7];
% - Line width -
heavy = 3;
light = 0.8;

height = 11;
aspect = 16/9;



%% Figure 2
figure('units','centimeter','Position',[10 10 1.75*height*aspect 1.75*height],'color','w')
% PANEL 1
subplot(2,2,1)
hold on;
h1a=plot(stepr,squeeze(rhomatmom(1,1,:)),'LineWidth',3,'Color','k');
h1b=plot(stepr,squeeze(rhomatmom(1,2:8,:)),'LineWidth',1,'Color',grayColor);
h1c=plot(stepr,squeeze(rhomatmom(1,10:11,:)),'LineWidth',1,'Color',grayColor);
xlabel('$\rho$','interpreter','latex');
ylabel('BGMS coefficient','interpreter','latex');
set(gca,'FontSize',16,'TickLabelInterpreter', 'latex');
title('Panel 1','interpreter','latex');
hold off;
%PANEL 2
subplot(2,2,2)
hold on;
h1a=plot(stepw,squeeze(sigwmatmom(1,1,:)),'LineWidth',3,'Color','k');
h1b=plot(stepw,squeeze(sigwmatmom(1,2:8,:)),'LineWidth',1,'Color',grayColor);
h1c=plot(stepw,squeeze(sigwmatmom(1,10:11,:)),'LineWidth',1,'Color',grayColor);
xlabel('$\sigma_w$','interpreter','latex');
set(gca,'FontSize',16,'TickLabelInterpreter', 'latex');
title('Panel 2','interpreter','latex');
hold off;
%PANEL 3
subplot(2,2,3)
hold on;
h1a=plot(stepe,squeeze(sigematmom(1,1,:)),'LineWidth',3,'Color','k');
h1b=plot(stepe,squeeze(sigematmom(1,2:8,:)),'LineWidth',1,'Color',grayColor);
h1c=plot(stepe,squeeze(sigematmom(1,10:11,:)),'LineWidth',1,'Color',grayColor);
xlabel('$\sigma_e$','interpreter','latex');
ylabel('BGMS coefficient','interpreter','latex');
set(gca,'FontSize',16,'TickLabelInterpreter', 'latex');
title('Panel 3','interpreter','latex');
hold off;
%PANEL 4
subplot(2,2,4)
hold on;
h1a=plot(stepv,squeeze(sigvmatmom(1,1,:)),'LineWidth',3,'Color','k');
h1b=plot(stepv,squeeze(sigvmatmom(1,2:8,:)),'LineWidth',1,'Color',grayColor);
h1c=plot(stepv,squeeze(sigvmatmom(1,10:11,:)),'LineWidth',1,'Color',grayColor);
xlabel('$\sigma_v$','interpreter','latex');
set(gca,'FontSize',16,'TickLabelInterpreter', 'latex');
title('Panel 4','interpreter','latex');
hold off;
saveas(gcf,fullfile(figure_folder,'figure2.eps'),'epsc');

% Figure 3
figure('units','centimeter','Position',[10 10 1.75*height*aspect 1.75*height],'color','w')
hold on;
h1a=plot(flip(stepp),squeeze(pqmatmom2(1,1,:)),'LineWidth',3,'Color','k');
h1b=plot(flip(stepp),squeeze(pqmatmom2(1,2:8,:)),'LineWidth',1,'Color',grayColor);
h1c=plot(flip(stepp),squeeze(pqmatmom2(1,10:11,:)),'LineWidth',1,'Color',grayColor);
xlabel('$q(1-p)$','interpreter','latex');
ylabel('BGMS coefficient','interpreter','latex');
set(gca,'FontSize',16,'TickLabelInterpreter', 'latex');
hold off;
saveas(gcf,fullfile(figure_folder,'figure3.eps'),'epsc');


