function summarize3StatePriorAndPosterior
%{
Summarize the priors and posteriors for both the USA
and the rest of the world Gibbs samplers from the 3-state model

Dick Startz
July 2016
modified September 2016
modified April 2017
modified May 2017
%}
tic;

nYears = 100;
hpdProb = 0.95;
nBins = 1000;

load('gibbsFrontier.mat');
load('PWT90Data.mat'); %temp
disp('Prior and posterior for US');
disp(['Prior mu0 = ',num2str(mu0)]);
disp(['Prior V0 = ',num2str(V0), ' prior sqrt(V0) = ',num2str(sqrt(V0))]);

disp(['Prior nu0 = ',num2str(nu0)]);
disp(['Prior sSqr0 = ',num2str(sSqr0)]);

disp(['Posterior median mu = ',num2str(median(mu))]);
disp(['Posterior mean mu = ',num2str(mean(mu))]);
disp(['Posterior std  mu = ',num2str(std(mu))]);

[lowerBound,upperBound] = hpd(mu,hpdProb,nBins);
disp([num2str(hpdProb),'% HPD for mu is [',...
    num2str(lowerBound),',',num2str(upperBound),']']);


disp(['Posterior median hUSA = ',num2str(median(hUSA))]);
disp(['Posterior mean hUSA = ',num2str(mean(hUSA))]);
disp(['Posterior std  hUSA = ',num2str(std(hUSA))]);
disp(['USA burnIns = ',num2str(burnInsUSA),' Gibbs keep = ',...
    num2str(gibbsKeepUSA)]);
[lowerBound,upperBound] = hpd(hUSA,hpdProb,nBins);
disp([num2str(hpdProb),'% HPD for hUSA is [',...
    num2str(lowerBound),',',num2str(upperBound),']']);
load(['C:\Users\startz\Documents\Next Hundred Years revision\'...
    'Version 2017\3-state programs pop weighted\'...
    'bayesianConvergenceEstimates3States.mat'],...
    'S');
SnoSplit= S;
clear S;



load('bayesianConvergenceEstimates3StatesSplit.mat');

disp(['Prior rho_c = U(',num2str(cLower_c),',',num2str(cUpper_c),')']);
disp(['Prior rho_d = U(',num2str(cLower_d),',',num2str(cUpper_d),')']);

% figure gap closing half life starting at initial gap
initialGap = weightedMean(yLast,popLast)/yStarLast;
halfLifeLower = gapHalfLife(cLower_c,initialGap);
halfLifeUpper = gapHalfLife(cUpper_c,initialGap);
disp(['Time to close half of initial gap of ',num2str(initialGap),...
    ' is ',num2str(halfLifeLower),' at lower limit and ',...
    num2str(halfLifeUpper),' at upper limit']);

rhoUtility(rho_n,'rho_n',cLower_n,cUpper_n,hpdProb,nBins,initialGap);
rhoUtility(rho_c,'rho_c',cLower_c,cUpper_c,hpdProb,nBins,initialGap);
rhoUtility(rho_d,'rho_d',cLower_d,cUpper_d,hpdProb,nBins,initialGap);

priorUtility(alpha_n,'alpha_n');
priorUtility(alpha_c,'alpha_c');
priorUtility(alpha_d,'alpha_d');

posteriorProbUtility(transitionMatrixEarly(1,:,:),...
    'Early from n',hpdProb,nBins);
posteriorProbUtility(transitionMatrixEarly(2,:,:),...
    'Early from c',hpdProb,nBins);
posteriorProbUtility(transitionMatrixEarly(3,:,:),...
    'Early from d',hpdProb,nBins);

posteriorProbUtility(transitionMatrix(1,:,:),'from n',hpdProb,nBins);
posteriorProbUtility(transitionMatrix(2,:,:),'from c',hpdProb,nBins);
posteriorProbUtility(transitionMatrix(3,:,:),'from d',hpdProb,nBins);

disp('Posterior mean early steady-state state probabilities');
disp(num2str(steadyStateProbabilities(mean(transitionMatrixEarly,3))));
disp('Posterior mean steady-state state probabilities');
disp(num2str(steadyStateProbabilities(mean(transitionMatrix,3))));

steadyStateHPD(transitionMatrixEarly,'early subperiod',hpdProb,nBins);
steadyStateHPD(transitionMatrix,'later subperiod',hpdProb,nBins);

