function reportOn3StateFutureSimulationsNewPop
%{

Create output needed for paper that comes from future simulations
Use Adrian's population projections

Dick Startz
July 2016
modified September 2016
modified April 2017
modified May 2017
modified January 2019
modified February 2019
%}
tic;
futureFileName = 'simulated3StateFuturesFrom2014Split';
load(futureFileName);
load('futurePopulation');
popWeight = futurePopulation./repmat(sum(futurePopulation,1),...
    size(futurePopulation,1),1);

% get initial world income
weightedROWIncomeLast = 0;
for iCountry = 1:nCountries
    weightedROWIncomeLast = weightedROWIncomeLast...
        + y(lastObForCountry(iCountry))*popWeightROWLast(iCountry);
end
disp(['Last observed ROW income = ', num2str(weightedROWIncomeLast)]);


weightedWorldIncomeLast = yStarLast*popWeightUSALast;
for iCountry = 1:nCountries
    weightedWorldIncomeLast = weightedWorldIncomeLast...
        + y(lastObForCountry(iCountry))*popWeightWorldLast(iCountry);
end
disp(['Last observed world income = ',...
    num2str(weightedWorldIncomeLast)]);
disp(['Last observed USA income = ', num2str(yStarLast)]);
disp(['Last observed gap = ', num2str(weightedROWIncomeLast/yStarLast)]);



medSimulatedFutureUSAIncome = median(exp(simulatedFutureLogUSA),2);
medSimulatedFutureUSAIncomeFinal = medSimulatedFutureUSAIncome(end);
disp(['Final USA Income = ',num2str(medSimulatedFutureUSAIncomeFinal)]);

% find median simulations for rest of world and entire world
% we do this by adding up weighted medians for each country

medFutureLevelIncome = nan(nCountries,nYears);
medWeightedROWIncome = zeros(nYears,1);
for iCountry = 1:nCountries
    medFutureLevelIncome(iCountry,:) =...
        median(exp(squeeze(simulatedFutureLogIncome(iCountry,:,:))),2);
    medWeightedROWIncome = medWeightedROWIncome...
        + (medFutureLevelIncome(iCountry,:).*popWeight(iCountry,:))';
end
medSimulatedROWIncomeFinal = medWeightedROWIncome(end);
disp(['Final ROW Income = ',num2str(medSimulatedROWIncomeFinal)]);

medWeightedWorldIncome...
    = medSimulatedFutureUSAIncome.*popWeight(end,:)';

