function reportOn3StateFutureSimulationsFromPriors
%{

Create output needed for paper that comes from future simulations
based on priors
Dick Startz
July 2016
modified September 2016
modified April 2017
modified May 2017
modified January 2019
%}
tic;
futureFileName = 'simulated3StateFuturesFrom2014SplitFromPriors';

load(futureFileName);

% 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,:)'*popWeightROWLast(iCountry);
end
medSimulatedROWIncomeFinal = medWeightedROWIncome(end);
disp(['Final ROW Income = ',num2str(medSimulatedROWIncomeFinal)]);

medWeightedWorldIncome...
    = medSimulatedFutureUSAIncome*popWeightUSALast;

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

% now without convergence
medWeightedWorldIncomeNoConvergence...
    = medSimulatedFutureUSAIncome*popWeightUSALast;

medIncomesNoConvergence = nan(nCountries,nYears);
for iCountry = 1:nCountries
    lastOb = lastObForCountry(iCountry);
    medIncomesNoConvergence(iCountry,:)...
        = exp(gap(lastOb))*medSimulatedFutureUSAIncome;
    medWeightedWorldIncomeNoConvergence =...
        medWeightedWorldIncomeNoConvergence + ...
        medIncomesNoConvergence(iCountry,:)'*popWeightWorldLast(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)'});
plot(lastYear,weightedWorldIncomeLast/1000,'og','linewidth',1.5);
legend([p1 p2 p3],'with estimated convergence','without convergence',...
    '2014 US income','location','Northwest');
hold off;
print('withAndWithoutConvergenceFromPriors','-depsc2','-tiff');
print('withAndWithoutConvergenceFromPriors','-dpng');



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),...
        populationWeights);
    weightedAverageSimPrc05(iYear) =...
        weightedMean(prc05SimulatedFutureLevelIncome(:,iYear),...
        populationWeights);
    weightedAverageSimPrc95(iYear) =...
        weightedMean(prc95SimulatedFutureLevelIncome(:,iYear),...
        populationWeights);
end

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

figure;
p1 = semilogy(yearsToPlot,weightedAverageSimMedian,'-r','linewidth',1.5);
hold on;


title({'Population weighted income: median and 90% HPD',...
    'thousands of 2011 dollars'});
semilogy(yearsToPlot,weightedAverageSimPrc05,':r','linewidth',1.5);
semilogy(yearsToPlot,weightedAverageSimPrc95,':r','linewidth',1.5);

realFutureFileName = 'simulated3StateFuturesFrom2014Split';
load(realFutureFileName);

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);

% 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,:)'*popWeightROWLast(iCountry);
end
medSimulatedROWIncomeFinal = medWeightedROWIncome(end);
disp(['Final ROW Income = ',num2str(medSimulatedROWIncomeFinal)]);

medWeightedWorldIncome...
    = medSimulatedFutureUSAIncome*popWeightUSALast;

for iCountry = 1:nCountries
    medWeightedWorldIncome = medWeightedWorldIncome...
        + median(exp(squeeze(simulatedFutureLogIncome(iCountry,:,:))),2)...
        *popWeightWorldLast(iCountry);
end
medSimulatedWorldIncomeFinal = medWeightedWorldIncome(end);

medSimulatedFutureLevelIncome = median(simulatedFutureLevelIncome,3);
prc05SimulatedFutureLevelIncome = prctile(simulatedFutureLevelIncome,5,3);
prc95SimulatedFutureLevelIncome = prctile(simulatedFutureLevelIncome,95,3);

weightedAverageSimMedianReal = zeros(nYears,1);
weightedAverageSimPrc05Real = zeros(nYears,1);
weightedAverageSimPrc95Real = zeros(nYears,1);

populationWeights = [pop(lastObForCountry);popUSALast];
for iYear = 1:nYears
    weightedAverageSimMedianReal(iYear) =...
        weightedMean(medSimulatedFutureLevelIncome(:,iYear),...
        populationWeights);
    weightedAverageSimPrc05Real(iYear) =...
        weightedMean(prc05SimulatedFutureLevelIncome(:,iYear),...
        populationWeights);
    weightedAverageSimPrc95Real(iYear) =...
        weightedMean(prc95SimulatedFutureLevelIncome(:,iYear),...
        populationWeights);
