%% Specify Gradient for HMC Sampler
% This example shows how to set up a Bayesian linear regression model for
% efficient posterior sampling using the Hamiltonian Monte Carlo (HMC)
% sampler. The prior distribution of the coefficients is a multivariate
% _t_, and the disturbance variance has an inverse gamma prior.
%%
% Consider the multiple linear regression model that predicts U.S. real
% gross national product (|GNPR|) using a linear combination of industrial
% production index (|IPI|), total employment (|E|), and real wages (|WR|). 
%
% $$\texttt{GNPR}_t=\beta_0+\beta_1\texttt{IPI}_t+\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 4-D _t_ distribution with 30
% degrees of freedom for each component, correlation matrix |C|, 
% location |ct|, and scale |st|.
%
% * $\sigma^2 \sim IG(10r_1,10r_2)$, with shape $10r_1$ and scale $10r_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 log of the joint prior distribution,
% $\pi\left(\beta,\sigma^2\right)$, given the values of $\beta$ and
% $\sigma^2$, and returns the gradient of the log prior density with respect to
% $\beta$ and $\gamma$.
%
% <include>priorMVTIGHMC.m</include>
%
%%
% Create an anonymous function that operates like |priorMVTIGHMC|, but accepts
% the parameter values only and holds the hyperparameter values fixed at
% arbitrarily chosen values.
rng(1); % For reproducibility
dof = 30;
V = rand(4,4);
Sigma = V'*V;
st = sqrt(diag(Sigma));
C = Sigma./(st*st');
ct = -10*rand(4,1);
a = 10*rand;
b = 10*rand;
logPDF = @(params)priorMVTIGHMC(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 |priorMVTIGHMC| and the variable names.
p = 3;
PriorMdl = bayeslm(p,'ModelType','custom','LogPDF',logPDF,...
    'VarNames',["IPI" "E" "WR"]);
%%
% |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. For numerical stability, convert the series to returns.
load Data_NelsonPlosser
X = price2ret(DataTable{:,PriorMdl.VarNames(2:end)});
y = price2ret(DataTable{:,'GNPR'});  
%%
% Estimate the marginal posterior distributions of $\beta$ and $\sigma^2$
% using the HMC sampler.  To monitor the progress of the sampler, specify a
% verbosity level of |2|.  For numerical stability, specify
% reparameterization of the disturbance variance.
options = sampleroptions('Sampler','hmc','VerbosityLevel',2)
PosteriorMdl = estimate(PriorMdl,X,y,'Options',options,'Reparameterize',true);
%%
% |PosteriorMdl| is an |empiricalblm| model object storing the draws from
% the posterior distributions.
%%
% View trace plots and autocorrelation function (ACF) plots of the draws
% from the posteriors.
for j = 1:4
    figure;
    subplot(2,1,1)
    plot(PosteriorMdl.BetaDraws(j,:));
    title(sprintf('Trace plot - Beta %d',j));
    xlabel('MCMC draw')
    ylabel('Simulation index')
    subplot(2,1,2)
    autocorr(PosteriorMdl.BetaDraws(j,:))
end

figure;
subplot(2,1,1)
plot(PosteriorMdl.Sigma2Draws);
title('Trace plot - Disturbance Variance');
xlabel('MCMC draw')
ylabel('Simulation index')
subplot(2,1,2)
autocorr(PosteriorMdl.Sigma2Draws)
%%
% The trace plots show adequate mixing and no transient behavior. The
% autocorrelation plots show low correlation among subsequent samples.


%% 
% Copyright 2012 The MathWorks, Inc.