%% Bayesian Stochastic Search Variable Selection
% This example shows how to implement stochastic search variable selection (SSVS), a Bayesian variable selection technique for linear regression models.
%% Introduction
% Consider this Bayesian linear regression model.
% y_t = \sum_k \beta_k x_{tk} + \varepsilon_t.
% The regression coefficients \beta_k|\sigma^2 \sim N(\mu_j,\sigma^2V_k).
% k = 0,...,p.
% The disturbances \varepsilon_t \sim N(0,\sigma^2).
% The disturbance variance \sigma^2 \sim IG(A,B), where IG(A,B) is the inverse gamma distribution with shape A and scale B.
% The goal of variable selection is to include only those predictors supported by data in the final regression model. One way to do this is to analyze the 2^p permutations of models, called regimes, where models differ by the coefficients that are included. If p is small, then you can fit all permutations of models to the data, and then compare the models by using performance measures, such as goodness-of-fit (for example, Akaike information criterion) or forecast mean squared error (MSE). However, for even moderate values of p, estimating all permutations of models is inefficient.
% A Bayesian view of variable selection is a coefficient, being excluded from a model, has a degenerate posterior distribution. That is, the excluded coefficient has a Dirac delta distribution, which has its probability mass concentrated on zero. To circumvent the complexity induced by degenerate variates, the prior for a coefficient being excluded is a Gaussian distribution with a mean of 0 and a small variance, for example N(0,0.1^2). Because the prior is concentrated around zero, the posterior must also be concentrated around zero.
% The prior for the coefficient being included can be N(\mu,V), where V is sufficiently away from zero and \mu is usually zero. This framework implies that the prior of each coefficient is a Gaussian mixture model.
% Consider the latent, binary random variables \gamma_k, k = 0,...,p, such that:
% \gamma_k = 1 indicates that \beta_k\sim N(0,\sigma^2V_{1k}) and that \beta_k is included in the model.
% \gamma_k = 0 indicates that \beta_k\sim N(0,\sigma^2V_{2k}) and that \beta_k is excluded from the model.
% \gamma_k\sim {\rm Bernoulli}(g_k).
% The sample space of \gamma_k has a cardinality of 2^{p+1}, and each element is a p+1-D vector of zeros or ones.
% V_{2k} is a small, positive number and V_{1k} > V_{2k}.
% Coefficients \beta_j and \beta_k, j \ne k are independent, a priori.
% One goal of SSVS is to estimate the posterior regime probabilities {g_k}, the estimates that determine whether corresponding coefficients should be included in the model. Given \beta_k, \gamma_k is conditionally independent of the data. Therefore, for k = 0,...,p, this equation represents the full conditional posterior distribution of the probability that variable k is included in the model:
% P(\gamma_k = 1\vert \beta,\sigma^2,\gamma_{\ne k})\propto g_k \phi\left(\beta_k;0,\sigma^2V_{1k}\right),
% where \phi\left(\mu,\sigma^2\right) is the pdf of the Gaussian distribution with scalar mean \mu and variance \sigma^2.
% Econometrics Toolbox™ has two Bayesian linear regression models that specify the prior distributions for SSVS: mixconjugateblm and mixsemiconjugateblm. The framework presented earlier describes the priors of the mixconjugateblm model. The difference between the models is that \beta and \sigma^2 are independent, a priori, for mixsemiconjugateblm models. Therefore, the prior variance of \beta_k is V_{k1} (\gamma_k =1) or V_{k2} (\gamma_k =0).
% After you decide which prior model to use, call bayeslm to create the model and specify hyperparameter values. Supported hyperparameters include:
% Intercept, a logical scalar specifying whether to include an intercept in the model.
% Mu, a (p + 1)-by-2 matrix specifying the prior Gaussian mixture means of \beta. The first column contains the means for the component corresponding to \gamma_k = 1, and the second column contains the means corresponding to \gamma_k = 0. By default, all means are 0, which specifies implementing SSVS.
% V, a (p + 1)-by-2 matrix specifying the prior Gaussian mixture variance factors (or variances) of \beta. Columns correspond to the columns of Mu. By default, the variance of the first component is 10 and the variance of the second component is 0.1.
% Correlation, a (p + 1)-by-(p + 1) positive definite matrix specifying the prior correlation matrix of \beta for both components. The default is the identity matrix, which implies that the regression coefficients are uncorrelated, a priori.
% Probability, a (p + 1)-D vector of prior probabilities of variable inclusion (g_k, k = 0,...,_p_) or a function handle to a custom function. \gamma_j and \gamma_k, j \ne k, are independent, a priori. However, using a function handle (@functionname), you can supply a custom prior distribution that specifies dependencies between \gamma_j and \gamma_k. For example, you can specify forcing x_2 out of the model if x_4 is included.
% After you create a model, pass it and the data to estimate. The estimate function uses a Gibbs sampler to sample from the full conditionals, and estimate characteristics of the posterior distributions of \beta and \sigma^2. Also, estimate returns posterior estimates of g_k.
% 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® 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 with 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 = gcf;
h.Position(3) = h.Position(3)*2.5;
title('Variable Box Plots');

