function quickEstimationCheck
%{
Do some quick summary figures of the Bayesian estimation

Dick Startz
September 2016
modified March 2017 to change in filename and state definition
modified April 2017 to add China plot
modified May 2017 for split
%}

load bayesianConvergenceEstimates3StatesSplit.mat;
[ meanFromN,stdFromN ] = dirichletMeanStd(alpha_n);
[ meanFromC,stdFromC ] = dirichletMeanStd(alpha_c);
[ meanFromD,stdFromD ] = dirichletMeanStd(alpha_d);
meanPriorTransitionMatrix = [meanFromN';meanFromC';meanFromD'];
stdPriorTransitionMatrix = [stdFromN';stdFromC';stdFromD'];

disp('Mean prior transition matrix');
disp(num2str(meanPriorTransitionMatrix));

disp('Std prior transition matrix');
disp(num2str(stdPriorTransitionMatrix));

disp('prior mean steady-state state probabilities');
disp(num2str(steadyStateProbabilities(meanPriorTransitionMatrix)));

disp('Mean transition matrix');
disp(num2str(mean(transitionMatrix,3)));
disp('lower 5% HPD transition matrix');
disp(num2str(prctile(transitionMatrix,5,3)));
disp('upper 95% HPD transition matrix');
disp(num2str(prctile(transitionMatrix,95,3)));

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



%{This doesn't work because probabilities do not add up
% now compute all steady states and find median
steadyStateProb = zeros(3,gibbsKeep);
for iDraw = 1:gibbsKeep
    steadyStateProb(:,iDraw) ...
        = steadyStateProbabilities(transitionMatrix(:,:,iDraw));
end
disp('Posterior median simulated steady-state state probabilities');
disp(num2str(median(steadyStateProb,2)));
%}

figure;
subplot(1,3,1);
histogram(rho_n);
title(['mean \rho_n = ',num2str(mean(rho_n))]);
subplot(1,3,2);
histogram(rho_c);
title(['mean \rho_c = ',num2str(mean(rho_c))]);
subplot(1,3,3);
histogram(rho_d);
title(['mean \rho_d = ',num2str(mean(rho_d))]);

for iCountry = 1:nCountries
    if strcmp(countryCodesCell{iCountry},'KOR')
        koreaCountryNumber = iCountry;
    elseif strcmp(countryCodesCell{iCountry},'ESP')
        spainCountryNumber = iCountry;
    elseif strcmp(countryCodesCell{iCountry},'CHN')
        chinaCountryNumber = iCountry;
    elseif strcmp(countryCodesCell{iCountry},'IND')
        indiaCountryNumber = iCountry;
    elseif strcmp(countryCodesCell{iCountry},'CAN')
        canadaCountryNumber = iCountry;
    end
end

% now plot Korea
plotUtility(koreaCountryNumber,countryNames,h,firstObForCountry,...
    lastObForCountry,S,yearOb,gap,popWeightWorldLast);
print -dpng 'growthKorea';
% now plot Spain
plotUtility(spainCountryNumber,countryNames,h,firstObForCountry,...
    lastObForCountry,S,yearOb,gap,popWeightWorldLast);
print -dpng 'growthSpain';

% now plot China
plotUtility(chinaCountryNumber,countryNames,h,firstObForCountry,...
    lastObForCountry,S,yearOb,gap,popWeightWorldLast);
print -dpng 'growthChina';
% now plot India
plotUtility(indiaCountryNumber,countryNames,h,firstObForCountry,...
    lastObForCountry,S,yearOb,gap,popWeightWorldLast);
print -dpng 'growthIndia';
% now plot Canada
plotUtility(canadaCountryNumber,countryNames,h,firstObForCountry,...
    lastObForCountry,S,yearOb,gap,popWeightWorldLast);
print -dpng 'growthCanada';

% Now compute US growth rates since 1990

disp('US growth rates to 2014')
for iYear = 64:-1:1
    g = (yStarLast/yStar(65-iYear)).^(1/iYear) - 1;
    disp([num2str(2014-iYear),'    ', num2str(g)]);
end
end
function plotUtility(countryNumber,countryNames,h,firstObForCountry,...
    lastObForCountry,S,yearOb,gap,popWeightWorldLast)
countryName = countryNames{countryNumber};
figure;
histogram(1./h(countryNumber,:));
title(['mean \sigma^2 = ',num2str(mean(1./h(countryNumber,:))),...
    ' ',countryName]);

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

year1990 = lastOb - 24; % bad hardcoding, bad.

disp(['1990 ratio for ',countryName,' = ',...
    num2str(exp(gap(year1990)))]);
disp(['Final ratio for ',countryName,' = ',...
    num2str(exp(gap(lastOb)))]);
disp(['Final population proportion = ',...
    num2str(popWeightWorldLast(countryNumber))]);
yyaxis left;
plot(yearOb(firstOb+1:lastOb),exp(gap(firstOb+1:lastOb)),'linewidth',1.5);
ylabel('ratio');
xlabel('year');

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

year1989 = find(yearOb(firstOb+1:lastOb) == 1989);
yearsToPlot = yearOb(firstOb+1:lastOb);
yearsToPlot(year1989) = nan;

plot(yearsToPlot,mean(SCountry==2,2),'linewidth',1.5);
hold on;
plot(yearsToPlot,mean(SCountry==3,2),'linewidth',1.5);
switch countryName
    case 'China'
        bestLocation = 'Southeast'
    case 'India'
        bestLocation = 'Northwest'
    otherwise
        bestLocation = 'best'
end
legend('relative income','prob S=c','prob S=d','location',bestLocation);
title(['Income relative to USA and estimated states',' ',countryName]);
hold off;
end

