function estimateBayesianConvergence3StatesSplit
%{

Gibbs sampler for Markov-switching convergence model
Dick Startz
September 2016
modified March 2017, states are now 1 for no converge 2 for converge
and 3 for diverge
modified April 17 to change priors
%}
tic;

reset(RandStream.getGlobalStream);
load PWT90Data;

burnIns =    10000;
gibbsKeep = 100000;

splitYear = 1990;

cLower_n = 0.995; % c_n not really used.
cUpper_n = 1;
cLower_c = 0.95;
cUpper_c = 0.98;
cLower_d = 1.02;
cUpper_d = 1.05;

%{
alpha_n = 40*[0.950;0.025;0.025]/10;
alpha_c=[10;2*995;.5]/10;
alpha_d = 40*[0.025;0.025;0.950]/10;
%}
alpha_n = [550;10;10]/14;
alpha_c=[10;550;5]/14;
alpha_d = [10;10;550]/14;


%{
We want to allow for heteroskedasticity
Note that Gamma prior for h, h = 1/sig2 -> Gamma(m,v)
                    E(h) = m and Var(h) = 2*(m^2)/v
(We do an internal translation to Matlab's parameterization)

Let's set m = 1/var(gap) and then set v so that the 
standard deviation is 4 times the mean
sqrt(2*(m^2)/v) =  4m => 2/16

%}
m = nan(nCountries,1);
for iCountry = 1:nCountries
    lastOb = lastObForCountry(iCountry);
    firstOb = firstObForCountry(iCountry);
    m(iCountry) = 1/var(gap(firstOb+1:lastOb));
end
v = 1/8;

[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,burnIns+gibbsKeep);

transitionMatrix = transitionMatrix(:,:,burnIns+1:end);
transitionMatrixEarly = transitionMatrixEarly(:,:,burnIns+1:end);
rho_n = rho_n(burnIns+1:end);
rho_c = rho_c(burnIns+1:end);
rho_d = rho_d(burnIns+1:end);
h = h(:,burnIns+1:end);
S = S(:,burnIns+1:end);

save bayesianConvergenceEstimates3StatesSplit -v7.3;
elapsedTime = toc;
disp(['Executed ',mfilename,' in ',num2str(elapsedTime),' seconds ',...
    datestr(now,'on mmmm-dd-yyyy at HH:MM:SS')]);
end

