Extracted main text — title through conclusion, appendix excluded. This is what our citation measures are computed over, published so the extraction can be checked by eye.
92,474 characters · 16 sections · 124 citation commands
Stochastic Variational Inference for GARCH Models
\email{[email removed]}
\email{[email removed]} \equalcont{These authors contributed equally.}
\email{[email removed]} \equalcont{These authors contributed equally.}
\email{[email removed]} \equalcont{These authors contributed equally.}
\affil*[1]{\orgdiv{School of Mathematics and Statistics}, \orgname{University of New South Wales}, \orgaddress{\street{Anita B. Lawrence Centre}, \city{Kensington}, \postcode{2052}, \state{New South Wales}, \country{Australia}}}
\affil[2]{\orgdiv{Research School of Finance, Actuarial Studies and Statistics}, \orgname{The Australian National University}, \orgaddress{\street{Building 26C, Kingsley Street}, \city{Canberra}, \postcode{2601}, \state{ACT}, \country{Australia}}}
\affil[3]{\orgdiv{UNSW Data Science Hub}, \orgname{University of New South Wales}, \orgaddress{\street{Anita B. Lawrence Centre}, \city{Kensington}, \postcode{2052}, \state{NSW}, \country{Australia}}}
\affil[4]{\orgdiv{School of Mathematics and Statistics}, \orgname{University of Sydney}, \orgaddress{\street{Carslaw Building}, \city{Camperdown}, \postcode{2006}, \state{NSW}, \country{Australia}}}
\affil[5]{\orgdiv{ARC Training Centre in Data Analytics for Resources and Environments}, \orgname{University of Sydney}, \orgaddress{\street{ Biomedical Building}, \city{South Eveleigh}, \postcode{2015}, \state{NSW}, \country{Australia}}}
\abstract{Stochastic variational inference algorithms are derived for fitting various heteroskedastic time series models. We examine Gaussian, t, and skew-t response GARCH models and fit these using Gaussian variational approximating densities. We implement efficient stochastic gradient ascent procedures based on the use of control variates or the reparameterization trick and demonstrate that the proposed implementations provide a fast and accurate alternative to Markov chain Monte Carlo sampling. Additionally, we present sequential updating versions of our variational algorithms, which are suitable for efficient portfolio construction and dynamic asset allocation. }
Most financial decisions are based on the trade-off between risk and return (10.2307/2975974, 10.2307/2975974). Modeling volatility, which is often used to measure risk, is therefore crucial in risk management. Traditional econometric models assume that the variance of asset returns is constant, but this assumption may be unrealistic as it overlooks the presence of heteroskedasticity in financial time series (engle2001garch, engle2001garch). To address this issue, engle1982autoregressive introduced auto-regressive conditional heteroskedastic (ARCH) models, which utilize lagged disturbance terms to explain a time-varying conditional variance process. ARCH models were extended to the generalized ARCH (GARCH) class by bollerslev1986generalized to allow for a more flexible lag structure where the current conditional variance depends not only on the past conditional variances but also on the past disturbance terms. Currently, GARCH models are widely employed to describe the dynamics of volatility processes; however, conventional ARCH and GARCH models assume Gaussian errors, which fail to capture common features of financial time series such as heavy-tailedness and skewness. To improve both goodness-of-fit and prediction accuracy, bollerslev1987conditionally and lambert2001modelling respectively studied the use of t and skew-t distributions for modeling error terms (see also degiannakis2004volatility, degiannakis2004volatility).
Maximum likelihood estimation is typically used to estimate the parameters of GARCH models gonzalez1999efficiency, as it guarantees ease of implementation and asymptotic efficiency silvennoinen2021consistency. However, uncertainty quantification of maximum likelihood estimates is usually performed by taking advantage of their asymptotic normality nakatsuma2000bayesian, which may lead to underestimating uncertainty and model risk as a consequence. Bayesian inference offers an attractive alternative as it naturally expresses the uncertainty around the estimation process by returning a posterior probability distribution for the model parameters. In a Bayesian context, nakatsuma2000bayesian proposed a Markov chain Monte Carlo (MCMC) method with Metropolis-Hasting steps for fitting GARCH models. ardia2008financial used MCMC to fit regime-switching GARCH models with t innovations. A fully automated MCMC procedure for Bayesian estimation in GARCH models with t innovations is available in the R package bayesGARCH ardia2010bayesian. li2021efficient propose a sequential Monte Carlo algorithm for Student t GARCH models. While several MCMC methods are available for GARCH models gerlach2006mcmc, henneke2011mcmc, virbickaite2015bayesian, contino2017bayesian, xaba2021performance, these may suffer from slow convergence and be computationally intensive, in particular for long time series shiferaw2019time, mozumder2021option,livingston2022bayesian.
Variational Bayes (VB) methods jordan1998introduction have emerged as a more computationally efficient alternative to MCMC. Early research on variational inference focused on deriving algorithms with closed-form updates for conjugate models -- see, for example, bishop2006pattern, ormerod2010explaining or maestrini2018variational for variational inference in skew-t regression through conjugate model representations. Over recent years, stochastic gradient methods for variational inference (e.g., paisley2012variational, paisley2012variational) have found wider interest, especially for non-conjugate models and as an alternative to avoid lengthy mathematical derivations.
The performance of stochastic optimization algorithms, including stochastic VB, greatly depends on the variance of gradients and this has led to the development of computational solutions to reduce it. ranganath2014black introduced a technique based on control variates to improve the convergence of stochastic VB. This method requires minimal analytic calculations and can be applied to a variety of models. kingma2013auto proposed an alternative variance reduction technique named reparametrization trick, which consists in rewriting an expectation so that the distribution with respect to which the gradient is taken is independent of the model parameters. Another way of efficiently using the gradient of the model joint probability density to increase computational speed is presented in titsias2014doubly and other approaches are discussed in tran2021practical.
Variational inference has found limited attention in the literature on GARCH models. This is partially due to the characteristics of time series, which represent an obstacle to scalable algorithm implementations. hoffman2013stochastic developed a scalable stochastic variational inference algorithm to analyze massive datasets for complex topic models by subsampling the data to form noisy estimates. However, there is no clear and unbiased way to subsample time series data for GARCH models politis2003impact, alonso2006comparison. tran2021variational proposed a manifold-based variational Bayes approach for GARCH models with Gaussian error terms that maps the variational parameters from the Euclidean space to a Riemannian manifold, but its extension to models with non-Gaussian errors is not immediate.
In recent years there have been several attempts to develop sequential variational inference approaches for time series models, i.e., implementations of variational algorithms that update inference as soon as new data points of the series are available. lambert2022recursive proposed recursive Gaussian variational approximations for models assuming observations are independent, which are therefore not suitable for GARCH models and time-dependent data. tomasetti2022updating developed an updating VB (UVB) method as a sequential extension of batch stochastic variational approximations for streaming data broderick2013streaming, foti2014stochastic. This method is based on deriving a sequence of posterior distributions that are used as priors for the following steps through updates resulting from Bayes' theorem. Each posterior update only requires data observed in the previous step. While very efficient, UVB targets the approximating pseudo-posterior rather than the true posterior, resulting in reduced accuracy. An alternative sequential implementation is the one proposed by gunawan2021variational, which targets the true posterior distribution instead. Their approach is named sequential stochastic VB (Seq-SVB) and it is used to perform prediction with factor stochastic volatility models. A potential drawback of Seq-SVB is that the likelihood function is computed for all the time series, affecting computational efficiency.
This work includes two main contributions. First, we develop stochastic variational Bayes (SVB) algorithms to approximate the posterior distributions of parameters in GARCH models with Gaussian, t and skew t errors (also known as innovations in the econometric literature). Our variational approximations are compared to an MCMC target. We assess the performance of variational inference using control variates or the reparametrization trick in terms of computational speed and estimation accuracy. Results demonstrate that our stochastic variational implementations are accurate while being faster than MCMC. Our second contribution is the adaptation of the stochastic VB approaches to allow for sequential updates. We implement the procedures developed by gunawan2021variational and tomasetti2022updating and compare them for the first time, in the specific case of GARCH models. Our results demonstrate that both UVB and Seq-SVB are faster than the base SVB and show the presence of a trade-off between accuracy and computational costs.
The remainder of the paper is organized as follows. Section (ref) reviews GARCH models, variational Bayes approximations and their stochastic versions. Section (ref) develops VB algorithms for GARCH models with different innovation terms. Section (ref) proposes sequential versions of the variational algorithms from the previous section for efficiently updating posterior distributions when new observations are available. Section (ref) presents the results of our simulation experiments and Section (ref) illustrates the use of VB for real-world financial time series. Section (ref) concludes the paper with a discussion.
Let $y_1,\ldots,y_T$ be a time series, for instance the log-difference of stock prices at two consecutive time points, i.e., $y_t = \log P_t - \log P_{t-1}$. The GARCH$(p,q)$ model for the $t$th observation, $t=1,\ldots,T$, takes the form
The time-varying volatility process for $\sigma_t^2$ follows a deterministic autoregressive structure and here we implicitly assume the dependence of this structure on the $\alpha_i$ and $\beta_j$ parameters writing $\sigma_t^2$ instead of $\sigma_t^2(\alpha_1,\ldots,\alpha_p,\beta_1,\ldots,\beta_q)$. The success of this class of models lies in the fact that the log-likelihood function $\ell(\theta) = \log L(\theta;y)$ of $\theta = (\omega,\alpha_1, \ldots, \alpha_p, \beta_1, \ldots, \beta_q)$ is tractable and has the following expression for model (ref):
In real applications it has been widely recognized that financial returns exhibit significant deviations from normality, and typically have heavier tails than the normal distribution shirvani2020stock. To address this issue, GARCH models have been extended to include t error terms bollerslev1987conditionally, baillie1989common, beine2002central. This is achieved by replacing the normal distribution assumption in (ref) with $\varepsilon_t \sim \mathcal{T}(0,1;\nu)$, where $\mathcal{T}(0,1;\nu)$ represents a standardized t distribution with degrees of freedom $\nu$. A random variable $x$ has a standardized t distribution $\mathcal{T}(0,1;\nu)$ with degrees of freedom $\nu>2$ if its probability density function is $$f(x)=\frac{\Gamma\Big(\frac{\nu+1}{2}\Big)}{\sqrt{\pi(\nu-2)}\Gamma\Big(\frac{\nu}{2}\Big)}\left(1+\frac{x^2}{\nu-2}\right)^{-\frac{\nu+1}{2}}.$$ The log-likelihood function for the t GARCH$(p,q)$ model defined by the parameter vector $\theta = (\nu,\omega,\alpha_1, \ldots, \alpha_p, \beta_1, \ldots, \beta_q)$ is
Skewness is a pervasive feature observed in empirical data which can be taken into account in financial time series modeling and inference jondeau2019average. lambert2001modelling incorporated the skew t distribution introduced by fernandez1998bayesian in GARCH models to account for skewness. Skew-t distributions are highly flexible and introduce asymmetry in the data model by adding only one parameter to the t distribution. In the corresponding model, the normal distribution assumption in (ref) is replaced by $\varepsilon_t \sim \mathcal{ST}(0,1;\nu,\xi)$, where the skewness parameter $\xi$ is a positive real number. A value of $\xi<1$ indicates negative skewness, while a value of $\xi>1$ indicates positive skewness. We follow the framework and definition provided by fernandez1998bayesian: a random variable $x$ has a skewed t distribution $\mathcal{ST}(0,1;\nu,\xi)$ with degrees of freedom $\nu>2$ and skewness parameter $\xi>0$ if its probability density function is
where
Note that the density function has been standardized to have a zero mean and unit variance. Defining the vector of parameters as $\theta = (\nu,\xi,\omega,\alpha_1, \ldots, \alpha_p, \beta_1, \ldots, \beta_q)$, the log-likelihood for the skewed t GARCH model with $T$ observations is given by
While the likelihood is still tractable, inference based on MCMC may be computationally costly in this case, and convergence slow iqbal2021bayesian.
Let $y$ denote a vector of observed data, and let $\theta$ be a vector of parameters. In a Bayesian setting, given a prior distribution $p(\theta)$ and a likelihood function $L(\theta; y)$, the prior beliefs about $\theta$ are updated after collecting the data to obtain the posterior distribution:
where $p({y})$ is the marginal likelihood, or evidence. In practice, the posterior distribution in Bayesian inference may be intractable, as its normalizing factor, or evidence, $p({y})$, is often not available in closed form. Consequently, simulation-based Markov chain Monte Carlo (MCMC) methods are typically employed to estimate the posterior distribution through sampling robert1999monte. However, MCMC sampling can be computationally expensive or impractical for complex models. An alternative approach is provided by variational Bayes (VB), which aims to approximate the intractable posterior distribution using a variational approximating density $q(\theta)$ selected from a family of densities $\mathcal{Q}$. In VB, the approximating density is chosen to minimize the Kullback-Leibler (KL) divergence between the variational approximating density and the true posterior density.
where $\lambda$ denotes a vector of parameters of the variational approximating density in the parameter space $\mathcal{M}$. Minimization of the KL divergence is equivalent to maximizing the evidence lower bound (ELBO) $\mathcal{L}(\lambda)$ (see, e.g., bishop2006pattern, bishop2006pattern or ranganath2014black, ranganath2014black), which is defined as
In general, $\mathcal{L}(\lambda)$ does not have a closed-form solution in non-conjugate models. Stochastic variational inference approaches based on stochastic gradient methods (paisley2012variational, paisley2012variational; nott2012regression, nott2012regression; hoffman2013stochastic, hoffman2013stochastic; kingma2013auto, kingma2013auto; ranganath2014black, ranganath2014black; rezende2014stochastic, rezende2014stochastic) allow for efficient optimization of $\mathcal{L}(\lambda)$, even when the ELBO cannot be evaluated analytically. In stochastic variational inference an initial value of the variational parameters $\lambda^{(0)}$ is stated and iteratively updated through
until a stopping rule is satisfied. Here $\widehat{\nabla_\lambda \mathcal{L}(\lambda^{(i)})}$ denotes an unbiased estimator of the gradient vector of $\mathcal{L}(\lambda)$ computed with respect to $\lambda$ and $\{a_i,i \geq 1\}$ is a sequence of vector-valued learning rates satisfying the Robbins-Monro conditions $\sum_i a_i = \infty$ and $\sum_i a^2_{i} < \infty$ robbins1951stochastic. The topic of adaptive learning rates has been widely discussed in the literature (e.g., polyak1964some, polyak1964some; duchi2011adaptive, duchi2011adaptive; tieleman2012lecture, tieleman2012lecture; zeiler2012adadelta, zeiler2012adadelta). In this work, we adopt ADAM adaptive learning rates kingma2014adam and a moving average stopping criterion (see tran2021practical, tran2021practical for more details).
It is possible to show that the gradient vector can be written in an expectation form computed with respect to the variational density $q_\lambda(\theta)$ as
where $h_{\lambda}(\theta)= \log\big\{ p(\theta)L(\theta;y)\big/q_{\lambda}(\theta)\big\}$ (e.g., see tran2021practical, tran2021practical). A simple estimator of the $j$-th element of the ELBO gradient is then given by
where $s=1,\ldots,S$ are samples from the variational distribution $q_\lambda(\theta)$.
In this work, we consider a popular form of VB known as Gaussian variational approximation, which assumes $q_\lambda(\theta)$ to be a multivariate Gaussian distribution with mean $\mu$ and covariance matrix $\Sigma$. Several authors have studied and discussed the parameterization of Gaussian variational approximations; see, for example, opper2009variational, opper2009variational; ormerod2012gaussian, ormerod2012gaussian; ong2018gaussian, ong2018gaussian; quiroz2022gaussian, quiroz2022gaussian. We adopt the Cholesky factorization suggested by titsias2014doubly and factorize the covariance matrix as $\Sigma = LL^\top$ where $L$ is a lower triangular matrix. Consequently, the true posterior distribution $p(\theta\vert y)$ can be approximated by a $\mathcal{N}_{d}(\mu, LL^\top)$ distribution; in this case $\lambda=(\mu^\top,\mbox{vech}(L)^\top)^\top$, where $\mbox{vech}(\cdot)$ is the half-vectorization operator. Introducing this factorization reduces the computational cost to reach convergence while preserving the positive definiteness of the covariance matrix $\Sigma$.
The performance of stochastic variational inference greatly depends on the variance of the gradient estimator liu2019variance. paisley2012variational and ranganath2014black adopt a variance reduction technique based on the use of control variates, where the gradient estimator function is replaced by another function with equivalent expectation but smaller variance. For any arbitrary number $c_j$, the $j$-th element of the estimator version of (ref) can be slightly modified to
which is still an unbiased estimator of the gradient of $\mathcal{L}(\lambda)$ (ranganath2014black, ranganath2014black; tran2021practical, tran2021practical). The optimal control variate $c_j$ that minimizes the variance can be calculated by taking the derivatives of the variance of the gradient estimator with respect to $c_j$:
where $j = 1, \ldots, d$.
The reparametrization trick (kingma2014adam, kingma2014adam; rezende2014stochastic, rezende2014stochastic; titsias2014doubly, titsias2014doubly) offers an attractive alternative to the control variates approach. Suppose samples are meant to be generated from $\theta \sim q_\lambda(\theta)$, then it is often possible to write $\theta=g_\lambda(\epsilon)$, where $\epsilon$ is an auxiliary variable with an independent marginal distribution $p(\epsilon)$ and $g_\lambda(\epsilon)$ is a deterministic function parameterized by $\lambda$. Using this deterministic transformation, it is possible to rewrite $\mathcal{L} (\lambda)$ in (ref) as an expectation with respect to the alternative measure $p(\epsilon)$, i.e. as
Then the gradient of the evidence lower bound $\nabla_\lambda\mathcal{L}(\lambda)$ can be calculated with respect to the measure $p(\epsilon)$
(see tran2021practical, tran2021practical for details). The reparametrization trick is used to express the gradient of the expectation as an expectation of a gradient. Provided that $g_{\lambda}$ is differentiable, then samples can be generated from $p(\varepsilon)$ to estimate $\nabla_{\lambda} \mathcal{L}(\lambda)$ via Monte Carlo.
When using the reparametrization trick and Gaussian variational approximations with a Cholesky factorization, the estimates of the $\mathcal{L}(\lambda)$ gradient depend on the gradient of the deterministic transformation function $\nabla_\lambda g_\lambda(\epsilon)$, as shown in (ref), and are computed with respect to the variational parameters $\mu$ and $L$. The gradient of $\mathcal{L}(\lambda)$ with respect to the mean vector $\mu$ can be estimated via Monte Carlo sampling from the alternative probability measure $p(\epsilon)$ by noting that
Let $\operatorname{vech}(A)$ denote the half-vectorization of a matrix $A$ of size $d\times d$, i.e. the vector of length $d(d+1)/2$ obtained by stacking the elements in the matrix lower triangle, including the elements along the diagonal of $A$, from left to right and top to bottom. The gradient of the half-vectorization of $L$ under $p(\epsilon)$ is
with $\theta=\mu+L\epsilon$.
In this section, we derive an algorithm for fitting GARCH models via SVB (Stochastic Variational Bayes). The algorithm provides approximate posterior densities for the parameters of interest in models with the distributional assumptions on the innovations of Section (ref).
We adopt the notation of wand2014fully and ong2018variational and start by deriving the algorithm using the control variates approach first, followed by the reparametrization trick.
To define an SVB (Stochastic Variational Bayes) algorithm, it is necessary to have an unbiased estimator of the gradient vector of $\mathcal{L}(\lambda)$ with respect to the variational parameter $\lambda$, denoted as $\widehat{\nabla_\lambda \mathcal{L}(\lambda)}$.
Considering Equation (ref), this estimate depends on:
where $\ell(\theta)$ is the log-likelihood function, $p(\theta)$ is the prior distribution and $q_{\lambda}(\theta)$ is the variational distribution. It is easy to notice that $\log q_{\lambda}(\theta)$ and, consequently, $\nabla_\lambda \log q_{\lambda}(\theta)$ depend only on the variational density $q_\lambda(\theta)$, and they do not rely on the specific model under examination. Here, we take $q_{\lambda}(\theta)$ to be a multivariate Gaussian distribution. The simplest approximation assumed the covariance matrix $\Sigma$ in the Gaussian variational density $q_\lambda(\theta)$ has a fully independent structure, meaning it is diagonal (mean-field approximation). However, this approximation lacks the ability to capture the dependence structure of the target posterior density and tends to underestimate the uncertainties of the model parameters ong2018gaussian. We will compare results using a dense or a diagonal covariance matrix.
According to tan2018gaussian, faster convergence in variational approximations for state-space models can be achieved by parametrizing the precision matrix instead of the covariance matrix. In our simulations, we have observed that, for GARCH models, variational algorithms using control variates converge faster when the approximation is based on a parameterization of the precision matrix. However, the opposite holds true when using the reparametrization trick. We define the precision matrix as $\Omega = \Sigma^{-1} = CC^\top$, where $C$ is a lower triangular matrix.
To implement an algorithm based on control variates, the variational parameters to be optimized are $\lambda=(\mu^\top,\operatorname{vech}(C)^\top)^\top$, where $\mu$ is a $d$-dimensional vector, and $\operatorname{vech}(C)$ is the half-vectorization of $C$, i.e. a vector of $d(d+1)/2$ elements obtained by vectorizing the lower triangular part of $C$. The log-density of the multivariate Gaussian distribution is
and the gradient of the log-multivariate Gaussian density with respect to $\mu$ and $\operatorname{vech}(C)$ is $$\nabla_\lambda \log q_\lambda(\theta) = (\nabla_\mu \log q_\lambda(\theta)^\top,\nabla_{\operatorname{vech}{(C)}} \log q_\lambda(\theta)^\top)^\top,$$ where
Here $\operatorname{diag}(1/C)$ denotes a diagonal matrix whose $i$-th entry is $1/C_{ii}$ and all the elements of $\mu$ and $\operatorname{vech}{(C)}$ are not restricted. To ensure that $C$ is positive-definite, we parameterize the optimization problem in terms of the lower triangular matrix $C'$, where $C'_{ij}=C_{ij}$ if $i\neq j$ and $C'_{ii}=\exp(C_{ii})$.
While producing an unbiased estimator of the ELBO with smaller variance, the use of control variates can slow down the variational algorithm implementation tran2021practical and the reparametrization trick can be used as an alternative to derive faster algorithms. This approach uses a reparametrized gradient estimator that is adapted to the specific model. From (ref), the generic expression for the gradient is
Finally, the gradient of the log multivariate Gaussian variational density with respect to $\theta$ is
For simplicity, we will restrict our attention to the GARCH(1,1) model. Extensions to the general GARCH(p,q) as given in (ref) are straightforward. The parameters of a GARCH(1,1) with Gaussian innovations are $\phi = (\omega,\alpha,\beta)^\top$, with $\omega > 0$, $\alpha>0$, $\beta>0$ and $\alpha+\beta <1$. tran2021variational suggest the following transformations to work with unconstrained parameters: $\alpha = \psi_1\psi_2$, $\beta = \psi_1(1-\psi_2)$ with $0<\psi_i<1$, $i=1,2$, $\theta_{\omega} =\log(\omega)$, $\theta_{\psi_1} =\log\{\psi_1/(1-\psi_1)\}$ and $\theta_{\psi_2} =\log\{\psi_2/(1-\psi_2)\}$, therefore $\theta_{\phi} = (\theta_\omega,\theta_{\psi_1},\theta_{\psi_2})^\top$, and $\theta=\theta_{\phi}$.
An inverse gamma prior distribution $IG(1,1)$ for $\omega$ and independent uniform prior distributions $U(0,1)$ for $\psi_1$ and $\psi_2$ are assumed. These prior distributions induces a negative log-gamma density $NLG(1,1)$ for $\theta_{\omega}$ and an independent standard logistic distribution for $\theta_{\psi_1}$ and $\theta_{\psi_2}$, i.e.
and
The log-likelihood function of the GARCH(1,1) model with Gaussian innovations is defined in (ref) and using (ref) this gives rise to
Using (ref), it is then possible to implement the control variates approach.
To implement a variational algorithm using the reparametrization trick, the gradients $\nabla_\theta \ell(\theta)$ and $\nabla_\theta \log p(\theta)= \left(\nabla_{\theta_\omega} \log p(\theta_\omega),\nabla_{\theta_{\psi_1}} \log p(\theta_{\psi_1}),\nabla_{\theta_{\psi_2}} \log p(\theta_{\psi_2})\right)^\top$ for the GARCH(1,1) model with Gaussian innovations are needed. The three components of the first gradient are
and
The gradient of the log-likelihood function for the Gaussian GARCH(1,1) model (see also fiorentini1996analytic, fiorentini1996analytic) in (ref) with respect to the unconstrained parameters is
and the partial derivatives of $\sigma_t^2$ with respect to the unconstrained parameters are given in Appendix (ref). Finally, $\nabla_\theta{h_{\lambda}(\theta)}$ defined in (ref) takes the following form:
where $\nabla_\theta \log p(\theta_{\phi})$ is stated in (ref)--(ref).
To adapt SVB with control variates to the GARCH model with t innovations, it is necessary to derive the model-specific $\log p(\theta)$ and $\ell(\theta)$. The log-likelihood function for the t model is given in Equation (ref). The prior distributions for $\omega$, $\alpha$, and $\beta$ in the time-varying variance equation for $\sigma_t^2$ remain the same as in the Gaussian case. As for the degrees of freedom parameter $\nu$, we choose a translated exponential prior density ardia2008financial.
Similarly to $\omega$, $\alpha$, and $\beta$, $\nu$ can also be mapped to an unrestricted space by setting $\nu=\log({\exp(\theta_\nu)+1})+2$. According to tan2018gaussian, this specific transformation ensures greater numerical stability compared to the transformation $\nu=\exp(\theta_\nu)+2$. The resulting log-prior density of $\theta_{\nu}$ is then
The parameters to be optimized are then the unconstrained ones, $\theta = (\theta_{\phi},\theta_{\nu})$, where $\theta_{\phi}=(\theta_{\omega},\theta_{\psi_1},\theta_{\psi_2})^T$, and $\phi=(\omega,\alpha,\beta)^T$.
The expression of $h_{\lambda}(\theta)$ defined in (ref) takes the following form for the t GARCH(1,1) model:
where $\log p(\theta)$ is given by the sum of (ref)--(ref) and (ref).
To implement SVB with the reparametrization trick, the gradients of the log-prior densities for the $\theta_{\phi}$ components remain the same as in Equations (ref)--(ref). Additionally, the gradient of the log-prior density for $\theta_\nu$ is given by
The derivational steps for the gradient of the log-likelihood function from the t GARCH model are provided in Appendix (ref) (also refer to levy2003analytic, levy2003analytic). Therefore, the expression for $\nabla_\theta{h_{\lambda}(\theta)}$, as defined in Equation (ref) for the t GARCH(1,1) model, takes the following form:
where the expression of $\nabla_\theta \ell(\theta)=\left(\frac{\partial \ell(\theta)}{\partial \theta_{\phi}},\frac{\partial \ell(\theta)}{\partial \theta_\nu}\right)^\top$ can be found in the appendix (ref) and $\nabla_\theta \log p(\theta)$ is stated in (ref)--(ref) and (ref).
The log-likelihood for the GARCH(1,1) model with standardized skewed t innovations is presented in Equation (ref). Additionally, the skewness parameter $\xi$ is assumed to follow an inverse gamma $IG(1,1)$ prior distribution, according to the same reasoning in the definition of the prior distribution for $\omega$:
Once again, the skewness parameter $\xi$ can be transformed into an unconstrained parameter using the transformation $\xi = \log{\exp(\theta_\xi)+1}$. As a result, the log-prior density of $\theta_\xi$ is given by
The original set of parameters $(\phi,\nu,\xi)^\top$ with $\phi = (\omega,\alpha,\beta)^\top$ is replaced by the unrestricted parameters $\theta=(\theta_{\phi},\theta_{\nu},\theta_{\xi})^\top$, where $\theta_{\phi} = (\theta_\omega,\theta_{\psi_1},\theta_{\psi_2})$. The prior densities for the constrained parameters $\omega$, $\alpha$, $\beta$, and $\nu$ remain the same as those of the t GARCH model described in Section (ref).
The expression for $h_{\lambda}(\theta)$ defined in Equation (ref) takes the following form for the skewed t GARCH(1,1) model:
where $\log p(\theta)$ is given by the sum of (ref)--(ref), (ref) and (ref).
To implement SVB with the reparametrization trick, the gradients of the log-prior densities with respect to $\theta_\omega, \theta_\alpha, \theta_\beta, \theta_\nu$ are required. These gradients can be found in Appendix (ref). Additionally, the gradient of the log-prior density for $\theta_\xi$ is
The gradient for the log-likelihood function of the skewed t GARCH model is given in Appendix (ref). From (ref),
where $\nabla_\theta \log p(y\vert \theta)$ is provided in Appendix (ref) and $\nabla_\theta \log p(\theta)$ is given by (ref)--(ref), (ref) and (ref).
Algorithms (ref) and (ref) outline SVB procedures for GARCH-type models with different types of innovations, namely Gaussian, t, and skewed t. These algorithms utilize the Gaussian Cholesky factorization and incorporate the use of control variates and the reparametrization trick, respectively.
In the previous section, we introduced SVB approximations for GARCH models. One characteristic of time series data, such as those described by GARCH models, is that data points may become available sequentially. Ideally, inference should be updated sequentially as new information becomes available, avoiding the inefficient re-evaluation of the likelihood function for the entire dataset. In this section, we present implementations of the variational approximations from the previous section that allow for sequential updating. We adapt the strategies of gunawan2021variational and tomasetti2022updating for the sequential updating of variational approximations for GARCH models for the first time in the literature and compare their performance.
In their work, tomasetti2022updating propose a procedure called updating variational Bayes (UVB), which enables sequential updates of posterior distributions without the need to re-compute the likelihood function at each step.
Consider a sequence of time points $T_1, T_2, \ldots$, and suppose we want to update the posterior density $p(\theta \vert y_{1:T_n})$ after observing new data points in the interval $[T_{n+1},T_{n+h}]$ to derive $p(\theta \vert y_{1:T_{n+h}})$. The posterior distribution computed up to time $T_n$ can then be used as the prior distribution for the subsequent observation interval. As a result, the likelihood function for the interval $[T_{n+1},T_{n+h}]$ is based solely on $h$ data points:
Hence, the posteriors for new observations can be updated sequentially by recursively applying Bayes' theorem.
We now define the UVB procedure for GARCH-type models. To initialize the UVB algorithm, we can employ the batch stochastic variational approach described in Section (ref) to compute the first variational approximation $q_{\lambda_n}(\theta\vert y_{1:T_n})$ of the posterior density $p(\theta\vert y_{1:T_n})$. After obtaining the initial variational approximation and observing new data points in the interval $[T_{n+1},T_{n+h}]$, we can approximate the posterior distribution defined in Equation (ref) with a pseudo-posterior obtained using the variational approximation of the posterior distribution after the first $T_n$ points:
Instead of minimizing the Kullback-Leibler divergence between the variational density $q_{\lambda_{n+h}}(\theta\vert y_{1:T_{n+h}})$ and the exact posterior density $p(\theta\vert y_{1:T_{n+h}})$, UVB modifies the objective function by substituting the exact posterior $p(\theta\vert y_{1:T_{n+h}})$ with the pseudo-posterior density $\Tilde{p}(\theta\vert y_{1:T_{n+h}})$. For each $n=1,2,\ldots$, the variational parameter $\lambda_{n+h}$ is optimized by maximizing the evidence lower bound (ELBO) using a batch stochastic variational approximation. This modification enables UVB to update the posterior density sequentially and efficiently. It is worth noting that in this article, we employ the faster reparametrization trick approach for UVB, whereas the original paper only considers the control variates case, which was found to be slower than the reparametrization trick approach in our numerical experiments. Algorithm (ref) outlines the steps involved in implementing the UVB approach.
The alternative approach proposed by gunawan2021variational, called sequential SVB (Seq-SVB), differs from UVB in that it targets the true posterior distribution instead of the pseudo-posterior distribution (ref).
Instead of sequentially using variational approximations to the posterior distribution, Seq-SVB optimizes the variational parameters $\lambda_n$ at time $T_n$ and employs them as initial values for the stochastic variational approximations of the variational parameter $\lambda_{n+1}$ at time $T_{n+1}$. Seq-SVB achieves greater accuracy than UVB as it directly targets the true posterior distribution. However, it comes at a higher computational cost since it requires evaluating the likelihood function at each time in the interval $[T_1, T_{n+1}]$.
Now, we define the Seq-SVB algorithm for GARCH models. In the initialization step, batch SVB, as outlined in Algorithm (ref) or (ref), is utilized for the data up to time $T_n$ to obtain the optimal variational parameters $\lambda$ up to time $T_n$. Subsequently, the first variational approximation $q_{\lambda_n}(\theta \vert y_{T_1:T_n})$ of the true posterior $p(\theta \vert y_{T_1:T_n})$ is obtained. After observing new data points up to time $T_{n+1}$, another SVB algorithm can be implemented using all the available data in the interval $[T_1,T_{n+1}]$, with the optimal $\lambda_n$ from the initialization serving as the initial value for the variational parameter. This procedure is repeated for the instances $(n+1), (n+2), \ldots$. The implementation of Seq-SVB is summarized in Algorithm (ref).
We first assess the accuracy of the batch SVB approximation discussed in Section (ref) by applying Algorithms (ref) and (ref) to GARCH models with the innovation types under examination. All computations were performed in R r2023 on the Katana computational cluster with 8-CPU cores and 4GB memory per core of the University of New South Wales (\url{https://docs.restech.unsw.edu.au/using_katana/about_katana/}).
We generated time series of lengths 1,000 and 5,000 from GARCH(1,1) models with Gaussian, t, and skew t innovations using the R package fGarch wuertz2022package. For each setting, we simulated 1,000 replicates. The GARCH models were generated using parameter values similar to those estimated from the S&P500 return index time series analyzed in Section (ref) with MCMC. Specifically, we used $\omega=0.1$, $\alpha=0.2$, and $\beta=0.75$ for Gaussian, t, and skew t innovations. We set $\nu=4$ for the t and skew t degrees of freedom, and $\xi=0.8$ as the skew t asymmetry parameter.
The algorithm parameters of Algorithms (ref) and (ref) were set to the default values employed in tran2021practical, with some slight adjustments made to improve convergence. Specifically, the maximum patience parameter was set to $P=100$, and the moving average window size was set to $t_W=25$. Both adaptive learning weights were set to $\beta_1=\beta_2=0.9$, and the adaptive learning rate was set to $\eta_0=0.02$. We used $\tau=1,000$ as the threshold. In order to investigate the impact of the number of Monte Carlo samples on the efficiency and computational time, we conducted experiments using different values of $S$ for both the control variates and reparametrization methods. Specifically, we fixed the number of MC samples at $S=5, 10,$ and $50$ for both methods. By comparing the results obtained with these different values of $S$, we can assess the influence of the number of MC samples on the estimation accuracy and computational efficiency of the two approaches. In this paper, we focus on presenting the results obtained using the control variates method with $S=10$, and the reparametrization trick method with $S=5$. These particular choices for the number of MC samples were selected based on similar computational time, ensuring a fair comparison between the two approaches (See table (ref)). The additional results obtained using different values of $S$ can be found in the appendix (ref), providing a comprehensive analysis of the sensitivity of the results to the choice of MC sample size.
Figure (ref) shows the ELBO values with respect to the number of iterations in implementations based on control variates and the reparametrization trick, using both the Cholesky factorization of the precision matrix $\Omega$ and the covariance matrix $\Sigma$. Interestingly, convergence appears to be faster when using the precision matrix with control variates, and when using the covariance matrix with the reparametrization trick. In particular, when using the covariance matrix with control variates, SVB may fail to reach convergence in some cases.
The accuracy of the approximation is analyzed by comparing it to the posterior distribution obtained through MCMC with 1,000,000 iterations, after discarding the first 100,000 samples as burn-in. MCMC was implemented using random walk Metropolis-Hastings with Gaussian proposal distributions. The mean and variance of the proposal distributions were set to the maximum likelihood estimates, following chan2016modeling. For comparison with the MCMC target, we computed the variational approximation accuracy with respect to the posterior density as
with 100% indicating perfect matching between the variational density $q_{\lambda}(\theta)$ and the posterior distribution $p(\theta \vert y)$ obtained via MCMC.
Figure (ref) compares the accuracy of the variational approximations obtained when using SVB (on the left) and mean-field SVB (on the right) with control variates for time series of length 1,000 from a GARCH(1,1) model with Gaussian innovations. As expected, the mean-field SVB approach, which uses a fully diagonal covariance matrix for the variational approximating density, results in significantly lower accuracy for all the parameters. On the other hand, the implementation based on a full covariance matrix allows to achieve high accuracy compared to the MCMC target.
Figure (ref) presents a comparison of the approximated posterior distributions obtained using SVB with a diagonal (mean-field) and non-diagonal covariance matrix, with control variates, or with the reparametrization trick and a non-diagonal covariance matrix, along with the MCMC posterior. On the right, the ELBO is illustrated as a function of the number of iterations. The results are based on time series of length 1,000 with Gaussian innovations. Once again, it is evident that the mean-field SVB approach tends to underestimate the posterior uncertainty. On the other hand, both control variates and the reparametrization trick produce highly accurate approximations. Control variates and the reparametrization trick lead to a higher average ELBO ($-781.57$ for the control variates approach and $-781.56$ for the reparametrization trick implementation) than the mean-field approximation ($ -782.26$), and the reparametrization trick algorithm achieved faster convergence in general.
Figure (ref) shows the efficiency of our SVB method in approximating the posterior densities of the parameters for a randomly selected T GARCH-type dataset. Both the control variates and reparametrization trick approaches exhibit an impressive level of agreement with the MCMC approximation, which utilized 1 million Monte Carlo samples. The overlapping of the posterior densities suggests that our SVB methods provide accurate and reliable estimates of the model parameters.
Figure (ref) illustrates the approximate posterior densities for SVB applied to one replicate of the simulation study generated from the GARCH model with skew t innovations. The figure also includes a plot showing the SVB lower bound as a function of the number of iterations for both the control variates and reparametrization trick implementations. All approximations appear to provide accurate results. The convergence rate suggests that the reparametrization trick reaches convergence faster compared to the algorithm using control variates.
Figure (ref) - (ref) (in appendix (ref)) present boxplots comparing the accuracy results of SVB algorithms for a GARCH(1,1) model with different innovations. For shorter time series, we did not observe a clear advantage in terms of accuracy between the control variates and reparametrization trick approaches. However, as the length of the time series increases, the reparametrization trick method consistently outperforms the control variates method in terms of accuracy. This is evident from the boxplots, where the accuracy results obtained using the control variates method are consistently lower and exhibit a wider spread compared to the reparametrization trick method.
Table (ref) provides a summary of the average accuracy levels of the posterior approximations for all the parameters in the three considered models. The table also includes a comparison between time series with 1,000 observations and the longer ones with 5,000 observations. It is evident that the mean-field approximation underestimates the uncertainty of the posterior distribution and exhibits the lowest accuracy. Across all datasets, the accuracy is consistently below 89%, with $\omega$ being the most challenging parameter to estimate. This issue persists across datasets of different lengths.
On the other hand, both the control variates and reparametrization trick methods achieve higher accuracy results when using an unrestricted covariance matrix. For both methods, the accuracy results for most parameters are consistently above 88% across different time series lengths except for the skewness parameter $\xi$ in the longer time series with the reparametrization trick approach. Interestingly, the performance of the control variates and reparametrization trick methods varies depending on the length of the time series. The control variates approach demonstrates higher accuracy for shorter time series, while the reparametrization trick method performs better for longer time series. This finding suggests that the choice of method may depend on the specific characteristics and length of the time series under consideration.
Table (ref) presents the computational running times of all the considered implementations, including a comparison with standard MCMC. The MCMC was stopped after 50,000 iterations, as suggested in tan2018gaussian and ong2018gaussian, and convergence was also checked using the Gelman-Rubin diagnostics and autocorrelations brooks1998general.
The table presents the average computational times for various approximation methods applied to 1,000 datasets with different numbers of Monte Carlo samples $S$. It is observed that the computational time increases as the number of MC samples increases. Comparing the same number of MC samples, the control variates approach exhibits faster computational times than the reparametrization trick method. However, it is important to note that the accuracy results of the control variates approach may not meet the same standards, particularly for longer time series. By selecting MC samples of $S=10$ for the control variates approach and $S=5$ for the reparametrization trick approach, we are able to achieve reasonable accuracy results while significantly reducing the computational time compared to the traditional MCMC sampling approximation method.
For time series with 1,000 observations, SVB with control variates is 4.37, 4.57, and 7.24 times faster than the traditional MCMC method for Gaussian, t, and skewed t GARCH models, respectively. The reparametrization trick method performs even faster, precisely 4.55, 4.64, and 7.93 times faster than MCMC. For the longer time series with 5,000 observations, SVB with control variates is 7.34, 7.40, and 9.51 times faster than MCMC in terms of average computational time, while SVB with the reparametrization trick is 4.86, 5.14, and 9.11 times faster for Gaussian, t, and skewed t GARCH models, respectively.
While offering a notable reduction in computational time compared to traditional MCMC methods, our SVB methods maintain reasonable accuracy in approximating the posterior densities in GARCH-type models. This combination of accuracy and computational speed makes SVB an attractive alternative for analyzing GARCH models in various applications.
We now investigate the performance of the sequential alternatives to SVB described in Section (ref) and compare the accuracy and computational times of UVB and Seq-SVB. We generate 1,000 datasets from the GARCH model with skewed t innovations and 1,000 observations and an equal number of replicates from the same model but with 5,000 observations. We perform this comparison focusing on the skewed t case, which comprises both Gaussian and t GARCH models as special cases.
For the sequential VB method, we exclusively employ the reparametrization trick approach with a full covariance matrix. This choice is motivated by the fact that the mean-field batch version offers significantly lower accuracy compared to the unrestricted case, as shown in the previous simulated experiments. The decision to choose the reparametrization trick approach over the control variates method for our sequential variational Bayes implementation was based on the availability of prior discussions and implementations as the control variates method has already been extensively discussed and implemented in tomasetti2022updating. On the other hand, the reparametrization trick for UVB has not been explored in previous work. Therefore, to expand the scope of our study and explore alternative approaches, we opted to implement the reparametrization trick for the sequential version.
For the time series with 1,000 observations, i.e., $T_N=1,000$, we choose ten different starting points $T_n = {500,550,600,\ldots,950}$ and perform sequential updates every $(T_N-T_n)/c$ observations, with $c={1,2,5,10}$. For example, if $T_N=1,000$, $T_n=600$ and $c=2$, the posterior distribution is first approximated using batch SVB on the first 600 observations, and then the posterior distribution is updated twice as new data in blocks of 200 observations arrive (see Figure (ref)).
Figure (ref) shows the average accuracy obtained across 1,000 simulated skew t GARCH datasets, each containing 1,000 observations. The accuracy is measured for different starting values and various numbers of sequential updates. Seq-SVB methods consistently outperform UVB in terms of accuracy. This can be attributed to the fact that Seq-SVB directly targets the true posteriors instead of aiming to approximate the pseudo-posterior distributions. Remarkably, Seq-SVB methods exhibit nearly identical accuracy regardless of the number of updates or the choice of starting points. This result is expected, as the only difference among Seq-SVB approaches is the initial value of the variational parameters for the SVB algorithm, which ultimately should converge to the same value. On the other hand, UVB approaches produce accuracy results that depend on the number of observations in the updates. As the number of sequential updates increases, the accuracy tends to decrease since the posterior approximation at each update is based on fewer observations, and each sequential update introduces a slight deviation from the true posterior. If the initial SVB approximation is based on more observations, the overall approximation is more accurate.
Figure (ref) shows the average computational time with different starting points for a single sequential update. Both UVB and Seq-SVB outperform batch SVB in terms of computational speed. Seq-SVB exhibits a small downward slope, indicating a slight decrease in computation time as the number of new data points decreases. UVB, on the other hand, shows a more noticeable downward slope, indicating a more significant reduction in computational time with fewer newly arrived data points.
In Seq-SVB, the initial values used in the optimization algorithm can have an impact on convergence and computational time. When fewer new data points are available, the initial values might already be relatively close to the final converged values at time $T_N$. As a result, the optimization algorithm requires fewer iterations to reach convergence, leading to a slight decrease in computational time.
Nevertheless, the significant reduction in computation time observed in the UVB approach can be attributed to the lower number of observations that need to be recomputed in the log-likelihood calculations. UVB only needs to incorporate the likelihood of the newly arrived data points rather than the entire dataset, resulting in faster computations.
Table (ref) presents the average computational time for a single update using both UVB and Seq-SVB with 1,000 and 5,000 observations. Regardless of the time series length, the trends remain consistent. UVB consistently outperforms Seq-SVB in terms of computational speed per update, and the computational gains are more significant when there are fewer newly arrived data points.
This advantage of UVB becomes more evident with longer time series lengths as Seq-SVB requires processing the entire data series for each update, while UVB only needs to incorporate the incremental time series.
We apply our batch and sequential variational Bayes methods to a real financial stock return dataset, specifically a time series of the Standard & Poor 500 (S&P500) stock index. After performing the data clearance procedure, which involved removing any missing or incomplete data, as well as calculating the log-returns of the time series, we obtained a dataset comprising 1,000 observations. These observations were recorded over a period spanning from 2nd January 2015 to 4th January 2019. The log-returns of the time series exhibit a range between -4.18% and 5.69%, with a median log return of 0.03%. The distribution of the log returns is negatively skewed, as evidenced by a skewness value of -0.33. Additionally, the kurtosis of the distribution is higher than that of a normal distribution, indicating excess kurtosis, with a kurtosis value of 4.59. We implemented the SVB approaches using Gaussian, T, and skewed T models on the provided dataset. Model selection criteria, namely the Akaike Information Criterion (AIC) and Bayesian Information Criterion (BIC), were used to assess the performance of the models. The AIC and BIC values for each model are reported in Table (ref).
The skewed T GARCH model returns the lowest values for both AIC and BIC. With an AIC value of 2165.44 and a BIC value of 2189.98, the skewed T GARCH model is selected as the best-fitting model for the dataset.
Figure (ref) showcases the estimated marginal posterior density plots for the parameters $\omega, \alpha, \beta, \nu,$ and $\xi$ of the skewed t GARCH model, using MCMC and batch SVB with control variates and the reparametrization trick. The density plots obtained from both the control variates method and the reparametrization trick closely resemble the estimates obtained from the MCMC method, indicating a high level of accuracy in approximating the posterior densities.
The accuracy results for the control variates method are 97.56% for $\omega$, 95.86% for $\alpha$, 97.42% for $\beta$, 93.23% for $\nu$ and 93.44% for $\xi$. The reparametrization trick implementation yields similar accuracy results: 94.06% for $\omega$, 98.13% for $\alpha$, 95.42% for $\beta$, 90.49% for $\nu$ and 92.80% for $\xi$. These accuracy percentages demonstrate the effectiveness of both implementations in providing reasonable approximations of the posterior densities from the skewed t GARCH model when applied to a real S&P500 dataset.
Figure (ref) illustrates the accuracy results obtained from applying the skewed t GARCH model to the real S&P 500 dataset using the Seq-SVB and UVB methods. Different starting values and varying numbers of sequential updates are considered.
Consistent with previous findings, the Seq-SVB approaches consistently outperform the UVB methods in terms of accuracy. The accuracy results in the UVB methods are significantly influenced by the choice of starting values and the number of updates. On the other hand, the Seq-SVB approach demonstrates robust accuracy regardless of the starting values and the number of updates.
The upward slope observed in the plot indicates that using fewer newly arrived data points leads to more accurate approximations. This suggests that incorporating a smaller portion of the newly observed data in each update results in improved accuracy. The Seq-SVB approach benefits from this property, as it only updates the variational parameters with a subset of the newly arrived data, leading to more accurate posterior approximations.
These results provide further evidence that the Seq-SVB methods are well-suited for sequential inference tasks, such as analyzing real financial datasets. They offer accurate approximations of the posterior distributions, and the accuracy is relatively insensitive to the choice of starting values and the number of updates.
Table (ref) and (ref) summarize the mean accuracy results for the parameters of the skewed t GARCH model using UVB and Seq-SVB approaches, compared to the batch SVB method. The accuracy results are obtained by considering ten different starting values ranging from $T_n = 500$ to $T_n = 950$.
The findings indicate that the accuracy of the UVB method decreases as the number of updates increases. This suggests that the UVB method may struggle to provide accurate results when a larger number of updates are applied. On the other hand, the Seq-SVB method exhibits more consistent accuracy results across different starting values, indicating its robustness and insensitivity to the number of updates.
However, it is important to note that the Seq-SVB method requires a longer computational time compared to UVB. The mean computational time for Seq-SVB, considering ten different starting values, is 47.62 seconds for a single update, while UVB only requires 17.11 seconds. Therefore, the choice between Seq-SVB and UVB should consider the trade-off between computational time and accuracy. Seq-SVB may provide higher accuracy but requires more computational resources, while UVB offers a more time-efficient solution at the expense of slightly lower accuracy.
In summary, the Seq-SVB method demonstrates robust accuracy across different starting values and update frequencies, while UVB is more sensitive to the number of updates but computationally faster. Researchers should consider these trade-offs based on their specific requirements and constraints when selecting between Seq-SVB and UVB for their applications.
This paper introduces and compares batch stochastic variational Bayes (SVB) algorithms for approximating posterior densities in GARCH-type models, including Gaussian, t, and skewed t GARCH models. The control variates and reparametrization trick methods are implemented and compared with traditional MCMC estimation. Both SVB methods demonstrate high accuracy in approximating the posterior densities of GARCH time series models.
For the models examined in this study, both the control variates and reparametrization trick approaches demonstrate notable improvements in computational efficiency compared to traditional MCMC methods. The computational time is influenced by the number of Monte Carlo samples $S$. In general, the control variates method exhibits faster computation times compared to the reparametrization trick method when the same number of MC samples is employed. However, the control variates method may struggle to maintain the same level of accuracy when using a small number of MC samples, especially for the longer time series.
We have provided sequential versions of the SVB approaches, namely UVB and Seq-SVB, to handle streaming data or scenarios with frequent updates. Both UVB and Seq-SVB outperform the batch SVB methods in terms of computational speed. UVB achieves significant time savings by only evaluating the likelihood of new observations in each update. However, UVB may sacrifice some accuracy compared to Seq-SVB, especially with a high number of sequential updates or a substantial number of new observations.
A promising approach is to run UVB and Seq-SVB in parallel, leveraging the strengths of both methods. Additionally, refreshing the pseudo-posterior in UVB using the posterior obtained from Seq-SVB can potentially improve accuracy while still benefiting from the computational efficiency of UVB.
In summary, stochastic variational inference methods, particularly the reparametrization trick approach, provide an effective and efficient framework for posterior approximation in GARCH time series models. The choice between UVB and Seq-SVB depends on the specific trade-off between computational speed and accuracy requirements, making them valuable tools for researchers and practitioners in various applications involving time series analysis.
\backmatter
\bmhead{Acknowledgments} This research includes computations using the computational cluster Katana supported by Research Technology Services at UNSW Sydney.