% %% Config

switch VAR_SET
    case 13
        % replace average size with average size based on job creation if
        % missing 
        vlog_age0_sz = log(vZit_age0_emp_msa ./ vZit_age0_firms_msa);
        vlog_age0_sz_alt = log(vZit_age0_job_ctn_births_msa ./ vZit_age0_firms_msa);
        vlog_age0_sz(isnan(vlog_age0_sz))=vlog_age0_sz_alt(isnan(vlog_age0_sz));
        
%         vlog_hpi=log(v_home_price_index);
%         dlog_hpi=[NaN(1,N); diff(vlog_hpi)];
        
        vlog_emp = log(vemp);
        vlog_pop = log(vpop);
        vfirm_entry_rate_lvl=100*[NaN(1,N); 2*vZit_age0_firms_msa(2:end,:)./(vZit_ageall_firms_msa(2:end,:)+vZit_ageall_firms_msa(1:end-1,:))];
        vfirm_entry_rate=100*[NaN(1,N); 2*diff(vZit_age0_firms_msa(1:end,:))./(vZit_ageall_firms_msa(2:end,:)+vZit_ageall_firms_msa(1:end-1,:))];
        vZit_age1_firmdeath_firms_msa_org=vZit_age1_firmdeath_firms_msa;
        for nn=1:size(vZit_age1_firmdeath_firms_msa,2)
            if any(isnan(vZit_age1_firmdeath_firms_msa(2:end,nn)))
                for tt=find(isnan(vZit_age1_firmdeath_firms_msa(2:end,nn)))
                    vZit_age1_firmdeath_firms_msa(1+tt,nn)=vZit_age1_firmdeath_firms_msa(1+tt-1,nn);
                end
            end
        end
        vfirm_exit_rate_lvl=100*[NaN(1,N); 2*vZit_age1_firmdeath_firms_msa(2:end,:)./(vZit_ageall_firms_msa(2:end,:)+vZit_ageall_firms_msa(1:end-1,:))]; 
        vfirm_exit_rate=100*[NaN(1,N); 2*diff(vZit_age1_firmdeath_firms_msa(1:end,:))./(vZit_ageall_firms_msa(2:end,:)+vZit_ageall_firms_msa(1:end-1,:))]; 
        vfirm_exit_rate_all_lvl=100*[NaN(1,N); 2*vZit_ageall_firmdeath_firms_msa(2:end,:)./(vZit_ageall_firms_msa(2:end,:)+vZit_ageall_firms_msa(1:end-1,:))]; 
        vfirm_exit_rate_all = 100*[NaN(1,N); 2*diff(vZit_ageall_firmdeath_firms_msa(1:end,:))./(vZit_ageall_firms_msa(2:end,:)+vZit_ageall_firms_msa(1:end-1,:))]; 
        vlog_emp_pop=vlog_emp-vlog_pop;
        dlog_firms=[NaN(1,N); diff(log(vZit_ageall_firms_msa))];
        dlog_pop=[NaN(1,N); diff(vlog_pop)];
        
       
        vjob_creation_rate_births_lvl=vjob_creation_rate_births;
        vjob_creation_rate_births=vjob_creation_growth_births;
        VarNames ={'vjob_creation_rate_births'; 'vlog_emp_pop';            'dlog_pop'; 'dlog_wage' };
        VarLabels={'Gross job creation (births)'        ; 'Employment/Pop (log, BDS)';       'Pop growth';'Wage growth'};
        VarLabelsShort={'Job creation (births)'        ; 'Employment/Pop';  'Pop growth'; 'Wage growth'};
        VarScale =[1 100 100 100];

        [~,jj]=find(isnan(v_migrant_rate_exm(vYear==1983,:)));
        v_migrant_rate_exm(vYear==1983,jj)=nanmean(v_migrant_rate_exm(vYear==1983,:));
        v_migrant_rate_exm=[NaN(1,N); v_migrant_rate_exm(1:end-1,:)];
        Mig_Pop_Corr=NaN(N,2);
        for mm=1:N
            Mig_Pop_Corr(mm,1)=corr(v_migrant_rate_exm(8:end-2,mm),dlog_pop(8:end-2,mm)*100);
            Mig_Pop_Corr(mm,2)=[1 0]*cov(v_migrant_rate_exm(8:end-2,mm),dlog_pop(8:end-2,mm)*100)*[0; 1]/var(dlog_pop(8:end-2,mm)*100);
        end
        disp('migration - population growth correlation (IQR + median)');
        display(prctile(Mig_Pop_Corr(:,1),[25 50 75]));
        disp('migration - population growth slope (IQR + median)');
        display(prctile(Mig_Pop_Corr(:,2),[25 50 75]));
        PerNames ={'vfirm_entry_rate'; ...
%                     'dlog_hpi';   ...
                    'vfirm_exit_rate';'v_migrant_rate_exm';'vfirm_exit_rate_all'; 'vlog_age0_sz'};
        PerLabels={'Firm entry rate' ; ...
%                     'House price growth'; ...
                    'Firm exit rate';'Net migration rate';'Firm exit rate (all)'; 'Startup size'};
        PerLabelsShort={'Firm entry rate' ; ...
%                     'House price growth'; ...
                    'Firm exit rate';'Net migration rate';'Firm exit rate (all)'; 'Startup size'};
        PerScale =[1 1 1 1 100];