for iCountry = 1:nCountries
    medWeightedWorldIncome = medWeightedWorldIncome...
        + (median(exp(squeeze(simulatedFutureLogIncome(iCountry,:,:))),2)...
        .*popWeight(iCountry,:)');
end
medSimulatedWorldIncomeFinal = medWeightedWorldIncome(end);
disp(['Final World Income = ',num2str(medSimulatedWorldIncomeFinal)]);

% now without convergence
medWeightedWorldIncomeNoConvergence...
    = medSimulatedFutureUSAIncome.*(popWeight(end,:)');

medIncomesNoConvergence = nan(nCountries,nYears);
for iCountry = 1:nCountries
    lastOb = lastObForCountry(iCountry);
    medIncomesNoConvergence(iCountry,:)...
        = exp(gap(lastOb))*medSimulatedFutureUSAIncome;
    medWeightedWorldIncomeNoConvergence =...
        medWeightedWorldIncomeNoConvergence + ...
        (medIncomesNoConvergence(iCountry,:).*popWeight(iCountry,:))';
end
medWeightedWorldIncomeNoConvergenceFinal...
    = medWeightedWorldIncomeNoConvergence(end);
disp(['Final world income no convergence = ',...
    num2str(medWeightedWorldIncomeNoConvergenceFinal)]);

disp(['Ratio of final world income to last USA income,',...
    ' no convergence = ',...
    num2str(medWeightedWorldIncomeNoConvergenceFinal/yStarLast)]);

disp(['Ratio of final world income to last USA income,',...
    ' with convergence = ',...
    num2str(medSimulatedWorldIncomeFinal/yStarLast)]);

disp(['Ratio of final USA income to last USA income = ',...
    num2str(medSimulatedFutureUSAIncomeFinal/yStarLast)]);

    

yearsToPlot = (lastYear+1):(lastYear+nYears);
figure;
p1 = plot(yearsToPlot,medWeightedWorldIncome/1000,'-b',...
    'linewidth',1.5);
hold on;
p2 = plot(yearsToPlot,medWeightedWorldIncomeNoConvergence/1000,...
    '--m','linewidth',1.5);
p3 = plot(yearsToPlot,yStarLast*ones(length(yearsToPlot),1)/1000,':k',...
    'linewidth',1.5);
title({'Estimated world GDP per capita with estimated convergence'...
    'and without convergence', '(thousands of 2011 dollars)'});
xlabel('Year');
ylabel('Output per capita');
plot(lastYear,weightedWorldIncomeLast/1000,'og','linewidth',1.5);
legend([p1 p2 p3],'with estimated convergence','without convergence',...
    '2014 US output per capita','location','Northwest');
hold off;
print('withAndWithoutConvergence new pop','-depsc2','-tiff');
print('withAndWithoutConvergence new pop','-dpng');
print -depsc -tiff 'Figure_2_color';
print -deps -tiff 'Figure_2_BW';



expectedGrowthRate =...
    (medSimulatedWorldIncomeFinal/weightedWorldIncomeLast)^(1/nYears)-1;
disp(['Expected world income growth rate over ',num2str(nYears),' = ',...
    num2str(expectedGrowthRate)]);
gapLevelFinal =...
    medSimulatedROWIncomeFinal/medSimulatedFutureUSAIncomeFinal;
disp(['Final gap over ',num2str(nYears),' years = ',...
    num2str(gapLevelFinal)]);
disp(['Ratio of final world income to current US income',...
    ' with convergence = ', num2str(medSimulatedWorldIncomeFinal/...
    medSimulatedFutureUSAIncomeFinal)]);

% now plot HPDs with and without parameter uncertainty
% add USA to end of list
simulatedFutureLevelIncome =...
    exp([simulatedFutureLogIncome;reshape(simulatedFutureLogUSA,...
    [1,nYears,nDraws])])/1000;
% now take the median and prctiles for each country,
% then weight them by population
medSimulatedFutureLevelIncome = median(simulatedFutureLevelIncome,3);
prc05SimulatedFutureLevelIncome = prctile(simulatedFutureLevelIncome,5,3);
prc95SimulatedFutureLevelIncome = prctile(simulatedFutureLevelIncome,95,3);

weightedAverageSimMedian = zeros(nYears,1);
weightedAverageSimPrc05 = zeros(nYears,1);
weightedAverageSimPrc95 = zeros(nYears,1);

%populationWeights = [pop(lastObForCountry);popUSALast];
for iYear = 1:nYears
    weightedAverageSimMedian(iYear) =...
        weightedMean(medSimulatedFutureLevelIncome(:,iYear),...
        popWeight(:,iYear));
    weightedAverageSimPrc05(iYear) =...
        weightedMean(prc05SimulatedFutureLevelIncome(:,iYear),...
        popWeight(:,iYear));
    weightedAverageSimPrc95(iYear) =...
        weightedMean(prc95SimulatedFutureLevelIncome(:,iYear),...
        popWeight(:,iYear));
end

figure;
plot(yearsToPlot,weightedAverageSimMedian,'-b','linewidth',1.5);
hold on;
title({'Population weighted median output per capita and 90% HPD',...
    'thousands of 2011 dollars'});
plot(yearsToPlot,weightedAverageSimPrc05,':b','linewidth',1.5);
plot(yearsToPlot,weightedAverageSimPrc95,':b','linewidth',1.5);
xlabel('Year');
hold off;

% now do the same thing with no parameter uncertainty
simulatedFutureCertainLevelIncome =...
    exp([simulatedFutureCertainLogIncome;reshape(simulatedFutureCertainLogUSA,...
    [1,nYears,nDraws])])/1000;
% now take the median and prctiles for each country,
% then weight them by population
medSimulatedFutureCertainLevelIncome =...
    median(simulatedFutureCertainLevelIncome,3);
prc05SimulatedFutureCertainLevelIncome =...
    prctile(simulatedFutureCertainLevelIncome,5,3);
prc95SimulatedFutureCertainLevelIncome =...
    prctile(simulatedFutureCertainLevelIncome,95,3);

weightedAverageSimMedianCertain = zeros(nYears,1);
weightedAverageSimPrc05Certain = zeros(nYears,1);
weightedAverageSimPrc95Certain = zeros(nYears,1);

for iYear = 1:nYears
    weightedAverageSimMedianCertain(iYear) =...
        weightedMean(medSimulatedFutureCertainLevelIncome(:,iYear),...
        popWeight(:,iYear));
    weightedAverageSimPrc05Certain(iYear) =...
        weightedMean(prc05SimulatedFutureCertainLevelIncome(:,iYear),...
        popWeight(:,iYear));
    weightedAverageSimPrc95Certain(iYear) =...
        weightedMean(prc95SimulatedFutureCertainLevelIncome(:,iYear),...
        popWeight(:,iYear));
end

% now same thing without shocks
simulatedFutureLevelIncomeNoShocks =...
    exp([simulatedFutureLogIncomeNoSchocks;...
    reshape(simulatedFutureLogUSANoSchock,...
    [1,nYears,nDraws])])/1000;
% now take the median and prctiles for each country,
% then weight them by population
medSimulatedFutureLevelIncomeNoShocks =...
    median(simulatedFutureLevelIncomeNoShocks,3);
prc05SimulatedFutureLevelIncomeNoShocks =...
    prctile(simulatedFutureLevelIncomeNoShocks,5,3);
prc95SimulatedFutureLevelIncomeNoShocks =...
    prctile(simulatedFutureLevelIncomeNoShocks,95,3);

weightedAverageSimMedianNoShocks = zeros(nYears,1);
weightedAverageSimPrc05NoShocks = zeros(nYears,1);
weightedAverageSimPrc95NoShocks = zeros(nYears,1);

for iYear = 1:nYears
    weightedAverageSimMedianNoShocks(iYear) =...
        weightedMean(medSimulatedFutureLevelIncomeNoShocks(:,iYear),...
        popWeight(:,iYear));
    weightedAverageSimPrc05NoShocks(iYear) =...
        weightedMean(prc05SimulatedFutureLevelIncomeNoShocks(:,iYear),...
        popWeight(:,iYear));
    weightedAverageSimPrc95NoShocks(iYear) =...
        weightedMean(prc95SimulatedFutureLevelIncomeNoShocks(:,iYear),...
        popWeight(:,iYear));
end



figure;
plot(yearsToPlot,weightedAverageSimMedianCertain,'-b','linewidth',1.5);
hold on;
title({'Population weighted median output per capita and 90% HPD',...
    'thousands of 2011 dollars, without parameter uncertainty'});
plot(yearsToPlot,weightedAverageSimPrc05Certain,':b','linewidth',1.5);
plot(yearsToPlot,weightedAverageSimPrc95Certain,':b','linewidth',1.5);
xlabel('Year');
hold off;

figure;
plot(yearsToPlot,weightedAverageSimMedianNoShocks,'-b','linewidth',1.5);
hold on;
title({'Population weighted median output per capita and 90% HPD',...
    'thousands of 2011 dollars, no shocks to frontier'});
plot(yearsToPlot,weightedAverageSimPrc05NoShocks,':b','linewidth',1.5);
plot(yearsToPlot,weightedAverageSimPrc95NoShocks,':b','linewidth',1.5);
xlabel('Year');
hold off;

figure;
p1 = semilogy(yearsToPlot,weightedAverageSimMedian,'-b','linewidth',1.5);
hold on;
p2 = semilogy(yearsToPlot,weightedAverageSimMedianCertain,'--r','linewidth',1.5);
p3 = semilogy(yearsToPlot,weightedAverageSimMedianNoShocks,':g','linewidth',1.5);

title({'Population weighted output per capita: median and 90% HPD',...
    'thousands of 2011 dollars'});
semilogy(yearsToPlot,weightedAverageSimPrc05,'-b','linewidth',1.5);
semilogy(yearsToPlot,weightedAverageSimPrc95,'-b','linewidth',1.5);
semilogy(yearsToPlot,weightedAverageSimPrc05Certain,'--r','linewidth',1.5);
semilogy(yearsToPlot,weightedAverageSimPrc95Certain,'--r','linewidth',1.5);
semilogy(yearsToPlot,weightedAverageSimPrc05NoShocks,':g','linewidth',1.5);
semilogy(yearsToPlot,weightedAverageSimPrc95NoShocks,':g','linewidth',1.5);

legend([p1 p2 p3],'Projection','Projection w/o parameter uncertainty',...
    'Projection w/o frontier shocks','location','Northwest');
%{
v = axis;
axis([v(1) v(2) 5 1000]);
%}
ax = gca;
ax.YAxis.TickLabels ={'10','100','1,000'};
ylabel('Output per capita (log scale)');
xlabel('Year');
hold off;
print('Projection with and without uncertainty new pop','-depsc2','-tiff');
print('Projection with and without uncertainty new pop','-dpng');
print -depsc -tiff 'Figure_5_color';
print -deps -tiff 'Figure_5_BW';

% find median and hpd gap and mean state probabilities
medWeightedROWGap = zeros(nYears,1);
prc05WeightedROWGap = zeros(nYears,1);
prc95WeightedROWGap = zeros(nYears,1);
prc10WeightedROWGap = zeros(nYears,1);
prc90WeightedROWGap = zeros(nYears,1);
meanWeightedROWStates = zeros(nYears,3); 
for iCountry = 1:nCountries
    gaps = exp(squeeze(simulatedCountryLogGaps(iCountry,:,:)));
    medWeightedROWGap = medWeightedROWGap...
        + median(gaps,2).*popWeight(iCountry,:)';
    prc05WeightedROWGap = prc05WeightedROWGap...
        + prctile(gaps,5,2).*popWeight(iCountry,:)';
    prc95WeightedROWGap = prc95WeightedROWGap...
        + prctile(gaps,95,2).*popWeight(iCountry,:)';
    prc10WeightedROWGap = prc10WeightedROWGap...
        + prctile(gaps,10,2).*popWeight(iCountry,:)';
    prc90WeightedROWGap = prc90WeightedROWGap...
        + prctile(gaps,90,2).*popWeight(iCountry,:)';
    for iProb = 1:3
        meanWeightedROWStates(:,iProb) = meanWeightedROWStates(:,iProb)...
           + mean(squeeze(simulatedFutureStates(iCountry,:,:))==iProb,2)...
            .*popWeight(iCountry,:)';
    end

end
disp(['Mean ROW weighted final states = ',...
    num2str(meanWeightedROWStates(end))]);
disp(['Median ROW weighted final Gap = ',...
    num2str(medWeightedROWGap(end))]);
disp('Median percentiles final gap');
disp([num2str(5),' ',num2str(prc05WeightedROWGap(end))]);
disp([num2str(10),' ',num2str(prc10WeightedROWGap(end))]);
disp([num2str(90),' ',num2str(prc90WeightedROWGap(end))]);
disp([num2str(95),' ',num2str(prc95WeightedROWGap(end))]);

% find Last gap and state
gapROWLast = 0;
meanStatesROWLast = zeros(3,1);
meanCountry1970FirstStates = nan(nCountries,3);
meanCountryLastStates = nan(nCountries,3);
lastGap = yLast/yStarLast;
year1970Gap = nan(nCountries,1);

for iCountry = 1:nCountries
    firstOb = firstObForCountry(iCountry);
    lastOb = lastObForCountry(iCountry);
    gapROWLast = gapROWLast + lastGap(iCountry)*popWeightROWLast(iCountry);
    for iProb = 1:3
        meanStatesROWLast(iProb) = meanStatesROWLast(iProb) +...
            mean(S(lastOb,:)==iProb)*popWeightROWLast(iCountry);
        if firstYearForCountry(iCountry) <=1970
            firstObToUse = firstOb + 1970 - firstYearForCountry(iCountry);
            meanCountry1970FirstStates(iCountry,iProb)...
                = mean(S(firstObToUse+1,:)==iProb);
            year1970Gap(iCountry) = exp(gap(firstObToUse));
        end
        meanCountryLastStates(iCountry,iProb) = mean(S(lastOb,:)==iProb);
    end
end

figure;
plot(yearsToPlot,meanWeightedROWStates,'linewidth',1.5);
hold on;
plot(lastYear,meanStatesROWLast,'og','linewidth',1.5);
title('Population weighted mean state');
xlabel('Year');
legend('n','c','d','location','best');

%v = axis;
%axis([v(1) v(2) 0 1]);
hold off;

figure;
plot(yearsToPlot,medWeightedROWGap,yearsToPlot,prc05WeightedROWGap,':m',...
    yearsToPlot,prc95WeightedROWGap,':m','linewidth',1.5);
hold on;
plot(lastYear,gapROWLast,'og','linewidth',1.5);
title('Median and 90% HPD for ratio of output per capita to frontier');
xlabel('Year');
hold off;


subplot(2,1,1);
plot(yearsToPlot,medWeightedROWGap,'linewidth',1.5);
hold on;
plot(lastYear,gapROWLast,'og','linewidth',1.5);

xlabel('Year');
ylabel('Median fraction of frontier');
title('Population weighted ratio of output per capita to frontier');
hold off;

subplot(2,1,2);
plot(yearsToPlot,meanWeightedROWStates(:,1),'--b','linewidth',1.5);
hold on;
plot(yearsToPlot,meanWeightedROWStates(:,2),'-r','linewidth',1.5);
plot(yearsToPlot,meanWeightedROWStates(:,3),':y','linewidth',1.5);

plot(lastYear,meanStatesROWLast,'og','linewidth',1.5);
ylabel('Probability of state');
title('Mean probability of state');
legend('n','c','d','location','best');
xlabel('Year');
hold off;

print('Future Gap and States 3-states new pop','-depsc2','-tiff');
print('Future Gap and States 3-states new pop','-dpng');
print -depsc -tiff 'Figure_4_color';
print -deps -tiff 'Figure_4_BW';

% now make Lorenz curves
% sort income, then sort populations the same way, then cumsum, include US
popWeightAllWorldLast = [popWeightWorldLast;popWeightUSALast];
incomeLast = [yLast;yStarLast];
[sortedIncomeLast,indexSortedIncomeLast] = sort(incomeLast);
sortedPopWeightWorldLast = popWeightAllWorldLast(indexSortedIncomeLast);

incomeFinal = medSimulatedFutureLevelIncome(:,end);
[sortedIncomeFinal,indexSortedIncomeFinal] = sort(incomeFinal);
popWeightAllWorldFinal = popWeight(:,end);
sortedPopWeightWorldFinal = popWeightAllWorldFinal(indexSortedIncomeFinal);

figure;
plot(cumsum(sortedPopWeightWorldLast),sortedIncomeLast/1000,...
    ':','linewidth',1.5);
hold on;
plot(cumsum(sortedPopWeightWorldFinal),sortedIncomeFinal,...
    'linewidth',1.5);

xlabel('Fraction of world population');
ylabel('Output per capita, thousands of 2011 dollars');
title('Cumulative output distribution versus population');

plot(cumsum(sortedPopWeightWorldLast),...
    ones(length(sortedPopWeightWorldLast),1)*yStarLast/1000,'--',...
    'linewidth',0.75);
plot(cumsum(sortedPopWeightWorldFinal),...
    ones(length(sortedPopWeightWorldFinal),1)...
    *medSimulatedFutureUSAIncomeFinal/1000,'--','linewidth',0.75);


legend(num2str(lastYear),num2str(yearsToPlot(end)),...
    [num2str(lastYear),' US'],[num2str(yearsToPlot(end)),' US'],...
    'location','northwest');
hold off;
print('Cumulative income distribution 3-states new pop','-depsc2','-tiff');
print('Cumulative income distribution 3-states new pop','-dpng');
print -depsc -tiff 'Figure_6_color';
print -deps -tiff 'Figure_6_BW';


giniLast = giniCoefficient(sortedIncomeLast,sortedPopWeightWorldLast);
disp(['Last year Gini coefficient = ',num2str(giniLast)]);

% now do Gini separately for each draw
simulatedFutureLogIncomeAll = simulatedFutureLogIncome;
simulatedFutureLogIncomeAll(end+1,:,:)=simulatedFutureLogUSA;
simulatedFutureLevelIncome = exp(simulatedFutureLogIncomeAll);
giniAll = nan(nYears,nDraws);
for iDraw = 1:nDraws
    for iYear = 1:nYears
        giniAll(iYear,iDraw) = ...
            giniCoefficient(...
            squeeze(simulatedFutureLevelIncome(:,iYear,iDraw)),...
            popWeight(:,iYear));
    end
end

figure;
medianGini = median(giniAll,2);
plot(yearsToPlot,medianGini,'linewidth',1.5);
hold on;
plot(yearsToPlot,prctile(giniAll,05,2),':b','linewidth',1.5);
plot(yearsToPlot,prctile(giniAll,95,2),':b','linewidth',1.5);
plot(lastYear,giniLast,'og','linewidth',1.5);
hold off;

title('Projected Gini Coefficients and 90% HPD');
xlabel('Year');
ylabel('Gini coefficient');
print('Gini over time improved new pop','-dpng');
print -depsc -tiff 'Figure_7_Color';
print -deps -tiff 'Figure_7_BW';


disp(['First simulated year median Gini coefficient = ',...
    num2str(medianGini(1))]);
disp(['Final year median Gini coefficient = ',num2str(medianGini(end))]);

% now make bubble plots
figure;
fudge = 30/.05;
scatter(lastGap,meanCountryLastStates(:,2),fudge*popWeightROWLast,...
    'b','o','linewidth',1);
xlabel('Ratio to frontier');
ylabel('Probabilty of convergent state');
title({'Probability of convergent state vs ratio to frontier',...
    '1970 vs 2014'});
hold on;

scatter(year1970Gap,meanCountry1970FirstStates(:,2),...
    fudge*popWeight(1:end-1,end),'g','d','linewidth',1);
legend('2014','1970');
v = axis;
axis([v(1) 1 v(3) v(4)]);
hold off;
print -depsc2 'bubble convergent 1970 new pop';
print -dpng 'bubble convergent 1970 new pop';
print -depsc -tiff 'Figure_A2_color';

figure;
scatter(lastGap,meanCountryLastStates(:,1),fudge*popWeightROWLast,'b','o');
xlabel('Ratio to frontier');
ylabel('Probabilty of non-convergent state');
title({'Probability of non-convergent state vs ratio to frontier',...
    'first observed year vs 2014'});
hold on;
scatter(year1970Gap,meanCountry1970FirstStates(:,1),fudge*popWeightROWLast,...
    'g','d');
legend('2014','first year');
v = axis;
axis([v(1) 1 v(3) v(4)]);
hold off;


elapsedTime = toc;

disp(['Executed ',mfilename,' in ',num2str(elapsedTime),' seconds on ',...
    datestr(datetime)]);
end

