function [transitionMatrix, transitionMatrixEarly,rho_n,rho_c,rho_d,h,S]...
    = allCountrySwitchingSampler3StatesWeightedSplit(gap,nCountries,...
    firstObForCountry, lastObForCountry,originalNobs,...
    alpha_n,alpha_c,alpha_d,cLower_n,cUpper_n,...
    cLower_c,cUpper_c,cLower_d,cUpper_d,m,v,popWeightROWLast,...
    firstYearForCountry,splitYear, gibbsDraws)
%{
Gibbs sampler for ar1 model
    This version fixes rho0 = 1
Input:
    gap               gap
    nCountries      number of countries
    firstObForCountry   first observation for country i
    lastObForCountry    last observation for country i
    originalNobs    length of gap, total number of valid observations
    
    alpha_n       Dirichlet priors
    alpha_c      
    alpha_d      
    cLower_n,cUpper_n uniform prior for rho_n
    cLower_c,cUpper_c uniform prior for rho_c
    cLower_d,cUpper_d uniform prior for rho_d
    m,v             Gamma prior for h, h = 1/sig2 -> Gamma(m,v)
                    E(h) = m and Var(h) = 2*(m^2)/v
    gibbsDraws  number of Gibbs draws to return

Output:
    transitionMatrix    state transition probabilities
    rho_n,rho_c,rho_d   AR1 coefficients
    h                   precision (inverse variance) of error term
    S                   states

Dick Startz
May 2016
modified September 2016
modified March 2017, 3 states, 1 for no converge 2 for converge
and 3 for diverge
modified April 2017 to put in country weights
%}
    
    transitionMatrix = nan(3,3,gibbsDraws);
    transitionMatrixEarly = nan(3,3,gibbsDraws);
    rho_n = nan(1,gibbsDraws);
    rho_c = nan(1,gibbsDraws);
    rho_d = nan(1,gibbsDraws);
    h = nan(nCountries,gibbsDraws);
    
    S = nan(originalNobs,gibbsDraws);   % stack states across countries in
    % location parallel to gap
    
    % do some initialization using prior means
    transitionMatrix(1,:,1) = dirichletMeanStd(alpha_n);
    transitionMatrix(2,:,1) = dirichletMeanStd(alpha_c);
    transitionMatrix(3,:,1) = dirichletMeanStd(alpha_d);

    transitionMatrixEarly(1,:,1) = dirichletMeanStd(alpha_n);
    transitionMatrixEarly(2,:,1) = dirichletMeanStd(alpha_c);
    transitionMatrixEarly(3,:,1) = dirichletMeanStd(alpha_d);
    %rho_n(1) = (cLower_n+cUpper_n)/2;
    rho_c(1) = (cLower_c+cUpper_c)/2;
    rho_d(1) = (cLower_d+cUpper_d)/2;
    h(:,1) = m;
    
    %   fix rho_n = 1
    rho_n(:) = 1;
    tic;
    for iDraw = 2:gibbsDraws
        for iCountry = 1:nCountries
            newSplitOb = firstObservationForYear(iCountry,...
                firstObForCountry,firstYearForCountry,splitYear);
            if isnan(newSplitOb)
                continue;
            end
            % do first subperiod
            lastOb = newSplitOb - 1;
            firstOb = firstObForCountry(iCountry);
            LHS = gap(firstOb+1:lastOb);
            RHS = gap(firstOb:lastOb-1);
            rhoVec = [rho_n(iDraw-1);rho_c(iDraw-1);rho_d(iDraw-1)];
            if ~(isempty(LHS) ||isempty(RHS))
            S(firstOb+1:lastOb,iDraw) = carterKohnAR1Restricted...
                (LHS,RHS, rhoVec, repmat(h(iCountry,iDraw-1),3,1),...
                transitionMatrixEarly(:,:,iDraw-1));
            end
        
            % do second subperiod
            if newSplitOb < firstOb - 1
                error(['second sample error ',num2str(iCountry)]);
                % saefty check
            end
            
            lastOb = lastObForCountry(iCountry);
            firstOb = newSplitOb;
            LHS = gap(firstOb:lastOb);
            RHS = gap(firstOb-1:lastOb-1);
            rhoVec = [rho_n(iDraw-1);rho_c(iDraw-1);rho_d(iDraw-1)];
            
            S(firstOb:lastOb,iDraw) = carterKohnAR1Restricted...
                (LHS,RHS, rhoVec, repmat(h(iCountry,iDraw-1),3,1),...
                transitionMatrix(:,:,iDraw-1));
            
            
            lastOb = lastObForCountry(iCountry); % whole sample for h
            firstOb = firstObForCountry(iCountry);
            LHS = gap(firstOb+1:lastOb);
            RHS = gap(firstOb:lastOb-1);
            
            h(iCountry,iDraw) = drawh(LHS,...
                [(S(firstOb+1:lastOb,iDraw)==1)...
                .*RHS,(S(firstOb+1:lastOb,iDraw)==2).*RHS,...
                (S(firstOb+1:lastOb,iDraw)==3).*RHS],...
                rhoVec,m(iCountry),v);
        end
        
            desiredState = 2;
            rho_c(iDraw) = drawSingleParamUPriorAllCountriesWeighted(gap,...
                firstObForCountry, lastObForCountry,nCountries,...
                S(:,iDraw),desiredState,h(:,iDraw),cLower_c,cUpper_c,...
                popWeightROWLast);
            
            desiredState = 3;
            rho_d(iDraw) = drawSingleParamUPriorAllCountriesWeighted(gap,...
                firstObForCountry, lastObForCountry,nCountries,...
                S(:,iDraw),desiredState,h(:,iDraw),cLower_d,cUpper_d,...
                popWeightROWLast);
        
        
            % do first subperiod
            transitionMatrixEarly(:,:,iDraw) = ...
                drawTransProbAllCountries3StatesWeightedSplit(S(:,iDraw),...
                firstObForCountry, lastObForCountry,nCountries,...
                alpha_n,alpha_c,alpha_d,popWeightROWLast,...
                firstYearForCountry,splitYear,true);
 
            % do second subperiod
            transitionMatrix(:,:,iDraw) = ...
                drawTransProbAllCountries3StatesWeightedSplit(S(:,iDraw),...
                firstObForCountry, lastObForCountry,nCountries,...
                alpha_n,alpha_c,alpha_d,popWeightROWLast,...
                firstYearForCountry,splitYear,false);
            
    end
end

