function summarizeData
%{
Generate useful summary statistics about growth

Dick Startz
September 2016
revised April 2017
revised February 2019
%}
newLine = char(10);
tic;
load PWT90Data;

disp(['Number of countries = ',num2str(nCountries)]);
disp(['Total number of observations = ',num2str(nobs)]);
dYStar = log(yStar(2:end))...
    - log(yStar(1:end-1));
meanUSAGrowth = mean(dYStar);
stdUSAMeanGrowth = std(dYStar)/sqrt(length(dYStar)-1); % for use below
disp(['Mean US growth = ',num2str(meanUSAGrowth)]);

obs1970 = find(USAyear==1970);
lastUSAOb = length(USApop); % to get a round Matlab confusion

dYStar70 = log(yStar(obs1970:lastUSAOb))...
    - log(yStar(obs1970-1:lastUSAOb-1));
meanUSAGrowth70 = mean(dYStar70);
disp(['Mean US growth since 1970 = ',num2str(meanUSAGrowth70)]);

disp(['First year of data = ',num2str(min(firstYearForCountry))]);
disp(['Last year of data = ',num2str(max(lastYearForCountry))]);
disp(['Last USA per capita GDP = ',num2str(yStarLast)]);

% now get mean growth and related statistics for each country
meanGrowth = nan(nCountries,1);
stdErrMeanGrowth = nan(nCountries,1);
pValueNegative = nan(nCountries,1);
pValueVersusUSA = nan(nCountries,1);

for iCountry = 1:nCountries
    lastOb = lastObForCountry(iCountry);
    firstOb = firstObForCountry(iCountry);
    nobs = lastOb-firstOb+1;
    annualGrowth = log(y(firstOb+1:lastOb)) - log(y(firstOb:lastOb-1));
    meanGrowth(iCountry) = mean(annualGrowth);
    stdErrMeanGrowth(iCountry) = std(annualGrowth)/sqrt(nobs-1);
    tStatNegative = meanGrowth(iCountry)/stdErrMeanGrowth(iCountry);
    pValueNegative(iCountry) = tcdf(tStatNegative,nobs-1);
    tVersusUSA = (meanGrowth(iCountry)-meanUSAGrowth)/...
        sqrt(stdUSAMeanGrowth^2+stdErrMeanGrowth(iCountry)^2);
    pValueVersusUSA(iCountry) = tcdf(tVersusUSA,nobs-1);
end
t = table(meanGrowth,pValueNegative,pValueVersusUSA,...
    'RowNames',countryNames)
countNegativeGrowth = sum(meanGrowth<0);
countSignificantNegativeGrowth = sum(pValueNegative<0.05);

countBelowFrontierGrowth = sum(meanGrowth<meanUSAGrowth);
countSignificantBelowFrontierGrowth = sum(pValueVersusUSA<0.05);

disp(['Number below zero = ',num2str(countNegativeGrowth)]);
disp(['Number significantly below zero = ',...
    num2str(countSignificantNegativeGrowth)]);
disp(['Number below frontier = ',num2str(countBelowFrontierGrowth)]);
disp(['Number significantly below frontier = ',...
    num2str(countSignificantBelowFrontierGrowth)]);

%{
    Compute average growth for those countries with growth below
    frontier growth
%}
growthDif = meanGrowth - meanUSAGrowth;
shortFall = mean(growthDif(growthDif<0));
averageFinalGap = mean(yLast-yStarLast);
disp(['Short fall in growth from frontier for those below frontier'...
    ' averages ',num2str(shortFall)]);
disp(['Of those, average final gap was ',num2str(averageFinalGap)]);
disp(['Suggesting rho_d = ',num2str(1-shortFall/averageFinalGap)]);


disp(['World GDP final year = ',num2str(worldGDPLast)]);
fractionGDPWithNegativeGrowth = sum(gdpLast(meanGrowth<0))/worldGDPLast;

fractionBelowFrontierGrowth = ...
    sum(gdpLast(meanGrowth<meanUSAGrowth))/worldGDPLast;

fractionGDPSignificantNegativeGrowth = ...
    sum(gdpLast(pValueNegative<0.05))/worldGDPLast;

fractionSignificantBelowFrontierGrowth = ...
    sum(gdpLast(pValueVersusUSA<0.05))/worldGDPLast;
disp(['Fraction of GDP below zero = ',num2str(fractionGDPWithNegativeGrowth)]);
disp(['Fraction GDP below frontier = ',...
    num2str(fractionBelowFrontierGrowth)]);

disp(['Fraction of GDP significantly below zero = ',...
    num2str(fractionGDPSignificantNegativeGrowth)]);
disp(['Fraction GDP significantly below frontier = ',...
    num2str(fractionSignificantBelowFrontierGrowth)]);
disp('Growth significantly below frontier and fraction of world GDP');
for iCountry = 1:nCountries
    if pValueVersusUSA(iCountry)<0.05
        disp(['GDP 2014 ',countryNames{iCountry},' = ',...
            num2str(gdpLast(iCountry)), ' fraction of world = ',...
            num2str(gdpLast(iCountry)/worldGDPLast)]);
    end
end

disp(' '); % leave a blank line

% find some summary statistics
disp(['Mean first year = ',num2str(mean(firstYearForCountry))]);
disp(['Min first year = ',num2str(min(firstYearForCountry))]);
disp(['Max first year = ',num2str(max(firstYearForCountry))]);

weightedFirstYear = weightedMean(firstYearForCountry,popLast);
disp(['weighted mean ROW first year = ',num2str(sum(weightedFirstYear))]);

disp(['Mean mean Growth = ',num2str(mean(meanGrowth))]);
disp(['Min mean Growth = ',num2str(min(meanGrowth))]);
disp(['Max mean Growth = ',num2str(max(meanGrowth))]);

