function runGibbsFrontier
%{
do Bayesian estimate for USA random walk

Dick Startz
July 2016
modified September 2016
%}

tic;
reset(RandStream.getGlobalStream);
load('PWT90Data.mat');

obs1970 = find(USAyear==1970);

yStar = yStar; %fool Matlab
dY = log(yStar(obs1970:end))...
    - log(yStar(obs1970-1:end-1));

nobs = length(yStar(obs1970:end));
meanUSAGrowth = mean(dY);
stdUSAGrowth = std(dY);
stdUSAGrowthMean = stdUSAGrowth/sqrt(nobs);

disp(['mean US growth since 1970 = ',num2str(meanUSAGrowth)]);
disp(['std = ',num2str(stdUSAGrowth)]);
disp(['std of mean = ', num2str(stdUSAGrowthMean)]);

burnInsUSA =    10000;
gibbsKeepUSA = 100000;
disp(['keep ',num2str(gibbsKeepUSA),' Gibbs draws after ',...
    num2str(burnInsUSA),' burn ins']);

mu0 = 0;
V0 = 4*stdUSAGrowth^2;
disp(['Prior mu0 = ',num2str(mu0)]);
disp(['Prior V0 = ',num2str(V0), ' prior sqrt(V0) = ',num2str(sqrt(V0))]);

nu0 = 1;
sSqr0 = stdUSAGrowth^2;
disp(['Prior nu0 = ',num2str(nu0)]);
disp(['Prior sSqr0 = ',num2str(sSqr0)]);

[mu, hUSA] = gibbsFrontier(dY,mu0,V0,nu0,sSqr0,burnInsUSA+gibbsKeepUSA);

mu = mu(burnInsUSA+1:end);
hUSA = hUSA(burnInsUSA+1:end);

disp(['Posterior mean mu = ',num2str(mean(mu))]);
disp(['Posterior std  mu = ',num2str(std(mu))]);

disp(['Posterior mean hUSA = ',num2str(mean(hUSA))]);
disp(['Posterior std  hUSA = ',num2str(std(hUSA))]);

histogram(mu);
title(['Posterior \mu, sample mean = ',num2str(meanUSAGrowth)]);
print('gibbsFrontierMu','-dpng');

figure;
histogram(hUSA);
title(['Posterior hUSA, sample inverse variance = ',...
    num2str(sqrt(stdUSAGrowth))]);
print('gibbsFrontierh','-dpng');

save gibbsFrontier.mat mu hUSA burnInsUSA gibbsKeepUSA mu0 nu0 V0 sSqr0;


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

