Get 20M+ Full-Text Papers For Less Than $1.50/day. Start a 14-Day Trial for You or Your Team.

Learn More →

Generalized Autoregressive Score Models in R: The GAS Package

Generalized Autoregressive Score Models in R: The GAS Package This paper presents the R package GAS for the analysis of time series under the Generalized Autoregressive Score (GAS) framework of Creal et al. (2013) and Harvey (2013). The distinctive feature of the GAS approach is the use of the score function as the driver of time{variation in the parameters of nonlinear models. The GAS package provides functions to simulate univariate and multivariate GAS processes, estimate the GAS parameters and to make time series forecasts. We illustrate the use of the GAS package with a detailed case study on estimating the time{varying conditional densities of a set of nancial assets. Keywords : GAS, time series models, score models, dynamic conditional score, R software. arXiv:1609.02354v1 [stat.CO] 8 Sep 2016 2 The R package GAS 1. Introduction Time{variation in the parameters describing a stochastic time series process is pervasive in al- most all applied scienti c elds. Early references to time series models include Kalman (1960) and Box and Jenkins (1970). In many settings, the model of interest is characterized by time{ varying parameters, for which the literature has proposed a myriad of possible speci cations. Recently, Creal et al. (2013) and Harvey (2013) note that many of the proposed models are either dicult to estimate (in particular, the class of stochastic volatility models reviewed in Shephard (2005)) and/or do not properly take the shape of the conditional distribution of the data into account. Creal et al. (2013) and Harvey (2013) therefore propose to use the score of the conditional density function as the main driver of time{variation in the param- eters of the time series process used to describe the data. A further advantage of using the conditional score as driver is that the estimation by Maximum Likelihood is straightforward. The resulting model is referred to as: Score{Driven model, Dynamic Conditional Score (DCS) model, or Generalized Autoregressive Score (GAS) model. In this article and accompanying R package, we use the GAS acronym for both. The R package GAS is conceived to be of relevance for the modelling of all types of time series data. It does not matter whether they are real{valued, integer{valued, (0,1){bounded or strictly positive, as long as there is a conditional density for which the score function and the Hessian are well{de ned. The practical relevance of the GAS framework has been illustrated in the case of nancial risk forecasting (see e.g., Harvey and Sucarrat (2014) for market risk, Oh and Patton (2013) for systematic risk, and Creal et al. (2014) for credit risk analysis), dependence modelling (see e.g., Harvey and Thiele (2015) and Janus et al. (2014)), and spatial econometrics (see e.g., Blasques et al. (2014d) and Catania and Bill e (2016)). For a more complete overview of the work on GAS models, we refer the reader to the GAS community page at http://www.gasmodel.com/. It is important to note that, even though the GAS framework has been developed by econome- tricians, it is exible enough to be used in all elds in which the use of time{varying parameter models is relevant. The main diculty in using GAS models is to derive the score and Hessian and implementing the Maximum Likelihood estimation of the resulting nonlinear models. The R package GAS answers these needs by proposing an integrated set of R functions to do time series analysis in the R statistical language (R Core Team 2016) under the GAS framework. The functionalities include: (i) estimation, (ii) prediction, (iii) simulation, (iv) backtesting, and (v) graphical representation of the results, implying that it is ready to use in real{life ap- plications. The user interface uses the R programming language, which has the advantage of being free and open source. However, most of the underlying routines are principally written in C++ exploiting the armadillo library (Sanderson 2010) and the R packages Rcpp (Ed- delbuettel and Fran cois 2011; Eddelbuettel et al. 2016a) and RcppArmadillo (Eddelbuettel and Sanderson 2014; Eddelbuettel et al. 2016b) to speed up the computations. Furthermore, since the package is written with the S4 methods, R users with basic programming knowledge will nd common functions such as coef(), plot() and show() to extract and analyze their A typical example is the class of (G)ARCH models in which the squared (demeaned) return is the driver of time{variation in the conditional variance, independently of the shape of the conditional distribution of the return. To see that this is counter{intuitive, consider the case of observing a 10% return when the conditional mean is 0% and the volatility is 3%. Under the assumption of a normal distribution, the 10% return is a strong signal of an increase in volatility, while under a fat{tailed Student{t distribution, the signal is weakened because of the higher probability that the extreme value is an observation from the tails. David Ardia, Kris Boudt, Leopoldo Catania 3 results. We believe that this aspect is of primary importance since it dramatically increases the number of potential users. The R package GAS is available from the CRAN repository at https://cran.r-project.org/package=GAS. Other codes available for speci c GAS models are available in the GAS community page at http://www.gasmodel.com/. For instance, the R package betategarch (Sucarrat 2013) allows us to estimate the beta{t{EGARCH model of Harvey (2013) and its skewed version introduced by Harvey and Sucarrat (2014). The outline of the paper is as follows: Section 2 reviews the GAS framework to de ne time{ varying parameter models, referring to the seminal works of Creal et al. (2013) and Harvey (2013). Section 3 introduces the R package GAS and illustrates how to to simulate, estimate and make predictions. Section 4 presents a real{life application to nancial data. Section 5 concludes. 2. The GAS framework to modeling time{varying parameters One of the most appealing characteristics of the GAS framework is its applicability to de ne time{varying parameter models in a large variety of univariate and multivariate time series settings. We try to be as general as possible in reviewing the GAS framework, and report in Appendix A the detailed equations for the speci c case of a conditionally Student{t distributed random variable. In this section, we rst introduce the notation and present the GAS model when the parameter space is unrestricted. We then show how a mapping function can be used to model the time{variation in the parameters when the parameter space is restricted. The section concludes by summarizing the Maximum Likelihood approach for GAS model estimation. 2.1. Model speci cation Let y 2 < be an N {dimensional random vector at time t with conditional distribution: y jy  p(y ;  ) ; (1) t 1:t1 t t 0 0 0 J where y  (y ; : : : ; y ) contains the past values of y up to time t 1 and  2   < 1:t1 1 t1 t t is a vector of time{varying parameters which fully characterizes p() and only depends on y and a set of static additional parameters , i.e.,   (y ; ) for all t. The main 1:t1 t 1:t1 feature of GAS models is that the evolution in the time{varying parameter vector  is driven by the score of the conditional distribution de ned in (1), together with an autoregressive component: + A s + B  ; (2) t+1 t t where, , A and B are matrices of coecients with proper dimensions collected in , and s is a vector which is proportional to the score of (1): s  S ( )r (y ;  ) : t t t t t t The matrix S is a J  J positive de nite scaling matrix known at time t and: @ log p(y ;  ) t t r (y ;  )  ; t t t t 4 The R package GAS is the score of (1) evaluated at  . Creal et al. (2013) suggest to set the scaling matrix S to t t a power > 0 of the inverse of the Information Matrix of  to account for the variance of r . More precisely: S ( )  I ( ) ; t t t t with: I ( )  E r (y ;  )r (y ;  ) ; (3) t t t1 t t t t t t where the expectation is taken with respect to the conditional distribution of y given y . t 1:t1 The additional parameter is xed by the user and usually takes value in the set f0; ; 1g. When = 0, S = I and there is no scaling. If = 1 ( = ), the conditional score r (y ;  ) t t t t is premultiplied by the inverse of (the square root of ) its covariance matrix I ( ). t t It is worth noting that, whatever the choice of , s is a Martingale Di erence (MD) with respect to the distribution of y given y , i.e., E [s ] = 0 for all t. Furthermore, when t 1:t1 t1 t = , the additional moment condition V [s ] = I can be easily derived. Due to the t1 t fact that s is a MD, if the spectral radius of B is less then one , the updating equation of  reported in (2) implies a mean reverting process for  through the long{term mean t t 1 1 (I B) , which means that the unconditional value of  is (I B) . It follows that, the J {valued vector  and the J J matrix B control for the level and the persistence of the process, respectively. The additional J  J matrix of coecients A, that premultiplies the scaled score s , controls for the impact of s to  . Speci cally, as detailed in Creal et al. (2013), the quantity t t+1 s indicates the direction to update the vector of parameters from  , to  , acting as t t t+1 a steepest ascent algorithm for improving the model local t given the current parameter position. Interestingly, this updating procedure resembles the well{known Newton{Raphson algorithm. Hence, A can be interpreted as the step of the update, and needs to be designed in a way to do not distort the signal coming from s ; see Section 2.3. 2.2. Reparametrization In (2) the parameter vector  has a linear speci cation and is thus unbounded. In practice, the parameter space of  is often restricted (  < ). For instance, when we model the scale parameter of a Student{t distribution, we need to ensure its positiveness. Even if this problem can be solved by imposing constraints on  (as is done in the GARCH model, see, Bollerslev 1986), the standard solution under the GAS framework is to use a (possibly J J e e nonlinear) link function () that maps  2 < into  and where  2 < has the linear t t t dynamic speci cation of (2). Speci cally, let  : < !  be a twice di erentiable vector{ valued mapping function such that ( ) =  . The updating equation for  is then given t t t by: ( ) t t (4) e e + Ae s + B ; t t t1 e e e e e e e where s  S ( )r (y ;  ) and r (y ;  ) represents the score of (1) with respect to  , and, t t t t t t t t t t e e e e e consequently, S ( ) can depend on the information matrix of  given by I ( ). Denote the t t t t t We denote by I the identity matrix of appropriate size. The spectral radius of a L  L matrix X is de ned as  (X)  max (j j; : : : ;j j), where  is the i{th 1 L i eigenvalue of X, for i = 1; : : : ; L. David Ardia, Kris Boudt, Leopoldo Catania 5 Jacobian matrix of () evaluated at  as follows: @( ) J ( )  : Then, the following relations hold: e e e r (y ;  ) = J ( ) r (y ;  ) t t t t t t t e e e e I ( ) = J ( ) I ( )J ( ) : t t t t t t This way, almost all the nonlinear constraints can be easily handled via the de nition of a proper mapping function () and its associated Jacobian matrix J (). The coecients to be estimated are gathered into   (; A; B) and estimated by numerically maximizing the (log-)likelihood function as detailed in Section 2.3. In Appendix B we discuss the choice of mapping function for GAS models in more details. 2.3. Maximum likelihood estimation A useful property of GAS models is that, given the past information and the static parameter vector , the vector of time{varying parameters,  , is perfectly predictable and the loglikeli- hood function can be easily evaluated via the prediction error decomposition. More precisely, for a sample of T realizations of y , collected in y , the vector of parameters  can be t 1:T estimated by Maximum Likelihood (ML) as the solution of: arg max L (; y ) ; (5) 1:T where: L (; y )  log p (y ;  ) + log p (y ;  ) ; 1:T 1 1 t t t=2 (I B) , and, for t > 1,   (y ; ). Note the dependence of  on  and y . 1 t 1:t1 t 1:t1 There are two important caveats in the ML evaluation of GAS models. The rst one is that, from a theoretical perspective, ML estimation of GAS models is an on{going research topic. General results are reported by Harvey (2013), Blasques et al. (2014a) and Blasques et al. (2014b), while results for speci c models have been derived by Andres (2014) and Blasques et al. (2014d). The second one is that, even when the ML estimator is consistent and asymptotically normal, the numerical maximization of the loglikelihood function in (5) can be challenging, because of the nonlinearities induced by  () and the way y enters the scaled score s . Consequently, t t when the optimizer is gradient{based, good starting values need to be selected for GAS models. In the R package GAS, starting values for the optimizer are chosen in the following way: (i) estimate the static version of the model (i.e., with A = 0 and B = 0) and set the initial value of  accordingly, and (ii) perform a grid search for the coecients contained in A and B. Further technical details are presented in Section 3.2. Clearly, the coecients , A and B in (4) are di erent from those of (2), however, for notational purposes, we continue to use the same notation. 6 The R package GAS Implementation of the models in the R package GAS follows the common approach in the GAS literature. First, matrices A and B are constrained to be diagonal. Second, in order to avoid an explosive pattern for  , the spectral radius of B is constrained to be less than one. Third, the positiveness of each element of A is imposed in order to do not distort the signal coming from the conditional score s . 3. The R package GAS The R package GAS o ers an integrated environment to deal with GAS models in R. Its structure is somehow similar to the R package rugarch (Ghalanos 2015b) for GARCH models, which is widely used by practitioners and academics. The similarities concern the steps the user has to do to perform her analysis as well as the type of functions she faces. Speci cally, the rst step is to specify the model, which means choosing: (i) the assumptions for the conditional distribution of the data, (ii) the set of parameters that have to vary over time and, (iii) the scaling mechanism for the conditional score. These steps are detailed in Section 3.1. Once the model is properly speci ed, the user can estimate the unknown parameters in  by numerical maximization of the log{likelihood function as detailed in Section 3.2. Finally, predictions according to the estimated model can be easily performed; see Section 3.3. Simulation of GAS models is presented in Section 3.4. Functions for: (i) speci cation, (ii) estimation, (iii) forecasting and (iv) simulation are avail- able for univariate and multivariate time series. The general nomenclature for the functions when we consider univariate time series is \UniGAS...()" and that for multivariate time series is \MultiGAS...()" . In the R package GAS, several datasets are also included for reproducibility purposes, such as: US in ation (cpichg), US unemployment rate (usunp), realized volatility of the S&P500 Index (sp500rv) and intraday bid and ask quotes for Citygroup corporation (tqdata). These datasets are freely available online; see the R documentation for references. In this section, we use the monthly US in ation measured as the logarithmic change in the CPI available from the Federal Reserve Bank of St. Louis website https://fred.stlouisfed.org/. This dataset can be easily loaded in the R workspace using: data("cpichg", package = "GAS"). 3.1. Speci cation Speci cation of GAS models is the rst step the user needs to undertake. This is achieved by using the UniGASSpec() and MultiGASSpec() functions, in the cases of univariate and multivariate models, respectively. Both functions accept three arguments and return an object of the class uGASSpec and mGASSpec, respectively. The three arguments are: - Dist: A character indicating the label of the conditional distribution assumed for the data. Available distributions can be displayed using the function DistInfo() and are reported in Table 1. By default Dist = "norm", i.e., the Gaussian distribution. - ScalingType: A character indicating the scaling mechanism for the conditional score, i.e., the value of the parameter in (3). Possible choices are "Identity" ( = 0), "Inv" ( = 1) and "InvSqrt" ( = ). Note that, for some distributions only ScalingType = "Identity" is supported; see function DistInfo() and Table 1. By David Ardia, Kris Boudt, Leopoldo Catania 7 Label Name Type Parameters # Scaling Type norm Gaussian univariate location, scale 2 Identity, Inv, InvSqrt std Student{t (i) univariate location, scale, shape 3 Identity, Inv, InvSqrt sstd Skew{Student{t (ii) univariate location, scale, skewness, shape 4 Identity ald Asymmetric Laplace (iii) univariate location, scale, skewness 3 Identity, Inv, InvSqrt poi Poisson (iv) univariate location 1 Identity, Inv, InvSqrt gamma Gamma univariate scale, shape 2 Identity, Inv, InvSqrt exp Exponential (v) univariate location 1 Identity, Inv, InvSqrt beta Beta (vi) univariate scale, shape 2 Identity, Inv, InvSqrt mvnorm Multivariate Gaussian multivariate location, scale, correlation 9 (vii) Identity mvt Multivariate Student{t multivariate location, scale, correlation, shape 10 (vii) Identity Table 1: Statistical distributions or which the R package GAS provides the functionality to simulate, estimate and forecast the time{variation in its parameters. The fth column, #, reports the number of parameters of the distribution. Note: (i) the usual Student{t distribution (not reparametrised in terms of the variance parameter), (ii) the reparametrised Skew{Student{t such that the location and scale parameters coincide with the mean and the standard deviation of the distribution as done in the rugarch package, (iii) the ald distribution used the ,  and  reparametrization, as speci ed in Kotz et al. (2001), (iv) for the Poisson distribution location means the usual intensity parameter, (v) for the Beta distribution shape means the usual parameter and scale means the usual parameter, (vi) for the Exponential distribution location means the usual rate parameter, (vii) for N = 3. default ScalingType = "Identity", i.e., no scaling occurs. - GASPar: A named list with boolean entries containing information about which pa- rameters of the conditional distribution have to be time{varying. Generally, each uni- variate distribution is identi ed by a series of maximum four parameters. These are indicated by location, scale, skewness and shape. Note that, for some distribu- tions, these labels are not strictly related to their literal statistical meaning. Indeed, for the Exponential distribution exp, the term location indicates the usual inten- sity rate parameter; see the DistInfo() function for more details. For multivariate distributions, the set of parameters is indicated by location, scale, correlation and shape . For example, in the case of a multivariate Student{t distribution with mean vector  , scale matrix   D R D , where D is the diagonal matrix of t t t t t scales and R is the correlation matrix, and  degrees of freedom, we have that: t t location refers to  , scale refers to D , correlation refers to R and shape refers t t to  . By default, GASPar = list(location = FALSE, scale = TRUE, skewness = FALSE, shape = FALSE) for the univariate case, and GASPar = list(location = FALSE, scale = TRUE, correlation = FALSE, shape = FALSE) for the multivariate case. The function MultiGASSpec() also accepts the additional boolean argument ScalarParameters controlling for the parametrization of A and B in (4). Setting ScalarParameters = TRUE (the default value), the coecients controlling the evolution of the location, scale and corre- lation parameters are constrained to be the same across each group. Speci cally, if y 2 < follows a GAS process with conditional multivariate Gaussian distribution, the vector of time{ In the R package GAS the information matrices and the scores are always computed using their analytical formulations. In the R documentation an extra parameter, shape2, is reported. This will be used for future extensions of the package. 8 The R package GAS varying parameters is  =  ;  ;  ;  ;  ;  ;  ;  ;  . If ScalarParameters t 1;t 2;t 3;t 1;t 2;t 3;t 21;t 31;t 32;t = TRUE, the matrix of coecients A is parameterized as: A  diag a ; a ; a ; a ; a ; a ; a ; a ; a ; while, if ScalarParameters = FALSE, the matrix of coecients A takes the form: A  diag a ; a ; a ; a ; a ; a ; a ; a ; a : 1 2 3 1 2 3 21 31 32 Hence, in the latter case, each element of  evolves heterogeneously with respect to the others. The same constraints are applied to B, which means that, if ScalarParameters = TRUE, for the general N case, the number of parameters decreases from 3N (N + 1) =2 to N (N + 1) =2 + 2. Additional constraints are introduced through the GASPar argument as in the univariate case; see help("MultiGASSpec"). As an illustration, assume that we want to specify a Student{t GAS model with time{varying conditional mean and scale parameters, but xed degree of freedom, i.e.,  = . This can be easily done with the following lines of code: R> GASSpec <- UniGASSpec(Dist = "std", ScalingType = "Identity", GASPar = list(location = TRUE, scale = TRUE, shape = FALSE)) Details about the object returned from UniGASSpec() are printed in the console by simply calling GASSpec: R> GASSpec ------------------------------------------------------- - Univariate GAS Specification - ------------------------------------------------------- Conditional distribution ------------------------------------------------------- Name: Student-t Label: std Type: univariate Parameters: location, scale, shape Number of Parameters: 3 References: ------------------------------------------------------- GAS specification ------------------------------------------------------- Score scaling type: Identity Time varying parameters: location, scale ------------------------------------------------------- Since the scaling matrix S is set to the identify matrix (i.e., ScalingType = "Identity") this model for the conditional Student{t distribution corresponds to the one described in Ap- pendix A. Multivariate GAS speci cations are analogously speci ed using the MultiGASSpec() function; see help("MultiGASSpec"). David Ardia, Kris Boudt, Leopoldo Catania 9 3.2. Estimation Similar to model speci cation, estimation is handled with two di erent functions for univariate and multivariate models: UniGASFit() and MultiGASFit(), respectively. These functions require only two arguments: the GAS speci cation object GASSpec and the data, and returns an object of the class uGASFit or mGASFit. As an example, let us estimate the GAS model previously speci ed using the US in ation data included in the R package GAS: R> data("cpichg", package = "GAS") R> Fit <- UniGASFit(GASSpec, cpichg) The computational time is less than one second on a modern computer. The optimizer used is the General Nonlinear Augmented Lagrange Multiplier method of Ye (1988) available in the R package Rsolnp (Ghalanos and Theussl 2016). Results can be inspected by calling the object Fit. R> Fit ------------------------------------------ - Univariate GAS Fit - ------------------------------------------ Model Specification: T = 276 Conditional distribution: std Score scaling type: Identity Time varying parameters: location, scale ------------------------------------------ Estimates: Estimate Std. Error t value Pr(>|t|) kappa1 0.03735 0.02991 1.249 1.059e-01 kappa2 -0.25994 0.13553 -1.918 2.755e-02 kappa3 0.92671 0.76127 1.217 1.117e-01 a1 0.07173 0.01780 4.030 2.787e-05 a2 0.45377 0.20828 2.179 1.468e-02 b1 0.94318 0.02600 36.272 0.000e+00 b2 0.85559 0.07185 11.908 0.000e+00 ------------------------------------------ Unconditional Parameters: location scale shape 0.6574 0.1653 6.5262 ------------------------------------------ Information Criteria: AIC BIC np llk 370.4 395.8 7.0 -178.2 ------------------------------------------ 10 The R package GAS Elapsed time: 0.01 mins The output printed in the console is divided into: (i) the summary of the model, (ii) the estimated coecients along with signi cance levels according to their asymptotic normal distribution, (iii) the unconditional level of the parameters, i.e., ((I B) b), (iv) AIC and BIC information criteria in addition to the number of estimated parameters (np), the log{likelihood (llk) evaluated at its optimum, and (v) the computation time. Concerning the estimated coecients, kappa1, kappa2 and kappa3 are the elements of vector in (9), i.e.,  ;  and  , respectively. Analogously, a1 and a2 are the estimates of a and a and b1 and b2 are estimates of b and b , where  refers to the scale parameter of the Student{t distribution; see Appendix A. Note that, since we have speci ed scale = FALSE in the UniGASSpec() function, coecients a3 and b3, corresponding to a and b are not reported (and constrained to zero during the optimization). The R package GAS provides several methods to extract the relevant estimated quantities for objects of the class uGASFit or mGASFit. They allow us to: (i) calculate several quan- tities of the estimated conditional distribution at each point in time, such as: quantiles, conditional moments and ltered parameters (see quantile(Fit), getMoments(Fit) and getFilteredParameters(Fit), respectively), (ii) extract the estimated coecients (coef(Fit)), (iii) generate a graphical representation of the results (plot(Fit)); see help("uGASFit") for details. 3.3. Forecasting Forecasting is a crucial aspect in applied time series analysis. Given the parametric assump- tion of GAS models, predictions are usually given in the form of density forecasts, i.e., the distribution of y jy for h  1. Knowing the predictive density, practitioners can extract T +h 1:T any relevant quantities such as future expected value E [y ] or (co{)variance V [y ]. T T +h T T +h For GAS models, the one{step ahead predictive distribution (h = 1) is analytically available while it needs to be estimated by simulation in the multi{step ahead case (h > 1). The R package GAS can handle both one{step and multi{step ahead forecasts. Consis- tent with previous nomenclature, functions for univariate and multivariate predictions are UniGASFor() and MultiGASFor(), respectively. These functions accept an object of the class uGASFit or mGASFit, created using the functions UniGASFit() and MultiGASFit(), and re- turn an object of the class uGASFor and mGASFor, respectively. Additional arguments are: • H: a numeric integer value representing the forecast horizon, i.e., h. By default H = 1. • B: a numeric integer value representing the number of draws to approximate the multi{ step ahead predictive distribution when h > 1. By default B = 1e4. • ReturnDraws: a boolean controlling if the simulated draws from y jy , y jy ; : : : , T +1 1:T T +2 1:T y jy have to be returned. By default ReturnDraws = FALSE. T +h 1:T Other arguments to perform rolling{type of forecasts are detailed in the documentation; see help("UniGASFor"). Practically, if we want to predict the next{year in ation (i.e., h = 12 with the monthly series cpichg), after having estimated the GAS model of Section 3.2, we can execute the following code: David Ardia, Kris Boudt, Leopoldo Catania 11 R> Forecast <- UniGASFor(Fit, H = 12) and inspect the results by calling the object Forecast: R> Forecast ------------------------------------------ - Univariate GAS Forecast - ------------------------------------------ Model Specification Conditional distribution: std Score scaling type: Identity Horizon: 12 Rolling forecast: FALSE ------------------------------------------ Parameters forecast: location scale shape T+1 0.10128 0.1524 6.526 T+2 0.09497 0.1737 6.526 T+3 0.09380 0.2151 6.526 T+4 0.09253 0.2577 6.526 T+5 0.08745 0.3020 6.526 .................... location scale shape T+8 0.08345 0.4219 6.526 T+9 0.07790 0.4575 6.526 T+10 0.07380 0.4900 6.526 T+11 0.07557 0.5199 6.526 T+12 0.07507 0.5465 6.526 which returns some model information and the predictions of future model parameters based on averages over B draws. Forecast is an object of the class uGASFor and comes with several methods to extract and visualize the results; see help("uGASFor"). As commonly done in time series analysis, predictions are generated from models tted to rolling windows. The R package GAS includes this functionality via UniGASRoll() and MultiGASRoll(). These functions accept several arguments that we brie y describe in the univariate case: • data: a vector of length T +ForecastLength containing all the observations. • GASSpec: an object of the class uGASSpec created with UniGASSpec(). • ForecastLength: a numeric integer which speci es the length of the out{of{sample. • RefitEvery: a numeric integer of periods before coecients are re{estimated. 12 The R package GAS • RefitWindow: a character for the type of the window. As in the R package rugarch, we de ne the options: RefitWindow = "recursive" and RefitWindow = "moving". If RefitWindow = "recursive" all past observations are used when the model is re{ estimated. If RefitWindow = "moving" initial observations are eliminated. Other arguments useful to tailor the forecasting procedure and to parallelize the code execu- tion are available and detailed in the R documentation; see help("UniGASRoll"). Suppose now we are interested in assessing the forecast performance of the GAS model with a Student{t conditional distribution and time{varying location and scale parameters, detailed in Appendix A, and speci ed in the object GASSpec in Section 3.1. We treat the last 150 observations of cpichg as out{of{sample and run a rolling{window forecast exercise using the following portion of code: Roll <- UniGASRoll(cpichg, GASSpec, ForecastLength = 150, RefitEvery = 3, RefitWindow = "moving") where model coecients are re{estimated quarterly (i.e., every three observations with monthly data) using a moving windows (RefitWindow = "moving"). The code automatically makes a series of one{step ahead rolling predictions according to the model estimated using only the past information. This way, the user can perform out{of{sample analysis with GAS models. The object Roll belongs to the class uGASRoll which, as uGASFit and uGASFor, comes with several methods to extract and represent the results; see help("uGASRoll"). 3.4. Simulation Simulation of univariate and multivariate GAS models is straightforward with the R package GAS. This can be easily done via UniGASSim() and MultiGASSim(); see the R documentation. Several examples, also investigating the nite sample properties of the ML estimator for GAS models, are reported in the inst/test/Simulation.R le included in the package tarball. Besides selecting the conditional distribution of the time series process, the user of course needs to specify also the static parameters  governing the dynamics in  . More precisely, for simulation of GAS models, the vector  and the matrices A and B, need to be speci ed. It is worth stressing that, the de nition of  can be tricky, especially for multivariate mod- els. The diculty emerges from the fact that,  determines the unconditional value of the reparametrised vector of parameters  , implying that, if the user wants to specify the model in terms of a target value  she needs to know the inverse of the mapping function  (). To address this problem, the functions UniUnmapParameters() and MultiUnmapParameters(), representing  (), are available for univariate and multivariate models, respectively. This way, the user can easily specify  such that   (I B) ( ). Table 2 lists the numerical bounds imposed for the univariate distributions, such that UniUnmapParameters() cannot takes values outside those ranges. For the multivariate case, please refer to the examples reported in the inst/test/SimulateGAS.R le included in the package tarball. Here we de ne the \target value" as the unconditional level of the parameters the user has in mind. This targeting approach requires the time{varying parameter model to be stationary, as explained in, e.g., Blasques et al. (2014c). David Ardia, Kris Boudt, Leopoldo Catania 13 Label location scale skewness shape norm < < { { std < < { (2:01; 50:0) snorm < < (0:1; 2:0) { sstd < < (0:1; 2:0) (2:01; 50:0) + + ald < < < { poi < { { { + + gamma < < { { exp < { { { + + beta < < { { Table 2: Overview of the restrictions on the allowed values for the parameters of the univariate distributions, for which the R package GAS provides the functionality to simulate, estimate and forecast the time variation in the parameters. When the parameter space is < , we use the exponential link function reported in (10) with c = 0, while when the space is of the type (a; b), we use the modi ed logistic link function reported in (11); see Appendix B. Suppose we want to simulate T = 1; 000 observations from the Student{t GAS model reported in Appendix A with time{varying location and scale, but constant shape parameters. Assume our target value for the parameters is  =  ;  ;  with  = 0:1;  = 1:5 and  = 7. The matrix A and B are de ned as: A = diag (0:1; 0:4; 0:0) B = diag (0:9; 0:95; 0:0) ; such that both the conditional mean and the conditional variance evolve quite persistently over time, while the shape parameter is constant. The implementation of UniUnmapParameters() and UniGASSim() proceeds as: A <- diag(c(0.1, 0.4, 0.0)) B <- diag(c(0.9, 0.95, 0.0)) ThetaStar <- c(0.1, 1.5, 7.0) Kappa <- (diag(3) - B) %*% UniUnmapParameters(ThetaStar, "std") Sim <- UniGASSim(T = 1000, Kappa, A, B, Dist = "std", ScalingType = "Identity") where Sim is an object of the class uGASSim and comes with several method such as show, plot, and getMoments, among others; see help("uGASSim"). 4. Applications to nancial data In order to illustrate how the R package GAS can be used in practical situations, we present an empirical application with univariate and multivariate time series of nancial returns. We consider daily log{returns (in percentage points) of the Dow Jones 30 constituents available in the dji30ret dataset. This dataset includes the closing value log{returns from March 3rd, 14 The R package GAS 1987 to February 3rd, 2009 for a total of 5,521 observations per series. The dataset can be easily loaded in the workspace using: R> library("GAS") R> data("dji30ret", package = "GAS") where dji30ret is a 5521  30 data.frame containing the daily log{returns. Our analysis is a typical out{of{sample exercise, meaning that: (i) we estimate the models using an in{ sample period, (ii) we do predictions of the conditional distribution for the observations in the out{of{sample period, (iii) and that we compare the models according to their out{of{sample performance. The models we consider are univariate/multivariate GAS models estimated with the R package GAS, and univariate/multivariate GARCH models estimated using the popular R packages rugarch (Ghalanos 2015b) and rmgarch (Ghalanos 2015a), respectively. The univariate spec- i cations we consider are: (i) the Skew{Student{t GAS model with only time{varying scale parameter (i.e., Dist = "sstd") and, (ii) the GARCH(1,1) model with Skew{Student{t dis- tributed error. For both models we employ the Skew{Student{t distribution of Fern andez and Steel (1998) reparametrised such that the location and scale parameters coincide with the mean and the standard deviation of the distribution as done in the rugarch package. For the multivariate speci cations, we consider: (i) the GAS model with conditional multi- variate Student{t distribution with time{varying scales and correlations used in Creal et al. (2011) and, (ii) the Dynamic Conditional Correlation (DCC) model of Engle (2002) with a conditional multivariate Student{t distribution. For simplicity, the multivariate analysis only considers three series of the whole dataset: Caterpillar Inc. (CAT), 3M (MMM) and P zer Inc. (PFE). The code used to specify the univariate and multivariate GAS models is: R> uGASSpec <- UniGASSpec(Dist = "sstd", ScalingType = "Identity", GASPar = list(scale = TRUE)) and: R> mGASSpec <- MultiGASSpec(Dist = "mvt", ScalingType = "Identity", GASPar = list(scale = TRUE, correlation = TRUE)) respectively. The last H = 3; 000 observations (from January 27th, 1991, to the end of the sample) com- pose the out{of{sample period. During the out{of{sample period, one{step ahead density predictions are constructed by the univariate and multivariate models. Models (and there- fore coecients) are re{estimated using a moving{window every hundredth observations, as detailed in Section 3.3. One{step ahead rolling prediction are then computed as: luGASRoll <- list() David Ardia, Kris Boudt, Leopoldo Catania 15 N <- ncol(dji30ret) for(i in 1:N){ luGASRoll[[i]] <- UniGASRoll(data = dji30ret[, i], GASSpec = uGASSpec, ForecastLength = 3000, RefitEvery = 100) names(luGASRoll) <- colnames(dji30ret) and: mGASRoll <- MultiGASRoll(data = dji30ret[, c("CAT", "MMM", "PFE")], GASSpec = mGASSpec, ForecastLength = 3000, RefitEvery = 100) for the univariate and multivariate cases, respectively. Let us now compare the ability of GAS and GARCH models in predicting the one{step ahead distribution using so{called scoring rules, which compare the predicted density with the ex{post realized value of the return and deliver a score which de nes a ranking across the alternative models at each point in time (Gneiting et al. 2007). Generally, we de ne b b S  S(y ; p(y ;  )) as the score at time t + 1 for having predicted p(y ;  ) t+1 t+1 t+1 t+1 t+1 t+1 when y has been realized. We consider two widely used scoring rules: t+1 • The average weighted Continuous Ranked Probability Score (wCRPS): T +H1 wCRPS  w (z) F z;  I dz; (6) t+1 fy <zg t+1 t=T where w (z) is a weight function that emphasizes regions of interest of the predictive distribution, such as the tails or the center. Similarly to Gneiting and Ranjan (2011), we consider the cases of: (i) a weighting that gives equal emphasis to all the parts of the distribution; w (z) = 1, (ii) a weighting that focuses on the center; w (z) =  (z); a;b (iii) a weighting that focuses on the tails; w (z) = 1  (z) = (0), (iv) a weighting a;b a;b that focuses on the right tail; w (z) =  (z), and (v) a weighting that focuses on the a;b left tail w (z) = 1  (z). The functions  (z) and  (z) are the pdf and cdf a;b a;b a;b of a Gaussian distribution with mean a and standard deviation b, respectively. The label uniform represents the case where equal emphasis is given to all the parts of the distribution. • The average Negative Log Score (NLS): T +H1 NLS  log p(y ;  ): (7) t+1 t+1 t=T 16 The R package GAS Consistent with Gneiting et al. (2007), we specify the Negative Log Score such that the \direction" between the two scoring rules is the same, i.e., forecasts with lower NLS and lower wCRPS are preferred. The two aforementioned scoring rules can be easily evaluated using the BacktestDensity() function available in the R package GAS. The BacktestDensity() function accepts an object of the class uGASRoll, and returns a list with two elements: (i) the averages negative LS and wCRPS as in (7) and (6), and (ii) their values at each point in time. Additional arguments are: • lower: numeric representing the lower bound used to approximate (6) by Monte Carlo integration as detailed in Gneiting and Ranjan (2011). • upper: numeric as lower but for the upper bound. • K: numeric integer representing the number of points used to discretize the integral in (6). By default K = 1000, plus the two numeric arguments, a and b, representing a and b in the weight functions, by default a = 0.0 and b = 1.0. In our case, in order to evaluate NLS and wCRPS for the rst asset we can simply run: R> DensityBacktest <- BacktestDensity(luGASRoll[[1]], lower = -1.0, upper = 1.0) R> DensityBacktest$average LS uniform center tails tail_r tail_l -2.389e+00 1.329e-02 5.300e-03 7.310e-06 6.645e-03 6.648e-03 where lower = -1.0 and upper = 1.0. Table 3 reports the test statistics for the Diebold and Mariano (1995) (DM) test of equal performance between the series of Negative Log Scores and weighed Continuous Ranked Probability Scores for univariate GAS and GARCH models across the out{of{sample pe- riod. Negative values indicate that GAS models generate more accurate predictions of the one{step ahead conditional distribution while positive values favour GARCH. We found that, for almost all the series, GAS outperforms GARCH at very high con dence levels according to both NLS and wCRPS. Interestingly, our results suggest that GAS delivers more accurate results whatever part of the conditional distribution the wCRPS emphasizes. For the multivariate analysis we only consider NLS. In this case, the DM test statistic is 4:15, which strongly favours the GAS model against the DCC speci cation. To further The two arguments lower and upper coincide with y and y in Equation 16 of Gneiting and Ranjan l u (2011), respectively. These are two numeric objects with no default value, i.e., the user have to de ne these values according to her research design. Equals to I in Equation 16 of Gneiting and Ranjan (2011). These values can be chosen in order to target some \optimal" prediction level, or to add more exibility and focus on speci c parts the predictive distribution; see Gneiting and Ranjan (2011). Chosen lower and upper values de ne a proper range for log{returns not in percentage points as the one considered here. David Ardia, Kris Boudt, Leopoldo Catania 17 Asset NLS Uniform Center Tails Tails{r Tails{l b a a c a b AA 1:99 2:39 2:39 1:45 2:46 2:32 a b b b b AXP 2:85 2:25 2:25 0:20 2:25 2:25 c b b b b BA 1:54 2:00 2:00 0:40 2:00 1:99 a c c b BAC 3:16 1:49 1:50 0:64 1:25 1:73 a a a a a C 4:07 3:14 3:14 0:29 3:18 3:09 a a a b a a CAT 5:47 5:38 5:38 1:94 5:32 5:43 CVX 0:71 1:13 1:13 1:27 1:12 1:13 DD 1:34 0:66 0:66 0:27 0:62 0:69 b b b b b DIS 2:17 1:76 1:76 0:40 1:78 1:74 a a a b a a GE 3:81 4:22 4:22 1:73 4:23 4:20 GM 0:09 0:44 0:45 0:45 0:51 0:37 a a a a a a HD 3:67 2:89 2:89 3:32 2:89 2:88 a a a a a a HPQ 4:38 4:44 4:44 3:42 4:42 4:45 a a a a a a IBM 3:92 4:03 4:03 2:95 4:00 4:06 a b b c b b INTC 3:16 1:72 1:72 1:48 1:80 1:65 a b b b b JNJ 3:81 2:08 2:08 0:41 2:01 2:16 b b b b b JPM 1:97 1:93 1:93 0:23 1:87 1:99 AIG 0:22 0:94 0:95 0:68 0:76 1:05 a a a a a KO 3:36 3:03 3:03 0:95 3:05 3:00 b b b b b MCD 2:04 1:96 1:96 1:02 1:96 1:95 a a a b a a MMM 4:22 4:62 4:62 2:26 4:62 4:62 a a a a a a MRK 3:82 5:10 5:10 3:70 5:18 5:02 a a a b a a MSFT 3:24 2:76 2:76 2:25 2:79 2:73 a a a a a a PFE 4:86 4:98 4:98 3:21 5:02 4:94 b b b b b b PG 1:97 1:99 1:99 2:07 2:01 1:97 T 0:31 0:16 0:16 0:46 0:14 0:18 c c c c c UTX 1:50 1:52 1:52 0:80 1:57 1:46 b b b a b b VZ 2:22 2:13 2:13 2:41 2:13 2:14 b b b b b WMT 2:15 1:83 1:83 0:13 1:88 1:78 XOM 0:17 0:39 0:39 1:26 0:42 0:37 Table 3: Test statistics for the Diebold and Mariano (1995) test of equal performance be- tween the series of negative Log Scores and weighed Continuous Ranked Probability Scores for univariate GAS and GARCH models across the out{of{sample logarithmic returns of Dow Jones 30 constituents. Negative values indicate that GAS models report more accurate pre- dictions of the one{step ahead conditional distribution while positive values favour GARCH. The apexes a; b and c represent rejection of the null hypothesis of Equal Predictive Ability at the 1%; 5% and 10% con dence levels, respectively. The out{of{sample period spans from January 27th, 1991, to February 3rd, 2009 for a total of 3,000 observations. investigate this result, we report in Figure 1 the Cumulative sum of the di erences between 18 The R package GAS 1997−03 1998−06 1999−09 2000−12 2002−03 2003−06 2004−09 2005−12 2007−03 2008−06 Figure 1: Cumulative out{of{sample Log Score di erences between the multivariate Student{t GAS and the DCC(1,1) model of Engle (2002) with multivariate Student{t errors. Periods when the plot line slopes upward represent periods in which GAS outperforms GARCH, while downward{sloping segments indicate periods when the GARCH forecast is more accurate. The blue shaded area represents periods of recession in the US economy according to the \USREC" series available from the Federal Reserve Bank of St. Louis web site at https: //fred.stlouisfed.org/series/USREC. the Log Scores (CLS) of GAS and DCC de ned as: t=T +l1 GAS DCC GASjDCC b b CLS  log p y ;  log p y ;  ; t+1 t+1 t+1 t+1 T :T +l t=T GAS DCC b b where p y ;  and p y ;  are the densities predicted from GAS and DCC t+1 t+1 t+1 t+1 evaluated in y , respectively. The series of Log Scores for the multivariate GAS models is t+1 available in the output of the BacktestDensity() function, or can be extracted using the LogScore method de ned for mGASRoll objects: R> LS_MGAS <- LogScore(mGASRoll) Looking at Figure 1, periods when the plot line slopes upward represent periods in which GAS outperforms DCC, while downward{sloping segments indicate periods when the DCC forecast is more accurate. From this plot, we clearly understand the output of the DM test. Interestingly, we found that GAS starts dominating DCC after 2000. David Ardia, Kris Boudt, Leopoldo Catania 19 5. Conclusion This article introduced the R package GAS for simulating, estimating and forecasting time{ varying parameter models under the Generalized Autoregressive Score framework. It allows practitioners in many scienti c areas to perform their applied research using GAS models in a user{friendly environment. We introduced the model speci cation in a general way and illustrated the package usage. In particular, we performed an empirical application using nancial data in which we compared the performance of univariate and multivariate GAS and GARCH models. Given the exibility of GAS models and the availability of several statistical distributions in the GAS package, a number of di erent applications can be easily handled, such as: (i) the analysis of integer valued time series using the Poisson GAS model (poi), (ii) the analysis of (0,1){bounded time series using the Beta GAS model (beta), (iii) the analysis of strictly positive time series with an inverse location/scale dependence using the Gamma GAS model (gamma). Finally, if you use R or GAS, please cite the software in publications. Computational details The results in this paper were obtained using R 3.2.3 (R Core Team 2016) with the packages: GAS version 0.1.4 (Catania et al. 2016), MASS version 7.3-45 and (Venables and Ripley 2002; Ripley 2015), Rcpp version 0.12.7 (Eddelbuettel and Fran cois 2011; Eddelbuettel et al. 2016a), RcppArmadillo version 0.7.400.2.0 (Eddelbuettel and Sanderson 2014; Eddelbuettel et al. 2016b), Rsolnp version 1.16 (Ghalanos and Theussl 2016), xts version 0.9-7 (Ryan and Ulrich 2015) and quantmod version 0.4-6 (Ryan 2015). R itself and all packages used are available from CRAN at http://CRAN.R-project.org/. The package GAS is available from the CRAN repository at https://cran.r-project.org/package=GAS. The version under development is available in GitHub at https://github.com/LeopoldoCatania/GAS. Computations were performed on a Genuine Intel® quad core CPU i7{3630QM 2.40Ghz processor. Code outputs were obtained using options(digits = 4, max.print = 40, prompt = "R> ") The folder inst/doc inside the GAS package tarball contains additional technical documen- tations. A step by step guide on how to add a new statistical distribution in the GAS package is reported in the le AddNewDistribution.pdf. Acknowledgments The authors acknowledge Google for nancial support via the Google Summer of Code 2016 project "GAS"; see https://summerofcode.withgoogle.com/projects/#4717600793690112. Any remaining errors or shortcomings are the authors' responsibility. References Andres P (2014). \Maximum Likelihood Estimates for Positive Valued Dynamic Score Models: The DySco Package." Computational Statistics & Data Analysis, 76, 34{42. doi:10.1016/ j.csda.2013.11.004. 20 The R package GAS Blasques F, Koopman SJ, Lucas A (2014a). \Maximum Likelihood Estimation for Correctly Speci ed Generalized Autoregressive Score Models: Feedback E ects, Contraction Condi- tions and Asymptotic Properties." techreport TI 14-074/III, Tinbergen Institute. URL http://www.tinbergen.nl/discussionpaper/?paper=2332. Blasques F, Koopman SJ, Lucas A (2014b). \Maximum Likelihood Estimation for Generalized Autoregressive Score Models." techreport TI 2014-029/III, Tinbergen Institute. URL http: //www.tinbergen.nl/discussionpaper/?paper=2286. Blasques F, Koopman SJ, Lucas A (2014c). \Stationarity and Ergodicity of Univariate Gen- eralized Autoregressive Score Processes." Electronic Journal of Statistics, 8(1), 1088{1112. doi:doi:10.1214/14-EJS924. Blasques F, Koopman SJ, Lucas A, Schaumburg J (2014d). \Spillover Dynamics for Systemic Risk Measurement using Spatial Financial Time Series Models." techreport TI 2014-103/III, Tinbergen Institute. URL http://www.tinbergen.nl/discussionpaper/?paper=2369. Bollerslev T (1986). \Generalized Autoregressive Conditional Heteroskedasticity." Journal of Econometrics, 31(3), 307{327. doi:10.1016/0304-4076(86)90063-1. Box GE, Jenkins GM (1970). Time Series Analysis: Forecasting and Control. Holden{Day, San Francisco. Catania L, Bill e AG (2016). \Dynamic Spatial Autoregressive Models with Autoregressive and Heteroskedastic Disturbances." Working paper, URL http://arxiv.org/abs/1602.02542. Catania L, Boudt K, Ardia D (2016). GAS: Generalised Autoregressive Score Models. R pack- age version 0.1.4, URL https://cran.r-project.org/package=GAS. Creal D, Koopman SJ, Lucas A (2011). \A Dynamic Multivariate Heavy{Tailed Model for Time{Varying Volatilities and Correlations." Journal of Business & Economic Statistics, 29(4), 552{563. doi:10.1198/jbes.2011.10070. Creal D, Koopman SJ, Lucas A (2013). \Generalized Autoregressive Score Models with Applications." Journal of Applied Econometrics, 28(5), 777{795. doi:10.1002/jae.1279. Creal D, Schwaab B, Koopman SJ, Lucas A (2014). \Observation{Driven Mixed{Measurement Dynamic Factor Models with an Application to Credit Risk." The Review of Economics and Statistics, 96(5), 898{915. doi:10.1162/REST_a_00393. Diebold FX, Mariano RS (1995). \Comparing Predictive Accuracy." Journal of Business & Economic Statistics, 13(3), 253{263. doi:10.1080/07350015.1995.10524599. Eddelbuettel D, Fran cois R (2011). \Rcpp: Seamless R and C++ Integration." Journal of Statistical Software, 40(8), 1{18. doi:10.18637/jss.v040.i08. Eddelbuettel D, Fran cois R, Allaire J, Ushey K, Kou Q, Bates D, Chambers J (2016a). Rcpp: Seamless R and C++ Integration. R package version 0.12.7, URL https://cran. r-project.org/package=Rcpp. Eddelbuettel D, Fran cois R, Bates D (2016b). RcppArmadillo: Rcpp Integration for the Armadillo Templated Linear Algebra Library. R package version 0.7.400.2.0, URL https: //cran.r-project.org/package=RcppArmadillo. David Ardia, Kris Boudt, Leopoldo Catania 21 Eddelbuettel D, Sanderson C (2014). \RcppArmadillo: Accelerating R with High{ Performance C++ Linear Algebra." Computational Statistics & Data Analysis, 71, 1054{ 1063. doi:10.1016/j.csda.2013.02.005. Engle RF (2002). \Dynamic Conditional Correlation: A Simple Class of Multivariate Gen- eralized Autoregressive Conditional Heteroskedasticity Models." Journal of Business & Economic Statistics, 20(3), 339{350. doi:10.1198/073500102288618487. Fern andez C, Steel MF (1998). \On Bayesian modeling of fat tails and skewness." Journal of the American Statistical Association, 93(441), 359{371. doi:10.1080/01621459.1998. Ghalanos A (2015a). rmgarch: Multivariate GARCH Models. R package version 1.3-0, URL https://cran.r-project.org/package=rmgarch. Ghalanos A (2015b). rugarch: Univariate GARCH Models. R package version 1.3-6, URL https://cran.r-project.org/package=rugarch. Ghalanos A, Theussl S (2016). Rsolnp: General Non{Linear Optimization using Augmented Lagrange Multiplier Method. R package version 1.16, URL https://cran.r-project.org/ package=Rsolnp. Gneiting T, Balabdaoui F, Raftery AE (2007). \Probabilistic Forecasts, Calibration and Sharpness." Journal of the Royal Statistical Society: Series B (Statistical Methodology), 69(2), 243{268. doi:10.1111/j.1467-9868.2007.00587.x. Gneiting T, Ranjan R (2011). \Comparing Density Forecasts using Threshold {and Quantile{ Weighted Scoring Rules." Journal of Business & Economic Statistics, 29(3), 411{422. doi: 10.1198/jbes.2010.08110. Harvey AC (2013). Dynamic Models for Volatility and Heavy Tails: With Applications to Financial and Economic Time Series. Cambridge University Press. Harvey AC, Luati A (2014). \Filtering with Heavy Tails." Journal of the American Statistical Association, 109(507), 1112{1122. doi:10.1080/01621459.2014.887011. Harvey AC, Sucarrat G (2014). \EGARCH Models with Fat Tails, Skewness and Leverage." Computational Statistics & Data Analysis, 76, 320{338. doi:10.1016/j.csda.2013.09. Harvey AC, Thiele S (2015). \Testing Against Changing Correlation." doi:10.1016/j. jempfin.2015.09.003. Forthcoming in the Journal of Empirical Finance. Jaeckel P, Rebonato R (1999). \The Most General Methodology for Creating a Valid Correlation Matrix for Risk Management and Option Pricing Purposes." Journal of Risk, 2(2), 17{28. URL http://www.risk.net/journal-of-risk/technical-paper/2161082/ the-methodology-creating-valid-correlation-matrix-risk-management-option-pricing-purposes. Janus P, Koopman SJ, Lucas A (2014). \Long Memory Dynamics for Multivariate Dependence under Heavy Tails." Journal of Empirical Finance, 29, 187{206. doi:10.1016/j.jempfin. 2014.09.007. 22 The R package GAS Kalman RE (1960). \A New Approach to Linear Filtering and Prediction Problems." Trans- actions of the ASME{Journal of Basic Engineering, 82(Series D), 35{45. Kotz S, Kozubowski T, Podgorski K (2001). The Laplace Distribution and Generalizations: A Revisit with Applications to Communications, Economics, Engineering, and Finance. Springer{Verlag. URL http://www.springer.com/us/book/9780817641665. Lucas A, Zhang X (2016). \Score{Driven Exponentially Weighted Moving Averages and Value{at{Risk Forecasting." International Journal of Forecasting, 32(2), 293{302. doi: 10.1016/j.ijforecast.2015.09.003. Oh DH, Patton AJ (2013). \Time{Varying Systemic Risk: Evidence from a Dynamic Cop- ula Model of CDS Spreads." doi:10.1080/07350015.2016.1177535. Forthcoming in the Journal of Business & Economic Statistics. Pinheiro JC, Bates DM (1996). \Unconstrained Parametrizations for Variance{Covariance Matrices." Statistics and Computing, 6(3), 289{296. doi:10.1007/BF00140873. Pourahmadi M, Wang X (2015). \Distribution of Random Correlation Matrices: Hyperspher- ical Parameterization of the Cholesky Factor." Statistics & Probability Letters, 106, 5{12. R Core Team (2016). R: A Language and Environment for Statistical Computing. R Foundation for Statistical Computing, Vienna, Austria. R version 3.2.3, URL https: //www.R-project.org/. Rapisarda F, Brigo D, Mercurio F (2007). \Parameterizing Correlations: A Geomet- ric Interpretation." IMA Journal of Management Mathematics, 18(1), 55{73. doi: 10.1016/j.spl.2015.06.015. Ripley B (2015). MASS: Support Functions and Datasets for Venables and Ripley's MASS. R package version 7.3-45, URL https://CRAN.R-project.org/package=MASS. Ryan JA (2015). quantmod: Quantitative Financial Modelling Framework. R package version 0.4-6, URL https://CRAN.R-project.org/package=quantmod. Ryan JA, Ulrich JM (2015). xts: Extensible Time Series. R package version 0.9-7, URL https://CRAN.R-project.org/package=xts. Sanderson C (2010). \Armadillo: An Open Source C++ Linear Algebra Library for Fast Prototyping and Computationally Intensive Experiments." Technical report, NICTA. URL http://arma.sourceforge.net/. Shephard N (2005). Stochastic Volatility: Selected Readings. Oxford University Press, Oxford. Sucarrat G (2013). \betategarch: Simulation, estimation and forecasting of Beta-Skew-t- EGARCH models." The R Journal, 5(2), 137{147. URL https://journal.r-project. org/archive/2013-2/sucarrat.pdf. Venables WN, Ripley BD (2002). Modern Applied Statistics with S. Fourth edition. Springer{ Verlag, New York. ISBN 0-387-95457-0, URL http://www.stats.ox.ac.uk/pub/MASS4. Ye Y (1988). Interior Algorithms for Linear, Quadratic, and Linearly Constrained Convex Programming. Ph.D. thesis, Stanford University. David Ardia, Kris Boudt, Leopoldo Catania 23 Zhu D, Galbraith JW (2010). \A Generalized Asymmetric Student-t Distribution with Ap- plication to Financial Econometrics." Journal of Econometrics, 157(2), 297{305. doi: 10.1016/j.jeconom.2010.01.013. 24 The R package GAS A. The GAS model with conditional Student{t distribution Let us consider the case where the distribution of the univariate random variable y 2 <, conditionally on y , is Student{t with location  , scale  and  degrees of freedom, i.e., 1:t1 t t t = ( ;  ;  ) and: t t t t ! t 2 2 2 (y  ) t t p(y ;  )   1 + : (8) t t t t 2 t t As will become clear, the score corresponding to the Student{t distribution has the advantage of dampening the e ect of extreme observations on the future volatility, when the Student{t has suciently fat tails. It has been used by Creal et al. (2013) and Lucas and Zhang (2016) under the name tGAS, and by Harvey (2013) and Harvey and Luati (2014) under the name Beta{t{EGARCH. Di erentiating the logarithm of (8) with respect to  leads to the score vector r (y ;  ) = t t t t (r ;r ;r ) , with: t t t ( + 1)(y  ) t t t (y  ) t t 1 + t t t t ( + 1) (y  ) 1 t t t t 2 2 (y  ) t t 2  1 + t t t t 1  + 1 1  1 t t 2 2 2 2 2 2 2 1 (y  ) ( + 1) (y  ) t t t t t log 1 + +   ; 2   2 (y  ) t t t t 2  1 + t t t t where () is the Digamma function. Without loss of generality, let us consider the case where = 0 with no reparametrization, i.e.,  =  . The results when 6= 0 and a t t mapping function () for  is introduced are qualitatively the same. Clearly, what controls for the response to extreme observations in the conditional score r (y ;  ) is the degree of t t t freedom parameter  . When  is small, say  = 3, the conditional distribution of y has t t t t high probability mass in the tails, which means that extreme observations, which would be considered outliers under the conditionally normal distribution, are likely to be observed. If we introduce the following mapping function for the unrestricted vector of parameter  = (e ;  ; e ) : t t t >  e t t e e ( )    exp( ) t t exp(e ) + c ; t t with c = 2 in order to ensure the existence of V [y ], then the GAS updating step for t1 t t when = 0 takes the form: ( ) t+1 t+1 (9) e e e + AJ ( ) r (y ;  ) + B ; t+1 t t t t t David Ardia, Kris Boudt, Leopoldo Catania 25 where   ( ;  ;  ) , A  diag(a ; a ; a ) and B  diag(b ; b ; b ). In this particular case, the Jacobian matrix J ( ) takes the form: 0 1 1 0 0 @ e A J ( ) = : 0 exp( ) 0 t t 0 0 exp(e ) Constraints on the evolution of the GAS parameters can be easily considered by xing the values of the A and B elements. For example, if the constraint  =  has to be imposed, we set a = b = 0 during the (log-)likelihood maximization. B. Mapping functions Now we brie y discuss the choice of the mapping function () for GAS models. We indicate e e the i{th element of  and  as  and  , respectively. Analogously, we refer to the i{th t t i;t i;t element of the vector{valued mapping function () as  (), such that  ( ) =  . i i i;t i;t Generally, there are three types of constraints we want to impose on  : i;t 1)  > c; c 2 < i;t 2)  2 (a; b), for a; b 2 < and b > a i;t 3)  2 (a; b)j 2  for a; b 2 < and b > a , i;t t the additional case when  2 <, and thus  =  , implicitly requires that  : < ! < is i;t i;t i;t i the identity function. The rst case,  > c, covers the situation where, for example,  is a scale parameter and, i;t i;t consequently, its positiveness has to be imposed (i.e., c = 0). In this case,  : < ! [c;1), and the exponential link function, de ned as: = exp( ) + c ; (10) i;t i;t can be employed. The second case,  2 (a; b), covers the situation where, for example, i;t p (;  ), is the asymmetric Student{t distribution of Zhu and Galbraith (2010), and  is its t i;t skew parameter de ned in (0; 1). In the more general case we have  : < ! (a; b), and thus, the modi ed logistic function: b a = a + ; (11) i;t 1 + exp( ) can be employed. The last case,  2 (a; b)j 2 , is more complicated and covers the i;t t situation where, for example, p (;  ) is a multivariate Gaussian distribution and  is one t i;t element of its correlation matrix R . Clearly, in this case   [1; 1], with the equivalence t i;t corresponding to the case N = 2. For the more general case N > 1, we need to ensure 0 N that R is positive de nite, i.e., x R x > 0;8x 2 < . Following Creal et al. (2011), we t t employ the hyperspherical coordinates transformation originally proposed by Pinheiro and The case  < c follows immediately. i;t 26 The R package GAS Bates (1996) and subsequently discussed in Jaeckel and Rebonato (1999), Rapisarda et al. (2007) and Pourahmadi and Wang (2015). We de ne the general (h; k){th lower diagonal element of R as  =  for h > k, h < N and e =  , for i = 1; : : : ; N (N 1) =2. t hk;t i;t hk;t i;t Pourahmadi and Wang (2015) show that: h1 m1 h1 X Y Y = c c + c c s s + c s s 1  h < k  N ; hk;t h1;t k1;t hm;t km;t hl;t kl;t hk;t hl;t kl;1 m=2 l=1 l=1 where c  cos e and s  sin e for all 1  h < k  N ensure that R hk;t hk;t hk;t hk;t t f g is a proper correlation matrix. ij;t i;j=1 These three speci cations for  () cover all the cases considered in this article and in the R package GAS. Additional information are reported in the R documentation. Fore details on () and  (); see help("UniMapParameters") and help("UniUnmapParameters") in the univariate case and help("MultiMapParameters") and help("MultiUnmapParameters") in the multivariate case. Aliation: David Ardia Institute of Financial Analysis University of Neuch^ atel, Switzerland and Department of Finance, Insurance and Real Estate Laval University, Canada E-mail: david.ardia@unine.ch Kris Boudt Vrije Universiteit Brussel, Belgium and Vrije Universiteit Amsterdam, The Netherlands E-mail: kris.boudt@vub.ac.be Leopoldo Catania (corresponding author) Department of Economics and Finance Faculty of Economics University of Rome, \Tor Vergata" Via Columbia, 2 00133 Rome, Italy E-mail: leopoldo.catania@uniroma2.it http://www.deepdyve.com/assets/images/DeepDyve-Logo-lg.png Statistics arXiv (Cornell University)

Generalized Autoregressive Score Models in R: The GAS Package

Statistics , Volume 2021 (1609) – Sep 8, 2016

Loading next page...
 
/lp/arxiv-cornell-university/generalized-autoregressive-score-models-in-r-the-gas-package-RxDUeseNbR

References

References for this paper are not available at this time. We will be adding them shortly, thank you for your patience.

eISSN
ARCH-3347
DOI
10.18637/jss.v088.i06
Publisher site
See Article on Publisher Site

Abstract

This paper presents the R package GAS for the analysis of time series under the Generalized Autoregressive Score (GAS) framework of Creal et al. (2013) and Harvey (2013). The distinctive feature of the GAS approach is the use of the score function as the driver of time{variation in the parameters of nonlinear models. The GAS package provides functions to simulate univariate and multivariate GAS processes, estimate the GAS parameters and to make time series forecasts. We illustrate the use of the GAS package with a detailed case study on estimating the time{varying conditional densities of a set of nancial assets. Keywords : GAS, time series models, score models, dynamic conditional score, R software. arXiv:1609.02354v1 [stat.CO] 8 Sep 2016 2 The R package GAS 1. Introduction Time{variation in the parameters describing a stochastic time series process is pervasive in al- most all applied scienti c elds. Early references to time series models include Kalman (1960) and Box and Jenkins (1970). In many settings, the model of interest is characterized by time{ varying parameters, for which the literature has proposed a myriad of possible speci cations. Recently, Creal et al. (2013) and Harvey (2013) note that many of the proposed models are either dicult to estimate (in particular, the class of stochastic volatility models reviewed in Shephard (2005)) and/or do not properly take the shape of the conditional distribution of the data into account. Creal et al. (2013) and Harvey (2013) therefore propose to use the score of the conditional density function as the main driver of time{variation in the param- eters of the time series process used to describe the data. A further advantage of using the conditional score as driver is that the estimation by Maximum Likelihood is straightforward. The resulting model is referred to as: Score{Driven model, Dynamic Conditional Score (DCS) model, or Generalized Autoregressive Score (GAS) model. In this article and accompanying R package, we use the GAS acronym for both. The R package GAS is conceived to be of relevance for the modelling of all types of time series data. It does not matter whether they are real{valued, integer{valued, (0,1){bounded or strictly positive, as long as there is a conditional density for which the score function and the Hessian are well{de ned. The practical relevance of the GAS framework has been illustrated in the case of nancial risk forecasting (see e.g., Harvey and Sucarrat (2014) for market risk, Oh and Patton (2013) for systematic risk, and Creal et al. (2014) for credit risk analysis), dependence modelling (see e.g., Harvey and Thiele (2015) and Janus et al. (2014)), and spatial econometrics (see e.g., Blasques et al. (2014d) and Catania and Bill e (2016)). For a more complete overview of the work on GAS models, we refer the reader to the GAS community page at http://www.gasmodel.com/. It is important to note that, even though the GAS framework has been developed by econome- tricians, it is exible enough to be used in all elds in which the use of time{varying parameter models is relevant. The main diculty in using GAS models is to derive the score and Hessian and implementing the Maximum Likelihood estimation of the resulting nonlinear models. The R package GAS answers these needs by proposing an integrated set of R functions to do time series analysis in the R statistical language (R Core Team 2016) under the GAS framework. The functionalities include: (i) estimation, (ii) prediction, (iii) simulation, (iv) backtesting, and (v) graphical representation of the results, implying that it is ready to use in real{life ap- plications. The user interface uses the R programming language, which has the advantage of being free and open source. However, most of the underlying routines are principally written in C++ exploiting the armadillo library (Sanderson 2010) and the R packages Rcpp (Ed- delbuettel and Fran cois 2011; Eddelbuettel et al. 2016a) and RcppArmadillo (Eddelbuettel and Sanderson 2014; Eddelbuettel et al. 2016b) to speed up the computations. Furthermore, since the package is written with the S4 methods, R users with basic programming knowledge will nd common functions such as coef(), plot() and show() to extract and analyze their A typical example is the class of (G)ARCH models in which the squared (demeaned) return is the driver of time{variation in the conditional variance, independently of the shape of the conditional distribution of the return. To see that this is counter{intuitive, consider the case of observing a 10% return when the conditional mean is 0% and the volatility is 3%. Under the assumption of a normal distribution, the 10% return is a strong signal of an increase in volatility, while under a fat{tailed Student{t distribution, the signal is weakened because of the higher probability that the extreme value is an observation from the tails. David Ardia, Kris Boudt, Leopoldo Catania 3 results. We believe that this aspect is of primary importance since it dramatically increases the number of potential users. The R package GAS is available from the CRAN repository at https://cran.r-project.org/package=GAS. Other codes available for speci c GAS models are available in the GAS community page at http://www.gasmodel.com/. For instance, the R package betategarch (Sucarrat 2013) allows us to estimate the beta{t{EGARCH model of Harvey (2013) and its skewed version introduced by Harvey and Sucarrat (2014). The outline of the paper is as follows: Section 2 reviews the GAS framework to de ne time{ varying parameter models, referring to the seminal works of Creal et al. (2013) and Harvey (2013). Section 3 introduces the R package GAS and illustrates how to to simulate, estimate and make predictions. Section 4 presents a real{life application to nancial data. Section 5 concludes. 2. The GAS framework to modeling time{varying parameters One of the most appealing characteristics of the GAS framework is its applicability to de ne time{varying parameter models in a large variety of univariate and multivariate time series settings. We try to be as general as possible in reviewing the GAS framework, and report in Appendix A the detailed equations for the speci c case of a conditionally Student{t distributed random variable. In this section, we rst introduce the notation and present the GAS model when the parameter space is unrestricted. We then show how a mapping function can be used to model the time{variation in the parameters when the parameter space is restricted. The section concludes by summarizing the Maximum Likelihood approach for GAS model estimation. 2.1. Model speci cation Let y 2 < be an N {dimensional random vector at time t with conditional distribution: y jy  p(y ;  ) ; (1) t 1:t1 t t 0 0 0 J where y  (y ; : : : ; y ) contains the past values of y up to time t 1 and  2   < 1:t1 1 t1 t t is a vector of time{varying parameters which fully characterizes p() and only depends on y and a set of static additional parameters , i.e.,   (y ; ) for all t. The main 1:t1 t 1:t1 feature of GAS models is that the evolution in the time{varying parameter vector  is driven by the score of the conditional distribution de ned in (1), together with an autoregressive component: + A s + B  ; (2) t+1 t t where, , A and B are matrices of coecients with proper dimensions collected in , and s is a vector which is proportional to the score of (1): s  S ( )r (y ;  ) : t t t t t t The matrix S is a J  J positive de nite scaling matrix known at time t and: @ log p(y ;  ) t t r (y ;  )  ; t t t t 4 The R package GAS is the score of (1) evaluated at  . Creal et al. (2013) suggest to set the scaling matrix S to t t a power > 0 of the inverse of the Information Matrix of  to account for the variance of r . More precisely: S ( )  I ( ) ; t t t t with: I ( )  E r (y ;  )r (y ;  ) ; (3) t t t1 t t t t t t where the expectation is taken with respect to the conditional distribution of y given y . t 1:t1 The additional parameter is xed by the user and usually takes value in the set f0; ; 1g. When = 0, S = I and there is no scaling. If = 1 ( = ), the conditional score r (y ;  ) t t t t is premultiplied by the inverse of (the square root of ) its covariance matrix I ( ). t t It is worth noting that, whatever the choice of , s is a Martingale Di erence (MD) with respect to the distribution of y given y , i.e., E [s ] = 0 for all t. Furthermore, when t 1:t1 t1 t = , the additional moment condition V [s ] = I can be easily derived. Due to the t1 t fact that s is a MD, if the spectral radius of B is less then one , the updating equation of  reported in (2) implies a mean reverting process for  through the long{term mean t t 1 1 (I B) , which means that the unconditional value of  is (I B) . It follows that, the J {valued vector  and the J J matrix B control for the level and the persistence of the process, respectively. The additional J  J matrix of coecients A, that premultiplies the scaled score s , controls for the impact of s to  . Speci cally, as detailed in Creal et al. (2013), the quantity t t+1 s indicates the direction to update the vector of parameters from  , to  , acting as t t t+1 a steepest ascent algorithm for improving the model local t given the current parameter position. Interestingly, this updating procedure resembles the well{known Newton{Raphson algorithm. Hence, A can be interpreted as the step of the update, and needs to be designed in a way to do not distort the signal coming from s ; see Section 2.3. 2.2. Reparametrization In (2) the parameter vector  has a linear speci cation and is thus unbounded. In practice, the parameter space of  is often restricted (  < ). For instance, when we model the scale parameter of a Student{t distribution, we need to ensure its positiveness. Even if this problem can be solved by imposing constraints on  (as is done in the GARCH model, see, Bollerslev 1986), the standard solution under the GAS framework is to use a (possibly J J e e nonlinear) link function () that maps  2 < into  and where  2 < has the linear t t t dynamic speci cation of (2). Speci cally, let  : < !  be a twice di erentiable vector{ valued mapping function such that ( ) =  . The updating equation for  is then given t t t by: ( ) t t (4) e e + Ae s + B ; t t t1 e e e e e e e where s  S ( )r (y ;  ) and r (y ;  ) represents the score of (1) with respect to  , and, t t t t t t t t t t e e e e e consequently, S ( ) can depend on the information matrix of  given by I ( ). Denote the t t t t t We denote by I the identity matrix of appropriate size. The spectral radius of a L  L matrix X is de ned as  (X)  max (j j; : : : ;j j), where  is the i{th 1 L i eigenvalue of X, for i = 1; : : : ; L. David Ardia, Kris Boudt, Leopoldo Catania 5 Jacobian matrix of () evaluated at  as follows: @( ) J ( )  : Then, the following relations hold: e e e r (y ;  ) = J ( ) r (y ;  ) t t t t t t t e e e e I ( ) = J ( ) I ( )J ( ) : t t t t t t This way, almost all the nonlinear constraints can be easily handled via the de nition of a proper mapping function () and its associated Jacobian matrix J (). The coecients to be estimated are gathered into   (; A; B) and estimated by numerically maximizing the (log-)likelihood function as detailed in Section 2.3. In Appendix B we discuss the choice of mapping function for GAS models in more details. 2.3. Maximum likelihood estimation A useful property of GAS models is that, given the past information and the static parameter vector , the vector of time{varying parameters,  , is perfectly predictable and the loglikeli- hood function can be easily evaluated via the prediction error decomposition. More precisely, for a sample of T realizations of y , collected in y , the vector of parameters  can be t 1:T estimated by Maximum Likelihood (ML) as the solution of: arg max L (; y ) ; (5) 1:T where: L (; y )  log p (y ;  ) + log p (y ;  ) ; 1:T 1 1 t t t=2 (I B) , and, for t > 1,   (y ; ). Note the dependence of  on  and y . 1 t 1:t1 t 1:t1 There are two important caveats in the ML evaluation of GAS models. The rst one is that, from a theoretical perspective, ML estimation of GAS models is an on{going research topic. General results are reported by Harvey (2013), Blasques et al. (2014a) and Blasques et al. (2014b), while results for speci c models have been derived by Andres (2014) and Blasques et al. (2014d). The second one is that, even when the ML estimator is consistent and asymptotically normal, the numerical maximization of the loglikelihood function in (5) can be challenging, because of the nonlinearities induced by  () and the way y enters the scaled score s . Consequently, t t when the optimizer is gradient{based, good starting values need to be selected for GAS models. In the R package GAS, starting values for the optimizer are chosen in the following way: (i) estimate the static version of the model (i.e., with A = 0 and B = 0) and set the initial value of  accordingly, and (ii) perform a grid search for the coecients contained in A and B. Further technical details are presented in Section 3.2. Clearly, the coecients , A and B in (4) are di erent from those of (2), however, for notational purposes, we continue to use the same notation. 6 The R package GAS Implementation of the models in the R package GAS follows the common approach in the GAS literature. First, matrices A and B are constrained to be diagonal. Second, in order to avoid an explosive pattern for  , the spectral radius of B is constrained to be less than one. Third, the positiveness of each element of A is imposed in order to do not distort the signal coming from the conditional score s . 3. The R package GAS The R package GAS o ers an integrated environment to deal with GAS models in R. Its structure is somehow similar to the R package rugarch (Ghalanos 2015b) for GARCH models, which is widely used by practitioners and academics. The similarities concern the steps the user has to do to perform her analysis as well as the type of functions she faces. Speci cally, the rst step is to specify the model, which means choosing: (i) the assumptions for the conditional distribution of the data, (ii) the set of parameters that have to vary over time and, (iii) the scaling mechanism for the conditional score. These steps are detailed in Section 3.1. Once the model is properly speci ed, the user can estimate the unknown parameters in  by numerical maximization of the log{likelihood function as detailed in Section 3.2. Finally, predictions according to the estimated model can be easily performed; see Section 3.3. Simulation of GAS models is presented in Section 3.4. Functions for: (i) speci cation, (ii) estimation, (iii) forecasting and (iv) simulation are avail- able for univariate and multivariate time series. The general nomenclature for the functions when we consider univariate time series is \UniGAS...()" and that for multivariate time series is \MultiGAS...()" . In the R package GAS, several datasets are also included for reproducibility purposes, such as: US in ation (cpichg), US unemployment rate (usunp), realized volatility of the S&P500 Index (sp500rv) and intraday bid and ask quotes for Citygroup corporation (tqdata). These datasets are freely available online; see the R documentation for references. In this section, we use the monthly US in ation measured as the logarithmic change in the CPI available from the Federal Reserve Bank of St. Louis website https://fred.stlouisfed.org/. This dataset can be easily loaded in the R workspace using: data("cpichg", package = "GAS"). 3.1. Speci cation Speci cation of GAS models is the rst step the user needs to undertake. This is achieved by using the UniGASSpec() and MultiGASSpec() functions, in the cases of univariate and multivariate models, respectively. Both functions accept three arguments and return an object of the class uGASSpec and mGASSpec, respectively. The three arguments are: - Dist: A character indicating the label of the conditional distribution assumed for the data. Available distributions can be displayed using the function DistInfo() and are reported in Table 1. By default Dist = "norm", i.e., the Gaussian distribution. - ScalingType: A character indicating the scaling mechanism for the conditional score, i.e., the value of the parameter in (3). Possible choices are "Identity" ( = 0), "Inv" ( = 1) and "InvSqrt" ( = ). Note that, for some distributions only ScalingType = "Identity" is supported; see function DistInfo() and Table 1. By David Ardia, Kris Boudt, Leopoldo Catania 7 Label Name Type Parameters # Scaling Type norm Gaussian univariate location, scale 2 Identity, Inv, InvSqrt std Student{t (i) univariate location, scale, shape 3 Identity, Inv, InvSqrt sstd Skew{Student{t (ii) univariate location, scale, skewness, shape 4 Identity ald Asymmetric Laplace (iii) univariate location, scale, skewness 3 Identity, Inv, InvSqrt poi Poisson (iv) univariate location 1 Identity, Inv, InvSqrt gamma Gamma univariate scale, shape 2 Identity, Inv, InvSqrt exp Exponential (v) univariate location 1 Identity, Inv, InvSqrt beta Beta (vi) univariate scale, shape 2 Identity, Inv, InvSqrt mvnorm Multivariate Gaussian multivariate location, scale, correlation 9 (vii) Identity mvt Multivariate Student{t multivariate location, scale, correlation, shape 10 (vii) Identity Table 1: Statistical distributions or which the R package GAS provides the functionality to simulate, estimate and forecast the time{variation in its parameters. The fth column, #, reports the number of parameters of the distribution. Note: (i) the usual Student{t distribution (not reparametrised in terms of the variance parameter), (ii) the reparametrised Skew{Student{t such that the location and scale parameters coincide with the mean and the standard deviation of the distribution as done in the rugarch package, (iii) the ald distribution used the ,  and  reparametrization, as speci ed in Kotz et al. (2001), (iv) for the Poisson distribution location means the usual intensity parameter, (v) for the Beta distribution shape means the usual parameter and scale means the usual parameter, (vi) for the Exponential distribution location means the usual rate parameter, (vii) for N = 3. default ScalingType = "Identity", i.e., no scaling occurs. - GASPar: A named list with boolean entries containing information about which pa- rameters of the conditional distribution have to be time{varying. Generally, each uni- variate distribution is identi ed by a series of maximum four parameters. These are indicated by location, scale, skewness and shape. Note that, for some distribu- tions, these labels are not strictly related to their literal statistical meaning. Indeed, for the Exponential distribution exp, the term location indicates the usual inten- sity rate parameter; see the DistInfo() function for more details. For multivariate distributions, the set of parameters is indicated by location, scale, correlation and shape . For example, in the case of a multivariate Student{t distribution with mean vector  , scale matrix   D R D , where D is the diagonal matrix of t t t t t scales and R is the correlation matrix, and  degrees of freedom, we have that: t t location refers to  , scale refers to D , correlation refers to R and shape refers t t to  . By default, GASPar = list(location = FALSE, scale = TRUE, skewness = FALSE, shape = FALSE) for the univariate case, and GASPar = list(location = FALSE, scale = TRUE, correlation = FALSE, shape = FALSE) for the multivariate case. The function MultiGASSpec() also accepts the additional boolean argument ScalarParameters controlling for the parametrization of A and B in (4). Setting ScalarParameters = TRUE (the default value), the coecients controlling the evolution of the location, scale and corre- lation parameters are constrained to be the same across each group. Speci cally, if y 2 < follows a GAS process with conditional multivariate Gaussian distribution, the vector of time{ In the R package GAS the information matrices and the scores are always computed using their analytical formulations. In the R documentation an extra parameter, shape2, is reported. This will be used for future extensions of the package. 8 The R package GAS varying parameters is  =  ;  ;  ;  ;  ;  ;  ;  ;  . If ScalarParameters t 1;t 2;t 3;t 1;t 2;t 3;t 21;t 31;t 32;t = TRUE, the matrix of coecients A is parameterized as: A  diag a ; a ; a ; a ; a ; a ; a ; a ; a ; while, if ScalarParameters = FALSE, the matrix of coecients A takes the form: A  diag a ; a ; a ; a ; a ; a ; a ; a ; a : 1 2 3 1 2 3 21 31 32 Hence, in the latter case, each element of  evolves heterogeneously with respect to the others. The same constraints are applied to B, which means that, if ScalarParameters = TRUE, for the general N case, the number of parameters decreases from 3N (N + 1) =2 to N (N + 1) =2 + 2. Additional constraints are introduced through the GASPar argument as in the univariate case; see help("MultiGASSpec"). As an illustration, assume that we want to specify a Student{t GAS model with time{varying conditional mean and scale parameters, but xed degree of freedom, i.e.,  = . This can be easily done with the following lines of code: R> GASSpec <- UniGASSpec(Dist = "std", ScalingType = "Identity", GASPar = list(location = TRUE, scale = TRUE, shape = FALSE)) Details about the object returned from UniGASSpec() are printed in the console by simply calling GASSpec: R> GASSpec ------------------------------------------------------- - Univariate GAS Specification - ------------------------------------------------------- Conditional distribution ------------------------------------------------------- Name: Student-t Label: std Type: univariate Parameters: location, scale, shape Number of Parameters: 3 References: ------------------------------------------------------- GAS specification ------------------------------------------------------- Score scaling type: Identity Time varying parameters: location, scale ------------------------------------------------------- Since the scaling matrix S is set to the identify matrix (i.e., ScalingType = "Identity") this model for the conditional Student{t distribution corresponds to the one described in Ap- pendix A. Multivariate GAS speci cations are analogously speci ed using the MultiGASSpec() function; see help("MultiGASSpec"). David Ardia, Kris Boudt, Leopoldo Catania 9 3.2. Estimation Similar to model speci cation, estimation is handled with two di erent functions for univariate and multivariate models: UniGASFit() and MultiGASFit(), respectively. These functions require only two arguments: the GAS speci cation object GASSpec and the data, and returns an object of the class uGASFit or mGASFit. As an example, let us estimate the GAS model previously speci ed using the US in ation data included in the R package GAS: R> data("cpichg", package = "GAS") R> Fit <- UniGASFit(GASSpec, cpichg) The computational time is less than one second on a modern computer. The optimizer used is the General Nonlinear Augmented Lagrange Multiplier method of Ye (1988) available in the R package Rsolnp (Ghalanos and Theussl 2016). Results can be inspected by calling the object Fit. R> Fit ------------------------------------------ - Univariate GAS Fit - ------------------------------------------ Model Specification: T = 276 Conditional distribution: std Score scaling type: Identity Time varying parameters: location, scale ------------------------------------------ Estimates: Estimate Std. Error t value Pr(>|t|) kappa1 0.03735 0.02991 1.249 1.059e-01 kappa2 -0.25994 0.13553 -1.918 2.755e-02 kappa3 0.92671 0.76127 1.217 1.117e-01 a1 0.07173 0.01780 4.030 2.787e-05 a2 0.45377 0.20828 2.179 1.468e-02 b1 0.94318 0.02600 36.272 0.000e+00 b2 0.85559 0.07185 11.908 0.000e+00 ------------------------------------------ Unconditional Parameters: location scale shape 0.6574 0.1653 6.5262 ------------------------------------------ Information Criteria: AIC BIC np llk 370.4 395.8 7.0 -178.2 ------------------------------------------ 10 The R package GAS Elapsed time: 0.01 mins The output printed in the console is divided into: (i) the summary of the model, (ii) the estimated coecients along with signi cance levels according to their asymptotic normal distribution, (iii) the unconditional level of the parameters, i.e., ((I B) b), (iv) AIC and BIC information criteria in addition to the number of estimated parameters (np), the log{likelihood (llk) evaluated at its optimum, and (v) the computation time. Concerning the estimated coecients, kappa1, kappa2 and kappa3 are the elements of vector in (9), i.e.,  ;  and  , respectively. Analogously, a1 and a2 are the estimates of a and a and b1 and b2 are estimates of b and b , where  refers to the scale parameter of the Student{t distribution; see Appendix A. Note that, since we have speci ed scale = FALSE in the UniGASSpec() function, coecients a3 and b3, corresponding to a and b are not reported (and constrained to zero during the optimization). The R package GAS provides several methods to extract the relevant estimated quantities for objects of the class uGASFit or mGASFit. They allow us to: (i) calculate several quan- tities of the estimated conditional distribution at each point in time, such as: quantiles, conditional moments and ltered parameters (see quantile(Fit), getMoments(Fit) and getFilteredParameters(Fit), respectively), (ii) extract the estimated coecients (coef(Fit)), (iii) generate a graphical representation of the results (plot(Fit)); see help("uGASFit") for details. 3.3. Forecasting Forecasting is a crucial aspect in applied time series analysis. Given the parametric assump- tion of GAS models, predictions are usually given in the form of density forecasts, i.e., the distribution of y jy for h  1. Knowing the predictive density, practitioners can extract T +h 1:T any relevant quantities such as future expected value E [y ] or (co{)variance V [y ]. T T +h T T +h For GAS models, the one{step ahead predictive distribution (h = 1) is analytically available while it needs to be estimated by simulation in the multi{step ahead case (h > 1). The R package GAS can handle both one{step and multi{step ahead forecasts. Consis- tent with previous nomenclature, functions for univariate and multivariate predictions are UniGASFor() and MultiGASFor(), respectively. These functions accept an object of the class uGASFit or mGASFit, created using the functions UniGASFit() and MultiGASFit(), and re- turn an object of the class uGASFor and mGASFor, respectively. Additional arguments are: • H: a numeric integer value representing the forecast horizon, i.e., h. By default H = 1. • B: a numeric integer value representing the number of draws to approximate the multi{ step ahead predictive distribution when h > 1. By default B = 1e4. • ReturnDraws: a boolean controlling if the simulated draws from y jy , y jy ; : : : , T +1 1:T T +2 1:T y jy have to be returned. By default ReturnDraws = FALSE. T +h 1:T Other arguments to perform rolling{type of forecasts are detailed in the documentation; see help("UniGASFor"). Practically, if we want to predict the next{year in ation (i.e., h = 12 with the monthly series cpichg), after having estimated the GAS model of Section 3.2, we can execute the following code: David Ardia, Kris Boudt, Leopoldo Catania 11 R> Forecast <- UniGASFor(Fit, H = 12) and inspect the results by calling the object Forecast: R> Forecast ------------------------------------------ - Univariate GAS Forecast - ------------------------------------------ Model Specification Conditional distribution: std Score scaling type: Identity Horizon: 12 Rolling forecast: FALSE ------------------------------------------ Parameters forecast: location scale shape T+1 0.10128 0.1524 6.526 T+2 0.09497 0.1737 6.526 T+3 0.09380 0.2151 6.526 T+4 0.09253 0.2577 6.526 T+5 0.08745 0.3020 6.526 .................... location scale shape T+8 0.08345 0.4219 6.526 T+9 0.07790 0.4575 6.526 T+10 0.07380 0.4900 6.526 T+11 0.07557 0.5199 6.526 T+12 0.07507 0.5465 6.526 which returns some model information and the predictions of future model parameters based on averages over B draws. Forecast is an object of the class uGASFor and comes with several methods to extract and visualize the results; see help("uGASFor"). As commonly done in time series analysis, predictions are generated from models tted to rolling windows. The R package GAS includes this functionality via UniGASRoll() and MultiGASRoll(). These functions accept several arguments that we brie y describe in the univariate case: • data: a vector of length T +ForecastLength containing all the observations. • GASSpec: an object of the class uGASSpec created with UniGASSpec(). • ForecastLength: a numeric integer which speci es the length of the out{of{sample. • RefitEvery: a numeric integer of periods before coecients are re{estimated. 12 The R package GAS • RefitWindow: a character for the type of the window. As in the R package rugarch, we de ne the options: RefitWindow = "recursive" and RefitWindow = "moving". If RefitWindow = "recursive" all past observations are used when the model is re{ estimated. If RefitWindow = "moving" initial observations are eliminated. Other arguments useful to tailor the forecasting procedure and to parallelize the code execu- tion are available and detailed in the R documentation; see help("UniGASRoll"). Suppose now we are interested in assessing the forecast performance of the GAS model with a Student{t conditional distribution and time{varying location and scale parameters, detailed in Appendix A, and speci ed in the object GASSpec in Section 3.1. We treat the last 150 observations of cpichg as out{of{sample and run a rolling{window forecast exercise using the following portion of code: Roll <- UniGASRoll(cpichg, GASSpec, ForecastLength = 150, RefitEvery = 3, RefitWindow = "moving") where model coecients are re{estimated quarterly (i.e., every three observations with monthly data) using a moving windows (RefitWindow = "moving"). The code automatically makes a series of one{step ahead rolling predictions according to the model estimated using only the past information. This way, the user can perform out{of{sample analysis with GAS models. The object Roll belongs to the class uGASRoll which, as uGASFit and uGASFor, comes with several methods to extract and represent the results; see help("uGASRoll"). 3.4. Simulation Simulation of univariate and multivariate GAS models is straightforward with the R package GAS. This can be easily done via UniGASSim() and MultiGASSim(); see the R documentation. Several examples, also investigating the nite sample properties of the ML estimator for GAS models, are reported in the inst/test/Simulation.R le included in the package tarball. Besides selecting the conditional distribution of the time series process, the user of course needs to specify also the static parameters  governing the dynamics in  . More precisely, for simulation of GAS models, the vector  and the matrices A and B, need to be speci ed. It is worth stressing that, the de nition of  can be tricky, especially for multivariate mod- els. The diculty emerges from the fact that,  determines the unconditional value of the reparametrised vector of parameters  , implying that, if the user wants to specify the model in terms of a target value  she needs to know the inverse of the mapping function  (). To address this problem, the functions UniUnmapParameters() and MultiUnmapParameters(), representing  (), are available for univariate and multivariate models, respectively. This way, the user can easily specify  such that   (I B) ( ). Table 2 lists the numerical bounds imposed for the univariate distributions, such that UniUnmapParameters() cannot takes values outside those ranges. For the multivariate case, please refer to the examples reported in the inst/test/SimulateGAS.R le included in the package tarball. Here we de ne the \target value" as the unconditional level of the parameters the user has in mind. This targeting approach requires the time{varying parameter model to be stationary, as explained in, e.g., Blasques et al. (2014c). David Ardia, Kris Boudt, Leopoldo Catania 13 Label location scale skewness shape norm < < { { std < < { (2:01; 50:0) snorm < < (0:1; 2:0) { sstd < < (0:1; 2:0) (2:01; 50:0) + + ald < < < { poi < { { { + + gamma < < { { exp < { { { + + beta < < { { Table 2: Overview of the restrictions on the allowed values for the parameters of the univariate distributions, for which the R package GAS provides the functionality to simulate, estimate and forecast the time variation in the parameters. When the parameter space is < , we use the exponential link function reported in (10) with c = 0, while when the space is of the type (a; b), we use the modi ed logistic link function reported in (11); see Appendix B. Suppose we want to simulate T = 1; 000 observations from the Student{t GAS model reported in Appendix A with time{varying location and scale, but constant shape parameters. Assume our target value for the parameters is  =  ;  ;  with  = 0:1;  = 1:5 and  = 7. The matrix A and B are de ned as: A = diag (0:1; 0:4; 0:0) B = diag (0:9; 0:95; 0:0) ; such that both the conditional mean and the conditional variance evolve quite persistently over time, while the shape parameter is constant. The implementation of UniUnmapParameters() and UniGASSim() proceeds as: A <- diag(c(0.1, 0.4, 0.0)) B <- diag(c(0.9, 0.95, 0.0)) ThetaStar <- c(0.1, 1.5, 7.0) Kappa <- (diag(3) - B) %*% UniUnmapParameters(ThetaStar, "std") Sim <- UniGASSim(T = 1000, Kappa, A, B, Dist = "std", ScalingType = "Identity") where Sim is an object of the class uGASSim and comes with several method such as show, plot, and getMoments, among others; see help("uGASSim"). 4. Applications to nancial data In order to illustrate how the R package GAS can be used in practical situations, we present an empirical application with univariate and multivariate time series of nancial returns. We consider daily log{returns (in percentage points) of the Dow Jones 30 constituents available in the dji30ret dataset. This dataset includes the closing value log{returns from March 3rd, 14 The R package GAS 1987 to February 3rd, 2009 for a total of 5,521 observations per series. The dataset can be easily loaded in the workspace using: R> library("GAS") R> data("dji30ret", package = "GAS") where dji30ret is a 5521  30 data.frame containing the daily log{returns. Our analysis is a typical out{of{sample exercise, meaning that: (i) we estimate the models using an in{ sample period, (ii) we do predictions of the conditional distribution for the observations in the out{of{sample period, (iii) and that we compare the models according to their out{of{sample performance. The models we consider are univariate/multivariate GAS models estimated with the R package GAS, and univariate/multivariate GARCH models estimated using the popular R packages rugarch (Ghalanos 2015b) and rmgarch (Ghalanos 2015a), respectively. The univariate spec- i cations we consider are: (i) the Skew{Student{t GAS model with only time{varying scale parameter (i.e., Dist = "sstd") and, (ii) the GARCH(1,1) model with Skew{Student{t dis- tributed error. For both models we employ the Skew{Student{t distribution of Fern andez and Steel (1998) reparametrised such that the location and scale parameters coincide with the mean and the standard deviation of the distribution as done in the rugarch package. For the multivariate speci cations, we consider: (i) the GAS model with conditional multi- variate Student{t distribution with time{varying scales and correlations used in Creal et al. (2011) and, (ii) the Dynamic Conditional Correlation (DCC) model of Engle (2002) with a conditional multivariate Student{t distribution. For simplicity, the multivariate analysis only considers three series of the whole dataset: Caterpillar Inc. (CAT), 3M (MMM) and P zer Inc. (PFE). The code used to specify the univariate and multivariate GAS models is: R> uGASSpec <- UniGASSpec(Dist = "sstd", ScalingType = "Identity", GASPar = list(scale = TRUE)) and: R> mGASSpec <- MultiGASSpec(Dist = "mvt", ScalingType = "Identity", GASPar = list(scale = TRUE, correlation = TRUE)) respectively. The last H = 3; 000 observations (from January 27th, 1991, to the end of the sample) com- pose the out{of{sample period. During the out{of{sample period, one{step ahead density predictions are constructed by the univariate and multivariate models. Models (and there- fore coecients) are re{estimated using a moving{window every hundredth observations, as detailed in Section 3.3. One{step ahead rolling prediction are then computed as: luGASRoll <- list() David Ardia, Kris Boudt, Leopoldo Catania 15 N <- ncol(dji30ret) for(i in 1:N){ luGASRoll[[i]] <- UniGASRoll(data = dji30ret[, i], GASSpec = uGASSpec, ForecastLength = 3000, RefitEvery = 100) names(luGASRoll) <- colnames(dji30ret) and: mGASRoll <- MultiGASRoll(data = dji30ret[, c("CAT", "MMM", "PFE")], GASSpec = mGASSpec, ForecastLength = 3000, RefitEvery = 100) for the univariate and multivariate cases, respectively. Let us now compare the ability of GAS and GARCH models in predicting the one{step ahead distribution using so{called scoring rules, which compare the predicted density with the ex{post realized value of the return and deliver a score which de nes a ranking across the alternative models at each point in time (Gneiting et al. 2007). Generally, we de ne b b S  S(y ; p(y ;  )) as the score at time t + 1 for having predicted p(y ;  ) t+1 t+1 t+1 t+1 t+1 t+1 when y has been realized. We consider two widely used scoring rules: t+1 • The average weighted Continuous Ranked Probability Score (wCRPS): T +H1 wCRPS  w (z) F z;  I dz; (6) t+1 fy <zg t+1 t=T where w (z) is a weight function that emphasizes regions of interest of the predictive distribution, such as the tails or the center. Similarly to Gneiting and Ranjan (2011), we consider the cases of: (i) a weighting that gives equal emphasis to all the parts of the distribution; w (z) = 1, (ii) a weighting that focuses on the center; w (z) =  (z); a;b (iii) a weighting that focuses on the tails; w (z) = 1  (z) = (0), (iv) a weighting a;b a;b that focuses on the right tail; w (z) =  (z), and (v) a weighting that focuses on the a;b left tail w (z) = 1  (z). The functions  (z) and  (z) are the pdf and cdf a;b a;b a;b of a Gaussian distribution with mean a and standard deviation b, respectively. The label uniform represents the case where equal emphasis is given to all the parts of the distribution. • The average Negative Log Score (NLS): T +H1 NLS  log p(y ;  ): (7) t+1 t+1 t=T 16 The R package GAS Consistent with Gneiting et al. (2007), we specify the Negative Log Score such that the \direction" between the two scoring rules is the same, i.e., forecasts with lower NLS and lower wCRPS are preferred. The two aforementioned scoring rules can be easily evaluated using the BacktestDensity() function available in the R package GAS. The BacktestDensity() function accepts an object of the class uGASRoll, and returns a list with two elements: (i) the averages negative LS and wCRPS as in (7) and (6), and (ii) their values at each point in time. Additional arguments are: • lower: numeric representing the lower bound used to approximate (6) by Monte Carlo integration as detailed in Gneiting and Ranjan (2011). • upper: numeric as lower but for the upper bound. • K: numeric integer representing the number of points used to discretize the integral in (6). By default K = 1000, plus the two numeric arguments, a and b, representing a and b in the weight functions, by default a = 0.0 and b = 1.0. In our case, in order to evaluate NLS and wCRPS for the rst asset we can simply run: R> DensityBacktest <- BacktestDensity(luGASRoll[[1]], lower = -1.0, upper = 1.0) R> DensityBacktest$average LS uniform center tails tail_r tail_l -2.389e+00 1.329e-02 5.300e-03 7.310e-06 6.645e-03 6.648e-03 where lower = -1.0 and upper = 1.0. Table 3 reports the test statistics for the Diebold and Mariano (1995) (DM) test of equal performance between the series of Negative Log Scores and weighed Continuous Ranked Probability Scores for univariate GAS and GARCH models across the out{of{sample pe- riod. Negative values indicate that GAS models generate more accurate predictions of the one{step ahead conditional distribution while positive values favour GARCH. We found that, for almost all the series, GAS outperforms GARCH at very high con dence levels according to both NLS and wCRPS. Interestingly, our results suggest that GAS delivers more accurate results whatever part of the conditional distribution the wCRPS emphasizes. For the multivariate analysis we only consider NLS. In this case, the DM test statistic is 4:15, which strongly favours the GAS model against the DCC speci cation. To further The two arguments lower and upper coincide with y and y in Equation 16 of Gneiting and Ranjan l u (2011), respectively. These are two numeric objects with no default value, i.e., the user have to de ne these values according to her research design. Equals to I in Equation 16 of Gneiting and Ranjan (2011). These values can be chosen in order to target some \optimal" prediction level, or to add more exibility and focus on speci c parts the predictive distribution; see Gneiting and Ranjan (2011). Chosen lower and upper values de ne a proper range for log{returns not in percentage points as the one considered here. David Ardia, Kris Boudt, Leopoldo Catania 17 Asset NLS Uniform Center Tails Tails{r Tails{l b a a c a b AA 1:99 2:39 2:39 1:45 2:46 2:32 a b b b b AXP 2:85 2:25 2:25 0:20 2:25 2:25 c b b b b BA 1:54 2:00 2:00 0:40 2:00 1:99 a c c b BAC 3:16 1:49 1:50 0:64 1:25 1:73 a a a a a C 4:07 3:14 3:14 0:29 3:18 3:09 a a a b a a CAT 5:47 5:38 5:38 1:94 5:32 5:43 CVX 0:71 1:13 1:13 1:27 1:12 1:13 DD 1:34 0:66 0:66 0:27 0:62 0:69 b b b b b DIS 2:17 1:76 1:76 0:40 1:78 1:74 a a a b a a GE 3:81 4:22 4:22 1:73 4:23 4:20 GM 0:09 0:44 0:45 0:45 0:51 0:37 a a a a a a HD 3:67 2:89 2:89 3:32 2:89 2:88 a a a a a a HPQ 4:38 4:44 4:44 3:42 4:42 4:45 a a a a a a IBM 3:92 4:03 4:03 2:95 4:00 4:06 a b b c b b INTC 3:16 1:72 1:72 1:48 1:80 1:65 a b b b b JNJ 3:81 2:08 2:08 0:41 2:01 2:16 b b b b b JPM 1:97 1:93 1:93 0:23 1:87 1:99 AIG 0:22 0:94 0:95 0:68 0:76 1:05 a a a a a KO 3:36 3:03 3:03 0:95 3:05 3:00 b b b b b MCD 2:04 1:96 1:96 1:02 1:96 1:95 a a a b a a MMM 4:22 4:62 4:62 2:26 4:62 4:62 a a a a a a MRK 3:82 5:10 5:10 3:70 5:18 5:02 a a a b a a MSFT 3:24 2:76 2:76 2:25 2:79 2:73 a a a a a a PFE 4:86 4:98 4:98 3:21 5:02 4:94 b b b b b b PG 1:97 1:99 1:99 2:07 2:01 1:97 T 0:31 0:16 0:16 0:46 0:14 0:18 c c c c c UTX 1:50 1:52 1:52 0:80 1:57 1:46 b b b a b b VZ 2:22 2:13 2:13 2:41 2:13 2:14 b b b b b WMT 2:15 1:83 1:83 0:13 1:88 1:78 XOM 0:17 0:39 0:39 1:26 0:42 0:37 Table 3: Test statistics for the Diebold and Mariano (1995) test of equal performance be- tween the series of negative Log Scores and weighed Continuous Ranked Probability Scores for univariate GAS and GARCH models across the out{of{sample logarithmic returns of Dow Jones 30 constituents. Negative values indicate that GAS models report more accurate pre- dictions of the one{step ahead conditional distribution while positive values favour GARCH. The apexes a; b and c represent rejection of the null hypothesis of Equal Predictive Ability at the 1%; 5% and 10% con dence levels, respectively. The out{of{sample period spans from January 27th, 1991, to February 3rd, 2009 for a total of 3,000 observations. investigate this result, we report in Figure 1 the Cumulative sum of the di erences between 18 The R package GAS 1997−03 1998−06 1999−09 2000−12 2002−03 2003−06 2004−09 2005−12 2007−03 2008−06 Figure 1: Cumulative out{of{sample Log Score di erences between the multivariate Student{t GAS and the DCC(1,1) model of Engle (2002) with multivariate Student{t errors. Periods when the plot line slopes upward represent periods in which GAS outperforms GARCH, while downward{sloping segments indicate periods when the GARCH forecast is more accurate. The blue shaded area represents periods of recession in the US economy according to the \USREC" series available from the Federal Reserve Bank of St. Louis web site at https: //fred.stlouisfed.org/series/USREC. the Log Scores (CLS) of GAS and DCC de ned as: t=T +l1 GAS DCC GASjDCC b b CLS  log p y ;  log p y ;  ; t+1 t+1 t+1 t+1 T :T +l t=T GAS DCC b b where p y ;  and p y ;  are the densities predicted from GAS and DCC t+1 t+1 t+1 t+1 evaluated in y , respectively. The series of Log Scores for the multivariate GAS models is t+1 available in the output of the BacktestDensity() function, or can be extracted using the LogScore method de ned for mGASRoll objects: R> LS_MGAS <- LogScore(mGASRoll) Looking at Figure 1, periods when the plot line slopes upward represent periods in which GAS outperforms DCC, while downward{sloping segments indicate periods when the DCC forecast is more accurate. From this plot, we clearly understand the output of the DM test. Interestingly, we found that GAS starts dominating DCC after 2000. David Ardia, Kris Boudt, Leopoldo Catania 19 5. Conclusion This article introduced the R package GAS for simulating, estimating and forecasting time{ varying parameter models under the Generalized Autoregressive Score framework. It allows practitioners in many scienti c areas to perform their applied research using GAS models in a user{friendly environment. We introduced the model speci cation in a general way and illustrated the package usage. In particular, we performed an empirical application using nancial data in which we compared the performance of univariate and multivariate GAS and GARCH models. Given the exibility of GAS models and the availability of several statistical distributions in the GAS package, a number of di erent applications can be easily handled, such as: (i) the analysis of integer valued time series using the Poisson GAS model (poi), (ii) the analysis of (0,1){bounded time series using the Beta GAS model (beta), (iii) the analysis of strictly positive time series with an inverse location/scale dependence using the Gamma GAS model (gamma). Finally, if you use R or GAS, please cite the software in publications. Computational details The results in this paper were obtained using R 3.2.3 (R Core Team 2016) with the packages: GAS version 0.1.4 (Catania et al. 2016), MASS version 7.3-45 and (Venables and Ripley 2002; Ripley 2015), Rcpp version 0.12.7 (Eddelbuettel and Fran cois 2011; Eddelbuettel et al. 2016a), RcppArmadillo version 0.7.400.2.0 (Eddelbuettel and Sanderson 2014; Eddelbuettel et al. 2016b), Rsolnp version 1.16 (Ghalanos and Theussl 2016), xts version 0.9-7 (Ryan and Ulrich 2015) and quantmod version 0.4-6 (Ryan 2015). R itself and all packages used are available from CRAN at http://CRAN.R-project.org/. The package GAS is available from the CRAN repository at https://cran.r-project.org/package=GAS. The version under development is available in GitHub at https://github.com/LeopoldoCatania/GAS. Computations were performed on a Genuine Intel® quad core CPU i7{3630QM 2.40Ghz processor. Code outputs were obtained using options(digits = 4, max.print = 40, prompt = "R> ") The folder inst/doc inside the GAS package tarball contains additional technical documen- tations. A step by step guide on how to add a new statistical distribution in the GAS package is reported in the le AddNewDistribution.pdf. Acknowledgments The authors acknowledge Google for nancial support via the Google Summer of Code 2016 project "GAS"; see https://summerofcode.withgoogle.com/projects/#4717600793690112. Any remaining errors or shortcomings are the authors' responsibility. References Andres P (2014). \Maximum Likelihood Estimates for Positive Valued Dynamic Score Models: The DySco Package." Computational Statistics & Data Analysis, 76, 34{42. doi:10.1016/ j.csda.2013.11.004. 20 The R package GAS Blasques F, Koopman SJ, Lucas A (2014a). \Maximum Likelihood Estimation for Correctly Speci ed Generalized Autoregressive Score Models: Feedback E ects, Contraction Condi- tions and Asymptotic Properties." techreport TI 14-074/III, Tinbergen Institute. URL http://www.tinbergen.nl/discussionpaper/?paper=2332. Blasques F, Koopman SJ, Lucas A (2014b). \Maximum Likelihood Estimation for Generalized Autoregressive Score Models." techreport TI 2014-029/III, Tinbergen Institute. URL http: //www.tinbergen.nl/discussionpaper/?paper=2286. Blasques F, Koopman SJ, Lucas A (2014c). \Stationarity and Ergodicity of Univariate Gen- eralized Autoregressive Score Processes." Electronic Journal of Statistics, 8(1), 1088{1112. doi:doi:10.1214/14-EJS924. Blasques F, Koopman SJ, Lucas A, Schaumburg J (2014d). \Spillover Dynamics for Systemic Risk Measurement using Spatial Financial Time Series Models." techreport TI 2014-103/III, Tinbergen Institute. URL http://www.tinbergen.nl/discussionpaper/?paper=2369. Bollerslev T (1986). \Generalized Autoregressive Conditional Heteroskedasticity." Journal of Econometrics, 31(3), 307{327. doi:10.1016/0304-4076(86)90063-1. Box GE, Jenkins GM (1970). Time Series Analysis: Forecasting and Control. Holden{Day, San Francisco. Catania L, Bill e AG (2016). \Dynamic Spatial Autoregressive Models with Autoregressive and Heteroskedastic Disturbances." Working paper, URL http://arxiv.org/abs/1602.02542. Catania L, Boudt K, Ardia D (2016). GAS: Generalised Autoregressive Score Models. R pack- age version 0.1.4, URL https://cran.r-project.org/package=GAS. Creal D, Koopman SJ, Lucas A (2011). \A Dynamic Multivariate Heavy{Tailed Model for Time{Varying Volatilities and Correlations." Journal of Business & Economic Statistics, 29(4), 552{563. doi:10.1198/jbes.2011.10070. Creal D, Koopman SJ, Lucas A (2013). \Generalized Autoregressive Score Models with Applications." Journal of Applied Econometrics, 28(5), 777{795. doi:10.1002/jae.1279. Creal D, Schwaab B, Koopman SJ, Lucas A (2014). \Observation{Driven Mixed{Measurement Dynamic Factor Models with an Application to Credit Risk." The Review of Economics and Statistics, 96(5), 898{915. doi:10.1162/REST_a_00393. Diebold FX, Mariano RS (1995). \Comparing Predictive Accuracy." Journal of Business & Economic Statistics, 13(3), 253{263. doi:10.1080/07350015.1995.10524599. Eddelbuettel D, Fran cois R (2011). \Rcpp: Seamless R and C++ Integration." Journal of Statistical Software, 40(8), 1{18. doi:10.18637/jss.v040.i08. Eddelbuettel D, Fran cois R, Allaire J, Ushey K, Kou Q, Bates D, Chambers J (2016a). Rcpp: Seamless R and C++ Integration. R package version 0.12.7, URL https://cran. r-project.org/package=Rcpp. Eddelbuettel D, Fran cois R, Bates D (2016b). RcppArmadillo: Rcpp Integration for the Armadillo Templated Linear Algebra Library. R package version 0.7.400.2.0, URL https: //cran.r-project.org/package=RcppArmadillo. David Ardia, Kris Boudt, Leopoldo Catania 21 Eddelbuettel D, Sanderson C (2014). \RcppArmadillo: Accelerating R with High{ Performance C++ Linear Algebra." Computational Statistics & Data Analysis, 71, 1054{ 1063. doi:10.1016/j.csda.2013.02.005. Engle RF (2002). \Dynamic Conditional Correlation: A Simple Class of Multivariate Gen- eralized Autoregressive Conditional Heteroskedasticity Models." Journal of Business & Economic Statistics, 20(3), 339{350. doi:10.1198/073500102288618487. Fern andez C, Steel MF (1998). \On Bayesian modeling of fat tails and skewness." Journal of the American Statistical Association, 93(441), 359{371. doi:10.1080/01621459.1998. Ghalanos A (2015a). rmgarch: Multivariate GARCH Models. R package version 1.3-0, URL https://cran.r-project.org/package=rmgarch. Ghalanos A (2015b). rugarch: Univariate GARCH Models. R package version 1.3-6, URL https://cran.r-project.org/package=rugarch. Ghalanos A, Theussl S (2016). Rsolnp: General Non{Linear Optimization using Augmented Lagrange Multiplier Method. R package version 1.16, URL https://cran.r-project.org/ package=Rsolnp. Gneiting T, Balabdaoui F, Raftery AE (2007). \Probabilistic Forecasts, Calibration and Sharpness." Journal of the Royal Statistical Society: Series B (Statistical Methodology), 69(2), 243{268. doi:10.1111/j.1467-9868.2007.00587.x. Gneiting T, Ranjan R (2011). \Comparing Density Forecasts using Threshold {and Quantile{ Weighted Scoring Rules." Journal of Business & Economic Statistics, 29(3), 411{422. doi: 10.1198/jbes.2010.08110. Harvey AC (2013). Dynamic Models for Volatility and Heavy Tails: With Applications to Financial and Economic Time Series. Cambridge University Press. Harvey AC, Luati A (2014). \Filtering with Heavy Tails." Journal of the American Statistical Association, 109(507), 1112{1122. doi:10.1080/01621459.2014.887011. Harvey AC, Sucarrat G (2014). \EGARCH Models with Fat Tails, Skewness and Leverage." Computational Statistics & Data Analysis, 76, 320{338. doi:10.1016/j.csda.2013.09. Harvey AC, Thiele S (2015). \Testing Against Changing Correlation." doi:10.1016/j. jempfin.2015.09.003. Forthcoming in the Journal of Empirical Finance. Jaeckel P, Rebonato R (1999). \The Most General Methodology for Creating a Valid Correlation Matrix for Risk Management and Option Pricing Purposes." Journal of Risk, 2(2), 17{28. URL http://www.risk.net/journal-of-risk/technical-paper/2161082/ the-methodology-creating-valid-correlation-matrix-risk-management-option-pricing-purposes. Janus P, Koopman SJ, Lucas A (2014). \Long Memory Dynamics for Multivariate Dependence under Heavy Tails." Journal of Empirical Finance, 29, 187{206. doi:10.1016/j.jempfin. 2014.09.007. 22 The R package GAS Kalman RE (1960). \A New Approach to Linear Filtering and Prediction Problems." Trans- actions of the ASME{Journal of Basic Engineering, 82(Series D), 35{45. Kotz S, Kozubowski T, Podgorski K (2001). The Laplace Distribution and Generalizations: A Revisit with Applications to Communications, Economics, Engineering, and Finance. Springer{Verlag. URL http://www.springer.com/us/book/9780817641665. Lucas A, Zhang X (2016). \Score{Driven Exponentially Weighted Moving Averages and Value{at{Risk Forecasting." International Journal of Forecasting, 32(2), 293{302. doi: 10.1016/j.ijforecast.2015.09.003. Oh DH, Patton AJ (2013). \Time{Varying Systemic Risk: Evidence from a Dynamic Cop- ula Model of CDS Spreads." doi:10.1080/07350015.2016.1177535. Forthcoming in the Journal of Business & Economic Statistics. Pinheiro JC, Bates DM (1996). \Unconstrained Parametrizations for Variance{Covariance Matrices." Statistics and Computing, 6(3), 289{296. doi:10.1007/BF00140873. Pourahmadi M, Wang X (2015). \Distribution of Random Correlation Matrices: Hyperspher- ical Parameterization of the Cholesky Factor." Statistics & Probability Letters, 106, 5{12. R Core Team (2016). R: A Language and Environment for Statistical Computing. R Foundation for Statistical Computing, Vienna, Austria. R version 3.2.3, URL https: //www.R-project.org/. Rapisarda F, Brigo D, Mercurio F (2007). \Parameterizing Correlations: A Geomet- ric Interpretation." IMA Journal of Management Mathematics, 18(1), 55{73. doi: 10.1016/j.spl.2015.06.015. Ripley B (2015). MASS: Support Functions and Datasets for Venables and Ripley's MASS. R package version 7.3-45, URL https://CRAN.R-project.org/package=MASS. Ryan JA (2015). quantmod: Quantitative Financial Modelling Framework. R package version 0.4-6, URL https://CRAN.R-project.org/package=quantmod. Ryan JA, Ulrich JM (2015). xts: Extensible Time Series. R package version 0.9-7, URL https://CRAN.R-project.org/package=xts. Sanderson C (2010). \Armadillo: An Open Source C++ Linear Algebra Library for Fast Prototyping and Computationally Intensive Experiments." Technical report, NICTA. URL http://arma.sourceforge.net/. Shephard N (2005). Stochastic Volatility: Selected Readings. Oxford University Press, Oxford. Sucarrat G (2013). \betategarch: Simulation, estimation and forecasting of Beta-Skew-t- EGARCH models." The R Journal, 5(2), 137{147. URL https://journal.r-project. org/archive/2013-2/sucarrat.pdf. Venables WN, Ripley BD (2002). Modern Applied Statistics with S. Fourth edition. Springer{ Verlag, New York. ISBN 0-387-95457-0, URL http://www.stats.ox.ac.uk/pub/MASS4. Ye Y (1988). Interior Algorithms for Linear, Quadratic, and Linearly Constrained Convex Programming. Ph.D. thesis, Stanford University. David Ardia, Kris Boudt, Leopoldo Catania 23 Zhu D, Galbraith JW (2010). \A Generalized Asymmetric Student-t Distribution with Ap- plication to Financial Econometrics." Journal of Econometrics, 157(2), 297{305. doi: 10.1016/j.jeconom.2010.01.013. 24 The R package GAS A. The GAS model with conditional Student{t distribution Let us consider the case where the distribution of the univariate random variable y 2 <, conditionally on y , is Student{t with location  , scale  and  degrees of freedom, i.e., 1:t1 t t t = ( ;  ;  ) and: t t t t ! t 2 2 2 (y  ) t t p(y ;  )   1 + : (8) t t t t 2 t t As will become clear, the score corresponding to the Student{t distribution has the advantage of dampening the e ect of extreme observations on the future volatility, when the Student{t has suciently fat tails. It has been used by Creal et al. (2013) and Lucas and Zhang (2016) under the name tGAS, and by Harvey (2013) and Harvey and Luati (2014) under the name Beta{t{EGARCH. Di erentiating the logarithm of (8) with respect to  leads to the score vector r (y ;  ) = t t t t (r ;r ;r ) , with: t t t ( + 1)(y  ) t t t (y  ) t t 1 + t t t t ( + 1) (y  ) 1 t t t t 2 2 (y  ) t t 2  1 + t t t t 1  + 1 1  1 t t 2 2 2 2 2 2 2 1 (y  ) ( + 1) (y  ) t t t t t log 1 + +   ; 2   2 (y  ) t t t t 2  1 + t t t t where () is the Digamma function. Without loss of generality, let us consider the case where = 0 with no reparametrization, i.e.,  =  . The results when 6= 0 and a t t mapping function () for  is introduced are qualitatively the same. Clearly, what controls for the response to extreme observations in the conditional score r (y ;  ) is the degree of t t t freedom parameter  . When  is small, say  = 3, the conditional distribution of y has t t t t high probability mass in the tails, which means that extreme observations, which would be considered outliers under the conditionally normal distribution, are likely to be observed. If we introduce the following mapping function for the unrestricted vector of parameter  = (e ;  ; e ) : t t t >  e t t e e ( )    exp( ) t t exp(e ) + c ; t t with c = 2 in order to ensure the existence of V [y ], then the GAS updating step for t1 t t when = 0 takes the form: ( ) t+1 t+1 (9) e e e + AJ ( ) r (y ;  ) + B ; t+1 t t t t t David Ardia, Kris Boudt, Leopoldo Catania 25 where   ( ;  ;  ) , A  diag(a ; a ; a ) and B  diag(b ; b ; b ). In this particular case, the Jacobian matrix J ( ) takes the form: 0 1 1 0 0 @ e A J ( ) = : 0 exp( ) 0 t t 0 0 exp(e ) Constraints on the evolution of the GAS parameters can be easily considered by xing the values of the A and B elements. For example, if the constraint  =  has to be imposed, we set a = b = 0 during the (log-)likelihood maximization. B. Mapping functions Now we brie y discuss the choice of the mapping function () for GAS models. We indicate e e the i{th element of  and  as  and  , respectively. Analogously, we refer to the i{th t t i;t i;t element of the vector{valued mapping function () as  (), such that  ( ) =  . i i i;t i;t Generally, there are three types of constraints we want to impose on  : i;t 1)  > c; c 2 < i;t 2)  2 (a; b), for a; b 2 < and b > a i;t 3)  2 (a; b)j 2  for a; b 2 < and b > a , i;t t the additional case when  2 <, and thus  =  , implicitly requires that  : < ! < is i;t i;t i;t i the identity function. The rst case,  > c, covers the situation where, for example,  is a scale parameter and, i;t i;t consequently, its positiveness has to be imposed (i.e., c = 0). In this case,  : < ! [c;1), and the exponential link function, de ned as: = exp( ) + c ; (10) i;t i;t can be employed. The second case,  2 (a; b), covers the situation where, for example, i;t p (;  ), is the asymmetric Student{t distribution of Zhu and Galbraith (2010), and  is its t i;t skew parameter de ned in (0; 1). In the more general case we have  : < ! (a; b), and thus, the modi ed logistic function: b a = a + ; (11) i;t 1 + exp( ) can be employed. The last case,  2 (a; b)j 2 , is more complicated and covers the i;t t situation where, for example, p (;  ) is a multivariate Gaussian distribution and  is one t i;t element of its correlation matrix R . Clearly, in this case   [1; 1], with the equivalence t i;t corresponding to the case N = 2. For the more general case N > 1, we need to ensure 0 N that R is positive de nite, i.e., x R x > 0;8x 2 < . Following Creal et al. (2011), we t t employ the hyperspherical coordinates transformation originally proposed by Pinheiro and The case  < c follows immediately. i;t 26 The R package GAS Bates (1996) and subsequently discussed in Jaeckel and Rebonato (1999), Rapisarda et al. (2007) and Pourahmadi and Wang (2015). We de ne the general (h; k){th lower diagonal element of R as  =  for h > k, h < N and e =  , for i = 1; : : : ; N (N 1) =2. t hk;t i;t hk;t i;t Pourahmadi and Wang (2015) show that: h1 m1 h1 X Y Y = c c + c c s s + c s s 1  h < k  N ; hk;t h1;t k1;t hm;t km;t hl;t kl;t hk;t hl;t kl;1 m=2 l=1 l=1 where c  cos e and s  sin e for all 1  h < k  N ensure that R hk;t hk;t hk;t hk;t t f g is a proper correlation matrix. ij;t i;j=1 These three speci cations for  () cover all the cases considered in this article and in the R package GAS. Additional information are reported in the R documentation. Fore details on () and  (); see help("UniMapParameters") and help("UniUnmapParameters") in the univariate case and help("MultiMapParameters") and help("MultiUnmapParameters") in the multivariate case. Aliation: David Ardia Institute of Financial Analysis University of Neuch^ atel, Switzerland and Department of Finance, Insurance and Real Estate Laval University, Canada E-mail: david.ardia@unine.ch Kris Boudt Vrije Universiteit Brussel, Belgium and Vrije Universiteit Amsterdam, The Netherlands E-mail: kris.boudt@vub.ac.be Leopoldo Catania (corresponding author) Department of Economics and Finance Faculty of Economics University of Rome, \Tor Vergata" Via Columbia, 2 00133 Rome, Italy E-mail: leopoldo.catania@uniroma2.it

Journal

StatisticsarXiv (Cornell University)

Published: Sep 8, 2016

There are no references for this article.