%% Select Variables Using Bayesian Lasso Regression
% Consider the multiple linear regression model that predicts US 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 the prior distributions are:
% For k = 0,...,3, \beta_k\vert\sigma^2 has a Laplace distribution with a mean of 0 and a scale of \sigma^2 /\lambda, where \lambda is the shrinkage parameter. The coefficients are conditionally independent.
% \sigma^2 \sim IG(A,B). A and B are the shape and scale, respectively, of an inverse gamma distribution.
% Create a prior model for Bayesian lasso regression. Specify the number of predictors, the prior model type, and variable names. Specify these shrinkages:
% 0.01 for the intercept
% 10 for IPI and WR
% 1e5 for E because it has a scale that is several orders of magnitude larger than the other variables
% The order of the shrinkages follows the order of the specified variable names, but the first element is the shrinkage of the intercept.

p = 3;
PriorMdl = bayeslm(p,'ModelType','lasso','Lambda',[0.01; 10; 1e5; 10],...
    'VarNames',["IPI" "E" "WR"]);

% PriorMdl is a lassoblm 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"};

%
% Perform Bayesian lasso regression by passing the prior model and data to estimate, that is, by estimating the posterior distribution of \beta and \sigma^2. Bayesian lasso regression uses Markov chain Monte Carlo (MCMC) to sample from the posterior. For reproducibility, set a random seed.

rng(1);
PosteriorMdl = estimate(PriorMdl,X,y);

% PosteriorMdl is an empiricalblm model object that stores draws from the posterior distributions of \beta and \sigma^2 given the data. estimate displays a summary of the marginal posterior distributions in the MATLAB® command line. Rows of the summary correspond to regression coefficients and the disturbance variance, and columns correspond 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 IPI is in [4.157, 4.799] 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.427.
% Plot the posterior distributions.

plot(PosteriorMdl)

% Given the shrinkages, the distribution of E is fairly dense around 0. Therefore, E might not be an important predictor.
% By default, estimate draws and discards a burn-in sample of size 5000. However, a good practice is to inspect a trace plot of the draws for adequate mixing and lack of transience. Plot a trace plot of the draws for each parameter. You can access the draws that compose the distribution (the properties BetaDraws and Sigma2Draws) using dot notation.

figure;
for j = 1:(p + 1)
    subplot(2,2,j);
    plot(PosteriorMdl.BetaDraws(j,:));
    title(sprintf('%s',PosteriorMdl.VarNames{j}));
end
figure;
plot(PosteriorMdl.Sigma2Draws);
title('Sigma2');

% The trace plots indicate that the draws seem to mix well. The plots show no detectable transience or serial correlation, and the draws do not jump between states.
% Copyright 2018 The MathWorks, Inc.
%
%