%% Bayesian Lasso Regression
% This example shows how to perform variable selection by using Bayesian
% lasso regression. 
%%
% Lasso regression is a linear regression technique that combines
% regularization and variable selection.  Regularization helps prevent
% overfitting by decreasing the magnitude of the regression coefficients.
% The frequentist view of lasso regression differs from that of other
% regularization techniques, such as, ridge regression, because lasso attributes a
% value of exactly 0 to regression coefficients corresponding to predictors
% that are insignificant or redundant.
%%
% Consider this multiple linear regression model: 
%
% $$y = \beta_0 1_n + X\beta + \varepsilon.$$
%
% * $y$ is a vector of _n_ responses.
% * $X$ is a matrix of _p_ corresponding observed predictor variables.
% * $\beta$ is a vector of _p_ regression coefficients.
% * $\beta_0$ is the intercept.
% * $1_n$ is a length _n_ column vector of ones.
% * $\varepsilon$ is a random vector of iid Gaussian disturbances with a mean of 0 and variance of $\sigma^2$.
% 
% The objective function of the frequentist view of lasso regression is
%
% $$f\left( {\beta_0,\beta;y,X} \right) = 0.5\left\| {y - \beta_0 1_n + X\beta } \right\|_2^2 + \lambda \left\| \beta  \right\|_1^2.$$ 
%
% $\lambda$ is the penalty (or shrinkage) term and, unlike other
% parameters, is not fit to the data.  You must specify a value for
% $\lambda$ before estimation, and a best practice is to tune it.  After
% you specify a value, $f$ is minimized with respect to the regression
% coefficients.  The resulting coefficients are the lasso estimates. For
% more details on the frequentist view of lasso regression, see
% <docid:econ_ug.bvmwdog>.
%%
% For this example, consider creating a predictive linear model for the US
% unemployment rate.  You want a model that generalizes well. In other
% words, you want to minimize the model complexity by removing all
% redundant predictors and all predictors that are uncorrelated with the
% unemployment rate. 
%% Load and Preprocess Data
% Load the US macroeconomic data set |Data_USEconModel.mat|. 
load Data_USEconModel
%%
% The data set includes the MATLAB(R) timetable |DataTable|, which contains 14
% variables measured from Q1 1947 through Q1 2009; |UNRATE| is the US
% unemployment rate.  For more details, enter |Description| at the
% command line.
%%
% Plot all series in the same figure, but in separate subplots. 
figure;
for j = 1:size(DataTable,2)
    subplot(4,4,j)
    plot(DataTable.Time,DataTable{:,j});
    title(DataTable.Properties.VariableNames(j));
end
%%
% All series except |FEDFUNDS|, |GS10|, |TB3MS|, and |UNRATE| appear to
% have an exponential trend.
%%
% Apply the log transform to those variables that have an exponential trend.
hasexpotrend = ~ismember(DataTable.Properties.VariableNames,...
    ["FEDFUNDS" "GD10" "TB3MS" "UNRATE"]);
DataTableLog = varfun(@log,DataTable,'InputVariables',...
    DataTable.Properties.VariableNames(hasexpotrend));
