%% Simulate Parameters from Posterior Distribution in Presence of Missing Values
% When simulate estimates the posterior distribution from which to draw parameters, it removes all rows in the data that contain at least one missing value (NaN). However, simsmooth uses Bayesian data augmentation to impute non-presample missing values during posterior estimation.
% Consider the 3-D VAR(4) model for the US inflation (INFL), unemployment (UNRATE), and federal funds (FEDFUNDS) rates.
% \left[ {\begin{array}{*{20}{l}}
% {\texttt{INFL}_t}\\
% {\texttt{UNRATE}_t}\\
% {\texttt{FEDFUNDS}_t}
% \end{array}} \right] = c + \sum\limits_{j = 1}^4 {{\Phi _j}\left[ {\begin{array}{*{20}{l}}
% {\texttt{INFL}_{t-j}\\
% {\texttt{UNRATE}_{t-j}\\
% {\texttt{FEDFUNDS}_{t-j}
% \end{array}} \right]} + \left[ {\begin{array}{*{20}{c}}
% {{\varepsilon _{1,t}}}\\
% {{\varepsilon _{2,t}}}\\
% {{\varepsilon _{3,t}}}
% \end{array}} \right].
% For all t, \varepsilon_t is a series of independent 3-D normal innovations with a mean of 0 and covariance \Sigma. Assume a weak conjugate prior distribution for the parameters \left({\left\lbrack \Phi_1 ,\ldotp \ldotp \ldotp ,\Phi_4 ,c\right\rbrack }^{\prime } ,\Sigma \right).
%% Load and Preprocess Data
% Load the US macroeconomic data set.  Compute the inflation rate, and stabilize the unemployment and federal funds rates.

load Data_USEconModel
seriesnames = ["INFL" "UNRATE" "FEDFUNDS"];
DataTable.INFL = 100*[NaN; price2ret(DataTable.CPIAUCSL)];
DataTable.DUNRATE = [NaN; diff(DataTable.UNRATE)];
DataTable.DFEDFUNDS = [NaN; diff(DataTable.FEDFUNDS)];
seriesnames(2:3) = "D" + seriesnames(2:3);

% Several series have leading NaN values because their measurements were not available at the same time as other measurements. Because leading NaN values can interfere with presample specification, remove all rows containing at least one leading missing value.

rmldDataTable = rmmissing(DataTable(:,seriesnames));

%
%% Create Prior Model
% Create a weak conjugate prior model by specifying large coefficient prior variances. Specify the response series names.

numseries = numel(seriesnames);
numlags = 4;

PriorMdl = conjugatebvarm(numseries,numlags,'SeriesNames',seriesnames);
numcoeffseqn = size(PriorMdl.V,1);
PriorMdl.V = 1e4*eye(numcoeffseqn);

%
%% Randomly Place Missing Values in Data
% To illustrate simulation in the presence of missing values, randomly place missing values in the data after the presample period.

rng(1)                          % For reproducibility
T = size(rmldDataTable,1);
idxpre = 1:PriorMdl.P;          % Presample period
idxest = (PriorMdl.P + 1):T;    % Estimation period 

nmissing = 20;                  % Simulate at most nmissing missing values
sidx = [randsample(idxest,nmissing,true); randsample(1:numseries,nmissing,true)];
lsidx = sub2ind([T,numseries],sidx(1,:),sidx(2,:));
MissData = table2array(rmldDataTable);
MissData(lsidx) = NaN;
MissDataTable = rmldDataTable;
MissDataTable{:,:} = MissData;

%
%% Simulate Parameters from Posterior
% Draw 1000 parameter sets from the posterior distribution by calling simsmooth, which estimates the posterior distribution of the parameters, and then forms the posterior predictive distribution.

[Coeff,Sigma] = simsmooth(PriorMdl,MissDataTable);

% Coeff is a 39-by-1000 matrix of randomly drawn coefficient vectors from the posterior. Sigma is a 3-by-3-by-1000 array of randomly drawn innovations covariance matrices.
% By default, simsmooth initializes the VAR model by using the first four observations in the data.
% To associate rows of Coeff to coefficients, obtain a summary of the prior distribution by using summarize. In a table, display the first set of randomly drawn coefficients with corresponding names.

Summary = summarize(PriorMdl,'off');
table(Coeff(:,1),'RowNames',Summary.CoeffMap)

% Alternatively, you can create an empirical model from the draws, and use summarize to display the model by specifying any display option.
% Display a summary of the posterior draws as an equation.

EmpMdl = empiricalbvarm(numseries,numlags,'SeriesNames',seriesnames,...
    'CoeffDraws',Coeff,'SigmaDraws',Sigma);
summarize(EmpMdl,'equation')

% Copyright 2020 The MathWorks, Inc.