disp('h is heteroskedastic');
disp('we set m to 1/var(gap) and v = 1/8 so std dev=4m');
disp(['mean prior m = ',num2str(mean(m))]);
disp(['mean prior v = ',num2str(mean(v))]);
disp(['mean posterior mean h = ',num2str(mean(mean(h,2)))]);
disp(['median posterior mean h = ',num2str(median(mean(h,2)))]);
disp(['mean posterior std h = ',num2str(mean(std(h,'',2)))]);
[lowerBound,upperBound] = hpd(mean(h,2),hpdProb,nBins);
disp([num2str(hpdProb),'% HPD for mean(h,2) is [',...
    num2str(lowerBound),',',num2str(upperBound),']']);


disp(['ROW burnIns = ',num2str(burnIns),' Gibbs keep = ',...
    num2str(gibbsKeep)]);

% now summarize state in each period
firstYearInSample = USAyear(1)+1; %remember first dif
lastYearInSample = USAyear(end);
nYearsInSample = lastYearInSample - firstYearInSample + 1;
meanStateCountryByYear = nan(nCountries,nYearsInSample,3);
meanStateCountryByYearnoSplit = nan(nCountries,nYearsInSample,3);
for iCountry = 1:nCountries
    firstOb = firstObForCountry(iCountry);
    lastOb = lastObForCountry(iCountry);
    firstYear = firstYearForCountry(iCountry) - (firstYearInSample-1) + 1;
    lastYear = lastYearInSample - (firstYearInSample-1);
    for iProb = 1:3
        meanStateCountryByYear(iCountry, firstYear:lastYear,iProb) = ...
            mean(S(firstOb+1:lastOb,:)==iProb,2);
        meanStateCountryByYearnoSplit(iCountry, firstYear:lastYear,iProb) = ...
            mean(SnoSplit(firstOb+1:lastOb,:)==iProb,2);
    end
end


meanStateByYear = zeros(nYearsInSample,3);
weightedMeanStateByYear = zeros(nYearsInSample,3);
meanStateByYearnoSplit = zeros(nYearsInSample,3);
weightedMeanStateByYearnoSplit = zeros(nYearsInSample,3);
for iYear = 1:nYearsInSample
    whichCountriesValid = ~isnan(meanStateCountryByYear(:,iYear,1));
    for iProb = 1:3
        meanStateByYear(iYear,iProb) =...
            mean(meanStateCountryByYear(whichCountriesValid,iYear,iProb));
        weightedMeanStateByYear(iYear,iProb) =...
            weightedMean(meanStateCountryByYear(whichCountriesValid,...
            iYear, iProb),...
            pop(lastObForCountry(whichCountriesValid)));
        meanStateByYearnoSplit(iYear,iProb) =...
            mean(meanStateCountryByYearnoSplit...
            (whichCountriesValid,iYear,iProb));
        weightedMeanStateByYearnoSplit(iYear,iProb) =...
            weightedMean(meanStateCountryByYearnoSplit...
            (whichCountriesValid, iYear, iProb),...
            pop(lastObForCountry(whichCountriesValid)));
    end
end
% use these for nicer plotting
yearsToPlot = 1951:2014;
obBeforeSplit = find(yearsToPlot == splitYear) - 1;
weightedMeanStateByYearToPlot = weightedMeanStateByYear;
weightedMeanStateByYearToPlot(obBeforeSplit,:) = nan;

terminalMeans = nan(nCountries,3);
for iCountry= 1:nCountries
    for iProb = 1:3
        terminalMeans(iCountry,iProb)...
            = meanStateCountryByYear(iCountry,end,iProb);
    end
end
for iProb = 1:3
    figure;
    plot(firstYearInSample:lastYearInSample,...
        weightedMeanStateByYear(:,iProb),...
        firstYearInSample:lastYearInSample,meanStateByYear(:,iProb),...
        'linewidth',1.5);
    ylabel('Means of country means');
    xlabel('Year');
    title(['Mean probability of state ',num2str(iProb)]);
    legend('Population weighted','Unweighted','location','NorthWest');
    print(['Mean State = ',num2str(iProb)],'-depsc2','-tiff');
    print(['Mean State = ',num2str(iProb)],'-dpng');
    
    disp(['Weighted probability of state ',num2str(iProb)]);
    table([firstYearInSample:lastYearInSample]',...
        weightedMeanStateByYear(:,iProb))
    
    % now get table of terminal probabilities

    [popsInOrder,sizeOrderOfCountries] = sort(pop(lastObForCountry),...
        'descend');
    disp(['for state = ',num2str(iProb)]);
    t = table(countryCodesCell(sizeOrderOfCountries),...
        terminalMeans(sizeOrderOfCountries,iProb),...
        popsInOrder,'VariableNames',{'countryCode','meanState','population'})
