%% Estimate Posterior Using Hamiltonian Monte Carlo Sampler
%%
% Consider the multiple linear regression model that predicts the US real
% gross national product (|GNPR|) using a linear combination of total
% employment (|E|) and real wages (|WR|).
%
% $$\texttt{GNPR}_t=\beta_0+\beta_2\texttt{E}_t+\beta_3\texttt{WR}_t+\varepsilon_t.$$
%
% For all $t$, $\varepsilon_t$ is a series of independent Gaussian
% disturbances with a mean of 0 and variance $\sigma^2$. Assume these 
% prior distributions:
%
% * $\beta_j\vert\sigma^2$ is a 3-D _t_ distribution with 10
% degrees of freedom for each component, correlation matrix |C|, 
% location |ct|, and scale |st|.
% * $\sigma^2 \sim IG(r_1,r_2)$, with shape $r_1$ and scale $r_2$.
%
% |bayeslm| treats these assumptions and the data likelihood as if the
% corresponding posterior is analytically intractable.
%%
% Declare a MATLAB(R) function that:
%
% * Accepts values of $\beta$ and $\sigma^2$ together in a column vector,
% and accepts values of the hyperparameters
% * Returns the value of the joint prior distribution,
% $\pi\left(\beta,\sigma^2\right)$, given the values of $\beta$ and
% $\sigma^2$
%
% <include>priorMVTIG.m</include>
%
%%
% Create an anonymous function that operates like |priorMVTIG|, but
% accepts the parameter values only, and holds the hyperparameter values
% fixed at arbitrarily chosen values.
prednames = ["E" "WR"];
p = numel(prednames);
numcoeff = p + 1;

rng(1); % For reproducibility
dof = 10;
V = rand(numcoeff);
Sigma = 0.5*(V + V') + numcoeff*eye(numcoeff);
st = sqrt(diag(Sigma));
C = diag(1./st)*Sigma*diag(1./st);
ct = rand(numcoeff,1);
a = 10*rand;
b = 10*rand;
logPDF = @(params)priorMVTIG(params,ct,st,dof,C,a,b);
%%
% Create a custom joint prior model for the linear regression parameters.
% Specify the number of predictors |p|. Also, specify the function handle
% for |priorMVTIG| and the variable names.
PriorMdl = bayeslm(p,'ModelType','custom','LogPDF',logPDF,...
    'VarNames',prednames);
%%
% |PriorMdl| is a |customblm| Bayesian linear regression model object
% representing the prior distribution of the regression coefficients and
% disturbance variance.
%%
% Load the Nelson-Plosser data set.  Create variables for the response and
% predictor series.
load Data_NelsonPlosser
X = DataTable{:,PriorMdl.VarNames(2:end)};
y = DataTable{:,"GNPR"};  
%%
% Estimate the marginal posterior distributions of $\beta$ and $\sigma^2$
% using the Hamiltonian Monte Carlo (HMC) sampler.  Specify drawing 10,000
% samples and a burn-in period of 1000 draws.
PosteriorMdl = estimate(PriorMdl,X,y,'Sampler','hmc','NumDraws',1e4,...
    'Burnin',1e3);
%%
% |PosteriorMdl| is an |empiricalblm| model object storing the draws from
% the posterior distributions.
%%
% View a trace plot and an ACF plot of the draws from the posterior of
% $\beta_1$ (for example) and the disturbance variance.  Do not plot the
% burn-in period.
figure;
subplot(2,1,1)
plot(PosteriorMdl.BetaDraws(2,1001:end));
title(['Trace Plot ' char(8212) ' \beta_1']);
xlabel('MCMC Draw')
ylabel('Simulation Index')
subplot(2,1,2)
autocorr(PosteriorMdl.BetaDraws(2,1001:end))

figure;
subplot(2,1,1)
plot(PosteriorMdl.Sigma2Draws(1001:end));
title(['Trace Plot ' char(8212) ' Disturbance Variance']);
xlabel('MCMC Draw')
ylabel('Simulation Index')
subplot(2,1,2)
autocorr(PosteriorMdl.Sigma2Draws(1001:end))
%%
% The MCMC sample of the disturbance variance appears to mix well.

%% 
% Copyright 2012 The MathWorks, Inc.