%         PerScale =[1 100 1 1 1 100];

        Suffix=[Suffix '_sm_wage_mig'];    
 case 130
        % replace average size with average size based on job creation if
        % missing 
        vlog_age0_sz = log(vZit_age0_emp_msa ./ vZit_age0_firms_msa);
        vlog_age0_sz_alt = log(vZit_age0_job_ctn_births_msa ./ vZit_age0_firms_msa);
        vlog_age0_sz(isnan(vlog_age0_sz))=vlog_age0_sz_alt(isnan(vlog_age0_sz));
        
%         vlog_hpi=log(v_home_price_index);
%         dlog_hpi=[NaN(1,N); diff(vlog_hpi)];
        
        vlog_emp = log(vemp);
        vlog_pop = log(vpop);
        fprintf(2,'Notice the different definition of the entry rate here! \n');
        vfirm_entry_rate_lvl=100*[NaN(1,N); 2*vZit_age0_firms_msa(2:end,:)./(vZit_ageall_firms_msa(2:end,:)+vZit_ageall_firms_msa(1:end-1,:))];
%         vjob_creation_growth_births = 100 * [NaN(1,size(vZit_age0_job_ctn_msa,2)); diff(vZit_age0_job_ctn_msa,1)] ./ vZit_ageall_denom_msa;
        vfirm_entry_rate=100*[NaN(1,N); 2*diff(vZit_age0_firms_msa(1:end,:))./(vZit_ageall_firms_msa(2:end,:)+vZit_ageall_firms_msa(1:end-1,:))];
        vZit_age1_firmdeath_firms_msa_org=vZit_age1_firmdeath_firms_msa;
        for nn=1:size(vZit_age1_firmdeath_firms_msa,2)
            if any(isnan(vZit_age1_firmdeath_firms_msa(2:end,nn)))
                for tt=find(isnan(vZit_age1_firmdeath_firms_msa(2:end,nn)))
                    vZit_age1_firmdeath_firms_msa(1+tt,nn)=vZit_age1_firmdeath_firms_msa(1+tt-1,nn);
                end
            end
        end
        fprintf(2,'Notice the different definition of the exit rate here! \n');
        vfirm_exit_rate_lvl=100*[NaN(1,N); 2*vZit_age1_firmdeath_firms_msa(2:end,:)./(vZit_ageall_firms_msa(2:end,:)+vZit_ageall_firms_msa(1:end-1,:))]; 
        vfirm_exit_rate=100*[NaN(1,N); 2*diff(vZit_age1_firmdeath_firms_msa(1:end,:))./(vZit_ageall_firms_msa(2:end,:)+vZit_ageall_firms_msa(1:end-1,:))]; 
        vfirm_exit_rate_all_lvl=100*[NaN(1,N); 2*vZit_ageall_firmdeath_firms_msa(2:end,:)./(vZit_ageall_firms_msa(2:end,:)+vZit_ageall_firms_msa(1:end-1,:))]; 
        vfirm_exit_rate_all = 100*[NaN(1,N); 2*diff(vZit_ageall_firmdeath_firms_msa(1:end,:))./(vZit_ageall_firms_msa(2:end,:)+vZit_ageall_firms_msa(1:end-1,:))]; 
        vlog_emp_pop=vlog_emp-vlog_pop;
        dlog_firms=[NaN(1,N); diff(log(vZit_ageall_firms_msa))];
        dlog_pop=[NaN(1,N); diff(vlog_pop)];
