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

p = 3;
PriorMdl = bayeslm(p,'ModelType','diffuse','VarNames',["IPI" "E" "WR"]);

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

PosteriorMdl = estimate(PriorMdl,X,y);

% Extract the posterior mean of \beta from the posterior model, and the posterior covariance of  \beta from the estimation summary returned by summarize.

estBeta = PosteriorMdl.Mu;
Summary = summarize(PosteriorMdl);
estBetaCov = Summary.Covariances{1:(end - 1),1:(end - 1)};

%
% 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 you can calculate probabilities directly, you can estimate the probability using Monte Carlo simulation instead.
% Draw 1e6 samples from the marginal posterior distribution of \beta.

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

% BetaSim is a 4-by- 1e6 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, 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.53.
%
% The marginal posterior distribution of the coefficient of WR is a t_{58}, but centered at 2.47 and scaled by 0.37. Directly compute the posterior probability that the coefficient of WR is less than 2.5.

center = estBeta(isWR);
stdBeta = sqrt(diag(estBetaCov));
scale = stdBeta(isWR);
t = (2.5 - center)/scale;
dof = 68;
directProb = tcdf(t,dof)

% The posterior probabilities are nearly identical.
%
% Copyright 2018 The MathWorks, Inc.