function generate3StateFutureSimulations
%{
generate 3-state simulated future data from Gibbs samples
Four scenarios
    actual simulations
    actual simulations assuming zero convergence
    simulations with no parameter uncertainty
    simulations with no shocks to frontier

save results

remember that "gap data" is in logs and needs to be
made into levels

Dick Startz
July 2016
modified September 2016
modified April 2017
%}
tic;
reset(RandStream.getGlobalStream);
inputGibbsFileName = 'bayesianConvergenceEstimates3StatesSplit';
inputUSAFileName = 'gibbsFrontier';
outputFileName = 'simulated3StateFuturesFrom2014Split';

load(inputGibbsFileName);
load(inputUSAFileName);

% now let's do a simulation from the posterior
% sample with replacement effectively
% we will save the results
nDraws = 10000;
nYears = 100;

finalGap = nan(nCountries,1);
for iCountry = 1:nCountries
    finalGap(iCountry) = gap(lastObForCountry(iCountry));
end

simulatedCountryLogGaps = nan(nCountries,nYears,nDraws);
simulatedFutureStates = nan(nCountries,nYears,nDraws);

for iDraw = 1:nDraws
    % first get common parameters
    gibbsDraw = randi(gibbsKeep);
    rho_n_Draw = rho_n(gibbsDraw);
    rho_c_Draw = rho_c(gibbsDraw);
    rho_d_Draw = rho_d(gibbsDraw);
    transitionMatrixDraw = squeeze(transitionMatrix(:,:,gibbsDraw));
    
    % now get country specific draws and forecast out
    for iCountry = 1:nCountries
        hDraw = h(iCountry,gibbsDraw);
        lastOb = lastObForCountry(iCountry);
        finalStateDraw = S(lastOb,gibbsDraw);
        
        epsilonDraw = randn/sqrt(hDraw);
        simulatedFutureStates(iCountry,1,iDraw)...
            = myRandomState(transitionMatrixDraw(finalStateDraw,:));
        %{
If gap is positive, then state can't be three
            %}
            switch finalStateDraw
                case 1
                    simulatedCountryLogGaps(iCountry,1,iDraw)...
                        = rho_n_Draw*finalGap(iCountry) + epsilonDraw;
                case 2
                    simulatedCountryLogGaps(iCountry,1,iDraw)...
                        = rho_c_Draw*finalGap(iCountry) + epsilonDraw;
                case 3
                    if finalGap(iCountry) <= 0
                        simulatedCountryLogGaps(iCountry,1,iDraw)...
                            = rho_d_Draw*finalGap(iCountry) + epsilonDraw;
                    else
                        simulatedCountryLogGaps(iCountry,1,iDraw)...
                            = rho_n_Draw*finalGap(iCountry) + epsilonDraw;
                    end
            end
            
            
            for futureYear = 2:nYears
                epsilonDraw = randn/sqrt(hDraw);
                simulatedFutureStates(iCountry,futureYear,iDraw)...
                    = myRandomState(transitionMatrixDraw(...
                    simulatedFutureStates(iCountry,futureYear-1,iDraw),:));
                switch simulatedFutureStates(iCountry,futureYear-1,iDraw)
                    case 1
                        simulatedCountryLogGaps(iCountry,futureYear,iDraw)...
                            = rho_n_Draw*simulatedCountryLogGaps...
                            (iCountry,futureYear-1,iDraw) + epsilonDraw;
                    case 2
                        simulatedCountryLogGaps(iCountry,futureYear,iDraw)...
                            = rho_c_Draw*simulatedCountryLogGaps...
                            (iCountry,futureYear-1,iDraw) + epsilonDraw;
                    case 3
                        if simulatedCountryLogGaps(iCountry,futureYear-1,...
                                iDraw) <= 0
                            simulatedCountryLogGaps(iCountry,futureYear,...
                                iDraw) = rho_d_Draw*simulatedCountryLogGaps...
                                (iCountry,futureYear-1,iDraw) + epsilonDraw;
                        else
                            simulatedCountryLogGaps(iCountry,futureYear,...
                                iDraw) = rho_n_Draw*simulatedCountryLogGaps...
                                (iCountry,futureYear-1,iDraw) + epsilonDraw;
                        end
                end
            end
    end
end

% Now generate forecasts for the frontier, do this one vectorized for no
% particularly good reason
% do it with and without shocks

simulatedFutureLogUSA = nan(nYears,nDraws);
simulatedFutureLogUSANoSchock = nan(nYears,nDraws);
frontierPosteriorDraws = randi(gibbsKeepUSA,nDraws,1);
simulatedMeans = mu(frontierPosteriorDraws);
simulatedStd = sqrt(1./hUSA(frontierPosteriorDraws));


simulatedShocks = normrnd(simulatedMeans',simulatedStd');
simulatedFutureLogUSA(1,:) = log(yStarLast) + simulatedShocks;
simulatedFutureLogUSANoSchock(1,:) = log(yStarLast) + simulatedMeans;

for iYear = 2:nYears
    simulatedShocks = normrnd(simulatedMeans',simulatedStd');
    simulatedFutureLogUSA(iYear,:) = simulatedFutureLogUSA(iYear-1,:) + ...
        simulatedShocks;
    simulatedFutureLogUSANoSchock(iYear,:) = ...
        simulatedFutureLogUSANoSchock(iYear-1,:) + simulatedMeans';
end