DataTableLog = [DataTableLog DataTable(:,DataTable.Properties.VariableNames(~hasexpotrend))];
%%
% |DataTableLog| is a timetable like |DataTable|, but those variables with
% an exponential trend are on the log scale.
%%
% Coefficients that have relatively large magnitudes tend to dominate the
% penalty in the lasso regression objective function.  Therefore, it is
% important that variables have a similar scale when you implement lasso
% regression. Compare the scales of the variables in |DataTableLog| by
% plotting their box plots on the same axis.
figure;
boxplot(DataTableLog.Variables,'Labels',DataTableLog.Properties.VariableNames);
h = gca;
h.XTickLabelRotation = 45;
title('Variable Box Plots');
%%
% The variables have fairly similar scales.
%%
% To determine how well time series models generalize, follow this
% procedure:
%
% # Partition the data into estimation and forecast samples.
% # Fit the models to the estimation sample.
% # Use the fitted models to forecast responses into the forecast horizon.
% # Estimate the forecast mean squared error (FMSE) for each model.
% # Choose the model with the lowest FMSE.
%
%% 
% Create estimation and forecast sample variables for the response and
% predictor data. Specify a forecast horizon of 4 years (16 quarters).
fh = 16;
y = DataTableLog.UNRATE(1:(end - fh));
yF = DataTableLog.UNRATE((end - fh + 1):end);
isresponse = DataTable.Properties.VariableNames == "UNRATE";
X = DataTableLog{1:(end - fh),~isresponse};
XF = DataTableLog{(end - fh + 1):end,~isresponse};
p = size(X,2); % Number of predictors
predictornames = DataTableLog.Properties.VariableNames(~isresponse);
%% Fit Simple Linear Regression Model
% Fit a simple linear regression model to the estimation sample. Specify
% the variable names (the name of the response variable must be the last
% element). Extract the _p_-values. Identify the insignificant coefficients
% by using a 5% level of significance.
SLRMdlFull = fitlm(X,y,'VarNames',DataTableLog.Properties.VariableNames)
slrPValue = SLRMdlFull.Coefficients.pValue;
isInSig = SLRMdlFull.CoefficientNames(slrPValue > 0.05)
%%
% |SLRMdlFull| is a |LinearModel| model object. The _p_-values suggest that, at
% a 5% level of significance, |CPIAUCSL|, |GS10|, |M2SL|, and 
% |FEDFUNDS| are insignificant variables.
%%
% Refit the linear model without the insignificant variables in the
% predictor data. Using the full and reduced models, predict the unemployment
% rate into the forecast horizon, then estimate the FMSEs.
idxreduced = ~ismember(predictornames,isInSig);
XReduced = X(:,idxreduced);
SLRMdlReduced = fitlm(XReduced,y,'VarNames',[predictornames(idxreduced) {'UNRATE'}]);

yFSLRF = predict(SLRMdlFull,XF);
yFSLRR = predict(SLRMdlReduced,XF(:,idxreduced));

fmseSLRF = sqrt(mean((yF - yFSLRF).^2))
fmseSLRR = sqrt(mean((yF - yFSLRR).^2))
%%
% The FMSE of the reduced model is less than the FMSE of the full model,
% which suggests that the reduced model generalizes better.
%% Fit Frequentist Lasso Regression Model
% Fit a lasso regression model to the data.  Use the default regularization
% path. |lasso| standardizes the data by default.  Because the variables
% are on similar scales, do not standardize them. Return
% information about the fits along the regularization path.
[LassoBetaEstimates,FitInfo] = lasso(X,y,'Standardize',false,...
    'PredictorNames',predictornames); 
