%% Simulate Regimes from SSVS Predictor Selection
% Consider the regression model in docid:econ_ug.bvj_ani-2.
% Assume these prior distributions for k = 0,...,3:
% \beta_k\vert\sigma^2,\gamma_k = \gamma_k\sigma\sqrt{V_{k1}}Z_1 + (1 - \gamma_k)\sigma\sqrt{V_{k2}}Z_2, where Z_1 and Z_2 are independent, standard normal random variables.  Therefore, the coefficients have a Gaussian mixture distribution. Assume all coefficients are conditionally independent, a priori, but they are dependent on the disturbance variance.
% \sigma^2 \sim IG(A,B). A and B are the shape and scale, respectively, of an inverse gamma distribution.
% \gamma_k \in \left\lbrace 0,1\right\rbraceand it represents the random variable-inclusion regime variable with a  discrete uniform distribution.
% Create a prior model for performing SSVS. Assume that \beta and \sigma^2 \;are dependent (a  conjugate mixture model). Specify the number of predictors p and the names of the regression coefficients.

p = 3;
PriorMdl = mixconjugateblm(p,'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'};

% Compute the number of possible regimes, that is, number of combinations that result from including and excluding variables in the model.

cardRegime = 2^(PriorMdl.Intercept + PriorMdl.NumPredictors)

%
% Simulate 10,000 regimes from the posterior distribution.

rng(1);
[~,~,RegimeSim] = simulate(PriorMdl,X,y,'NumDraws',10000);

% RegimeSim is a 4-by-1000 logical matrix. Rows correspond to the variables in Mdl.VarNames, and columns correspond to draws from the posterior distribution.
% Plot a histogram of the regimes visited. Recode the regimes so that they are readable. Specifically, for each regime, create a string that identifies the variables in the model, and separate the variables with dots.

cRegime = num2cell(RegimeSim,1);
cRegime = categorical(cellfun(@(c)join(PriorMdl.VarNames(c),"."),cRegime));
cRegime(ismissing(cRegime)) = "NoCoefficients";

histogram(cRegime);
title('Variables Included in Models')
ylabel('Frequency');

% Compute the marginal posterior probability of variable inclusion.

table(mean(RegimeSim,2),'RowNames',PriorMdl.VarNames,...
    'VariableNames',"Regime")

%
% Copyright 2018 The MathWorks, Inc.