%now do the same with no parameter uncertainty
simulatedFutureCertainLogUSA = nan(nYears,nDraws);
simulatedMeans = median(mu);
simulatedStd = median(sqrt(1./hUSA));


simulatedShocks = normrnd(simulatedMeans,simulatedStd,nDraws,1);
simulatedFutureCertainLogUSA(1,:) = log(yStarLast) + simulatedShocks;

for iYear = 2:nYears
    simulatedShocks = normrnd(simulatedMeans,simulatedStd,nDraws,1)';
    simulatedFutureCertainLogUSA(iYear,:)...
        = simulatedFutureCertainLogUSA(iYear-1,:) + simulatedShocks;
end

% now generate forecasts for each country by the frontier plus gap
simulatedFutureLogIncome = nan(nCountries,nYears,nDraws);
for iCountry = 1:nCountries
    simulatedFutureLogIncome(iCountry,:,:)...
        = simulatedFutureLogUSA...
        + squeeze(simulatedCountryLogGaps(iCountry,:,:));
end

% now the same thing, only assuming the gap stays at the last observed gap
simulatedNoConvergenceFutureLogIncome = nan(nCountries,nYears,nDraws);
for iCountry = 1:nCountries
    simulatedNoConvergenceFutureLogIncome(iCountry,:,:)...
        = simulatedFutureLogUSA + finalGap(iCountry);
end

% now do this without parameter uncertainty
% like first simulation but with constant parameters except states
simulatedCountryCertainLogGaps = nan(nCountries,nYears,nDraws);
simulatedFutureCertainStates = nan(nCountries,nYears,nDraws);

% first get common parameters

rho_n_Draw = median(rho_n);
rho_c_Draw = median(rho_c);
rho_d_Draw = median(rho_d);
transitionMatrixMedian = median(transitionMatrix,3);
hDrawForCountries = median(h,2);

for iDraw = 1:nDraws
    gibbsDraw = randi(gibbsKeep);
    
    % now get country specific draws and forecast out
    for iCountry = 1:nCountries
        hDraw = hDrawForCountries(iCountry);
        lastOb = lastObForCountry(iCountry);
        finalStateDraw = S(lastOb,gibbsDraw);
        
        epsilonDraw = randn/sqrt(hDraw);
        simulatedFutureCertainStates(iCountry,1,iDraw)...
            = myRandomState(transitionMatrixMedian(finalStateDraw,:));
        switch finalStateDraw
            case 1
                simulatedCountryCertainLogGaps(iCountry,1,iDraw)...
                    = rho_n_Draw*finalGap(iCountry) + epsilonDraw;
            case 2
                simulatedCountryCertainLogGaps(iCountry,1,iDraw)...
                    = rho_c_Draw*finalGap(iCountry) + epsilonDraw;
            case 3
                if finalGap(iCountry) <= 0
                   simulatedCountryCertainLogGaps(iCountry,1,iDraw)...
                    = rho_d_Draw*finalGap(iCountry) + epsilonDraw;
                else
                    simulatedCountryCertainLogGaps(iCountry,1,iDraw)...
                    = rho_n_Draw*finalGap(iCountry) + epsilonDraw;
                end
        end
        
        
        for futureYear = 2:nYears
            epsilonDraw = randn/sqrt(hDraw);
            simulatedFutureCertainStates(iCountry,futureYear,iDraw)...
                = myRandomState(transitionMatrixMedian(...
                simulatedFutureCertainStates(iCountry,...
                futureYear-1,iDraw),:));
            switch simulatedFutureCertainStates(iCountry,futureYear-1,iDraw)
                case 1
                    simulatedCountryCertainLogGaps(iCountry,futureYear,...
                        iDraw)...
                        = rho_n_Draw*...
                        simulatedCountryCertainLogGaps(iCountry,...
                        futureYear-1,iDraw) + epsilonDraw;
                case 2
                    simulatedCountryCertainLogGaps(iCountry,futureYear,...
                        iDraw)...
                        = rho_c_Draw*...
                        simulatedCountryCertainLogGaps(iCountry,...
                        futureYear-1,iDraw) + epsilonDraw;
                case 3
                    if simulatedCountryCertainLogGaps(iCountry,...
                            futureYear-1,iDraw) <= 0
                        simulatedCountryCertainLogGaps(iCountry,...
                            futureYear,iDraw)...
                        = rho_d_Draw*...
                        simulatedCountryCertainLogGaps(iCountry,...
                        futureYear-1,iDraw) + epsilonDraw;
                    else
                         simulatedCountryCertainLogGaps(iCountry,...
                            futureYear,iDraw)...
                        = rho_n_Draw*...
                        simulatedCountryCertainLogGaps(iCountry,...
                        futureYear-1,iDraw) + epsilonDraw;
                    end
                    
            end
        end
    end
end

simulatedFutureCertainLogIncome = nan(nCountries,nYears,nDraws);
for iCountry = 1:nCountries
    simulatedFutureCertainLogIncome(iCountry,:,:)...
        = simulatedFutureCertainLogUSA...
        + squeeze(simulatedCountryCertainLogGaps(iCountry,:,:));
end

% now make simulations with no shock to frontier
simulatedFutureLogIncomeNoSchocks = nan(nCountries,nYears,nDraws);
for iCountry = 1:nCountries
    simulatedFutureLogIncomeNoSchocks(iCountry,:,:)...
        = simulatedFutureLogUSANoSchock...
        + squeeze(simulatedCountryLogGaps(iCountry,:,:));
end


save(outputFileName,'-v7.3');
elapsedTime = toc;

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