% %         
       
        fprintf(2,'Notice the different definition of the job creation rate here! \n');
        vjob_creation_rate_births_lvl=vjob_creation_rate_births;
        vjob_creation_rate_births=vjob_creation_growth_births;
        VarNames ={'vjob_creation_rate_births'; 'vlog_emp_pop';            'dlog_pop'; 'dlog_wage' ;'vfirm_entry_rate';'vfirm_exit_rate_all'};
        VarLabels={'Gross job creation (births)'        ; 'Employment/Pop (log, BDS)';       'Pop growth';'Wage growth'; 'Firm entry rate' ;'Firm exit rate (all)'};
        VarLabelsShort={'Job creation (births)'        ; 'Employment/Pop';  'Pop growth'; 'Wage growth'; 'Firm entry rate' ;'Firm exit rate (all)'};
        VarScale =[1 100 100 100 1 1];

        [~,jj]=find(isnan(v_migrant_rate_exm(vYear==1983,:)));
        v_migrant_rate_exm(vYear==1983,jj)=nanmean(v_migrant_rate_exm(vYear==1983,:));
        v_migrant_rate_exm=[NaN(1,N); v_migrant_rate_exm(1:end-1,:)];
        Mig_Pop_Corr=NaN(N,2);
        for mm=1:N
            Mig_Pop_Corr(mm,1)=corr(v_migrant_rate_exm(8:end-2,mm),dlog_pop(8:end-2,mm)*100);
            Mig_Pop_Corr(mm,2)=[1 0]*cov(v_migrant_rate_exm(8:end-2,mm),dlog_pop(8:end-2,mm)*100)*[0; 1]/var(dlog_pop(8:end-2,mm)*100);
        end
        disp('migration - population growth correlation (IQR + median)');
        display(prctile(Mig_Pop_Corr(:,1),[25 50 75]));
        disp('migration - population growth slope (IQR + median)');
        display(prctile(Mig_Pop_Corr(:,2),[25 50 75]));
        PerNames ={ ...
%                     'dlog_hpi';  ...
                    'vfirm_exit_rate';'v_migrant_rate_exm';'vlog_age0_sz'};
        PerLabels={ ...
%                     'House price growth'; ...
                    'Firm exit rate';'Net migration rate';'Startup size'};
        PerLabelsShort={ ...
%                     'House price growth'; ...
                    'Firm exit rate';'Net migration rate';'Startup size'};
        PerScale =[ ...
%                     100 ...
                    1 1 100];

        Suffix=[Suffix '_lg_wage_mig'];    
    otherwise
        error('Choice of variables not set up')
end


if ~exist('VarLabelsShort', 'var')
    VarLabelsShort=VarLabels;
end
p=length(VarNames);
pp=length(PerNames);

TargetMSA_ID=37980;
TargetMSA_Name=cMSA{vMSA==TargetMSA_ID}(1:min(strfind(cMSA{vMSA==TargetMSA_ID},'-'),strfind(cMSA{vMSA==TargetMSA_ID},','))-1);
figure;
for p_idx=1:p+pp
    subplot(ceil((p+pp)/2),2,p_idx);
    if p_idx<=p
        eval(['plot(vYear(vYear<=EndYear), (' VarNames{p_idx} '(vYear<=EndYear, vMSA==TargetMSA_ID)-mean(' VarNames{p_idx}  '(vYear<=EndYear,:),2))*VarScale(p_idx));']);
        title([VarLabels{p_idx} ' in ' TargetMSA_Name ' vs the nation'])
    else
        eval(['plot(vYear(vYear<=EndYear), (' PerNames{p_idx-p} '(vYear<=EndYear, vMSA==TargetMSA_ID)-mean(' PerNames{p_idx-p}  '(vYear<=EndYear,:),2))*PerScale(p_idx-p));']);
        title([PerLabels{p_idx-p} ' in ' TargetMSA_Name ' vs the nation'])
    end
    hold all;
    grid on;
end
 
