%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% "The Distributional Effects of COVID-19 and Optimal Mitigation
%  Policies" (Sewon Hur, International Economic Review, 2023)
%
% run_01.m — reads the Fortran model output in results/ and produces:
%   * model figures (saved to figures/)
%   * the model inputs for Tables 3, 4, and 5 (printed to the console)
%   * results/model<jobnum>.xls, used by run_02_figures.do
%
% Settings:
%   jobnum   selects the experiment (1 = baseline; see list below)
%   projects selects which policy configurations to load; presets for
%            each table/figure are given below
%   nprint   1 = save figures as PDF, 0 = display only
%
% Run this from the package root (the folder containing results/).
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

clear;
close all;
addpath('auxiliary')

nprint = 1;
% 0: do not save figures
% 1: save figures as .pdf

runtran = 1;

jobnum = 1;
jobnumstr = sprintf('job%1i',jobnum);
if jobnum==0
    rootdir = 'results/nopandemic/';
elseif jobnum==1
    rootdir = 'results/baseline/';
elseif jobnum==2
    rootdir = 'results/2-taus/';
elseif jobnum==3
    rootdir = 'results/3-vsl/';
elseif jobnum==4
    rootdir = 'results/4-eta/';
elseif jobnum==5
    rootdir = 'results/5-kappa/';
elseif jobnum==6
    rootdir = 'results/6-pigou/';
elseif jobnum==7
    rootdir = 'results/7-gamma/';
elseif jobnum==8
    rootdir = 'results/8-utilde/';
elseif jobnum==9
    rootdir = 'results/9-r0/';
end

if(nprint == 0)
    set(0,'DefaultFigureWindowStyle','docked')
else
    set(0,'DefaultFigureWindowStyle','default')
end

% projects = [2 3 4 5]; % Tables 3 and 4
% projects = [2 9 6 7 8]; % Table 5
% projects = [2];       % Figure 4 panels
projects = [2 3];       % Figure 3 model series (writes results/model1.xls); also makes Figure 4 panels 
% 1: SIR (exogenous)
% 2: SIR Macro (no mitigation)
% 3: 2 + US policy
% 4: 2 +    US subsidy only
% 5: 2 +    US lockdown only
% 6: 2 + constrained optimum
% 7: 2 +    optimal subsidy only
% 8: 2 +    optimal tax only
% 9: 2 + output maximizing
% 10: 2 +    maximizing subsidy only
% 11: 2 +    maximizing tax only
Np = length(projects);

myPaperSize = [3.333333333333333 2.604166666666667];
lowerLeft = [2 2];
options.units = 'inches';
options.FontMode = 'Fixed';
options.FixedFontSize = 12;
options.format = 'pdf';
options.FontName = 'arial';
options.Renderer = 'painters';

