%% Inspect Augmented Data
% Consider the 3-D VAR(4) model of docid:econ_ug.mw_93ab52c2-232a-4877-8152-03c57b96825f.
% 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);

% Remove all rows that contain leading missing values.

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

% Create a weak conjugate prior model. 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 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;

% Draw 1000 parameter sets from the posterior distribution by calling simsmooth. Return the values that the simulation smoother imputes for the missing observations.

[~,~,NaNDraws] = simsmooth(PriorMdl,MissDataTable);

% NaNDraws is a 19-by-1000 matrix of randomly drawn response vectors from the posterior predictive distribution. Elements correspond to the missing values in the data ordered by a column-wise search. For example, NaNDraws(3,1) is the first randomly drawn imputed response of the third missing value in the data.  Find the corresponding value in the data.

[idxi,idxj] = find(ismissing(MissDataTable),3);
responsename = seriesnames(idxj(end))
observationtime = MissDataTable.Time(idxi(end))

% Plot the empirical distribution of the imputed values of the inflation rate during Q3-65.

histogram(NaNDraws(3,:))
title('Q3-65 Inflation Rate Empirical Distribution')

% Copyright 2020 The MathWorks, Inc.