% The variables have fairly similar scales.
% To tune the prior Gaussian mixture variance factors, follow this procedure:
% Partition the data into estimation and forecast samples.
% Fit the models to the estimation sample and specify, for all k, V_{1k} = \{10,\;50,\;100\} and V_{2k} = \{0.05,\;0.1,\;0.5\}.
% Use the fitted models to forecast responses into the forecast horizon.
% Estimate the forecast MSE for each model.
% Choose the model with the lowest forecast MSE.
% George and McCulloch suggest another way to tune the prior variances of \beta in docid:econ_ug.mw_932b44e0-518f-4376-a7a4-22f208097e83.
%
% 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);

%
%% Create Prior Bayesian Linear Regression Models
% Create prior Bayesian linear regression models for SSVS by calling bayeslm and specifying the number of predictors, model type, predictor names, and component variance factors. Assume that \beta and \sigma^2 are dependent, a priori (mixconjugateblm model).

V1 = [10 50 100];
V2 = [0.05 0.1 0.5];
numv1 = numel(V1);
numv2 = numel(V2);

PriorMdl = cell(numv1,numv2); % Preallocate

for k = 1:numv2
    for j = 1:numv1
        V = [V1(j)*ones(p + 1,1) V2(k)*ones(p + 1,1)];
        PriorMdl{j,k} = bayeslm(p,'ModelType','mixconjugateblm',...
            'VarNames',predictornames,'V',V);
    end
end

% PriorMdl is a 3-by-3 cell array, and each cell contains a mixconjugateblm model object.
%
% Plot the prior distribution of log_GDP for the models in which V2 is 0.5.

for j = 1:numv1
    [~,~,~,h] = plot(PriorMdl{j,3},'VarNames',"log_GDP");
    title(sprintf("Log GDP, V1 = %g, V2 = %g",V1(j),V2(3)));
    h.Tag = strcat("fig",num2str(V1(j)),num2str(V2(3)));
end

% The prior distributions of \beta have the spike-and-slab shape. When V1 is low, more of the distribution is concentrated around 0, which makes it more difficult for the algorithm to attribute a high value for beta. However, variables the algorithm identifies as important are regularized, in that the algorithm does not attribute a high magnitude to the corresponding coefficients.
% When V1 is high, more density occurs well away from zero, which makes it easier for the algorithm to attribute non-zero coefficients to important predictors. However, if V1 is too high, then important predictors can have inflated coefficients.
%
%% Perform SSVS Variable Selection
% To perform SSVS, estimate the posterior distributions by using estimate. Use the default options for the Gibbs sampler.

PosteriorMdl = cell(numv1,numv2);
PosteriorSummary = cell(numv1,numv2);

rng(1); % For reproducibility
for k = 1:numv2
    for j = 1:numv1
        [PosteriorMdl{j,k},PosteriorSummary{j,k}] = estimate(PriorMdl{j,k},X,y,...
            'Display',false);
    end
end

% Each cell in PosteriorMdl contains an empiricalblm model object storing the full conditional posterior draws from the Gibbs sampler. Each cell in PosteriorSummary contains a table of posterior estimates. The Regime table variable represents the posterior probability of variable inclusion (g_k).
%
% Display a table of posterior estimates of g_k.

RegimeTbl = table(zeros(p + 2,1),'RowNames',PosteriorSummary{1}.Properties.RowNames);
for k = 1:numv2
    for j = 1:numv1
        vname = strcat("V1_",num2str(V1(j)),"__","V2_",num2str(V2(k)));
        vname = replace(vname,".","p");
        tmp = table(PosteriorSummary{j,k}.Regime,'VariableNames',vname);
        RegimeTbl = [RegimeTbl tmp];
    end
end
RegimeTbl.Var1 = [];
RegimeTbl

% Using an arbitrary threshold of 0.10, all models agree that FEDFUNDS is an insignificant or redundant predictor. When V1 is high, TB3MS borders on being insignificant.
%
% Forecast responses and compute forecast MSEs using the estimated models.

yhat = zeros(fh,numv1*numv2);
fmse = zeros(numv1,numv2);

for k = 1:numv2
    for j = 1:numv1
        idx = ((k - 1)*numv1 + j); 
        yhat(:,idx) = forecast(PosteriorMdl{j,k},XF);
        fmse(j,k) = sqrt(mean((yF - yhat(:,idx)).^2));
    end
end

%
% Identify the variance factor settings that yield the minimum forecast MSE.

minfmse = min(fmse,[],'all');
[idxminr,idxminc] = find(abs(minfmse - fmse) < eps);
bestv1 = V1(idxminr)
bestv2 = V2(idxminc)

%
% Estimate an SSVS model using the entire data set and the variance factor settings that yield the minimum forecast MSE.

XFull = [X; XF];
yFull = [y; yF];
EstMdl = estimate(PriorMdl{idxminr,idxminc},XFull,yFull);

% EstMdl is an empiricalblm model representing the result of performing SSVS. You can use EstMdl to forecast the unemployment rate given future predictor data , for example.
%
% Copyright 2018 The MathWorks, Inc.