%% Estimate Posterior Probability Using Monte Carlo Simulation
% Consider the linear regression model in Estimate Marginal Posterior Distributions.
% Create a prior model for the regression coefficients and disturbance variance, then estimate the marginal posterior distributions. Turn the estimation display off.

dof = 50;
C = eye(4);
ct = [-25; 4; 0; 3];
st = ones(4,1);
a = 3;
b = 1;
logPDF = @(params)priorMVTIG(params,ct,st,dof,C,a,b);

p = 3;
PriorMdl = bayeslm(p,'ModelType','custom','LogPDF',logPDF,...
    'VarNames',["IPI" "E" "WR"]);

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

width = [20,0.5,0.01,1,20];
thin = 10;
numDraws = 1e5/thin;

rng(1) % For reproducibility
PosteriorMdl = estimate(PriorMdl,X,y,'Width',width,'Thin',thin,...
    'NumDraws',numDraws,'Display',false);

% Estimate posterior distribution summary statistics for \beta by using the draws from the posterior distribution stored in posterior model.

estBeta = mean(PosteriorMdl.BetaDraws,2);
EstBetaCov = cov(PosteriorMdl.BetaDraws');

%
% Suppose that if the coefficient of real wages (WR)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;
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 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 WR is less than 2.5 is about 0.4430.
%
% Copyright 2012 The MathWorks, Inc.