function projectTypicalGap3States
%{
Project out the typical gap starting in either a converged or
nonconverged state


Dick Startz
July 2016
modified September 2016
modified April 2017
modified May 2017
%}
tic;
reset(RandStream.getGlobalStream);

fileInUSA = 'bayesianConvergenceEstimates3StatesSplit';
load(fileInUSA);

nYears = 100;
yearsToPlot = (lastYear+1):(lastYear+nYears);

medianRho_c = median(rho_c,2);
medianRho_d = median(rho_d,2);
medianTransitionMatrix = median(transitionMatrix,3);
disp('Median transition matrix');
disp(num2str(medianTransitionMatrix));

nDraws = 10000;
weightedGDPPerCapitaLast = weightedMean(yLast,popLast);
lastGap = weightedGDPPerCapitaLast/yStarLast;
lastLogGap = log(lastGap);

simulatedLogGapsFromConvergent = nan(nYears,nDraws);
simulatedLogGapsFromNonConvergent = nan(nYears,nDraws);
SDrawNonCon = nan(nYears,nDraws);
for iDraw = 1:nDraws
    SDrawCon = myRandomState(medianTransitionMatrix(2,:));
    switch SDrawCon
        case 1
            simulatedLogGapsFromConvergent(1,iDraw) = lastLogGap;
        case 2
            simulatedLogGapsFromConvergent(1,iDraw)...
                = medianRho_c*lastLogGap;
        case 3
            simulatedLogGapsFromConvergent(1,iDraw)...
                = medianRho_d*lastLogGap;
    end
    
    for iYear = 2:nYears
        SDrawCon = myRandomState(medianTransitionMatrix(SDrawCon,:));
        switch SDrawCon
            case 1
                simulatedLogGapsFromConvergent(iYear,iDraw)...
                    = simulatedLogGapsFromConvergent(iYear-1,iDraw);
            case 2
                simulatedLogGapsFromConvergent(iYear,iDraw)...
                    = medianRho_c*simulatedLogGapsFromConvergent(iYear-1,iDraw);
            case 3
                simulatedLogGapsFromConvergent(iYear,iDraw)...
                    = medianRho_d*simulatedLogGapsFromConvergent(iYear-1,iDraw);
        end
    end
    
    SDrawNonCon(1,iDraw) = myRandomState(medianTransitionMatrix(1,:));
    switch SDrawNonCon(1,iDraw)
        case 1
            simulatedLogGapsFromNonConvergent(1,iDraw) = lastLogGap;
        case 2
            simulatedLogGapsFromNonConvergent(1,iDraw)...
                = medianRho_c*lastLogGap;
        case 3
            simulatedLogGapsFromNonConvergent(1,iDraw)...
                = medianRho_d*lastLogGap;
    end
    
    for iYear = 2:nYears
        SDrawNonCon(iYear,iDraw) = myRandomState(medianTransitionMatrix(SDrawNonCon(iYear-1,iDraw),:));
        switch SDrawNonCon(iYear,iDraw)
            case 1
                simulatedLogGapsFromNonConvergent(iYear,iDraw)...
                    = simulatedLogGapsFromNonConvergent(iYear-1,iDraw);
            case 2
                simulatedLogGapsFromNonConvergent(iYear,iDraw)...
                    = medianRho_c*simulatedLogGapsFromNonConvergent(iYear-1,iDraw);
            case 3
                simulatedLogGapsFromNonConvergent(iYear,iDraw)...
                    = medianRho_d*simulatedLogGapsFromNonConvergent(iYear-1,iDraw);
        end
    end
end
figure;
plot(yearsToPlot,mean(SDrawNonCon==2,2),'linewidth',1.5);
title('SDrawNonCon');

figure;

plot(yearsToPlot,mean(exp(simulatedLogGapsFromConvergent),2),...
    'linewidth',1.5);
hold on;
plot(yearsToPlot,mean(exp(simulatedLogGapsFromNonConvergent),2),...
    'linewidth',1.5);
legend('Initial convergent','Initial nonconvergent',...
    'location','northwest');
title({'Simulated ratio income-to-frontier'...
    ' at median posterior parameters'});
xlabel('Year');
hold off;
print('simulated Convergence 3-state','-depsc2','-tiff');
print('simulated Convergence 3-state','-dpng');

% now just look at probability of states
probCAtYear = nan(3,nYears);
probCAtYear(:,1) = medianTransitionMatrix*[0;1;0];
for iYear = 2:nYears
    probCAtYear(:,iYear) = medianTransitionMatrix*probCAtYear(:,iYear-1);
end

probNAtYear = nan(3,nYears);
probNAtYear(:,1) = medianTransitionMatrix*[1;0;0];
for iYear = 2:nYears
    probNAtYear(:,iYear) = medianTransitionMatrix*probNAtYear(:,iYear-1);
end

figure;
plot(yearsToPlot,probCAtYear(2,:),yearsToPlot,probNAtYear(2,:),...
    'linewidth',1.5);
legend('start converge','start nonconverge');
title('Probability of convergent state starting');
xlabel('year');

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

