%% Perform Variable Selection Using Default Lasso Shrinkage
% Consider the linear regression model in docid:econ_ug.mw_3c5ee5cc-859b-42d1-af6b-c8e525af7948.
% Create a prior model for performing Bayesian lasso regression. Specify the number of predictors p and the names of the regression coefficients.

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

% PriorMdl stores the shrinkage values for all predictors in its Lambda property. shrinkage(1) is the shrinkage for the intercept, and the elements of shrinkage(2:end) correspond to the coefficients of the predictors in Mdl.VarNames. The default shrinkage for the intercept is 0.01, and the default is 1 for all other coefficients.
%
% Load the Nelson-Plosser data set. Create variables for the response and predictor series. Because lasso is sensitive to variable scales, standardize all variables.

load Data_NelsonPlosser
X = DataTable{:,PriorMdl.VarNames(2:end)};
y = DataTable{:,'GNPR'};

X = (X - nanmean(X))./nanstd(X);
y = (y - nanmean(y))/nanstd(y);

%
% Although this example standardizes variables, you can specify different shrinkage values for each coefficient by setting the Lambda property of PriorMdl to a numeric vector of shrinkage values.
% Implement Bayesian lasso regression by estimating the marginal posterior distributions of \beta and \sigma^2. Because Bayesian lasso regression uses Markov chain Monte Carlo  (MCMC) for estimation, set a random number seed to reproduce the results.

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 at the 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 E  (standardized) is in [-0.110, 0.365] 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.
% 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 plot shows no detectable transience or serial correlation, and the draws do not jump between states.
%
% Plot the posterior distributions of the coefficients and disturbance variance.

figure;
plot(PosteriorMdl)

% E and WR might not be important predictors because 0 is within the region of high density in their posterior distributions.
% Copyright 2018 The MathWorks, Inc.