Extracted main text — title through conclusion, appendix excluded. This is what our citation measures are computed over, published so the extraction can be checked by eye.
63,238 characters · 17 sections · 73 citation commands
Shrinkage in the Time-Varying Parameter Model Framework Using the R Package shrinkTVP
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 parameters 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 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 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. fru-wag:sto employ a spike and slab prior, while continuous shrinkage priors have been utilised as a regularization alternative in, e.g., bel-etal:hie_tv, bit-fru:ach and cad-etal:tri. For an excellent review of shrinkage priors, with a particular focus on high dimensional regression, the reader is directed to bha-etal:las.
In this paper, we describe the \proglang{R} package \pkg{shrinkTVP} 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 bit-fru:ach and cad-etal:tri. The computationally intensive Markov chain Monte Carlo (MCMC) algorithms in the package are written in \proglang{C++} and interfaced with \proglang{R} R via the \pkg{Rcpp} edd-bal:ext and the \pkg{RcppArmadillo} 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++}.
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} zei:zoo and \pkg{xts} xts as well as time series formats like \code{ts} are supported. Estimation output is compatible with the popular \proglang{R} package \pkg{coda} 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.
\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} 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 pet-etal:dyn introduces the methodology and many \proglang{R} code examples. As of now, priors are not designed to encourage shrinkage and \pkg{shrinkTVP} complements \pkg{dlm} in this regard.
The \proglang{R} package \pkg{bvarsv} 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 pri:tim. We refer to del-pro:tim2 for details on the MCMC algorithm and a later correction of the original scheme. 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 pri:tim as example data set, \code{usmacro.update}, which we use in our predictive exercise in Section (ref) to showcast the effect of introducing shrinkage priors on time-varying parameters.
Additional packages emerged very recently. The \proglang{R} package \pkg{tvReg} 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. For an illustrative application, a TVP-VAR(4) model is fitted to the \code{usmacro.update} data set mentioned above, using the function \code{tvVAR}. The \proglang{R} package \pkg{walker} walker facilitates the estimation of DLMs and generalized DLMs using MCMC algorithms provided by \proglang{Stan} carpenter2017stan. For inference, the importance sampling method of ISMCMC is implemented within a Hamiltonian Monte Carlo framework. The \proglang{R} package \pkg{bsts} 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}.
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.
Among others, \pkg{ncvreg} ncvreg is useful for fitting standard penalized regression estimators, \pkg{glmnet} sim-etal:reg allows elastic-net regularization for a variety of models, \pkg{horseshoe} pas:hor implements the horseshoe prior, while \pkg{shrink} dun-etal:shr provides various shrinkage methods for linear, generalized linear, and Cox regression models. \pkg{biglasso} zen:big aims at very fast lasso-type regularization for high-dimensional linear regression. Recent \proglang{R} packages include \pkg{NormalBetaPrime} bai:nor for Bayesian univariate and \pkg{MBSP} 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} 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.
The remainder of the paper is organized as follows. Section (ref) briefly introduces TVP models and normal-gamma-gamma shrinkage priors, and describes the MCMC algorithms for posterior simulation. The package \pkg{shrinkTVP} is introduced in Section (ref). 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) 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. Finally, Section (ref) concludes the paper.
Let us recall the state space form of a TVP model. For $t = 1, \ldots, T$, we have that
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 $\mathbf{Q} = \text{Diag}(\theta_1, \ldots, \theta_d)$ 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, \mathbf{Q})$, with initial mean $\bm \beta = (\beta_1, \ldots, \beta_d)$. Model ((ref)) can be rewritten equivalently in the non-centered parametrization as
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 (SV) specification. In the latter case, the log-volatility $h_t = \log \sigma^2_t$ follows an AR(1) model jac-etal:bayJBES,kas-fru:anc, kas:dea. More specifically,
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 nak:tim, sim:com by capturing some of the variability in the error term.
We place conditionally independent normal-gamma-gamma (NGG) priors cad-etal:tri,gri-bro:hie, 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 triple gamma prior on the innovation variances $\theta_j$, for $j = 1, \ldots, d$. The NGG can be represented as a conditionally normal distribution, where the component specific variance is itself a compound probability distribution resulting from two gamma distributions. In this representation, it looks as follows
Letting $c^\xi$ and $c^\tau$ go to infinity results in a normal-gamma (NG) prior gri-bro:inf on the $\sqrt{\theta}_j$'s and $\beta_j$'s. It has a representation as a conditionally normal distribution, with the component specific variance following a gamma distribution, that is
From here, letting $a^\xi$ and $a^\tau$ go to infinity yields a normal prior with fixed variance, also known as ridge regression:
We refer to $a^\xi$ and $a^\tau$ as the pole parameters, as marginally more mass is placed around zero as they become smaller. $c^\xi$ and $c^\tau$ are referred to as the tail parameters, as they control the amount of mass in the tails of the distribution, with smaller values equating to heavier tails. Finally, the parameters $\kappa_B^2$ and $\lambda_B^2$ are dubbed the {global shrinkage parameters}, as they influence how strongly {all} parameters are pulled to zero. The larger $\kappa_B^2$ and $\lambda_B^2$, the stronger this effect.
One of the key benefits of the NGG prior is that many interesting shrinkage priors are contained within it as special or limiting cases. Beyond the NG prior mentioned above, two such cases are the horseshoe prior car-etal:han and the Bayesian Lasso par-cas:bay. The former results from an NGG prior with the pole and tail parameters equal to $0.5$, while the latter is a special case of the NG prior with a pole parameter fixed to one. As the connection between the NGG prior and the horseshoe prior may not be entirely obvious from the parameterization presented here, the interested reader is referred to cad-etal:tri for details.
The parameters $a^\xi$, $a^\tau$, $c^\xi$, $c^\tau$, $\kappa_B^2$ and $\lambda_B^2$ can be learned from the data through appropriate prior distributions. Results from cad-etal:tri motivate the use of different distributions for these parameters under the NGG and NG prior. In the NGG case, the scaled global shrinkage parameters conditionally follow F distributions, depending on their respective pole and tail parameters:
The scaled tail and pole parameters, in turn, follow beta distributions:
These priors are chosen as they imply a uniform prior on a suitably defined model size, see cad-etal:tri for details. In the NG case the global shrinkage parameters follow independent gamma distributions:
In order to learn the pole parameters in the NG case, we generalize the approach taken in bit-fru:ach and place the following gamma distributions as priors:
which correspond to the exponential priors used in bit-fru:ach when $\alpha_{a^\xi}=1$ and $\alpha_{a^\tau}=1$. The parameters $\alpha_{a^\xi}$ and $\alpha_{a^\tau}$ act as degrees of freedom and allow the prior to be bounded away from zero.
In the homoscedastic case we employ a hierarchical prior, where the scale of an inverse gamma prior for $\sigma^2$ follows a gamma distribution, that is,
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 (ref) are chosen as in kas-fru:anc, that is
with hyperparameters $b_\mu, B_\mu, a_\phi, b_\phi,$ and $B_\sigma$.
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 bit-fru:ach and cad-etal:tri for further details.
Step 4 presents a fork in the algorithm, as different parameterizations are used in the NGG and NG case, as to improve mixing. For details on the exact parameterization used in the NGG case, see cad-etal:tri. Additionally, not all sampling steps are performed in all prior setups. If, for example, the user has defined that $\kappa_B^2$ should not be learned from the data, then this step is not executed.
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). We employ the procedure described in mcc-etal:sim which exploits the sparse, block tri-diagonal structure of the precision matrix of the full conditional distribution of $\tilde \boldsymbol{\beta} =( \tilde \boldsymbol{\beta}_0, \ldots, \tilde \boldsymbol{\beta}_T)$, to speed up computations.
Moreover, as described in bit-fru:ach, in step 3 we make use of the ancillarity-sufficiency interweaving strategy (ASIS) introduced by 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 kas-fru:anc, multivariate factor SV models kas-etal:eff and dynamic linear state space models sim-etal:int.
\paragraph{Adaptive Metropolis-within-Gibbs} For the pole and tail parameters, no full conditionals exist and a Metropolis-Hastings step has to be performed. To improve mixing, \code{shrinkTVP} supports adaptive Metropolis-within-Gibbs as in rob-ros:exa. The algorithm works as follows. For each parameter $i$ that is being learned from the data, let $s_i$ represent the standard deviation of the proposal distribution. After the $n^{th}_i$ batch of $m_i$ iterations, update $s_i$ according to the following rule:
The starting value of $s_i$, $m_i$, $c_i$ and $d_i$ can all be set by the user. Additionally, if adaptive Metropolis-within-Gibbs is not desired, it can be switched off and a simple Metropolis-Hastings step will be performed.
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 Equation (ref) is estimated in a Bayesian fashion with the NG prior defined in equations (ref), (ref), (ref) and (ref) with the following choice for the hyperparameters: $d_1 = d_2 = e_1 = e_2 = 0.001$, $\alpha_{a^\xi}=\alpha_{a^\tau}=5$ and $\beta_{a^\xi}=\beta_{a^\tau}=10$, implying a prior mean of $\mbox{\rm E}(a^\xi)= \mbox{\rm E}(a^\tau)=0.1$. The error is assumed to be homoscedastic, with prior defined in Equation (ref) 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} 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.
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 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}.
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.
More granular control over the prior setup can be exercised by passing additional arguments to \code{shrinkTVP}. The most important argument in this regard is \code{mod_type}, which is used to specify whether the normal-gamma-gamma (\code{mod_type = "triple"}), the normal-gamma (\code{mod_type = "double"}) or ridge regression (\code{mod_type = "ridge"}) is used. Beyond this, the user can specify the hyperparameters given in Section (ref) and has the possibility to fix one or both of the values of the global shrinkage parameters ($\kappa_B^2$, $\lambda_B^2$) and the pole and tail parameters ($a^\tau$, $a^\xi$, $c^\tau$, $c^\xi$). By default, these parameters are learned from the data. The benefit of this flexibility 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. Under an NGG prior, for example, setting the pole and tail parameters equal to $1/2$ results in a horseshoe prior on the $\sqrt\theta_j$'s and the $\beta_j$'s, respectively. If the user desires a higher degree of sparsity, this can be achieved by setting the pole parameters to a value closer to zero. Table (ref) gives an overview of different model specifications. Note that different hyperparameter values can be chosen for the variances and the means of the initial values.
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. If the learning of a parameter is deactivated and no specific fixed value is provided, \code{shrinkTVP} will resort to default values. These equate to 0.1 for the pole and tail parameters and 20 for the global shrinkage parameters. 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 pole parameters}
It is possible to set the pole parameter $a^\xi$($a^\tau$) 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:
\paragraph{Fixing the global shrinkage parameters} The user can choose to fix the value of $\kappa_B^2$($\lambda_B^2$) by specifying the argument \code{kappa2_B} (\code{lambda2_B}), after setting \code{learn_kappa2_B} (\code{learn_lambda2_B}) 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$:
\paragraph{Changing the prior type} To change the model type, the input argument \code{mod_type} has to be supplied. It has to be a string equal to either \code{"triple"}, \code{"double"} or \code{"ridge"}. As an example, we fit a hierarchical NGG prior, both on the $\sqrt{\theta_j}$'s and on the $\beta_j$'s:
\paragraph{Fixing the tail parameters} Much like the pole parameters, the tail parameter $c^\xi$ ($c^\tau$) can also be fixed to a value. This is done by setting \code{learn_c_xi} (\code{learn_c_tau}) to \code{FALSE} and then supplying the input parameter \code{c_xi} (\code{c_tau}). As an example, the code below fits a non-hierarchical horseshoe prior, both on the $\sqrt{\theta_j}$'s and on the $\beta_j$'s:
The stochastic volatility specification defined in Equation (ref) 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 with an NG prior in which all the parameters are learned and the observation equation errors are modeled through stochastic volatility:
The priors on the SV parameters are the ones defined in Equation (ref), with hyperparameters fixed to $b_\mu = 0$ , $B_\mu = 1$, $a_\phi = 5$, $b_\phi = 1.5$ , and $B_\sigma = 1$.
Beyond simply switching off parts of the hierarchical structure of the prior setup, users can also modify the hyperparameters governing the hyperprior distributions. This can be done through the arguments \code{hyperprior_param} and \code{sv_param}, which both have to be named lists. Hyperparameters not specified by the user will be set to default values, which can be found in the help file of the \code{shrinkTVP} function. Note, however, that the dependence structure (e.g. $\kappa^2_B$ depends on $a^\xi$ and $c^\xi$ in the NGG specification) can not be changed. As such, if the user desires to change the hyperparameters of a prior that depends on other parameters, this can only be achieved by deactivating the learning of the parameters higher up in the hierarchy and fixing them to specific values. To demonstrate how to change specific hyperparameters, the code below modifies those governing the prior on $a^\xi$:
The Metropolis-Hastings algorithm discussed in Section (ref) can be tuned via the argument \code{MH_tuning}. Similar to \code{hyperprior_param} and \code{sv_param}, it is a named list where values that are not supplied are replaced by standard values. By default, adaptive Metropolis-within-Gibbs is activated for all parameters learned from the data that requrire a Metropolis-Hastings step. Below is an example where the adaptive Metropolis is deactivated for one of the pole parameters and slightly tuned for the other:
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 NG prior, the values are:
When some parameters are fixed by the user, the corresponding output value is omitted. Additionally, increasing or decreasing the amount of levels in the hierarchy of the prior also changes which values are returned. For example, if \code{mod_type} is changed to \code{"triple"} and the learning of the tail parameters $c^\xi$ and $c^\tau$ is not deactivated, then the output will also contain the respective parameter draws in \code{c_xi} and \code{c_tau}. 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 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% highest posterior density region and effective sample size (ESS) for the non time-varying parameters.
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) for the corresponding output.
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 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).
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 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 method from this package that users are familiar with (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), where \code{pars = "theta_sr"}, changes 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), one can run
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 log of the posterior 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
where $\bm \psi$ is the set of model parameters 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})$. 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)). However, bit-fru:ach propose a more efficient approximation of the predictive density, 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$.homoscedastic error case, given $M$ samples from the posterior distribution of the parameters and the latent variables up to $t_0$, Monte Carlo integration of the resulting predictive density yields following mixture approximation,
where the conditional predictive densities are Gaussian and the conditional moments depend on the MCMC draws. The mean $\hat{y} _{t_0+1} ^{(m)}$ and the variance $S _ {t_0+1} ^{(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 $\boldsymbol{\Sigma}_{t_0}$ of the posterior distribution of $\tilde {\bm \beta}_{t_0}$. These quantities can be obtained by iteratively calculating $\boldsymbol{\Sigma}_{t}$ and $\bm m_{t}$ up to time $t_0$, as described in mcc-etal:sim:
The quantities $\bm c_t$, $\boldsymbol \Omega_{tt}$ and $\boldsymbol \Omega_{t-1,t}$ for $t=1, \ldots, t_0$ are given in Appendix (ref).
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 (ref), and used to calculate the 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:
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).
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}.
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.
Thanks to its vectorised nature, \code{eval_pred_dens} 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} R package. The following snippet of code exploits this behaviour to plot the posterior predictive density. The result can be seen in Figure (ref).
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 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, that is $T=250$. The same dataset up to 2001:Q3 was used by 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:
In the snippet of code below, we estimate a TVP model with a fully hierarchical NG prior for $60000$ iterations, with a thinning of $10$ and a burn-in of $10000$, hence keeping $5000$ posterior draws.
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) shows the paths of $\bm \beta_t$ evolving over time, and Figure (ref) 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.
It appears clear by looking at Figure (ref) 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). 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 $p(\sqrt{\theta}_j|{\mathbf y})$ of the standard deviations $\sqrt{\theta}_j$ is a strong indication of time variability in the associated parameter $\beta_{{j}{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 eleven different prior choices: (1) the full hierarchical NGG prior, (2) the hierarchical NGG prior with fixed $a^\xi = a^\tau = c^\xi = c^\tau = 0.1$, (3) the NGG prior with $a^\xi = a^\tau = c^\xi = c^\tau = 0.1$ and $\kappa_B^2 = \lambda_B^2 = 20$, (4) the hierarchical horseshoe prior, (5) the horseshoe prior $\kappa_B^2 = \lambda_B^2 = 20$, (6) the full hierarchical NG prior, (7) the hierarchical NG prior with fixed $a^\xi = a^\tau = 0.1$, (8) the NG prior with $a^\xi = a^\tau = 0.1$ and $\kappa_B^2 = \lambda_B^2 = 20$, (9) the hierarchical Bayesian Lasso, and (10) the Bayesian Lasso with $\kappa_B^2 = \lambda_B^2 = 20$ and (11) ridge regression with $\kappa_B^2 = \lambda_B^2 = 20$. Figure (ref) shows the cumulative LPDSs for the last 50 quarters of the \code{usmacro.update} dataset. The default prior, the fully hierarchical NG prior on both the $\beta_j$'s and the $\sqrt{\theta_j}$'s, performs the best in terms of prediction. In Appendix (ref) we show how to obtain LPDSs for different models and points in time, using the packages \pkg{foreach} wes:for and \pkg{doParallel} wes:doP.
The goal of this paper was to introduce the reader to the functionality of the \proglang{R} package \pkg{shrinkTVP} kna-etal:shr. This \proglang{R} package provides a fully Bayesian approach for statistical inference in TVP models with shrinkage priors. On the one hand, the package provides an easy entry point for users who 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 NG 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 or the horseshoe prior 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. 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 bit-fru:ach, exploiting a representation of this model as a system of independent TVP models with univariate responses.