for i_proj = 1:Np
    
    if projects(i_proj) == 1
        projectname{i_proj} = 'SIR';
        projectdir = strcat(rootdir,'exogenous/');
    elseif projects(i_proj) == 2
        if projects(1)==1
            projectname{i_proj} = 'SIR macro';
        else
            projectname{i_proj} = 'No public mitigation';
        end
        projectdir = strcat(rootdir,'no_mitigation/');
    elseif projects(i_proj) == 3
        projectname{i_proj} = 'US policy';
        projectdir = strcat(rootdir,'us_policy/');
    elseif projects(i_proj) == 4
        projectname{i_proj} = 'Subsidy only';
        projectdir = strcat(rootdir,'us_subsidy/');
    elseif projects(i_proj) == 5
        projectname{i_proj} = 'Lockdown only';
        projectdir = strcat(rootdir,'us_lockdown/');
    elseif projects(i_proj) == 6
        projectname{i_proj} = 'Constrained optimum';
        projectdir = strcat(rootdir,'optimum/');
    elseif projects(i_proj) == 7
        projectname{i_proj} = 'Optimal subsidy only';
        projectdir = strcat(rootdir,'optimum_sub/');
    elseif projects(i_proj) == 8
        projectname{i_proj} = 'Optimal tax only';
        projectdir = strcat(rootdir,'optimum_tax/');
    elseif projects(i_proj) == 9
        projectname{i_proj} = 'Output maximizing';
        projectdir = strcat(rootdir,'maximizing/');
    elseif projects(i_proj) == 10
        projectname{i_proj} = 'Maximizing subsidy only';
        projectdir = strcat(rootdir,'max_sub/');
    elseif projects(i_proj) == 11
        projectname{i_proj} = 'Maximizing tax only';
        projectdir = strcat(rootdir,'max_tax/');
    end
    if runtran==0
        projectname{i_proj} = 'US policy';
        projectdir = strcat(rootdir,'us_policy/');
    end
    
    %% READ DATA %%
    filename = strcat(projectdir,'PARAMS_MATLAB.out');
    params = importdata(filename,' ',1);
    days = params.data(1);
    N_k = params.data(2);
    N_z = params.data(3);
    jmax = params.data(4);
    beta_adj = params.data(5);
    beta_e{i_proj} = params.data(6);
    beta_c{i_proj} = params.data(7);
    beta_l{i_proj} = params.data(8);
    pi_x = params.data(9);
    hoscap = params.data(10);
    psi_j = params.data(11);
    initk = params.data(12);
    hikfrac = params.data(13);
    gam = params.data(19);
    clear params
    
    if i_proj==1
        
        filename = strcat(projectdir,'equilibrium_init.out');
        agg_init = importdata(filename,' ',1);
        wss = agg_init.data(1);
        rss = agg_init.data(2);
        aggk = agg_init.data(3);
        aggc = agg_init.data(4);
        aggy = agg_init.data(5);
        aggl = agg_init.data(6);
        aggh = agg_init.data(7);
        ss = agg_init.data(8);
        tau = agg_init.data(9);
        clear agg_init        
        
        
        % optimal saving kprime(j,z,k)
        filename = strcat(projectdir,'kss.out');
        formatSpec = '%18f%18f%18f%18f%18f%18f%18f%18f%18f%18f%18f%18f%18f%18f%18f%18f%18f%18f%18f%18f%18f%f%[^\n\r]';
        fileID = fopen(filename,'r');
        dataArray = textscan(fileID, formatSpec, 'Delimiter', '', 'WhiteSpace', '', 'TextType', 'string',  'ReturnOnError', false);
        fclose(fileID);
        vals = [dataArray{1:end-1}];
        
        kss(1,:,:)=reshape(vals(:,2:2+N_z-1)',1,N_z,N_k);
        kss(2,:,:)=reshape(vals(:,2+N_z:2+2*N_z-1)',1,N_z,N_k);
        kss(3,:,:)=reshape(vals(:,2+2*N_z:2+3*N_z-1)',1,N_z,N_k);
        
        clearvars filename formatSpec fileID dataArray ans vals;
        
        
        % consumption(j,z,k)
        filename = strcat(projectdir,'css.out');
        formatSpec = '%18f%18f%18f%18f%18f%18f%18f%18f%18f%18f%18f%18f%18f%18f%18f%18f%18f%18f%18f%18f%18f%f%[^\n\r]';
        fileID = fopen(filename,'r');
        dataArray = textscan(fileID, formatSpec, 'Delimiter', '', 'WhiteSpace', '', 'TextType', 'string',  'ReturnOnError', false);
        fclose(fileID);
        vals = [dataArray{1:end-1}];
        
        css(1,:,:)=reshape(vals(:,2:2+N_z-1)',1,N_z,N_k);
        css(2,:,:)=reshape(vals(:,2+N_z:2+2*N_z-1)',1,N_z,N_k);
        css(3,:,:)=reshape(vals(:,2+2*N_z:2+3*N_z-1)',1,N_z,N_k);
        
        clearvars filename formatSpec fileID dataArray ans vals;
        
        
        % hours(j,z,k)
        filename = strcat(projectdir,'lss.out');
        formatSpec = '%18f%18f%18f%18f%18f%18f%18f%18f%18f%18f%18f%18f%18f%18f%18f%18f%18f%18f%18f%18f%18f%f%[^\n\r]';
        fileID = fopen(filename,'r');
        dataArray = textscan(fileID, formatSpec, 'Delimiter', '', 'WhiteSpace', '', 'TextType', 'string',  'ReturnOnError', false);
        fclose(fileID);
        vals = [dataArray{1:end-1}];
        
        lss(1,:,:)=reshape(vals(:,2:2+N_z-1)',1,N_z,N_k);
        lss(2,:,:)=reshape(vals(:,2+N_z:2+2*N_z-1)',1,N_z,N_k);
        lss(3,:,:)=reshape(vals(:,2+2*N_z:2+3*N_z-1)',1,N_z,N_k);
        
        clearvars filename formatSpec fileID dataArray ans vals;
        
        
        % Savings rate
        filename = strcat(projectdir,'sss.out');
        formatSpec = '%18f%18f%18f%18f%18f%18f%18f%18f%18f%18f%18f%18f%18f%18f%18f%18f%18f%18f%18f%18f%18f%f%[^\n\r]';
        fileID = fopen(filename,'r');
        dataArray = textscan(fileID, formatSpec, 'Delimiter', '', 'WhiteSpace', '', 'TextType', 'string',  'ReturnOnError', false);
        fclose(fileID);
        vals = [dataArray{1:end-1}];
        
        sss(1,:,:)=reshape(vals(:,2:2+N_z-1)',1,N_z,N_k);
        sss(2,:,:)=reshape(vals(:,2+N_z:2+2*N_z-1)',1,N_z,N_k);
        sss(3,:,:)=reshape(vals(:,2+2*N_z:2+3*N_z-1)',1,N_z,N_k);
        
        clearvars filename formatSpec fileID dataArray ans vals;
        
        % wealth distribution(j,z,k)
        filename = strcat(projectdir,'muss.out');
        formatSpec = '%18f%18f%18f%18f%18f%18f%18f%18f%18f%18f%18f%18f%18f%18f%18f%18f%18f%18f%18f%18f%18f%f%[^\n\r]';
        fileID = fopen(filename,'r');
        dataArray = textscan(fileID, formatSpec, 'Delimiter', '', 'WhiteSpace', '', 'TextType', 'string',  'ReturnOnError', false);
        fclose(fileID);
        vals = [dataArray{1:end-1}];
        
        muss(1,:,:)=reshape(vals(:,2:2+N_z-1)',1,N_z,N_k);
        muss(2,:,:)=reshape(vals(:,2+N_z:2+2*N_z-1)',1,N_z,N_k);
        muss(3,:,:)=reshape(vals(:,2+2*N_z:2+3*N_z-1)',1,N_z,N_k);
        
        clearvars filename formatSpec fileID dataArray ans vals;
        
        % vsl(j,z,k)
        filename = strcat(projectdir,'v_vsl.out');
        formatSpec = '%18f%18f%18f%18f%18f%18f%18f%18f%18f%18f%18f%18f%18f%18f%18f%18f%18f%18f%18f%18f%18f%f%[^\n\r]';
        fileID = fopen(filename,'r');
        dataArray = textscan(fileID, formatSpec, 'Delimiter', '', 'WhiteSpace', '', 'TextType', 'string',  'ReturnOnError', false);
        fclose(fileID);
        vals = [dataArray{1:end-1}];
        
        vsl(1,:,:)=reshape(vals(:,2:2+N_z-1)',1,N_z,N_k);
        vsl(2,:,:)=reshape(vals(:,2+N_z:2+2*N_z-1)',1,N_z,N_k);
        vsl(3,:,:)=reshape(vals(:,2+2*N_z:2+3*N_z-1)',1,N_z,N_k);
        
        clearvars filename formatSpec fileID dataArray ans vals;
        
        
        % k_grid
        filename = strcat(projectdir,'k_grid.out');
        k_grid=dlmread(filename);
        
        [i,j]=size(k_grid);
        k_grid=reshape(k_grid',i*j,1);      % transpose array into column vector
        k_grid=k_grid(1:N_k);               % store only data (drop empty terms)
        [~,ik_zero] = min(abs(k_grid));
        
        
        % z_grid
        filename = strcat(projectdir,'z_grid.out');
        z_grid=dlmread(filename);
        [i,j]=size(z_grid);
        z_grid=reshape(z_grid',i*j,1);      % transpose array into column vector
        z_grid=z_grid(1:N_z);               % store only data (drop empty terms)
        
        
        % age-income profile eta_j
        filename = strcat(projectdir,'eta_j.out');
        eta_j=dlmread(filename);
        [i,j]=size(eta_j);
        eta_j=reshape(eta_j',i*j,1);        % transpose array into column vector
        eta_j=eta_j(1:jmax);                % store only data (drop empty terms)
        
        
        % ptran(z,z)
        filename = strcat(projectdir,'ptran.out');
        ptran=dlmread(filename);
        [i,j]=size(ptran);
        ptran=reshape(ptran',i*j,1);        % transpose array into column vector
        ptran=ptran(1:N_z*N_z);             % store only data (drop empty terms)
        ptran=reshape(ptran,N_z,N_z);       % reshape vector into matrix
        
        
        % Erd(z)
        filename = strcat(projectdir,'Erd.out');
        Erd=dlmread(filename);
        [i,j]=size(Erd);
        Erd=reshape(Erd',i*j,1);            % transpose array into column vector
        Erd=Erd(1:N_z);                     % store only data (drop empty terms)
        
        
        % wfh_z(z)
        filename = strcat(projectdir,'wfh_z.out');
        wfh_z=dlmread(filename);
        [i,j]=size(wfh_z);
        wfh_z=reshape(wfh_z',i*j,1);        % transpose array into column vector
        wfh_z=wfh_z(1:(jmax-1)*N_z);	    % store only data (drop empty terms)
        wfh_z=reshape(wfh_z,(jmax-1),N_z); 

        % labor income per working age person
        lincomewap = wss*aggl/sum(sum(sum(muss(1:2,:,:))));
        lincomewap_ann = lincomewap*(365/days);

        
        %% Calibration statistics        
        
        % Get some ginis (individuals, ages 25+)
        wealth = reshape(repmat(reshape(k_grid,[1 1 N_k]),[jmax N_z 1]),jmax*N_z*N_k,1);
        pop = reshape(muss,jmax*N_z*N_k,1);
        consump = reshape(css,jmax*N_z*N_k,1);
        earnings = NaN(jmax,N_z,N_k);
        earnings(1,:,:) = wss*(1-tau)*eta_j(1)*repmat(z_grid,[1,N_k]).*squeeze(lss(1,:,:));
        earnings(2,:,:) = wss*(1-tau)*eta_j(2)*repmat(z_grid,[1,N_k]).*squeeze(lss(2,:,:));
        earnings(3,:,:) = ss;
        
        gini_k = gini(pop,wealth);
        gini_c = gini(pop,consump); 
        gini_y = gini(pop,earnings);        
               
        cumpop = squeeze(sum(muss,[1,2]));
        cumpop = cumsum(cumpop);
        for i=2:N_k
            if cumpop(i)<=cumpop(i-1)
                cumpop(i) = cumpop(i-1)+eps;
            end
        end
        ik25 = interp1(cumpop,1:1:N_k,.25);
        ik50 = interp1(cumpop,1:1:N_k,.5);
        ik75 = interp1(cumpop,1:1:N_k,.75);
        k25 = k_grid(floor(ik25))*(ceil(ik25)-ik25) + k_grid(ceil(ik25))*(ik25-floor(ik25));
        k75 = k_grid(floor(ik75))*(ceil(ik75)-ik75) + k_grid(ceil(ik75))*(ik75-floor(ik75));
        k7525 = k75/k25;     
        
        fprintf('\n<<inputs for TABLE 3>>')
        fprintf('\nCapital / annual GDP: %3.1f', aggk/(aggy*365/days) );
        fprintf('\nHours per worker: %2.1f', aggh/sum(sum(sum(muss(1:2,:,:))))*112 );
        fprintf('\nretirement income (percent of earnings/worker): %3.1f', 100*ss/lincomewap )
        fprintf('\ndisposable earnings gini: %3.2f', gini_y );
        fprintf('\nconsumption gini: %3.2f', gini_c );
        fprintf('\nwealth gini: %3.2f', gini_k );
        fprintf('\nwealth p75/p25: %3.1f\n', k7525 );

        clear consump earnings wealth pop   

        ik50lo = floor(ik50);
        ik50hi = ik50lo+1;
        
        vsl_lzlk = (sum(sum(sum( ...
            muss(1,1:4,1:ik50lo).*vsl(1,1:4,1:ik50lo) )))...
            + sum(sum(sum( ...
            muss(2,1:3,1:ik50lo).*vsl(2,1:3,1:ik50lo) ))) ) ...
            /(sum(sum(sum(muss(1,1:4,1:ik50lo) ))) ...
            +sum(sum(sum(muss(2,1:3,1:ik50lo) ))) );
        vsl_hzhk = (sum(sum(sum( ...
            muss(1,5:7,ik50hi:N_k).*vsl(1,5:7,ik50hi:N_k) )))...
            + sum(sum(sum( ...
            + muss(2,4:7,ik50hi:N_k).*vsl(2,4:7,ik50hi:N_k) ))) ) ...
            /(sum(sum(sum(muss(1,5:7,ik50hi:N_k) ))) ...
            +sum(sum(sum(muss(2,4:7,ik50hi:N_k) ))) );
        
        fprintf('\nvsl (low wage/low wealth): %3.0f', vsl_lzlk/(365/14) );
        fprintf('\nvsl (high wage/high wealth): %3.0f', vsl_hzhk/(365/14) );

        
    end
    
    %% TRANSITION %%
    if runtran
        
        filename = strcat(projectdir,'PARAMS_MATLAB.out');
        params = importdata(filename,' ',1);
        if length(params.data)>13
            sui_cbar{i_proj} = params.data(14);
            hcut{i_proj} = params.data(15);
            mitdur{i_proj} = params.data(16);
            lockdown{i_proj} = params.data(17);
            hcap{i_proj} = params.data(18);    
        end
        clear params  
        
        
        filename = strcat(projectdir,'equilibrium_tran.out');
        startRow = 2;
        formatSpec = '%5f%15f%15f%15f%15f%15f%15f%15f%15f%15f%15f%15f%15f%15f%15f%15f%15f%15f%15f%15f%15f%15f%15f%15f%15f%15f%15f%15f%15f%15f%15f%15f%15f%15f%15f%15f%15f%15f%15f%15f%15f%15f%15f%15f%15f%15f%15f%15f%15f%15f%15f%15f%15f%15f%15f%15f%15f%15f%15f%15f%15f%15f%15f%15f%15f%15f%15f%15f%15f%15f%15f%15f%15f%f%[^\n\r]';
        fileID = fopen(filename,'r');
        dataArray = textscan(fileID, formatSpec, 'Delimiter', '', 'WhiteSpace', '', 'TextType', 'string', 'HeaderLines' ,startRow-1, 'ReturnOnError', false, 'EndOfLine', '\r\n');
        fclose(fileID);
        vals = [dataArray{1:end-1}];
        
        wpath{i_proj} = vals(:,2);
        rpath{i_proj} = vals(:,3);
        kpath{i_proj} = vals(:,4);
        cpath{i_proj} = vals(:,5);
        ypath{i_proj} = vals(:,6);
        lpath{i_proj} = vals(:,7);
        hourspath{i_proj} = vals(:,8);
        bigspath{i_proj} = vals(:,9);
        bigipath{i_proj} = vals(:,10);
        bigrpath{i_proj} = vals(:,11);
        s1path{i_proj} = vals(:,12);
        i1path{i_proj} = vals(:,13);
        r1path{i_proj} = vals(:,14);
        s2path{i_proj} = vals(:,15);
        i2path{i_proj} = vals(:,16);
        r2path{i_proj} = vals(:,17);
        s3path{i_proj} = vals(:,18);
        i3path{i_proj} = vals(:,19);
        r3path{i_proj} = vals(:,20);
        c1path{i_proj} = vals(:,21);
        h1path{i_proj} = vals(:,22);
        c2path{i_proj} = vals(:,23);
        h2path{i_proj} = vals(:,24);
        c3path{i_proj} = vals(:,25);
        h3path{i_proj} = vals(:,26);
        c1spath{i_proj} = vals(:,27);
        c2spath{i_proj} = vals(:,28);
        c3spath{i_proj} = vals(:,29);
        c1ipath{i_proj} = vals(:,30);
        c2ipath{i_proj} = vals(:,31);
        c3ipath{i_proj} = vals(:,32);
        c1rpath{i_proj} = vals(:,33);
        c2rpath{i_proj} = vals(:,34);
        c3rpath{i_proj} = vals(:,35);
        oh1spath{i_proj} = vals(:,36);
        oh2spath{i_proj} = vals(:,37);
        oh3spath{i_proj} = vals(:,38);
        oh1ipath{i_proj} = vals(:,39);
        oh2ipath{i_proj} = vals(:,40);
        oh3ipath{i_proj} = vals(:,41);
        oh1rpath{i_proj} = vals(:,42);
        oh2rpath{i_proj} = vals(:,43);
        oh3rpath{i_proj} = vals(:,44);
        htot1spath{i_proj} = vals(:,45);
        htot2spath{i_proj} = vals(:,46);
        htot3spath{i_proj} = vals(:,47);
        htot1ipath{i_proj} = vals(:,48);
        htot2ipath{i_proj} = vals(:,49);
        htot3ipath{i_proj} = vals(:,50);
        htot1rpath{i_proj} = vals(:,51);
        htot2rpath{i_proj} = vals(:,52);
        htot3rpath{i_proj} = vals(:,53);
        dpath{i_proj} = vals(:,54);
        aggtaxpath{i_proj} = vals(:,55);
        aggsuipath{i_proj} = vals(:,56);
        pvdef{i_proj} = vals(1,57);
        oc1spath{i_proj} = vals(:,58);
        oc2spath{i_proj} = vals(:,59);
        oc3spath{i_proj} = vals(:,60);
        oc1ipath{i_proj} = vals(:,61);
        oc2ipath{i_proj} = vals(:,62);
        oc3ipath{i_proj} = vals(:,63);
        oc1rpath{i_proj} = vals(:,64);
        oc2rpath{i_proj} = vals(:,65);
        oc3rpath{i_proj} = vals(:,66);
        takeuppath{i_proj} = vals(:,67);
        betaepath{i_proj} = vals(:,68);
        d1path{i_proj} = vals(:,69);
        d2path{i_proj} = vals(:,70);
        d3path{i_proj} = vals(:,71);
        clearvars filename startRow formatSpec fileID dataArray ans vals;
                
        
        if aggtaxpath{i_proj}(2)>0
            tauc{i_proj}=aggtaxpath{i_proj}(1);
        else
            tauc{i_proj}=0;
        end
        
        if aggsuipath{i_proj}(2)>0
            taus{i_proj}=aggsuipath{i_proj}(1);
        else
            taus{i_proj}=0;
        end
        
        htot1path{i_proj} = htot1spath{i_proj} ...
            + htot1ipath{i_proj} ...
            + htot1rpath{i_proj};
        htot2path{i_proj} = htot2spath{i_proj} ...
            + htot2ipath{i_proj} ...
            + htot2rpath{i_proj};
        hourstotpath{i_proj} = htot1path{i_proj} + htot2path{i_proj};
        
        oh1path{i_proj} = oh1spath{i_proj} ...
            + oh1ipath{i_proj} ...
            + oh1rpath{i_proj};
        oh2path{i_proj} = oh2spath{i_proj} ...
            + oh2ipath{i_proj} ...
            + oh2rpath{i_proj};
        if projects(i_proj) == 1
            ohpath{i_proj} = hourstotpath{i_proj};
        else
            ohpath{i_proj} = oh1path{i_proj} + oh2path{i_proj};
        end
        
        oc1path{i_proj} = oc1spath{i_proj} ...
            + oc1ipath{i_proj} ...
            + oc1rpath{i_proj};
        oc2path{i_proj} = oc2spath{i_proj} ...
            + oc2ipath{i_proj} ...
            + oc2rpath{i_proj};
        oc3path{i_proj} = oc3spath{i_proj} ...
            + oc3ipath{i_proj} ...
            + oc3rpath{i_proj};
        if projects(i_proj) == 1
            ocpath{i_proj} = cpath{i_proj};
        else
            ocpath{i_proj} = oc1path{i_proj} + oc2path{i_proj} + oc3path{i_proj};
        end
        
        repnum{i_proj} = NaN*bigipath{i_proj};
        if projects(i_proj)==1 % exogenous
            repnum{i_proj}(1) = beta_e{i_proj}*bigspath{i_proj}(1)/pi_x;
            for t = 1:length(repnum{i_proj})
                if bigipath{i_proj}(t)>0
                    repnum{i_proj}(t) = beta_e{i_proj}*bigspath{i_proj}(t)/pi_x;
                end
            end
        else % SIR Macro
            repnum{i_proj}(1) = beta_c{i_proj}*(oc1spath{i_proj}(1)+oc2spath{i_proj}(1)+oc3spath{i_proj}(1))^2/pi_x ...
                + beta_l{i_proj}*(oh1spath{i_proj}(1)+oh2spath{i_proj}(1))^2/pi_x ...
                + beta_e{i_proj}*bigspath{i_proj}(1)/pi_x;
            for t = 1:length(repnum{i_proj})
                if bigipath{i_proj}(t)>=0.001
                    repnum{i_proj}(t) = ...
                        beta_c{i_proj}*(oc1spath{i_proj}(t)+oc2spath{i_proj}(t)+oc3spath{i_proj}(t)).*(oc1ipath{i_proj}(t)+oc2ipath{i_proj}(t)+oc3ipath{i_proj}(t))/pi_x./bigipath{i_proj}(t) ...
                        + beta_l{i_proj}*(oh1spath{i_proj}(t)+oh2spath{i_proj}(t)).*(oh1ipath{i_proj}(t)+oh2ipath{i_proj}(t))/pi_x./bigipath{i_proj}(t) ...
                        + betaepath{i_proj}(t)*bigspath{i_proj}(t)/pi_x;
                end
            end
        end
        
        tcut = floor(2*365/days);
        yloss{i_proj} = 1-sum(ypath{i_proj}(1:tcut))/(aggy*tcut);
        yindex{i_proj} = 100*sum(ypath{i_proj}(1:tcut))/(aggy*tcut);
        
        ocipath{i_proj} = oc1ipath{i_proj}+oc2ipath{i_proj}+oc3ipath{i_proj};
        
        c1spath{i_proj} = c1spath{i_proj}./s1path{i_proj};
        c1ipath{i_proj} = c1ipath{i_proj}./i1path{i_proj};
        c1rpath{i_proj} = c1rpath{i_proj}./r1path{i_proj};
        c2spath{i_proj} = c2spath{i_proj}./s2path{i_proj};
        c2ipath{i_proj} = c2ipath{i_proj}./i2path{i_proj};
        c2rpath{i_proj} = c2rpath{i_proj}./r2path{i_proj};
        c3spath{i_proj} = c3spath{i_proj}./s3path{i_proj};
        c3ipath{i_proj} = c3ipath{i_proj}./i3path{i_proj};
        c3rpath{i_proj} = c3rpath{i_proj}./r3path{i_proj};
        
        oh1spath{i_proj} = oh1spath{i_proj}./s1path{i_proj};
        oh1ipath{i_proj} = oh1ipath{i_proj}./i1path{i_proj};
        oh1rpath{i_proj} = oh1rpath{i_proj}./r1path{i_proj};
        oh2spath{i_proj} = oh2spath{i_proj}./s2path{i_proj};
        oh2ipath{i_proj} = oh2ipath{i_proj}./i2path{i_proj};
        oh2rpath{i_proj} = oh2rpath{i_proj}./r2path{i_proj};
        oh3spath{i_proj} = oh3spath{i_proj}./s3path{i_proj};
        oh3ipath{i_proj} = oh3ipath{i_proj}./i3path{i_proj};
        oh3rpath{i_proj} = oh3rpath{i_proj}./r3path{i_proj};
        
        htot1spath{i_proj} = htot1spath{i_proj}./s1path{i_proj};
        htot1ipath{i_proj} = htot1ipath{i_proj}./i1path{i_proj};
        htot1rpath{i_proj} = htot1rpath{i_proj}./r1path{i_proj};
        htot2spath{i_proj} = htot2spath{i_proj}./s2path{i_proj};
        htot2ipath{i_proj} = htot2ipath{i_proj}./i2path{i_proj};
        htot2rpath{i_proj} = htot2rpath{i_proj}./r2path{i_proj};
        htot3spath{i_proj} = htot3spath{i_proj}./s3path{i_proj};
        htot3ipath{i_proj} = htot3ipath{i_proj}./i3path{i_proj};
        htot3rpath{i_proj} = htot3rpath{i_proj}./r3path{i_proj};
        
        poppath{i_proj} = bigipath{i_proj}+bigspath{i_proj}+bigrpath{i_proj};
        cumdpath{i_proj} = cumsum(dpath{i_proj});

        filename = strcat(projectdir,'equilibrium_tran_2.out');
        startRow = 2;
        formatSpec = '%5f%15f%15f%15f%15f%15f%15f%15f%15f%15f%15f%15f%15f%15f%15f%15f%15f%15f%f%[^\n\r]';
        fileID = fopen(filename,'r');
        dataArray = textscan(fileID, formatSpec, 'Delimiter', '', 'WhiteSpace', '', 'TextType', 'string', 'HeaderLines' ,startRow-1, 'ReturnOnError', false, 'EndOfLine', '\r\n');
        fclose(fileID);
        vals = [dataArray{1:end-1}];
        
        oh_y_lzlkpath{i_proj} = vals(:,2);
        oh_y_hzlkpath{i_proj} = vals(:,3);
        oh_y_lzhkpath{i_proj} = vals(:,4);
        oh_y_hzhkpath{i_proj} = vals(:,5);
        oh_m_lzlkpath{i_proj} = vals(:,6);
        oh_m_hzlkpath{i_proj} = vals(:,7);
        oh_m_lzhkpath{i_proj} = vals(:,8);
        oh_m_hzhkpath{i_proj} = vals(:,9);
        oc_y_lzlkpath{i_proj} = vals(:,10);
        oc_y_hzlkpath{i_proj} = vals(:,11);
        oc_y_lzhkpath{i_proj} = vals(:,12);
        oc_y_hzhkpath{i_proj} = vals(:,13);
        oc_m_lzlkpath{i_proj} = vals(:,14);
        oc_m_hzlkpath{i_proj} = vals(:,15);
        oc_m_lzhkpath{i_proj} = vals(:,16);
        oc_m_hzhkpath{i_proj} = vals(:,17);
        oc_r_lkpath{i_proj} = vals(:,18);
        oc_r_hkpath{i_proj} = vals(:,19);
        clearvars filename startRow formatSpec fileID dataArray ans vals;

        % welf(j,z,k)
        filename = strcat(projectdir,'welf.out');
        formatSpec = '%18f%18f%18f%18f%18f%18f%18f%18f%18f%18f%18f%18f%18f%18f%18f%18f%18f%18f%18f%18f%18f%f%[^\n\r]';
        fileID = fopen(filename,'r');
        dataArray = textscan(fileID, formatSpec, 'Delimiter', '', 'WhiteSpace', '', 'TextType', 'string',  'ReturnOnError', false);
        fclose(fileID);
        vals = [dataArray{1:end-1}];
        
        welf{i_proj}(1,:,:)=reshape(vals(:,2:2+N_z-1)',1,N_z,N_k);
        welf{i_proj}(2,:,:)=reshape(vals(:,2+N_z:2+2*N_z-1)',1,N_z,N_k);
        welf{i_proj}(3,:,:)=reshape(vals(:,2+2*N_z:2+3*N_z-1)',1,N_z,N_k);
        
        clearvars filename formatSpec fileID dataArray ans vals;
        if isfile(strcat(projectdir,'vt_s.out'))
            
            % vt_s(j,z,k) (susceptible)
            filename = strcat(projectdir,'vt_s.out');
            formatSpec = '%5f%18f%18f%18f%18f%18f%18f%18f%18f%18f%18f%18f%18f%18f%18f%18f%18f%18f%18f%18f%18f%18f%f%[^\n\r]';
            fileID = fopen(filename,'r');
            dataArray = textscan(fileID, formatSpec, 'Delimiter', '', 'WhiteSpace', '', 'TextType', 'string',  'ReturnOnError', false);
            fclose(fileID);
            vals = [dataArray{1:end-1}];

            [i,~] =size(vals);
            Tend = i/N_k;
            for t=1:Tend
                for idx_k=1:N_k
                    i = (t-1)*N_k+idx_k;
                    vt_s{i_proj}(1,:,idx_k,t) = vals(i,3:3+N_z-1);
                    vt_s{i_proj}(2,:,idx_k,t) = vals(i,3+N_z:3+2*N_z-1);
                    vt_s{i_proj}(3,:,idx_k,t) = vals(i,3+2*N_z:3+3*N_z-1);
                end
            end

            clearvars filename formatSpec fileID dataArray ans vals;
            

            % vt_i(j,z,k) (infected)
            filename = strcat(projectdir,'vt_i.out');
            formatSpec = '%5f%18f%18f%18f%18f%18f%18f%18f%18f%18f%18f%18f%18f%18f%18f%18f%18f%18f%18f%18f%18f%18f%f%[^\n\r]';
            fileID = fopen(filename,'r');
            dataArray = textscan(fileID, formatSpec, 'Delimiter', '', 'WhiteSpace', '', 'TextType', 'string',  'ReturnOnError', false);
            fclose(fileID);
            vals = [dataArray{1:end-1}];

            [i,~] =size(vals);
            Tend = i/N_k;
            for t=1:Tend
                for idx_k=1:N_k
                    i = (t-1)*N_k+idx_k;
                    vt_i{i_proj}(1,:,idx_k,t) = vals(i,3:3+N_z-1);
                    vt_i{i_proj}(2,:,idx_k,t) = vals(i,3+N_z:3+2*N_z-1);
                    vt_i{i_proj}(3,:,idx_k,t) = vals(i,3+2*N_z:3+3*N_z-1);
                end
            end

            clearvars filename formatSpec fileID dataArray ans vals;
            
        
            % consumption(j,z,k,susceptible,t)
            filename = strcat(projectdir,'ct_s.out');
            formatSpec = '%5f%18f%18f%18f%18f%18f%18f%18f%18f%18f%18f%18f%18f%18f%18f%18f%18f%18f%18f%18f%18f%18f%f%[^\n\r]';
            fileID = fopen(filename,'r');
            dataArray = textscan(fileID, formatSpec, 'Delimiter', '', 'WhiteSpace', '', 'TextType', 'string',  'ReturnOnError', false);
            fclose(fileID);
            vals = [dataArray{1:end-1}];

            [i,~] =size(vals);
            Tend = i/N_k;
            for t=1:Tend
                for idx_k=1:N_k
                    i = (t-1)*N_k+idx_k;
                    ct_s{i_proj}(1,:,idx_k,t) = vals(i,3:3+N_z-1);
                    ct_s{i_proj}(2,:,idx_k,t) = vals(i,3+N_z:3+2*N_z-1);
                    ct_s{i_proj}(3,:,idx_k,t) = vals(i,3+2*N_z:3+3*N_z-1);
                end
            end
            clearvars filename formatSpec fileID dataArray ans vals;
            

            % cfht(j,z,k,susceptible,t)
            filename = strcat(projectdir,'cfht_s.out');
            formatSpec = '%5f%18f%18f%18f%18f%18f%18f%18f%18f%18f%18f%18f%18f%18f%18f%18f%18f%18f%18f%18f%18f%18f%f%[^\n\r]';
            fileID = fopen(filename,'r');
            dataArray = textscan(fileID, formatSpec, 'Delimiter', '', 'WhiteSpace', '', 'TextType', 'string',  'ReturnOnError', false);
            fclose(fileID);
            vals = [dataArray{1:end-1}];

            [i,~] =size(vals);
            Tend = i/N_k;
            for t=1:Tend
                for idx_k=1:N_k
                    i = (t-1)*N_k+idx_k;
                    cfht_s{i_proj}(1,:,idx_k,t) = vals(i,3:3+N_z-1);
                    cfht_s{i_proj}(2,:,idx_k,t) = vals(i,3+N_z:3+2*N_z-1);
                    cfht_s{i_proj}(3,:,idx_k,t) = vals(i,3+2*N_z:3+3*N_z-1);
                end
            end
            clearvars filename formatSpec fileID dataArray ans vals;
            

            % consumption(j,z,k,infected,t)
            filename = strcat(projectdir,'ct_i.out');
            formatSpec = '%5f%18f%18f%18f%18f%18f%18f%18f%18f%18f%18f%18f%18f%18f%18f%18f%18f%18f%18f%18f%18f%18f%f%[^\n\r]';
            fileID = fopen(filename,'r');
            dataArray = textscan(fileID, formatSpec, 'Delimiter', '', 'WhiteSpace', '', 'TextType', 'string',  'ReturnOnError', false);
            fclose(fileID);
            vals = [dataArray{1:end-1}];

            [i,~] =size(vals);
            Tend = i/N_k;
            for t=1:Tend
                for idx_k=1:N_k
                    i = (t-1)*N_k+idx_k;
                    ct_i{i_proj}(1,:,idx_k,t) = vals(i,3:3+N_z-1);
                    ct_i{i_proj}(2,:,idx_k,t) = vals(i,3+N_z:3+2*N_z-1);
                    ct_i{i_proj}(3,:,idx_k,t) = vals(i,3+2*N_z:3+3*N_z-1);
                end
            end
            clearvars filename formatSpec fileID dataArray ans vals;

            
            % hours(j,z,k,susceptible,t)
            filename = strcat(projectdir,'lt_s.out');
            formatSpec = '%5f%18f%18f%18f%18f%18f%18f%18f%18f%18f%18f%18f%18f%18f%18f%18f%18f%18f%18f%18f%18f%18f%f%[^\n\r]';
            fileID = fopen(filename,'r');
            dataArray = textscan(fileID, formatSpec, 'Delimiter', '', 'WhiteSpace', '', 'TextType', 'string',  'ReturnOnError', false);
            fclose(fileID);
            vals = [dataArray{1:end-1}];

            [i,~] =size(vals);
            Tend = i/N_k;
            for t=1:Tend
                for idx_k=1:N_k
                    i = (t-1)*N_k+idx_k;
                    lt_s{i_proj}(1,:,idx_k,t) = vals(i,3:3+N_z-1);
                    lt_s{i_proj}(2,:,idx_k,t) = vals(i,3+N_z:3+2*N_z-1);
                    lt_s{i_proj}(3,:,idx_k,t) = vals(i,3+2*N_z:3+3*N_z-1);
                end
            end
            clearvars filename formatSpec fileID dataArray ans vals;

            
            % hours(j,z,k,infected,t)
            filename = strcat(projectdir,'lt_i.out');
            formatSpec = '%5f%18f%18f%18f%18f%18f%18f%18f%18f%18f%18f%18f%18f%18f%18f%18f%18f%18f%18f%18f%18f%18f%f%[^\n\r]';
            fileID = fopen(filename,'r');
            dataArray = textscan(fileID, formatSpec, 'Delimiter', '', 'WhiteSpace', '', 'TextType', 'string',  'ReturnOnError', false);
            fclose(fileID);
            vals = [dataArray{1:end-1}];

            [i,~] =size(vals);
            Tend = i/N_k;
            for t=1:Tend
                for idx_k=1:N_k
                    i = (t-1)*N_k+idx_k;
                    lt_i{i_proj}(1,:,idx_k,t) = vals(i,3:3+N_z-1);
                    lt_i{i_proj}(2,:,idx_k,t) = vals(i,3+N_z:3+2*N_z-1);
                    lt_i{i_proj}(3,:,idx_k,t) = vals(i,3+2*N_z:3+3*N_z-1);
                end
            end
            clearvars filename formatSpec fileID dataArray ans vals;

            
            % hours(j,z,k,recovered,t)
            filename = strcat(projectdir,'lt_r.out');
            formatSpec = '%5f%18f%18f%18f%18f%18f%18f%18f%18f%18f%18f%18f%18f%18f%18f%18f%18f%18f%18f%18f%18f%18f%f%[^\n\r]';
            fileID = fopen(filename,'r');
            dataArray = textscan(fileID, formatSpec, 'Delimiter', '', 'WhiteSpace', '', 'TextType', 'string',  'ReturnOnError', false);
            fclose(fileID);
            vals = [dataArray{1:end-1}];

            [i,~] =size(vals);
            Tend = i/N_k;
            for t=1:Tend
                for idx_k=1:N_k
                    i = (t-1)*N_k+idx_k;
                    lt_r{i_proj}(1,:,idx_k,t) = vals(i,3:3+N_z-1);
                    lt_r{i_proj}(2,:,idx_k,t) = vals(i,3+N_z:3+2*N_z-1);
                    lt_r{i_proj}(3,:,idx_k,t) = vals(i,3+2*N_z:3+3*N_z-1);
                end
            end
            clearvars filename formatSpec fileID dataArray ans vals;


            hcappath{i_proj} = 112*hcap{i_proj}*ones(Tend,1);
            for t=4:Tend
                hcappath{i_proj}(t) = 112*hcap{i_proj}* ...
                    ((365*2-t*14)/(365*2-42))^2 ...
                    + 112*(1-((365*2-t*14)/(365*2-42))^2);
            end

            lt_max = max(max(lt_s{i_proj},lt_i{i_proj}),lt_r{i_proj});
            [hcappath{i_proj},...
                112*squeeze(max(max(max(...
                lt_max(1:2,:,:,:).*repmat(1-wfh_z,[1,1,N_k]),[],1),[],2),[],3))];
            
        end % if isfile(strcat(projectdir,'vt_s.out'))
        
        
        ik50lo = floor(ik50);
        ik50hi = ik50lo+1;
        
        old_lk{i_proj} = 100*sum(sum(sum(muss(3,1,1:ik50lo)...
            .*welf{i_proj}(3,1,1:ik50lo))))/sum(sum(sum(muss(3,1,1:ik50lo))));
        old_hk{i_proj} = 100*sum(sum(sum(muss(3,1,ik50hi:N_k)...
            .*welf{i_proj}(3,1,ik50hi:N_k))))/sum(sum(sum(muss(3,1,ik50hi:N_k))));
        
        working_age{i_proj} = 100*sum(sum(sum(muss(1:2,:,:)...
            .*welf{i_proj}(1:2,:,:))))/sum(sum(sum(muss(1:2,:,:))));      
        retired{i_proj} = 100*sum(sum(muss(3,:,:)...
            .*welf{i_proj}(3,:,:)))/sum(sum(muss(3,:,:)));
        low_z{i_proj} = 100*(sum(sum( ...
            muss(1,1:4,:).*welf{i_proj}(1,1:4,:))) ...
            + sum(sum( ...
            muss(2,1:3,:).*welf{i_proj}(2,1:3,:))) ) ...
            /(sum(sum(muss(1,1:4,:))) + sum(sum(muss(2,1:3,:))) );
        high_z{i_proj} = 100*(sum(sum( ...
            muss(1,5:7,:).*welf{i_proj}(1,5:7,:))) ...
            + sum(sum( ...
            muss(2,4:7,:).*welf{i_proj}(2,4:7,:))) ) ...
            /(sum(sum(muss(1,5:7,:))) + sum(sum(muss(2,4:7,:))) );
        low_k{i_proj} = 100*sum(sum(sum(muss(:,:,1:ik50lo)...
            .*welf{i_proj}(:,:,1:ik50lo))))/sum(sum(sum(muss(:,:,1:ik50lo))));
        high_k{i_proj} = 100*sum(sum(sum(muss(:,:,ik50hi:N_k)...
            .*welf{i_proj}(:,:,ik50hi:N_k))))/sum(sum(sum(muss(:,:,ik50hi:N_k))));
        
        avgwelf{i_proj} = 100*sum(sum(sum(muss.*welf{i_proj})))/sum(sum(sum(muss)));
        if projects(1) == 2
            poswelf{i_proj} = 100*sum(sum(sum(muss.*(welf{i_proj}-welf{1}>0))))/sum(sum(sum(muss)));
        end
        
        popadj = 217229155/329135084;
        peak_inf{i_proj} = max(bigipath{i_proj});
        deaths{i_proj} = 217229155*max(cumdpath{i_proj});
        death_per_10k{i_proj} = max(cumdpath{i_proj})*10000*popadj;
		
		dfrac_young{i_proj} = sum(d1path{i_proj})/max(cumdpath{i_proj})*100;
		dfrac_mid{i_proj} = sum(d2path{i_proj})/max(cumdpath{i_proj})*100;

        
        tstart = 2020 + 72/365;
        tgrid = [tstart:days/365:tstart+days/365*(length(kpath{i_proj})-1)];
         
        %% Figure 4 (response to pandemic, no public mitigation) %%
        
        iz10 = 2;
        iz25 = 3;
        iz75 = 5;
        iz90 = 6;

        if isfile(strcat(projectdir,'vt_s.out'))
           
            % (dead plotting code removed for the replication package)
%             xlabel('time')
%             xlim([2020.15 2022])
%             ylim([0 120])
%             xticks([2020:0.25:2022])
%             xticklabels({'','2020q2','','2020q4','','2021q2','','2021q4'})
%             ylabel('index (steady state=100)')
%             legend('low-wage/low-wealth','high-wage/low-wealth','low-wage/high-wealth','high-wage/high-wealth','location','best','Interpreter','latex')
%             if nprint
%                 box on
%                 legend off
%                 if jobnum==0
%                     filename = sprintf('figures/oh_policy_m_s_no_%1i.pdf',projects(i_proj));
%                 elseif jobnum<=5
%                     filename = sprintf('figures/oh_policy_m_s_base_%1i.pdf',projects(i_proj));
%                 else
%                     filename = sprintf('figures/oh_policy_m_s_%2i_%1i.pdf',jobnum,projects(i_proj));
%                 end
%                 set(gcf,'PaperUnits','inches');
%                 set(gcf,'PaperPositionMode','manual');
%                 set(gcf,'PaperSize',myPaperSize);
%                 set(gcf,'Units','inches');
%                 set(gcf,'Position',[lowerLeft myPaperSize]);
%                 export_pdf(gcf,filename);
%             else
%                 title('Outside Hours policy (middle age+susceptible)')
%             end

            %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

            figure;
            hold on
            plot([2020,tgrid(1:Tend+1)],[100; ...
                100*oc_y_lzlkpath{i_proj}(1:Tend+1) ...
                /oc_y_lzlkpath{i_proj}(1)],'b','LineWidth',2);
            plot([2020,tgrid(1:Tend+1)],[100; ...
                100*oc_y_hzlkpath{i_proj}(1:Tend+1) ...
                /oc_y_hzlkpath{i_proj}(1)],'r--','LineWidth',2);
            plot([2020,tgrid(1:Tend+1)],[100; ...
                100*oc_y_lzhkpath{i_proj}(1:Tend+1) ...
                /oc_y_lzhkpath{i_proj}(1)],'k-.','LineWidth',2);
            plot([2020,tgrid(1:Tend+1)],[100; ...
                100*oc_y_hzhkpath{i_proj}(1:Tend+1) ...
                /oc_y_hzhkpath{i_proj}(1)],'m:','LineWidth',2);
            hold off
            %xlabel('time')
            %xlim([2020.15 2022])
            ylim([0 105])
            %xticks([2020:0.25:2022])
            %xticklabels({'','2020q2','','2020q4','','2021q2','','2021q4'})
            xlim([2020.15 2021.5])
            ax = gca;
            ax.XTick = [2020:0.25:2022];
            ax.XTickLabel = '';
            myLabels = {'',' q2 ',' q3 ',' q4 ','q1',' q2 ','q3',' q4 '; ...
                '','(2020)','','','(2021)','','',''};
            for i = 1:length(myLabels)
                text(ax.XTick(i), ax.YLim(1), sprintf('%s\n%s', myLabels{:,i}), ...
                    'horizontalalignment','center','verticalalignment','top');
            end
            ax.XLabel.String = sprintf('\n','');
            ylabel('index (steady state=100)')
            legend('low-wage/low-wealth','high-wage/low-wealth','low-wage/high-wealth','high-wage/high-wealth','location','south','Interpreter','latex','FontSize',11)
            if nprint
                box on
                if jobnum==0
                    filename = sprintf('figures/oc_y_no_%1i.pdf',projects(i_proj));
                elseif jobnum<=5
                    filename = sprintf('figures/oc_y_base_%1i.pdf',projects(i_proj));
                else
                    filename = sprintf('figures/oc_y_%2i_%1i.pdf',jobnum,projects(i_proj));
                end
                set(gcf,'PaperUnits','inches');
                set(gcf,'PaperPositionMode','manual');
                set(gcf,'PaperSize',myPaperSize);
                set(gcf,'Units','inches');
                set(gcf,'Position',[lowerLeft myPaperSize]);
                export_pdf(gcf,filename);
            else
                title('Outside Consumption (young)')
            end

            figure;
            hold on
            plot([2020,tgrid(1:Tend+1)],[100; ...
                100*oc_m_lzlkpath{i_proj}(1:Tend+1) ...
                /oc_m_lzlkpath{i_proj}(1)],'b','LineWidth',2);
            plot([2020,tgrid(1:Tend+1)],[100; ...
                100*oc_m_hzlkpath{i_proj}(1:Tend+1) ...
                /oc_m_hzlkpath{i_proj}(1)],'r--','LineWidth',2);
            plot([2020,tgrid(1:Tend+1)],[100; ...
                100*oc_m_lzhkpath{i_proj}(1:Tend+1) ...
                /oc_m_lzhkpath{i_proj}(1)],'k-.','LineWidth',2);
            plot([2020,tgrid(1:Tend+1)],[100; ...
                100*oc_m_hzhkpath{i_proj}(1:Tend+1) ...
                /oc_m_hzhkpath{i_proj}(1)],'m:','LineWidth',2);
            hold off
            %xlabel('time')
            ylabel('index (steady state=100)')
            ylim([0 105])
            %xlim([2020.15 2022])
            %xticks([2020:0.25:2022])
            %xticklabels({'','2020q2','','2020q4','','2021q2','','2021q4'})
            xlim([2020.15 2021.5])
            ax = gca;
            ax.XTick = [2020:0.25:2022];
            ax.XTickLabel = '';
            myLabels = {'',' q2 ',' q3 ',' q4 ','q1',' q2 ','q3',' q4 '; ...
                '','(2020)','','','(2021)','','',''};
            for i = 1:length(myLabels)
                text(ax.XTick(i), ax.YLim(1), sprintf('%s\n%s', myLabels{:,i}), ...
                    'horizontalalignment','center','verticalalignment','top');
            end
            ax.XLabel.String = sprintf('\n','');
            legend('low-wage/low-wealth','high-wage/low-wealth','low-wage/high-wealth','high-wage/high-wealth','location','best','Interpreter','latex')
            if nprint
                box on
                legend off
                if jobnum==0
                    filename = sprintf('figures/oc_m_no_%1i.pdf',projects(i_proj));
                elseif jobnum<=5
                    filename = sprintf('figures/oc_m_base_%1i.pdf',projects(i_proj));
                else
                    filename = sprintf('figures/oc_m_%2i_%1i.pdf',jobnum,projects(i_proj));
                end
                set(gcf,'PaperUnits','inches');
                set(gcf,'PaperPositionMode','manual');
                set(gcf,'PaperSize',myPaperSize);
                set(gcf,'Units','inches');
                set(gcf,'Position',[lowerLeft myPaperSize]);
                export_pdf(gcf,filename);
            else
                title('Outside consumption (middle age)')
            end

            figure;
            hold on
            plot([2020,tgrid(1:Tend+1)],[100; ...
                100*oc_r_lkpath{i_proj}(1:Tend+1) ...
                /oc_r_lkpath{i_proj}(1)],'b','LineWidth',2);
            plot([2020,tgrid(1:Tend+1)],[100; ...
                100*oc_r_hkpath{i_proj}(1:Tend+1) ...
                /oc_r_hkpath{i_proj}(1)],'m:','LineWidth',2);
            hold off
            %xlabel('time')
            %xlim([2020.15 2022])
            ylim([0 105])
            %xticks([2020:0.25:2022])
            %xticklabels({'','2020q2','','2020q4','','2021q2','','2021q4'})
            xlim([2020.15 2021.5])
            ax = gca;
            ax.XTick = [2020:0.25:2022];
            ax.XTickLabel = '';
            myLabels = {'',' q2 ',' q3 ',' q4 ','q1',' q2 ','q3',' q4 '; ...
                '','(2020)','','','(2021)','','',''};
            for i = 1:length(myLabels)
                text(ax.XTick(i), ax.YLim(1), sprintf('%s\n%s', myLabels{:,i}), ...
                    'horizontalalignment','center','verticalalignment','top');
            end
            ax.XLabel.String = sprintf('\n','');
            ylabel('index (steady state=100)')
            legend('low-wealth','high-wealth','location','best','Interpreter','latex','FontSize',11)
            if nprint
                box on
                %legend off
                if jobnum==0
                    filename = sprintf('figures/oc_o_s_no_%1i.pdf',projects(i_proj));
                elseif jobnum<=5
                    filename = sprintf('figures/oc_o_base_%1i.pdf',projects(i_proj));
                else
                    filename = sprintf('figures/oc_o_%2i_%1i.pdf',jobnum,projects(i_proj));
                end
                set(gcf,'PaperUnits','inches');
                set(gcf,'PaperPositionMode','manual');
                set(gcf,'PaperSize',myPaperSize);
                set(gcf,'Units','inches');
                set(gcf,'Position',[lowerLeft myPaperSize]);
                export_pdf(gcf,filename);
            else
                title('Outside consumption (old)')
            end            

            figure;
            hold on
            plot([2020,tgrid(1:Tend+1)],[100; ...
                100*oh_y_lzlkpath{i_proj}(1:Tend+1) ...
                /oh_y_lzlkpath{i_proj}(1)],'b','LineWidth',2);
            plot([2020,tgrid(1:Tend+1)],[100; ...
                100*oh_y_hzlkpath{i_proj}(1:Tend+1) ...
                /oh_y_hzlkpath{i_proj}(1)],'r--','LineWidth',2);
            plot([2020,tgrid(1:Tend+1)],[100; ...
                100*oh_y_lzhkpath{i_proj}(1:Tend+1) ...
                /oh_y_lzhkpath{i_proj}(1)],'k-.','LineWidth',2);
            plot([2020,tgrid(1:Tend+1)],[100; ...
                100*oh_y_hzhkpath{i_proj}(1:Tend+1) ...
                /oh_y_hzhkpath{i_proj}(1)],'m:','LineWidth',2);
            hold off
            %xlabel('time')
            %xlim([2020.15 2022])
            ylim([0 105])
            %xticks([2020:0.25:2022])
            %xticklabels({'',['2020';' q2 '],'','2020q4','','2021q2','','2021q4'})
            xlim([2020.15 2021.5])
            ax = gca;
            ax.XTick = [2020:0.25:2022];
            ax.XTickLabel = '';
            myLabels = {'',' q2 ',' q3 ',' q4 ','q1',' q2 ','q3',' q4 '; ...
                '','(2020)','','','(2021)','','',''};
            for i = 1:length(myLabels)
                text(ax.XTick(i), ax.YLim(1), sprintf('%s\n%s', myLabels{:,i}), ...
                    'horizontalalignment','center','verticalalignment','top');
            end
            ax.XLabel.String = sprintf('\n','');
            ylabel('index (steady state=100)')
            legend('low-wage/low-wealth','high-wage/low-wealth','low-wage/high-wealth','high-wage/high-wealth','location','best','Interpreter','latex')
            if nprint
                box on
                legend off
                if jobnum==0
                    filename = sprintf('figures/oh_y_no_%1i.pdf',projects(i_proj));
                elseif jobnum<=5
                    filename = sprintf('figures/oh_y_base_%1i.pdf',projects(i_proj));
                else
                    filename = sprintf('figures/oh_y_%2i_%1i.pdf',jobnum,projects(i_proj));
                end
                set(gcf,'PaperUnits','inches');
                set(gcf,'PaperPositionMode','manual');
                set(gcf,'PaperSize',myPaperSize);
                set(gcf,'Units','inches');
                set(gcf,'Position',[lowerLeft myPaperSize]);
                export_pdf(gcf,filename);
            else
                title('Outside Hours (young)')
            end


            figure;
            hold on
            plot([2020,tgrid(1:Tend+1)],[100; ...
                100*oh_m_lzlkpath{i_proj}(1:Tend+1) ...
                /oh_m_lzlkpath{i_proj}(1)],'b','LineWidth',2);
            plot([2020,tgrid(1:Tend+1)],[100; ...
                100*oh_m_hzlkpath{i_proj}(1:Tend+1) ...
                /oh_m_hzlkpath{i_proj}(1)],'r--','LineWidth',2);
            plot([2020,tgrid(1:Tend+1)],[100; ...
                100*oh_m_lzhkpath{i_proj}(1:Tend+1) ...
                /oh_m_lzhkpath{i_proj}(1)],'k-.','LineWidth',2);
            plot([2020,tgrid(1:Tend+1)],[100; ...
                100*oh_m_hzhkpath{i_proj}(1:Tend+1) ...
                /oh_m_hzhkpath{i_proj}(1)],'m:','LineWidth',2);
            hold off
            %xlabel('time')
            %xlim([2020.15 2022])
            ylim([0 105])
            %xticks([2020:0.25:2022])
            %xticklabels({'','2020q2','','2020q4','','2021q2','','2021q4'})
            xlim([2020.15 2021.5])
            ax = gca;
            ax.XTick = [2020:0.25:2022];
            ax.XTickLabel = '';
            myLabels = {'',' q2 ',' q3 ',' q4 ','q1',' q2 ','q3',' q4 '; ...
                '','(2020)','','','(2021)','','',''};
            for i = 1:length(myLabels)
                text(ax.XTick(i), ax.YLim(1), sprintf('%s\n%s', myLabels{:,i}), ...
                    'horizontalalignment','center','verticalalignment','top');
            end
            ax.XLabel.String = sprintf('\n','');
            ylabel('index (steady state=100)')
            legend('low-wage/low-wealth','high-wage/low-wealth','low-wage/high-wealth','high-wage/high-wealth','location','best','Interpreter','latex')
            if nprint
                box on
                legend off
                if jobnum==0
                    filename = sprintf('figures/oh_m_no_%1i.pdf',projects(i_proj));
                elseif jobnum<=5
                    filename = sprintf('figures/oh_m_base_%1i.pdf',projects(i_proj));
                else
                    filename = sprintf('figures/oh_m_%2i_%1i.pdf',jobnum,projects(i_proj));
                end
                set(gcf,'PaperUnits','inches');
                set(gcf,'PaperPositionMode','manual');
                set(gcf,'PaperSize',myPaperSize);
                set(gcf,'Units','inches');
                set(gcf,'Position',[lowerLeft myPaperSize]);
                export_pdf(gcf,filename);
            else
                title('Outside Hours (middle age)')
            end

        end
        
        
        figure;
        hold on
        plot(tgrid,100*i1path{i_proj}/sum(sum(muss(1,:,:))),'b','LineWidth',2);
        plot(tgrid,100*i2path{i_proj}/sum(sum(muss(2,:,:))),'r-','LineWidth',2);
        plot(tgrid,100*i3path{i_proj}/sum(sum(muss(3,:,:))),'k-.','LineWidth',2);
        hold off
        %xlabel('time')
        %xlim([2020.15 2022])
        %xticks([2020:0.25:2022])
        %xticklabels({'','2020q2','','2020q4','','2021q2','','2021q4'})
        xlim([2020.15 2021.5])
        ax = gca;
        ax.XTick = [2020:0.25:2022];
        ax.XTickLabel = '';
        myLabels = {'',' q2 ',' q3 ',' q4 ','q1',' q2 ','q3',' q4 '; ...
            '','(2020)','','','(2021)','','',''};
        for i = 1:length(myLabels)
            text(ax.XTick(i), ax.YLim(1), sprintf('%s\n%s', myLabels{:,i}), ...
                'horizontalalignment','center','verticalalignment','top');
        end
        ax.XLabel.String = sprintf('\n','');
        ylabel('index (steady state=100)')
        ylabel('percent')
        ymax = max(0.1,ceil(1000*max(i1path{i_proj})*3)/10);
        if ymax>1
            ymax = ceil(100*max(i1path{i_proj})*3);
        end
        ylim([0 ymax])
        if ymax<1
            yticks([0:0.1:ymax])
        elseif ymax<2
            yticks([0:0.5:ymax])
        elseif ymax<8
            yticks([0:1:ymax])
        else
            yticks([0:5:ymax])
        end
        legend('young','middle','old','Location','best','Interpreter','latex','FontSize',11)
        if nprint
            box on
            if jobnum==0
                filename = sprintf('figures/tran_disease_age_no_%1i.pdf',projects(i_proj));
            elseif jobnum<=5
                filename = sprintf('figures/tran_disease_age_base_%1i.pdf',projects(i_proj));
            else
                filename = sprintf('figures/tran_disease_age_%2i_%1i.pdf',jobnum,projects(i_proj));
            end
            set(gcf,'PaperUnits','inches');
            set(gcf,'PaperPositionMode','manual');
            set(gcf,'PaperSize',myPaperSize);
            set(gcf,'Units','inches');
            set(gcf,'Position',[lowerLeft myPaperSize]);
            export_pdf(gcf,filename);
        else
            title('Disease progression')
        end        
        
    end % if runtran
    
    
end % for i_proj

if runtran

            
    fprintf('\n<<inputs for TABLES 4 and 5>>')
    fprintf('\nwelfare:       %6.1f  %6.1f  %6.1f  %6.1f  %6.1f', avgwelf{:} );
    fprintf('\n working age:  %6.1f  %6.1f  %6.1f  %6.1f  %6.1f', working_age{:} );
    fprintf('\n retired:      %6.1f  %6.1f  %6.1f  %6.1f  %6.1f', retired{:} );
    fprintf('\n low wage:     %6.1f  %6.1f  %6.1f  %6.1f  %6.1f', low_z{:} );
    fprintf('\n high wage:    %6.1f  %6.1f  %6.1f  %6.1f  %6.1f', high_z{:} );
    fprintf('\n low wealth:   %6.1f  %6.1f  %6.1f  %6.1f  %6.1f', low_k{:} );
    fprintf('\n high wealth:  %6.1f  %6.1f  %6.1f  %6.1f  %6.1f', high_k{:} );
    fprintf('\n-----------------------------------------')
    if projects(1)==2 && Np>1
    fprintf('\npolicy support:%6.1f  %6.1f  %6.1f  %6.1f  %6.1f', poswelf{:} );
    end

    fprintf('\n2-year output: %6.1f  %6.1f  %6.1f  %6.1f  %6.1f', yindex{:} );
    fprintf('\ndeaths/10k:    %6.1f  %6.1f  %6.1f  %6.1f  %6.1f\n', death_per_10k{:} );

    if projects(i_proj) == 3
    xx = [(0:1:53)' ypath{i_proj}(1:54) cpath{i_proj}(1:54) ...
        dpath{i_proj}(1:54) ocpath{i_proj}(1:54) ...
        hourstotpath{i_proj}(1:54) ohpath{i_proj}(1:54) ...
        i1path{i_proj}(1:54)+i2path{i_proj}(1:54)+i3path{i_proj}(1:54) ...
        takeuppath{i_proj}(1:54) betaepath{i_proj}(1:54) ...
		d1path{i_proj}(1:54) d2path{i_proj}(1:54) d3path{i_proj}(1:54) ...
        oh_y_lzlkpath{i_proj}(1:54) oh_y_hzlkpath{i_proj}(1:54) ...
        oh_y_lzhkpath{i_proj}(1:54) oh_y_hzhkpath{i_proj}(1:54) ...
        oh_m_lzlkpath{i_proj}(1:54) oh_m_hzlkpath{i_proj}(1:54) ...
        oh_m_lzhkpath{i_proj}(1:54) oh_m_hzhkpath{i_proj}(1:54) ...
        oc_y_lzlkpath{i_proj}(1:54) oc_y_hzlkpath{i_proj}(1:54) ...
        oc_y_lzhkpath{i_proj}(1:54) oc_y_hzhkpath{i_proj}(1:54) ...
        oc_m_lzlkpath{i_proj}(1:54) oc_m_hzlkpath{i_proj}(1:54) ...
        oc_m_lzhkpath{i_proj}(1:54) oc_m_hzhkpath{i_proj}(1:54) ...
        oc_r_lkpath{i_proj}(1:54) oc_r_hkpath{i_proj}(1:54) ];
    filename = sprintf('results/model%i.xls',jobnum);
    writematrix(xx,filename,'Range','A1');
    end
end

if length(projects)>1
    x = 100*[max(cumdpath{1})];
    y = 100*[yloss{1}];
    for i = 2:Np
        x = [x;100*max(cumdpath{i})];
        y = [y;100*yloss{i}];
    end
    if Np==2
        g = {projectname{1};projectname{2}};
    elseif Np==3
        g = {projectname{1};projectname{2};projectname{3}};
    elseif Np==4
        g = {projectname{1};projectname{2};projectname{3};projectname{4}};
    elseif Np==5
        g = {projectname{1};projectname{2};projectname{3};projectname{4};projectname{5}};
    elseif Np==6
        g = {projectname{1};projectname{2};projectname{3};projectname{4};projectname{5};projectname{6}};
    elseif Np==7
        g = {projectname{1};projectname{2};projectname{3};projectname{4};projectname{5};projectname{6};projectname{7}};
    end
    

    
    
    figure;
    hold on
    for i=1:Np
        if i==1
            plot(tgrid,repnum{i},'k','LineWidth',2);
        elseif i==2
            plot(tgrid,repnum{i},'r--','LineWidth',2);
        elseif i==3
            plot(tgrid,repnum{i},'b-.','LineWidth',2);
        elseif i==4
            plot(tgrid,repnum{i},'m:','LineWidth',2);
        else
            plot(tgrid,repnum{i},'LineWidth',2);
        end
    end
    hold off
    %xlabel('time')
    xlim([2020.15 2021.5])
    %xticks([2020:0.25:2022])
    ylim([0.5,2.5])
    yticks([0:1:3])
    %xticklabels({'','2020q2','','2020q4','','2021q2','','2021q4'})
    ax = gca;
    ax.XTick = [2020:0.25:2022];
    ax.XTickLabel = '';
    myLabels = {'',' q2 ',' q3 ',' q4 ','q1',' q2 ','q3',' q4 '; ...
        '','(2020)','','','(2021)','','',''};
    for i = 1:length(myLabels)
        text(ax.XTick(i), ax.YLim(1), sprintf('%s\n%s', myLabels{:,i}), ...
            'horizontalalignment','center','verticalalignment','top');
    end
    ax.XLabel.String = sprintf('\n','');
    legend(g,'location','northeast','Interpreter','latex','FontSize',11)
    
    if nprint
        box on
        if jobnum==0
            filename = sprintf('figures/tran_repnum_no_%1i%1i%1i%1i%1i%1i%1i%1i.pdf',projects(1:Np));
        elseif jobnum<=5
            filename = sprintf('figures/tran_repnum_base_%1i%1i%1i%1i%1i%1i%1i%1i.pdf',projects(1:Np));
        else
            filename = sprintf('figures/tran_repnum_%2i_%1i%1i%1i%1i%1i%1i%1i%1i.pdf',jobnum,projects(1:Np));
        end
        set(gcf,'PaperUnits','inches');
        set(gcf,'PaperPositionMode','manual');
        set(gcf,'PaperSize',myPaperSize);
        set(gcf,'Units','inches');
        set(gcf,'Position',[lowerLeft myPaperSize]);
        export_pdf(gcf,filename);
    else
        title('Reproduction number')
    end
    
    figure;
    hold on
    for i=1:Np
        if i==1
            plot(tgrid,100*bigipath{i},'k','LineWidth',2);
        elseif i==2
            plot(tgrid,100*bigipath{i},'r--','LineWidth',2);
        elseif i==3
            plot(tgrid,100*bigipath{i},'b-.','LineWidth',2);
        elseif i==4
            plot(tgrid,100*bigipath{i},'m:','LineWidth',2);
        else
            plot(tgrid,100*bigipath{i},'LineWidth',2);
        end
    end
    hold off
    %xlabel('time')
    %xlim([2020.15 2022])
    %xticks([2020:0.25:2022])
    %xticklabels({'','2020q2','','2020q4','','2021q2','','2021q4'})
    ylabel('percent')
    ymax = ceil(max(100*bigipath{1}));
    for i=1:Np
        ymax = max(ymax,ceil(max(100*bigipath{i})));
    end
    ylim([0,ymax*1.05])
    if ymax<10
        yticks([0:1:ymax])
    else
        yticks([0:5:ymax])
    end
    xlim([2020.15 2021.5])
    ax = gca;
    ax.XTick = [2020:0.25:2022];
    ax.XTickLabel = '';
    myLabels = {'',' q2 ',' q3 ',' q4 ','q1',' q2 ','q3',' q4 '; ...
        '','(2020)','','','(2021)','','',''};
    for i = 1:length(myLabels)
        text(ax.XTick(i), ax.YLim(1), sprintf('%s\n%s', myLabels{:,i}), ...
            'horizontalalignment','center','verticalalignment','top');
    end
    ax.XLabel.String = sprintf('\n','');
    legend(g,'location','northwest','Interpreter','latex','FontSize',11)
    if nprint
        box on
        %legend off
        if jobnum==0
            filename = sprintf('figures/tran_disease_no_%1i%1i%1i%1i%1i%1i%1i%1i.pdf',projects(1:Np));
        elseif jobnum<=5
            filename = sprintf('figures/tran_disease_base_%1i%1i%1i%1i%1i%1i%1i%1i.pdf',projects(1:Np));
        else
            filename = sprintf('figures/tran_disease_%2i_%1i%1i%1i%1i%1i%1i%1i%1i.pdf',jobnum,projects(1:Np));
        end
        set(gcf,'PaperUnits','inches');
        set(gcf,'PaperPositionMode','manual');
        set(gcf,'PaperSize',myPaperSize);
        set(gcf,'Units','inches');
        set(gcf,'Position',[lowerLeft myPaperSize]);
        export_pdf(gcf,filename);
    else
        title('Infected rate')
    end
    
    figure;
    hold on
    for i=1:Np
        if i==1
            plot(tgrid,100*cumdpath{i},'k','LineWidth',2);
        elseif i==2
            plot(tgrid,100*cumdpath{i},'r--','LineWidth',2);
        elseif i==3
            plot(tgrid,100*cumdpath{i},'b-.','LineWidth',2);
        elseif i==4
            plot(tgrid,100*cumdpath{i},'m:','LineWidth',2);
        else
            plot(tgrid,100*cumdpath{i},'LineWidth',2);
        end
    end
    hold off
    %xlabel('time')
    ylabel('percent')
    xlim([2020.15 2021.25])
    %xticks([2020:0.25:2022])
    %xticklabels({'','2020q2','','2020q4','','2021q2','','2021q4'})
    ylabel('percent')
    ymax = ceil(max(100*cumdpath{1}));
    for i=1:Np
        ymax = max(ymax,ceil(max(100*cumdpath{i})));
    end
    ylim([0,ymax])
    if ymax<2
        yticks([0:0.5:ymax])
    elseif ymax<10
        yticks([0:1:ymax])
    else
        yticks([0:5:ymax])
    end
    xlim([2020.15 2021.5])
    ax = gca;
    ax.XTick = [2020:0.25:2022];
    ax.XTickLabel = '';
    myLabels = {'',' q2 ',' q3 ',' q4 ','q1',' q2 ','q3',' q4 '; ...
        '','(2020)','','','(2021)','','',''};
    for i = 1:length(myLabels)
        text(ax.XTick(i), ax.YLim(1), sprintf('%s\n%s', myLabels{:,i}), ...
            'horizontalalignment','center','verticalalignment','top');
    end
    ax.XLabel.String = sprintf('\n','');
    legend(g,'location','best','Interpreter','latex')
    if nprint
        box on
        legend off
        if jobnum==0
            filename = sprintf('figures/tran_deaths_no_%1i%1i%1i%1i%1i%1i%1i%1i.pdf',projects(1:Np));
        elseif jobnum<=5
            filename = sprintf('figures/tran_deaths_base_%1i%1i%1i%1i%1i%1i%1i%1i.pdf',projects(1:Np));
        else
            filename = sprintf('figures/tran_deaths_%2i_%1i%1i%1i%1i%1i%1i%1i%1i.pdf',jobnum,projects(1:Np));
        end
        set(gcf,'PaperUnits','inches');
        set(gcf,'PaperPositionMode','manual');
        set(gcf,'PaperSize',myPaperSize);
        set(gcf,'Units','inches');
        set(gcf,'Position',[lowerLeft myPaperSize]);
        export_pdf(gcf,filename);
    else
        title('Death rate (cumulative)')
    end    

    
    figure;
    hold on
    for i=1:length(projects)
        if i==1
            plot([2020,tgrid],[100;100*ohpath{i}/hourstotpath{i}(1)],'k','LineWidth',2);
        elseif i==2
            plot([2020,tgrid],[100;100*ohpath{i}/hourstotpath{i}(1)],'r--','LineWidth',2);
        elseif i==3
            plot([2020,tgrid],[100;100*ohpath{i}/hourstotpath{i}(1)],'b-.','LineWidth',2);
        elseif i==4
            plot([2020,tgrid],[100;100*ohpath{i}/hourstotpath{i}(1)],'m:','LineWidth',2);
        else
            plot([2020,tgrid],[100;100*ohpath{i}/hourstotpath{i}(1)],'LineWidth',2);
        end
    end
    hold off
    xlabel('time')
    ylabel('index (steady state=100)')
    %xticks([2020:0.25:2022])
    %xticklabels({'','2020q2','','2020q4','','2021q2','','2021q4'})
    ylim([0,105])
    yticks([0:20:120])
    xlim([2020.15 2021.5])
    ax = gca;
    ax.XTick = [2020:0.25:2022];
    ax.XTickLabel = '';
    myLabels = {'',' q2 ',' q3 ',' q4 ','q1',' q2 ','q3',' q4 '; ...
        '','(2020)','','','(2021)','','',''};
    for i = 1:length(myLabels)
        text(ax.XTick(i), ax.YLim(1), sprintf('%s\n%s', myLabels{:,i}), ...
            'horizontalalignment','center','verticalalignment','top');
    end
    ax.XLabel.String = sprintf('\n','');
    legend(g,'location','best','Interpreter','latex')
    if nprint
        box on
        legend off
        if jobnum==0
            filename = sprintf('figures/tran_oh_no_%1i%1i%1i%1i%1i%1i%1i%1i.pdf',projects(1:Np));
        elseif jobnum<=5
            filename = sprintf('figures/tran_oh_base_%1i%1i%1i%1i%1i%1i%1i%1i.pdf',projects(1:Np));
        else
            filename = sprintf('figures/tran_oh_%2i_%1i%1i%1i%1i%1i%1i%1i%1i.pdf',jobnum,projects(1:Np));
        end
        set(gcf,'PaperUnits','inches');
        set(gcf,'PaperPositionMode','manual');
        set(gcf,'PaperSize',myPaperSize);
        set(gcf,'Units','inches');
        set(gcf,'Position',[lowerLeft myPaperSize]);
        export_pdf(gcf,filename);
    else
        title('Aggregate outside hours')
    end    
    
    figure;
    hold on
    for i=1:length(projects)
        if i==1
            plot([2020,tgrid],[100;100*ocpath{i}/ocpath{i}(1)],'k','LineWidth',2);
        elseif i==2
            plot([2020,tgrid],[100;100*ocpath{i}/ocpath{i}(1)],'r--','LineWidth',2);
        elseif i==3
            plot([2020,tgrid],[100;100*ocpath{i}/ocpath{i}(1)],'b-.','LineWidth',2);
        elseif i==4
            plot([2020,tgrid],[100;100*ocpath{i}/ocpath{i}(1)],'m:','LineWidth',2);
        else
            plot([2020,tgrid],[100;100*ocpath{i}/ocpath{i}(1)],'LineWidth',2);
        end
    end
    hold off
    %xlabel('time')
    ylabel('index (steady state=100)')
    %xlim([2020.15 2022])
    %xticks([2020:0.25:2022])
    %xticklabels({'','2020q2','','2020q4','','2021q2','','2021q4'})
    ylim([0,105])
    yticks([0:20:120])
    xlim([2020.15 2021.5])
    ax = gca;
    ax.XTick = [2020:0.25:2022];
    ax.XTickLabel = '';
    myLabels = {'',' q2 ',' q3 ',' q4 ','q1',' q2 ','q3',' q4 '; ...
        '','(2020)','','','(2021)','','',''};
    for i = 1:length(myLabels)
        text(ax.XTick(i), ax.YLim(1), sprintf('%s\n%s', myLabels{:,i}), ...
            'horizontalalignment','center','verticalalignment','top');
    end
    ax.XLabel.String = sprintf('\n','');
    legend(g,'location','best','Interpreter','latex')
    if nprint
        box on
        legend off
        if jobnum==0
            filename = sprintf('figures/tran_oc_no_%1i%1i%1i%1i%1i%1i%1i%1i.pdf',projects(1:Np));
        elseif jobnum<=5
            filename = sprintf('figures/tran_oc_base_%1i%1i%1i%1i%1i%1i%1i%1i.pdf',projects(1:Np));
        else
            filename = sprintf('figures/tran_oc_%2i_%1i%1i%1i%1i%1i%1i%1i%1i.pdf',jobnum,projects(1:Np));
        end
        set(gcf,'PaperUnits','inches');
        set(gcf,'PaperPositionMode','manual');
        set(gcf,'PaperSize',myPaperSize);
        set(gcf,'Units','inches');
        set(gcf,'Position',[lowerLeft myPaperSize]);
        export_pdf(gcf,filename);
    else
        title('Aggregate outside consumption')
    end
    
    
    figure;
    hold on
    for i=1:length(projects)
        if i==1
            plot([2020,tgrid],[100;100*ypath{i}/aggy],'k','LineWidth',2);
        elseif i==2
            plot([2020,tgrid],[100;100*ypath{i}/aggy],'r--','LineWidth',2);
        elseif i==3
            plot([2020,tgrid],[100;100*ypath{i}/aggy],'b-.','LineWidth',2);
        elseif i==4
            plot([2020,tgrid],[100;100*ypath{i}/aggy],'m:','LineWidth',2);
        else
            plot([2020,tgrid],[100;100*ypath{i}/aggy],'LineWidth',2);
        end
    end
    hold off
    %xlabel('time')
    ylabel('index (steady state=100)')
    %xlim([2020.15 2022])
    %xticks([2020:0.25:2022])
    %xticklabels({'','2020q2','','2020q4','','2021q2','','2021q4'})
    ylim([0,105])
    yticks([0:20:120])
    xlim([2020.15 2021.5])
    ax = gca;
    ax.XTick = [2020:0.25:2022];
    ax.XTickLabel = '';
    myLabels = {'',' q2 ',' q3 ',' q4 ','q1',' q2 ','q3',' q4 '; ...
        '','(2020)','','','(2021)','','',''};
    for i = 1:length(myLabels)
        text(ax.XTick(i), ax.YLim(1), sprintf('%s\n%s', myLabels{:,i}), ...
            'horizontalalignment','center','verticalalignment','top');
    end
    ax.XLabel.String = sprintf('\n','');
    legend(g,'location','best','Interpreter','latex')
    if nprint
        box on
        legend off
        if jobnum==0
            filename = sprintf('figures/tran_y_no_%1i%1i%1i%1i%1i%1i%1i%1i.pdf',projects(1:Np));
        elseif jobnum<=5
            filename = sprintf('figures/tran_y_base_%1i%1i%1i%1i%1i%1i%1i%1i.pdf',projects(1:Np));
        else
            filename = sprintf('figures/tran_y_%2i_%1i%1i%1i%1i%1i%1i%1i%1i.pdf',jobnum,projects(1:Np));
        end
        set(gcf,'PaperUnits','inches');
        set(gcf,'PaperPositionMode','manual');
        set(gcf,'PaperSize',myPaperSize);
        set(gcf,'Units','inches');
        set(gcf,'Position',[lowerLeft myPaperSize]);
        export_pdf(gcf,filename);
    else
        title('Aggregate output')
    end

end % if projecttype==0










% Helper: PDF export matching the original hgexport behavior
% (fixed 12pt fonts; appends .pdf when the name lacks an extension)
function export_pdf(fig, filename)
    if ~endsWith(filename,'.pdf'); filename = [filename '.pdf']; end
    set(findall(fig,'-property','FontSize'),'FontSize',12);
    exportgraphics(fig, filename, 'ContentType','vector');
end