%%
% By default, |lasso| fits the lasso regression model 100 times by using a
% sequence of values for $\lambda$.  Therefore, |LassoBetaEstimates| is a
% 13-by-100 matrix of regression coefficient estimates.  Rows correspond to
% predictor variables in |X|, and columns correspond to values of $\lambda$.
% |FitInfo| is a structure that includes the values of $\lambda$
% (|FitInfo.Lambda|), the mean squared error for each fit (|FitInfo.MSE|),
% and the estimated intercepts (|FitInfo.Intercept|).
%%
% Compute the FMSE of each model returned by |lasso|.
yFLasso = FitInfo.Intercept + XF*LassoBetaEstimates;
fmseLasso = sqrt(mean((yF - yFLasso).^2,1));
%%
% Plot the magnitude of the regression coefficients with respect to the
% shrinkage value.
hax = lassoPlot(LassoBetaEstimates,FitInfo);
L1Vals = hax.Children.XData;
yyaxis right
h = plot(L1Vals,fmseLasso,'LineWidth',2,'LineStyle','--');
legend(h,'FMSE','Location','SW');
ylabel('FMSE');
title('Frequentist Lasso')
%%
% A model with 7 predictors (df = 7) seems to balance minimal FMSE and
% model complexity well.  Obtain the coefficients corresponding to the
% model containing 7 predictors and yielding minimal FMSE.
fmsebestlasso = min(fmseLasso(FitInfo.DF == 7));
idx = fmseLasso == fmsebestlasso;
bestLasso = [FitInfo.Intercept(idx); LassoBetaEstimates(:,idx)];
table(bestLasso,'RowNames',["Intercept" predictornames])
%%
% The frequentist lasso analysis suggests that the variables |CPIAUCSL|,
% |GCE|, |GDP|, |GPDI|, |PCEC|, and |FEDFUNDS| are either insignificant or
% redundant.
%% Fit Bayesian Lasso Regression Model
% In the Bayesian view of lasso regression, the prior distribution of the
% regression coefficients is Laplace (double exponential), with mean 0 and
% scale $\sigma/\lambda$, where $\lambda$ is the fixed shrinkage parameter
% and $\sigma^2\sim IG(A,B)$. For more details, see
% <docid:econ_ug.object_lassoblm>. 
%%
% As with the frequentist view of lasso regression, if you increase
% $\lambda$, then the sparsity of the resulting model increases
% monotonically. However, unlike frequentist lasso, Bayesian lasso has
% insignificant or redundant coefficient modes that are not exactly 0.  Instead,
% estimates have a posterior distribution and are insignificant or
% redundant if the region around 0 has high density.
%%
% When you implement Bayesian lasso regression in MATLAB(R), be aware of
% several differences between the Statistics and Machine Learning
% Toolbox(TM) function <docid:stats_ug.bs25w54-1 |lasso|> and the
% Econometrics Toolbox(TM) object <docid:econ_ug.object_lassoblm
% |lassoblm|> and its associated functions.
% 
% * |lassoblm| is part of an object framework, whereas |lasso| is a
% function.  The object framework streamlines econometric workflows. 
% * Unlike |lasso|, |lassoblm| does not standardize the predictor data.
% However, you can supply different shrinkage values for each coefficient.
% This feature helps maintain the interpretability of the coefficient
% estimates.
% * |lassoblm| applies one shrinkage value to each coefficient; it
% does not accept a regularization path like |lasso|.
%
%%
% Because the |lassoblm| framework is suited to performing Bayesian lasso
% regression for one shrinkage value per coefficient, you can use a for
% loop to perform lasso over a regularization path.
%%
% Create a Bayesian lasso regression prior model by using |bayeslm|.
% Specify the number of predictor variables and the variable names. Display
% the default shrinkage value for each coefficient stored in the |Lambda|
% property of the model.
PriorMdl = bayeslm(p,'ModelType','lasso','VarNames',predictornames);
table(PriorMdl.Lambda,'RowNames',PriorMdl.VarNames)
%%
% |PriorMdl| is a |lassoblm| model object.  |lassoblm| attributes
% a shrinkage of |1| for each coefficient except the intercept, which
% has a shrinkage of |0.01|.  These default values might not be useful in
% most applications; a best practice is to experiment with many values.
%%
% Consider the regularization path returned by |lasso|.  Inflate the
% shrinkage values by a factor of $n/\sqrt(MSE_j)$, where $MSE_j$ is the MSE
% of lasso run $j$, $j$ = 1 through the number of shrinkage values.
ismissing = any(isnan(X),2);
n = sum(~ismissing); % Effective sample size
lambda = FitInfo.Lambda*n./sqrt(FitInfo.MSE);

%%
% Consider the regularization path returned by |lasso|.  Loop through the
% shrinkage values at each iteration:
%
% # Estimate the posterior distributions of the regression coefficients and
% the disturbance variance, given the shrinkage and data.  Because the
% scales of the predictors are close, attribute the same shrinkage to each
% predictor, and attribute a shrinkage of |0.01| to the intercept. Store
% the posterior means and 95% credible intervals.
% # For the lasso plot, if a 95% credible interval includes 0, then set the
% corresponding posterior mean to zero.
% # Forecast into the forecast horizon.
% # Estimate the FMSE.
%
numlambda = numel(lambda);

% Preallocate
BayesLassoCoefficients = zeros(p+1,numlambda);
BayesLassoCI95 = zeros(p+1,2,numlambda);
fmseBayesLasso = zeros(numlambda,1);
BLCPlot = zeros(p+1,numlambda);

