%\VignetteIndexEntry{Dynamic Model Averaging for Practitioners in Economics and Finance: The eDMA Package}
%\VignetteKeywords{eDMA}
%\VignettePackage{eDMA}
\documentclass[nojss,shortnames]{jss}
\usepackage{amsfonts}
\usepackage{amsmath}
\usepackage{bbm}
\usepackage{rotating}
\usepackage{booktabs}
\usepackage{multirow}
%:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
% INPUT DEFINITIONS
% :::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
%\input{DEF.tex}
\def\qmo{``}
\def\qmc{''}
\def\qmcsp{'' }

\def\btheta{\mbox{\boldmath $\theta$}}
\def\beeta{\mbox{\boldmath $\eta$}}
\def\bF{\mathbf{F}}
\def\bW{\mathbf{W}}
\def\bR{\mathbf{R}}
\def\bC{\mathbf{C}}
\newcommand{\indep}{\rotatebox[origin=c]{90}{$\models$}}

\newcommand{\insertfloat}[1]{%
	\begin{center}
		[Insert #1 about here]%
	\end{center}%
}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% declarations for jss.cls %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%% almost as usual
\author{Leopoldo Catania\\University of Rome, \qmo Tor Vergata\qmcsp \And
        Nima Nonejad\\ Aalborg University and CREATES}
\title{Dynamic Model Averaging for Practitioners in Economics and Finance: The \pkg{eDMA} Package}

%% for pretty printing and a nice hypersummary also set:
\Plainauthor{Leopoldo Catania, Nima Nonejad} %% comma-separated
\Plaintitle{Dynamic Model Averaging for Practitioners in Economics and Finance: The eDMA Package} %% without formatting
\Shorttitle{\pkg{eDMA}: Efficient Dynamic Model Averaging} %% a short title (if necessary)

%% an abstract and keywords
\Abstract{
\cite{raftery_etal.2010} introduce an estimation technique referred to as Dynamic
Model Averaging (DMA). In their application, DMA is used to the problem
of predicting the output strip thickness for a cold rolling mill,
where the output is measured with a time delay. Recently, DMA has also shown to be
very useful in macroeconomic and financial applications. In this paper, we present the \pkg{eDMA} package for DMA estimation in R, which is especially suited for practitioners in economics
and finance. Our implementation proves to be up to 133 times faster then a standard implementation on a single--core CPU. With the help of this package, practitioners are able perform DMA on a standard
PC without resorting to large clusters, which are not easily available to all researchers.
We demonstrate the usefulness of this package through simulation
experiments and an empirical application using quarterly U.S. inflation data.
}
\Keywords{Dynamic model averaging, Multi core CPU, Parallel computing, \proglang{R}, \code{OpenMP}}
\Plainkeywords{Dynamic model averaging, Multi core CPU, Parallel computing, R, OpenMP} %% without formatting
%% at least one keyword must be supplied

%% publication information
%% NOTE: Typically, this can be left commented and will be filled out by the technical editor
%% \Volume{50}
%% \Issue{9}
%% \Month{June}
%% \Year{2012}
%% \Submitdate{2012-06-04}
%% \Acceptdate{2012-06-04}

%% The address of (at least) one author should be given
%% in the following format:
\Address{
  Leopoldo Catania\\
  Department of Economics and Finance\\
  Faculty of Economics\\
  University of Rome, \qmo Tor Vergata\qmcsp\\
  Via Columbia, 2\\
  00133 Rome, Italy\\
  E-mail: \email{leopoldo.catania@uniroma2.it}\\
  %

  Nima Nonejad\\
  Department of Mathematical Sciences\\
  Aalborg University and CREATES\\
  Fredrik Bajers Vej 7G\\
  9220, Aalborg East, Denmark\\
  E-mail: \email{nimanonejad@gmail.com}\\
  %URL: \url{http://eeecon.uibk.ac.at/~zeileis/}
}


%% It is also possible to add a telephone and fax number
%% before the e-mail in the following format:
%% Telephone: +43/512/507-7103
%% Fax: +43/512/507-2851

%% for those who use Sweave please include the following line (with % symbols):
%% need no \usepackage{Sweave.sty}

%% end of declarations %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


\begin{document}

%% include your article here, just as usual
%% Note that you should use the \pkg{}, \proglang{} and \code{} commands.
%
\section[Introduction]{Introduction}
\label{sec: Introduction}
%
Modeling and forecasting economic variables such as real GDP, inflation and equity premium
is of clear importance for researchers in economics and finance. For instance, forecasting inflation is
crucial for central banks with regards to conducting optimal monetary policy. Similarly, understanding and predicting equity premium is one of the most widely important topics discussed in financial economics as it has great implications on portfolio choice and risk management, see for instance \cite{dangl_halling.2012} among many others.

In order to obtain the best forecast as possible, practitioners often
try to take advantage of the many potential predictors available and
seek to combine the information from these predictors in an optimal
way, see \cite{stock_watson.1999}, \cite{stock_watson.2008} and \cite{groen_etal.2013} just to mention a few references. In the context of economic applications, \cite{koop_korobilis.2011} and \cite{koop_korobilis.2012}, implement a technique developed by \cite{raftery_etal.2010}, referred to as Dynamic
Model Averaging (DMA). The original intent of the mentioned article is to predict the output strip thickness for a cold rolling
mill, where the output is measured with a time delay. Typically, DMA consists of thousands of models based on all possible combinations of the predictors available. Furthermore, besides allowing for time--variation in the regression coefficients of each individual model, DMA also allows the relevant model set to change with time as well by introducing a forgetting factor. \cite{koop_korobilis.2011} and \cite{koop_korobilis.2012} argue (very elegantly) that by slightly adjusting the original framework, DMA can be very useful in economic contexts, especially inflation forecasting.\footnote{Specifically, \cite{koop_korobilis.2012} change the conditional volatility formula of \cite{raftery_etal.2010} arguing that their original formula is not suited in the context of economic data. Their suggestion proves more compatible with the behavior of economic data such as inflation, see also \cite{dangl_halling.2012} for a similar approach.} \cite{dangl_halling.2012} provide further suggestions on how to
improve DMA such that it better suits the nature of economic and financial data.
%
The aforementioned authors, also provide a very useful variance decomposition scheme using the output from the estimation procedure to understand which source of uncertainty explains the degree of variation in the observed variable. \cite{byrne_etal.2014}, among others, use the modifications proposed in \cite{dangl_halling.2012} to model currency exchange--rate behavior. We must also emphasize that DMA is not solely limited to these series and can be used in a wide range of economic applications such as: Forecasting realized volatility as well as house, oil and commodity prices.

However, designing an efficient DMA algorithm remains a challenging issue.
As previously mentioned, DMA considers all possible combinations from
a set of predictors at each point in time. Typically, many candidate variables
are available and, as a consequence, it poses a limit given the computational facilities at hand, which for many practitioners typically consists of a standard 8 core CPU. Thus, while handling a relatively small number of model combinations, usually between 1000 to 3000, allows one to perform DMA using standard loops and software, handling larger number of combinations becomes very burdensome, especially in the context of memory consumption, see also \cite{koop_korobilis.2012}.

In order to deal with this issue,
\cite{onorante_raftery.2016} suggest a strategy that considers not the
whole model space, but rather a subset of models and dynamically optimizes
the choice of models at each point in time. However, \cite{onorante_raftery.2016} have to assume that models
do not change too fast over time, which is not an ideal assumption when dealing with financial and in some cases monthly economic data. Furthermore, the approach of \cite{onorante_raftery.2016} is not able to clearly state why certain predictors can be considered as less important than others or vice versa. Finally, incorporating the modifications suggested in \cite{dangl_halling.2012}, which are very important with regards to interpretation of results, is not possible using the approach of \cite{onorante_raftery.2016}.

% in this paper
In this paper, we introduce the \pkg{eDMA} package for \proglang{R} \citep{R.2015}, which efficiently implements a DMA procedure based on \cite{raftery_etal.2010} and \cite{dangl_halling.2012}. The routines in the \pkg{eDMA} package are principally written in \proglang{C++} using the \code{armadillo} library of \cite{sanderson.2010} and then made available in \proglang{R} exploiting the \pkg{Rcpp} and \pkg{RcppArmadillo} packages of \cite{Rcpp} and \cite{RcppArmadillo}, respectively. Furthermore, the \code{OpenMP} API \citep{openmp.2008} is used to speedup the computations when a shared memory multiple processors hardware is available, which, nowadays, is standard for the majority of commercial laptops. However, if the hardware does not have multiple processors, the \pkg{eDMA} package can still be used with the classical sequential CPU implementation.

Our aim is to provide a package that can be used by a broad audience from different academic fields who are interested in implementing DMA in their research. Furthermore, our package enables practitioners, to perform DMA using a large number of predictors without needing to understand and possibly implement complex programming concepts such as \qmo how to efficiently allocate memory\qmc, or \qmo how to efficiently parallelize the computations\qmc.

It is worth noting that, within the \proglang{R} environment, the package \pkg{dma} of \cite{dma} downloadable from CRAN can be used to perform the DMA of \cite{raftery_etal.2010}. However, \pkg{dma} has several weaknesses such as (i): Does not allow for the extensions of \cite{dangl_halling.2012}, which are important in the context of interpreting the amount of time--variation in the regression coefficients and performing a variance decomposition analysis, (ii): It is remarkably slow compared to our alternative, (iii): Requires a very large amount of RAM when executed for moderately large applications, (iv): Does not allow for parallel computing. We refer the reader interested in these aspects to Section 4, where we report a comparative analysis between \pkg{dma} and \pkg{eDMA} using simulated data. \pkg{eDMA} permits us to also perform Bayesian Model Averaging (BMA) and Bayesian Model Selection (BMS) for linear regression models with constant coefficients implemented, for example, in the \proglang{R} packages \pkg{BMA} \citep{BMA} and \pkg{BMS} \citep{zeugner_feldkircher.2015}.

The empirical application in Section \ref{sec:EmpiricalApplication} presents different features of this package, provides illustrations on how practitioners can implement DMA and more importantly use the output from the estimation procedure to interpret results.

The structure of this paper is as follows: Sections \ref{sec:TheModel} and \ref{sec:ModifiedDMA} briefly introduce
DMA and its extensions. Section \ref{sec:ThePackage} presents
the technical aspects. Section \ref{sec:computational} provides an intuitive description of the challenges
that DMA posses from a computational point of view and proposes solutions. Section \ref{sec:EmpiricalApplication} provides an empirical application
to demonstrate the advantages of \pkg{eDMA} from a practical
point of view. Therefore, practitioners who are solely interested
on how to implement DMA using the \pkg{eDMA} package can skip Sections \ref{sec:TheModel} and \ref{sec:ModifiedDMA}. Finally, Section \ref{sec:conclusion}
concludes.

\section[The model]{The model}
\label{sec:TheModel}

In this section, we briefly introduce the DMA algorithm of \cite{raftery_etal.2010}.
We then pinpoint some of the main challenges associated with applying
DMA in economic and financial applications. Thereafter, we propose
solutions in Section \ref{sec:computational}.

Let $y_t$ denote the dependent variable at time $t$. Our goal is to model $y_t$ using a pool of $i = 1,\dots,k$, candidate Dynamic Linear Models (DLM) of the types introduced in \cite{west_harrison.1999} and \cite{raftery_etal.2010},

\begin{align}
y_{t} & =  \bF_{t}^{\left(i\right)\prime}\btheta_{t}^{\left(i\right)}+\varepsilon_{t}^{\left(i\right)},\quad\varepsilon_{t}^{\left(i\right)}\sim N\left(0,V_{t}^{\left(i\right)}\right)\label{eq:2.1}\\
\btheta_{t}^{\left(i\right)} & =  \btheta_{t-1}^{\left(i\right)}+\beeta_{t}^{\left(i\right)},\quad\beeta_{t}^{\left(i\right)}\sim N\left(0,\bW_{t}^{\left(i\right)}\right),\label{eq:2.2}
\end{align}

where $\bF^{(i)}_t$ is a subset from the total $n$ predictors. Let $p$ denote the numbers of
predictors in $\bF_{t}^{\left(i\right)}$ for model $i$. Then, $\btheta_{t}^{\left(i\right)}$
is a $p\times1$ vector of time--varying regression coefficients, which
evolve according to (\ref{eq:2.2}) and determine the impact of $\bF^{(i)}_t$ on $y_t$. Note, we do not assume any
systematic movements in $\btheta_{t}^{\left(i\right)}$. On the contrary,
we consider changes in $\btheta_{t}^{\left(i\right)}$ as unpredictable.\footnote{See \cite{dangl_halling.2012} and \cite{koop_korobilis.2012} for a similar
model specification.}

The conditional variances, $V^{\left(i\right)}_{t}$ and $\bW^{\left(i\right)}_{t}$, for model $i$ are the unknown quantities associated with the observational equation, \eqref{eq:2.1}, and the state equation, \eqref{eq:2.2}. Obviously, when $\bW^{\left(i\right)}_{t}=\boldsymbol{0}$ for $t=1,\dots,T$, then $\btheta^{\left(i\right)}_{t}$ is constant over time. Thus, our model can nest the specification of constant regression coefficients. As $\bW^{\left(i\right)}_{t}$ increases, the variation in $\btheta^{\left(i\right)}_{t}$ varies according to Equation~\ref{eq:2.2}. However, this does not mean that $\btheta^{(i)}_t$ needs to change at every time period. For instance, we can simply have periods where $\bW^{(i)}_t=0$ in which $\btheta^{(i)}_t=\btheta^{(i)}_{t-1}$ and other periods when  $\btheta^{(i)}_t$ is allowed to change according to Equation~\ref{eq:2.2}.\footnote{As indicated below, we model time--variation in $\btheta_t^{(i)}$ through a forgetting factor, $\delta$. The way we update $\delta$, which determines the magnitude of the shocks that hit $\btheta^{(i)}_t$ at each $t$, avoids any unreasonable behavior of $\btheta_t^{(i)}$. Thus, we do not need to put any restrictions on $\btheta_t^{(i)}$ itself, see \cite{dangl_halling.2012} for a similar approach.}

In DMA, we consider a total of $k=2^{n}-1$ possible combinations of the predictors at each point in time while contemporaneously assuming that $\btheta^{(i)}_t$ changes according to Equation~\ref{eq:2.2}.\footnote{The model $y_t=\varepsilon_t$ is \textbf{not} considered in the universe of models, see also \cite{dangl_halling.2012}.} DMA then averages forecasts/predictions across these combinations using a recursive updating scheme based on the predictive likelihood. The predictive likelihood measures the ability of a model to predict $y_t$, thus making it
the central quantity of interest for model evaluation. Apparently, models containing important combinations of predictors receive high predictive likelihood values, which means that these models obtain higher weights in the averagning process. Besides averaging, we can also use the prediction/forecasts of the model receiving the highest probability among all model combinations considered at each point in time. In this case, we are performing Dynamic Model Selection (DMS), see also \cite{koop_korobilis.2012}.

In the context of estimation, DMA avoids the difficult
task of specifying $\bW^{\left(i\right)}_{t}$ for each individual model. Instead DMA relies on using a forgetting
factor, $0<\delta\leq1$. This parameter describes the loss of information
through time, which simplifies things greatly from a practical point of view. The forgetting factor avoids the need to estimate $\bW^{(i)}_t$ for each individual model. Particularly, using the same notation as the Appendix in \cite{dangl_halling.2012}, we define the variables in the Kalman recursions of the $i$--t model as: (i): $\bR^{\left(i\right)}_{t}$, the unconditional variance of $\btheta^{\left(i\right)}_t$ (see Equation 14 in the Appendix of \cite{dangl_halling.2012}), (ii): $\bC^{\left(i\right)}_{t}$, the estimator for the variance of $\btheta^{\left(i\right)}_t$, (see Equation 20 in the Appendix of \cite{dangl_halling.2012}), and (iii): $S^{\left(i\right)}_{t}$, the estimator of the observational variance (see Equation 17 in the Appendix of \cite{dangl_halling.2012}).
%
Then, using $\delta$, we can rewrite $\bR^{\left(i\right)}_{t}=\delta^{-1}\bC^{\left(i\right)}_{t-1}$, indicating that there is a relationship between $\bW^{\left(i\right)}_{t}$ and $\delta$, which is given as $\bW^{\left(i\right)}_{t}=\left(1-\delta\right)/\delta \bC^{\left(i\right)}_{t-1}$. In other words, the loss of information is proportional to the covariance of the state parameters.
This way, we can control the magnitude of the shocks that affect $\btheta^{\left(i\right)}_{t}$
by adjusting $\delta$ instead of directly estimating $\bW^{(i)}_t$. Accordingly, $\delta=1$ corresponds to
$\bW^{\left(i\right)}_{t}=\boldsymbol{0}$, which means that $\btheta_{t}^{(i)}$ is constant over time. For
$\delta<1$, we introduce time--variation in $\btheta^{\left(i\right)}_{t}$. For instance, when $\delta=0.99$, in the context quarterly data, observations five years ago receive approximately $80\%$ as much weight as last period's observation, which corresponds to gradual time--variation in $\btheta^{(i)}_t$. When $\delta=0.95$, observations $20$ periods ago receive only about $35\%$ as much weight as last period's observations, suggesting relatively more drastic change in $\btheta^{(i)}_t$ at each point in time. Evidently, while this renders the model more
flexible to adapt to changes in $y_{t}$, the increased variability
in $\btheta^{\left(i\right)}_{t}$ translates into high prediction variance. Thus, finding the correct value of $\delta$ at each point in time is extremely important with regards to the ability of the model to capture the correct magnitude of variations in $\btheta^{(i)}_t$ an thus produce optimal forecasts/predictions. The same
line of argumentation holds for the parameter, $\alpha$, which controls the forgetting for the entire model set, see \cite{raftery_etal.2010} for more details. For instance, $\alpha=0.95$ means that past model performance is less important than when $\alpha=0.99$. We must also determine a way to model the evolution of $V_{t}^{\left(i\right)}$. Here, we have
several choices, which we go into more details below, see point (c).

To summarize, DMA depends on:

\begin{itemize}
\item[(a)] The number of predictors to consider. Typically, in economic applications, the set of predictors contains
exogenous predictors as well as lagged values of $y_{t}$. For instance,
in the context of forecasting quarterly inflation, besides considering
predictors such as unemployment rate and T--bill rates, \cite{koop_korobilis.2012} also consider the first three lags of $y_{t}$ as explanatory variables. Of course, how to distinguish useful variables from noise remains an important issue. For instance, following a recession or revisions of data, the set of predictors to use probably changes.

%
\item[(b)] The choice of the forgetting factors, $\alpha$ and $\delta$. In many applications
$\alpha\in\left\{ 0.98,0.99,1\right\} $ work very well and results do not change drastically across different values of $\alpha$. On the other
hand, as previously mentioned, we often find that the choice of $\delta$ is more important.
\cite{koop_korobilis.2012} fix $\delta$ at $\{0.95,0.98,0.99,1.00\}$ and run DMA using each of these values. They find that results differ considerably in terms
of out--of--sample forecasts, see \cite{koop_korobilis.2012} for more details. Evidently, in many economic applications,
it is very plausible that $\delta$ would indeed be time--varying.
For instance, it is plausible to expect that $\delta$ is relatively low in
recessions or periods of market turmoil (as there is considerable
time--variation in $\btheta^{\left(i\right)}_{t}$ in these periods). Conversely, $\delta$ is expected
to be close to $1.00$ during tranquil and periods of low volatility. \cite{dangl_halling.2012} propose a very elegant solution to this problem by considering a grid of values for $\delta$ and incorporate this in the DMA setting by averaging over all possible combinations of the predictors as well as over a grid value for $\delta$. Furthermore, this procedure can also be used to obtain more information from the data through a variance decomposition scheme, see below for more details.
%
\item[(c)] Modeling $V_{t}^{\left(i\right)}$: In this paper, we make
things easy for conjugate analysis by using the following assumptions, see also \cite{dangl_halling.2012} and \cite{byrne_etal.2014} for the same approach:
We assume that $V_{t}^{\left(i\right)}=V^{\left(i\right)}$ for all
$t$. For time $t=0$, we specify a Normal prior on $\btheta_{0}^{\left(i\right)}$
and a Gamma prior on $\phi^{\left(i\right)}=V^{-1\left(i\right)}$.
The time $t$ estimate of the variance of the error term in (\ref{eq:2.1})
is then given as Equation 17 in the Appendix of \cite{dangl_halling.2012}. More importantly, by using these assumptions,
we find that, when we integrate the conditional density of $y_{t}$
over the values of $\btheta_{t}^{\left(i\right)}$ and $V^{\left(i\right)}$
to obtain the predictive density, $p\left(y_{t}\vert\mathcal{F}_{t-1}\right)$, the corresponding density has a closed--form solution, which is
given as $p\left(y_{t}\vert\mathcal{F}_{t-1}\right)\sim t_{n^{\left(i\right)}_{t}}\left(\hat y^{\left(i\right)}_{t},Q^{\left(i\right)}_{t}\right)$,
where $t_{n_{t}^{\left(i\right)}}$ stands for the Student--t distribution
with $n_{t}^{\left(i\right)}$ degrees--of--freedom and mean and scale
given by $\hat y^{\left(i\right)}_{t}$ and $Q_{t}^{\left(i\right)}$, see Equations 13 and 16 of \cite{dangl_halling.2012},
respectively.
\end{itemize}

We can also follow \cite{koop_korobilis.2012}
and use an Exponentially Weighted Moving Average (EWMA) estimate of
$V_{t}^{\left(i\right)}$. However, this requires the practitioner to also consider an additional parameter, namely, $\kappa$, which increases the computation burden. Furthermore, by experimenting with EWMA on a smaller model, we find that $\delta$ and $\kappa$ in many ways are intertwined, in the sense that we obtain the same magnitudes of variation in $\btheta^{(i)}_t$ as for (c) when we allow $\kappa$ and $\delta$ to vary over time. In this case, compared to (c), we obtain higher estimates of $\delta$ over time whereas we estimate $\kappa$ close to $0.96$, which is the value suggested in \cite{koop_korobilis.2012}. We observe the same phenomena when we allow $\alpha$ to vary with $\delta$. Overall, our conclusion is that it is best to use (c) and fix $\alpha$ close to $0.99$ for monthly and quarterly data. This way, we maintain a very parsimonious model structure and are never in doubt with regards to under (over) estimating the true magnitude of variation in $\delta$ ($\alpha$).

\section[Modified DMA]{Modified DMA}
\label{sec:ModifiedDMA}
%
Below, we present the DMA algorithm modified to incorporate the extensions mentioned above. We refer the reader to \cite{dangl_halling.2012} for more details.
%These Equations serve as a basis for the calculations in Section \ref{sec:FOO}.
Let $M_{i}$ denote a model containing a specific set of predictors
chosen from a set of $k=2^n-1$ candidates and $\delta_{j}\in\left\{ \delta_{1},...,\delta_{d}\right\} $
denote a specific choice of the degree of time--variation in the regression coefficients
at time $t$. The total posterior density of model $M_{i}$ and forgetting
$\delta_{j}$ at time $t$, $p\left(M_{i},\delta_{j}\vert\mathcal{F}_{t}\right)$,
is given as
\begin{align}
p\left(M_{i},\delta_{j}\vert\mathcal{F}_{t}\right) & =  p\left(M_{i}\vert\delta_{j},\mathcal{F}_{t}\right)p\left(\delta_{j}\vert\mathcal{F}_{t}\right).\nonumber
\end{align}
In order to obtain $p\left(M_{i}\vert\mathcal{F}_{t}\right)$ and $p\left(\delta_{j}\vert\mathcal{F}_{t}\right)$,
we can use that
\begin{align}
p\left(M_{i}\vert\mathcal{F}_{t}\right) & =  \sum_{j=1}^{d}p\left(M_{i}\vert\delta_{j},\mathcal{F}_{t}\right)p\left(\delta_{j}\vert\mathcal{F}_{t}\right).\label{eq:3.2}
\end{align}
The term, $p\left(M_{i}\vert\delta_{j},\mathcal{F}_{t}\right)$, in Equation~\ref{eq:3.2} is given as
\begin{align}
p\left(M_{i}\vert\delta_{j},\mathcal{F}_{t}\right) & =  \frac{p\left(y_{t}\vert M_{i},\delta_{j},\mathcal{F}_{t-1}\right)p\left(M_{i}\vert\delta_{j},\mathcal{F}_{t-1}\right)}{\sum_{l=1}^{k}p\left(y_{t}\vert M_{l},\delta_{j},\mathcal{F}_{t-1}\right)p\left(M_{l}\vert\delta_{j},\mathcal{F}_{t-1}\right)}\label{eq:3.3}
\end{align}
where
\begin{align}
p\left(M_{i}\vert\delta_{j},\mathcal{F}_{t-1}\right) & =  \frac{p\left(M_{i}\vert\delta_{j},\mathcal{F}_{t-1}\right)^{\alpha}}{\sum_{l=1}^{k}p\left(M_{l}\vert\delta_{j},\mathcal{F}_{t-1}\right)^{\alpha}}.\label{eq:3.4}
\end{align}
The second term on the right--hand side of Equation~\ref{eq:3.2} is given as
\begin{align}\label{eq:delta}
p\left(\delta_{j}\vert\mathcal{F}_{t}\right) & =  \frac{p\left(y_{t}\vert\delta_{j},\mathcal{F}_{t-1}\right)p\left(\delta_{j}\vert\mathcal{F}_{t-1}\right)}{\sum_{l=1}^{d}p\left(y_{t}\vert\delta_{l},\mathcal{F}_{t-1}\right)p\left(\delta_{l}\vert\mathcal{F}_{t-1}\right)}.
\end{align}
where
\begin{align}
p\left(\delta_{j}\vert\mathcal{F}_{t-1}\right) & =  \frac{p\left(\delta_{j}\vert\mathcal{F}_{t-1}\right)^{\alpha}}{\sum_{l=1}^{d}p\left(\delta_{l}\vert\mathcal{F}_{t-1}\right)^{\alpha}}.\nonumber
\end{align}
As previously mentioned, $p\left(y_{t}\vert M_{i},\delta_{j},\mathcal{F}_{t-1}\right)\sim t_{n_{t}}\left(\hat y_{i,t}^{\left(j\right)},Q_{i,t}^{(j)}\right)$,
where $\hat y_{i,t}^{\left(j\right)}$ and $Q_{i,t}^{(j)}$ are the quantities of model
$M_{i}$, $i=1,...,k$, conditional on $\delta_{j}$, $j=1,...,d$, and $\alpha$.
Typically, $p\left(M_{i},\delta_{j}\vert\mathcal{F}_{0}\right)=1/(d\cdot k)$
such that, initially, all model combinations and degrees of time--variation
are equally likely. Thereafter, as a new observation arrives, model probabilities are updated using
the above recursions.
%
\subsection[Using the output from DMA]{Using the output from DMA}
%
For practitioners, the most interesting output from DMA are:
\begin{itemize}
\item[(i)] The predictive mean of $y_{t+1}$ conditional on $\mathcal{F}_{t}, \hat y_{t+1}$.
This is simply an average of each of the individual model predictive
means. That is
\begin{align}
\hat y_{t+1} & =  \sum_{j=1}^{d}\E\left[y_{t+1}^{\left(j\right)}\vert\mathcal{F}_{t}\right]p\left(\delta_{j}\vert\mathcal{F}_{t}\right),\label{eq:4.5}
\end{align}
where
\begin{align}
\E\left[y_{t+1}^{\left(j\right)}\vert\mathcal{F}_{t}\right] & =  \sum_{i=1}^{k}\E\left[y_{i,t+1}^{\left(j\right)}\vert\mathcal{F}_{t}\right]p\left(M_{i}\vert\delta_{j},\mathcal{F}_{t}\right).\nonumber
\end{align}
%
The formulas for the predictive density are very similar. We only replace the predictive mean with the predictive density inside Equation~\ref{eq:4.5}, see \cite{dangl_halling.2012}. Besides averaging over the individual predictive means/densities, we can simply choose the predictive mean/density associated with the model with the highest posterior probability, see Equation~\ref{eq:3.2}. Henceforth, we label this as Dynamic Model Selection (DMS). Furthermore, as mentioned in \cite{koop_korobilis.2012}, for DMA  and DMS, with $\delta$ and $\alpha$ fixed at 1, we have Bayesian model averaging (BMA) and Bayesian Model Selection (BMS) based on exact predictive likelihood, see for instance \cite{zeugner_feldkircher.2015}.\footnote{\cite{zeugner_feldkircher.2015} also implement BMA using the MC$^3$ algorithm relying on Markov Chain Monte Carlo (MCMC) techniques. However, their framework does not allow for time--variation in the regression coefficients nor model size.}
%
\item[(ii)] Quantities such as the expected size of the predictor, $\E\left[Size_{t}\right]=\Sigma_{i=1}^{k}Size^{(i)}p\left(M_{i}\vert\mathcal{F}_{t}\right)$,
where $Size^{(i)}$ be the number of predictors in model
$i$. This quantity reveals the average number of predictors in the DMA, see \cite{koop_korobilis.2012}. Similarly, we can compute the number of predictors for the model with the highest posterior probability, \eqref{eq:3.2}, at each point in time, which give the optimal model size at time $t$.
%
\item[(iii)] Posterior inclusion probabilities for the predictor. That is,
at each $t$ we calculate $\sum_{i=1}^{k}1_{\left(i\subset m\right)}p\left(M_{i}\vert\mathcal{F}_{t}\right)$,
where $1_{\left(i\subset m\right)}$ is an indicator function taking
the value of either $0$ or $1$ and $m$ is the $m$th predictor. We can also report the highest posterior model probability or the sum of the top $10\%$ model probabilities among all model combinations after the effect of $\delta$ is integrated out. This information can be used to determine if there is a group or an individual model that obtains relatively high posterior probability.
\item[(iv)] Posterior weighted average of $\delta$ at each point in time, see Equation~\ref{eq:delta}.
%
\item[(v)] Posterior weighted average estimates of $\btheta_{t}$ for DMA,
\begin{align}\label{eq:theta_post}
\E\left[\btheta_{t}\vert\mathcal{F}_{t}\right] & =  \sum_{j=1}^{d}\E\left[\btheta_{t}^{\left(j\right)}\vert\mathcal{F}_{t}\right]p\left(\delta_{j}\vert\mathcal{F}_{t}\right),
\end{align}
where
\begin{align}
\E\left[\btheta_{t}^{\left(j\right)}\vert\mathcal{F}_{t}\right] & =  \sum_{i=1}^{k}\E\left[\theta_{i,t}^{\left(j\right)}\vert\mathcal{F}_{t}\right]p\left(M_{i}\vert\delta_{j},\mathcal{F}_{t}\right).\nonumber
\end{align}

\item[(vi)] Variance decomposition of the data, $\textrm{Var}\left(y_{t+1}\vert\mathcal{F}_t\right)$,
 decomposed into

\begin{align}
\VAR\left(y_{t+1}\vert\mathcal{F}_t\right) & =  \sum_{j=1}^{d}\left[\sum_{i=1}^{k}\left(S_{t}\vert M_{i},\delta_{j},\mathcal{F}_{t}\right)p\left(M_{i}\vert\delta_{j},\mathcal{F}_{t}\right)\right]p\left(\delta_{j}\vert\mathcal{F}_{t}\right)\nonumber\\
 & + \sum_{j=1}^{d}\left[\sum_{i=1}^{k}\left(\bF_{t}^\prime\bR_{t}\bF_{t}\vert M_{i},\delta_{j},\mathcal{F}_{t}\right)p\left(M_{i}\vert\delta_{j},\mathcal{F}_{t}\right)\right]p\left(\delta_{j}\vert\mathcal{F}_{t}\right)\nonumber\\
 & + \sum_{j=1}^{d}\left[\sum_{i=1}^{k}\left(\hat{y}_{t+1,i}^{\left(j\right)}-\hat{y}_{t+1}^{\left(j\right)}\right)^{2}p\left(M_{i}\vert\delta_{j},\mathcal{F}_{t}\right)\right]p\left(\delta_{j}\vert\mathcal{F}_{t}\right)\nonumber\\
 & + \sum_{j=1}^{d}\left(\hat{y}_{t+1}^{\left(j\right)}-\hat{y}_{t+1}\right)^{2}p\left(\delta_{j}\vert\mathcal{F}_{t}\right).\label{eq:varDec}
\end{align}

The first term is the observational variance, Obs. The remaining terms
are: Variance due to errors in the estimation of the coefficients,
Coeff, variance due to uncertainty with respect to the choice
of the predictors, Mod, and variance due to uncertainty with
respect to the choice of the degree of time--variation in the regression
coefficients, TVP, see \cite{dangl_halling.2012} for more details.
\end{itemize}

\section[The eDMA package for R]{The \pkg{eDMA} package for \proglang{R}}

\label{sec:ThePackage}
The \pkg{eDMA} package for \proglang{R} offers an integrated environment for practitioners in economics and finance to perform our DMA algorithm. It is principally written in \proglang{C++} exploiting the \pkg{armadillo} library of \cite{sanderson.2010} to speed up computations. The relevant functions are then made available in \proglang{R} through the \pkg{Rcpp} and \pkg{RcppArmadillo} packages of \cite{Rcpp} and \cite{RcppArmadillo}, respectively. It also makes use of the \code{OpenMP} API \citep{openmp.2008} to parallelize part of the routines needed to perform DMA. Furthermore, multiple processors are automatically used if supported by the hardware, however, as will be discussed later, the user is free to manage the level of resources used by the program.

The \pkg{eDMA} package is written using the S4 object oriented language, meaning that classes and methods are available in the code. Specifically, \proglang{R} users will find common methods such as \code{plot()}, \code{show()}, \code{as.data.frame()},\code{coef()} and \code{residuals()}, among others, in order to visualise the output of DMA and extract estimated quantities.

Package \pkg{eDMA} is available from GitHub at \url{https://github.com/LeopoldoCatania/eDMA} and can be installed using the command:\footnote{Note that, to build the \pkg{eDMA} package under Windows it is advised to install Rtools from \url{https://cran.r-project.org/bin/windows/Rtools/} before preceding to the next steps.}

\begin{CodeChunk}
\begin{CodeInput}
R> library("devtools")
R> install_github("LeopoldoCatania/eDMA")
\end{CodeInput}
\end{CodeChunk}

Once the package is correctly installed and loaded, the user faces one function named \code{DMA()} to perform DMA over a series of possible models. The \code{DMA()} function accepts a series of arguments and returns an object of the class \code{DMA} which comes with several methods, see Section \ref{sec:methods}. The arguments the \code{DMA()} function accepts are:
%
\begin{itemize}
  \item \code{formula}: An object of class \code{formula} (or one that can be coerced to that class): A symbolic description of the model to be fitted. The formula should include all the predictors one chooses to include. The inclusion of the constant term follows the usual \proglang{R} practice, \emph{i.e.}, it is included by default and can be removed if necessary. For instance, in order to model \code{y ~ x}, however, without the constant, we can write for example, \code{y ~ x - 1}, see \code{help(formula)}. This implementation follows the common practice for \proglang{R} users, see \emph{e.g.}, the \pkg{plm} package of \cite{plm}.
  %
  \item \code{data}: A \code{data.frame} (or object coercible by \code{as.data.frame()} to a \code{data.frame}) containing the variables in the model. If \code{data} is an object of the class \code{ts}, \code{zoo} or \code{xts}, then the time information is used in the graphical representation of the results as well as for the estimated quantities. The dimension of \code{data} is $T\times \left(1 + n\right)$, containing at each row, the dependent variables $y_t$ and the predictors $\bF_t$, that is $\left(y_t, \bF_t^\prime\right)$, for all $t = 1,\dots, T$.
  %
  \item \code{vDelta}: A $d\times 1$ numeric vector representing a grid of $\delta$. Typically we choose the following grid: $\{0.90,~0.91,~\dots,~1.00\}$. By default \code{vDelta = c(0.90, 0.95, 0.99)}.
  %
  \item \code{dAlpha}: A numeric variable representing $\alpha$ in Equation~\ref{eq:3.4}. By default \code{dAlpha = 0.99}.
  %
  \item \code{vKeep}: A numeric vector of indices representing the predictors that must be always included in the models. The models that do not include the variables declared in \code{vKeep} are automatically discarded. The indices must be consistent with the model description given in \code{formula}. For instance, if the first and fourth variables always have to be included, then we must set vKeep=c(1, 4). Notice that, the intercept (if not removed from \code{formula}) is always in the first position. \code{vKeep} can also be a character vector indicating the names of the predictors if these are consistent with the provided \code{formula}. Furthermore, if \code{vKeep = "KS"} the \qmo Kitchen Sink\qmcsp formulation is adopted, \emph{i.e.}, all the predictors are always included, see, \emph{e.g.}, \cite{Paye.2012}. By default all the combinations are considered, \code{vKeep = NULL}.
      %
  \item \code{bZellnerPrior}: A boolean variable indicating whether the Zellner's prior \citep[see][]{dangl_halling.2012} should be used for the coefficients at time $t=0$. By default \code{bZellnerPrior = FALSE}.
  %
  \item \code{dG}: A numeric variable equal to $100$ by default. If \code{bZellnerPrior = TRUE}, this represents the prior hyperparameter value, $g$, in Equation 4 of \cite{dangl_halling.2012}. Otherwise, if \code{bZellnerPrior = FALSE}, it represents the scaling factor for the variance covariance matrix of the Normal prior for $\btheta_0$, \emph{i.e.}, $\btheta_0\sim N(0,dG\times \mathbb{I})$, where $\mathbb{I}$ is the identity matrix. We generally recommend practitioners to use the default prior, especially in the context of quarterly data, where we typically have $200$ to $300$ observations. For longer time--series, the differences between priors is of relatively less importance.
%
  \item \code{bParallelize}: A boolean variable indicating wether to use multiple processors to speed up the computations. By default \code{bParallelize = TRUE}. Since the use of multiple processors is basically effortless for the user, we suggest to not change this value. Furthermore, if the hardware does not permit parallel computations, the program will automatically adapt to run on a single core.
  %
  \item \code{iCores}: An integer indicating the number of cores to use if \code{bParallelize = TRUE}. By default, all but one cores are used. The number of cores is guessed using the
        \code{detectCores()} function from the \pkg{parallel} package. The choice of the number of cores depends somehow from the specific application, namely the length of the time--series $T$ and the number of the predictors $n$. However, as detailed in \cite{chapman_etal.2008}, the level of parallelization of the code should be traded off with the increase in computational time due to threads communications. Consequently, the user can fine tune its application depending on its hardware changing this parameter.
\end{itemize}

The \code{DMA()} function returns an object of the \textit{formal} class \code{DMA}.\footnote{see, \code{help("class")} and \code{help("DMA-class")}.} This object contains model information and the estimated quantities. It is organized in three slots: \code{model},  \code{Est}, \code{data}. The slot, \code{model}, contains information about the specification used to perform DMA. Examples are: The number of considered models and the computational time in seconds. The slot, \code{Est}, contains the estimated quantities such as: Point forecasts, Predictive likelihood, Posterior inclusion probabilities of the predictors, Filtered estimates of the regression coefficients, $\btheta_t$, and so on. Finally, the slot, \code{data}, includes the data passed to the \code{DMA()} function, organised in the vector of responses \code{vY} and a design matrix \code{mF}.
%For reproducibility purposes, the data that we use in the empirical application in this paper are included in the \pkg{eDMA} package, see Section \ref{sec:EmpiricalApplication} for more details.

\subsection[Using eDMA]{Using \pkg{eDMA}}

After having installed \pkg{eDMA}, it can be easily loaded using:

\begin{CodeChunk}
\begin{CodeInput}
R> library("eDMA")
\end{CodeInput}
\end{CodeChunk}

Thereafter, model estimation can be performed using the R commands reported below.

In order to illustrate how \pkg{eDMA} works in practice, we present the following guidelines using simulated data. We also provide an application using quarterly inflation data in Section \ref{sec:EmpiricalApplication}.

We simulate a time--series of $T=500$ observations from

\begin{equation}\label{eq:dgp}
  y_t = \bF_t^\prime\btheta_t + \sqrt{0.1}\varepsilon_t,\quad \varepsilon_t\overset{iid}{\sim}\mathcal{N}\left(0, 1\right).
\end{equation}

The first four elements of $\btheta_t$ vary according to random--walks, whereas the remaining elements in $\btheta_t$ are equal to zero at all time periods. In other words, $\btheta_t = \left(\theta_{1,t}, \theta_{2,t}, \theta_{3,t}, \theta_{4,t}, \theta_{5,t},\theta_{6,t}\right)^\prime$ with

\begin{equation}\label{eq:theta_rw}
  \theta_{k,t} = \theta_{k,t-1} + \sqrt{0.01}\eta_{k,t},\quad\eta_{k,t}\overset{iid}{\sim}\mathcal{N}\left(0, 1\right),
\end{equation}

for $k = 1,2,3,4$, and $\eta_{k,t}\indep\eta_{j,t}$, for all $k\neq j$. The last two elements of $\btheta_t$ are equal to zero, that is, $\theta_{5,t}=\theta_{6,t}=0$ for $t=1,\dots, T$. The first element of the $6\times 1$ vector, $\bF_t$, is one, representing the constant term. The remaining elements are generated from a standard Gaussian distribution, \emph{i.e.}, $\bF_t = \left(1.0, x_{2,t}, x_{3,t}, x_{4,t}, x_{5,t}, x_{6,t}\right)^\prime$, where $x_{k,t}\overset{iid}{\sim}\mathcal{N}\left(0, 1\right)$ and $x_{k,t}\indep x_{j,t}$ for all $k\neq j$.
We simulate the data in the following way (that is $\theta_{5,t}=\theta_{6,t}=0$) to illustrate that DMA is indeed able to identify the correct variables. In other words, the inclusion probabilities of the last two predictors ought to be zero as they do not impact $y_t$ through $\bF_t$. Conversely, inclusion probabilities of the first four predictors ought to converge to 1.

This data is simulated using the \code{SimulateDLM()} function available in \pkg{eDMA}, details are reported in the \proglang{R} documentation, see \code{help("SimulateDLM")}. We organize the data in a \code{data.frame} named \code{SimData}, which is included in \pkg{eDMA} and can be loaded into the workspace by executing:

\begin{CodeChunk}
\begin{CodeInput}
R> data("SimData", package = "eDMA")
\end{CodeInput}
\end{CodeChunk}

DMA is then performed using the function \code{DMA()} as:

\begin{CodeChunk}
\begin{CodeInput}
R> Fit <- DMA(y ~ x2 + x3 + x4 + x5 + x6 , data = SimData,
          vDelta = seq(0.9, 1.0, 0.01))
\end{CodeInput}
\end{CodeChunk}

Information on the DMA procedure is available by typing:

\begin{CodeChunk}
\begin{CodeInput}
R> Fit

------------------------------------------
-        Dynamic Model Ageraging         -
------------------------------------------

Model Specification	
T     = 500
n     = 6
d     = 11
Alpha = 0.99
Model combinations = 63
Model combinations including averaging over delta = 693
------------------------------------------
Prior : Multivariate Gaussian with mean vector 0
        and covariance matrix equal to: 100 x diag(6)
------------------------------------------
The grid for delta:

Delta =  0.90, 0.91, 0.92, 0.93, 0.94, 0.95,
         0.96, 0.97, 0.98, 0.99, 1.00
------------------------------------------

Elapsed time	: 0.57 secs
\end{CodeInput}
\end{CodeChunk}

Note, we specify a grid of eleven equally spaced values for $\delta$ ($d=11$) ranging from $0.90$ to $1.00$. Furthermore, since we do not specify any value for \code{bZellnerPrior} and \code{bParallelize}, their default values, \code{bZellnerPrior = FALSE} and \code{bParallelize = TRUE} have been used.

In order to extract the quantities estimated by DMA, the user can relay on the \code{as.data.frame()} method. \code{as.data.frame()} accepts two arguments: (i): An object of the class \code{DMA} and (ii): A character string, \code{which}, indicating the quantity to extract. Possible values for \code{which} are:
\begin{itemize}
  \item \code{"vyhat"}: Point forecasts of DMA, see Equation~\ref{eq:4.5}. \code{"vyhat_DMS"} for point forecast according to DMS.
  \item \code{"mincpmt"}: Posterior inclusion probabilities of the predictors at each point in time, see \cite{koop_korobilis.2012} for more details.
  \item \code{"vsize"}: Expected number of predictors (average size), see \cite{koop_korobilis.2012} and point (ii) at page 7.
  \item \code{"vsize_DMS"}: Number of predictors in the model with the highest posterior model probability, at each point in time, see Equation~\ref{eq:3.2}.
  \item \code{"mtheta"}: Filtered estimates of the regression coefficients for DMA, see Equation~\ref{eq:theta_post}.
  \item \code{"mpmt"}: Posterior probability of the degrees of instability, see Equation~\ref{eq:delta}.
  \item \code{"vLpdfhat"}: Predictive (log--)likelihood of DMA, see \cite{dangl_halling.2012}.
  \item \code{"vLpdfhat_DMS"}: Predictive (log--)likelihood according to DMS.
  \item \code{"vdeltahat"}: Posterior weighted average of $\delta$, that is $\hat\delta_t = \sum_{j=1}^{d}\delta_j p\left(\delta_j\vert\mathcal{F}_t\right)$.
  \item \code{"mvdec"}: Individual components of Equation~\ref{eq:varDec}, see point (vi) in page 8 and \cite{dangl_halling.2012} for more details. The function returns a $T\times 5$ matrix whose columns contain the variables.
  \begin{itemize}
    \item \code{vobs}: Observational variance, Obs.
    \item \code{vcoeff}: Variance due to errors in the estimation of the coefficients, Coeff.
    \item \code{vmod}: Variance due to model uncertainty, Mod.
    \item \code{vtvp}: Variance due to uncertainty with respect to the choice of the degrees of time--variation in the regression coefficients, TVP.
       \item \code{vtotal}: Total variance, that is \code{vtotal} = \code{vobs} + \code{vcoeff} + \code{vmod} + \code{vtvp} .
  \end{itemize}
  \item \code{"vhighmp_DMS"}: Highest posterior model probability, \emph{i.e.}, $\underset{i}{\max}~P\left(M_i\vert\mathcal{F}_t\right),\quad t =1,\dots,T$.
  \item \code{"vhighmpTop01_DMS"}: Sum of the $10\%$ highest posterior model probabilities.

\end{itemize}

The additional \code{numeric} argument, \code{iBurnPeriod}, determines the length of the burn--in period, \emph{i.e.}, results before $t$=\code{iBurnPeriod} are discarded. By default, \code{iBurnPeriod = NULL}, meaning that no burn--in period is considered. For instance, in order to extract the posterior inclusion probabilities of the predictors, with a burn--in period of 50 observations, we can easily run the following command

\begin{CodeChunk}
\begin{CodeInput}
R> PostProb = as.data.frame(Fit, which = "mincpmt", iBurnPeriod = 50)
\end{CodeInput}
\end{CodeChunk}

which returns a $(T-$\code{iBurnPeriod}$)\times 7$ matrix of inclusion probabilities for the predictors at each point in time. Final values of \code{PostProb} are printed as:

\begin{CodeChunk}
\begin{CodeInput}
R> round(tail(PostProb), 2)
\end{CodeInput}
\begin{CodeOutput}
       (Intercept) x2 x3 x4   x5   x6
[445,]           1  1  1  1 0.06 0.03
[446,]           1  1  1  1 0.06 0.03
[447,]           1  1  1  1 0.07 0.03
[448,]           1  1  1  1 0.07 0.03
[449,]           1  1  1  1 0.07 0.03
[450,]           1  1  1  1 0.09 0.04
\end{CodeOutput}
\end{CodeChunk}

Furthermore, if the supplied data is a \code{ts}, \code{zoo} or \code{xts} object, the class membership is automatically transferred to the output of the \code{as.data.frame()} method.

The \code{plot()} method is also available for the class \code{DMA}. Specifically, this method prints an interactive menu in the console permitting the user to chose between a series of interesting graphical representation of the estimated quantities. It can be straightforwardly executed running:

\begin{CodeChunk}
\begin{CodeInput}
R> plot(Fit)
\end{CodeInput}
\begin{CodeOutput}
Print 1-11 or 0 to exit
 1: Point forecasts
 2: Predictive likelihood
 3: Posterior weighted average of delta
 4: Posterior inclusion probabilities of the predictors
 5: Posterior probability of the degree of instability
 6: Filtered estimates of the regression coefficients, theta
 7: Observational variance
 8: Variance due to errors in the estimation of the coefficients
 9: Variance due to model uncertainty
10: Variance due to uncertainty with respect to the choice of
    the degrees of time-variation in the regression coefficients
11: Expected number of predictors (average size)
12: Number of predictors (highest posterior model probability) (DMS)
13: Highest posterior model probability (DMS)
14: Point forecasts (highest posterior model probability) (DMS)
15: Predictive likelihood (highest posterior model probability) (DMS)
\end{CodeOutput}
\end{CodeChunk}
%
and selecting the desiderated options. The additional character argument, \code{which}, can be supplied in order to directly plot one particular quantity. Possible values for \code{which} are the same of the \code{as.data.frame()} method. Similar to \code{as.data.frame()}, the additional \code{numeric} argument \code{iBurnPeriod} determines the length of the burn--in period. Typically, it takes around $30$ to $50$ for the model to adapt to the time--series given the prior. Therefore, in almost all applications, the first $30$ to $50$ observations are discarded.

The code:

\begin{CodeChunk}
\begin{CodeInput}
R> plot(Fit, which = "mincpmt", iBurnPeriod = 50)
\end{CodeInput}
\end{CodeChunk}

plots the inclusion probabilities for the predictors discarding the first 50 observations. The outcome is reported in Figure~\ref{fig:inclusionprob_simulation}. As expected, $x_1$ to $x_4$ quickly converge to $1$ after few observations. Conversely, the inclusion probabilities of the last two predictors with loading factor equal to zero, quickly converge to $0$.

\begin{figure}[!t]
\centering
\includegraphics[width=1\textwidth]{Simulation_InclusionProbabilities.pdf}
\caption{Posterior inclusion probabilities of the predictors using simulated data.}
%%%
\label{fig:inclusionprob_simulation}
\end{figure}

\subsection[Additional methods for the DMA class]{Additional methods for the \code{DMA} class}
\label{sec:methods}

The \code{DMA} class comes with several methods for extracting and representing estimated quantities. The \code{plot()}, \code{as.data.frame()} and \code{show()} methods have been previously introduced, additional methods are: \code{summary()}, \code{coef()}, \code{residuals()}, \code{inclusion.prob()}, and \code{pred.like()}.

For instance, the \code{summary} method prints a summary of the estimated model directly in the console. The code:

\begin{CodeChunk}
\begin{CodeInput}
R> summary(Fit, iBurnPeriod = 50)
\end{CodeInput}
\end{CodeChunk}

produces the output:

\begin{CodeChunk}
\begin{CodeInput}
Call:
 DMA(formula =  vY ~ x1 + x2 + x3 + x4 + x5 )

Residuals:
    Min      1Q  Median      3Q     Max
-2.0408 -0.3833  0.0421  0.4427  2.3140

Coefficients:
            E[theta_t] SD[theta_t] E[P(theta_t)] SD[P(theta_t)]
(Intercept)       0.51        0.68          1.00           0.00
x1               -0.64        0.65          0.90           0.29
x2                2.11        1.74          0.93           0.20
x3               -1.43        1.02          0.99           0.02
x4                0.01        0.03          0.07           0.07
x5                0.00        0.01          0.06           0.04

Variance contribution (in percentage points):
  vobs vcoeff   vmod   vtvp
  3.18   1.73   0.08  95.01

Top 10% included predictors:	(Intercept)

Forecast Performance:
                         DMA      DMS
MSE                    0.486    0.482
MAD                    0.537    0.532
Predictive Likehood -462.360 -463.083
\end{CodeInput}
\end{CodeChunk}

where the quantities, \code{E[theta_t], SD[theta_t], E[P(theta_t)]} and \code{SD[P(theta_t)]} represent the means and standard deviations across the time dimension of the filtered estimates of $\btheta_t^{(i)}$, and the inclusion probabilities after burn-in.

The last part of the summary, (\code{Forecast Performance}), prints the output of the \code{BacktestDMA()} function implemented in \pkg{eDMA}. \code{BacktestDMA()} accepts a \code{DMA} object and returns a \code{matrix} with out--of--sample Mean Squared Error (MSE), Mean Absolute Deviation (MAD) and Predictive Likelihood, computed according to DMA and DMS, see \code{help("BacktestDMA")}.

The additional methods: \code{coef()}, \code{residuals()}, \code{inclusion.prob()}, and \code{pred.like()} are wrapper to the \code{as.data.frame()} method and focus on particular estimated quantities, for instance:

\begin{itemize}
  \item[-] \code{coef()}: Returns a $T\times n$ \code{matrix} with the filtered regressor coefficients, $\btheta_t,\quad t = 1,\dots,T$.
  \item[-] \code{residuals()}: Extract the residuals of the model, \emph{i.e.}, $y_t - \hat{y_t},\quad t = 1,\dots, T$. The additional \code{boolean} argument \code{standardize} controls if the standardize residuals should be returned. By default \code{standardize = FALSE}. The additional \code{character} argument, \code{Type}, permits to choose between residuals evaluated using DMA (\code{"DMA"}) or DMS (\code{"DMS"}). By default \code{Type = "DMA"}.
  \item[-] \code{inclusion.prob()}: Extract the inclusion probabilities of the predictors. Analogous to \code{as.data.frame(object, which = "mincpmt", iBurnPeriod)}.
  \item[-] \code{pred.like()}: Extract the predictive log likelihood series. The additional argument \code{Type} permits to choose between predictive likelihoods evaluated using DMA and DMS. By default \code{Type = "DMA"}. Similar to the above variables, \code{pred.like()} accepts \code{iBurnPeriod}.
\end{itemize}

\section[Computational challenges]{Computational challenges}
%
\label{sec:computational}
%
Although estimation of DMA does not require resorting to simulation
methods, in many economic applications, performing DMA can become computationally
very cumbersome.
%
%
As it can be seen from the set of recursions from
the Section \ref{sec:ModifiedDMA}, DMA consists of considering a large number of model combinations. In many cases, DMA tends to occupy a large chunk of Random--Access Memory (RAM). Often on a standard PC, the system basically runs out of memory due to the large number of combinations and the amount of information that must be saved. Therefore, it limits the use of DMA to middle--sized data sets. For instance, in their seminal paper, \cite{koop_korobilis.2012} use DMA to forecast quarterly inflation. Thus, $y_{t}$ in Equation~\ref{eq:2.1} is the percentage changes in the quarterly U.S.
GDP price deflator and $\bF_{t}$ consists of $14$ exogenous predictors and three lags of $y_{t}$ for a total of $17$ variables. However, handling $2^{17}$
combinations reveals to be very burdensome in their programming framework. Therefore, \cite{koop_korobilis.2012}
choose to include three lags of inflation in all model combinations and thus reduce the model space to $2^{14}$ model combinations.

\begin{figure}
\centering
\includegraphics[width=1\textwidth]{ComputationalTime.pdf}
\caption{Computational time for \code{DMA()} using simulated data. Each panel represents computation time in minutes for \code{DMA()} using different sample sizes, $T$, number of predictors, $n$, and values of $d$, the number of points in the grid of $\delta$. The values for $d$ range between 2 and 10, the line at the bottom of each subfigure is for $d=2$, the one immediately above is $d=3$ and so on until the last which is for $d=10$. Computations are performed on a standard Intel Core i7--4790 processor with 8 threads and 8 GB of RAM with Ubuntu 12.04 server edition.}
%%%
\label{fig:comptime}
%\end{sidewaysfigure}
\end{figure}

We can argue that DMA can impose a very substantial challenge for the practitioner when dealing with a large number of predictors,
namely that, besides dealing with the task of transforming mathematical
equations from paper to codes, handling data and estimation
issues, a practitioners also has to overcome \qmo technical/computer science\qmcsp challenges
such as how to deal with extensive memory consumption and how to use multiple cores instead of a single core to speed up computation time.
Although one can always improve the computational procedure by \qmo coding smarter\qmcsp or discovering ways to optimize memory allocation,
it seems unreasonable to expect that practitioners in economics
should have extensive knowledge of computer science concepts
such as those stated above.

In this paper, we provide practical solutions to these problems. First, reduction in computation time is implemented by writing all the code in \proglang{C++} using the \code{armadillo} library of \cite{sanderson.2010}.
%
Second, we  exploit multiple processors through the \code{OpenMP} API whenever the hardware is suited for that.
%
The combination of \proglang{C++} routines and parallel processing permits to dramatically speed up the computations over the same code written in plain \proglang{R}.

In order to provide an intuitive example of the advantages of our package, we report a comparison between our code and the available \pkg{dma} package of \cite{dma}. For this experiment, since the \pkg{dma} package cannot operate over a grid value of $\delta$, we fix $\delta$ at $0.95$. We simulate $T=\{100,500,1000\}$ observations from a DLM with $n = \{4,6,8,10,12,14,16\}$ predictors and evaluate the differences in the computational time of the \code{dma()} function in the \pkg{dma} package and the \code{DMA()} function in the presented \pkg{eDMA} package. The experiment is performed on a standard Intel Core i7--4790 processor with 8 threads and Ubuntu 12.04 server edition.

Table~\ref{tab:RafComparison} reports the ratio of the CPU time for different values of $T$ and $n$ between \code{dma()} and \code{DMA()}. As one can note, the decrease in computational time in favor of our package is huge. For example, in the case $T = 500$ and $n = 16$, \code{dma()} takes 37.5 minutes while \code{DMA()} only 1.8.
%
It is also worth stressing that, the benefit of using \pkg{eDMA} does not only concern the possibility of running moderately large applications in a reasonable time using a commercial hardware, but also enables practitioners to run application with a large number of exogenous variables. To give an idea of the computational time a \pkg{eDMA} user faces, we report a second simulation study. We simulate from a DLM with $T = \{100,200,...,900,1000\}$, $n = \{2,3,..,18\}$ and run \code{DMA()} using a grid of values for $\delta$ between $0.9$ and $1.0$ with different spaces $d$, namely $d=\{2,3,\dots,10\}$. Figure~\ref{fig:comptime} displays the computational time in minutes for all the combination of $T,n,d$. The lines reported in each subfigure represent the computational time for a specific choice of $d$. The line at the bottom of each subfigure is for $d=2$,\footnote{In this case $\delta$ can takes values $\delta=0.9$ and $\delta = 1.0$.} the one immediately above is for $d=3$ and so on until $d=10$. From the Figure, we can see that, when $T\le400$, even for $n=18$ and $d=10$, the computational time is less then 15 minutes. Such sample sizes are relatively common in economic applications. When $T$ increases, computational time increases linearly. For example, when $T=800, n = 18$ and $d=10$, computational time is 30 minutes, which is the double of the same case with $T=400$.

The other relevant problem with DMA is the RAM usage. Specifically, if we want to store the quantities defined in Equations \ref{eq:2.1} and \ref{eq:3.3}, we need to define two arrays of dimension $T\times d\times\left(2^n-1\right)$. These kind of objects are not present in the \pkg{eDMA} package since we rely on the markovian nature of the model clearly evident from Equation~\ref{eq:2.1}. In this respect, we keep track of the quantities coming from Equation~\ref{eq:3.3} and $p\left(y_{t}\vert M_{i},\delta_{j},\mathcal{F}_{t-1}\right)$ only for two consecutive periods during the loop over $T$. %However, arrays of size $2\times d\times\left(2^n-1\right)$ cannot be eliminated at all in our framework.
RAM usage is still efficiently performed in the \pkg{eDMA} package. Indeed, the computer where we run all our simulations has only $8$GB of RAM.

\begin{table}[!t]
\centering
%\resizebox{1\columnwidth}{!}{%
\setlength\tabcolsep{17pt}
\begin{tabular}{lccccccc}
\toprule
$T/n$ & 4 & 6 & 8 & 10 & 12 & 14 & 16\\
\cmidrule(lr){1-1}\cmidrule(lr){2-8}
100 & 10.9 & 92.6 & 133.2 & 88.4 & 69.4 & 70.5 & 58.4 \\
500 & 37.5 & 29.4 & 31.5 & 30.7 & 25.7 & 26.5 & 25.4 \\
1000 & 13.0 & 15.0 & 13.8 & 12.9 & 12.7 & 13.5 & 13.8 \\
\bottomrule
\end{tabular}
%}
\caption{\footnotesize{Ratio of computation time between the \code{dma()} function from the \pkg{dma} package of \cite{dma} and the \code{DMA()} function of the \pkg{eDMA} package for different values of $T$ and $n$.}}
\label{tab:RafComparison}
\end{table}

\section[A DMA example: Inflation data]{A DMA example: Inflation data}
\label{sec:EmpiricalApplication}

This section shows how to run DMA using the \proglang{R} \citep{R.2015} package
\pkg{DMA} \citep{eDMA}. We use a time--series of quarterly U.S. inflation rate for illustration, show how to obtain
posterior output and briefly interpret results. The example can be thought
of as a typical assignment for a researcher at a central bank who
is interested in forecasting inflation several--quarters ahead and understand the relationship between inflation and business cycles.

\subsection[Data]{Data}

In this analysis, we use the data of \cite{groen_etal.2013}.\footnote{The data is downloadable from \url{http://www.tandfonline.com/doi/suppl/10.1080/07350015.2012.727718}.}
Particularly, as a measure of inflation, $y_{t}$, we consider quarterly
log changes in the Gross Domestic Product implicit price deflator (GDPDEF) observed
at a quarterly frequency ranging from 1960q1 to 2011q2. The number
of exogenous predictors are fifteen. This number is in accordance
with typical \qmo real--world\qmcsp application, see \cite{dangl_halling.2012}, \cite{koop_korobilis.2012} and \cite{groen_etal.2013}.

The predictors are: Real GDP in volume terms (ROUTP), real
durable personal consumption expenditures in volume terms (RCONS), real residential investment in
volume terms (RINVR), the import deflator (PIMP), the unemployment
ratio (UNEMP), non--farm payrolls data on employment (NFPR), housing
starts (HSTS), the real spot price of oil (OIL), the real food commodities
price index (FOOD) the real raw material commodities price index
(RAW), and the M2 monetary aggregate (M2), which can reflect information
on the current stance of monetary policy and liquidity in the economy
as well as spending in households. In addition, we also use data on
the term structure of interest rates approximated by means of:
The level factor (YL), the slope factor (TS) and curvature factor
(CS). Finally, we proxy inflation expectations through the one--year
ahead inflation expectations that come from the Reuters/Michigan Survey
of Consumers (MS). We include the data in the \pkg{eDMA} package as a \code{xts} object of dimension $206\times 16$ named \code{USData}.

For most series, we use the percentage change of the original series
in order to remove possible stochastic and deterministic trends. Exceptions
are HSTS, for which we use the logarithm of the respective levels,
as well as UNEMP, YL,TS, CS and MS, where we use the \qmo raw\qmcsp levels,
see also \cite{groen_etal.2013} for more details. Finally, since inflation
very persistence, besides these $15$ predictors, we follow \cite{groen_etal.2013} and also include four inflation lags,
$y_{t-1},...,y_{t-4}$, as predictors. In \pkg{eDMA}, we implement the function, \code{Lag()}, which allows us to lag variables delivered in the form of vector or matrices. For instance, to lag the \code{numeric} vector \code{GDPDEF} of length $T$ by one period, we simply run

\begin{CodeChunk}
\begin{CodeInput}
R> Lag(GDPDEF, 1)
\end{CodeInput}
\end{CodeChunk}

which returns a \code{numeric} vector of length $T$ containing the lagged inflation. Values that are not available (\emph{i.e.}, the value $y_0$ in this example) are replaced by \code{NA}.

\subsection[Model estimation]{Model estimation}

We have a total of $2^{19} = 524288$ model combinations.\footnote{The model which include only the constant is also considered. Indeed, when \code{vKeep = NULL}, the number of models is $2^n - 1$, however, when \code{vKeep != NULL}, the number of models is $2^b$, where \code{b = n - length(vKeep)}.} Furthermore, we let $\delta=\left\{ 0.9,0.91,...,1\right\} $ such
that we have a total of $\left(2^{19}\right)\cdot11 = 5767168$ combinations. We
set $\alpha=0.99$ and specify a noninformative prior over the model
combinations, $p\left(M_{s}\mid\mathcal{F}_{0}\right)=1/\left(d\cdot k\right),\quad s = 1,\dots,d\cdot k$,
such that initially, all models are equally likely. We then update these model
probabilities as new information arrives. As previously mentioned,
we include the constant in all models, see also \cite{groen_etal.2013}.

In terms of implementation, we start by loading the \pkg{eDMA} package and the data set by typing:

\begin{CodeChunk}
\begin{CodeInput}
R> library("eDMA")
R> data("USData")
\end{CodeInput}
\end{CodeChunk}

In order to perform DMA using the \code{DMA()} function, we write:

\begin{CodeChunk}
\begin{CodeInput}
R> Fit <- DMA(GDPDEF ~  Lag(GDPDEF, 1) + Lag(GDPDEF, 2) +
                        Lag(GDPDEF, 3) + Lag(GDPDEF, 4) +
                        Lag(ROUTP, 1)  + Lag(RCONS, 1)  +
                        Lag(RINVR, 1)  + Lag(PIMP, 1)   +
                        Lag(UNEMP, 1)  + Lag(NFPR, 1)   +
                        Lag(HSTS, 1)   + Lag(M2, 1)     +
                        Lag(OIL, 1)    + Lag(RAW, 1)    +
                        Lag(FOOD, 1)   + Lag(YL, 1)     +
                        Lag(TS, 1)     + Lag(CS, 1)     +
                        Lag(MS, 1), data = USData,
                        vDelta = seq(0.90, 1.00, 0.01), vKeep = 1)
\end{CodeInput}
\end{CodeChunk}

In terms of implementation, we suggest using the non--informative prior,
\code{bZellnerPrior = FALSE}, which is the default, see Section \ref{sec:sensitivity} for a prior sensitivity analysis. More details on the model can be made available by typing \code{Fit}:

\begin{CodeChunk}
\begin{CodeInput}
R> Fit
\end{CodeInput}
\begin{CodeOutput}
------------------------------------------
-        Dynamic Model Ageraging         -
------------------------------------------

Model Specification	
T     = 202
n     = 20
d     = 11
Alpha = 0.99
Model combinations = 524288
Model combinations including averaging over delta = 5767168
------------------------------------------
Prior : Multivariate Gaussian with mean vector 0
        and covariance matrix equal to: 100 x diag(20)

Variables always included : (Intercept)
------------------------------------------
The grid for delta:

Delta =  0.90, 0.91, 0.92, 0.93, 0.94, 0.95,
         0.96, 0.97, 0.98, 0.99, 1.00
------------------------------------------

Elapsed time : 1686.14 secs
\end{CodeOutput}
\end{CodeChunk}

As it can be seen, the total estimation time of our DMA is 1686.14 seconds corresponding to around
28 minutes on an Intel Core i7-3630QM processor. A complete summary of the estimation is available as:

\begin{CodeChunk}
\begin{CodeInput}
R> summary(Fit, iBurnPeriod = 32)
\end{CodeInput}
\begin{CodeOutput}
Call:
 DMA(formula =  Lag(GDPDEF, 1) + Lag(GDPDEF, 2) +
                Lag(GDPDEF, 3) + Lag(GDPDEF, 4) +
                Lag(ROUTP, 1)  + Lag(RCONS, 1)  +
                Lag(RINVR, 1)  + Lag(PIMP, 1)   +
                Lag(UNEMP, 1)  + Lag(NFPR, 1)   +
                Lag(HSTS, 1)   + Lag(M2, 1)     +
                Lag(OIL, 1)    + Lag(RAW, 1)    +
                Lag(FOOD, 1)   + Lag(YL, 1)     +
                Lag(TS, 1)     + Lag(CS, 1)     +
                Lag(MS, 1) )

Residuals:
    Min      1Q  Median      3Q     Max
-1.3494 -0.2990 -0.0138  0.2218  1.5653

Coefficients:
               E[theta_t] SD[theta_t] E[P(theta_t)] SD[P(theta_t)]
(Intercept)          0.11        0.16          1.00           0.00
Lag(GDPDEF, 1)       0.41        0.18          0.83           0.30
Lag(GDPDEF, 2)       0.02        0.02          0.19           0.11
Lag(GDPDEF, 3)       0.09        0.06          0.38           0.22
Lag(GDPDEF, 4)       0.12        0.06          0.52           0.23
Lag(ROUTP, 1)        0.00        0.01          0.15           0.09
Lag(RCONS, 1)        0.00        0.00          0.14           0.07
Lag(RINVR, 1)        0.02        0.02          0.21           0.11
Lag(PIMP, 1)         0.21        0.09          0.83           0.27
Lag(UNEMP, 1)       -0.03        0.06          0.22           0.12
Lag(NFPR, 1)         0.02        0.01          0.22           0.14
Lag(HSTS, 1)         0.02        0.03          0.20           0.09
Lag(M2, 1)           0.01        0.01          0.16           0.06
Lag(OIL, 1)         -0.03        0.09          0.36           0.22
Lag(RAW, 1)          0.00        0.01          0.16           0.07
Lag(FOOD, 1)         0.01        0.01          0.20           0.12
Lag(YL, 1)           0.27        0.44          0.36           0.32
Lag(TS, 1)           0.01        0.04          0.15           0.07
Lag(CS, 1)          -0.03        0.07          0.19           0.12
Lag(MS, 1)           0.02        0.03          0.18           0.09

Variance contribution (in percentage points):
  vobs vcoeff   vmod   vtvp
 43.59  11.71  12.84  31.86

Top 10% included predictors:	(Intercept), Lag(PIMP, 1)

Forecast Performance:
                         DMA      DMS
MSE                    0.235    0.251
MAD                    0.361    0.369
Predictive Likehood -106.352 -122.943
\end{CodeOutput}
\end{CodeChunk}

Below, we go into more details
with regards to how to use the output from the estimation procedure to interpret our results.

\subsection[Using the output from eDMA]{Using the output from \pkg{eDMA}}

The output can be divided into two main parts: (a): Full--sample
analysis, (b): Out--of--sample analysis. With regards to (a), the most
interesting quantities are:~\code{mincpmt},~\code{vsize}, ~\code{mtheta},~
\code{vdeltahat}, and \code{mvdec}, see Section \ref{sec:ThePackage}.

For instance, the inclusion probabilities of the predictors for the last part of the sample can be
printed by:

\begin{CodeChunk}
\begin{CodeInput}
R> InclusionProb <- inclusion.prob(Fit, iBurnPeriod = 32)
R> tail(round(InclusionProb[, 1:4], 2))
\end{CodeInput}
\begin{CodeOutput}
           (Intercept) Lag(GDPDEF, 1) Lag(GDPDEF, 2) Lag(GDPDEF, 3)
2010-01-01           1           0.99           0.43           0.71
2010-04-01           1           1.00           0.44           0.72
2010-07-01           1           1.00           0.45           0.73
2010-10-01           1           1.00           0.45           0.73
2011-01-01           1           1.00           0.45           0.73
2011-04-01           1           1.00           0.45           0.73
\end{CodeOutput}
\end{CodeChunk}

\begin{figure}[!t]
\centering
\includegraphics[width=1\textwidth]{InclusionProbabilities.pdf}
\caption{Posterior inclusion probabilities for the most important predictors of DMA. Panels (a), (b) and (c): First, third and fourth lags of inflation. Panel (d): Import deflator (PIMP). Panel (e): Inflation expectations (MS). Panel (f): M2 monetary aggregate (M2). Panel (g): Real spot price of oil (OIL). Panel (h): Level factor of the term structure (YL). We refer the reader to \cite{groen_etal.2013} for more details regarding the variables. The gray vertical bars indicate business cycle peaks, \emph{i.e.}, the point at which an economic expansion transitions to a recession, based on National Bureau of Economic Research (NBER) business cycle dating.}
%%%
\label{fig:incprob}
\end{figure}

The above matrix shows the inclusion probabilities of: The constant and
$y_{t-1},...,y_{t-3}$, from $2010q1$ to $2011q2$. Notice that, the inclusion probabilities of the
constant term, \code{(Intercept)}, are always equal to $1$ as every model contains this
term (since we set \code{vKeep = 1}), see (iii) in page 7 of this paper.

\begin{figure}[!t]
\centering
\includegraphics[width=1\textwidth]{FilteredBeta.pdf}
\caption{Filtered estimates of the regression coefficients for the most important predictors of DMA. Panels (a), (b) and (c): First, third and fourth lags of inflation. Panel (d): Import deflator (PIMP). Panel (e): Inflation expectations (MS). Panel (f): M2 monetary aggregate (M2). Panel (g): Real spot price of oil (OIL). Panel (h): Level factor for the terms structure (YL). We refer the reader to \cite{groen_etal.2013} for more details regarding the variables. The gray vertical bars indicate business cycle peaks, \emph{i.e.}, the point at which an economic expansion transitions to a recession, based on National Bureau of Economic Research (NBER) business cycle dating.}
%%%
\label{fig:beta}
\end{figure}

In Figure~\ref{fig:incprob}, we report the inclusion probabilities for predictors
that are important at least one point in time. To be precise, any
predictor where the inclusion probability is never above 0.2 is excluded.
In these plots, we also make evident NBER recorded recessions
(shaded gray bars). Overall, we observe a good amount of time--variation
these plots. The lags of inflation, except for $y_{t-2}$ all seem
important. The import deflator (PIMP) also receives high posterior
probability throughout the sample. Inflation expectation (MS) and
M2 receive higher probabilities towards the end of the sample. However,
the inclusion probability of these predictors at any point in time
is relatively low. Real spot price of oil (OIL) receives high inclusion probability
during the post Great Moderation era, whereas we observe the opposite
trend for YL. In addition to the inclusion probabilities, we also report filtered
estimates of the regression coefficients for these predictors in Figure~\ref{fig:beta}. These quantities are extracted from \code{Fit} simply using

\begin{CodeChunk}
\begin{CodeInput}
R> mTheta <- coef(Fit, iBurnPeriod = 32)
\end{CodeInput}
\end{CodeChunk}

Besides these variables, the output from DMA can be used to analyze:

\begin{figure}[!t]
\centering
\includegraphics[width=1\textwidth]{CombinedPlots.pdf}
\caption{Posterior output for DMA. Panel (a): Posterior weighted average estimate of $\delta$. Panel (b):  Number of predictors for the model with the highest posterior probability. Panel (c):  Sum of top $10\%$ inclusion probabilities. Panel (d): Observational variance. Panel (e): Variance due to errors in the estimation of the coefficients. Panel (f) Variance due to model uncertainty (Mod, \textit{solid}) and variance due to uncertainty with respect to the choice of the degrees of time-variation in the regression coefficients (TVP, \textit{red-dotted}). The gray vertical bars indicate business cycle peaks, \emph{i.e.}, the point at which an economic expansion transitions to a recession, based on National Bureau of Economic Research (NBER) business cycle dating.}
%%%
\label{fig:combined}
\end{figure}

%\begin{itemize}
%\item
The magnitude of time--variation in the model parameters, \code{"vdeltahat"},
which is the posterior weighted average of $\delta$ at each point in
time. We report this estimate in panel (a) of Figure~\ref{fig:combined}. The analogous plot in \proglang{R} can be obtained using:

\begin{CodeChunk}
\begin{CodeInput}
R> plot(Fit, which = "vdeltahat", iBurnPeriod = 32)
\end{CodeInput}
\end{CodeChunk}

Except for the recessions in $1990$ and $2001$, there
is a very intuitive relationship between $\delta$ and recession periods.
Typically, $\delta$ falls during recessions, which fares well with
the notion that $\btheta_{t}$ changes rapidly as the model needs to
adopt to the changes in the data. Conversely, $\delta$ remains high
and close to $1$ during the Great Moderation. We can also use \code{inclusion.prob()} to extract the posterior probability of each value of $\delta$. We can print these quantities for the final part of the sample using:

\begin{CodeChunk}
\begin{CodeInput}
R> InclusionProbDelta <- inclusion.prob(Fit, iBurnPeriod = 32)
R> round(tail(InclusionProbDelta), 2)
\end{CodeInput}
\begin{CodeOutput}
           0.9 0.91 0.92 0.93 0.94 0.95 0.96 0.97 0.98 0.99    1
2010-01-01   0    0    0    0 0.01 0.02 0.06 0.16 0.30 0.28 0.17
2010-04-01   0    0    0    0 0.01 0.02 0.06 0.15 0.29 0.29 0.17
2010-07-01   0    0    0    0 0.01 0.02 0.05 0.15 0.30 0.29 0.18
2010-10-01   0    0    0    0 0.01 0.02 0.06 0.16 0.30 0.29 0.16
2011-01-01   0    0    0    0 0.01 0.02 0.06 0.16 0.30 0.29 0.16
2011-04-01   0    0    0    0 0.01 0.02 0.06 0.17 0.32 0.28 0.14
\end{CodeOutput}
\end{CodeChunk}

where the column names are the values of $\delta$.

%\item
In panel (b), we report the number of predictors contained in the model with the
highest posterior probability , $p\left(M_i\vert \mathcal{F}_t\right)$,at each point in time. This can be achieved by:

\begin{CodeChunk}
\begin{CodeInput}
R> plot(Fit, which = "vsize_DMS", iBurnPeriod = 32)
\end{CodeInput}
\end{CodeChunk}

Alternatively, we can also plot the expected number of predictors at each point in
time replacing \code{which = "vsize_DMS"} by \code{which = "vsize"}. A very interesting result from
panel (b) is that, although we have $19$ predictors, at each point
in time the best model contains only a few predictors. Furthermore,
there is also evidence of time--variation in the number of predictors as we have as many as eight and as few as two predictors over the sample.

%\item
We can also use posterior model probabilities to obtain an idea
of how important is model averaging. In panel (c), we report the sum of
the posterior inclusion probabilities for the $10\%$ of models (\code{which = "vhighmpTop01_DMS"}). If this number is high, then it means that
relatively few models dominate, obtaining relatively high posterior probabilities. Conversely,
if this number is low, then no individual (or group of) models receive
high probabilities, which provides evidence in favor of averaging over predictors instead of selecting.

%\item
Finally, in panels (d), (e) and (f), we report the variance decomposition
analysis (\code{which = "mvdec"}). Evidently, the dominant source of uncertainty is the observational
variance. This is not surprising as random fluctuation are expected to dominate uncertainty.
However, uncertainty about the correct magnitude of time--variation in $\btheta_t$ is very
high during recessions such as the Great recession of $2008$, see the red--dashed line in panel (f). We
believe that practitioners can build on the variance decomposition
results to better understand the sources of variation in inflation. Indeed, as our preliminary results indicate, model uncertainty and estimation errors in the coefficients play a minor role compared to uncertainty due to time--variation in $\btheta_t$.
%\end{itemize}

\subsection[Forecasts]{Forecasts}

\begin{table}[!t]
\setlength\tabcolsep{3.0pt}{\footnotesize{}}%
\begin{tabular}{cl}
\hline
Model & Description\\
\hline\\
$\mathcal{M}_{0}$ & \parbox{14cm}{Plain AR(4) model: The constant term and $y_{t-1},\dots,y_{t-4}$ are always included. We set $\alpha=1$ and $\delta=1$.}\\[15pt]

$\mathcal{M}_{1}$ & \parbox{14cm}{Time-varying AR(4) model: The constant term and $y_{t-1},...,y_{t-4}$ are always included. We set $\alpha=0.99$ and average over $\delta_{1},\dots,\delta_{d}$.}\\[15pt]

$\mathcal{M}_{2}$ & \parbox{14cm}{DMA using $y_{t-1},\dots,y_{t-4}$: The constant term is always included. We set $\alpha=0.99$ and average over the combinations of $y_{t-1},...,y_{t-4}$ and $\delta_{1},\dots,\delta_{d}$.}\\[15pt]

$\mathcal{M}_{3}$ & \parbox{14cm}{DMA using $y_{t-1},\dots,y_{t-4}$ and the exogenous predictors: The constant term is always included. We set $\alpha=0.99$ and average over the combinations of predictors as well as $\delta_{1},\dots,\delta_{d}$.}\\[18pt]

$\mathcal{M}_{4}$ & \parbox{14cm}{DMS using $y_{t-1},\dots,y_{t-4}$ and the exogenous predictors: The constant term is always included. We set $\alpha=0.99$ and select the model with the highest posterior probability at each $t$ and use it to forecasts.}\\[18pt]

$\mathcal{M}_{5}$ & BMA: DMA with $\alpha=1$ and $\delta=1$.\\[10pt]
$\mathcal{M}_{6}$ & BMS: DMS with $\alpha=1$ and $\delta=1$.\\[10pt]

$\mathcal{M}_{7}$ & \parbox{14cm}{Kitchen Sink: The constant term, $y_{t-1},\dots,y_{t-4}$ and all exogenous predictors are always included.  We set $\alpha=0.99$ and average only over $\delta_{1},\dots,\delta_{d}$.}\\
\hline
\end{tabular}
\caption{Model specifications. The first column is the model index. The second column provides a brief description of each individual model.}
\label{tab:models}
\end{table}

A very important feature of DMA is out--of--sample forecasting,
see \cite{koop_korobilis.2011} and \cite{koop_korobilis.2012}. In this
section, we illustrate how our package can be used to perform forecasting.

In Table~\ref{tab:models}, we provide an overview of several alternative models.
Notice that, all models can be estimate using our package. For instance, the plain AR(4) model, ($\mathcal{M}_0$), can be estimated by setting $\delta=1.0$, $\alpha=1.0$, using the code:

\begin{CodeChunk}
\begin{CodeInput}
R> Fit_M0 <- DMA(GDPDEF ~  Lag(GDPDEF, 1) + Lag(GDPDEF, 2) +
                           Lag(GDPDEF, 3) + Lag(GDPDEF, 4),
                           data = USData, vDelta = 1.00,
                           dAlpha = 1.00, vKeep = c(1, 2, 3, 4, 5))
\end{CodeInput}
\end{CodeChunk}

The same holds for Bayesian Model Averaging (BMA, $\mathcal{M}_5$) and Bayesian Model Selection (BMS, $\mathcal{M}_6$) by setting
$\delta=1.0$, $\alpha=1.0$. Thus, \pkg{eDMA} also relates to the \pkg{BMS} package of \cite{zeugner_feldkircher.2015}.

We use the models to obtain one and five quarter ahead forecasts through
direct forecasting, see \cite{marcellino_etal.2006}.

Table~\ref{tab:back} reports the Mean Squared Error (MSE) and the Predictive Likelihood Difference (PLD) of
$\mathcal{M}_{i}$, $i=1,\dots,7$, over $\mathcal{M}_{0}$ (the benchmark)
at $h=1$ and $h=5$. Compared
to the benchmark, $\mathcal{M}_{1}$ provides relatively small improvements
at $h=1$, whereas at $h=5$, we obtain reductions in MSE by about
$17\%$ over the benchmark. The improvement in PLD at $h=5$ is also
notable. By averaging over $y_{t-1},...y_{t-4}$, we obtain more gains,
which indicates that allowing for time--variation in the number of
predictors is important. DMA using lags of inflation as well as $15$
additional predictors is the top performer, regardless of $h$. Conversely,
the improvements of DMS over $\mathcal{M}_0$ are much more modest. This
result is understandable as panel (c) in Figure~\ref{fig:combined} demonstrates that
there is no individual model that performs overwhelmingly better than
other specifications. As a result, DMS is outperformed by DMA.

As previously mentioned, DMA (DMS) with $\alpha=\delta=1$ correspond
to BMA (BMS). Overall, compared to
the benchmark model, BMA provides some improvements in density forecasts,
whereas we do not observe any improvements in terms of point forecasts.

\begin{table}[!t]
  \centering
    %\resizebox{1\columnwidth}{!}{%
  \setlength\tabcolsep{28pt}
  \begin{tabular}{lcccc}
  \toprule
Model  & \multicolumn{2}{c}{$h = 1$} & \multicolumn{2}{c}{$h = 5$}\\
  \cmidrule(lr){1-1}\cmidrule(lr){2-3}\cmidrule(lr){4-5}
 & MSE & PLD & MSE & PLD\\
\cmidrule(lr){1-1}\cmidrule(lr){2-3}\cmidrule(lr){4-5}
$\mathcal{M}_{1}$ & 1.020 & 0.244 & 0.835 & 20.522 \\
$\mathcal{M}_{2}$ & 0.982 & 1.544 & 0.694 & 37.192 \\
$\mathcal{M}_{3}$ & 0.969 & 11.700 & 0.645 & 73.819 \\
$\mathcal{M}_{4}$ & 1.038 & -4.890 & 0.780 & 26.996 \\
$\mathcal{M}_{5}$ & 0.976 & 6.809 & 1.127 & 9.234 \\
$\mathcal{M}_{6}$ & 1.187 & -18.358 & 1.310 & -33.352 \\
$\mathcal{M}_{7}$ & 1.762 & -7.125 & 1.242 & 32.293 \\
\bottomrule
  \end{tabular}
  %}
    \caption{Mean Squared Error (MSE) and Predictive Likelihood Difference (PLD) of $\mathcal{M}_{i},\quad i=1,\dots,7$ compared to $\mathcal{M}_{0}$ for $h = 1$ and $h = 5$ quarters ahead out--of--sample forecasts.}
        \label{tab:back}
\end{table}

Finally, as an alternative to these models, we can consider the Kitchen
Sink model (the model with all predictors, $\mathcal{M}_7$) where we only average over $\delta$. Compared to $\mathcal{M}_{0}$, the kitchen
sink model does not provide any improvements at $h=1$. At $h=5$,
we observe improvements in density forecasts. However, the kitchen
sink model is outperformed by DMA.

\subsection[Why does DMA performs so well ?]{Why does DMA performs so well ?}

To investigate how quickly our techniques adapt to changes in data,
we report the accumulated log--PLD for several models over the benchmark
in panels (a)--(d) of Figure~\ref{fig:pld_size}. These can be obtained using the \code{pred.like()} method available for \code{DMA} objects. For instance, we create the two vectors \code{vPL_M0} and \code{vPL_DMA} containing the log--Predictive Likelihood of $\mathcal{M}_0$ and $\mathcal{M}_3$ using:

\begin{figure}[!t]
\centering
\includegraphics[width=1\textwidth]{PriorAndESize.pdf}
\caption{Accumulated PDL and the optimal number of predictors for prior sensitivity analysis. Panel (a): $\mathcal{M}_1$ over $\mathcal{M}_0$, Panel (b): $\mathcal{M}_3$ over $\mathcal{M}_0$, Panel (c): $\mathcal{M}_5$ over $\mathcal{M}_0$, Panel (d): $\mathcal{M}_7$ over $\mathcal{M}_0$. Panels (e)--(h): Number of predictors for the model with the highest posterior probability using Zellenr's prior with $g = 0.5, 20, 100, T$. The gray vertical bars indicate business cycle peaks, \emph{i.e.}, the point at which an economic expansion transitions to a recession, based on National Bureau of Economic Research (NBER) business cycle dating.}
%%%
\label{fig:pld_size}
\end{figure}

\begin{CodeChunk}
\begin{CodeInput}
R> vPL_M0 <- pred.like(Fit_M0, iBurnPeriod = 32)
R> vPL_M3 <- pred.like(Fit, iBurnPeriod = 32)
\end{CodeInput}
\end{CodeChunk}

and compute the accumulated log-PLD of $\mathcal{M}_3$ over $\mathcal{M}_0$ as:
\begin{CodeChunk}
\begin{CodeInput}
R> vPLD_M3.M0  <- vPL_M3 - vPL_M0
\end{CodeInput}
\end{CodeChunk}

which is reported in panel (b) of Figure~\ref{fig:pld_size}.

In panels (a), (b), (c) and (d) of Figure~\ref{fig:pld_size} a value of zero corresponds to equal support of both
models, positive values are in support of the model of choice over
$\mathcal{M}_{0}$ and negative values show support of $\mathcal{M}_{0}$
over the model of choice at time $t$. In these panels, we decompose the effects
of (i): Allowing for time--variation in the regression coefficients,
(ii): Allowing for model uncertainty but no time--variation in the
regression coefficients and (iii): Allowing for time--variation in
the regression coefficients and model uncertainty.

In panel (a), we see that the time--varying AR(4) ($\mathcal{M}_1$) model outperforms
the benchmark during the recession in the early $1970$s. On the other
hand, both models provide basically same results. The Great Recession
is associated with a small improvements in favor of $\mathcal{M}_{1}$.
Compared to the plain AR(4) model, it takes about $30$ observations
to provide compelling evidence for DMA. These improvements increase
till the end of the sample. Furthermore, we also see that DMA improves
model performance in recession as well as expansion periods. Compared
to BMA, the improvements of DMA are mostly concentrated around recession
periods, which indicates the importance of allowing for time--variation of $\btheta_t$ in periods of economic turmoil.

\subsection[Prior sensitivity analysis]{Prior sensitivity analysis}
\label{sec:sensitivity}

In this section, we evaluate the robustness of our results with regards
to Zellner's prior shrinkage parameter, $g$, see (4) and (5) in \cite{dangl_halling.2012}. Intuitively, a smaller value
of $g$ means more shrinkage around the prior mean of $\btheta_0$, \emph{i.e.},
$\boldsymbol{0}$. Furthermore, if $g=0$, then $p\left(M_{i}\mid\mathcal{F}_{t}\right)$
is equal for all models. The larger is $g$, the more we are willing
to move away from the model priors in response to what we observe
in the data. Thus, the choice of $g$ can have important implications
with respect to out--of--sample forecasting results.

In order to shed light on this issue, we re--estimate DMA with $g$ equals to $0.1$,
$20$, $100$ and $T$ (using \code{bZellnerPrior = TRUE}) and see to which extent it influences out--of--sample
results, see Table~\ref{tab:priorsensitivity} and Figure~\ref{fig:pld_size}. Compared to our default priors, we see that very informative values such as $g=0.1$
and $20$ worsen results considerably as DMA tends to favor models with higher number of predictors, which obviously  does not improve forecast performance, see panels (e) and (f) of Figure~\ref{fig:pld_size}. As we increase $g$, we obtain closer
results to the default prior. Therefore, we generally recommend practitioners to use the default prior.

\begin{table}[!t]
  \centering
    %\resizebox{1\columnwidth}{!}{%
  \setlength\tabcolsep{55pt}
  \begin{tabular}{lcc}
  \toprule
Prior & MSE & PLD \\
\cmidrule(lr){1-1}\cmidrule(lr){2-3}
$g = 0.1$ & 4.875 & -175.145 \\
$g = 20$ & 1.255 & -27.141 \\
$g = 100$ & 1.072 & -0.573 \\
$g = T$ & 1.014 & 6.920 \\
\bottomrule
  \end{tabular}
  %}
    \caption{Mean Squared Error (MSE) and Predictive Likelihood Difference (PLD) of DMA using the following values of $g$: $0.1,20,100,T$ and $\mathcal{M}_{0}$ for $h = 1$.}
     \label{tab:priorsensitivity}
\end{table}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section[Conclusion]{Conclusion}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\label{sec:conclusion}
%
In this paper, we present the \pkg{eDMA} package for \proglang{R}. %This package allows to perform DMA as in \cite{dangl_halling.2012} accounting for different levels of parameter instability through time.
%
%DMA has been recently introduced in the economic and financial literature by \cite{koop_korobilis.2011} and \cite{koop_korobilis.2012} and is becoming widely used by practitioners and researchers in universities and central banks.
The purpose of \pkg{eDMA} is to offer an integrated environment to perform DMA using the available \code{DMA()} function, which enables practitioners to perform DMA exploiting multiple processors without major efforts. Furthermore, \proglang{R} users will find common methods to represent and extract estimated quantities such as \code{plot()}, \code{as.data.frame()}, \code{coef()} and \code{residuals()}.

Overall, \pkg{eDMA} is able to: (i): Incorporate the extensions introduced in \cite{dangl_halling.2012}, which is particularly relevant for economic and financial applications, (ii): Compared to other approaches, our package is much faster, (iii): It requires a smaller amount of RAM even in cases of moderately large applications, and (iv): It allows for parallel computing.

In Section \ref{sec:computational}, we also detail the expected time the program takes to perform DMA under different sample sizes, number of predictors and number of grid points. For typical economic applications, estimation time is around 30 min using a commercial laptop. Very large applications can still benefit from the use of \pkg{eDMA} when performed on desktop or even clusters, without additional effort from the user. %In \pkg{eDMA} the

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section*{Computational details}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

The results in this paper were obtained using \proglang{R} 3.2.3 \citep{R.2015} with the
packages: \pkg{eDMA} version 1.0 \citep{eDMA},
\pkg{Rcpp} version 0.12.5 \citep{Rcpp.2011,Rcpp}, \pkg{RcppArmadillo} version 0.7.100.3.1 \citep{RcppArmadillo.2014,RcppArmadillo}, \pkg{xts} version 0.9-7 \citep{xts} and \pkg{devtools} version 1.1.1 \citep{devtools}.
\proglang{R} itself and all packages used are available
from \proglang{CRAN} at \url{http://CRAN.R-project.org/}. The package \pkg{eDMA} is available from GitHub at \url{https://github.com/LeopoldoCatania/eDMA}.
Computations were performed on
a Genuine Intel\textregistered{} quad core CPU i7--3630QM 2.40Ghz processor.

\clearpage

%:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
\begin{thebibliography}{27}
\newcommand{\enquote}[1]{``#1''}
\providecommand{\natexlab}[1]{#1}
\providecommand{\url}[1]{\texttt{#1}}
\providecommand{\urlprefix}{URL }
\expandafter\ifx\csname urlstyle\endcsname\relax
  \providecommand{\doi}[1]{doi:\discretionary{}{}{}#1}\else
  \providecommand{\doi}{doi:\discretionary{}{}{}\begingroup
  \urlstyle{rm}\Url}\fi
\providecommand{\eprint}[2][]{\url{#2}}

\bibitem[{Byrne \emph{et~al.}(2014)Byrne, Korobilis, and
  Ribeiro}]{byrne_etal.2014}
Byrne JP, Korobilis D, Ribeiro PJ (2014).
\newblock \enquote{On the Sources of Uncertainty in Exchange Rate
  Predictability.}
\newblock \emph{Available at SSRN 2502586}.
\newblock
  \urlprefix\url{http://papers.ssrn.com/sol3/papers.cfm?abstract_id=2502586}.

\bibitem[{Catania and Nonejad(2016)}]{eDMA}
Catania L, Nonejad N (2016).
\newblock \emph{\pkg{eDMA}: Efficient Dynamic Model Averaging}.
\newblock \proglang{R} package version 1.0.0,
  \urlprefix\url{https://github.com/LeopoldoCatania/eDMA}.

\bibitem[{Chapman \emph{et~al.}(2008)Chapman, Jost, and Van
  Der~Pas}]{chapman_etal.2008}
Chapman B, Jost G, Van Der~Pas R (2008).
\newblock \emph{Using \code{OpenMP}: Portable Shared Memory Parallel
  Programming}, volume~10.
\newblock MIT press, Cambridge, US.

\bibitem[{Croissant and Millo(2008)}]{plm}
Croissant Y, Millo G (2008).
\newblock \enquote{Panel Data Econometrics in {R}: The \pkg{plm} Package.}
\newblock \emph{Journal of Statistical Software}, \textbf{27}(2).
\newblock \urlprefix\url{http://www.jstatsoft.org/v27/i02/}.

\bibitem[{Dangl and Halling(2012)}]{dangl_halling.2012}
Dangl T, Halling M (2012).
\newblock \enquote{Predictive Regressions with Time--Varying Coefficients.}
\newblock \emph{Journal of Financial Economics}, \textbf{106}(1), 157--181.
\newblock \doi{doi:10.1016/j.jfineco.2012.04.003}.

\bibitem[{Eddelbuettel and Fran\c{c}ois(2011)}]{Rcpp.2011}
Eddelbuettel D, Fran\c{c}ois R (2011).
\newblock \enquote{\pkg{Rcpp}: {S}eamless \proglang{R} and \proglang{C++}
  Integration.}
\newblock \emph{Journal of Statistical Software}, \textbf{40}(8), 1--18.
\newblock \doi{10.18637/jss.v040.i08}.

\bibitem[{Eddelbuettel \emph{et~al.}(2016{\natexlab{a}})Eddelbuettel,
  Fran\c{c}ois, Allaire, Ushey, Kou, Bates, and Chambers}]{Rcpp}
Eddelbuettel D, Fran\c{c}ois R, Allaire J, Ushey K, Kou Q, Bates D, Chambers J
  (2016{\natexlab{a}}).
\newblock \emph{\pkg{Rcpp}: {S}eamless \proglang{R} and \proglang{C++}
  Integration}.
\newblock \proglang{R} package version 0.12.5,
  \urlprefix\url{https://cran.r-project.org/package=Rcpp}.

\bibitem[{Eddelbuettel \emph{et~al.}(2016{\natexlab{b}})Eddelbuettel,
  Fran\c{c}ois, and Bates}]{RcppArmadillo}
Eddelbuettel D, Fran\c{c}ois R, Bates D (2016{\natexlab{b}}).
\newblock \emph{\pkg{RcppArmadillo}: \pkg{Rcpp} Integration for the
  \pkg{Armadillo} Templated Linear Algebra Library}.
\newblock \proglang{R} package version 0.7.100.3.1,
  \urlprefix\url{https://cran.r-project.org/package=RcppArmadillo}.

\bibitem[{Eddelbuettel and Sanderson(2014)}]{RcppArmadillo.2014}
Eddelbuettel D, Sanderson C (2014).
\newblock \enquote{\pkg{RcppArmadillo}: {A}ccelerating \proglang{R} with
  High--Performance \proglang{C++} Linear Algebra.}
\newblock \emph{Computational Statistics \& Data Analysis}, \textbf{71},
  1054--1063.
\newblock \doi{10.1016/j.csda.2013.02.005}.

\bibitem[{Groen \emph{et~al.}(2013)Groen, Paap, and
  Ravazzolo}]{groen_etal.2013}
Groen JJ, Paap R, Ravazzolo F (2013).
\newblock \enquote{Real--Time Inflation Forecasting in a Changing World.}
\newblock \emph{Journal of Business \& Economic Statistics}, \textbf{31}(1),
  29--44.
\newblock \doi{10.1080/07350015.2012.727718}.

\bibitem[{Koop and Korobilis(2011)}]{koop_korobilis.2011}
Koop G, Korobilis D (2011).
\newblock \enquote{UK Macroeconomic Forecasting with Many Predictors: Which
  Models Forecast Best and When Do They Do So?}
\newblock \emph{Economic Modelling}, \textbf{28}(5), 2307--2318.
\newblock \doi{10.1016/j.econmod.2011.04.008}.

\bibitem[{Koop and Korobilis(2012)}]{koop_korobilis.2012}
Koop G, Korobilis D (2012).
\newblock \enquote{Forecasting Inflation Using Dynamic Model Averaging.}
\newblock \emph{International Economic Review}, \textbf{53}(3), 867--886.
\newblock \doi{10.1111/j.1468-2354.2012.00704.x}.

\bibitem[{Marcellino \emph{et~al.}(2006)Marcellino, Stock, and
  Watson}]{marcellino_etal.2006}
Marcellino M, Stock JH, Watson MW (2006).
\newblock \enquote{A comparison of direct and iterated multistep AR methods for
  forecasting macroeconomic time series.}
\newblock \emph{Journal of Econometrics}, \textbf{135}(1-2), 499 -- 526.
\newblock ISSN 0304-4076.
\newblock \doi{10.1016/j.jeconom.2005.07.020}.
\newblock
  \urlprefix\url{http://www.sciencedirect.com/science/article/pii/S030440760500165X}.

\bibitem[{McCormick \emph{et~al.}(2016)McCormick, Raftery, and Madigan}]{dma}
McCormick TH, Raftery A, Madigan D (2016).
\newblock \emph{\pkg{dma}: Dynamic Model Averaging}.
\newblock R package version 1.2-3,
  \urlprefix\url{https://CRAN.R-project.org/package=dma}.

\bibitem[{Onorante and Raftery(2016)}]{onorante_raftery.2016}
Onorante L, Raftery AE (2016).
\newblock \enquote{Dynamic Model Averaging in Large Model Spaces Using Dynamic
  Occam's Window.}
\newblock \emph{European Economic Review}, \textbf{81}, 2--14.
\newblock \doi{10.1016/j.euroecorev.2015.07.013}.

\bibitem[{OpenMP(2008)}]{openmp.2008}
OpenMP A (2008).
\newblock \emph{\code{OpenMP} Application Program Interface, v. 3.0}.
\newblock \urlprefix\url{http://www.openmp.org/mp-documents/spec30.pdf}.

\bibitem[{Paye(2012)}]{Paye.2012}
Paye BS (2012).
\newblock \enquote{`D\'{e}j\`{a} vol': Predictive Regressions for Aggregate
  Stock Market Volatility Using Macroeconomic Variables.}
\newblock \emph{Journal of Financial Economics}, \textbf{106}(3), 527 -- 546.
\newblock ISSN 0304-405X.
\newblock \doi{10.1016/j.jfineco.2012.06.005}.
\newblock
  \urlprefix\url{http://www.sciencedirect.com/science/article/pii/S0304405X12001316}.

\bibitem[{\proglang{R} {Core Team}(2016)}]{R.2015}
\proglang{R} {Core Team} (2016).
\newblock \emph{\proglang{R}: {A} Language and Environment for Statistical
  Computing}.
\newblock \proglang{R} Foundation for Statistical Computing, Vienna, Austria.
\newblock \proglang{R} version 3.2.3,
  \urlprefix\url{https://www.R-project.org/}.

\bibitem[{Raftery \emph{et~al.}(2015)Raftery, Hoeting, Volinsky, Painter, and
  Yeung}]{BMA}
Raftery A, Hoeting J, Volinsky C, Painter I, Yeung KY (2015).
\newblock \emph{\pkg{BMA}: Bayesian Model Averaging}.
\newblock R package version 3.18.6,
  \urlprefix\url{https://CRAN.R-project.org/package=BMA}.

\bibitem[{Raftery \emph{et~al.}(2010)Raftery, K{\'a}rn{\`y}, and
  Ettler}]{raftery_etal.2010}
Raftery AE, K{\'a}rn{\`y} M, Ettler P (2010).
\newblock \enquote{Online Prediction Under Model Uncertainty via Dynamic Model
  Averaging: Application to a Cold Rolling Mill.}
\newblock \emph{Technometrics}, \textbf{52}(1), 52--66.
\newblock \doi{10.1198/TECH.2009.08104}.

\bibitem[{Ryan and Ulrich(2015)}]{xts}
Ryan JA, Ulrich JM (2015).
\newblock \emph{\pkg{xts}: {E}xtensible Time Series}.
\newblock \proglang{R}~package version 0.9-7,
  \urlprefix\url{https://CRAN.R-project.org/package=xts}.

\bibitem[{Sanderson(2010)}]{sanderson.2010}
Sanderson C (2010).
\newblock \enquote{\pkg{Armadillo}: {A}n Open Source \proglang{C++} Linear
  Algebra Library for Fast Prototyping and Computationally Intensive
  Experiments.}
\newblock \emph{Technical report}, NICTA.
\newblock \urlprefix\url{http://arma.sourceforge.net/}.

\bibitem[{Stock and Watson(1999)}]{stock_watson.1999}
Stock JH, Watson MW (1999).
\newblock \enquote{Forecasting Inflation.}
\newblock \emph{Journal of Monetary Economics}, \textbf{44}(2), 293--335.
\newblock \doi{10.1016/S0304-3932(99)00027-6}.

\bibitem[{Stock and Watson(2008)}]{stock_watson.2008}
Stock JH, Watson MW (2008).
\newblock \enquote{Phillips Curve Inflation Forecasts.}
\newblock \emph{Technical report}, National Bureau of Economic Research.
\newblock \urlprefix\url{http://www.nber.org/papers/w14322}.

\bibitem[{West and Harrison(1999)}]{west_harrison.1999}
West M, Harrison J (1999).
\newblock \emph{Bayesian Forecasting \& Dynamic Models}.
\newblock Springer--Verlag, Berlin.

\bibitem[{Wickham and Chang(2016)}]{devtools}
Wickham H, Chang W (2016).
\newblock \emph{\pkg{devtools}: Tools to Make Developing R Packages Easier}.
\newblock R package version 1.11.1,
  \urlprefix\url{https://CRAN.R-project.org/package=devtools}.

\bibitem[{Zeugner and Feldkircher(2015)}]{zeugner_feldkircher.2015}
Zeugner S, Feldkircher M (2015).
\newblock \enquote{Bayesian Model Averaging Employing Fixed and Flexible
  Priors: The BMS Package for R.}
\newblock \emph{Journal of Statistical Software}, \textbf{68}(1), 1--37.
\newblock \doi{10.18637/jss.v068.i04}.

\end{thebibliography}

\clearpage

\end{document}
