%% Estimate Posterior Probability Using Monte Carlo Simulation
% Consider the linear regression model in docid:econ_ug.bvlxeoy-1.
% Create a prior model for the regression coefficients and disturbance variance, then estimate the marginal posterior distributions.

load Data_NelsonPlosser
VarNames = {'IPI'; 'E'; 'WR'};
X = DataTable{:,VarNames};
y = DataTable{:,'GNPR'};  

Mdl0 = fitlm(X,y);
coeff = Mdl0.Coefficients.Estimate;
CoeffCov = Mdl0.CoefficientCovariance;
coeffSE = Mdl0.Coefficients.SE;
CoeffCorr = diag(1./coeffSE)*CoeffCov*diag(1./coeffSE);
rng(1); % For reproducibility
n = 1e6;
dof = 1;
a = 1;
b = 10;
BetaDraws = mvtrnd(CoeffCorr,dof,n)'.*coeffSE + coeff;
Sigma2Draws = 1./gamrnd(a,1/b,1,n);

p = 3;
PriorMdl = bayeslm(p,'ModelType','empirical','BetaDraws',BetaDraws,...
    'Sigma2Draws',Sigma2Draws,'VarNames',VarNames);
PosteriorMdl = estimate(PriorMdl,X,y);

% Suppose that if the coefficient of real wages is below 2.5, then a policy is enacted. Although the posterior distribution of WR is known, and so you can calculate probabilities directly, you can estimate the probability using Monte Carlo simulation instead.
%
% Draw 1,000,000 samples from the marginal posterior distribution of \beta.

NumDraws = 1e6;
BetaSim = simulate(PosteriorMdl,'NumDraws',NumDraws);

% BetaSim is a 4-by-1,000,000 matrix containing the draws. Rows correspond to the regression coefficient and columns to successive draws.
%
% Isolate the draws corresponding to the coefficient of real wages (WR), and then identify which draws are less than 2.5.

isWR = PosteriorMdl.VarNames == "WR";
wrSim = BetaSim(isWR,:);
isWRLT2p5 = wrSim < 2.5;

%
% Find the marginal posterior probability that the regression coefficient of WR is below 2.5 by computing the proportion of draws that are less than 2.5.

probWRLT2p5 = mean(isWRLT2p5)

% The posterior probability that the coefficient of real wages is less than 2.5 is about 0.56.
%
% Copyright 2018 The MathWorks, Inc.