function reportOnWhichCountriesMatter
%{

Figure out which of the countries account for most of closing
the gap.

Basically figure out predicted median growth for a country versus
what it would be if at last observed gap. Then sort these from largest
to smallest.

Dick Startz
July 2016
modified September 2016
modified April 2017 to read in right file
modifed May 2017
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...
        + yLast(iCountry)*popWeightROWLast(iCountry);
end


weightedWorldIncomeLast = yStarLast*popWeightUSALast;
for iCountry = 1:nCountries
    weightedWorldIncomeLast = weightedWorldIncomeLast...
        + yLast(iCountry)*popWeightWorldLast(iCountry);
end


medSimulatedFutureUSAIncome = median(exp(simulatedFutureLogUSA),2);
medSimulatedFutureUSAIncomeFinal = medSimulatedFutureUSAIncome(end);

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

% now for each country figure out the difference between its forecast
% and non-convergence terminal value

diffLast = medFutureLevelIncome(:,end) - medIncomesNoConvergence(:,end);
% now weighted difference
weightedDiffLast = diffLast.*popWeightWorldLast;

disp(['Total difference between converged and non-converged = ',...
    num2str(sum(weightedDiffLast))]);

[sortedweightedDiffLast,indexSorted] = sort(weightedDiffLast,'descend');
percentageLost = sortedweightedDiffLast/medSimulatedWorldIncomeFinal;
lastWeight = popWeightWorldLast(indexSorted);
finalWeight = popWeight(indexSorted);
t1 = table(percentageLost, lastWeight,finalWeight,...
    'RowNames',countryNames(indexSorted))



elapsedTime = toc;

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