end
p4 = semilogy(yearsToPlot,weightedAverageSimMedianReal,'-b','linewidth',1.5);
semilogy(yearsToPlot,weightedAverageSimPrc05Real,':b','linewidth',1.5);
semilogy(yearsToPlot,weightedAverageSimPrc95Real,':b','linewidth',1.5);

legend([p1, p4],'From Prior','From Posterior','location','Northwest');
%{
v = axis;
axis([v(1) v(2) 5 1000]);
%}
ax = gca;
%ax.YAxis.TickLabels ={'10','100','1,000'};
ylabel('Income (log scale)');
hold off;
print('Projection From Priors','-depsc2','-tiff');
print('Projection From Priors','-dpng');

figure;
p1 = plot(yearsToPlot,weightedAverageSimMedian,'-r','linewidth',1.5);
hold on;


title({'Population weighted income: median and 90% HPD',...
    'thousands of 2011 dollars'});
plot(yearsToPlot,weightedAverageSimPrc05,':r','linewidth',1.5);
plot(yearsToPlot,weightedAverageSimPrc95,':r','linewidth',1.5);

p4 = plot(yearsToPlot,weightedAverageSimMedianReal,'-b','linewidth',1.5);
plot(yearsToPlot,weightedAverageSimPrc05Real,':b','linewidth',1.5);
plot(yearsToPlot,weightedAverageSimPrc95Real,':b','linewidth',1.5);

legend([p1, p4],'From Prior','From Posterior','location','Northwest');
%{
v = axis;
axis([v(1) v(2) 5 1000]);
%}
ax = gca;
%ax.YAxis.TickLabels ={'10','100','1,000'};
ylabel('Income');
hold off;
print('Projection From Priors - linear','-depsc2','-tiff');
print('Projection From Priors - linear','-dpng');

return;
% 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)*popWeightROWLast(iCountry);
    prc05WeightedROWGap = prc05WeightedROWGap...
        + prctile(gaps,5,2)*popWeightROWLast(iCountry);
    prc95WeightedROWGap = prc95WeightedROWGap...
        + prctile(gaps,95,2)*popWeightROWLast(iCountry);
    prc10WeightedROWGap = prc10WeightedROWGap...
        + prctile(gaps,10,2)*popWeightROWLast(iCountry);
    prc90WeightedROWGap = prc90WeightedROWGap...
        + prctile(gaps,90,2)*popWeightROWLast(iCountry);
    for iProb = 1:3
        meanWeightedROWStates(:,iProb) = meanWeightedROWStates(:,iProb)...
           + mean(squeeze(simulatedFutureStates(iCountry,:,:))==iProb,2)...
            *popWeightROWLast(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');

%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 income 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('Years');
ylabel('Median fraction of frontier');
title('Population weighted ratio of income to frontier');
hold off;

subplot(2,1,2);
plot(yearsToPlot,meanWeightedROWStates,'linewidth',1.5);
hold on;
plot(lastYear,meanStatesROWLast,'og','linewidth',1.5);
ylabel('Probability of state');
title('Mean probability of state');
legend('n','c','d');
hold off;

print('Future Gap and States 3-states From Priors','-depsc2','-tiff');
print('Future Gap and States 3-states From Priors','-dpng');

% 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 = [medFutureLevelIncome(:,end);...
    medSimulatedFutureUSAIncomeFinal];
[sortedIncomeFinal,indexSortedIncomeFinal] = sort(incomeFinal);
sortedPopWeightWorldFinal = popWeightAllWorldLast(indexSortedIncomeFinal);

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

xlabel('Fraction of world population');
ylabel('Income, thousands of 2011 dollars');
title('Cumulative income 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 From Priors','-depsc2','-tiff');
print('Cumulative income distribution 3-states From Priors','-dpng');

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

% now make bubble plots
figure;
fudge = 20/.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*popWeightROWLast,'g','d','linewidth',1);
legend('2014','1970');
v = axis;
axis([v(1) 1 v(3) v(4)]);
hold off;
print -depsc2 'bubble convergent 1970 From Priors';
print -dpng 'bubble convergent 1970 From Priors';

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