% Estimate and forecast
rng(10); % For reproducibility
for j = 1:numlambda
    PriorMdl.Lambda = lambda(j);
    [EstMdl,Summary] = estimate(PriorMdl,X,y,'Display',false);
    BayesLassoCoefficients(:,j) = Summary.Mean(1:(end - 1));
    BLCPlot(:,j) = Summary.Mean(1:(end - 1));
    BayesLassoCI95(:,:,j) = Summary.CI95(1:(end - 1),:);
    idx = BayesLassoCI95(:,2,j) > 0 & BayesLassoCI95(:,1,j) <= 0;
    BLCPlot(idx,j) = 0;
    yFBayesLasso = forecast(EstMdl,XF);
    fmseBayesLasso(j) = sqrt(mean((yF - yFBayesLasso).^2));
end

%%
% At each iteration, |estimate| runs a Gibbs sampler to sample from the
% full conditionals (see <docid:econ_ug.bvj1dzd>) and returns the
% |empiricalblm| model |EstMdl|, which contains the draws and the
% estimation summary table |Summary|. You can specify parameters for the
% Gibbs sampler, such as the number of draws or thinning scheme. 
%%
% A good practice is to diagnose the posterior sample by inspecting trace
% plots.  However, because 100 distributions were drawn, this example
% proceeds without performing that step.
%%
% Plot the posterior means of the coefficients with respect to the
% normalized L1 penalties (sum of the magnitudes of the coefficients except
% the intercept). On the same plot, but in the right y axis, plot the
% FMSEs.
L1Vals = sum(abs(BLCPlot(2:end,:)),1)/max(sum(abs(BLCPlot(2:end,:)),1));

figure;
plot(L1Vals,BLCPlot(2:end,:))
xlabel('L1');
ylabel('Coefficient Estimates');
yyaxis right
h = plot(L1Vals,fmseBayesLasso,'LineWidth',2,'LineStyle','--');
legend(h,'FMSE','Location','SW');
ylabel('FMSE');
title('Bayesian Lasso')
%%
% The model tends to generalize better as the normalized L1 penalty
% increases past 0.1.  To balance minimal FMSE and model complexity, choose
% the simplest model with FMSE closest to 0.68.
idx = find(fmseBayesLasso > 0.68,1);
fmsebestbayeslasso = fmseBayesLasso(idx);
bestBayesLasso = BayesLassoCoefficients(:,idx);
bestCI95 = BayesLassoCI95(:,:,idx);
table(bestBayesLasso,bestCI95,'RowNames',["Intercept" predictornames])
%%
% One way to determine if a variable is insignificant or redundant is
% by checking whether 0 is in its corresponding coefficient 95% credible
% interval. Using this criterion, |CPIAUCSL|, |M2SL|, |PCEC|, and
% |FEDFUNDS| are insignificant.
%%
% Display all estimates in a table for comparison.
SLRR = zeros(p + 1,1);
SLRR([true idxreduced]) = SLRMdlReduced.Coefficients.Estimate;
table([SLRMdlFull.Coefficients.Estimate; fmseSLRR],...
    [SLRR; fmseSLRR],...
    [bestLasso; fmsebestlasso],...
    [bestBayesLasso; fmsebestbayeslasso],'VariableNames',...
    {'SLRFull' 'SLRReduced' 'Lasso' 'BayesLasso'},...
    'RowNames',[PriorMdl.VarNames; 'MSE'])
%%
% After you chose a model, re-estimate it using the entire data set.
% For example, to create a predictive Bayesian lasso regression model,
% create a prior model and specify the shrinkage yielding the simplest
% model with the minimal FMSE, then estimate it using the entire data set.
bestLambda = lambda(idx);
PriorMdl = bayeslm(p,'ModelType','lasso','Lambda',bestLambda,...
    'VarNames',predictornames);
PosteriorMdl = estimate(PriorMdl,[X; XF],[y; yF]);
%%
% Copyright 2018 The MathWorks, Inc.