weightedMeanGrowth = weightedMean(meanGrowth,popLast);
disp(['weighted mean ROW mean Growth = ',num2str(sum(weightedMeanGrowth))]);

disp(['USA GDP Per capita final year = ',num2str(yStarLast)]);

disp(['Mean first year income = ',num2str(mean(yFirst))]);
disp(['Min first year income = ',num2str(min(yFirst))]);
disp(['Max first year income = ',num2str(max(yFirst))]);

disp(['Mean last year income = ',num2str(mean(yLast))]);
disp(['Min last year income = ',num2str(min(yLast))]);
disp(['Max last year income = ',num2str(max(yLast))]);

weightedGDPPerCapitaLast...
    = weightedMean(yLast,popLast);
disp(['weighted mean last year income = ',...
    num2str(weightedGDPPerCapitaLast)]);

% now get ratios
disp(['Mean unweighted first year ratio = ',...
    num2str(mean(yFirst/yStarFirst))]);
disp(['Min first year ratio = ',...
    num2str(min(yFirst/yStarFirst))]);
disp(['Max first year ratio = ',...
    num2str(max(yFirst/yStarFirst))]);

weightedRatioFirst = weightedFirstYear/yStarFirst(1);
disp(['weighted mean first year ratio = ',...
    num2str(weightedRatioFirst)]);

disp(['Mean unweighted last year ratio = ',...
    num2str(mean(yLast/yStarLast))]);
disp(['Min last year ratio = ',...
    num2str(min(yLast/yStarLast))]);
disp(['Max last year ratio = ',...
    num2str(max(yLast/yStarLast))]);

meanWeightedRatioLast = mean(weightedGDPPerCapitaLast/yStarLast);
disp(['Mean weighted last year ratio = ', num2str(meanWeightedRatioLast)]);
% now let's see who those high-income countries are
[yLastSorted,indexSorted] = sort(yLast,'descend');
disp('Last year income');
table(countryNames(indexSorted),yLastSorted,yLastSorted/yStarLast,...
    'VariableNames',{'Country','Last_Income','Last_income_ratio'})

disp('');
% now figure growth rate to converge in T years
T = [20, 50, 100];
for iT = 1:length(T)
    requiredGrowth ...
        = (1/meanWeightedRatioLast)^(1/T(iT))*(1+meanUSAGrowth70) - 1;
    disp(['Growth rate to converge in ',num2str(T(iT)),' years = ',...
        num2str(requiredGrowth)]);
end

% make figure showing weighted fraction of population observed 
% and weighted average gap
nYearsInSample = lastYear - firstYear + 1; 
countryInSampleInYear = nan(nCountries,nYearsInSample);
gapInSampleInYear = nan(nCountries,nYearsInSample);

for iCountry = 1:nCountries
    firstOb = firstObForCountry(iCountry);
    lastOb = lastObForCountry(iCountry);
    firstYearCountry = firstYearForCountry(iCountry);
    lastYearCountry = lastYearForCountry(iCountry); 
    rangeInSample = ...
        firstYearCountry-firstYear+1:lastYearCountry-firstYear+1;
    countryInSampleInYear(iCountry,rangeInSample) = 1;
    gapInSampleInYear(iCountry, rangeInSample) = exp(gap(firstOb:lastOb));
end
weightedFractionPopulation = nan(nYearsInSample,1);
weightedGap = nan(nYearsInSample,1);
unweightedMeanGap = nan(nYearsInSample,1); 
for iYear = 1:nYearsInSample
    whichCountriesValid = ~isnan(countryInSampleInYear(:,iYear));
    weightedFractionPopulation(iYear) = ...
        sum(popLast(whichCountriesValid))/popROWLast;
    weightedGap(iYear) = weightedMean...
        (gapInSampleInYear(whichCountriesValid,iYear),...
        popLast(whichCountriesValid));
    unweightedMeanGap(iYear) = mean...
        (gapInSampleInYear(whichCountriesValid,iYear));
end
figure;
yyaxis left;
plot(firstYear:lastYear, weightedFractionPopulation, 'linewidth',1.5);
v = axis;
axis([firstYear,lastYear,v(3),v(4)]);
ylabel(['Coverage of world population by ',num2str(lastYear),...
    ' weights']);
xlabel('Year');

yyaxis right;
plot(firstYear:lastYear,weightedGap,...
    'linewidth',1.5);
yLabelHandle = ylabel('Population weighted average output per capita ratio','rot',-90);
yLabelHandle.Position(1)=yLabelHandle.Position(1)+2; % make it look nicer
title({'Sample coverage of population and average ratio'...
    'output per capita-to-frontier excluding USA',...
    ['Weighted by ',num2str(lastYear),' population']});
print('Population and income ratio','-depsc2','-tiff');
print('Population and income ratio','-dpng');

table((firstYear:lastYear)',weightedFractionPopulation,...
    weightedGap, unweightedMeanGap)

% now do back-of-the-envelope
% current gap will be cut to a quarter with 2 chances in 3.

finalGapWouldBe = (1/3)*(1-meanWeightedRatioLast)...
    + (2/3)*(1-meanWeightedRatioLast)/4;
finalRatioWouldBe = 1 - finalGapWouldBe;

extraGrowth = (finalRatioWouldBe/meanWeightedRatioLast)^.01;
disp('Back of the envelope calculation of 2/3rds chance of reducing');
disp(['gap to 1/4 gives 100 year annual growth rate of ',...
    num2str(extraGrowth)]);

save PWT90Data;
elapsedTime = toc;
disp(['Executed ',mfilename,' in ',num2str(elapsedTime),' seconds ',...
    datestr(now,'on mmmm-dd-yyyy at HH:MM:SS')]);
end