TargetMSA_ID=37980;
TargetMSA_Name=cMSA{vMSA==TargetMSA_ID}(1:min(strfind(cMSA{vMSA==TargetMSA_ID},'-'),strfind(cMSA{vMSA==TargetMSA_ID},','))-1);
figure;
for p_idx=1:p+pp
    subplot(ceil((p+pp)/2),2,p_idx);
    if exist('old_vjob_creation_rate_births', 'var') && p_idx==1
        display('old_vjob_creation_rate_births');
        display([nanmean(nanmean(old_vjob_creation_rate_births(vYear<=EndYear, :),2)), nanmean(nanstd(old_vjob_creation_rate_births(vYear<=EndYear, :) ,[],2)), nanstd(nanmean(old_vjob_creation_rate_births(vYear<=EndYear, :),2))]);
    end
    if p_idx<=p
        display(VarNames{p_idx});
        display([eval(['nanmean(nanmean(' VarNames{p_idx} '(vYear<=EndYear, :)*VarScale(p_idx),2));']), eval(['nanmean(nanstd(' VarNames{p_idx} '(vYear<=EndYear, :)*VarScale(p_idx),[],2));']), eval(['nanstd(nanmean(' VarNames{p_idx} '(vYear<=EndYear, :)*VarScale(p_idx),2));'])]);
        eval(['plot(vYear(vYear<=EndYear), (' VarNames{p_idx} '(vYear<=EndYear, vMSA==TargetMSA_ID)-mean(' VarNames{p_idx}  '(vYear<=EndYear,:),2))*VarScale(p_idx));']);
        hold all;
%         eval(['hndl=plot(vYear(vYear<=EndYear), prctile(' VarNames{p_idx} '(vYear<=EndYear, :)-mean(' VarNames{p_idx}  '(:,:),2), [5 25 75 95],2)''*VarScale(p_idx));']);
        title([VarLabelsShort{p_idx} ' in ' TargetMSA_Name ' rel. to U.S.'])
%         set(hndl, {'Color'}, {[0.5 0 0];[0.5 0 0];[0.5 0 0];[0.5 0 0]},{'LineStyle'}, {':';'--';'--';':'});
    else
        display(PerNames{p_idx-p});
        display([eval(['nanmean(nanmean(' PerNames{p_idx-p} '(vYear<=EndYear, :)*PerScale(p_idx-p),2));']), eval(['nanmean(nanstd(' PerNames{p_idx-p} '(vYear<=EndYear, :)*PerScale(p_idx-p),[],2));']), eval(['nanstd(nanmean(' PerNames{p_idx-p} '(vYear<=EndYear, :)*PerScale(p_idx-p),2));'])]);
        eval(['plot(vYear(vYear<=EndYear), (' PerNames{p_idx-p} '(vYear<=EndYear, vMSA==TargetMSA_ID)-mean(' PerNames{p_idx-p}  '(vYear<=EndYear,:),2))*PerScale(p_idx-p));']);
        hold all;
%         eval(['hndl=plot(vYear(vYear<=EndYear), prctile(' PerNames{p_idx-p} '(vYear<=EndYear, :)-mean(' PerNames{p_idx-p}  '(:,:),2), [5 25 75 95],2)''*PerScale(p_idx-p));']);
        title([PerLabelsShort{p_idx-p} ' in ' TargetMSA_Name ' rel. to U.S.'])
%         set(hndl, {'Color'}, {[0.5 0 0];[0.5 0 0];[0.5 0 0];[0.5 0 0]},{'LineStyle'}, {':';'--';'--';':'});
    end
    grid on;   
end
 
%% Regional variation plots
MyFontSize=18;

figure;
for p_idx=1:p+pp
    subplot(ceil((p+pp)/2),2,p_idx);
    % Core VAR
    if p_idx<=p