end

% plot three states together
figure;
p1 = plot(firstYearInSample:lastYearInSample,...
    weightedMeanStateByYearToPlot(:,1),...
    firstYearInSample:lastYearInSample,...
    weightedMeanStateByYearToPlot(:,2),...
    firstYearInSample:lastYearInSample,...
    weightedMeanStateByYearToPlot(:,3),...
    'linewidth',1.5);
hold on
ax = gca;
ax.ColorOrderIndex = 2;
%{
plot(firstYearInSample:lastYearInSample,...
    weightedMeanStateByYearnoSplit(:,1),'--',...
    firstYearInSample:lastYearInSample,...
    weightedMeanStateByYearnoSplit(:,2),'--',...
    firstYearInSample:lastYearInSample,...
    weightedMeanStateByYearnoSplit(:,3),'--',...
    'linewidth',1.5);
%}
plot(firstYearInSample:lastYearInSample,...
    weightedMeanStateByYearnoSplit(:,2),'--',...
    'linewidth',1.5);
ylabel('Means of country means');
xlabel('Year');
title({'Population weighted mean probability of state',...
    'subperiods solid line, whole period dashed line'});
legend(p1,'n','c','d','location','West');
print('Mean States','-depsc2','-tiff');
print('Mean States','-dpng');

% plot convergent state only
figure;
plot(firstYearInSample:lastYearInSample,...
    weightedMeanStateByYearToPlot(:,2), 'linewidth',1.5);
hold on

plot(firstYearInSample:lastYearInSample,...
    weightedMeanStateByYearnoSplit(:,2),'--',...
    'linewidth',1.5);
ylabel('Means of country means');
xlabel('Year');
title('Population weighted mean probability of convergent state');
legend('subperiods','whole sample','location','West');
print('Mean Convergent State Split','-depsc2','-tiff');
print('Mean Convergent State Split','-dpng');

% plot gap and probability
figure;
yyaxis left;
plot(yearsToPlot,weightedGap(2:end),...
    'linewidth',1.5);
ylabel('Output per capita ratio');
xlabel('Year');

title({'Output per capita relative to USA and',...
    'population-weighted average probability of convergent state'});

hold on;
plotUtility(chinaCountryNumber,countryNames,firstObForCountry,...
    lastObForCountry,S,yearOb,gap);
plotUtility(indiaCountryNumber,countryNames,firstObForCountry,...
    lastObForCountry,S,yearOb,gap);

yyaxis right;
%leave out 1989 so as to not get an ugly line

plot(yearsToPlot,weightedMeanStateByYearToPlot(:,2),'linewidth',1.5);

ylabel('Estimated probability of convergent state');
legend('World','China','India',...
    'Prob of conv. state', 'location','northwest');
hold off;
print -dpng 'IncomeHistory';
print -depsc2 'IncomeHistory';

% now same thing for seminar leaving out India and China
figure;
yyaxis left;
plot(yearsToPlot,weightedGap(2:end),...
    'linewidth',1.5);
ylabel('Output per capita ratio');

title({'Output per capita relative to USA and',...
    'population-weighted average probability of convergent state'});
xlabel('Year');

hold on;
yyaxis right;
%leave out 1989 so as to not get an ugly line

plot(yearsToPlot,weightedMeanStateByYearToPlot(:,2),'linewidth',1.5);

ylabel('Estimated probability of convergent state');
legend('World output per capita','Prob of conv. state', 'location','northwest');
hold off;
print -dpng 'IncomeHistoryWorld';
print -depsc2 'IncomeHistoryWorld';


disp(['Last year gap = ',num2str(weightedGap(end))]);
% What would growth rate be if there is one third chance of
% gap falling by 7/8
gapAtFinal = (2/3)*weightedGap(end) + (1/3)*(1-weightedGap(end)*(7/8));
disp(['Predicted final gap = ',num2str(gapAtFinal)]);
requiredGrowth = (gapAtFinal/weightedGap(end)).^(1/100) - 1 ;
disp(['Required growth over 100 years = ',num2str(requiredGrowth)]);

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

function plotUtility(countryNumber,countryNames,firstObForCountry,...
    lastObForCountry,S,yearOb,gap)
    countryName = countryNames{countryNumber};

