%% Estimate Marginal Posterior Distributions
% Consider the linear regression model in Create Custom Multivariate t Prior Model for Coefficients.
% Create an anonymous function that operates like priorMVTIG, but accepts the parameter values only and holds the hyperparameter values fixed at their values.

dof = 50;
C = eye(4);
ct = [-25; 4; 0; 3];
st = ones(4,1);
a = 3;
b = 1;
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.

p = 3;
PriorMdl = bayeslm(p,'ModelType','custom','LogPDF',logPDF,...
    'VarNames',["IPI" "E" "WR"])

%
% 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. Specify a width for the slice sampler that is close to the posterior standard deviation of the parameters assuming a diffuse prior model. Reduce serial correlation by specifying a thinning factor of 10, and reduce the effective default number of draws by a factor of 10.

width = [20,0.5,0.01,1,20];
thin = 10;
numDraws = 1e5/thin;
rng(1) % For reproducibility
PosteriorMdl = estimate(PriorMdl,X,y,'Width',width,'Thin',thin,...
    'NumDraws',numDraws);

% PosteriorMdl is an empiricalblm model object storing draws from the posterior distributions of \beta and \sigma^2 given the data.  estimate displays a summary of the marginal posterior distributions to the command window. Rows of the summary correspond to regression coefficients and the disturbance variance, and columns to characteristics of the posterior distribution. The characteristics include:
% CI95, which contains the 95% Bayesian equitailed credible intervals for the parameters. For example, the posterior probability that the regression coefficient of WR is in [1.939, 3.222] is 0.95.
% Positive, which contains the posterior probability that the parameter is greater than 0. For example, the probability that the intercept is greater than 0 is 0.
% estimate derives the posterior characteristics from draws from the posterior distributions, which MATLAB® stores as matrices in the properties BetaDraws and Sigma2Draws.
%
% To monitor mixing and convergence of the MCMC sample, construct trace plots. In the BetaDraws property, draws correspond to columns and parameters to rows.

figure;
for j = 1:4
    subplot(2,2,j)
    plot(PosteriorMdl.BetaDraws(j,:))
    title(sprintf('Trace Plot of %s',PosteriorMdl.VarNames{j}));
end

figure;
plot(PosteriorMdl.Sigma2Draws)
title('Trace Plot of Sigma2');

% The trace plots indicate adequate mixing and convergence, and there are no transient effects to remove.
%
% Copyright 2012 The MathWorks, Inc.