%\VignetteIndexEntry{Shrinkage in the Time-Varying Parameter Model Framework Using the R Package shrinkTVP}
%\VignetteEngine{R.rsp::tex}

\documentclass[article, nojss]{jss}

%% -- LaTeX packages and custom commands ---------------------------------------

%% recommended packages
\usepackage{caption}
\usepackage{thumbpdf,lmodern}
\usepackage{framed}
\usepackage{bm}
\usepackage{natbib}
\usepackage{color}
\usepackage[plain,noend]{algorithm2e}
\usepackage{graphicx}
\usepackage{amsmath,amssymb}
\usepackage[utf8]{inputenc}
\usepackage{multirow}
\usepackage{booktabs}
\usepackage{framed}
%\usepackage{lipsum}

\usepackage[usenames,dvipsnames]{xcolor}
%\newcommand{\comment}[1]{\textcolor{Green}{#1}}  % changes by Sylvia
\newcommand{\comment}[1]{#1}  % changes by Sylvia
\newcommand{\commentred}[1]{#1}  % changes by Sylvia
%% changes in the notation by Sylvia

\newcommand{\yhat}{\hat{y}}
\newcommand{\yS}{S}
\newcommand{\QQ}{\comment{\mathbf{Q}}} % Sylvia
%\renewcommand{\QQ}{Q} % AAP
\newcommand{\aphi}{\comment{a_\phi}} % Sylvia
%\newcommand{\aphi}{a_0^{sv}}  % AAP
\newcommand{\bphi}{\comment{b_\phi}} % Sylvia
%\newcommand{\aphi}{b_0^{sv}}  % AAP
\newcommand{\Bsv}{\comment{B_\sigma}} % Sylvia
%\newcommand{\Bsv}{B_\sigma^{sv} % AAP
\newcommand{\OOmega}{\comment{\Omegav}} % Sylvia
%\newcommand{\OOmega}{\Omega} % AAP
\newcommand{\SSigma}{\comment{\Sigmam}} % Sylvia
%\newcommand{\SSigma}{\Sigma} % AAP
 %% new commands by Sylvia
\newcommand{\bct}[2]{\comment{\beta_{{#1}{#2}}}} % j first, t second
\newcommand{\bvt}{\comment{\tilde{\betav}}}
\newcommand{\xv}{\comment{\bm x}}
\newcommand{\cv}{\comment{\bm c}}
\newcommand{\imm}[1]{^{\comment{(#1)}}}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\newcommand{\atauprold}{\lambda^{[a^{\tau}]}}  % ANGELA
\newcommand{\axiprold}{\lambda^{[a^{\xi}]}}    % ANGELA
\newcommand{\ataupr}{b^{\tau}}  % Sylvia
\newcommand{\axipr}{b^{\xi}}    % Sylvia
\newcommand{\gammav}{\boldsymbol{\gamma}}

\newcommand{\e}{\mbox{\rm e}}
\newcommand{\dimy}{r}                     % dimension of y_t
\newcommand{\tr}{{\tiny \mbox{\rm tr}}}
\newcommand{\ytr}{\ym^{\tr}}
\newcommand{\LPSo}[1]{\LPS^{\star}_{#1}}

\newcommand{\weg}[1]{}
\newcommand{\dx}{\,\mathrm{d}x}
\newcommand{\V}{\mathrm{V}}
\newcommand{\B}{\mathrm{B}}
\newcommand{\Exp}{\mathcal{E}}
\newcommand{\Cm}{{\mathbf C}}
\newcommand{\Bm}{{\mathbf B}}
\newcommand{\Am}{{\mathbf A}}
\newcommand{\Dm}{{\mathbf D}}
\newcommand{\cm}{{\mathbf c}}
\newcommand{\Fm}{{\mathbf F}}
\newcommand{\ystar}{y^\star}
\renewcommand{\arraystretch}{1.2}
\renewcommand{\vfill}{\vspace*{\fill}}
\newcommand{\nb}{\small}
\newcommand{\spc}{\phantom{$-$}}
% general parameters
\newcommand{\rvYm}{\mathbf{Y}}
\newcommand{\EVfs}{\mathrm{EV}}
\newcommand{\Gamfun}[1]{\Gamma (#1)}
\newcommand{\kfRc}{  R}  % symbol for the cov of the obs noise
\newcommand{\kfR}{{\mathbf \kfRc}}  % symbol for the cov of the obs noise
\newcommand{\pitrue}{\pi ^{\rm true}}
\newcommand{\thc}{\vartheta} %notation used to summarize  univariate unknown model parameters
\newcommand{\thmod}{{\mathbf{\boldsymbol{\thc}}}} %notation used to summarize  multivaraite   unknown vector
\newcommand{\thetav}{{\mathbf{\boldsymbol{\theta}}}} %notation used to summarize  multivaraite   unknown vector
%%%% Specific new commands for this paper
\newcommand{\iid}{\mbox{\rm i.i.d.}}
\newcommand{\bfz}{{\mathbf{0}}}
\newcommand{\ypro}{u}
\newcommand{\ypros}{y^{\star}}
\newcommand{\yprov}{\ym^{u}}
\newcommand{\yprodiff}{z}
\newcommand{\yprodiffv}{\mathbf{\yprodiff}}
\newcommand{\tm}{\mathbf{t}}
\newcommand{\errordiff}{\varepsilon}
\newcommand{\errordiffv}{\boldsymbol{\varepsilon}}
\newcommand{\psiv}{\boldsymbol{\psi}}
\newcommand{\tauv}{\boldsymbol{\tau}}
\newcommand{\xiv}{\boldsymbol{\xi}}
\newcommand{\Omegav}{\boldsymbol \Omega}
\newcommand{\sigmaerr}{R}
\newcommand{\alphad}{r} %sdimension of regression coefficient
\newcommand{\scale}{\omega} % scaling factor for variance heterogeneity
%\newcommand{\scalev}{{\mathbf{\scale}}}% fettdruck geht nicht
\newcommand{\scalev}{\boldsymbol{\scale}}
\newcommand{\im}[1]{^{(#1)}}
\newcommand{\labset}{L}
\newcommand{\Rv}{\mathbf{R}} % all mixture indicators
\newcommand{\Pm}{\mathbf{P}}
\newcommand{\Qm}{\mathbf{Q}}
\newcommand{\Sm}{\mathbf{S}}
\newcommand{\mv}{\mathbf{m}}
\newcommand{\sa}{Q}
\newcommand{\ev}{\mathbf{e}}
\newcommand{\z}{\phantom{0}}
\newcommand{\zz}{\phantom{00}}
\newcommand{\zzz}{\phantom{000}}
\newcommand{\Probsym}{\mbox{\rm Pr}}

\newcommand{\betac}{\beta} %symbol for regression parameter
\newcommand{\betad}{d} %sdimension of regression coefficient
\newcommand{\betar}{\boldsymbol{\betac}} % vector of regression coefficients
\newcommand{\betav}{\betar} % vector of regression coefficients - notation used in some chapter
\newcommand{\betai}[1]{\indiv{\betav}_{#1}} %notation for a subject specific parameter
\newcommand{\betavi}{\indiv{\betav}} %notation for a subject specific parameter
\newcommand{\betavitilde}{\indiv{\tilde{\betav}}}
\newcommand{\betainoncen}{b}
\newcommand{\zetai}[1]{\indiv{\betainoncen}_{#1}} %notation for a subject specific parameter
\newcommand{\zetaitilde}[1]{\indiv{\tilde{\betainoncen}}_{#1}} %notation for a subject specific parameter
\newcommand{\zetavi}[1]{\indiv{\mathbf{\betainoncen}}_{#1}} %notation for a subject specific parameter
\newcommand{\zetavitilde}[1]{\indiv{\tilde{\mathbf{\betainoncen}}}_{#1}} %notation for a subject specific parameter
\newcommand{\Psim}{\boldsymbol{\Psi}}
\newcommand{\Xbeta}{{  \mathbf \Xz}}
\newcommand{\Xz}{x}% symbol used for the row in a regressor matrix.
% \newcommand{\Normal}[1]{ \mathcal{N}\left(#1\right)}
\newcommand{\Normult}[2]{ \mathcal{N} _{#1}\left(#2\right)}
\newcommand{\Mulnom}[1]{\mbox{\rm MulNom}\left(#1\right)}
\newcommand{\bs}{b} %mean of prior and posterior of mu
\newcommand{\Bs}{B} %variance- prior and posterior of variance
\newcommand{\br}{{\mathbf{\bs}}} %mean of prior and posterior of beta
\newcommand{\Br}{{\mathbf{\Bs}}} %variance-covariance of prior and posterior of beta (independence prior)
\newcommand{\yc}{y}
\newcommand{\ydens}{\yc}  % argument of a univariate rv in a density
\newcommand{\ym}{{\mathbf \yc}} % y multivariate observation
\newcommand{\yv}{\ym}
\newcommand{\constant}{\mbox{\rm constant}}
\newcommand{\pdf}[3]{f_{ {\footnotesize #1}}(#2;#3)}
\newcommand{\Normalpdfa}[2]{\pdf{\mathcal{N}}{#1}{#2}}
\newcommand{\Normal}[1]{ \mathcal{N}\left(#1\right)}
\newcommand{\NoGo}[1]{ \mathcal{NG}\left(#1\right)}
\newcommand{\Normalpdf}{\varphi}
\newcommand{\Gammad}[1]{ \mathcal{G}\left(#1\right)}
\newcommand{\IGammad}[1]{ \mathcal{IG}\left(#1\right)}
\newcommand{\GIG}[3]{ \mathcal{GIG}\left(#1,#2,#3\right)}
\newcommand{\Lap}[3]{ Lap\left(#1,#2\right)}
\newcommand{\InvGau}[2]{ \mathcal{InvGau}\left(#1,#2\right)}
\newcommand{\xm}{\mathbf{x}} % vector of regression coefficients

\newcommand{\LPS}{{\mbox{\rm LPDS}}}

\newcommand{\Fd}[1]{\mbox{\rm F}\left(#1\right)}
\newcommand{\Td}[2]{\mbox{\rm T}_{#1}\left(#2\right)}
\newcommand{\indic}[1]{I\{#1\}}
\newcommand{\trans}[1]{#1 ^{'}} % transposed sign in linear algebra
% \newcommand{\Exp}[1]{\mathcal{E}\left(#1\right)}
\newcommand{\EV}{\mathcal{EV}}
\newcommand{\KS}{\mathcal{KS}}
\newcommand{\Kv}{{\mathbf  K}}
\newcommand{\Logistic}{\mathcal{LO}}
\newcommand{\GenLogistic}{\mathcal{LG}}
\newcommand{\Real}{\Re}
\newcommand{\Betadis}[1]{\mathcal{B}\left(#1\right)}
\newcommand{\Betafun}[1]{B (#1)}
\newcommand{\new}{^\mathrm{new}}
\newcommand{\old}{^\mathrm{old}}
\newcommand{\minusindex}[1]{ - #1 }
\newcommand{\Chi}{\chi^2}
\newcommand{\Chisqu}[1]{\chi^2 _{#1}}
\newcommand{\Xbetar}{\Xbeta^r} % design matrix for the random effects - vector
\newcommand{\Xbetaf}{\Xbeta^f} % design matrix for the fixed effects - vector
\newcommand{\alphaf}{\alpha} % notation for fixed effects
\newcommand{\alphav}{\boldsymbol{\alphaf}} % vector of regression coefficients
\newcommand{\Sigmam}{\boldsymbol{\Sigma}} % vector of regression coefficients
\newcommand{\omegav}{\boldsymbol{\omega}} % vector of regression coefficients
\newcommand{\Edis}[1]{\mbox{\rm E}\left(#1\right)} % Expectation of a rv in displaymode
\newcommand{\Bincoef}[2]{\left( \begin{array}{c}   #1 \\#2 \end{array}\right)}
\newcommand{\indivsymb}{s}
\newcommand{\indiv}[1]{ #1 ^{\indivsymb}} %notation for a subject specific parameter
\newcommand{\kfwc}{ w}
\newcommand{\kfw}{{\mathbf{\kfwc}}}  % symbol for the state noise
\newcommand{\wt}[1]{\kfw_{#1}}
\newcommand{\Qrcm}{{\mathbf{\kfQ}}} % covariance in random effects model
\newcommand{\mum}{\boldsymbol{\mu}} % multivariate mean of a normal distribtuion-
\newcommand{\Cholrcm}{{\mathbf{C}}} % group mean  in mixture model
\newcommand{\kfQc}{ Q }  % symbol for the cov of the state noise
\newcommand{\kfQ}{{\mathbf{\kfQc}}}  % symbol for the cov of the state noise
\newcommand{\Xbetamat}{{  \mathbf X}}
\newcommand{\Xbetarmat}{\Xbetamat^r} % design matrix for the random effects - matrix
\newcommand{\Xbetafmat}{\Xbetamat^f} % design matrix for the fixed effects - matrix
\newcommand{\error}{\varepsilon} % symbol for error
\newcommand{\errorv}{\boldsymbol{\error}}
\newcommand{\Ferror}{F_{\error}}
\newcommand{\ferror}{f_{\error}}
\newcommand{\Vrcm}{{\mathbf{V}}} % error variance in the marginal model
\newcommand{\Diag}[1]{\mbox{\rm Diag}\left(#1\right)}
\newcommand{\dimmat}[2]{(#1\times #2)} % Dimension of a matrix
\newcommand{\Bino}[1]{\mbox{\rm BiNom}\left(#1\right)}
\newcommand{\zs}{z}% symbol used for explanatory variable
\newcommand{\zv}{\mathbf{\zs}}% symbol used for explanatory variable
\newcommand{\betaci}[1]{\indiv{\betac}_{#1}} %notation for a subject specific parameter
\newcommand{\Gammainv}[1]{\mathcal{G}^{-1} \left(#1\right)}


\newcommand{\Uniform}[1]{\mathcal{U}\left[#1\right]}
\newcommand{\Ew}[1]{\mbox{\rm E}(#1)}   % Expectation of a rv
\newcommand{\Vw}[1]{\mbox{\rm V}(#1)}   % Expectation of a rv
\newcommand{\pl}{\pi}  %notation used for the probabilities in a logit model
\newcommand{\plv}{{\mathbf{\boldsymbol{\pl}}}}  %notation used for the probabilities in a logit model
\newcommand{\plt}[1]{\pl_{#1}} %probability with index for a binary model
\newcommand{\plsw}[2]{\pl_{#1,#2}} %swithcing probability, single component; first index:class, second
\newcommand{\identm}{{\mathbf I}}
\newcommand{\identy}[1]{{\identm}_{#1}}
\newcommand{\kfW}{{\mathbf{W}}}  % symbol for the weight matrix
\newcommand{\dm}{{\mathbf{d}}}  % symbol for the weight matrix
\newcommand{\minusmcmc}[1]{_{\minusindex{#1}}} % exclude MCMC
\newcommand{\Studmult}[2]{t _{#1}\left(#2\right)}
\newcommand{\rv}{{\mathbf{r}}}

\newcommand{\Lihood}[2]{L \left(#1 |#2 \right) }
\newcommand{\expND}[0]{-\frac{1}{2} \sum_{i=1}^{n}\left(\xm_{i} - \mu \right)' \rho \left(\xm_{i} - \mu   \right)}
\newcommand{\expNDtwo}[0]{         \sum_{i=1}^{n}\left(\xm_{i} - \mu \right)' \rho \left(\xm_{i} - \mu    \right)}
\newcommand{\expNDtwominusrho}[0]{ \sum_{i=1}^{n}\left(\xm_{i} - \mu \right)' \left(\xm_{i} - \mu         \right)}
\newcommand{\expNDthree}[0]{\sum_{i=1}^{n}\left(\xm_{i} - \bar \xm \right)' \rho \left(\xm_{i} - \bar \xm \right)}
\newcommand{\expNDfour}[0]{\sum_{i=1}^{n}\left(\xm_{i} - \bar \xm \right)'  \left(\xm_{i} - \bar \xm      \right)}

%% new custom commands
\newcommand{\class}[1]{`\code{#1}'}
\newcommand{\fct}[1]{\code{#1()}}
\newtheorem{alg}{Algorithm}

%% -- Article metainformation (author, title, ...) -----------------------------

%% - \author{} with primary affiliation
%% - \Plainauthor{} without affiliations
%% - Separate authors by \And or \AND (in \author) or by comma (in \Plainauthor).
%% - \AND starts a new line, \And does not.
\Plainauthor{Angela Bitto-Nemling, Annalisa Cadonna, Sylvia Fr\"uhwirth-Schnatter, Peter Knaus}

\author{Angela Bitto-Nemling
\And Annalisa Cadonna\\WU Vienna
\And Sylvia Fr\"uhwirth-Schnatter\\WU Vienna
\AND Peter Knaus\\WU Vienna}

%% - \title{} in title case
%% - \Plaintitle{} without LaTeX markup (if any)
%% - \Shorttitle{} with LaTeX markup (if any), used as running title
\title{Shrinkage in the Time-Varying Parameter Model Framework Using the \proglang{R} Package \pkg{shrinkTVP}}
\Plaintitle{Shrinkage in the Time-Varying Parameter Model Framework Using the R Package \pkg{shrinkTVP}}
\Shorttitle{Shrinkage for TVP Models Using \pkg{shrinkTVP}}


%% - \Abstract{} almost as usual
%\Abstract{
%The  \proglang{R} package \pkg{shrinkTVP}
%provides a fully Bayesian implementation of time-varying parameters (TVP) models with shrinkage priors both on the variances and the initial values means, with the aim to automatically reduce time-varying parameters to static ones, and shrink the static parameters to zero. The package employs a Markov chain Monte
%Carlo (MCMC) sampler to obtain draws from the posterior
%distribution of the parameters and provides summary and plot methods to extract the main information from the posterior draws. The package can straightforwardly utilized by a beginner R user, while allowing fine tuning to more experienced partitioner.
%The goal of this paper is present the functionalities of \pkg{shrinkTVP}. in addition, it provides an brief introduction to TVP models and variance shrinkage, and a few working examples.
%}
\Abstract{
Time-varying parameter (TVP) models are widely used in time series analysis to
flexibly deal with processes which gradually change over time. However, the risk of
overfitting in TVP models is well known. This issue can be dealt with using appropriate
global-local shrinkage priors, which pull time-varying parameters towards static ones. In this paper, we introduce the  \proglang{R} package \pkg{shrinkTVP} \citep{kna-etal:shr},  which
provides a fully Bayesian implementation of shrinkage priors for TVP models, taking advantage of recent developments in the literature, in particular that of \cite{bit-fru:ach}.
The package \pkg{shrinkTVP}  allows for posterior simulation of the parameters through an efficient Markov Chain Monte Carlo
scheme. Moreover, summary and visualization methods, as well as the possibility of assessing predictive performance through log predictive density scores, are provided.
The computationally intensive tasks have been implemented in \proglang{C++} and interfaced with \proglang{R}.  The paper includes a brief overview
of the models and shrinkage priors implemented in the package. Furthermore, core functionalities are illustrated, both with simulated and real data.
}




%% - \Keywords{} with LaTeX markup, at least one required
%% - \Plainkeywords{} without LaTeX markup (if necessary)
%% - Should be comma-separated and in sentence case.
%% - \Keywords{} with LaTeX markup, at least one required
%% - \Plainkeywords{} without LaTeX markup (if necessary)
%% - Should be comma-separated and in sentence case.
\Keywords{Bayesian inference, Gibbs sampler, Markov chain Monte Carlo (MCMC),
 normal-gamma prior,  time-varying parameter (TVP) models, log predictive density scores}

%% - \Address{} of at least one author
%% - May contain multiple affiliations for each author
%%   (in extra lines, separated by \{and}\\).
%% - May contain multiple authors for the same affiliation
%%   (in the same first line, separated by comma).

%% - \Address{} of at least one author
%% - May contain multiple affiliations for each author
%%   (in extra lines, separated by \emph{and}\\).
%% - May contain multiple authors for the same affiliation
%%   (in the same first line, separated by comma).
\Address{
Peter Knaus\\
Institute for Statistics and Mathematics\\
Department of Finance, Accounting and Statistics\\
WU Vienna University of Economics and Business\\
Welthandelsplatz 1, Building D4, Entrance A, 3rd Floor\\
1020 Vienna, Austria\\
E-mail: \email{peter.knaus@wu.ac.at}
URL: \url{https://www.wu.ac.at/statmath/faculty-staff/faculty/peter-knaus}
}





\begin{document}
%% -- Introduction -------------------------------------------------------------

%% - In principle "as usual".
%% - But should typically have some discussion of both _software_ and _methods_.
%% - Use \proglang{}, \pkg{}, and \code{} markup throughout the manuscript.
%% - If such markup is in (sub)section titles, a plain text version has to be
%%   added as well.
%% - All software mentioned should be properly \cite-d.
%% - All abbreviations should be introduced.
%% - Unless the expansions of abbreviations are proper names (like "Journal
%%   of Statistical Software" above) they should be in sentence case (like
%%   "generalized linear models" below).

\section[Introduction]{Introduction} \label{sec:intro}

 %\textcolor{green}{(To discuss: abstract submitted by me (Angela) in december or different version -if yes, why? McCausland details and LPDS. Algorithm details; To Do for ABN: check bib references; are we consistent with notation especially hyperparameters etc in the documentation and formulas;  )
 %}.

Time-varying parameter (TVP) models are
widely used in time series analysis, because of their flexibility and ability to capture gradual changes in the model \comment{parameters} %variables
over time. The popularity of TVP models in macroeconomics and finance is based on the fact that, in most applications, the influence of certain predictors on the outcome variables varies over time \citep{pri:tim, dan-hal:pre, bel-etal:hie_tv}. TVP models, while capable of reproducing salient features of the data in a very effective way, present a concrete risk of overfitting, as often only a small subset of the parameters are time-varying. Hence, in the last decade, there has been a growing need for models and methods able to discriminate between time-varying and static parameters in TVP models.
A key contribution in this direction was the introduction of the non-centered parameterization of TVP models in  \cite{fru-wag:sto}, which recasts the problem of variance selection and shrinkage in terms of variable selection, thus allowing any tool used to this end in multiple regression models to be used to perform selection or shrinkage of variances.
\cite{fru-wag:sto} employ a spike and slab prior, while continuous shrinkage priors have been utilised as a regularization alternative in, e.g., \cite{bel-etal:hie_tv} and
\cite{bit-fru:ach}. For an excellent review of shrinkage priors, with a particular focus on high dimensional regression, the reader is directed to  \cite{bha-etal:las}.

In this paper, we describe the \proglang{R}  package  \pkg{shrinkTVP}  \citep{kna-etal:shr} for Bayesian TVP models with shrinkage. The package is available under the general public license (GPL $\geq$ 2) from the Comprehensive R
Archive Network (CRAN) at \url{https://cran.r-project.org/web/packages/shrinkTVP}.
The package efficiently implements recent developments in the Bayesian literature, in particular the ones
presented in \cite{bit-fru:ach}. The  computationally intensive Markov chain Monte Carlo (MCMC) algorithms
in the package are written in \proglang{C++} and interfaced with \proglang{R} \citep{R}  via the \pkg{Rcpp} \citep{edd-bal:ext} and the \pkg{RcppArmadillo} \citep{edd-san:acc} packages.
This approach combines the ease-of-use of R and its underlying functional programming paradigm with the computational speed of \proglang{C++}.

\commentred{The  package \pkg{shrinkTVP} is designed to provide an easy entry point for fitting TVP models with shrinkage priors, while also giving more experienced users the option to adapt the model to their needs. This is achieved by providing a robust baseline model that can be estimated by only passing the data, while also allowing the user to specify more advanced options.
Additionally, the \pkg{shrinkTVP} package is  designed to
ensure compatibility with well-known times series formats and to complement other packages. As input objects,  time series from the  \proglang{R} packages  \pkg{zoo}  \citep{zei:zoo} and \pkg{xts} \citep{xts} as well as time series formats like \code{ts} are supported.
Estimation output   is compatible  with the popular  \proglang{R} package  \pkg{coda}  \citep{plu:cod}   which can  be easily applied for  convergence diagnostic tests, among others.
Coupled with intuitive summary and plot methods, \pkg{shrinkTVP} is  a package that is  easy to use while remaining highly flexible. }

\commentred{\pkg{shrinkTVP} is, to our knowledge, the only \proglang{R} package that combines TVP models with shrinkage priors on the time-varying components  in a Bayesian framework.
Several  \proglang{R} packages deal  with statistical inference for  various specific classes of  state space models, of which TVP models are a special case.
The most popular \proglang{R} package in this field is \pkg{dlm} \citep{pet:anr},  a comprehensive package providing routines for maximum likelihood estimation, Kalman filtering and smoothing, and Bayesian analysis for dynamic linear models  (DLMs).
The  accompanying  book \citep{pet-etal:dyn} introduces the methodology and many \proglang{R} code examples.
%The prior proposed in (INCLUDE PAGE OF BOOK) does not serve as shrinkage prior as discussed in \citep{Sylvia and Helga}.
As of now, priors are not designed to encourage shrinkage and   \pkg{shrinkTVP}  complements   \pkg{dlm} in this regard.}
%Needless to say, it is possible to manually extend parts of the \pkg{dlm} package with appropriate routines to receive a TVP shrinkage model. Certainly, we recommend the user to facilitate her workload by usinginstead.
% which also provides parameter settings to reproduce findings in the \pkg{dlm} package %ramework.
%WHAT I WANNA SAY IS: you can either extend \pkg{dlm} to shrinkage or reduce shrinkTVP to %\pkg{dlm} only, as the editors ask for complementing other packages.

\commentred{The \proglang{R} package \pkg{bvarsv} \citep{kru:bva} implements  Bayesian inference for  vector autoregressive (VAR) models with time-varying parameters  (TVP-VAR)
and stochastic volatility for multivariate time series  as introduced by  \citep{pri:tim}. %, including noted by \citep{del-pro:tim2}.
 We refer to \citep{del-pro:tim2} for details on  the MCMC algorithm  and a later correction of the original scheme.
 %Thein \pkg{bvarsv} requires  at each iteration multiple calls of a recursive Kalman filter type algorithm,
 In addition to the very user friendly estimation function \code{bvar.sv.tvp},  \pkg{bvarsv} provides posterior predictive distributions and  enables impulse response analysis.  The package  includes the macroeconomic data set analysed  in \citep{pri:tim}  as   example data set,  \code{usmacro.update},  which we use in our predictive exercise in Section~\ref{sec:usmacro} to showcast the effect of introducing shrinkage priors on time-varying parameters.}


%\comment{ No normal-gamma shrinkage packages exist which can be easily linked or attached to any time-varying parameter model package.
%The reason for this absence is clearly due to the structure of a TVP plus shrinkage model. Its estimation demands a sophisticated MCMC scheme as developed in \citep{bit-fru:ach}.
%But \pkg{shrinkTVP} is, to our knowledge, the only \proglang{R} package that combines time-varying parameter models with shrinkage priors. }

%In the TVP models context, the most popular \proglang{R} package is \pkg{dlm} %\citep{pet:anr}, which provides routines for maximum likelihood \comment{estimation}, %Kalman filtering and smoothing, and Bayesian analysis for dynamic linear models %(DLMs), of which TVP models are a subset.
%\comment{The package \pkg{bsts}  \citep{sco:bst} performs  Bayesian analysis for the %closely related class of  structural time series models.} Moreover, a number of %\proglang{R} packages providing regularization and shrinkage methods are available.  %For example, \pkg{shrink} \citep{dun-etal:shr} implements various shrinkage methods %for linear, generalized linear, or Cox regressions, \pkg{biglasso} \citep{zen:big} %aims at very fast lasso-type models for high-dimensional linear regression and %\pkg{glmnet} \citep{sim-etal:reg} provides efficient procedures for implementing %elastic-net regularization for a variety of models. With regards to Bayesian %shrinkage, the normal-beta prime shrinkage prior is implemented in the package %\pkg{NormalBetaPrime} \citep{bai:nor} and the popular horseshoe prior in the package %\pkg{horseshoe} \citep{pas:hor}. Both of these packages focus on high-dimensional %regression models, and do not provide shrinkage for TVP models.

\commentred{Additional packages emerged very recently. %  in this highly active field.
The \proglang{R} package \pkg{tvReg} \citep{tvreg}  presents a user friendly compendium of  many common linear  TVP models,  including  standard linear regression as well as autoregressive,  seemingly unrelated equation  and  VAR models.
Estimation is based on kernel smoothing techniques. % One specific feature  of this  package are its customizable plotting and summary functions.
For an illustrative  application,  a TVP-VAR(4) model  is fitted  % and time-varying impulse response functions  are
  to the \code{usmacro.update} data set mentioned above,   using the   function  \code{tvVAR}.}
  \commentred{The  \proglang{R} package \pkg{walker} \citep{walker}  facilitates the estimation of DLMs and generalized  DLMs using MCMC  algorithms provided by \proglang{Stan} \citep{carpenter2017stan}.  For inference, the importance sampling method of \citep{ISMCMC} is implemented within a Hamiltonian Monte Carlo framework.}
\commentred{The \proglang{R} package \pkg{bsts} \citep{sco:bst} performs Bayesian analysis for  structural time series models, a highly  relevant  class of  state space models   including DLMs.  \pkg{bsts} is a very powerful package that  allows
shrinkage  for static  regression coefficients using spike and slab priors.
However, as for any other  packages mentioned above, variation of the dynamic  components is not regularized in   \pkg{bsts}.}

 \commentred{A main contribution of  the  package \pkg{shrinkTVP}  is  bridging the  active field of  \proglang{R} packages for state space models with the even more active field of \proglang{R} packages that  provide regularization and shrinkage methods  for common regression type models.
 % such as linear regression or generalized linear regression,  typically in a high-dimensional setting.
Among others,  \pkg{ncvreg} \citep{ncvreg} is useful for fitting standard penalized regression estimators,
 \pkg{glmnet} \citep{sim-etal:reg} allows  elastic-net regularization for a variety of models,
   \pkg{horseshoe} \citep{pas:hor}  implements the  horseshoe prior, while  \pkg{shrink} \citep{dun-etal:shr} provides various shrinkage methods for linear, generalized linear, and Cox regression models. \pkg{biglasso} \citep{zen:big} aims at very fast lasso-type regularization  for high-dimensional linear regression.
Recent \proglang{R} packages include   \pkg{NormalBetaPrime} \citep{bai:nor}
 for  Bayesian univariate
and \pkg{MBSP} \citep{MBSP} for  Bayesian multivariate linear regression analysis using, respectively,  the normal-beta prime  and the three parameter beta normal family for inducing shrinkage. The \proglang{R} package  \pkg{monomvn} \citep{monomvn}  employs a  normal-gamma prior in the specific situation
of  Bayesian inference for  multivariate normal and Student-$t$ data with a monotone  pattern of missing data.
}
%Again time varying parameter models are not covered in this package.
% All of these packages focus on high-dimensional regression models, and do not provide shrinkage for TVP models.
% }



%{\color{magenta}


%The \proglang{R} package \pkg{bayesreg} \citep{bayesreg} (SHALL WE INCLUDE THIS ONE? ) is again an %\proglang{R} package offering Bayesian linear regression with Gaussian or heavy-tailed error models and %Bayesian logistic regression with ridge, lasso, horseshoe and horseshoe+ estimators. Time varying parameter %models are not covered in this package.
% }




 \commentred{The remainder of the paper is organized as follows.
Section~\ref{sec:model} briefly introduces TVP \comment{models} and normal-gamma shrinkage priors, and describes the MCMC algorithm for posterior simulation. The package \pkg{shrinkTVP} is introduced in Section~\ref{sec:pkgshrinkTVP}. In particular, we illustrate how to run the MCMC sampler using the main function \code{shrinkTVP}, how to choose a specific model, and how to conduct posterior inference using the return object of \code{shrinkTVP}. Section~\ref{sec:LPDS} explains how to assess model   performance   by calculating log predictive density scores (LPDSs), and how to use LPDSs to compare the predictive performances of different priors. This is illustrated using the \code{usmacro.update} data set. % from the \pkg{bvarsv} \citep{kru:bva} package.
Finally, Section~\ref{sec:conclusions} concludes the paper.}


\section{Model specification and estimation } \label{sec:model}

\subsection{TVP models}
%In particular, we focus on the non-centered parametrization of state space models, which  recasts the variance selection problem as a variable selection problem \citep{fru-wag:sto}.
Let us recall the state space form of a TVP model. For $t = 1, \ldots, T$, we have that
\begin{equation}
\label{eq:centeredpar}
\begin{aligned}
&y_{t} =   \bm x_t \bm {\beta_{t}}  +  \epsilon_{t} , \qquad
\epsilon_{t} \sim \mathcal N (0,\sigma^2_t), \\
& \bm {\beta}_{t} =  \bm {\beta}_{t-1} + \bm w_{t}, \qquad   \bm w_{t}  \sim \mathcal N_d (0, \QQ),
\end{aligned}
\end{equation}
where $y_t$ is a univariate response variable and $\bm x_t = (x_{t 1}, x_{t 2}, \ldots, x_{t  d})$ is a $d$-dimensional row vector containing the regressors at time $t$, with $x_{t 1}$ corresponding to the intercept.
For simplicity, we assume here that $\QQ = \text{Diag}(\theta_1, \ldots, \theta_d)$ \comment{is a diagonal matrix}, implying that the state innovations are conditionally independent.
Moreover, we assume the initial value follows a normal distribution, i.e., $\bm \beta_{0} \sim \mathcal N_d (\bm \beta, \QQ)$\comment{, with initial mean
$\bm \beta = (\beta_1, \ldots, \beta_d)$}.
Model (\ref{eq:centeredpar}) can be rewritten equivalently in the non-centered parametrization as
\begin{equation}
\label{eq:noncenteredpar}
\begin{aligned}
&y_t= \bm x_t \bm \beta +   \bm x_t \text{Diag}(\sqrt{\theta}_1, \ldots, \sqrt{\theta}_d)
\tilde{\bm \beta}_{t} +  \epsilon_t, \quad  \epsilon_t \sim \mathcal N (0,\sigma^2_t),\\
&\tilde{\bm \beta}_{t} =\tilde {\bm \beta}_{t-1} + \tilde{\bm u}_{t}, \qquad \tilde{\bm u}_{t} \sim  \mathcal N_d (0, I_d),
\end{aligned}
\end{equation}
with $ \tilde{\bm \beta}_{0} \sim \mathcal{N}_d (\bm 0, I_d)$, where $I_d$ is the $d$-dimensional identity matrix.

\pkg{shrinkTVP} is capable of modelling the observation error both homoscedastically, i.e., $\sigma^2_t \equiv \sigma^2$ for all $t = 1, \ldots, T$ and heteroscedastically, via a stochastic volatility \comment{(SV)}  specification. In the latter case, the log-volatility $h_t = \log \sigma^2_t$ follows an AR(1) model \citep{jac-etal:bayJBES,kas-fru:anc, kas:dea}. More specifically,
\begin{eqnarray}
	\label{eq:svht}
	h_t | h_{t-1}, \mu, \phi, \sigma_\eta^2 \sim \mathcal{N} \left(\mu +\phi ( h_{t-1}-\mu),\sigma^2_\eta \right),
\end{eqnarray}
with initial state $h_0 \sim \mathcal N \left(\mu, \sigma_\eta^2/(1-\phi^2) \right)$.
The stochastic volatility model on the errors can prevent the detection of spurious variations in the TVP coefficients \citep{nak:tim, sim:com} by capturing some of the variability in the error term.


\subsection{Prior Specification} \label{sec:priors}

\subsubsection{Shrinkage priors on variances and model parameters}

We place conditionally independent normal-gamma priors \citep{gri-bro:inf} both on the standard deviations of the innovations, that is the $\sqrt{\theta_j}$'s,
and on the means of the initial value $\beta_j$, for $j = 1, \ldots, d$. Note that, in the case of the standard deviations,
this can equivalently be seen as a double gamma prior on the innovation variances $\theta_j$, for $j = 1, \ldots, d$.
This prior can be represented as a conditionally normal distribution, with the component specific variance following a gamma distribution, that is
\begin{eqnarray} \label{eq:normal}
\sqrt{\theta}_j|\xi^{2}_j  \sim \Normal{0,\xi_j^{2}}, \qquad \xi_j^2|a^\xi,\kappa^2  \sim  \Gammad{a^\xi,\frac{a^\xi \kappa^2}{2}}, \\
\beta_{j}|\tau^2_j \sim  \Normal{0,\tau^2_j}, \qquad \tau_j^2|a^\tau ,\lambda^2 \sim  \Gammad{a^\tau,\frac{a^\tau \lambda^2}{2}}.
\end{eqnarray}
The   \comment{prior variances} % component specific parameters
$\xi^{2}_j$ and $\tau^2_j$ are referred to as {local shrinkage parameters}, in that they control the strength with which each individual parameter \comment{$\sqrt{\theta}_j$ and $\beta_j$} is pulled \comment{toward} zero. The parameters $\kappa^2$ and $\lambda^2$ are dubbed {global shrinkage parameters}, as they determine how strongly {all} parameters are pulled to zero. \comment{Since $\Ew{\theta_j|a^\xi,\kappa^2}=2/\kappa^2$ and $\Ew{\beta^2_{j}|a^\tau ,\lambda^2}=2/\lambda^2$, the}  %. The
 larger  $\kappa^2$ and $\lambda^2$, the stronger  \comment{this effect}. % the effteco.
Finally, we refer to $a^\xi$ and $a^\tau$ as {shrinkage adaption parameters}. As $a^\xi$ and $a^\tau$ decrease, marginally more mass is placed around zero and jointly more mass is put on sparse specifications of the model. In particular, setting the local adaption parameters, $a^\xi$ and $a^\tau$, equal to one results in a Bayesian Lasso \citep{par-cas:bay} prior on the \comment{$\sqrt\theta_j$'s} and the \comment{$\beta_j$'s}, respectively.

The parameters $\kappa^2$, $\lambda^2$, $a^\xi$, $a^\tau$ can be learned from the data through appropriate prior distributions.
As priors for the global shrinkage parameters, we use
\begin{align} \label{eq:equNG03}
	\kappa^2 \sim \mathcal G (d_1, d_2), \qquad \lambda^2 \sim \mathcal G (e_1, e_2).
\end{align}
Moreover, in order to learn the shrinkage adaption parameters, we generalize the approach taken in \cite{bit-fru:ach} and place the following gamma distributions as priors:
\begin{align} \label{eq:equNG02}
a^\xi\sim \mathcal G(\nu^\xi, \nu^\xi b^\xi), \qquad
a^\tau \sim \mathcal G(\nu^\tau, \nu^\tau b^\tau),
\end{align}
which correspond to the exponential priors used in \cite{bit-fru:ach} when $\nu^\xi=1$ and $\nu^\tau=1$. The parameters $\nu^\xi$ and $\nu^\tau$ act as degrees of freedom and allow the prior to be bounded away from zero.




\subsubsection{Prior on the volatility parameter}


In the homoscedastic case we employ a hierarchical prior, \comment{where the scale of an inverse gamma prior  for $\sigma^2$ follows a gamma distribution}, that is,
\begin{eqnarray} \label{eq:priorsigma}
	\sigma^2|C_0 \sim \Gammainv{c_0,C_0}, \qquad  C_0 \sim \Gammad{g_0,G_0},
\end{eqnarray}
with hyperparameters $c_0$, $g_0$, and $G_0$.

In the case of stochastic volatility, the priors on the parameters $\mu$, $\phi$ and $\sigma^2_\eta$ in Equation~\eqref{eq:svht} are chosen as in \citet{kas-fru:anc}, that is
\begin{eqnarray} \label{eq:volpriors}
	\mu  \sim \mathcal{N}( b_\mu, B_\mu ), \quad \dfrac{\phi +1 }{2} \sim \mathcal{B}(\aphi, \bphi), \quad  \sigma^2_\eta \sim \mathcal{G}(1/2, 1/2 \Bsv ),
\end{eqnarray}
\comment{with hyperparameters $b_\mu, B_\mu, \aphi, \bphi,$ and $\Bsv$.}

\subsection{MCMC sampling algorithm}

\label{sec:MCMC}
The package \pkg{shrinkTVP} implements an MCMC Gibbs sampling algorithm with Metropolis-Hastings steps to obtain draws from the posterior distribution of the model parameters.
Here, we roughly sketch the sampling algorithm and refer the interested reader to \cite{bit-fru:ach} for further details.

\newpage

\commentred{\begin{alg} \label{facsvalg}
\mbox{\rm Gibbs Sampling Algorithm}
 \begin{enumerate} \itemsep 0mm
\item[\mbox{\rm 1.}] \mbox{\rm Sample the latent states} $\tilde \betav =( \tilde \betav_0, \ldots, \tilde \betav_T)$ \mbox{\rm  in the non-centered parametrization from  a}\\
   \mbox{\rm  multivariate  normal  distribution;}
%
\item[\mbox{\rm 2.}] \mbox{\rm Sample jointly} $\beta_1, \dots,\beta_d,$ \mbox{\rm  and} $\sqrt{\theta_1},\dots,\sqrt{\theta_d}$ \mbox{\rm  in the non-centered parametrization from}\\
    \mbox{\rm  a multivariate  normal  distribution;}
%
\item[\mbox{\rm 3.}] \mbox{\rm Perform an ancillarity-sufficiency interweaving step and redraw  each} $\beta_1, \dots,\beta_d$ \mbox{\rm from a}\\
\mbox{\rm   normal  distribution and each} ${\theta_1},\dots,{\theta_d}$ \mbox{\rm from a generalized inverse Gaussian distribution}\\
\mbox{\rm  using \pkg{GIGrvg} \citep{hoe-ley:gig};}
%
\item[\mbox{\rm 4.}]
\begin{enumerate}
	\item[\mbox{\rm (a)}] \mbox{\rm Sample the variance shrinkage adaption parameter} $ a^\xi$ \mbox{\rm  using a random walk}\\
\mbox{\rm Metro\-polis-Hastings step;}
	\item[\mbox{\rm (b)}] \mbox{\rm Sample the parameter shrinkage adaption parameter} $a^\tau$ \mbox{\rm   using a random walk}\\
\mbox{\rm   Metro\-polis-Hastings step;}
\end{enumerate}
%
\item[\mbox{\rm 5.}]
\begin{enumerate}
\item[\mbox{\rm (a)}] \mbox{\rm Sample the local shrinkage parameters} $\xi_j^2$, \mbox{\rm for} $j = 1, \ldots, d$,  \mbox{\rm from conditionally}\\ \mbox{\rm  independent  generalized inverse Gaussian distributions;}
%
\item[\mbox{\rm (b)}]  \mbox{\rm Sample the local shrinkage parameters} $\tau_j^2$, \mbox{\rm for} $j = 1, \ldots, d$,  \mbox{\rm from conditionally}\\ \mbox{\rm  independent  generalized inverse Gaussian distributions;}
\end{enumerate}
%
\item[\mbox{\rm 6.}]  \mbox{\rm Sample the error variance} $\sigma^2$ \mbox{\rm from an inverse gamma distribution in the homoscedastic}
 \mbox{\rm case or, in the SV case, sample the level} $\mu$, \mbox{\rm the persistence} $\phi$, % and
    \mbox{\rm  the volatility of the vola-}\\ \mbox{\rm tility} $\sigma^2_{\eta}$  \mbox{\rm and the log-volatilities} $\bm h= (h_0, \ldots, h_T)$
     \mbox{\rm  using \pkg{stochvol} \citep{kas:dea}.}
\end{enumerate}
\end{alg}}

When fitting the  %hierarchical
model \comment{under the full hierarchical shrinkage prior defined in Equations~\eqref{eq:normal}--\eqref{eq:equNG02}}, all steps in Algorithm~\ref{facsvalg} are performed, while steps 4(a), 4(b), 5(a) and 5(b) are skipped in certain prior setups.

One key feature of the algorithm is the joint sampling of the time-varying parameters $\tilde{\bm \beta}_t$, for $t=0, \ldots, T$ in step 1 of Algorithm~\ref{facsvalg}. We employ the procedure described in
\cite{mcc-etal:sim} which exploits the sparse, block tri-diagonal structure of the precision matrix of the full conditional distribution of
 \comment{$\tilde \betav =( \tilde \betav_0, \ldots, \tilde \betav_T)$},  %$\tilde{\bm \beta}_t$, for $t=0, \ldots, T$,
 to speed up computations.

Moreover, as described in \cite{bit-fru:ach}, in step 3 we make use of the  ancillarity-sufficiency interweaving strategy (ASIS) introduced by  \citet{yu-men:cen}. ASIS is well known to improve mixing by sampling certain parameters both in the centered and non-centered parameterization.
This strategy has been successfully applied to univariate SV models \citep{kas-fru:anc}, multivariate factor SV models \citep{kas-etal:eff}  and  dynamic linear state space models \citep{sim-etal:int}.



\section[The shrinkTVP package]{The \pkg{shrinkTVP} package}
\label{sec:pkgshrinkTVP}




\subsection{Running the model}
The core function of the package \pkg{shrinkTVP} is the function \code{shrinkTVP}, which serves as an R-wrapper for the actual sampler coded in \proglang{C++}. The function works out-of-the-box, meaning that estimation can be performed with minimal user input. With default settings, the TVP model in \comment{Equation}~\eqref{eq:centeredpar} is estimated in a Bayesian fashion with priors~\eqref{eq:normal} to~\eqref{eq:equNG02} and the following choice for the hyperparameters: $d_1 = d_2 = e_1 = e_2 = 0.001$, $\nu^\xi=\nu^\tau=5$ and $b^\xi=b^\tau=10$, \comment{implying a prior mean of  $\Ew{a^\xi}= \Ew{a^\tau}=0.1$}.
The error is assumed to be homoscedastic, with prior defined in \comment{Equation}~\eqref{eq:priorsigma} and hyperparameters $c_0 = 2.5$, $g_0 = 5$, and $G_0 = g_0/(c_0 - 1)$.

The only compulsory argument is an object of class \ ``formula'', which most users will be familiar with (see, for example, the use in the function \code{lm} in the package \pkg{stats} \citep{R}). The second argument is an optional data frame, containing the response variable and the covariates. Exemplary usage of this function is given in the code snippet below, along with the default output.
All code was on run on a personal computer with an Intel i5-8350U CPU.

\begin{CodeChunk}
\begin{CodeInput}
R> library("shrinkTVP")
R>
R> # Baseline model
R> set.seed(123)
R> sim <- simTVP(theta = c(0.2, 0, 0), beta_mean = c(1.5, -0.3, 0))
R> data <- sim$data
R> res <- shrinkTVP(y ~ x1 + x2, data = data)
\end{CodeInput}
\begin{CodeOutput}
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
Timing (elapsed): 4.647 seconds.
3228 iterations per second.

Converting results to coda objects and summarizing draws... Done!
\end{CodeOutput}
\end{CodeChunk}

Note that the data in the example is generated by the function \code{simTVP}, which can create synthetic datasets of varying sizes for illustrative purposes. The inputs \code{theta} and \code{beta} can be used to specify the true $\theta_1, \ldots, \theta_d$ and $\beta_1, \ldots, \beta_d$ used in the data generating process, in order to evaluate how well \code{shrinkTVP} recaptures these true values. The values correspond to the ones used in the synthetic example of \cite{bit-fru:ach}.

The user can specify the following MCMC algorithm parameters:  \code{niter}, which determines the number of MCMC iterations including the burn-in,
\code{nburn}, which equals the number of MCMC iterations discarded as burn-in, and
\code{nthin}, indicating the thinning parameter, meaning that every nthin-th draw is kept and returned.
The default values are \code{niter = 10000}, \code{nburn = round(niter/2)} and \code{nthin = 1}.

\commentred{The user is strongly encouraged to check convergence of the  produced Markov chain, especially for a large number of covariates. The output is made \pkg{coda} compatible, so that the user can utilize the tools provided by the excellent \proglang{R} package  to assess convergence.}





\subsection{Specifying the priors}  \label{sec3:priors}


More granular control over the prior setup can be exercised by specifying all or a subset of the parameters. In addition to changing the hyperparameters given in Section~\ref{sec:priors}, the user can fix  \comment{one or both} %the
values of the global shrinkage \comment{parameters} ($\kappa^2$, $\lambda^2$) and the shrinkage adaption parameters ($a^\tau$, $a^\xi$), instead of learning them from the data as done in the default setting. The benefit of this is twofold: on the one hand, desired degrees of sparsity and global shrinkage can be achieved through fixing the hyperparameters; on the other hand, interesting special cases arise from setting certain values of hyperparameters. For example, setting the local adaption parameters, $a^\xi$ and $a^\tau$, equal to one results in a Bayesian Lasso \citep{par-cas:bay} prior on the \comment{$\sqrt\theta_j$'s} and the \comment{$\beta_j$'s}, respectively. If the user desires a higher degree of sparsity, this can be achieved by setting the shrinkage adaption parameters to a value \comment{closer} to zero.  Table~\ref{tab:tablepriors} gives an overview of different model specifications. Note that separate prior choices can be made for the variances and the means of the initial values.


\begin{table}[]
	\centering
	\begin{tabular}{@{}lllll@{}}
	\toprule
									& \multicolumn{2}{l}{Shrinkage on $\sqrt{\theta_j}$}                          & \multicolumn{2}{l}{Shrinkage on $\beta_j$}                           \\ \midrule
									& $a^\xi$                               & $\kappa^2$              & $a^\tau$                                 & $\lambda^2$             \\ \cmidrule (r{4pt}){2-3} \cmidrule (l){4-5}
	Full hierarchical  \comment{shrinkage prior} % model
& $\mathcal{G}(\nu^\xi, \nu^\xi b^\xi)$ & $\mathcal{G}(d_1, d_2)$ & $\mathcal{G}(\nu^\tau, \nu^\tau b^\tau)$ & $\mathcal{G}(e_1, e_2)$ \\
	Hierarchical normal-gamma prior & fixed                                 & $\mathcal{G}(d_1, d_2)$ & fixed                                    & $\mathcal{G}(e_1, e_2)$ \\
	Normal-gamma prior              & fixed                                 & fixed                   & fixed                                    & fixed                   \\
\comment{Bayesian Lasso} & \comment{fixed at 1} &  fixed   & \comment{fixed at 1} & fixed \\
\bottomrule
	\end{tabular}
	\caption{Overview of different possible model specifications.}
	\label{tab:tablepriors}
\end{table}

In the following, we give some examples of models that can be estimated with the \pkg{shrinkTVP} package. In particular, we demonstrate how certain combinations of input arguments correspond to different model specifications. Note that in the following snippets of code, the argument \code{display_progress} is always set to \code{FALSE}, in order to suppress the progress bar and other outputs.

\paragraph{Fixing the shrinkage adaption parameters}

It is possible \comment{to set} %fix
the shrinkage adaption parameter $a^\xi$($a^\tau$) \comment{to a fixed value} through the input argument \code{a_xi} (\code{a_tau}), after setting \code{learn_a_xi} (\code{learn_a_tau}) to \code{FALSE}. As an example, we show how to fit a hierarchical Bayesian Lasso, both on the $\sqrt{\theta_j}$'s and on the $\beta_j$'s:

\begin{CodeInput}
R> res_hierlasso <- shrinkTVP(y ~ x1 + x2, data = data,
+    learn_a_xi = FALSE,  learn_a_tau = FALSE,
+    a_xi = 1, a_tau = 1, display_progress = FALSE)
\end{CodeInput}


\paragraph{Fixing the global shrinkage parameters}
The user can choose to fix the value of $\kappa^2$($\lambda^2$) by specifying the argument \code{kappa2} (\code{lambda2}), after setting \code{learn_k2} (\code{learn_lambda2}) to \code{FALSE}. In the code below, we give an example on how to fit a (non-hierarchical) Bayesian Lasso on both $\sqrt{\theta_j}$'s and $\beta_j$'s, with corresponding global shrinkage parameters fixed both to $100$:

\begin{CodeInput}
R> res_lasso <- shrinkTVP(y ~ x1 + x2, data = data,
+    learn_a_xi = FALSE, learn_a_tau = FALSE, a_xi = 1, a_tau = 1,
+    learn_kappa2 = FALSE, learn_lambda2 = FALSE, kappa2 = 100, lambda2 = 100,
+    display_progress = FALSE)
\end{CodeInput}



\subsection{Stochastic volatility specification}

The stochastic volatility specification defined in Equation~\eqref{eq:svht} can be used by setting the option \code{sv} to \code{TRUE}. This is made possible by a call to the \code{update_sv} function exposed by the \pkg{stochvol} package.
The code below fits a model in which all the parameters are learned and the observation equation errors are modeled through stochastic volatility:

\begin{CodeInput}
R> res_sv <- shrinkTVP(y ~ x1 + x2, data = data, sv = TRUE,
+    display_progress = FALSE)
\end{CodeInput}

The priors on the SV parameters are the ones defined in Equation~\eqref{eq:volpriors}, with hyperparameters fixed to
$b_\mu = 0$ , $B_\mu = 1$, $\aphi = 5$, $\bphi = 1.5$ , and $\Bsv = 1$.

\subsection{Specifying the hyperparameters}

Beyond simply switching off parts of the hierarchical structure of the prior setup, users can also modify the hyperparameters governing the prior distributions. This can be done through the arguments \code{hyperprior_param} and \code{sv_param}, which both have to be named lists.

In addition to user defined hyperparameters, unspecified parameters will be set to default values, as defined in Section~\ref{sec3:priors}.



\begin{CodeInput}
R> res_hyp <- shrinkTVP(y ~ x1 + x2, data = data,
+    hyperprior_param = list(b_xi = 5, nu_xi = 10),
+    display_progress = FALSE)
\end{CodeInput}



\subsection{Posterior inference: Summarize and visualize the posterior distribution}

The return value of \code{shrinkTVP} is an object of type \code{shrinkTVP}, which is a named list containing a variable number of elements, depending on the prior specification. For the default model, the values are:

\begin{enumerate}
	\itemsep0em
	\item the parameter draws of $\sigma^2$  in \code{sigma2},
	\item the parameter draws of $(\sqrt{\theta_1}, \dots, \sqrt{\theta_d})$ in \code{theta_sr},
	\item the parameter draws of $\bm \beta = (\beta_1, \dots, \beta_d)$ in \code{beta_mean},
	\item a list holding $d$ \code{mcmc.tvp} objects (one for each $\comment{\bm{\beta}_j=(\bct{j}{0}, \ldots, \bct{j}{T})}$)
containing the parameter draws in \code{beta},
	\item the parameter draws of $\xi_1^2,  \ldots,\xi_d^2,$ in \code{xi2},
	\item the parameter draws of $a^{\xi}$ in \code{a_xi},
	\item some acceptance statistics for the Metropolis-Hastings step for $a^{\xi}$ in \code{a_xi_acceptance},
	\item the parameter draws of $\tau_1^2, \ldots, \tau_d^2$ in \code{tau2},
	\item the parameter draws of $a^{\tau}$ in \code{a_tau},
	\item some acceptance statistics for the Metropolis-Hastings step for $a^{\tau}$ in \code{a_tau_acceptance},
	\item the parameter draws for $\kappa^2$ in \code{kappa2},
	\item the parameter draws for $\lambda^2$ in \code{lambda2},
	\item the parameter draws of $C_0$ in \code{C0},
	\item the prior hyperparameters in \code{priorvals},
	\item the design matrix, the response and the formula in \code{model}, %\footcomment{ad 15:Data matrix?}
	\item summary statistics for the parameter draws in \code{summaries}
	 and
	%\textcolor{blue}{
objects required for the \code{LPDS} function in  \code{LPDS_comp}.
\end{enumerate}
When some parameters are fixed by the user, the corresponding output value is omitted. %\footcomment{Do the fixed parameters become part of  the priorvals?}
In the SV case, the draws for the parameters of the SV model on the errors are contained in \code{sv_mu}, \code{sv_phi} and \code{sv_sigma}.  For details, see \cite{kas:dea}.

The two main tools for summarizing the output of \code{shrinkTVP} are the \code{summary} and \code{plot} methods implemented for \code{shrinkTVP} objects. \code{summary} has two arguments beyond the \code{shrinkTVP} object itself, namely \code{digits} and \code{showprior}, which control the output displayed. \code{digits} indicates the number of decimal places to round the posterior summary statistics to, while \code{showprior} determines whether or not to show the prior distributions resulting from the user input. In the example below, the default \code{digits} value of 3 is used, while the prior specification is omitted. The output of \code{summary} consists of the mean, standard deviation, median, 95\% \comment{highest} posterior density region and effective sample size (ESS) for the non time-varying parameters.


\begin{CodeChunk}
\begin{CodeInput}
R> summary(res, showprior = FALSE)
\end{CodeInput}
\begin{CodeOutput}
Summary of 5000 MCMC draws after burn-in of 5000.

Statistics of posterior draws of parameters (thinning = 1):

 param                   mean    sd        median HPD 2.5% HPD 97.5% ESS
 sigma2                  0.994   0.125     0.985  0.75     1.238     1340.972

 abs(theta_sr_Intercept) 0.418   0.063     0.411  0.298    0.538     365.77
 abs(theta_sr_x1)        0.014   0.024     0.001  0        0.065     126.673
 abs(theta_sr_x2)        0.003   0.007     0      0        0.015     759.047

 beta_mean_Intercept     0.128   0.399     0      -0.355   1.211     327.235
 beta_mean_x1            -0.242  0.159     -0.268 -0.485   0.008     111.541
 beta_mean_x2            -0.001  0.034     0      -0.085   0.07      4311.249

 xi2_Intercept           97.956  4131.527  0.243  0.005    9.821     5000
 xi2_x1                  0.578   28.674    0      0        0.061     5000
 xi2_x2                  0.019   0.53      0      0        0.004     3737.874

 a_xi                    0.084   0.038     0.077  0.021    0.16      231.749

 tau2_Intercept          31.334  1645.674  0      0        3.025     4062.968
 tau2_x1                 727.492 50802.992 0.087  0        4.855     5000
 tau2_x2                 0.112   1.704     0      0        0.069     3790.064

 a_tau                   0.084   0.039     0.078  0.019    0.16      127.561

 kappa2                  31.075  89.542    1.929  0        168.95    3950.702

 lambda2                 61.895  212.079   2.761  0        314.075   961.691

 C0                      1.726   0.632     1.648  0.635    2.988     4971.06
\end{CodeOutput}
\end{CodeChunk}

%\textcolor{blue}{
The \code{plot} method can be used to visualize the posterior distribution estimated by \code{shrinkTVP}. Aside from a \code{shrinkTVP} object, its main argument is \code{pars}, a character vector containing the names of the parameters to visualize. \code{plot} will call either \code{plot.mcmc.tvp} from the \pkg{shrinkTVP} package if the parameter is time-varying or \code{plot.mcmc} from the \pkg{coda} package, if the parameter is non time-varying.  The default value of \code{pars} is \code{c("beta")}, leading to \code{plot.mcmc.tvp} being called on each of the $\beta_{jt}$, for $j =1, \ldots, d$. See the code below for an example and Figure~\ref{fig:beta} for the corresponding output.
%}
\begin{CodeInput}
R> plot(res)
\end{CodeInput}

\begin{figure}[h!]
	\centering
	\includegraphics[width=\linewidth]{./figs/plot_beta.pdf}
	\caption{Visualization of the  evolution of the \comment{time-varying  parameter  $\bm{\beta}_j=(\bct{j}{0}, \ldots, \bct{j}{T}), j=1, 2, 3,$
	over time  $t=0,\ldots,T$}, as provided by the \code{plot} method. \code{plot} is in turn calling \code{plot.mcmc.tvp} on the individual \code{mcmc.tvp} objects. The median is displayed as a black line, and the shaded areas indicate the \comment{pointwise} 95\% and 50\% posterior credible intervals.}
\label{fig:beta}
\end{figure}

%\textcolor{blue}{
The \code{plot.mcmc.tvp} method displays empirical posterior credible intervals of a time-varying parameter over time, i.e., $\beta_{jt}$, for $j =1, \ldots, d$ and $\sigma^2_t$ in the case of stochastic volatility. By default, the \comment{pointwise} 95\% and 50\% posterior credible intervals are displayed as shaded areas layered on top of one another, with the median represented by a black line, with an additional grey, dashed line at zero. To ensure that users have flexiblity in the plots created, a host of options are implemented for customisation. The bounds of the credible intervals can be modified through the \code{probs} input, allowing for different levels of uncertainty visualization. The arguments \code{quantlines} and \code{shaded} take boolean vectors as inputs, and determine if the corresponding credible intervals will be displayed through shading and/or lines. The shaded areas can be customised via the arguments \code{shadecol} and \code{shadealpha}, which determine the color and the degree of transparency of the shaded areas. The lines representing the quantiles can be adjusted through \code{quantlty}, \code{quantcol} and \code{quantlwd}, which modify the line type, color and line width, respectively. In the spirit of \proglang{R}, all of these arguments are vectorised and the supplied vectors are recycled in the typical \proglang{R} fashion if necessary. The first element of these vectors is always applied to the outermost credible interval, the second to the second outermost and so forth. The horizontal line at zero can be similarly adjusted through \code{zerolty}, \code{zerolwd} and \code{zerocol} or entirely turned off by setting \code{drawzero} equal to \code{FALSE}. All further arguments are passed on to the standard \code{plot} method, allowing for changes to the line representing the median and other plot modifications that users of \proglang{R} are familiar with. An example of possible customisation can be seen in the code below, with the corresponding output being Figure~\ref{fig:beta2}.
%}

\begin{CodeInput}
R> plot(res, pars = "beta", xlim = c(100, 200),
+    probs = seq(0.1, 0.9, by = 0.1),
+    quantlines = T, quantcol = color[c(1:4, 4:1)], quantlty = 1,
+    quantlwd = 3, col = color[5], lwd = 3, shadecol = "gold1")
\end{CodeInput}
%In the code below, this behavior is exploited to add a title to the plot, resulting in Figure~\ref{fig:sigma2}.

%Aside from the \code{shrinkTVP} object, its arguments are the %\textcolor{red}{mandatory} \code{pars}, a character vector containing the names of the %parameters to visualize, and in the case of time-varying parameters \code{nplot}, %which controls the number of plots to display per page. The names supplied in %\code{pars} have to coincide with the names of the list elements of the %\code{shrinkTVP} object. The default value is \code{c("beta")}, i.e., empirical %quantiles of the posterior of $\bm \beta_t$ over time are shown. If there are too many %plots to fit on one page, \code{plot} will cycle through all parameters to display. It %will call either \code{plot.mcmc} from the \pkg{coda} package, if the parameter is non %time-varying, or \code{plot.mcmc.tvp} from the \pkg{shrinkTVP} package for %time-varying parameters, passing all additional arguments specified via \code{...} to %the respective plotting functions.
%See the code below for an example and  Figure~\ref{fig:beta} for the corresponding output.


%\footcomment{According to our definition of $\xm_t$, {\tt x1} corresponds to $x_{t2}$ and {\tt x2} corresponds to $x_{t3}$ in Figure~1-3. x-label in Figure~1 should be $t$.}
%\textcolor{blue}{
%In the case of large dimension $d$, the argument \code{nplot} controls the number of plots to display per page before a new page is generated. It should be noted that this only affects the plots of time varying parameters, as \code{nplot} does not change how \pkg{coda} handles multiple pages of plots.}
%\textcolor{blue}{
%The \code{plot} method will pass any additional arguments on to the respective plotting functions, allowing for some flexibility in the displayed plots. As mentioned,
%\code{plot.mcmc.tvp} defaults to plots using colored shadings. This can be toggled via the boolean vectors \code{quantlines} and  \code{shaded}, which determine if the corresponding quantiles will be displayed through shading and/or lines.
%Moreover, the user can customise the shaded areas and the pointwise quantiles lines though a variety of arguments. Two examples of customized plots, for $\bm \beta_{t}$ and $\sigma^2_t$, are demonstrated in the code below and the resulting Figure~\ref{fig:beta2} and Figure~\ref{fig:sigma2}, respectively.
%}


\begin{figure}[h!]
\centering
\includegraphics[width=\linewidth]{./figs/plot_beta_2.pdf}
\caption{Visualization of the  evolution of  the \comment{time-varying  parameter  $\bm \beta_{t}$ over time  $t=100,\ldots,200$ for $j=1, \ldots, 3$}. In this example, the x-axis of the plot was restricted with \code{xlim}, the color of the shaded areas was changed to yellow and colored solid lines have been added to delimit the credible intervals. The colored lines represent the median and the \comment{pointwise}	10\%, 20\%, 30\% 40\%, 60\%, 70\%, 80\%, and 90\% quantiles.}
\label{fig:beta2}
\end{figure}


To visualize other parameters via the \code{plot} method, the user has to change the \code{pars} argument. \code{pars} can either be set to a single character object or to a vector of characters containing the names of the parameter draws to display. In the latter case, the \code{plot} method will display groups of plots at a time, prompting the user to move on \comment{to} the next series of plots, similarly to how \pkg{coda} handles long plot outputs. Naturally, as all parameter draws are converted to \pkg{coda} objects, any
\comment{method  from this package  that users are familiar with} %methods that users are familiar with from this package
(e.g., to check convergence) can be applied to the parameter draws contained in a \code{shrinkTVP} object. An example of this can be seen in Figure~\ref{fig:theta}, where \code{pars = "theta_sr"},  \comment{changes} % changing
the output to a graphical summary of the parameter draws of $\sqrt{\theta_1}, \dots, \sqrt{\theta_d}$, using \pkg{coda}'s \code{plot.mcmc} function. To obtain Figure~\ref{fig:theta}, one can run
%
\begin{CodeInput}
R> plot(res, pars = "theta_sr")
\end{CodeInput}



\begin{figure}[h!]
\centering
\includegraphics[width=\textwidth,height=\textheight,keepaspectratio]{./figs/plot_theta.pdf}
\caption{\comment{Trace plots (left column) and kernel density estimates of the posterior density (right column) for the parameters
$\sqrt{\theta}_1, \dots, \sqrt{\theta}_3$}, as provided by the \code{plot} method. \code{plot} is in turn calling \pkg{coda}'s \code{plot.mcmc}.}
\label{fig:theta}
\end{figure}



\section{Predictive performances and model comparison}
\label{sec:LPDS}

Within a Bayesian framework, a natural way to predict a future observation is through its posterior predictive density.  For this reason, log-predictive density scores (LPDSs) provide a means of assessing how well the model performs in terms of prediction on real data. The log-predictive density score for time $t_0 +1$
is obtained by evaluating at $y_{t_0 +1}$ the
\comment{log of the posterior predictive density} % obtained by fitting the model to the previous  data points $y_{1}, \ldots, y_{t_0}$.}
%log-predictive density
obtained by fitting the model to the previous $t_0$ data points.
Given the data up to time $t_0$, the posterior predictive density at time $t_0 + 1$ is given by
\begin{align} \label{MCpred}
p(y_{t_0 + 1}| y_{1}, \ldots, y_{t_0}, \bm x_{t_0 +1} ) =  \int p(y_{t_0 + 1}| \bm x_{t_0 +1}, \bm \psi )  p (\bm \psi| y_{1}, \ldots, y_{t_0} )d \bm \psi,
\end{align}
where $\bm \psi$ is the set of model parameters \comment{and latent variables up to $t_0+1$.  For a TVP model with homoscedastic  errors, $\bm \psi=  (\tilde {\bm \beta}_0,\ldots \tilde {\bm \beta}_{t_0 +1}, \sqrt \theta_1, \ldots, \sqrt{\theta_d}, \beta_1, \ldots, \beta_d, \sigma^2)$,
whereas for a TVP model with  SV errors,
 $\bm \psi=  (\tilde {\bm \beta}_0,\ldots \tilde {\bm \beta}_{t_0 +1}, \sqrt \theta_1, \ldots, \sqrt{\theta_d}, \beta_1, \ldots, \beta_d,
 \sigma^2_1, \ldots,\sigma^2_{t_0 +1})$.}
\comment{Given $M$ samples from the posterior distribution of the parameters and latent variables, $p (\bm \psi| y_{1}, \ldots, y_{t_0} )$, Monte Carlo integration
could be applied immediately to approximate (\ref{MCpred}).
However, \cite{bit-fru:ach} propose a more efficient approximation of the predictive density,   \comment{the so-called conditionally optimal Kalman mixture approximation which is}
obtained by analytically integrating out $\tilde {\bm \beta}_{t_0+1}$ from the likelihood at time $t_0 +1$.}

% In our model   $\bm \psi=  (\tilde {\bm \beta}_0,\ldots \tilde {\bm \beta}_T, \sqrt \theta_1, \ldots, \sqrt{\theta_d}, \beta_1, \ldots, \beta_d, \sigma^2)$ and $\bm \psi=  (\tilde {\bm \beta}_0,\ldots \tilde {\bm \beta}_T, \sqrt \theta_1, \ldots, \sqrt{\theta_d}, \beta_1, \ldots, \beta_d, \sigma^2_1, \ldots,\sigma^2_T)$ in the SV case. However, in a TVP model, the likelihood depends on the one-step ahead value of the states, $\tilde {\bm \beta}_{t_0+1}$, which are neither observed nor estimated.  As a solution to this, \cite{bit-fru:ach} propose an approximation of the predictive density,



 In the homoscedastic error case, given $M$ samples from the posterior distribution of the parameters  \comment{and the latent variables up to $t_0$},
 \commentred{Monte Carlo integration  of  the resulting predictive density yields following   mixture approximation,
\begin{eqnarray} \label{eq:approx_mixture}
\hspace{-0.15cm}
 p(y_{t_0+1}| y_{1},\ldots, y_{t_0}, \xv_{t_0+1}) &\approx& \dfrac{1}{M} \sum_{m=1}^M
\Normalpdfa{y_{t_0+1}}{\yhat _{t_0+1} \imm{m}, \yS _{t_0+1} \imm{m} }, \\
 \yhat \imm{m} _{t_0+1}& =& \bm x_{t_0+1} \bm \beta^{(m)} + \bm F_{t_0+1}\imm{m}  \bm m_{t_0} \imm{m}, \nonumber  \\
 \yS \imm{m} _{t_0+1} &= & \bm F_{t_0+1} \imm{m} (\SSigma_{t_0} \imm{m} + I_d) ( \bm F_{t_0+1}\imm{m}) ^\top + (\sigma^{2})\imm{m} , \nonumber
%%&\bm F_{t_0+1}\imm{m}  \bm m_{t_0} \imm{m}
%% &p(y_{t_0+1}| \bm x_{t_0+1} \bm \beta^{(m)} + \\
%%&\bm F_{t_0+1}\imm{m}  \bm m_{t_0} \imm{m}, \bm F_{t_0+1} \imm{m} (\SSigma_{t_0} \imm{m} + I_d)
% %( \bm F_{t_0+1}\imm{m}) ^\top + (\sigma^{2})\imm{m} ),
\end{eqnarray}
 where  the conditional predictive  densities  are  Gaussian and   the  conditional  moments   depend on the MCMC draws.
 The   mean  $\yhat _{t_0+1} \imm{m}$  and  the variance  $\yS _ {t_0+1} \imm{m} $}
are computed for  the  $m$th MCMC iteration from
 $\bm F_{t_0+1}=  \bm x_{t_0+1} \text{Diag}(\sqrt{\theta_1}, \ldots, \sqrt{\theta_d)}$ and
 the mean  $\bm m_{t_0}$ and  the covariance matrix  $\SSigma_{t_0}$  of the posterior distribution of $\tilde {\bm \beta}_{t_0}$.
These quantities can be obtained by iteratively calculating $\SSigma_{t}$ and $\bm m_{t}$ up to time $t_0$, as described in \cite{mcc-etal:sim}:
\begin{align*}
&\SSigma_1 = (\OOmega_{11})^{-1}, \qquad \bm m_1 = \SSigma_1 \bm c_1,\\
&\SSigma_t = (\OOmega_{tt} - \OOmega_{t-1,t}^{\top} \SSigma_{t-1} \OOmega_{t-1,t})^{-1}, \qquad \bm m_t = \SSigma_t (\bm c_t - \OOmega_{t-1,t}^{\top} \bm m_{t-1} ).
\end{align*}
The quantities $\cv_t$,  $\OOmega_{tt}$ and  $\OOmega_{t-1,t}$ for $t=1,  \ldots, t_0$ are given in Appendix~\ref{sec:mccau}.

For the SV case,  it is still  possible to analytically integrate  out $\tilde {\bm \beta}_{t_0+1}$ from the likelihood at time $t_0 +1$ conditional on
a known value of $\sigma^2_{t_0+1}$, however it is not possible to integrate  the likelihood with respect to both  latent variables $\tilde {\bm \beta}_{t_0+1}$  and $\sigma^2_{t_0+1}$.
Hence, at each MCMC iteration a draw is taken from the predictive distribution of  $\sigma^2_{t_0+1}=\exp (h_{t_0+1})$, derived from Equation~\eqref{eq:svht},  and used to calculate the \comment{conditional predictive density of   $y_{t_0+1}$}.
The approximation of the one-step ahead predictive density can then be obtained through the following steps:
\begin{enumerate}
	\item    for each MCMC draw of $(\mu, \phi,\sigma_{\eta}^2) \imm{m}$ and $h_{t_0} \imm{m}$,  obtain a  draw of $(\sigma^{2}_{t_0+1})\imm{m}$; 	\item calculate the   conditionally optimal Kalman   mixture approximation as \commentred{in (\ref{eq:approx_mixture}) with following
slightly different values $\yS _ {t_0+1} \imm{m} $:}
\commentred{	\begin{align*}
	%\hspace{-1cm}
&  \yS _ {t_0+1} \imm{m} =  \bm F_{t_0+1}  \imm{m} (\SSigma_{t_0} \imm{m} + I_d) ( \bm F_{t_0+1} \imm{m}) ^\top +  (\sigma^{2}_{t_0+1})\imm{m}, 		\end{align*}
	where $\bm F_{t_0+1}$ and $\SSigma_{t_0}$ are the same as defined above.}
\end{enumerate}


%\textcolor{blue}{
These calculations can be performed by the \code{LPDS} function, based on a fitted TVP model resulting from a call to \code{shrinkTVP}. The function's arguments are an object of class \code{shrinkTVP} and  \code{data_test},  a data frame with one row, containing covariates and response at time $t_0 + 1$. The following snippet of code fits a \code{shrinkTVP} model to synthetic data up to $T - 1$, and then calculates the LPDS at time $T$. The obtained LPDS score is then displayed.
For an example on how to calculate LPDSs for $k$ points in time, please see~Section~\ref{sec:usmacro}. %}

\begin{figure}[t!]
	\centering
	\includegraphics[width=0.6\linewidth]{./figs/plot_dens.pdf}
	\caption{{One-step ahead predictive density $p(y_{t_0+1}| y_{1},\ldots, y_{t_0}, \xv_{t_0+1}) $ for a synthetic data set. The black vertical line represents the true realisation of $y_{{t_0}+1}$}.}
	\label{fig:pred_dens}
\end{figure}

\begin{CodeChunk}
\begin{CodeInput}
R> res_LPDS <- shrinkTVP(y ~ x1 + x2, data = data[1:(nrow(data) - 1),],
+    display_progress = FALSE)
R> LPDS(res_LPDS, data[nrow(data), ])
\end{CodeInput}
\begin{CodeOutput}
[1] -1.234404
\end{CodeOutput}
\end{CodeChunk}


% \textcolor{blue}{
An additional functionality provided by the package \pkg{shrinkTVP} is the evaluation of the one-step ahead predictive density through the function
%}
\code{eval_pred_dens}.
% \textcolor{blue}{
It takes as inputs an object of class \code{shrinkTVP}, a one row data frame containing $\bm x_{t_0+1}$ and a point, or vector of points, at which the predictive density is to be evaluated. It returns a vector of the same length, containing the value of the density at the points the user supplied. An example of this can be seen in the code below. %}

\begin{CodeChunk}
\begin{CodeInput}
R> eval_pred_dens(1:3, res_LPDS, data[nrow(data), ])
\end{CodeInput}
\begin{CodeOutput}
[1] 0.0003977231 0.0043372607 0.0283156058
\end{CodeOutput}
\end{CodeChunk}


%\textcolor{blue}{
Thanks to its vectorised nature, \code{eval_pred_dens}  %\textcolor{blue}{
can be plugged directly into functions that expect an expression that evaluates to the length of the input, such as the \code{curve} function from the \pkg{graphics} \citep{R} package. The following snippet of code exploits this behaviour to plot the posterior predictive density. The result can be seen in Figure~\ref{fig:pred_dens}. %}


\begin{CodeChunk}
\begin{CodeInput}
R> curve(eval_pred_dens(x, res_LPDS, data[nrow(data), ]), to = 12,
+    ylab = bquote("p(" * y[t[0]+1] * "\uff5c" * y[1] * ","
+    * ldots * "," ~ y[t[0]] * "," ~ x[t[0]+1] * ")"),
+    xlab = expression(y[t[0]+1]), lwd = 2.5, col = "skyblue", axes = FALSE)
R> abline(v = data$y[nrow(data)])
R> axis(1)
R> axis(2)
\end{CodeInput}
\end{CodeChunk}




\section{Predictive exercise: usmacro dataset}
\label{sec:usmacro}


In the following, we provide a brief demonstration on how to use the \pkg{shrinkTVP} package on real data and compare different prior specifications via LPDSs. Specifically, we consider the \code{usmacro.update} dataset from the \pkg{bvarsv} package \citep{kru:bva}. The dataset \code{usmacro.update} contains the inflation rate, unemployment rate and treasury bill interest rate for the United States, from 1953:Q1 to 2015:Q2, \comment{that is $T=250$}. The same dataset up to 2001:Q3 was used by \cite{pri:tim}. The response variable is the inflation rate \code{inf}, while the predictors are the lagged inflation rate \code{inf_lag}, the lagged unemployed rate \code{une_lag} and the lagged treasury bill interest \code{tbi_lag}. We construct our dataset as follows:


\begin{CodeInput}
R> library("bvarsv")
R> data("usmacro.update")
R>
R> # Create matrix of lags and create final data set
R> lags <- usmacro.update[1:(nrow(usmacro.update) - 1), ]
R> colnames(lags) <- paste0(colnames(lags), "_lag")
R> us_data <- data.frame(inf = usmacro.update[2:nrow(usmacro.update), "inf"],
+    lags)
\end{CodeInput}



In the snippet of code below, we run the default \comment{TVP model  with the full hierarchical shrinkage prior} %full hierarchical model
for $60000$ iterations, with a thinning of $10$ and a burn-in of $10000$, hence keeping $5000$ posterior draws.

\begin{CodeInput}
R> us_res <- shrinkTVP(inf ~ inf_lag + une_lag + tbi_lag, us_data,
+    niter = 60000, nburn = 10000, nthin = 10,
+    display_progress = FALSE)
\end{CodeInput}



Once we have fit the model, we can perform posterior inference by using the \code{summary} and \code{plot} methods. The summary is shown below, while Figure~\ref{fig:beta_us} shows the paths of $\bm \beta_t$ evolving over time, and Figure~\ref{fig:theta_us} displays the trace plots (left column) and posterior densities (right column) of $\sqrt{\theta_1}, \ldots, \sqrt{\theta_4}$ obtained via the \code{plot} method.


\begin{CodeChunk}
\begin{CodeInput}
R> summary(us_res, showprior = FALSE)
\end{CodeInput}
\begin{CodeOutput}
Summary of 50000 MCMC draws after burn-in of 10000.

Statistics of posterior draws of parameters (thinning = 10):

 param                   mean    sd       median HPD 2.5% HPD 97.5% ESS
 sigma2                  0.019   0.006    0.018  0.008    0.03      1500.19

 abs(theta_sr_Intercept) 0.141   0.026    0.142  0.09     0.192     832.707
 abs(theta_sr_inf_lag)   0.043   0.006    0.043  0.031    0.056     2328.585
 abs(theta_sr_une_lag)   0.005   0.006    0.002  0        0.017     85.479
 abs(theta_sr_tbi_lag)   0.001   0.002    0      0        0.007     414.397

 beta_mean_Intercept     0.329   0.413    0.126  -0.123   1.192     527.28
 beta_mean_inf_lag       0.748   0.183    0.758  0.384    1.092     1039.996
 beta_mean_une_lag       -0.12   0.073    -0.135 -0.236   0.003     150.987
 beta_mean_tbi_lag       0.009   0.024    0      -0.024   0.074     636.097

 xi2_Intercept           1.125   22.568   0.03   0.001    0.87      5000
 xi2_inf_lag             1.312   75.481   0.004  0        0.22      5000
 xi2_une_lag             0.082   3.564    0      0        0.009     5000
 xi2_tbi_lag             0.01    0.456    0      0        0.001     5000

 a_xi                    0.097   0.043    0.092  0.022    0.181     403.786

 tau2_Intercept          8.933   238.096  0.054  0        9.116     5000
 tau2_inf_lag            40.038  1536.691 0.741  0        22.044    5000
 tau2_une_lag            2.432   55.47    0.032  0        2.527     5000
 tau2_tbi_lag            0.678   28.873   0      0        0.09      5000

 a_tau                   0.1     0.043    0.095  0.02     0.182     483.085

 kappa2                  114.933 254.208  22.776 0        523.803   5000

 lambda2                 9.153   32.301   0.924  0        40.918    5000

 C0                      0.132   0.063    0.122  0.03     0.258     2499.323
\end{CodeOutput}
\end{CodeChunk}

\begin{figure}[t!]
	\centering
	\includegraphics[width=\linewidth]{./figs/plot_usmacro_beta.pdf}
	\caption{Visualization of the evolution of  \comment{the time-varying  parameter  $\bm{\beta}_j=(\bct{j}{0}, \ldots, \bct{j}{T})$ over time  $t=0,\ldots,T$  for $j=1, \ldots, 4$} %$\bm \beta_t$ over time
for the \code{usmacro.update} dataset. The median is displayed as a black line, and the shaded areas indicate the \comment{pointwise} 95\% and 50\% posterior credible intervals.}
	\label{fig:beta_us}
\end{figure}

It appears clear by looking at Figure~\ref{fig:beta_us} that the intercept and the parameter associated with the lagged inflation rate are time-varying, while  the parameters associated with the lagged treasury bill interest rate and the lagged unemployment rate are relatively constant. This can be confirmed by looking at the posterior distributions of the corresponding standard deviations, displayed in Figure~\ref{fig:theta_us}. The posterior densities of the standard deviations associated with the intercept and the lagged inflation are bimodal, with very little mass around zero. This bimodality results from the non-identifiability of the sign of the standard deviation. As a convenient side effect, noticeable bimodality in the density plots of the posterior distribution
  \comment{$p(\sqrt{\theta}_j|\ym)$ of the standard deviations  $\sqrt{\theta}_j$}   % of the standard deviations
  is a strong indication of time variability in the associated parameter  \comment{$\bct{j}{t}$.} % over time $t$.}
  Conversely, the posterior densities of the standard deviations associated with the lagged unemployment and the lagged treasury bill interest rate have a pronounced spike at zero, indicating strong model evidence in favor of constant parameters. Moreover, the path of the parameter of the treasury bill interest rate is centered at zero, indicating that this parameter is neither time-varying nor significant.





In order to compare the predictive performances of different shrinkage priors, we calculate one-step ahead LPDSs for the last 50  points in time for five different prior choices: (1) the full hierarchical \comment{shrinkage} prior, (2) the hierarchical normal-gamma \comment{prior}  with fixed $a^\xi = a^\tau = 0.1$, (3) the normal-gamma \comment{prior} with $a^\xi = a^\tau = 0.1$ and $\kappa^2 = \lambda^2 = 20$, (4) the hierarchical Bayesian Lasso, and (5) the Bayesian Lasso with $\kappa^2 = \lambda^2 = 20$. Figure~\ref{fig:LPDS} shows the cumulative LPDSs for the last 50 quarters of the \code{usmacro.update} dataset. The default prior, that is the fully hierarchical \comment{shrinkage}  prior on both the $\beta_j$'s and the \comment{$\sqrt{\theta_j}$'s}, performs the best in terms of prediction amongst the five fitted models. In Appendix~\ref{sec:multicore} we show how to obtain LPDSs for different models and points in time, using the packages \pkg{foreach} \citep{wes:for} and \pkg{doParallel} \citep{wes:doP}.



\begin{figure}[t!]
	\centering
	\includegraphics[width=\linewidth]{./figs/plot_usmacro_theta.pdf}
	\caption{\comment{Trace plots (left column) and kernel density estimates of the posterior density (right column) for the parameters} $\sqrt{\theta_1}, \ldots, \sqrt{\theta_4}$ for the \code{usmacro.update} dataset.}
	\label{fig:theta_us}
\end{figure}


\begin{figure}[h!]
	\centering
	\includegraphics[width=\linewidth]{./figs/LPDS_macro.pdf}
	\captionof{figure}{Cumulative LPDSs scores for the last 50 quarters of the \code{usmacro.update} dataset, for five different shrinkage priors: (1) the full hierarchical \comment{shrinkage} prior, (2) the hierarchical normal-gamma \comment{prior} with fixed $a^\xi = a^\tau = 0.1$, (3) the normal-gamma
 \comment{prior}  with $a^\xi = a^\tau = 0.1$ and $\kappa^2 = \lambda^2 = 20$, (4) the hierarchical Bayesian Lasso, and (5) the Bayesian Lasso with $\kappa^2 = \lambda^2 = 20$.}
	\label{fig:LPDS}
\end{figure}




\section{Conclusions}
\label{sec:conclusions}

The goal of this  paper was to introduce the reader to the functionality of the \proglang{R} package \pkg{shrinkTVP} \citep{kna-etal:shr}. This \proglang{R} package provides a fully Bayesian approach for statistical inference in TVP models with shrinkage priors.
\comment{On the one hand, the package provides an easy entry point for users
who %aim at applying TVP models and
want  to pass on only their data in a first step of exploring  TVP models for  their specific application context.
Running the function  \code{shrinkTVP} under  the default model with a fully hierarchical shrinkage prior with predefined hyperparameters,  estimation of a TVP model becomes as easy  as using the well-known function \code{lm} for a standard linear regression model.
  On the other hand, exploiting numerous advanced options of the package,
  the more experienced user can also explore alternative model specifications such as the Bayesian Lasso and use log-predictive density scores to compare
  various model specifications.}


Various examples of the usage of \pkg{shrinkTVP} were given, and the \code{summary} and \code{plot} methods for
straightforward posterior inference were illustrated.
Furthermore, a predictive exercise with the dataset \code{usmacro.updade} from the package \pkg{bvarsv} was performed, with a focus on the calculation of LPDSs using \code{shrinkTVP}. The default model in \pkg{shrinkTVP} showed better performance than its competitors in terms of cumulative LPDSs.
\comment{While these  examples  were confined to univariate responses,  the package can also be applied in a multivariate context, for instance to the
sparse TVP Cholesky SV model considered in \citet{bit-fru:ach}, exploiting a representation of this model as a system of independent  TVP models
with univariate responses.}

\comment{A plan for further versions of the package is to implement additional shrinkage priors for TVP models such as the well-known horseshoe prior
\citep{bha-etal:las}}.


%\clearpage\newpage

\appendix

\section{Appendix: Full conditional distribution of the latent states}
\label{sec:mccau}
Let $y^\star_t = y_t - \bm x_t \betav$  and  $\Fm_t=\bm x_{t}\Diag{\sqrt{\theta_{1}},\dots,\sqrt{\theta_{d}}}$ for $t=1,\dots,T$. Conditional on all other variables,  the joint density for the state process
   $ \tilde \betav =( \tilde \betav_0, \tilde \betav_1, \dots, \tilde \betav_T)$ is multivariate normal.
This distribution can be written in
   terms of the tri-diagonal precision matrix $\OOmega$  and the mean vector $\cm$ \citep{mcc-etal:sim}:
\begin{eqnarray} \label{postbetav}
    \tilde \betav  | \betav, \QQ ,  \sigma^2_1, \ldots,  \sigma^2_T, y^\star_1, \ldots  y^\star_T
   \sim \Normult{(T+1)  d }{\OOmega^{-1}\cm,\OOmega^{-1} }
\end{eqnarray}
   where:
   \begin{equation*}
\OOmega  =   \begin{bmatrix}
\OOmega_{00}&\OOmega_{01}  & 0 & & \\
\OOmega^{\top}_{01} &  \OOmega_{11}&\OOmega_{12} &0&0&\\
0 &  \OOmega^{\top}_{12}&\OOmega_{22}&\OOmega_{23}&\ddots&\vdots\\
&  0&\OOmega^{\top}_{23}&\ddots&\ddots&0 \\
&  \vdots & \ddots & \ddots &\OOmega_{T-1,T-1}&\OOmega_{T-1,T} \\
&  0&\ldots &0 &\OOmega^{\top}_{T-1,T}&\OOmega_{TT} \\
  \end{bmatrix},  \quad
 \cm  =  \begin{bmatrix}
  \cm_0\\
  \cm_1\\
 \cm_2 \\
  \vdots \\
  \cm_T
  \end{bmatrix} .
\end{equation*}
In this representation, each submatrix $\OOmega_{ts}$ is a matrix of dimension $d \times d$  defined as
\begin{eqnarray*}
 \OOmega_{00} &=&  2 I_d,\\
%  \OOmega_{11} &\equiv& \mathbf{v}^{\top}_1 A_{1}  \mathbf{v}_1 + I_d + P^{-1}_0,\\
  \OOmega_{tt} &=& \Fm^{\top}_t   \Fm _t /\sigma^2_t +   2 I_d,\quad t= 1,\dots, T-1,\\
%  \OOmega_{tt} &\equiv& \mathbf{v}^{\top}_t  A_{11,t} \mathbf{v}_t +  A_{22,t} + A_{22,t-1},\quad t= 1,\dots, T-1,\\
   % \OOmega_{TT} &\equiv& \mathbf{v}^{\top}_T  A_{11,T}   \mathbf{v}_T  + A_{22,T},\\
      \OOmega_{TT} &=& \Fm^{\top}_T   \Fm_T /\sigma^2_T  + I_d,\\
   \OOmega_{t-1,t} &=& -  I_d, \quad t=1, \dots, T,
\end{eqnarray*}
where $I_d$ is the $d\times d$ identity matrix and
 $\cm_{t}$ is a  column vector of  dimension $d \times 1$, defined as
\begin{eqnarray*}
\cm_{0}  =  \bfz, \qquad
   \cm_{t}  =  (\Fm^{\top}_t /\sigma^2_t) y_t^{\star},	\quad t = 1,\dots, T .
\end{eqnarray*}
In the homoscedastic case, $\sigma^2_1=\ldots = \sigma^2_T=\sigma^2$.




\section{Multicore LPDS calculation}
\label{sec:multicore}

In the code below, the following \proglang{R} packages are used:
\pkg{doParallel} \citep{wes:doP}, \pkg{foreach} \citep{wes:for},
\pkg{zoo} \citep{zei:zoo}, and \pkg{RhpcBLASctl} \citep{Nak:Rhp}.

%\textcolor{green}{ABN: i have not found a zoo manual to cite?}



\begin{CodeInput}
R> # Calculate LPDS in multicore
R> # Load libraries for multicore computations
R> library("doParallel")
R> library("foreach")
R>
R> # For manipulating dates
R> library("zoo")
R>
R> # Load library for controlling number of BLAS threads
R> library("RhpcBLASctl")
R>
R> # Define how many periods to calculate LPDS for
R> Tmax <- nrow(us_data) - 1
R> T0 <- Tmax - 49
R>
R> # Determine number of cores to be used and register parallel backend
R> ncores <- 4
R> cl <- makeCluster(ncores)
R> registerDoParallel(cl)
R>
R> lpds <- foreach(t = T0:Tmax, .combine = "cbind",
+    .packages = c("RhpcBLASctl", "shrinkTVP")) %dopar% {
+
+      niter <- 60000
+      nburn <- 30000
+      nthin <- 10
+    # Set number of BLAS threads, so they don't interfere with each other
+    blas_set_num_threads(1)
+
+    # Create data_t from all data up to time t and
+    # y_test and x_test from data at time t+1
+    data_test <- us_data[t+1,]
+    data_t <- us_data[1:t,]
+
+    # Run MCMC to calculate all LPDS
+    res_base <- shrinkTVP(inf ~ inf_lag + une_lag + tbi_lag, data = data_t,
+      niter = niter, nburn = nburn, nthin = nthin,
+      hyperprior_param = list(nu_tau = 1, nu_xi = 1))
+
+    res_las_hier <- shrinkTVP(inf ~ inf_lag + une_lag + tbi_lag, data = data_t,
+      niter = niter, nburn = nburn, nthin = nthin,
+      learn_a_xi = FALSE, learn_a_tau = FALSE,
+      a_xi = 1, a_tau = 1)
+
+    res_las <- shrinkTVP(inf ~ inf_lag + une_lag + tbi_lag, data = data_t,
+      niter = niter, nburn = nburn, nthin = nthin,
+      learn_a_xi = FALSE, learn_a_tau = FALSE,
+      a_xi = 1, a_tau = 1,
+      learn_kappa2 = FALSE, learn_lambda2 = FALSE,
+      kappa2 = 20, lambda2 = 20)
+
+      res_ng_hier <- shrinkTVP(inf ~ inf_lag + une_lag + tbi_lag, data = data_t,
+        niter = niter, nburn = nburn, nthin = nthin,
+        learn_a_xi = FALSE, learn_a_tau = FALSE,
+        a_xi = 0.1, a_tau = 0.1)
+
+      res_ng <- shrinkTVP(inf ~ inf_lag + une_lag + tbi_lag, data = data_t,
+        niter = niter, nburn = nburn, nthin = nthin,
+        learn_a_xi = FALSE, learn_a_tau = FALSE,
+        a_xi = 0.1, a_tau = 0.1,
+        learn_kappa2 = FALSE, learn_lambda2 = FALSE)
+
+      lpds_res <- c(LPDS(res_base, data_test),
+                    LPDS(res_ng_hier, data_test),
+                    LPDS(res_ng, data_test),
+                    LPDS(res_las_hier, data_test),
+                    LPDS(res_las, data_test))
+
+      rm("res_base", "res_ng_hier",  "res_ng", "res_las_hier", "res_las")
+
+      return(lpds_res)
+  }
R> stopCluster(cl)
R>
R> cumu_lpds <- apply(lpds, 1, cumsum)
R> color <- brewer.pal(5, "Set2")
R> # Plot results
R> par(mar=c(6,4,1,1))
R> colnames(cumu_lpds) <- c("Default", "Hierarchical Normal Gamma",
+    "Non-hierarchical Normal Gamma", "Hierarchical Lasso",
+    "Non-hierarchical Lasso")
R>
R> matplot(cumu_lpds, type = "l", ylab = "Cumulative LPDS",
+    xaxt = "n", xlab = "", col = color, lty=1, lwd=2 )
R>
R> # Extract labels from time series
R> labs = as.yearmon(time(usmacro.update))[T0:Tmax + 1][c(FALSE, TRUE)]
R>
R> # Create custom axis labels
R> axis(1, at = (1:length(T0:Tmax))[c(FALSE, TRUE)], labels = FALSE)
R> text(x=(1:length(T0:Tmax))[c(FALSE, TRUE)],
+    y=par()$usr[3]-0.05*(par()$usr[4]-par()$usr[3]),
+    labels=labs, srt=45, adj=1, xpd=TRUE)
R>
R> # Add legend
R> legend("topright", colnames(cumu_lpds), col = color,
+    lty = 1,lwd = 2, bty = "n", cex = 0.8)
\end{CodeInput}


\bibliography{shrinkTVP}
%% -----------------------------------------------------------------------------

\end{document}