firstOb = firstObForCountry(countryNumber);
lastOb = lastObForCountry(countryNumber);
SCountry = S(firstOb+1:lastOb,:);

%yyaxis left;
plot(yearOb(firstOb+1:lastOb),exp(gap(firstOb+1:lastOb)),'linewidth',1.5);
%{
ylabel('ratio');
xlabel('year');
yyaxis right;
plot(yearOb(firstOb+1:lastOb),mean(SCountry==2,2));

title([countryName]);
hold off;
%}
end

function t = gapHalfLife(rho,initialGap)
%{
model is
log G_t = rho^t*log(G_0)
find t that sets
G_t = 1 - .5(1-G_0)

log(1 - .5(1-G_0)) = rho^t*log(G_0)

%}
finalRatio = log(1-(1-initialGap)/2)/log(initialGap);
t = log(finalRatio)/log(rho);

end

function rhoUtility(rho,rhoName,cLower,cUpper,hpdProb,nBins,initialGap)
disp(['Posterior median ',rhoName,' = ',num2str(median(rho))]);
meanrho = mean(rho);
disp(['Posterior mean ',rhoName,' = ',num2str(meanrho)]);
stdrho = std(rho);
disp(['Posterior std ',rhoName,' = ',num2str(stdrho)]);
disp(['Posterior mean ',rhoName,' is ',...
    num2str(abs(cLower-meanrho)/stdrho),...
    ' stds above lower limit and ',...
    num2str(abs(cUpper-meanrho)/stdrho),' stds below upper limit']);
[lowerBound,upperBound] = hpd(rho,hpdProb,nBins);
disp([num2str(hpdProb),'% HPD for ',rhoName,' is [',...
    num2str(lowerBound),',',num2str(upperBound),']']);
figure;
hrho = histogram(rho,nBins);
title(['Posterior ',rhoName]);
disp(['Fraction of ',rhoName,' posterior at upper bound = ',...
    num2str(mean(rho==cUpper))]);
disp(['Fraction of ',rhoName,' posterior in upper bin = ',...
    num2str(hrho.Values(end)/sum(hrho.Values)),' last edge = ',...
    num2str(hrho.BinEdges(end),12)]);
disp(['Fraction of ',rhoName,' posterior above .975 = ',...
    num2str(mean(rho>=.975))]);
disp(['Fraction of ',rhoName,' posterior above .975 = ',...
    num2str(mean(rho>=.975))]);

halfLifeMean = gapHalfLife(meanrho,initialGap);
disp(['Time to close half of initial gap of ',num2str(initialGap),...
    ' is ',num2str(halfLifeMean),' at posterior mean']);
end

function priorUtility(alpha,alphaName)
disp(['Prior ',alphaName])
disp(num2str(alpha'));
[ meanFrom,stdFrom ] = dirichletMeanStd(alpha);
disp(['mean = ',num2str(meanFrom')]);
disp(['std = ',num2str(stdFrom')]);
end

function posteriorProbUtility(transitionMatrixRow,fromName,hpdProb,nBins)
disp(['Row ',fromName]);
disp(['Posterior median = ',num2str(median(transitionMatrixRow,3))]);
meanP = mean(transitionMatrixRow,3);
disp(['Posterior mean = ',num2str(meanP)]);
disp(['Posterior std = ',num2str(std(transitionMatrixRow,0,3))]);

s = size(transitionMatrixRow);
for iProb = 1:s(2)
    [lowerBound,upperBound] = hpd(squeeze(transitionMatrixRow(:,iProb,:))...
        ,hpdProb,nBins);
    disp([num2str(hpdProb),'% HPD for probability ',num2str(iProb),' is [',...
        num2str(lowerBound),',',num2str(upperBound),']']);
end
end

function steadyStateHPD(transitionMatrix,fromName,hpdProb,nBins)
gibbsDraws = size(transitionMatrix,3);
steadyStateProbs = nan(3,gibbsDraws);
for iDraw = 1:gibbsDraws
    steadyStateProbs(:,iDraw) = ...
        steadyStateProbabilities(squeeze(transitionMatrix(:,:,iDraw)));
end
for iProb = 1:3
    [lowerBound,upperBound] = hpd(steadyStateProbs(iProb,:),...
        hpdProb,nBins);
    disp([num2str(hpdProb),'% HPD for ',fromName,' steady-state ',...
        num2str(iProb),...
        ' is [',num2str(lowerBound),',',num2str(upperBound),']']);
end
end