function [mu, h] = gibbsFrontier(deltaY,mu0,V0,nu0,sSqr0,gibbsDraws)
%{

Gibbs sampler for frontier growth assuming Normal Independent Gamma
prior

Not that there is anything special, it's really just estimating the mean
See Koop page 60ff for algorithm

deltaY      data
mu0         mean normal prior
V0          variance normal prior
nu0         degrees of freedom gamma prior
sSqr0       inverse prior mean gamma prior

mu, h       posterior sample

Dick Startz
July 2016
%}

n = length(deltaY);
nuPost = n + nu0;
sumDeltaY = sum(deltaY);

mu = nan(gibbsDraws,1);
h = nan(gibbsDraws,1);

h(1) = 1/sSqr0; % arbitary initialization
mu(1) = mu0;

for iDraw = 2:gibbsDraws
    VPost = 1/((1/V0) + h(iDraw-1)*n);
    betaPost = VPost*((1/V0)*mu0 + h(iDraw-1)*sumDeltaY);
    mu(iDraw) = normrnd(betaPost,sqrt(VPost));
    
    ssr = sum((deltaY - mu(iDraw)).^2);
    sSqrPost = (ssr + nu0*sSqr0)/nuPost;
    
    % Matlab uses different conventions than Koop for Gamma
    h(iDraw) = gamrnd(nuPost/2,(1/sSqrPost)/(nuPost/2));
end

end