%         eval(['plot(vYear(vYear<=EndYear), (' VarNames{p_idx} '(vYear<=EndYear, vMSA==TargetMSA_ID)-mean(' VarNames{p_idx}  '(vYear<=EndYear,:),2))*VarScale(p_idx));']);
%         hold all;
        CurrentData = eval([VarNames{p_idx} '(vYear<=EndYear, :);']);
        YearFE = nanmean(CurrentData,2);
        MSAFE = nanmean(CurrentData-YearFE*ones(1,size(CurrentData,2)),1);
        PredictedData_MSA = ones(size(CurrentData,1),1)*MSAFE ;
        PredictedData_Year = YearFE*ones(1,size(CurrentData,2));
        ResidualData = CurrentData - PredictedData_MSA - PredictedData_Year;
        DescriptiveCorr = [corr(ResidualData(~isnan(ResidualData(:))), CurrentData(~isnan(ResidualData(:)))) , ...
                            corr(PredictedData_MSA(~isnan(ResidualData(:))), CurrentData(~isnan(ResidualData(:)))), corr(PredictedData_Year(~isnan(ResidualData(:))), CurrentData(~isnan(ResidualData(:))))];

        
        eval(['hndl=plot(vYear(vYear<=EndYear), prctile(' VarNames{p_idx} '(vYear<=EndYear, :), [5 25 50 75 95],2)''*VarScale(p_idx), ''LineWidth'', 2);']);
        title(['Distribution of ' VarLabelsShort{p_idx} ])
        set(hndl, {'Color'}, {[0.5 0 0];[0.5 0 0];[0 0 0]; [0.5 0 0];[0.5 0 0]},{'LineStyle'}, {':';'--';'-';'--';':'});
        figure;
        eval(['hndl=plot(vYear(vYear<=EndYear), prctile(' VarNames{p_idx} '(vYear<=EndYear, :), [5 25 50 75 95],2)''*VarScale(p_idx), ''LineWidth'', 2);']);
%         title({['Distribution of ' VarLabelsShort{p_idx} '.']; [ 'MSA FE R-squared=' num2str( DescriptiveCorr(2)^2, '%1.2f') '. Year FE R-squared=' num2str( DescriptiveCorr(3)^2, '%1.2f') '.']})
        title({['MSA FE R-squared=' num2str( DescriptiveCorr(2)^2, '%1.2f') '.'];['Year FE R-squared=' num2str( DescriptiveCorr(3)^2, '%1.2f') '.']})
        set(hndl, {'Color'}, {[0.5 0 0];[0.5 0 0];[0 0 0]; [0.5 0 0];[0.5 0 0]},{'LineStyle'}, {':';'--';'-';'--';':'});
        axis([min(vYear) EndYear ylim]);
        old_ylim=ylim;
        ylabel('%');
        grid on;
        set(findall(gcf,'type','axes'),'fontsize',MyFontSize, 'fontweight', 'normal');
        set(findall(gcf,'type','text'),'fontSize',MyFontSize, 'fontweight', 'normal');                
        axis([min(vYear) EndYear old_ylim]);
        print('-depsc', ['../Graphs/Descriptive' VarNames{p_idx}  Suffix '.eps']);
        close(gcf);
        
        disp([VarLabelsShort{p_idx} num2str( DescriptiveCorr(1)^2, ': residual share=%1.2f :: ') num2str( DescriptiveCorr(2)^2, 'FE share=%1.2f')])
    % Periphery VAR
    else
        CurrentData = eval([PerNames{p_idx-p} '(vYear<=EndYear, :);']);
        YearFE = nanmean(CurrentData,2);
        MSAFE = nanmean(CurrentData-YearFE*ones(1,size(CurrentData,2)),1);
        PredictedData_MSA = ones(size(CurrentData,1),1)*MSAFE ;
        PredictedData_Year = YearFE*ones(1,size(CurrentData,2));
        ResidualData = CurrentData - PredictedData_MSA - PredictedData_Year;
        DescriptiveCorr = [corr(ResidualData(~isnan(ResidualData(:))), CurrentData(~isnan(ResidualData(:)))) , ...
                            corr(PredictedData_MSA(~isnan(ResidualData(:))), CurrentData(~isnan(ResidualData(:)))), corr(PredictedData_Year(~isnan(ResidualData(:))), CurrentData(~isnan(ResidualData(:))))];

        eval(['hndl=plot(vYear(vYear<=EndYear), prctile(' PerNames{p_idx-p} '(vYear<=EndYear, :), [5 25 50 75 95],2)''*PerScale(p_idx-p), ''LineWidth'', 2);']);
        title(['Distribution of ' PerLabelsShort{p_idx-p} ])
        set(hndl, {'Color'}, {[0.5 0 0];[0.5 0 0];[0 0 0]; [0.5 0 0];[0.5 0 0]},{'LineStyle'}, {':';'--';'-';'--';':'});
        figure;
        eval(['hndl=plot(vYear(vYear<=EndYear), prctile(' PerNames{p_idx-p} '(vYear<=EndYear, :), [5 25 50 75 95],2)''*PerScale(p_idx-p), ''LineWidth'', 2);']);
