function generate3StateFutureSimulationsFromPriors
%{
generate 3-state simulated future data from priors
Two scenarios
    actual simulations
    actual simulations assuming zero convergence

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
modified January 2019
%}
tic;
reset(RandStream.getGlobalStream);
inputGibbsFileName = 'bayesianConvergenceEstimates3StatesSplit';
inputUSAFileName = 'gibbsFrontier';
outputFileName = 'simulated3StateFuturesFrom2014SplitFromPriors';

load(inputGibbsFileName); %We need this for some auxiliary stuff
load(inputUSAFileName);

% now let's do a simulation from the posterior
% sample with replacement effectively
% we will save the results
%nDraws = 10000;
nDraws = 1000;
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);

transitionMatrix = nan(3,3);
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));
    %}
    rho_n_Draw = 1;
    rho_c_Draw = unifrnd(cLower_c,cUpper_c);
    rho_d_Draw = unifrnd(cLower_d,cUpper_d);
    transitionMatrix(1,:) = drchrnd(alpha_n');
    transitionMatrix(2,:) = drchrnd(alpha_c');
    transitionMatrix(3,:) = drchrnd(alpha_d');
    SS_Prob = steadyStateProbabilities(transitionMatrix);
    
    % now get country specific draws and forecast out
    for iCountry = 1:nCountries
        lastOb = lastObForCountry(iCountry); % whole sample for h
        firstOb = firstObForCountry(iCountry);
        LHS = gap(firstOb+1:lastOb);
        hDraw = drawh(LHS,ones(length(LHS),1),mean(LHS),m(iCountry),v);
        lastOb = lastObForCountry(iCountry);
        % now draw the final state from the steady-state transition probs
        kludgeVector = [1;2;3];
        finalStateDraw = kludgeVector(logical(mnrnd(1,drchrnd(SS_Prob')')));
        
        epsilonDraw = randn/sqrt(hDraw);
        simulatedFutureStates(iCountry,1,iDraw)...
            = myRandomState(transitionMatrix(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(transitionMatrix(...
                    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);
simulatedMeans = repmat(mu0,nDraws,1);
simulatedStd = repmat(sqrt(V0),nDraws,1);


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

for iYear = 2:nYears
    simulatedShocks = normrnd(simulatedMeans',simulatedStd');
    simulatedFutureLogUSA(iYear,:) = simulatedFutureLogUSA(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




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

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

