https://github.com/cran/GAS
Tip revision: 4553011d1e8e2444581a5d50a8936c9b05603eb4 authored by Leopoldo Catania on 23 August 2016, 14:57:49 UTC
version 0.1.2
version 0.1.2
Tip revision: 4553011
GAS_vignette_0.1.1.Rnw
%\VignetteIndexEntry{Generalized Autoregressive Score Models in R: The GAS Package}
%\VignetteKeywords{GAS}
%\VignettePackage{GAS}
\documentclass[nojss,shortnames]{jss}
\usepackage{amsmath}
\usepackage{amsfonts}
\usepackage{graphicx}
\usepackage{caption}
\usepackage{subcaption}
\usepackage{booktabs}
\usepackage{xspace}
\usepackage{subdepth}
\usepackage{float}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% DECLARATIONS
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% DEF LEO
\def\btheta{\mbox{\boldmath $\theta$}}
\def\bxi{\mbox{\boldmath $\xi$}}
\def\bkappa{\mbox{\boldmath $\kappa$}}
\def\bmu{\mbox{\boldmath $\mu$}}
\def\bSigma{\mathbf{\Sigma}}
\def\by{\mathbf{y}}
\def\bs{\mathbf{s}}
\def\bA{\mathbf{A}}
\def\bB{\mathbf{B}}
\def\bI{\mathbf{I}}
\def\bX{\mathbf{X}}
\def\bD{\mathbf{D}}
\def\bR{\mathbf{R}}
\def\bx{\mathbf{x}}
\def\diag{\mbox{diag}}
\def\qmo{``}
\def\qmc{''}
\def\qmcsp{'' }
\newcommand\bbone{\ensuremath{\mathbbm{1}}}
\newcommand{\suchthat}{\;\ifnum\currentgrouptype=16 \middle\fi|\;}
\newcommand\independent{\protect\mathpalette{\protect\independenT}{\perp}}
\def\independenT#1#2{\mathrel{\rlap{$#1#2$}\mkern2mu{#1#2}}}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% DEF DAVE
\newcommand{\dave}[1]{\textbf{\textcolor{red}{#1}}}
\newcommand{\kris}[1]{\textbf{\textcolor{blue}{#1}}}
\newcommand{\leo}[1]{\textbf{\textcolor{green}{#1}}}
\providecommand\Student{\ensuremath{\text{Student--}t}\xspace}
\providecommand\SStudent{\ensuremath{\text{Skew--Student--}t}\xspace}
\providecommand{\V}{\mathsf{V}}
\providecommand{\bydef}{\ensuremath{\equiv}\xspace}
\newcommand{\insertfloat}[1]{%
\begin{center}
[Insert #1 about here]%
\end{center}%
}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% TITLE PAGE
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\author{%
David Ardia\\University of Neuch\^atel\\Laval University \And
Kris Boudt\\ Vrije Universiteit Brussel\\Vrije Universiteit Amsterdam\\\And
Leopoldo Catania\\University of Rome\\ \qmo Tor Vergata\qmcsp}
\Plainauthor{David Ardia, Kris Boudt, Leopoldo Catania} %% comma-separated
\title{Generalized Autoregressive Score Models in \proglang{R}:\\ The \pkg{GAS} Package}
\Plaintitle{Generalized Autoregressive Score Models in R: The GAS Package} %% without formatting
\Shorttitle{The \proglang{R} package \pkg{GAS}} %% a short title (if necessary)
\Abstract{%
This paper presents the \proglang{R} package \pkg{GAS} for the analysis of time series under the Generalized Autoregressive Score (GAS) framework of
\citet{harvey.2013} and \citet{creal_etal.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 \pkg{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 \pkg{GAS} package with a detailed case study on estimating the time--varying risk of a financial portfolio.
}
\Keywords{GAS, time series models, score models, dynamic conditional score, \proglang{R} software}
\Plainkeywords{GAS, time series models, score models, dynamic conditional score, R software}
\Address{
David Ardia\\
Institute of Financial Analysis\\
University of Neuch\^atel, Switzerland\\
and\\
Department of Finance, Insurance and Real Estate\\
Laval University, Canada\\
E-mail: \email{david.ardia@unine.ch}\\
\\
Kris Boudt\\
Vrije Universiteit Brussel, Belgium\\
and\\
Vrije Universiteit Amsterdam, The Netherlands\\
E-mail: \email{kris.boudt@vub.ac.be}\\
\\
Leopoldo Catania (corresponding author)\\
Department of Economics and Finance\\
Faculty of Economics\\
University of Rome, \qmo Tor Vergata\qmcsp\\
Via Columbia, 2\\
00133 Rome, Italy\\
E-mail: \email{leopoldo.catania@uniroma2.it}\\
%
}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% MAIN
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{document}
\newpage
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section[introduction]{Introduction}\label{sec:intro}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Time--variation in the parameters of a model describing time series process is pervasive in almost all applied scientific fields. Early references to time series models include \citet{kalman.1960} and \citet{box_jenkins.1970}. In many settings, the model of interest is characterized by time--varying parameters, for which the literature on time--varying parameter models has proposed a myriad of possible specifications. Recently, \citet{harvey.2013} and \citet{creal_etal.2013} note that many of the proposed models are either difficult to estimate (in particular, the class of stochastic volatility models reviewed in \cite{shephard.2005}) and/or do not properly take the shape of the conditional distribution of the data into account.\footnote{A typical example is the class of (G)ARCH models in which the squared (de-meaned) 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 distribution, the signal is weakened because of the possibility that the extreme value is an observation from the tails.} \citet{harvey.2013} and \citet{creal_etal.2013} therefore propose to use the score of the conditional density function as the main driver of time--variation in the parameters 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 generally referred to as the Dynamic Conditional Score (DCS) model or the Generalized Autoregressive Score (GAS) model, . In this article and accompanying \proglang{R} package, we use the GAS acronym for both.
The \proglang{R} package \pkg{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-defined. The practical relevance of the GAS framework has been illustrated in the case of financial risk forecasting (see e.g. \citet{harvey_sucarrat.2014} for market risk, \citet{oh_patton.2013} for systematic risk, and \citet{creal_etal.2014} for credit risk analysis), dependence modelling (see \emph{e.g.}, \citet{harvey_thiele.2014} and \citet{janus_etal.2014}), and spatial econometrics (see \emph{e.g.}, \citet{blasques_etal.2014a} and \citet{catania_bille.2016}). For a more complete overview of the work on GAS models, we refer the reader to the GAS community page at \href{http://www.gasmodel.com/}{http://www.gasmodel.com/}.
It is important to note that, even though the GAS framework has been developed by econometricians, it is flexible enough to be used in all fields in which the use of time--varying parameter models is relevant. The main difficulty in using GAS models is to derive the score and Hessian, and implementing the maximum likelihood estimation of the resulting nonlinear models, as well as the generation of time series forecasts. The \proglang{R} package \pkg{GAS} answers these needs by proposing an integrated set of \proglang{R} functions to do time series analysis in the \proglang{R} statistical language \citep{R} under the GAS framework. The functionality includes: (i) estimation, (ii) prediction, (iii) simulation, (iv) backtesting, and (v) graphical representation of the results, and is ready to use in real--life applications. The user interface uses the \proglang{R} programming language, which has the advantage of being free and open source. However, most of the underlying routines are principally written in \proglang{C++} exploiting the \pkg{armadillo} library \citep{sanderson.2010} and the \proglang{R} packages \pkg{Rcpp} \citep{Rcpp.2011,Rcpp} and \pkg{RcppArmadillo} \citep{RcppArmadillo.2014,RcppArmadillo} to speed up the computations. Furthermore, since the package is written with the \proglang{S4} methods, \proglang{R} users with basic programming knowledge will find common functions such as \code{coef()}, \code{plot()} and \code{show()} to extract and analyze their results. We believe that this aspect is of primary importance since it dramatically increases the number of potential users. The \proglang{R} package \pkg{GAS} is available from the CRAN repository at \url{https://cran.r-project.org/package=GAS}.
% paper schedule
The outline of the paper is as follows: Section~\ref{sec:models} reviews the GAS framework to define time--varying parameter models, referring to the seminal works of \citet{harvey.2013} and \citet{creal_etal.2013}. Section~\ref{sec:package} introduces the \proglang{R} package \pkg{GAS} and illustrates how to to simulate, estimate and make predictions. Section~\ref{sec:application} presents a real--life application to financial data. Section~\ref{sec:conclusion} concludes.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section[The GAS framework to modeling time--varying parameters]{The GAS framework to modeling time--varying parameters}\label{sec:models}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
One of the most appealing characteristics of the GAS framework is its applicability to define 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 the detailed equations for a specific model in Appendix~\ref{sec:gas_t}. First, we 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 to estimating GAS models.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsection[Model specification]{Model specification}\label{sec:notation}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Let $\by_t\in\Re^N$ be an $N$--dimensional random vector at
time $t$ with conditional distribution:
%
\begin{equation}~\label{eq:dist}
\by_t\vert\by_{1:t-1} \sim p(\by_t;\btheta_t) \,,
\end{equation}
%
where $\by_{1:t-1}\bydef(\by_1',\dots,\by_{t-1}')'$ contains the past values of $\by_t$ up to time $t-1$ and $\btheta_t\in\Theta\subseteq\Re^J$ is a vector of time--varying parameters which fully characterizes $p(\cdot)$ and only depends on $\by_{1:t-1}$ and a set of static additional parameters $\bxi$, \emph{i.e.}, $\btheta_t \bydef \btheta(\by_{1:t-1},\bxi)$ for all $t$. The main feature of GAS models is that the evolution in the time--varying parameter vector $\btheta_t$ is driven by the score of the conditional distribution defined in~\eqref{eq:dist}, together with an autoregressive component:
%
\begin{equation}~\label{eq:updating}
\btheta_{t+1} \bydef \bkappa + \bA \, \bs_t + \bB\,\btheta_t \,,
\end{equation}
%
where, $\bkappa$, $\bA$ and $\bB$ are matrices of coefficients with proper dimensions collected in $\bxi$, and $\bs_t$ is a vector which is proportional to the score of~\eqref{eq:dist}:
%
\begin{equation}~\label{eq:scaled_score}
\bs_t \bydef S_t(\btheta_t) \, \nabla_t (\by_t,\btheta_t) \,,
\end{equation}
%
where $S_t$ is a $J\times J$ positive definite scaling matrix known at time $t$ and:
%
\begin{equation}~\label{eq:score}
\nabla_t(\by_t,\btheta_t) \bydef \frac{\partial\log p(\by_t;\btheta_t)}{\partial\btheta_t} \,,
\end{equation}
%
is the score of~\eqref{eq:dist} evaluated at $\btheta_t$. \citet{creal_etal.2013} suggest to set the scaling matrix to a power $\gamma>0$ of the inverse of the Information Matrix of $\btheta_t$ to account for the variance of $\nabla_t$. More precisely:
%
\begin{equation*}
S_t \bydef \mathcal{I}_t(\btheta_t)^{-\gamma} \,,
\end{equation*}
%
with:
%
\begin{equation}~\label{eq:information_matrix}
\mathcal{I}_t(\btheta_t) \bydef \E_{t-1}\left[\nabla_t(\by_t,\btheta_t)\nabla_t(\by_t,\btheta_t)'\right] \,,
\end{equation}
%
where the expectation is taken with respect to the conditional distribution of $\by_t$ given $\by_{1:t-1}$. The additional parameter $\gamma$ is fixed by the user and usually takes value in the set $\{0,\tfrac{1}{2},1\}$. When $\gamma=0$, $S_t=\bI$ and there is no scaling. If $\gamma=1$ ($\gamma=\tfrac{1}{2}$), the conditional score $\nabla_t(\by_t,\btheta_t)$ is premultiplied by the inverse of (the square root of) its covariance matrix $\mathcal{I}_t(\btheta_t)$. It is worth noting that, whatever the choice of $\gamma$, $\bs_t$ is a Martingale Difference (MD) with respect to the distribution of $\by_t$ given $\by_{1:t-1}$, \emph{i.e.}, $\E_{t-1}\left[\bs_t\right]=\boldsymbol{0}$ for all $t$. Furthermore, when $\gamma = \tfrac{1}{2}$, the additional moment condition $\V_{t-1}\left[\bs_t\right]=\bI$ can be easily derived. Due to the fact that $\bs_t$ is a MD, if the spectral radius of $\bB$ is less then one\footnote{The spectral radius of a $L\times L$ matrix $\bX$ is defined as $\tau\left(\bX\right)\bydef\max\left(\vert\tau_1\vert, \dots, \vert\tau_L\vert\right)$, where $\tau_i$ is the $i$--th eigenvalue of $\bX$, for $i=1,\dots,L$.}, the updating equation of $\btheta_t$ reported in~\eqref{eq:updating} implies a mean reverting process for $\btheta_t$ through the long--term mean $\left(\bI - \bB\right)^{-1}\bkappa$, which means that the unconditional value of $\btheta_t$ is $\left(\bI - \bB\right)^{-1}\bkappa$. It follows that, the $J$--valued vector $\bkappa$ and the $J\times J$ matrix $\bB$, control for the level and the persistence of the process, respectively. The additional $J\times J$ matrix of coefficients $\bA$, that premultiplies the scaled score $\bs_t$, controls for the impact of $\bs_t$ to $\btheta_{t+1}$. Specifically, as detailed in \cite{creal_etal.2013}, the quantity $\bs_t$ indicates the direction to update the vector of parameters from $\btheta_t$, to $\btheta_{t+1}$, acting as a steepest ascent algorithm for improving the model local fit given the current parameter position. Interestingly, this updating procedure resembles the well--known Newton--Raphson algorithm. Hence, $\bA$ can be interpreted as the step of the update, and needs to be designed in a way to not distort the signal coming from $\bs_t$; see Section~\ref{sec:ml_gas}.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsection[Reparametrization]{Reparametrization}\label{sec:reparametrisation}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
In~\eqref{eq:updating} the parameter vector $\btheta_t$ has a linear specification and is thus unbounded. In practice, the parameter space of $\btheta_t$ is often restricted ($\Theta\subset\Re^J$). For instance, when we model the scale parameter of a \Student distribution, we need to ensure its positiveness. Even if this problem can be solved by imposing constraints on $\bxi$ (as is done in the GARCH model, see \cite{bollerslev.1987}), the standard solution under the GAS framework is to use a (possibly nonlinear) link function $\Lambda(\cdot)$ that maps $\widetilde\btheta_t\in\Re^J$ into $\btheta_t$ and where $\widetilde\btheta_t\in\Re^J$ has the linear dynamic specification of~\eqref{eq:updating}.
Specifically, let $\Lambda:\Re^J\to\Theta$ be a twice differentiable vector--valued mapping function such that $\Lambda(\widetilde\btheta_t) = \btheta_t$. The updating equation for $\btheta_t$ is then given by:
%
\begin{align}\label{eq:updating_rep}
\btheta_t &= \Lambda(\widetilde\btheta_t)\\
\widetilde\btheta_t &\bydef \bkappa + \bA\widetilde\bs_t + \bB\widetilde\btheta_{t-1} \,,
\end{align}
%
where $\widetilde\bs_t \bydef \widetilde S_t(\widetilde\btheta_t) \,\widetilde\nabla_t(\by_t,\widetilde\btheta_t)$ and $\widetilde\nabla_t(\by_t,\widetilde\btheta_t)$ represents the score of~\eqref{eq:dist} with respect to $\widetilde\btheta_t$, and, consequently, $\widetilde S_t(\widetilde\btheta_t)$ can depend on the information matrix of $\widetilde\btheta_t$ given by $\mathcal{\widetilde I}_t(\widetilde\btheta_t)$. Nicely, the following relations hold:
%
\begin{align}~\label{eq:relations}
\widetilde\nabla_t(\by_t,\widetilde\btheta_t) &= \mathcal{J}(\widetilde\btheta_t)'\nabla_t(\by_t,\btheta_t) \\
\mathcal{\widetilde I}_t(\widetilde\btheta_t) &= \mathcal{J}(\widetilde\btheta_t)'\mathcal{I}_t(\btheta_t)\mathcal{J}(\widetilde\btheta_t) \,,
\end{align}
%
where:
%
\begin{equation}~\label{eq:jacobian}
\mathcal{J}(\widetilde\btheta_t) \bydef \left.\frac{\partial\Lambda(\widetilde\btheta_t)}{\partial\widetilde\btheta_t}\right.\,,
\end{equation}
%
is the Jacobian matrix of $\Lambda(\cdot)$ evaluated at $\widetilde\btheta_t$.\footnote{The new scaled score $\widetilde\bs_t$ is still a MD with $\V_{t-1}\left[\widetilde\bs_t\right] = \mathcal{\widetilde I}_t(\widetilde\btheta_t)$.} This way, almost all the nonlinear constraints can be easily handled via the definition of a proper mapping function $\Lambda(\cdot)$ and its associated Jacobian matrix $\mathcal{J}(\cdot)$. The coefficients to be estimated are gathered into $\bxi \bydef (\bkappa, \bA, \bB)$ and estimated by numerically maximizing the (log-)likelihood function as detailed in Section~\ref{sec:ml_gas}.\footnote{Clearly, the coefficients $\bkappa$, $\bA$ and $\bB$ in~\eqref{eq:updating_rep} are different from those of~\eqref{eq:updating}, however, for notational purposes, we continue to use the same notation.} In Appendix~\ref{sec:mapping} we discuss the choice of mapping function for GAS models in more details.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsection[Maximum likelihood estimation]{Maximum likelihood estimation}\label{sec:ml_gas}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
A useful property of GAS models is that, given the past information and the static parameter vector $\bxi$, the time--varying parameters are perfectly predictable and the loglikelihood function can be easily evaluated via the prediction error decomposition. More precisely, for a sample of $T$ realizations of $\by_t$, collected in $\by_{1:T}$, the vector of parameters $\bxi$ can be estimated by Maximum Likelihood (ML) as the solution of:
\begin{equation}\label{eq:llk_max}
\widehat{\bxi} \bydef \underset{\bxi}{\arg\max}~\mathcal{L}\left(\bxi,\by_{1:T}\right),
\end{equation}
where:
\begin{equation}\label{eq:llk}
\mathcal{L}\left(\bxi,\by_{1:T}\right) \bydef \log p\left(\by_{1};\btheta_1\right) + \sum_{t=2}^{T}\log p\left(\by_{t};\btheta_t(\bxi, \by_{1:t-1})\right)\,,
\end{equation}
and $\btheta_1 \bydef (\bI - \bB)^{-1}\bkappa$. Note the dependence of $\btheta_t$ from $\bxi$ and $\by_{1:t-1}$.
There are two important caveats in the ML evaluation of GAS models. The first one is that, from a theoretical perspective, ML estimation of GAS models is an on--going research topic. General results are reported by \citet{blasques_etal.2014c}, \cite{blasques_etal.2014b} and \citet{harvey.2013}, while results for specific models have been derived by \citet{blasques_etal.2014a} and \citet{andres.2014}.
The second one is that, even when the ML estimator is consistent and asymptotically normal, the numerical maximization~\eqref{eq:llk_max} can be challenging, because of the nonlinearities induced by $\Lambda\left(\cdot\right)$ and by the way $\by_t$ enters the scaled score $\bs_t$. Consequently, when the optimizer is gradient--based, good starting values need to be selected for GAS models. In the \proglang{R} package \pkg{GAS}, starting values for the optimizer are chosen in the following way: (i) estimate the static version of the model (\emph{i.e.}, with $\bA=\boldsymbol{0}$ and $\bB=\boldsymbol{0}$) and set the initial value of $\bkappa$ accordingly, and (ii) perform a grid search for the coefficients contained in $\bA$ and $\bB$. Further technical details are presented in Section~\ref{sec:gas_est}.
Implementation of the models in the \proglang{R} package \pkg{GAS} follows the common approach in the GAS literature. First, matrices $\bB$ and $\bA$ are constrained to be diagonal. Second, in order to avoid an explosive pattern for $\widetilde\btheta_t$, the spectral radius of $\bB$ is constrained to be less then one. Third, the positiveness of each element of $\bA$ is imposed in order to do not distort the signal coming from the conditional score $\bs_t$.% (which provides the direction for parameters update as in the Newton--Raphson algorithm).
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section[The R package GAS]{The \proglang{R} package \pkg{GAS}}\label{sec:package}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
The \proglang{R} package \pkg{GAS} offers an integrated environment to deal with GAS models in \proglang{R}. Its structure is somehow similar to the \proglang{R} package \pkg{rugarch} \citep{rugarch} 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. Specifically, the first 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~\ref{sec:gas_spec}. Once the model is properly specified, the user can estimate the unknown parameters in $\bxi$ by numerical maximization of the log--likelihood function as detailed in Section~\ref{sec:gas_est}. Finally, predictions according to the estimated model can be easily performed; see Section~\ref{sec:gas_fore}. Simulation of GAS models is presented in Section~\ref{sec:gas_sim}.
Functions for: (i) specification, (ii) estimation, (iii) forecasting and (iv) simulation are available for univariate and multivariate time series. The general nomenclature for the functions when we consider univariate time series is \qmo\code{UniGAS\dots()}\qmcsp and that for multivariate time series is \qmo\code{MultiGAS\dots()}\qmcsp.
In the \proglang{R} package \pkg{GAS}, several datasets are also included for reproducibility
purposes, such as: US inflation (\code{cpichg}), US unemployment rate (\code{usunp}), Realized volatility of the S\&P500 Index (\code{sp500rv}), intraday bid and ask quotes for Citygroup corporation (\code{tqdata}); these datasets are freely available online; see the \proglang{R} documentation for references. In this Section, we use the monthly US inflation measured as the logarithmic change in the CPI available from the Federal Reserve Bank of St. Louis website \href{https://fred.stlouisfed.org/}{https://fred.stlouisfed.org/}. This dataset can be easily loaded in the \proglang{R} workspace using \code{data("cpichg")}.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsection[Specification]{Specification}\label{sec:gas_spec}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Specification of GAS models is the first step the user needs to undertake. This is achieved by using the \code{UniGASSpec()} and \code{MultiGASSpec()} functions, in the cases of univariate and multivariate models, respectively. Both functions accept three arguments and return an object of the class \code{uGASSpec} and \code{mGASSpec}, respectively. The three arguments are:
\begin{itemize}
\item[-] \code{Dist}: A \code{character} indicating the label of the conditional distributions assumed for the data. Available distributions can be displayed using the function \code{DistInfo()} and are reported in Table~\ref{tab:dist}. By default \code{Dist = "norm"}.
%
\item[-] \code{ScalingType}: A \code{character} indicating the scaling mechanism for the conditional score, \emph{i.e.}, the value of the $\gamma$ parameter in~\eqref{eq:information_matrix}. Possible choices are \code{"Identity"} ($\gamma = 0$), \code{"Inv"} ($\gamma = 1$) and \code{"InvSqrt"} ($\gamma = \tfrac{1}{2}$). Note that, for some distributions only \code{ScalingType = "Identity"} is supported; see function \code{DistInfo()} and Table~\ref{tab:dist}. By default \code{ScalingType = "Identity"}, \emph{i.e.}, no scaling occurs.\footnote{In the \pkg{GAS} package the information matrices and the scores are always computed using their analytical formulations.}
%
\item[-] \code{GASPar}: A named \code{list} with \code{boolean} entries containing information about which parameters of the conditional distribution have to be time--varying. Generally, each univariate distribution is identified by a series of maximum four parameters. These are indicated by \code{location, scale, skewness} and \code{shape}. Note that, for same distributions, these labels are not strictly related to their literal statistical meaning. Indeed, for the Exponential distribution \code{exp}, the term \code{location} indicates the usual intensity rate parameter; see the \code{DistInfo()} function for more details. For multivariate distributions, the set of parameters is indicated by \code{location, scale, correlation} and \code{shape}. For example, in the case of a multivariate \Student distribution with mean vector $\bmu_t$, scale matrix $\bSigma_t \bydef \bD_t\bR_t\bD_t$, where $\bD_t$ is the diagonal matrix of scales and $\bR_t$ is the correlation matrix, and $\nu_t$ degrees of freedom, we have that: \code{location} refers to $\bmu_t$, \code{scale} refers to $\bD_t$, \code{correlation} refers to $\bR_t$ and \code{shape} refers to $\nu_t$. By default, \code{GASPar = list(location = FALSE, scale = TRUE, skewness = FALSE, shape = FALSE)}.
\end{itemize}
\begin{table}[!t]
\centering
\resizebox{1.0\columnwidth}{!}{%
\begin{tabular}{ccccccc}
\toprule
Label & Name & Type & Parameters & \# & Scaling Type \\
\hline
\code{norm} & Gaussian & univariate & \code{location, scale} & 2 & \code{Identity, Inv, InvSqrt} \\
\code{std} & \Student (i) & univariate & \code{location, scale, shape} & 3 & \code{Identity, Inv, InvSqrt} \\
%\code{ast} & Asymmetric Student-t with two tail decay parameters & univariate & \code{location, scale, skewness, shape, shape2} & 5 & \code{Identity, Inv, InvSqrt} \\
%\code{ast1} & Asymmetric Student-t with two one decay parameter & univariate & \code{location, scale, skewness, shape} & 4 & \code{Identity, Inv, InvSqrt} \\
\code{sstd} & \SStudent (ii) & univariate & \code{location, scale, skewness, shape} & 4 & \code{Identity} \\
\code{ald} & Asymmetric Laplace (iii) & univariate & \code{location, scale, skewness} & 3 & \code{Identity, Inv, InvSqrt} \\
\code{poi} & Poisson (iv)& univariate & \code{location} & 1 & \code{Identity, Inv, InvSqrt} \\
\code{ber} & Bernoulli & univariate & \code{location} & 1 & \code{Identity, Inv, InvSqrt} \\
\code{gamma} & Gamma & univariate & \code{scale, shape} & 2 & \code{Identity, Inv, InvSqrt} \\
\code{exp} & Exponential (v)& univariate & \code{location} & 1 & \code{Identity, Inv, InvSqrt} \\
\code{beta} & Beta (vi)& univariate & \code{scale, shape} & 2 & \code{Identity, Inv, InvSqrt} \\
\code{mvnorm} & Multivariate Gaussian & multivariate & \code{location, scale, correlation} & 9 (vii) & \code{Identity} \\
\code{mvt} & Multivariate Student-t & multivariate & \code{location, scale, correlation, shape} & 10 (vii) & \code{Identity} \\
\bottomrule
\end{tabular}
}
\caption{Statistical distributions included in the \proglang{R} package \pkg{GAS}. The fifth column, \#, reports the number of parameters of the distribution. Note: (i) The usual \Student distribution (not reparametrised in terms of the variance parameter),
%(ii) \code{ast1} is the Constraint version of ast see \citet{zhu_galbraith.2010},
(ii) reparametrised such that the location and scale parameters coincide with the mean and the standard deviation of the distribution as done in the \pkg{rugarch} package. (iii) for the \code{ald} distribution the $\theta$, $\sigma$, $\kappa$ reparametrization is used; see \citet{kotz_etal.2001}, (iv) For the Poisson distribution \code{location} means the usual intensity parameter, (v) For the Beta distribution \code{shape} means the usual alpha parameter and \code{scale} means the usual beta parameter, (vi) For the Exponential distribution 'location' means the usual rate parameter. (vii) For $N=3$.}
\label{tab:dist}
\end{table}
The function \code{MultiGASSpec()} also accepts the additional boolean argument \code{ScalarParameters} controlling for the parametrization of $\bA$ and $\bB$ in~\eqref{eq:updating_rep}. Setting \code{ScalarParameters = TRUE} (the default value), the coefficients controlling the evolution of the location, scale and correlation parameters are constrained to be the same across each group. Specifically, if $\by_t\in\Re^3$ follows a GAS process with conditional multivariate Gaussian distribution, the set of parameters is $\btheta_t = \left(\mu_{1,t}, \mu_{2,t}, \mu_{3,t}, \sigma_{1,t}, \sigma_{2,t}, \sigma_{3,t}, \rho_{21,t}, \rho_{31,t},\rho_{32,t}\right)^\prime$, and the matrix of coefficients $\bA$ is parameterized as:
\begin{equation}\label{eq:scalarparameters}
\bA \bydef \diag\left(a_\mu, a_\mu, a_\mu, a_\sigma, a_\sigma, a_\sigma, a_\rho, a_\rho, a_\rho\right),
\end{equation}
while, if \code{ScalarParameters = FALSE}, the matrix of coefficients $\bA$ takes the form:
\begin{equation}\label{eq:scalarparameters}
\bA \bydef \diag\left(a_{\mu_1}, a_{\mu_2}, a_{\mu_3}, a_{\sigma_1}, a_{\sigma_2}, a_{\sigma_3}, a_{\rho_{21}}, a_{\rho_{31}}, a_{\rho_{32}}\right).
\end{equation}
Hence, in the latter case, each element of $\btheta_t$ evolves heterogeneously with respect to the others. The same constraints are applied to $\bB$, which means that, if \code{ScalarParameters = TRUE}, for the general $N$ case, the number of parameters decreases from $3N\left(N+1\right)/2$ to $N\left(N+1\right)/2 + 2$. Additional constraints are introduced through the \code{GASPar} argument as in the univariate case; see \code{help("MultiGASSpec")}.
As an illustration, assume that we want to specify a \Student GAS model with time--varying conditional mean and scale parameters, but fixed degree of freedom, \emph{i.e.}, $\nu_t = \nu$. This can be easily done with the following lines of code:
%
\begin{CodeChunk}
\begin{CodeInput}
R> GASSpec <- UniGASSpec(Dist = "std", ScalingType = "Identity",
GASPar = list(location = TRUE, scale = TRUE,
shape = FALSE))
\end{CodeInput}
\end{CodeChunk}
%
Details about the object returned from \code{UniGASSpec()} are printed in the console by simply calling \code{GASSpec}:
%
\begin{CodeChunk}
\begin{CodeInput}
R> GASSpec
\end{CodeInput}
\begin{CodeOutput}
-------------------------------------------------------
- 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
-------------------------------------------------------
\end{CodeOutput}
\end{CodeChunk}
%
As we have not specified a scaling (\code{ScalingType = "Identity"}) this model is the same as the one reported in Appendix~\ref{sec:gas_t}. Multivariate GAS specifications are analogously specified using the \code{MultiGASSpec()} function; see \code{help(MultiGASSpec)}.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsection[Estimation]{Estimation}\label{sec:gas_est}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Similar to model specification, estimation is handled with two different functions for univariate and multivariate models, \code{UniGASFit()} and \code{MultiGASFit()}, respectively. These functions require only two arguments: The GAS specification object \code{GASSpec} and the data, and returns an object of the class \code{uGASFit} or \code{mGASFit}. As an example, let us estimate the GAS model previously specified using the US inflation data included in the \proglang{R} package \pkg{GAS}:
%
\begin{CodeChunk}
\begin{CodeInput}
R> data("cpichg", package = "GAS")
R> Fit <- UniGASFit(GASSpec, cpichg)
\end{CodeInput}
\end{CodeChunk}
%
The computational time is less than one second on a modern computer. The optimizer used is the General Nonlinear Augmented Lagrange Multiplier method of \cite{ye.1988} available in the \proglang{R} package \pkg{Rsolnp} \citep{Rsolnp}. Results can be inspected by calling the object \code{Fit}.
%
\begin{CodeChunk}
\begin{CodeInput}
R> Fit
\end{CodeInput}
\begin{CodeOutput}
------------------------------------------
- Univariate GAS Fit -
------------------------------------------
Model Specification :
T = 276
Conditional distribution : std
Score scaling type : Identity
Time varying parameres : 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
\end{CodeOutput}
\end{CodeChunk}
%
The output printed in the console is divided into: (i) the summary of the model, (ii) the estimated coefficients along with significance levels according to their asymptotic normal distribution, (iii) the unconditional level of the parameters, \emph{i.e.}, $\Lambda\left((\bI - \bB)^{-1}\bkappa\right)$, (iv) AIC and BIC information criteria in addition to the the log--likelihood evaluated at its optimum, and (v) the computation time.
Concerning the estimated coefficients, \code{kappa1}, \code{kappa2} and \code{kappa3} are the elements of vector $\bkappa$ in~\eqref{eq:update_studentt}, \emph{i.e.}, $\kappa_\mu, \kappa_\phi$ and $\kappa_\nu$, respectively. Analogously, \code{a1} and \code{a2} are the estimates of $a_\mu$ and $a_\phi$ and \code{b1} and \code{b2} are estimates of $b_\mu$ and $b_\phi$. Note that, since we have specified \code{scale = FALSE} in the \code{UniGASSpec()} function, coefficients \code{a3} and \code{b3}, corresponding to $a_\nu$ and $b_\nu$ are not reported (and constrained to zero during the log--likelihood optimization).
Several methods to extract estimated quantities are available for objects of the class \code{uGASFit} or \code{mGASFit}. They allow us to: (i) calculate several quantities of the estimated conditional distribution at each point in time, such as: quantiles, conditional moments and filtered parameters (see \code{quantile(Fit)}, \code{getMoments(Fit)} and \code{getFilteredParameters(Fit)}, respectively), (ii) extract the estimated coefficients (\code{coef(Fit)}), (iii) generate a graphical representation of the results (\code{plot(Fit)}); see \code{help(uGASFit)} for details.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsection[Forecasting]{Forecasting}\label{sec:gas_fore}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Forecasting is a crucial aspect in applied time series analysis. Given the parametric assumption of GAS models, predictions are usually given in the form of density forecasts, \emph{i.e.}, the distribution of $\by_{T+h}\vert\by_{1:T}$ for $h\geq1$. Knowing the predictive density, practitioners can extract any relevant quantities such as future expected value $\E_T\left[\by_{T+h}\right]$ or (co--)variance $\V_T\left[\by_{T+h}\right]$. 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 \proglang{R} package \pkg{GAS} can handle both one--step and multi--step ahead forecasts. Consistent with previous nomenclature, functions for univariate and multivariate predictions are \code{UniGASFor()} and \code{MultiGASFor()}, respectively. These
functions accept an object of the class \code{uGASFit} or \code{mGASFit}, created using the functions \code{UniGASFit()} and \code{MultiGASFit()}, and return an object of the class \code{uGASFor} and \code{mGASFor}, respectively. Additional arguments are:
%
\begin{itemize}
\item \code{H}: a \code{numeric} integer value representing the forecast horizon, \emph{i.e.}, $h$. By default \code{iH = 1}.
\item \code{B}: a \code{numeric} integer value representing the number of draws to approximate the multi--step ahead predictive distribution when $h>1$. By default, \code{B = 1e4}.
\item \code{ReturnDraws}: a \code{boolean} controlling if the simulated draws from $\by_{T+1}\vert \by_{1:T}$, $\by_{T+2}\vert \by_{1:T},\dots$, $\by_{T+h}\vert \by_{1:T}$ have to be returned. By default \code{ReturnDraws = FALSE}.
\end{itemize}
%
Other arguments to perform rolling--type of forecasts are detailed in the documentation; see \code{help("UniGASFor")}. Practically, if we want to predict the next year inflation (\emph{i.e.}, $h=12$ with the monthly series \code{cpichg}), after having estimated the GAS model of Section~\ref{sec:gas_est}, we can execute the following code:
%
\begin{CodeChunk}
\begin{CodeInput}
R> Forecast <- UniGASFor(Fit, H = 12)
\end{CodeInput}
\end{CodeChunk}
%
and inspect the results by calling the object \code{Forecast}:
%
\begin{CodeChunk}
\begin{CodeInput}
R> Forecast
\end{CodeInput}
\begin{CodeOutput}
------------------------------------------
- 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.10219 0.1753 6.526
T+3 0.09718 0.2188 6.526
T+4 0.10207 0.2586 6.526
T+5 0.09407 0.2987 6.526
....................
location scale shape
T+8 0.09603 0.4207 6.526
T+9 0.08959 0.4458 6.526
T+10 0.08804 0.4901 6.526
T+11 0.07340 0.5159 6.526
T+12 0.07060 0.5405 6.526
\end{CodeOutput}
\end{CodeChunk}
%
which returns some model information and the predictions of future model parameters based on averages over \code{B} draws. \code{Forecast} is an object of the class \code{uGASFor} and comes with several methods to extract and visualize the results; see \code{help("uGASFor")}.
As commonly done in (financial) econometrics, predictions are generated from models fitted to rolling windows. The \proglang{R} package \pkg{GAS} includes this functionality via \code{UniGASRoll()} and \code{MultiGASRoll()}. These functions accept several arguments that we briefly describe in the univariate case:
%
\begin{itemize}
\item \code{data}: a vector of length $T + $\code{ForecastLength} containing all the observations.
\item \code{GASSpec}: an object of the class \code{uGASSpec} created with \code{UniGASSpec()}.
\item \code{ForecastLength}: a \code{numeric} integer which specifies the length of the out--of--sample.
\item \code{RefitEvery}: a \code{numeric} integer of periods before coefficients are re--estimated.
\item \code{RefitWindow}: a \code{character} for the type of the window. As in the \proglang{R} package \pkg{rugarch}, we define the options: \code{RefitWindow = "recursive"} and \code{RefitWindow = "moving"}. If \code{RefitWindow = "recursive"} all past observations are used when the model is re--estimated. If \code{RefitWindow = "moving"} initial observations are eliminated.
\end{itemize}
%
Other arguments useful to tailor the forecasting procedure and to parallelize the code execution are available and detailed in the \proglang{R} documentation; see \code{help("UniGASRoll")}.
Suppose now we are interested in assessing the forecast ability of the GAS model with \Student conditional distribution and time varying location and scale parameters, detailed in Section~\ref{sec:gas_t}, and specified in the object \code{GASSpec} in Appendix~\ref{sec:gas_spec}. We treat the last 150 observations of \code{cpichg} as out--of--sample and run a rolling--window forecast exercise using the following portion of code:
%
\begin{CodeChunk}
\begin{CodeInput}
Roll <- UniGASRoll(cpichg, GASSpec, ForecastLength = 150,
RefitEvery = 3, RefitWindow = "moving")
\end{CodeInput}
\end{CodeChunk}
%
where model coefficients are re--estimated quarterly (\emph{i.e.}, every three observations). 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 \code{Roll} belongs to the class \code{uGASRoll} which, as \code{uGASFit} and \code{uGASFor}, comes with several methods to extract and represent
the results; see \code{help("uGASRoll")}.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsection[Simulation]{Simulation}\label{sec:gas_sim}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Simulation of univariate and multivariate GAS models is straightforward with the \proglang{R} package \pkg{GAS}. This can be easily done via \code{UniGASSim()} and \code{MultiGASSim()}; see the \proglang{R} documentation. Several examples, also investigating the finite sample properties of the ML estimator for GAS models, are reported in the \code{inst/test/Simulation.R} file included in the package tarball. For simulation of GAS models, the vector $\bkappa$ and the matrices $\bB$ and $\bA$, need to be specified. It is worth stressing that, the definition of $\bkappa$ can be tricky, especially for multivariate models. The difficulty emerges from the fact that, $\bkappa$ determines the unconditional value of the reparametrised vector of parameters $\widetilde\btheta_t$, implying that, if the user wants to specify the model in terms of a target value $\btheta^*$ she needs to know the inverse of the mapping function $\Lambda\left(\cdot\right)$.\footnote{Here we define the \qmo target value\qmcsp as the unconditional level of the parameters the user has in mind, however, in order to have $\E\left[\btheta_t\right] = \btheta^*$, the model has to be stationary, see \emph{e.g.}, \cite{blasques_etal.2014e}.} To address this problem, the functions \code{UniUnmapParameters()} and \code{MultiUnmapParameters()}, representing $\Lambda^{-1}\left(\cdot\right)$, are available for univariate and multivariate models, respectively. This way, the user can easily specify $\bkappa$ such that $\bkappa = \left(\bI - \bB\right)\Lambda^{-1}\left(\btheta^*\right)$. Table~\ref{tab:bounds} lists the numerical bounds imposed for the univariate distributions, such that \code{UniUnmapParameters()} cannot takes values outside those ranges. For the multivariate case, please refer to the examples reported in the \code{inst/test/SimulateGAS.R} file included in the package tarball.
Suppose we want to simulate $T = 1,000$ observations from the \Student GAS model reported in Appendix~\ref{sec:gas_t} with time--varying location and scale, but constant shape parameters. Assume our target value for the parameters is $\btheta^* = \left(\mu^*, \sigma^*, \nu^*\right)^\prime$ with $\mu^* = 0.1, \sigma^* = 1.5$ and $\nu^* = 7$. The matrix $\bA$ and $\bB$ are defined as:
\begin{align}\label{eq:sim_a_b}
\bA &= \diag\left(0.1, 0.4, 0.0\right)\\
\bB &= \diag\left(0.9, 0.95, 0.0\right),
\end{align}
such that both the conditional mean and the conditional variance evolve quite persistently over time, while the shape parameter is constant. The implementation of \code{UniUnmapParameters()} and \code{UniGASSim()} proceed as:
\begin{CodeChunk}
\begin{CodeInput}
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(iT = 1000, Kappa, A, B,
Dist = "std", ScalingType = "Identity")
\end{CodeInput}
\end{CodeChunk}
\begin{table}[!t]
\centering
\setlength\tabcolsep{25pt}
\begin{tabular}{ccccc}
\toprule
Label & \code{location} & \code{scale} & \code{skewness} & \code{shape} \\
\hline
\code{norm} & $\Re^{~~}$ & $\Re^+$ & -- & -- \\
\code{std} & $\Re^{~~}$ & $\Re^+$ & -- & $\left(4.0,50.0\right)$ \\
\code{snorm} & $\Re^{~~}$ & $\Re^+$ & $\left(0.1,2.0\right)$ & -- \\
\code{sstd} & $\Re^{~~}$ & $\Re^+$ & $\left(0.1,2.0\right)$ & $\left(4.0,50.0\right)$ \\
\code{ald} & $\Re^{~~}$ & $\Re^+$ & $\Re^+$ & -- \\
\code{poi} & $\Re^+$ & -- & -- & -- \\
\code{gamma} & $\Re^+$ & $\Re^+$ & -- & -- \\
\code{exp} & $\Re^+$ & -- & -- & -- \\
\code{beta} & $\Re^+$ & $\Re^+$ & -- & -- \\
\code{norm} & $\Re^{~~}$ & $\Re^+$ & -- & -- \\
\bottomrule
\end{tabular}
\caption{Numerical bounds for univariate distributions in the \proglang{R} package \pkg{GAS}. For the case when the space is $\Re^+$ we use the exponential link function reported in Equation~\ref{eq:exponential_link} with $c=0$, while in the case when the space is of the type $(a, b)$ we use the modified logistic link function reported in Equation~\ref{eq:logistic_link}.}
\label{tab:bounds}
\end{table}
\vspace{-1.5cm}
where \code{Sim} is an object of the class \code{uGASSim} and comes with several method such as \code{show}, \code{plot}, \code{getMoments}, among others; see \code{help("uGASSim")}.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section[Applications to financial data]{Applications to financial data}\label{sec:application}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
In order to illustrate how the \proglang{R} package \pkg{GAS} can be used in practical situations, we report an empirical application with univariate and multivariate time series of financial returns. We consider daily log--returns (in percentage points) of the Dow Jones 30 constituents available in the \code{dji30ret} dataset originally included in the \pkg{rugarch} package and made available also in the \pkg{GAS} package. This dataset includes the closing value log--returns from March 3rd, 1987 to February 3rd, 2009 for a total of 5,521 observations per series. The dataset can be easily loaded in the workspace using:
\begin{CodeChunk}
\begin{CodeInput}
R> library("GAS")
R> data("dji30ret", package = "GAS")
\end{CodeInput}
\end{CodeChunk}
where \code{dji30ret} is a $5521 \times 30$ \code{data.frame} containing the daily logarithmic 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 during the out--of--sample period, (iii) we compare the models according to their out--of--sample performance.
The models we consider are univariate/multivariate GAS models estimated with the \proglang{R} package \pkg{GAS}, and univariate/multivariate GARCH models estimated using the popular \proglang{R} packages \pkg{rugarch} \citep{rugarch} and \pkg{rmgarch} \citep{rmgarch}, respectively. The univariate specifications we consider are: (i) the \SStudent GAS model with only time--varying scale parameter (\emph{i.e.}, \code{Dist = "sstd"}) and, (ii) the GARCH(1,1) model with \SStudent distributed error. For both models we employ the \SStudent distribution of \cite{fernandez_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 \pkg{rugarch} packages.
For the multivariate specifications, we consider: (i) the GAS model with conditional multivariate \Student distribution with time--varying scales and correlations used in \citet{creal_etal.2011} and, (ii) the Dynamic Conditional Correlation (DCC) model of \citet{engle.2002} with a conditional multivariate \Student distribution. For simplicity, the multivariate analysis only considers three series of the whole dataset, these are: Caterpillar Inc. (CAT), 3M (MMM) and Pfizer Inc. (PFE).
The code used to specify the univariate and multivariate GAS models is:
\begin{CodeChunk}
\begin{CodeInput}
R> uGASSpec <- UniGASSpec(Dist = "sstd",
ScalingType = "Identity",
GASPar = list(scale = TRUE))
\end{CodeInput}
\end{CodeChunk}
and:
\begin{CodeChunk}
\begin{CodeInput}
R> mGASSpec <- MultiGASSpec(Dist = "mvt",
ScalingType = "Identity",
GASPar = list(scale = TRUE,
correlation = TRUE))
\end{CodeInput}
\end{CodeChunk}
respectively.
% Experiment design
The last $H = 3,000$ observations (from January 27th, 1991, to the end of the sample) compose 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 therefore coefficients) are re--estimated using a moving--window every hundredth observations, as detailed in Section~\ref{sec:gas_fore}. One--step ahead rolling prediction are then computed as:
\begin{CodeChunk}
\begin{CodeInput}
luGASRoll <- list()
N <- ncol(dji30ret)
for(i in 1:N){
luGASRoll[[i]] <- UniGASRoll(data = dji30ret[, i],
GASSpec = uGASSpec,
ForecastLength = 3000,
RefitEvery = 100)
}
names(luGASRoll) <- colnames(dji30ret)
\end{CodeInput}
\end{CodeChunk}
and:
\begin{CodeChunk}
\begin{CodeInput}
mGASRoll <- MultiGASRoll(data = dji30ret[, c("CAT", "MMM", "PFE")],
GASSpec = mGASSpec,
ForecastLength = 3000,
RefitEvery = 100)
\end{CodeInput}
\end{CodeChunk}
for the univariate and multivariate case, respectively.
Let us now compare the ability of GAS and GARCH models in predicting the one--step ahead distribution using two \textit{proper} scoring rules \citep{gneiting_etal.2007}. Scoring rules compare the predicted density with the ex--post realized value of the return and deliver a score which defines a ranking across the alternative models at each point in time. Generally, we define $S_{t+1}\bydef S(y_{t+1}, p(y_{t+1};\widehat\btheta_{t+1}))$ as the score at time $t+1$ for having predicted $p(y_{t+1};\widehat\btheta_{t+1})$ when $y_{t+1}$ has been realized. We consider two widely used scoring rules:
\begin{itemize}
\item The average Negative Log Score (NLS):
\begin{equation}\label{eq:logscore}
\overline{NLS} \bydef -\frac{1}{H}\sum_{t = T}^{T+H-1} \log p(y_{t+1};\widehat\btheta_{t+1}).
\end{equation}
\item The average weighted Continuous Ranked Probability Score (wCRPS):
\begin{equation}\label{eq:crsp}
\overline{wCRPS} \bydef \frac{1}{H}\sum_{t = T}^{T+H-1}\int_{-\infty}^{\infty}w\left(z\right)\left(F\left(z;\widehat\btheta_{t+1}\right) - \mathbb{I}_{\{y_{t+1}<z\}}\right)^2\mathrm{d}z,
\end{equation}
where $w\left(z\right)$ is a weight function that emphasizes regions of interest of the predictive distribution, such as the tails or the center. Models with lower $\overline{NLS}$ and $\overline{wCRPS}$ are preferred.\footnote{Consistent with \citet{gneiting_etal.2007}, we specify the Negative Log Score such that the \qmo direction\qmcsp between the two scoring rules is the same.} Similarly to \citet{gneiting_ranjan.2012}, we consider the cases:
\begin{itemize}
\item[(i)] $w\left(z\right) = 1$: Uniform,
\item[(ii)] $w\left(z\right) = \phi_{a,b}\left(z\right)$ : Center,
\item[(iii)] $w\left(z\right) = 1 - \phi_{a,b}\left(z\right)/\phi_{a,b}\left(0\right)$: Tails,
\item[(iv)] $w\left(z\right) = \Phi_{a,b}\left(z\right)$: Right tail,
\item[(v)] $w\left(z\right) = 1 - \Phi_{a,b}\left(z\right)$: Left tail,
\end{itemize}
where $\phi_{a,b}\left(z\right)$ and $\Phi_{a,b}\left(z\right)$ are the \textit{pdf} and \textit{cdf} of a Gaussian distribution with mean $a$ and standard deviation $b$, respectively. The label \qmo Uniform\qmcsp represents the case where equal emphasis is given to all the parts of the distribution.
\end{itemize}
The two aforementioned scoring rules can be easily evaluated using the \code{BacktestDensity()} function available in the \proglang{R} package \pkg{GAS}. The \code{BacktestDensity()} function accepts an object of the class \code{uGASRoll}, and returns a list with two elements: (i) the averages negative LS and wCRPS as in~\eqref{eq:logscore} and~\eqref{eq:crsp}, and (ii) their values at each point in time. Additional arguments are:
\begin{itemize}
\item \code{lower}: \code{numeric} representing the lower bound used to approximate Equation~\ref{eq:crsp} by Monte Carlo integration as detailed in \cite{gneiting_ranjan.2012}.
\item \code{upper}: \code{numeric} as \code{lower} but for the upper bound.\footnote{The two arguments \code{lower} and \code{upper} coincide with $y_l$ and $y_u$ in Equation 16 of \cite{gneiting_ranjan.2012}, respectively. These are two numeric objects with no default value, \emph{i.e.}, the user have to define these values according to her research design.}
\item \code{K}: \code{numeric} integer representing the number of points used to discretize the integral in Equation~\ref{eq:crsp}.\footnote{Equals to $I$ in Equation 16 of \cite{gneiting_ranjan.2012}.} By default, \code{K = 1000},
\end{itemize}
plus the two \code{numeric} arguments, \code{a} and \code{b}, representing $a$ and $b$ in the weight functions, by default \code{a = 0.0} and \code{b = 1.0}.\footnote{These values can be chosen in order to target some \qmo optimal\qmcsp prediction level, or to add more flexibility and focus on specific part the predictive distribution; see \cite{gneiting_ranjan.2012}.}
\begin{table}[!t]
\centering
\setlength\tabcolsep{15pt}
%\resizebox{1\columnwidth}{!}{%
\begin{tabular}{lrrrrrr}
\toprule
Asset & LS & uniform & center & tails & tails--r & tails--l\\
\cmidrule(lr){1-1}\cmidrule(lr){2-7}
AA & $-2.12^b$ & $-2.18^b$ & $-2.18^b$ & $-0.69^{~}$ & $-2.26^b$ & $-2.09^b$ \\
AXP & $-3.06^a$ & $-2.54^a$ & $-2.54^a$ & $-0.81^{~}$ & $-2.59^a$ & $-2.48^a$ \\
BA & $-1.62^c$ & $-2.06^b$ & $-2.07^b$ & $-0.39^{~}$ & $-2.07^b$ & $-2.06^b$ \\
BAC & $-2.85^a$ & $-1.05^{~}$ & $-1.06^{~}$ & $0.64^{~}$ & $-0.94^{~}$ & $-1.15^{~}$ \\
C & $-4.35^a$ & $-2.83^a$ & $-2.84^a$ & $-0.28^{~}$ & $-3.02^a$ & $-2.63^a$ \\
CAT & $-5.84^a$ & $-5.96^a$ & $-5.96^a$ & $-2.86^a$ & $-5.91^a$ & $-6.00^a$ \\
CVX & $0.57^{~}$ & $1.14^{~}$ & $1.14^{~}$ & $1.33^c$ & $1.13^{~}$ & $1.14^{~}$ \\
DD & $-2.40^a$ & $-1.86^b$ & $-1.86^b$ & $-1.57^c$ & $-1.84^b$ & $-1.89^b$ \\
DIS & $-2.58^a$ & $-2.30^b$ & $-2.30^b$ & $-1.19^{~}$ & $-2.31^b$ & $-2.28^b$ \\
GE & $-3.82^a$ & $-4.64^a$ & $-4.64^a$ & $-2.50^a$ & $-4.69^a$ & $-4.58^a$ \\
GM & $0.50^{~}$ & $0.11^{~}$ & $0.10^{~}$ & $0.90^{~}$ & $0.03^{~}$ & $0.19^{~}$ \\
HD & $-3.84^a$ & $-3.05^a$ & $-3.05^a$ & $-3.94^a$ & $-3.07^a$ & $-3.03^a$ \\
HPQ & $-4.73^a$ & $-4.90^a$ & $-4.90^a$ & $-3.84^a$ & $-4.88^a$ & $-4.92^a$ \\
IBM & $-4.71^a$ & $-4.67^a$ & $-4.67^a$ & $-3.28^a$ & $-4.64^a$ & $-4.69^a$ \\
INTC & $-3.05^a$ & $-1.72^b$ & $-1.72^b$ & $-1.23^{~}$ & $-1.79^b$ & $-1.64^c$ \\
JNJ & $-3.95^a$ & $-2.26^b$ & $-2.26^b$ & $-0.46^{~}$ & $-2.19^b$ & $-2.34^a$ \\
JPM & $-2.21^b$ & $-2.20^b$ & $-2.20^b$ & $-0.67^{~}$ & $-2.14^b$ & $-2.25^b$ \\
AIG & $-0.38^{~}$ & $0.89^{~}$ & $0.89^{~}$ & $0.64^{~}$ & $0.68^{~}$ & $1.02^{~}$ \\
KO & $-3.75^a$ & $-3.72^a$ & $-3.72^a$ & $-2.02^b$ & $-3.74^a$ & $-3.69^a$ \\
MCD & $-2.20^b$ & $-1.94^b$ & $-1.94^b$ & $-1.15^{~}$ & $-1.94^b$ & $-1.93^b$ \\
MMM & $-4.69^a$ & $-5.24^a$ & $-5.24^a$ & $-2.97^a$ & $-5.25^a$ & $-5.24^a$ \\
MRK & $-3.83^a$ & $-4.90^a$ & $-4.90^a$ & $-3.63^a$ & $-4.98^a$ & $-4.81^a$ \\
MSFT & $-3.95^a$ & $-3.23^a$ & $-3.23^a$ & $-2.44^a$ & $-3.26^a$ & $-3.19^a$ \\
PFE & $-5.02^a$ & $-5.22^a$ & $-5.22^a$ & $-3.56^a$ & $-5.27^a$ & $-5.18^a$ \\
PG & $-2.24^b$ & $-2.13^b$ & $-2.13^b$ & $-2.21^b$ & $-2.15^b$ & $-2.11^b$ \\
T & $-0.09^{~}$ & $-0.01^{~}$ & $-0.01^{~}$ & $0.42^{~}$ & $0.00^{~}$ & $-0.03^{~}$ \\
UTX & $-1.80^b$ & $-1.73^b$ & $-1.73^b$ & $-0.95^{~}$ & $-1.78^b$ & $-1.67^b$ \\
VZ & $-2.58^a$ & $-2.48^a$ & $-2.48^a$ & $-3.35^a$ & $-2.47^a$ & $-2.49^a$ \\
WMT & $-2.74^a$ & $-2.47^a$ & $-2.47^a$ & $-0.36^{~}$ & $-2.52^a$ & $-2.43^a$ \\
XOM & $0.03^{~}$ & $0.33^{~}$ & $0.33^{~}$ & $1.27^{~}$ & $0.35^{~}$ & $0.30^{~}$ \\
\bottomrule
\end{tabular}
%}
\caption{Test statistics for the \cite{diebold_mariano.1995} test between 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 predictions of the one--step ahead conditional distribution while positive values favour GARCH. The apexes $a,b$ and $c$ represents rejection of the null hypothesis of Equal Predictive Ability at the $1\%, 5\%$ and $10\%$ confidence levels, respectively. The out--of--sample period spans from January 27th, 1991, to February 3rd, 2009 for a total of 3,000 observations.}
\label{tab:BackDist}
\end{table}
In our case, in order to evaluate $\overline{NLS}$ and $\overline{wCRPS}$ for the first asset we can simply run:
\begin{CodeChunk}
\begin{CodeInput}
R> DensityBacktest <- BacktestDensity(luGASRoll[[1]],
lower = -1.0, upper = 1.0)
R> DensityBacktest$average
\end{CodeInput}
\begin{CodeOutput}
LS uniform center tails tail_r tail_l
-2.387e+00 1.331e-02 5.307e-03 7.353e-06 6.654e-03 6.658e-03
\end{CodeOutput}
\end{CodeChunk}
where \code{lower = -1.0} and \code{upper = 1.0} works well for log--returns not in percentage points as the one considered here.
Table~\ref{tab:BackDist} reports the test statistics for the \citet{diebold_mariano.1995} (DM) test between the series of Negative Log Scores and weighed Continuous Ranked Probability Scores for univariate GAS and GARCH models across the out--of--sample period. Negative values indicate that GAS models generate more accurate predictions of the one--step ahead conditional distribution while positive values favours GARCH. We found that, for almost all the series, GAS outperforms GARCH at very high confidence 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 the average Negative Log Score, $\overline{NLS}$. In this case, the DM test statistic is $-0.34$, which still favours the GAS model against the DCC specification, however without being statistically different from zero (the associated $p$--value is about $0.36$). To further investigate this result, we report in Figure~\ref{fig:cumlogscore} the Cumulative sum of the differences between the Log Scores (CLS) of GAS and DCC defined as:
\begin{equation}\label{eq:diff}
CLS_{T:T+l}^{GAS\vert DCC} \bydef \sum_{t = T}^{t = T+l-1}\log p\left(\by_{t+1};\widehat{\btheta}_{t+1}^{GAS}\right) - \log p\left(\by_{t+1};\widehat{\btheta}_{t+1}^{DCC}\right)\,,
\end{equation}
where $p\left(\by_{t+1};\widehat{\btheta}_{t+1}^{GAS}\right)$ and $p\left(\by_{t+1};\widehat{\btheta}_{t+1}^{DCC}\right)$ are the densities predicted from GAS and DCC evaluated in $\by_{t+1}$, respectively. The series of Log Scores for the multivariate GAS models is available in the output of the \code{BacktestDensity()} function, or can be extracted using the \code{LogScore} method defined for \code{mGASRoll} objects:
\begin{CodeChunk}
\begin{CodeInput}
R> LS_MGAS <- LogScore(mGASRoll)
\end{CodeInput}
\end{CodeChunk}
Looking at Figure~\ref{fig:cumlogscore}, 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 dominates DCC for the first part of the out--of--sample period after the end of the early 2000's dot--com bubble, while the converse holds for the last part of the sample. This result suggests that possible model combinations between GAS and DCC can improve the predictions for multivariate returns series. We left this as a further research topic.
\begin{figure}[!t]
\centering
\includegraphics[width=1\textwidth]{CumLogScore_Multi}
\caption{Cumulative out--of--sample Log Score difference between the multivariate Student--t GAS and the DCC(1,1) model of \citet{engle.2002} with multivariate \Student distributed 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 for US according to the \qmo USREC\qmcsp series available from the Federal Reserve Bank of St. Louis web site: \url{https://fred.stlouisfed.org/series/USREC}.}
%%%
\label{fig:cumlogscore}
\end{figure}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section[Conclusion]{Conclusion}\label{sec:conclusion}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% it is important
This article introduced the \proglang{R} package \pkg{GAS} for simulating, estimating and forecasting with Generalized Autoregressive Score models in \proglang{R}. We briefly reviewed GAS models and the fast growing literature that relies on this class of models in applied works. To our knowledge, the \proglang{R} package \pkg{GAS} is the first complete framework for GAS models allowing practitioners in many scientific areas to perform their applied research in an user--friendly environment.
% the empirical application
We introduced the model specification in a general way and illustrated the package's usage. In particular, we performed and empirical application to financial data in which we compared the performance of univariate and multivariate GAS and GARCH models.
% what you can do with it
Given the flexibility of GAS models and the availability of several statistical distributions in the \pkg{GAS} package, a number of different applications can be easily handled, such as: (i) the analysis of integer valued time series using the Poisson GAS model (\code{poi}), (ii) the analysis of (0,1)--bounded time series using the Beta GAS model (\code{beta}), (iii) the analysis of strictly positive time series with an inverse location/scale dependence using the Gamma GAS model (\code{gamma}). We
are convinced that GAS models are relevant in many fields and hope that the \proglang{R} package \pkg{GAS} will be fruitful for many
researchers like econometricians or applied statisticians.
Finally, if you use \proglang{R} or \pkg{GAS}, please cite the software in publications.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section*{Computational details}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
The results in this paper were obtained using \proglang{R} 3.2.3 \citep{R} with the
packages: \pkg{GAS} version 0.1.1 \citep{GAS}, \pkg{MASS} version 7.3-45 and \citep{MASS.2002,MASS},
\pkg{Rcpp} version 0.12.5 \citep{Rcpp.2011,Rcpp}, \pkg{RcppArmadillo} version 0.7.100.3.1 \citep{RcppArmadillo.2014,RcppArmadillo}, \pkg{Rsolnp} version 1.15 \citep{Rsolnp}, \pkg{xts} version 0.9-7 \citep{xts} and \pkg{quantmod} version 0.4-5 \citep{quantmod}.
\proglang{R} itself and all packages used are available
from \proglang{CRAN} at \url{http://CRAN.R-project.org/}. The
package \pkg{GAS} is available from the CRAN repository at \url{https://cran.r-project.org/package=GAS}. The version under development is available in GitHub at \url{https://github.com/LeopoldoCatania/GAS}.
Computations were performed on
a Genuine Intel\textregistered{} quad core CPU i7--3630QM 2.40Ghz processor. Code outputs were
obtained using \code{options(digits = 4, max.print = 40, prompt = "R> ")}
The folder \code{inst/doc} inside the \pkg{GAS} package tarball contains additional technical documentations. A step by step guide on how to add a new statistical distribution in the \pkg{GAS} package is reported in the file \code{AddNewDistribution.pdf}.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section*{Acknowledgments}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
The authors acknowledge Google for financial support via the
Google Summer of Code 2016 project
"GAS"; see \url{https://summerofcode.withgoogle.com/projects/#4717600793690112}.
Any remaining errors or shortcomings are the authors' responsibility.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% JSS has its own style
\begin{thebibliography}{46}
\newcommand{\enquote}[1]{``#1''}
\providecommand{\natexlab}[1]{#1}
\providecommand{\url}[1]{\texttt{#1}}
\providecommand{\urlprefix}{URL }
\expandafter\ifx\csname urlstyle\endcsname\relax
\providecommand{\doi}[1]{doi:\discretionary{}{}{}#1}\else
\providecommand{\doi}{doi:\discretionary{}{}{}\begingroup
\urlstyle{rm}\Url}\fi
\providecommand{\eprint}[2][]{\url{#2}}
\bibitem[{Andres(2014)}]{andres.2014}
Andres P (2014).
\newblock \enquote{Maximum Likelihood Estimates for Positive Valued Dynamic
Score Models: {T}he \pkg{DySco} Package.}
\newblock \emph{Computational Statistics \& Data Analysis}, \textbf{76},
34--42.
\newblock \doi{10.1016/j.csda.2013.11.004}.
\bibitem[{Blasques \emph{et~al.}(2014{\natexlab{a}})Blasques, Koopman, and
Lucas}]{blasques_etal.2014b}
Blasques F, Koopman SJ, Lucas A (2014{\natexlab{a}}).
\newblock \enquote{Maximum Likelihood Estimation for Correctly Specified
Generalized Autoregressive Score Models: {F}eedback Effects, Contraction
Conditions and Asymptotic Properties.}
\newblock \emph{techreport TI 14-074/III}, Tinbergen Institute.
\newblock \urlprefix\url{http://www.tinbergen.nl/discussionpaper/?paper=2332}.
\bibitem[{Blasques \emph{et~al.}(2014{\natexlab{b}})Blasques, Koopman, and
Lucas}]{blasques_etal.2014c}
Blasques F, Koopman SJ, Lucas A (2014{\natexlab{b}}).
\newblock \enquote{Maximum Likelihood Estimation for Generalized Autoregressive
Score Models.}
\newblock \emph{techreport TI 2014-029/III}, Tinbergen Institute.
\newblock \urlprefix\url{http://www.tinbergen.nl/discussionpaper/?paper=2286}.
\bibitem[{Blasques \emph{et~al.}(2014{\natexlab{c}})Blasques, Koopman, Lucas,
and Schaumburg}]{blasques_etal.2014a}
Blasques F, Koopman SJ, Lucas A, Schaumburg J (2014{\natexlab{c}}).
\newblock \enquote{Spillover Dynamics for Systemic Risk Measurement using
Spatial Financial Time Series Models.}
\newblock \emph{techreport TI 2014-103/III}, Tinbergen Institute.
\newblock \urlprefix\url{http://www.tinbergen.nl/discussionpaper/?paper=2369}.
\bibitem[{Blasques \emph{et~al.}(2014{\natexlab{d}})Blasques, Koopman, Lucas
\emph{et~al.}}]{blasques_etal.2014e}
Blasques F, Koopman SJ, Lucas A, \emph{et~al.} (2014{\natexlab{d}}).
\newblock \enquote{Stationarity and Ergodicity of Univariate Generalized
Autoregressive Score Processes.}
\newblock \emph{Electronic Journal of Statistics}, \textbf{8}(1), 1088--1112.
\newblock \doi{doi:10.1214/14-EJS924}.
\bibitem[{Bollerslev(1987)}]{bollerslev.1987}
Bollerslev T (1987).
\newblock \enquote{A Conditionally Heteroskedastic Time Series Model for
Speculative Prices and Rates of Return.}
\newblock \emph{The Review of Economics and Statistics}, \textbf{69}(3),
542--547.
\newblock \doi{10.2307/1925546}.
\bibitem[{Box and Jenkins(1970)}]{box_jenkins.1970}
Box GE, Jenkins GM (1970).
\newblock \emph{Time Series Analysis: {F}orecasting and Control}.
\newblock Holden--Day, San Francisco.
\bibitem[{Catania and {Bill{\'e}}(2016)}]{catania_bille.2016}
Catania L, {Bill{\'e}} AG (2016).
\newblock \enquote{Dynamic Spatial Autoregressive Models with Autoregressive
and Heteroskedastic Disturbances.}
\newblock Working paper, \urlprefix\url{http://arxiv.org/abs/1602.02542}.
\bibitem[{Catania \emph{et~al.}(2016)Catania, Boudt, and Ardia}]{GAS}
Catania L, Boudt K, Ardia D (2016).
\newblock \emph{\pkg{GAS}: {G}eneralised Autoregressive Score Models}.
\newblock \proglang{R}~package version 0.1.1,
\urlprefix\url{https://cran.r-project.org/package=GAS}.
\bibitem[{Creal \emph{et~al.}(2011)Creal, Koopman, and Lucas}]{creal_etal.2011}
Creal D, Koopman SJ, Lucas A (2011).
\newblock \enquote{A Dynamic Multivariate Heavy--Tailed Model for Time--Varying
Volatilities and Correlations.}
\newblock \emph{Journal of Business \& Economic Statistics}, \textbf{29}(4),
552--563.
\newblock \doi{10.1198/jbes.2011.10070}.
\bibitem[{Creal \emph{et~al.}(2013)Creal, Koopman, and Lucas}]{creal_etal.2013}
Creal D, Koopman SJ, Lucas A (2013).
\newblock \enquote{Generalized Autoregressive Score Models with Applications.}
\newblock \emph{Journal of Applied Econometrics}, \textbf{28}(5), 777--795.
\newblock \doi{10.1002/jae.1279}.
\bibitem[{Creal \emph{et~al.}(2014)Creal, Schwaab, Koopman, and
Lucas}]{creal_etal.2014}
Creal D, Schwaab B, Koopman SJ, Lucas A (2014).
\newblock \enquote{Observation--Driven Mixed--Measurement Dynamic Factor Models
with an Application to Credit Risk.}
\newblock \emph{The Review of Economics and Statistics}, \textbf{96}(5),
898--915.
\newblock \doi{10.1162/REST_a_00393}.
\bibitem[{Diebold and Mariano(1995)}]{diebold_mariano.1995}
Diebold FX, Mariano RS (1995).
\newblock \enquote{Comparing Predictive Accuracy.}
\newblock \emph{Journal of Business \& Economic Statistics}, \textbf{13}(3),
253--263.
\newblock \doi{10.1080/07350015.1995.10524599}.
\bibitem[{Eddelbuettel and Fran\c{c}ois(2011)}]{Rcpp.2011}
Eddelbuettel D, Fran\c{c}ois R (2011).
\newblock \enquote{\pkg{Rcpp}: {S}eamless \proglang{R} and \proglang{C++}
Integration.}
\newblock \emph{Journal of Statistical Software}, \textbf{40}(8), 1--18.
\newblock \doi{10.18637/jss.v040.i08}.
\bibitem[{Eddelbuettel \emph{et~al.}(2016{\natexlab{a}})Eddelbuettel,
Fran\c{c}ois, Allaire, Ushey, Kou, Bates, and Chambers}]{Rcpp}
Eddelbuettel D, Fran\c{c}ois R, Allaire J, Ushey K, Kou Q, Bates D, Chambers J
(2016{\natexlab{a}}).
\newblock \emph{\pkg{Rcpp}: {S}eamless \proglang{R} and \proglang{C++}
Integration}.
\newblock \proglang{R}~package version 0.12.5,
\urlprefix\url{https://cran.r-project.org/package=Rcpp}.
\bibitem[{Eddelbuettel \emph{et~al.}(2016{\natexlab{b}})Eddelbuettel,
Fran\c{c}ois, and Bates}]{RcppArmadillo}
Eddelbuettel D, Fran\c{c}ois R, Bates D (2016{\natexlab{b}}).
\newblock \emph{\pkg{RcppArmadillo}: \pkg{Rcpp} Integration for the
\pkg{Armadillo} Templated Linear Algebra Library}.
\newblock \proglang{R}~package version 0.7.100.3.1,
\urlprefix\url{https://cran.r-project.org/package=RcppArmadillo}.
\bibitem[{Eddelbuettel and Sanderson(2014)}]{RcppArmadillo.2014}
Eddelbuettel D, Sanderson C (2014).
\newblock \enquote{\pkg{RcppArmadillo}: {A}ccelerating \proglang{R} with
High--Performance \proglang{C++} Linear Algebra.}
\newblock \emph{Computational Statistics \& Data Analysis}, \textbf{71},
1054--1063.
\newblock \doi{10.1016/j.csda.2013.02.005}.
\bibitem[{Engle(2002)}]{engle.2002}
Engle RF (2002).
\newblock \enquote{Dynamic Conditional Correlation: {A} Simple Class of
Multivariate Generalized Autoregressive Conditional Heteroskedasticity
Models.}
\newblock \emph{Journal of Business \& Economic Statistics}, \textbf{20}(3),
339--350.
\newblock \doi{10.1198/073500102288618487}.
\bibitem[{Fern{\'a}ndez and Steel(1998)}]{fernandez_steel.1998}
Fern{\'a}ndez C, Steel MF (1998).
\newblock \enquote{On Bayesian modeling of fat tails and skewness.}
\newblock \emph{Journal of the American Statistical Association},
\textbf{93}(441), 359--371.
\newblock \doi{10.1080/01621459.1998.10474117}.
\bibitem[{Ghalanos(2015{\natexlab{a}})}]{rmgarch}
Ghalanos A (2015{\natexlab{a}}).
\newblock \emph{\pkg{rmgarch}: {M}ultivariate {GARCH} Models}.
\newblock \proglang{R}~package version 1.3-0,
\urlprefix\url{https://cran.r-project.org/package=rmgarch}.
\bibitem[{Ghalanos(2015{\natexlab{b}})}]{rugarch}
Ghalanos A (2015{\natexlab{b}}).
\newblock \emph{\pkg{rugarch}: {U}nivariate {GARCH} Models}.
\newblock \proglang{R}~package version 1.3-6,
\urlprefix\url{https://cran.r-project.org/package=rugarch}.
\bibitem[{Ghalanos and Theussl(2016)}]{Rsolnp}
Ghalanos A, Theussl S (2016).
\newblock \emph{\pkg{Rsolnp}: {G}eneral Non--Linear Optimization using
Augmented {L}agrange Multiplier Method}.
\newblock \proglang{R}~package version 1.16,
\urlprefix\url{https://cran.r-project.org/package=Rsolnp}.
\bibitem[{Gneiting \emph{et~al.}(2007)Gneiting, Balabdaoui, and
Raftery}]{gneiting_etal.2007}
Gneiting T, Balabdaoui F, Raftery AE (2007).
\newblock \enquote{Probabilistic forecasts, calibration and sharpness.}
\newblock \emph{Journal of the Royal Statistical Society: Series B (Statistical
Methodology)}, \textbf{69}(2), 243--268.
\newblock \doi{10.1111/j.1467-9868.2007.00587.x}.
\bibitem[{Gneiting and Ranjan(2012)}]{gneiting_ranjan.2012}
Gneiting T, Ranjan R (2012).
\newblock \enquote{Comparing density forecasts using threshold-- and
quantile--weighted scoring rules.}
\newblock \emph{Journal of Business \& Economic Statistics}.
\newblock \doi{10.1198/jbes.2010.08110}.
\bibitem[{Harvey(2013)}]{harvey.2013}
Harvey AC (2013).
\newblock \emph{Dynamic Models for Volatility and Heavy Tails: {W}ith
Applications to Financial and Economic Time Series}.
\newblock Cambridge University Press.
\bibitem[{Harvey and Luati(2014)}]{harvey_luati.2014}
Harvey AC, Luati A (2014).
\newblock \enquote{Filtering with Heavy Tails.}
\newblock \emph{Journal of the American Statistical Association},
\textbf{109}(507), 1112--1122.
\newblock \doi{10.1080/01621459.2014.887011}.
\bibitem[{Harvey and Sucarrat(2014)}]{harvey_sucarrat.2014}
Harvey AC, Sucarrat G (2014).
\newblock \enquote{EGARCH Models with Fat Tails, Skewness and Leverage.}
\newblock \emph{Computational Statistics \& Data Analysis}, \textbf{76},
320--338.
\newblock \doi{10.1016/j.csda.2013.09.022}.
\bibitem[{Harvey and Thiele(2015)}]{harvey_thiele.2014}
Harvey AC, Thiele S (2015).
\newblock \enquote{Testing Against Changing Correlation.}
\newblock \doi{10.1016/j.jempfin.2015.09.003}.
\newblock Forthcoming in Journal of Empirical Finance.
\bibitem[{Jaeckel and Rebonato(1999)}]{jaeckel_rebonato.1999}
Jaeckel P, Rebonato R (1999).
\newblock \enquote{The Most General Methodology for Creating a Valid
Correlation Matrix for Risk Management and Option Pricing Purposes.}
\newblock \emph{Journal of Risk}, \textbf{2}(2), 17--28.
\newblock
\urlprefix\url{http://www.risk.net/journal-of-risk/technical-paper/2161082/the-methodology-creating-valid-correlation-matrix-risk-management-option-pricing-purposes}.
\bibitem[{Janus \emph{et~al.}(2014)Janus, Koopman, and Lucas}]{janus_etal.2014}
Janus P, Koopman SJ, Lucas A (2014).
\newblock \enquote{Long Memory Dynamics for Multivariate Dependence under Heavy
Tails.}
\newblock \emph{Journal of Empirical Finance}, \textbf{29}, 187--206.
\newblock \doi{10.1016/j.jempfin.2014.09.007}.
\bibitem[{Kalman(1960)}]{kalman.1960}
Kalman RE (1960).
\newblock \enquote{A New Approach to Linear Filtering and Prediction Problems.}
\newblock \emph{Transactions of the ASME--Journal of Basic Engineering},
\textbf{82}(Series D), 35--45.
\bibitem[{Kotz \emph{et~al.}(2001)Kotz, Kozubowski, and
Podgorski}]{kotz_etal.2001}
Kotz S, Kozubowski T, Podgorski K (2001).
\newblock \emph{The {L}aplace Distribution and Generalizations: {A} Revisit
with Applications to Communications, Economics, Engineering, and Finance}.
\newblock Springer--Verlag.
\newblock \urlprefix\url{http://www.springer.com/us/book/9780817641665}.
\bibitem[{Lucas and Zhang(2016)}]{lucas_zhang.2016}
Lucas A, Zhang X (2016).
\newblock \enquote{Score--Driven Exponentially Weighted Moving Averages and
{V}alue--at--{R}isk Forecasting.}
\newblock \emph{International Journal of Forecasting}, \textbf{32}(2),
293--302.
\newblock \doi{10.1016/j.ijforecast.2015.09.003}.
\bibitem[{Oh and Patton(2013)}]{oh_patton.2013}
Oh DH, Patton AJ (2013).
\newblock \enquote{Time--Varying Systemic Risk: {E}vidence from a Dynamic
Copula Model of {CDS} Spreads.}
\newblock \doi{10.1080/07350015.2016.1177535}.
\newblock Forthcoming in Journal of Business \& Economic Statistics.
\bibitem[{Pinheiro and Bates(1996)}]{pinheiro_bates.1996}
Pinheiro JC, Bates DM (1996).
\newblock \enquote{Unconstrained Parametrizations for Variance--Covariance
Matrices.}
\newblock \emph{Statistics and Computing}, \textbf{6}(3), 289--296.
\newblock \doi{10.1007/BF00140873}.
\bibitem[{Pourahmadi and Wang(2015)}]{pourahmadi_wang.2015}
Pourahmadi M, Wang X (2015).
\newblock \enquote{Distribution of Random Correlation Matrices:
{H}yperspherical Parameterization of the {C}holesky Factor.}
\newblock \emph{Statistics \& Probability Letters}, \textbf{106}, 5--12.
\bibitem[{\proglang{R} {Core Team}(2016)}]{R}
\proglang{R} {Core Team} (2016).
\newblock \emph{\proglang{R}: {A} Language and Environment for Statistical
Computing}.
\newblock \proglang{R} Foundation for Statistical Computing, Vienna, Austria.
\newblock \proglang{R}~version 3.2.3,
\urlprefix\url{https://www.R-project.org/}.
\bibitem[{Rapisarda \emph{et~al.}(2007)Rapisarda, Brigo, and
Mercurio}]{rapisarda_etal.2007}
Rapisarda F, Brigo D, Mercurio F (2007).
\newblock \enquote{Parameterizing Correlations: {A} Geometric Interpretation.}
\newblock \emph{IMA Journal of Management Mathematics}, \textbf{18}(1), 55--73.
\newblock \doi{10.1016/j.spl.2015.06.015}.
\bibitem[{Ripley(2015)}]{MASS}
Ripley B (2015).
\newblock \emph{\pkg{MASS}: {S}upport Functions and Datasets for {V}enables and
{R}ipley's {MASS}}.
\newblock \proglang{R}~package version 7.3-45,
\urlprefix\url{https://CRAN.R-project.org/package=MASS}.
\bibitem[{Ryan(2015)}]{quantmod}
Ryan JA (2015).
\newblock \emph{\pkg{quantmod}: {Q}uantitative Financial Modelling Framework}.
\newblock \proglang{R}~package version 0.4-5,
\urlprefix\url{https://CRAN.R-project.org/package=quantmod}.
\bibitem[{Ryan and Ulrich(2015)}]{xts}
Ryan JA, Ulrich JM (2015).
\newblock \emph{\pkg{xts}: {E}xtensible Time Series}.
\newblock \proglang{R}~package version 0.9-7,
\urlprefix\url{https://CRAN.R-project.org/package=xts}.
\bibitem[{Sanderson(2010)}]{sanderson.2010}
Sanderson C (2010).
\newblock \enquote{\pkg{Armadillo}: {A}n Open Source \proglang{C++} Linear
Algebra Library for Fast Prototyping and Computationally Intensive
Experiments.}
\newblock \emph{Technical report}, NICTA.
\newblock \urlprefix\url{http://arma.sourceforge.net/}.
\bibitem[{Shephard(2005)}]{shephard.2005}
Shephard N (2005).
\newblock \emph{Stochastic volatility: selected readings}.
\newblock Oxford University Press, Oxford.
\bibitem[{Venables and Ripley(2002)}]{MASS.2002}
Venables WN, Ripley BD (2002).
\newblock \emph{Modern Applied Statistics with {S}}.
\newblock Fourth edition. Springer--Verlag, New York.
\newblock ISBN 0-387-95457-0,
\urlprefix\url{http://www.stats.ox.ac.uk/pub/MASS4}.
\bibitem[{Ye(1988)}]{ye.1988}
Ye Y (1988).
\newblock \emph{Interior Algorithms for Finear, Quadratic, and Linearly
Constrained Convex Programming}.
\newblock Ph.D. thesis, Stanford University.
\bibitem[{Zhu and Galbraith(2010)}]{zhu_galbraith.2010}
Zhu D, Galbraith JW (2010).
\newblock \enquote{A Generalized Asymmetric {Student-t} Distribution with
Application to Financial Econometrics.}
\newblock \emph{Journal of Econometrics}, \textbf{157}(2), 297--305.
\newblock \doi{10.1016/j.jeconom.2010.01.013}.
\end{thebibliography}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
\clearpage
\appendix
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section[The GAS model with conditional Student-t distribution]{The GAS model with conditional \Student distribution}\label{sec:gas_t}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Let's consider the case when the distribution of the univariate random variable $y_t\in\Re$, conditionally on $\by_{1:t-1}$, is \Student with location $\mu_t$, scale $\phi_t$ and $\nu_t$ degrees of freedom, \emph{i.e.}, $\btheta_t = (\mu_t,\phi_t,\nu_t)^\prime$ and:
\begin{equation}~\label{eq:std}
p(y_t;\btheta_t) \bydef \frac{\Gamma\left(\frac{\nu_t + 1}{2}\right)}{\Gamma\left(\frac{\nu_t}{2}\right)\phi_t\sqrt{\pi\nu_t}}\left(1 + \frac{\left(y_t - \mu_t\right)^2}{\nu_t\phi_t^2}\right)^{-\frac{\nu_t + 1}{2}} \,.
\end{equation}
This model illustrates the benefits of the filter defined in~\eqref{eq:updating} when the data are contaminated by outliers or naturally exhibits extreme values and has been used by \citet{creal_etal.2013} and \citet{lucas_zhang.2016} under the name tGAS, and by \citet{harvey.2013} and \citet{harvey_luati.2014} under the name Beta--t--EGARCH.
Differentiating the logarithm of~\eqref{eq:std} with respect to $\btheta_t$ leads to the score vector $\nabla_t(y_t,\btheta_t)=(\nabla_t^\mu,\nabla_t^\phi,\nabla_t^\nu)'$, with:
%
\begin{align}~\label{eq:score_std}
\begin{split}
\nabla_t^\mu &\bydef \frac{(\nu_t + 1)(y_t - \mu_t)}{\nu_t\phi_t\left(1 + \frac{\left(y_t - \mu_t\right)^2}{\nu_t\phi_t}\right)} \\
\nabla_t^\phi &\bydef \frac{\left(\nu_t + 1\right)\left(y_t - \mu_t\right)^2}{2\nu_t\phi_t^2\left(1 + \frac{\left(y_t - \mu_t\right)^2}{\nu_t\phi_t}\right)} - \frac{1}{\phi_t}\\
\nabla_t^\nu &\bydef \frac{1}{2}\psi\left(\frac{\nu_t + 1}{2}\right) - \frac{1}{2}\psi\left(\frac{\nu_t}{2}\right) - \frac{1}{2\nu_t}\\
&- \frac{1}{2}\log\left(1 + \frac{\left(y_t - \mu_t\right)^2}{\nu_t\phi_t}\right) + \frac{\left(\nu_t + 1\right)\left(y_t - \mu_t\right)^2}{2\nu_t^2\phi_t\left(1 + \frac{\left(y_t - \mu_t\right)^2}{\nu_t\phi_t}\right)} \,,
\end{split}
\end{align}
%
where $\psi(\cdot)$ is the Digamma function. Without loss of generality, let us consider the case when $\gamma = 0$ with no reparametrization, \emph{i.e.}, $\btheta_t = \widetilde\btheta_t$. The results when $\gamma \neq 0$ and a mapping function $\Lambda(\cdot)$ for $\btheta_t$ is introduced are qualitatively the same. Clearly, what controls for the response to extreme observations in the conditional score $\nabla_t(y_t,\btheta_t)$ is the degree of freedom parameter $\nu_t$. Specifically, when $\nu_t$ is small, say $\nu_t = 3$, the conditional distribution of $y_t$ has high probability mass in the tails, which means that extreme observations, which would be considered outliers by the Gaussian yardstick, are likely to be observed.
If we introduce the following mapping function for the unrestricted vector of parameter $\widetilde\btheta_t = (\widetilde\mu_t,\widetilde\phi_t,\widetilde\nu_t)'$:
%
\begin{equation}~\label{eq:mapping}
\Lambda(\widetilde\btheta_t) \bydef \begin{cases}
\mu_t =& \!\!\widetilde\mu_t \\
\phi_t =& \!\!\exp(\widetilde\phi_t) \\
\nu_t =& \!\!\exp(\widetilde\nu_t) + c\,,
\end{cases}
\end{equation}
%
with $c=2$ in order to ensure the existence of $V_{t-1}\left[y_t\right]$, then the GAS updating step for $\btheta_t$ when $\gamma = 0$ takes the form:
%
\begin{align}~\label{eq:update_studentt}
\btheta_{t+1} &= \Lambda(\widetilde\btheta_{t+1})\\
\widetilde\btheta_{t+1} &= \bkappa + \bA\mathcal{J}(\widetilde\btheta_t)'\nabla_t(y_t,\btheta_t) + \bB\widetilde\btheta_t \,,
\end{align}
%
where $\bkappa\bydef(\kappa_\mu, \kappa_\phi, \kappa_\nu)'$, $\bA \bydef \diag(a_\mu, a_\phi, a_\nu)$ and $\bB \bydef \diag(b_\mu, b_\phi, b_\nu)$. In this particular case, the Jacobian matrix $\mathcal{J}(\widetilde\btheta_t)$ takes the form:
%
\begin{equation}
\mathcal{J}(\widetilde\btheta_t) = \begin{pmatrix}
1 & 0 & 0 \\
0 & \exp(\widetilde\phi_t) & 0 \\
0 & 0 & \exp(\widetilde\nu_t)
\end{pmatrix}\,.
\end{equation}
Constraints on the evolution of the GAS parameters can be easily considered fixing the values of the $\bA$ and $\bB$ elements. For example, if the constraint $\nu_t=\nu$ has to be imposed, we set $a_\nu = b_\nu = 0$ during the (log-)likelihood maximization.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section[Mapping functions]{Mapping functions}\label{sec:mapping}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Now we briefly discuss the choice of the mapping function $\Lambda(\cdot)$ for GAS models. We indicate the $i$--th element of $\btheta_t$ and $\widetilde\btheta_t$ as $\theta_{i,t}$ and $\widetilde\theta_{i,t}$, respectively. Analogously, we refer to the $i$--th element of the vector--valued mapping function $\Lambda(\cdot)$ as $\lambda_{i}(\cdot)$, such that $\lambda_i(\widetilde\theta_{i,t})=\theta_{i,t}$.
Generally, there are three types of constraints the time--series analyst wants to impose to $\theta_{i,t}$, namely:
\begin{itemize}
\item[1)] $\theta_{i,t}>c,\quad c\in\Re$
\item[2)] $\theta_{i,t}\in\left(a,b\right)$, for $a,b\in\Re$ and $b>a$
\item[3)] $\theta_{i,t}\in\left(a,b\right)\vert \btheta_t \in\Theta$ for $a,b\in\Re$ and $b>a$~,
\end{itemize}
the additional case when $\theta_{i,t}\in\Re$, and thus $\widetilde\theta_{i,t} = \theta_{i,t}$, implicitly requires that $\lambda_i:\Re\to\Re$ is the identity function.
The first case, $\theta_{i,t}>c$,\footnote{The case $\theta_{i,t}<c$ follows immediately.} covers the situation where, for example, $\theta_{i,t}$ is a scale parameter and, consequently, its positiveness has to be imposed (\emph{i.e.}, $c=0$). In this case, $\lambda_i:\Re\to\left[c,\infty\right)$, and the exponential link function, defined as:
\begin{equation}\label{eq:exponential_link}
\theta_{i,t} = \exp(\widetilde\theta_{i,t}) + c\,,
\end{equation}
can be employed. The second case, $\theta_{i,t}\in\left(a,b\right)$, covers the situation where, for example, $p\left(\cdot;\btheta_t\right)$, is the Asymmetric Student--t Distribution of \cite{zhu_galbraith.2010}, and $\theta_{i,t}$ is its skew parameter defined in $\left(0,1\right)$. In the more general case we have $\lambda_i:\Re\to\left(a, b\right)$, and thus, the modified logistic function:
\begin{equation}\label{eq:logistic_link}
\theta_{i,t} = a + \frac{b - a}{1 + \exp(-\widetilde\theta_t)} \,,
\end{equation}
can be employed. The last case, $\theta_{i,t}\in\left(a,b\right)\vert \btheta_t \in\Theta$, is more complicated and covers the situation where, for example, $p\left(\cdot;\btheta_t\right)$ is a multivariate Gaussian distribution and $\theta_{i,t}$ is one element of its correlation matrix $\bR_{t}$. Clearly, in this case $\theta_{i,t}\in\left[-1,1\right]$, with the equivalence corresponding to the case $N=2$. For the more general case $N>1$, we need to ensure that $\bR_t$ is positive definite, \emph{i.e.}, $\bx^\prime\bR_t\bx>0,\forall\bx\in\Re^N$. Following \cite{creal_etal.2011}, we employ the hyperspherical coordinates transformation originally proposed by \cite{pinheiro_bates.1996} and subsequently discussed in \cite{jaeckel_rebonato.1999}, \cite{rapisarda_etal.2007} and \cite{pourahmadi_wang.2015}. We define the general $(h,k)$--th lower diagonal element of $\bR_t$ as $\rho_{hk,t}=\theta_{i,t}$ for $h>k$, $h<N$ and $\widetilde\rho_{hk,t}=\widetilde\theta_{i,t}$, for $i=1,\dots,N\left(N-1\right)/2$. \cite{pourahmadi_wang.2015} show that:
\begin{equation}\label{eq:hyper_link}
\rho_{hk,t} = c_{h1,t}c_{k1,t} + \sum_{m=2}^{h-1}c_{hm,t}c_{km,t}\prod_{l=1}^{m-1}s_{hl,t}s_{kl,t} + c_{hk,t}\prod_{l=1}^{h-1}s_{hl,t}s_{kl,1}\quad 1\le h<k\le N\,,
\end{equation}
where $c_{hk,t} \bydef \cos\left(\widetilde\rho_{hk,t}\right)$ and $s_{hk,t} \bydef \sin\left(\widetilde\rho_{hk,t}\right)$ for all $1\le h<k\le N$ ensure that $\bR_t \bydef \{\rho_{ij,t}\}_{i,j=1}^N$ is a proper correlation matrix. %Furthermore, $\Lambda\left(\cdot\right)$ is bijective if the additional constraints $\tilde\rho_{hk,t}\in\left[-\pi,\pi\right]$ for all $1\le h<k\le N$ are imposed, however, as discussed in \cite{creal_etal.2011}, these constraints can be removed in practical situations.
These three specifications for $\lambda_i\left(\cdot\right)$ cover all the cases considered in this article and in the \proglang{R} package \pkg{GAS}. Additional information are reported in the \proglang{R} documentation; see \code{help("UniMapParameters")} (\code{help("MultiMapParameters")}) and \code{help("UniUnmapParameters")} (\code{help("MultiUnmapParameters")}), for details on $\Lambda\left(\cdot\right)$ and $\Lambda^{-1}\left(\cdot\right)$ in the univariate (multivariate) case, respectively.
\end{document}