%         title(['Distribution of ' PerLabelsShort{p_idx-p} '. Fixed effect R-squared=' num2str( DescriptiveCorr(2)^2, '%1.2f')])
%         title({['Distribution of ' PerLabelsShort{p_idx-p} '.']; {['MSA FE R-squared=' num2str( DescriptiveCorr(2)^2, '%1.2f') '.'];['Year FE R-squared=' num2str( DescriptiveCorr(3)^2, '%1.2f') '.']}})
        title({['MSA FE R-squared=' num2str( DescriptiveCorr(2)^2, '%1.2f') '.'];['Year FE R-squared=' num2str( DescriptiveCorr(3)^2, '%1.2f') '.']})
            
        set(hndl, {'Color'}, {[0.5 0 0];[0.5 0 0];[0 0 0]; [0.5 0 0];[0.5 0 0]},{'LineStyle'}, {':';'--';'-';'--';':'});
        old_ylim=ylim;
        ylabel('%');
        grid on;
        set(findall(gcf,'type','axes'),'fontsize',MyFontSize, 'fontweight', 'normal');
        set(findall(gcf,'type','text'),'fontSize',MyFontSize, 'fontweight', 'normal');           
        axis([min(vYear) EndYear old_ylim]);
        print('-depsc', ['../Graphs/Descriptive' PerNames{p_idx-p}  Suffix '.eps']);
        close(gcf);

        disp([PerLabelsShort{p_idx-p} num2str( DescriptiveCorr(1)^2, ': residual share=%1.2f :: ') num2str( DescriptiveCorr(2)^2, 'FE share=%1.2f')])

    end
    axis([min(vYear) EndYear ylim]);
    hold all;
    grid on;
end


if exist('vjob_creation_rate_births_lvl', 'var')
    CurrentData = eval(['vjob_creation_rate_births_lvl(vYear<=EndYear, :);']);
    YearFE = nanmean(CurrentData,2);
    MSAFE = nanmean(CurrentData-YearFE*ones(1,size(CurrentData,2)),1);

    PredictedData_MSA = ones(size(CurrentData,1),1)*MSAFE ;
    PredictedData_Year = YearFE*ones(1,size(CurrentData,2));
    ResidualData = CurrentData - PredictedData_MSA - PredictedData_Year;
    DescriptiveCorr = [corr(ResidualData(~isnan(ResidualData(:))), CurrentData(~isnan(ResidualData(:)))) , ...
                        corr(PredictedData_MSA(~isnan(ResidualData(:))), CurrentData(~isnan(ResidualData(:)))), corr(PredictedData_Year(~isnan(ResidualData(:))), CurrentData(~isnan(ResidualData(:))))];

    disp(['JC births (lvl)' num2str( DescriptiveCorr(1)^2, ': residual share=%1.2f :: ') num2str( DescriptiveCorr(2)^2, 'FE share=%1.2f :: ') num2str( DescriptiveCorr(3)^2, 'Year FE share=%1.2f')])

    figure;
    eval(['hndl=plot(vYear(vYear<=EndYear), prctile(vjob_creation_rate_births_lvl(vYear<=EndYear, :), [5 25 50 75 95],2)''*1, ''LineWidth'', 2);']);
    title({['MSA FE R-squared=' num2str( DescriptiveCorr(2)^2, '%1.2f') '.'];['Year FE R-squared=' num2str( DescriptiveCorr(3)^2, '%1.2f') '.']})
    set(hndl, {'Color'}, {[0.5 0 0];[0.5 0 0];[0 0 0]; [0.5 0 0];[0.5 0 0]},{'LineStyle'}, {':';'--';'-';'--';':'});
        old_ylim=ylim;
        ylabel('%');
        grid on;
        set(findall(gcf,'type','axes'),'fontsize',MyFontSize, 'fontweight', 'normal');
        set(findall(gcf,'type','text'),'fontSize',MyFontSize, 'fontweight', 'normal');              
        axis([min(vYear) EndYear old_ylim]);
    print('-depsc', ['../Graphs/Descriptivevjob_creation_rate_births_lvl' Suffix '.eps']);
end
