Semiparametric GARCH via Bayesian model averaging
Wilson Ye Chen Richard H. Gerlach


1]Wilson Ye Chen
2]Richard H. Gerlach
[1]School of Mathematical and Physical Sciences, University of Technology Sydney
[1]Australian Research Council Centre for Excellence in Mathematical and Statistical Frontiers
[2]Discipline of Business Analytics, The University of Sydney

Semiparametric GARCH via Bayesian model averaging
    [
    August 24, 2017
=================================================




As the dynamic structure of the financial markets is subject to dramatic changes, a model capable of providing consistently accurate volatility estimates must not make strong assumptions on how prices change over time. Most volatility models impose a particular parametric functional form that relates an observed price change to a volatility forecast (news impact function). We propose a new class of functional coefficient semiparametric volatility models where the news impact function is allowed to be any smooth function, and study its ability to estimate volatilities compared to the well known parametric proposals, in both a simulation study and an empirical study with real financial data. We estimate the news impact function using a Bayesian model averaging approach, implemented via a carefully developed Markov chain Monte Carlo (MCMC) sampling algorithm. Using simulations we show that our flexible semiparametric model is able to learn the shape of the news impact function from the observed data. When applied to real financial time series, our new model suggests that the news impact functions are significantly different in shapes for different asset types, but are similar for the assets of the same type.



Keywords: volatility, news impact function, Markov chain Monte Carlo, functional coefficient, regression spline, heavy tail





§ INTRODUCTION


Extensions of the GARCH model have been actively developed in the literature for the past three decades since the seminal studies by <cit.> and <cit.>. As the normality assumption of the innovation distribution of the GARCH model is often inadequate for describing the excess kurtosis of the daily returns, one stream of studies focuses on GARCH models with flexible parametric innovation distributions. Examples include the Student's t distributions <cit.>, generalised error distribution (GED) <cit.>, skewed-t distribution <cit.>, normal inverse Gaussian (NIG) distribution <cit.>, exponential generalized beta distribution of the second kind (EGB2) <cit.>, asymmetric Laplace distribution <cit.>, and the two-sided Weibull distribution <cit.>. Another line of research focuses on developing flexible conditional volatility dynamics. Examples include the EGARCH <cit.>, GJR-GARCH <cit.>, threshold-GARCH <cit.>, NAGARCH <cit.>, smooth transition GARCH <cit.>, GARCH model with slowly-varying volatility component <cit.>, Beta-t-GARCH <cit.>, and the more general family of score driven models <cit.>.

If either the innovation distribution or the conditional volatility dynamics is allowed to be of unknown form and estimated nonparametrically, the resulting model is referred to as a semiparametric GARCH model. To avoid efficiency loss or inconsistency due to misspecification, several studies propose to estimate the innovation density with a nonparametric estimator, such as the kernel density estimator. For example, see <cit.>, <cit.>, <cit.>, and <cit.> for details. Let us introduce a few notations before proceeding with the semiparametric GARCH models of the latter kind, which is the focus of this paper.

Consider a discrete-time real-valued stochastic process denoted by {y_t  t ∈ℤ} (The abbreviated notation {y_t} is used hereafter.), such that y_t = σ_tϵ_t, where {ϵ_t} is an i.i.d. sequence with (ϵ_t) = 0 and (ϵ_t) = 1. By construction, {y_t} is a martingale difference sequence. It follows that σ^2_t is the conditional variance of y_t; i.e., (y_tℱ_t-1) = σ^2_t, where ℱ_t-1 = σ({y_s  s ≤ t-1}) denotes the smallest σ-algebra containing the past observations of the process, and represents the available information at t-1. Assume that σ^2_t is the value of a function of possibly infinitely many past observations,

    σ^2_t = h(y_t-1,  y_t-2, … ad inf),

where h(·) is referred to as the volatility function.

Most extensions of the GARCH model, including those with nonparametrically estimated densities, make particular parametric assumptions about the form of the volatility function h(·), which plays a key role in the modelling of the conditional volatility process. With the hope of not having the need to choose amongst the various alternative parametric proposals, several studies focused on the nonparametric generalisations of the volatility function. Studies such as <cit.> and <cit.> assumed that h(·) is a smooth but unknown J-dimensional function for some J < ∞,

    σ^2_t = h(y_t-s_1 y_t-s_J),

with the lags {s_1 s_J} selected a priori, and estimated h(·) nonparametrically using a J-dimensional smoother. For example, notice that (<ref>) can be written in the form

    y^2_t = h(y_t-s_1 y_t-s_J) + u_t,

where u_t = σ^2_t(ϵ^2_t - 1) and {u_t} is a martingale difference sequence, which suggests that one can estimate h(·) using a kernel regression method by regressing y^2_t on the lagged observations y_t-s_1 y_t-s_J. In practice, the model in (<ref>) can only be estimated for a small J as it suffers from the so called “curse of dimensionality"; the optimal rate of convergence of a consistent nonparametric estimator of h(·) decreases rapidly as the dimensionality J increases; see for example <cit.> and <cit.>. As the conditional variance processes estimated from many empirical datasets are highly persistent, a large value of J is typically needed for the model in (<ref>) to adequately capture the empirically observed high persistence. In order to improve the best achievable rate of convergence, studies such as <cit.>, <cit.>, <cit.>, and <cit.> focused on additive models of the form

    σ^2_t = ∑_i=1^Jψ_i(θ)m(y_t-i),

for some J < ∞, where the coefficient functions {ψ_1(θ) ψ_J(θ)} are assumed known up to a parameter vector θ, and the smooth but unknown univariate function m(·) is to be estimated nonparametrically. Studies such as <cit.>, <cit.>, and <cit.> considered a more general case where J = ∞ in (<ref>). Their models include the popular GARCH(1,1) model as a special case, while (<ref>) is a semiparametric ARCH(J) model. Notice that in (<ref>) if we let J = ∞, θ = β, and ψ_i(β) = β^i-1, with β∈ (0,1), then the volatility function in (<ref>) can be written in the form σ^2_t = βσ^2_t-1 + m(y_t-1). The function m(·) in (<ref>) is often referred to as the “news impact function" in the literature, as it updates the conditional volatility based on the new information contained in the observation y_t-1. Another alternative strategy for modelling the volatility function h(·) was investigated in <cit.> and <cit.> where h(·) was assumed to be a smooth but unknown bivariate function

    σ^2_t = h(y_t-1, σ^2_t-1),

and was estimated nonparametrically using a bivariate smoothing method.

In this paper, we propose a novel semiparametric GARCH model, where the sequence of unobserved conditional variances {σ^2_t} follows a functional coefficient autoregressive (FAR) process. We model the coefficient function by a quadratic regression spline, where both the number and the positions of the knots are selected using a Bayesian model averaging (BMA) approach. A crucial step in BMA is to compute the expectations of random variables of interest with respect to the joint posterior distribution of the model parameters and the knot placements. For this purpose, we develop an adaptive Markov chain Monte Carlo (MCMC) algorithm to simulate from such joint posterior distribution. Similar to the additive model in (<ref>), the proposed approach only requires the estimation of an unknown univariate function, however, unlike (<ref>), our model nests the GARCH(1,1) model. Compared with a kernel regression method with a global bandwidth, a regression spline with strategically placed knots has the advantage of being locally adaptive <cit.>, and thus can accommodate varying degrees of smoothness across regions of the domain.

The rest of the paper is structured as follows. Section <ref> presents our semiparametric GARCH model, shows how various parametric models can be written as special cases of the proposed model, and defines the spline functional coefficient. Section <ref> describes the model in a Bayesian context, and explains how the scalar parameters and the functional coefficient are estimated using a MCMC sampling algorithm. Section <ref> presents a simulation study comparing the performance of the proposed model to the popular parametric alternatives. Section <ref> demonstrates the proposed model by applying it to ten real financial time series, and presents a model comparison study based on the averaged DIC – a slight generalisation of the usual single model DIC.



§ THE PROPOSED MODEL




 §.§ Functional coefficient semiparametric GARCH


Continuing with the notations introduced in Section <ref>, let {r_t  t ∈ℤ} denote a real-valued discrete-time stochastic process, and ℱ_t = σ({r_s  s ≤ t}) denote the natural filtration of {r_t}. Our proposed functional coefficient semiparametric GARCH (denoted by SP-GARCH hereafter) model is defined as follows.

    r_t   = μ + y_t, 
    
    y_t   = σ_tϵ_t, 
    σ^2_t   = ω + g(ϵ_t-1)σ^2_t-1,

where {ϵ_t} is an i.i.d. sequence of innovations with (ϵ_t) = 0 and (ϵ_t) = 1. It follows that μ and σ^2_t are, respectively, the conditional mean and conditional variance of r_t; i.e., (r_tℱ_t-1) = μ and (r_tℱ_t-1) = σ^2_t. In (<ref>), we make the assumption that the conditional mean of r_t is constant over time. The constant mean assumption is often reasonable in empirical applications such as when {r_t} is a model for a sequence of daily returns of a financial asset <cit.>. The distribution of the innovation ϵ_t can either be chosen from a wide range of parametric families with finite variance such as those mentioned in Section <ref> or modelled nonparametrically. In this paper, we take a parametric approach and assume that ϵ_t∼ F_sdt,ν, for ν > 2, where F_sdt,ν denotes a standardised Student-t distribution with ν degrees of freedom. (A random variable z follows F_sdt,ν if it is given by z = x[(ν - 2) / ν]^1/2, where x follows a Student-t distribution with ν degrees of freedom. It follows that (z) = 0 and (z) = 1.) The typical assumption that ϵ_t follows a heavy-tailed distribution is motivated by the stylised fact that the distribution of daily returns, standardised using an estimated volatility model, is heavy-tailed and approximately symmetric for many financial assets <cit.>.

The proposed model in (<ref>) admits the interpretation that the sequence of unobserved conditional variances {σ^2_t} follows an autoregressive process with a functional coefficient g(·). It can be seen that the function g(·) plays a similar role to that of the news impact function m(·) in (<ref>), as it determines how σ^2_t-1 is updated given the value of the innovation ϵ_t-1, which can be considered as news. Notice that the argument of g(·) is the lagged innovation ϵ_t-1, as opposed to the lagged centred observation y_t-1. We consider this as an advantage of the new model over the previous semiparametric proposals. For example, from a modeller's perspective, since the distribution of ϵ_t is either fully known or at least known up to its third moment (by construction), and is usually not time dependent, one can utilise such prior knowledge to improve the statistical properties of an estimator by imposing a certain structure on g(·). Furthermore, because the scale of ϵ_t is the same across different datasets, the estimates of g(·) can be directly compared in empirical applications.

Notice that g(ϵ_t) is an i.i.d. random variable. Let us suppose that g(ϵ_t) is replaced by a constant γ and the process {σ^2_t} has an initial value of σ^2_0 > 0, so that {σ^2_t} becomes a fully deterministic process. For γ∈ (0,1), we can see that γ controls the rate at which σ^2_t approaches its lower bound ω / (1 - γ) as t →∞. Since g(ϵ_t) is a random variable, the quantity [g(ϵ_t)] controls the average rate of decay of {σ^2_t}, and is referred to as the persistence of the conditional variance process. The process {r_t} defined in (<ref>) is covariance stationary if

    [g(ϵ_t)] ∈ [0, 1).

If {r_t} is covariance stationary, its unconditional variance is given by

    (r_t) = (σ_t^2) = ω/1 - [g(ϵ_t)].


The GARCH(1,1) model and many of its popular parametric extensions can be written as special cases of the new SP-GARCH model in (<ref>) with particular forms of the coefficient function g(·). For example, Table <ref> lists the forms of the coefficient function when the first order GARCH, GJR-GARCH, NAGARCH, and Beta-t-GARCH models are written as special cases of (<ref>), where β, α, α_1, α_2, c, and ν are parameters, taking their values typically in a subset of ℝ, and I_(-∞,0)(ϵ) takes the value 1 if ϵ < 0 and 0 otherwise. Figure <ref> shows the corresponding plots of these functions in the interval [-4, 4] and how their shapes change for the various values of their parameters.

Notice that all the models except for the GARCH allow a negative innovation to contribute differently to the conditional variance from a positive innovation of the same magnitude. For the GJR-GARCH and the Beta-t-GARCH models, the asymmetry of g(·) is created by scaling an even function differently on each side of zero. The scaling factor is α_1 + α_2 for a negative argument, and is α_1 for a positive argument. Moreover, the NAGARCH model is the only one that allows the minimum of g(·) to occur away from 0. For the Beta-t-GARCH model, the shape of g(·) is directly linked to the innovation distribution, which is assumed to be F_sdt,ν.

The Beta-t-GARCH model can be considered as a specific member of the family of score driven models known as generalised autoregressive score (GAS) models in <cit.> and dynamic conditional score (DCS) models in <cit.>. In a score driven volatility model, g(·) is a transformed gradient function (i.e., score) of the conditional distribution of r_t with respect to σ_t^2. The gradient function of the Student-t distribution dampens the effect of extreme observations on conditional variance when the degree of freedom parameter is small. Notice that u →ϵ^2 as ν→∞. Hence, the Beta-t-GARCH model is reduced to GJR-GARCH when the innovation distribution is Gaussian.

It is straightforward to show that the stationarity condition and the unconditional variance of each of the parametric models is consistent with those of the new model in (<ref>) and (<ref>).









 §.§ Spline coefficient function


In this paper, we generalise the various parametric proposals by assuming that the coefficient function g(·) in (<ref>) is an unknown but smooth univariate function. We model g(·) nonparametrically by a polynomial regression spline. A pth degree polynomial spline, for p ∈ℕ_0 = {0, ℕ}, is a function constructed piece-wise where each piece is a pth degree polynomial function. The set of points in the domain at which the pieces are joint together are referred to as knots. It is worth mentioning at this point that we consider a so called regression spline approach where there are significantly fewer number of knots than the number of observations, and the knot positions are selected as part of the modelling task. For example, a simple strategy is to place the knots either equally spaced or at specific sample quantiles. This is different from the smoothing spline approach (See <cit.>, for example.), where the knots are the observations themselves and a roughness penalty is added explicitly to reduce the effective degrees of freedom.

Specifically, the coefficient function in our model is given by

    g(ϵ) = ∑_i=0^p b_iϵ^i + ∑_i=1^Kβ_i (ϵ - k_i)^p_+,

where p ∈ℕ_0 is the degree of the spline, K ∈ℕ_0 is the number of knots, {k_1 k_K} is the sequence of knots satisfying k_1 < ⋯ < k_K, and (ϵ - k_i)^p_+ is known as a pth degree truncated power function defined by

    (ϵ - k_i)^p_+ = (ϵ - k_i)^pI_[k_i,∞)(ϵ),

with I_[k_i,∞)(ϵ) taking the value 1 if ϵ≥ k_i and 0 otherwise. In (<ref>), the pth degree spline is represented as a linear combination of the functions

    ϵ^0ϵ^p,  (ϵ - k_1)^p_+ (ϵ - k_K)^p_+,

where b_0 b_p, β_1β_K are the coefficients in ℝ. It can be shown that the set of functions in (<ref>) is a basis for the space of pth degree splines with knots k_1 k_K (See <cit.> and <cit.>, for example.). Furthermore, any linear combination of the functions in (<ref>) has p - 1 continuous derivatives. The particular basis whose members are given by (<ref>) is known as a truncated power basis. In this paper, we restrict our attention to quadratic splines where p = 2 and

    g(ϵ) = b_0 + b_1ϵ + b_2ϵ^2 + ∑_i=1^Kβ_i (ϵ - k_i)^2_+.

It can be shown that, when p = 2, several well known parametric models can be expressed exactly in terms of the spline coefficients. For example, following from the definitions in Table <ref>, the first order GARCH corresponds to the case where b_1 = 0 and K = 0; GJR-GARCH corresponds to the case where b_0 = β, b_1 = 0, b_2 = α_1 + α_2, β_1 = -α_2, K = 1, and k_1 = 0; NAGARCH corresponds to the case where b_0 = β + α c^2, b_1 = -2 α c, b_2 = α, and K = 0.

Assume that the distribution of ϵ_t is both standardised and symmetric. Using the definition in (<ref>) and the distributional assumption of ϵ_t, we can derive the persistence of the SP-GARCH conditional variance process as follows.

    [g(ϵ_t)] = ∫_-∞^∞g(ϵ_t)f(ϵ_t)dϵ_t = b_0 + b_2 + ∑_i=1^Kβ_i c_i,

where f(·) is the density function of the innovation distribution, and c_i is given by

    c_i = ∫_k_i^∞ (ϵ_t - k_i)^2 f(ϵ_t) d ϵ_t.

The integral for c_i in (<ref>) can be obtained in closed form if f(·) is the density function of the standard normal distribution.

    c_i = - 1/√(2 π)exp(-k_i^2/2) k + 1/2 (1 + k_i^2) ( k_i/√(2)),

where (x) = 1 - (x) denotes the complementary error function. If f(·) is instead the density function of F_sdt,ν, as per our assumption in Section <ref>, the expression for c_i is rather complicated, and the evaluation of such expression is subject to numerical difficulties. In such case, as a practical alternative, the value of c_i can be approximated with high accuracy using a numerical integration algorithm, such as the univariate adaptive quadrature. Figure <ref> shows plots of c_1 c_9 computed using the adaptive quadrature for F_sdt,ν with ν∈{2.2,  3,  8}, and using the closed-form expression in (<ref>) for the standard normal distribution, where the corresponding knots k_1 k_9 are placed at the quantile levels 0.1  0.9 of F_sdt,8. It can be seen that, even at ν = 8, the values of {c_i} are already very close to those computed analytically for the standard normal distribution, which suggests that, for the t distribution, equation (<ref>) can be used to approximate the integral for c_i when ν is moderate or large (for example, ν≥ 8).






§ BAYESIAN INFERENCE




 §.§ Bayesian model averaging


Our goal in Section <ref> is to estimate the degrees of freedom ν, the mean μ, the constant parameter ω of the volatility function, and the coefficient function g(·) of the SP-GARCH model defined in (<ref>), where g(·) is modelled as a quadratic regression spline with a sequence of knots {k_1 k_K} defined in (<ref>). The choice of the knot sequence is important in estimating functions by regression splines, as it imposes a certain structure on the estimates. Ideally, we would want a small number of well positioned knots, so that we have a parsimonious model that can also account for the important details in the data. We employ a Bayesian variable selection approach where model averaging occurs naturally, and in doing so, we account for two types of uncertainty: the uncertainty in the knot sequence {k_1 k_K} and the uncertainty in the parameter values conditional on a particular knot sequence. For a detailed treatment on Bayesian model selection and model averaging, see <cit.>, <cit.>, and the references therein.

Suppose that there are M ∈ℕ candidate models, indexed by τ∈{1  M}. Let θ_τ∈Θ_τ denote the parameter vector associated with the τth model, where the set Θ_τ defines the parameter space, and let 𝐫 = (r_1 r_T) denote a vector of observations. Further suppose that we are interested in estimating a quantity ϕ, such that conditional on 𝐫, ϕ = η(θ_τ, τ), for some function η(·). For example, ϕ may be the one-day-ahead forecast of the conditional volatility, or the value of g(ϵ_t) at some specific value of ϵ_t. Under the assumption of model uncertainty, the posterior mean of ϕ can be taken as a Bayesian point estimate, which is obtained by averaging over all the candidate models;

    ϕ̂_Bayes = (ϕ𝐫)    = ∑_τ=1^M[∫_Θ_τη(θ_τ, τ)p(θ_τ, τ𝐫)dθ_τ] 
       = ∑_τ=1^M[∫_Θ_τη(θ_τ, τ)p(θ_ττ, 𝐫)dθ_τ]p(τ𝐫) 
       = ∑_τ=1^M(ϕτ, 𝐫)p(τ𝐫).

While for some problems one can obtain the closed-form expressions for both the conditional expectation (ϕτ, 𝐫) and the posterior model probabilities p(τ𝐫) (For example, see <cit.>.), for many realistic models they are analytically intractable, and approximation methods must be used. Even in the case of <cit.>, the summation in (<ref>) is infeasible to compute due to a large number of candidate models. We develop an adaptive MCMC algorithm to simulate from the joint posterior distribution p(θ_τ, τ𝐫) defined on the space

    ⋃_τ∈{1  M}Θ_τ×{τ}.

If {(θ_τ^[i], τ^[i])  i ∈{1  N}} denotes a sample from p(θ_τ, τ𝐫) generated by a MCMC sampling algorithm, an approximation of ϕ̂_Bayes can then be computed as

    ϕ̂_MCMC = 1/N∑_i=1^Nη(θ_τ^[i], τ^[i]).

It can be seen from (<ref>) that model averaging is achieved without additional efforts via a single run of the MCMC sampler.

An advantage of the basis representation in (<ref>) is that candidate models can be represented as submodels (or restricted models) nested within a supermodel. If we consider a model with a large number of knots (i.e., a large K) as a supermodel, then a candidate model with a specific knot sequence can be considered as a submodel with a set of knot coefficients {β_i} restricted to zero for some i ∈{1  K}. In a Bayesian approach, the knot selection can be achieved by allowing a positive prior probability that each knot coefficient β_i is exactly zero.



 §.§ Joint posterior


In this section, we define the the joint posterior distribution of the model parameters and the knot sequences. We adopt the paradigm where there is a large number of potential knots, and a subset of which is assigned to a candidate model. Let θ = (ν, μ, ω,  b_0,  b_1,  b_2, β_1β_K) denote the parameter vector of the full model, and 𝐤 = (k_1 k_K) denote the vector of the potential knots. Each candidate model is indexed by a binary vector 𝐦 = (m_1 m_K) such that

    m_i =
    
    0    if β_i = 0,
    
    1    if β_i 0.

The joint posterior distribution of the pair (θ, 𝐦) given data is then proportional to the likelihood multiplied by the joint prior;

    p(θ, 𝐦𝐫) ∝ p(𝐫θ, 𝐦) p(θ, 𝐦).


Following from the model definition in (<ref>) and (<ref>), the likelihood is given by

    p(𝐫θ, 𝐦) = p(𝐫θ) = I_Θ(θ) ∏_t=1^Tp(r_tℱ_t-1, θ),
 
where the conditional distribution of each observation r_t is a Student-t distribution with ν degrees of freedom, whose mean is μ and whose variance is σ^2_t, and I_Θ(θ) is equal to one if θ∈Θ and zero otherwise. It should be made clear in (<ref>) that 𝐫 and 𝐦 are independent conditional on θ; the vector θ encodes all the information needed to generate 𝐫, as a result of the nested model set-up. We assume that ℱ_0 = σ(∅), and set the initial conditional variance σ^2_1 to the sample variance of the realisation of 𝐫. We use the indicator function I_Θ(·) to restrict the parameter space, by setting the likelihood to zero when θ∉Θ. The restricted parameter space Θ is defined as

    Θ = {θ 2 < ν≤ 200, 
    ω > 0, 
    0 < b_0 + b_2 + ∑_i=1^Kβ_ic_i < 1, 
    σ^2_t > 0, for t ∈{1  T+1}},

so that when θ∈Θ, the process {r_t} is covariance stationary with a well defined unconditional variance given by (<ref>), and the process {σ^2_t} is positive for t ∈{1  T+1}. Notice that the positivity restriction on the conditional variance process in (<ref>) is data-dependent, as σ^2_t is a function of both θ and 𝐫. This type of implicit restrictions was used in <cit.> to replace the sufficient conditions, which are explicitly available in terms of the model parameters, but more restrictive. The constants {c_i} are given by (<ref>). Since the density function f(·) in (<ref>) is of F_sdt,ν, the integral is evaluated numerically using adaptive quadrature. From an implementation perspective, these constants can be computed ex-ante on a fine grid for a range of degrees of freedom ν, and stored in a lookup table. During the evaluation of the indicator function I_{θ∈Θ}(θ), a table-lookup operation is computationally less expensive than running a quadrature algorithm.

The joint prior can be constructed as

    p(θ, 𝐦) = p(θ𝐦) p(𝐦).

When 𝐦 is known, we define the conditional prior for θ as

    p(θ𝐦) ∝ν^-2∏_i=1^K[ m_if_Gauss(β_i;  0, σ^2_β) + (1 - m_i)δ(β_i)],

where f_Gauss(β_i;  0, σ^2_β) is a Gaussian density, with mean zero and variance σ^2_β, the delta function δ(β_i) represents a point probability mass at β_i = 0. Conditional on 𝐦, the parameters are assumed to be independent a priori. The conditional prior in (<ref>) implies that β_i is known to be exactly zero a priori when m_i = 0. As the prior is centred at zero for each β_i, the value of σ^2_β controls the level of shrinkage of the nonzero knot coefficients. We treat σ^2_β as a tuning parameter whose value is set by the user. From (<ref>), we can see that the prior for ν behaves similarly to a half-Cauchy prior with a scale parameter equal to one. It can be shown that this prior is flat on ν^-1. Similar priors for ν were employed in studies such as <cit.> and <cit.> in the context of the GARCH models. For all the parameters other than ν and β_1β_K, the conditional prior in (<ref>) is flat. Finally, a flat prior is used for the model probabilities,

    p(𝐦) = 0.5^K∝ 1,

so that each model is given an equal prior probability. More generally, one can use the independent Bernoulli priors, p(𝐦) = ∏_i=1^K[π_i^m_i(1 - π_i)^(1 - m_i)], with different values of π_1π_K to express a prior preference for certain knots over the others. The main feature of the the joint prior specified in (<ref>) to (<ref>) is that each knot coefficient β_i is given a positive probability a priori to be exactly zero. This type of prior for β_i is known in the literature as the spike-and-slab prior <cit.>.



 §.§ Adaptive MCMC algorithm


As the closed-form expressions are unavailable for the exact Bayesian estimates of the type in (<ref>) based on the joint-posterior distribution of our model, we develop an adaptive MCMC sampling algorithm to simulate from p(θ, 𝐦𝐫) defined in (<ref>) to (<ref>), and compute the approximate estimates using (<ref>). The aim is to construct an ergodic Markov chain whose stationary distribution is p(θ, 𝐦𝐫), and then generate a sample path {(θ^[i], 𝐦^[i])  i ∈{1  N}} of the chain. In general, our sampling algorithm involves the following steps:


  * Let 𝐦^[1] be an arbitrary binary vector of length K, and i = 1.

  * Generate a proposal (θ^*, 𝐦^*) from the joint distribution q(θ^*, 𝐦^*𝐦^[i]) by


  * generating 𝐦^* from q(𝐦^*𝐦^[i]),

  * generating θ^* from q(θ^*𝐦^*, 𝐦^[i]).


  * Let (θ^[i+1], 𝐦^[i+1]) = (θ^*, 𝐦^*) with probability u, or (θ^[i+1], 𝐦^[i+1]) = (θ^[i], 𝐦^[i]) with probability 1 - u.

  * Increment i by one, and repeat from step 2 until i = N.


In step 2, the joint proposal distribution is constructed in a similar way to that of the prior distribution;

    q(θ^*, 𝐦^*𝐦^[i]) = q(θ^*𝐦^*, 𝐦^[i])q(𝐦^*𝐦^[i]),

where θ^* = (ν^*, μ^*, ω^*,  b_0^*,  b_1^*,  b_2^*, β_1^*β_K^*) and 𝐦^* = (m_1^* m_K^*). We can generate a realisation (θ^*, 𝐦^*) from (<ref>) by first generating a knot configuration 𝐦^*, followed by generating a parameter vector θ^* conditional on 𝐦^*, as in the steps 2(a) and 2(b). In step 2(a), conditional on 𝐦^[i], we generate a knot configuration 𝐦^* from the following Bernoulli mixture proposal density.

    q(𝐦^*𝐦^[i]) = ∏_j=1^K[m^[i]_jf_Bern(m^*_j;  1-π) + (1-m^[i]_j)f_Bern(m^*_j; π)],

where f_Bern(m^*_j; π) = π^m^*_j(1-π)^(1-m^*_j). By employing the proposal density in (<ref>), the resulting chain {𝐦^[i]} is a random walk in the K-dimensional binary vector space, where the parameter π controls the step-size. When π < 0.5, the realisations of m^[i]_j and m^[i+1]_j are more likely to be the same, and when π > 0.5, the opposite is true.

Let θ_𝐦^* denote a (6 + ∑_i=1^Km_i^*)-dimensional vector containing only those elements of θ^* that are known a priori to be nonzero (I.e., θ_𝐦^* is the parameter vector of the submodel implied by 𝐦^*), and θ_𝐦^* denote a (∑_i=1^K(1 - m_i^*))-dimensional vector containing only those knot coefficients that correspond to the zero elements of 𝐦^*. In step 2(b), when 𝐦^* is known, we write the proposal density for θ^* as

    q(θ^*𝐦^*, 𝐦^[i]) = q(θ^*𝐦^*) = q(θ_𝐦^*𝐦^*)q(θ_𝐦^*𝐦^*),

where q(θ_𝐦^*𝐦^*) = δ(θ_𝐦^*) is a multivariate delta function representing a point mass at θ_𝐦^* = (0  0). Notice that θ^* and 𝐦^[i] are conditionally independent given 𝐦^*. We generate θ_𝐦^* from a mixture of multivariate Gaussian distributions with a different scale for each component,

    q(θ_𝐦^*𝐦^*) = ∑_j=1^N_mix w_j f_mvn(θ_𝐦^*; θ̂_𝐦, ς_jΣ̂_𝐦),

where ∑_j=1^N_mix w_j = 1 and f_mvn(θ_𝐦^*; θ̂_𝐦, ς_jΣ̂_𝐦) is a multivariate Gaussian density with mean vector θ̂_𝐦 and covariance matrix ς_jΣ̂_𝐦. The number of components N_mix, the vector of weights 𝐰 = (w_1 w_N_mix), and the vector of scales ς = (ς_1ς_N_mix) are tuning parameters that can be adjusted to help the sampling algorithm efficiently explore the posterior distribution.

The mean vector θ̂_𝐦 and the base covariance matrix Σ̂_𝐦 must be set carefully to ensure that the proposal density for θ_𝐦^* in (<ref>) is a good approximation of the posterior density of the submodel

    p(θ_𝐦𝐦^*, 𝐫) ∝ p(𝐫θ_𝐦, 𝐦^*)p(θ_𝐦𝐦^*).

We first generate a sample path from a pilot MCMC run, targeting the submodel posterior in (<ref>), which is easily obtained by imposing zero restrictions on the conditional posterior of the full model p(θ𝐦^*, 𝐫). For this pilot run, we employ the random walk Metropolis (RMW) kernel and a multivariate Gaussian proposal distribution, whose covariance matrix is adaptively estimated, up to a scaling factor, using the realised states of the chain, while the scaling factor is adaptively tuned to target an acceptance rate of 0.234, as in <cit.>. After discarding an initial burn-in sample, we set θ̂_𝐦 and Σ̂_𝐦 to be the sample mean and sample covariance matrix of the pilot chain. It is still acceptable if the stationary distribution of this pilot chain is not exactly (<ref>) due to adaptation <cit.>, since the sole purpose of the pilot chain is to learn good parameter values of θ̂_𝐦 and Σ̂_𝐦. To reduce computational cost, we save the values of θ̂_𝐦 and Σ̂_𝐦 for each newly visited knot configuration so that the rerunning of the RWM sampler is not required when the same configuration is revisited.

Note that conditional on 𝐦^*, the proposal distribution is fixed, and the sampling algorithm effectively becomes an independence sampler <cit.> targeting (<ref>). Similar to the method of importance sampling, the independence sampler is known to be efficient if the proposal density is a good approximation of the target while being heavier-tailed <cit.>. We set N_mix = 3, 𝐰 = (0.85,  0.1,  0.05), and ς = (1,  10,  100), with the hope that Gaussian mixture proposal density is heavier-tailed than the target posterior in (<ref>), while we approximate the location, scale, and orientation of target by choosing θ̂_𝐦 and Σ̂_𝐦 using the pilot chain.

In step 3, as (<ref>) is a density with respect to the natural measure on the product space ℝ^6+K×{0, 1}^K, the standard Metropolis-Hastings acceptance probability is a valid option for u <cit.>, and is given by

    u    = min[p(θ^*, 𝐦^*𝐫) q(θ^[i], 𝐦^[i]θ^*, 𝐦^*)/p(θ^[i], 𝐦^[i]𝐫) q(θ^*, 𝐦^*θ^[i], 𝐦^[i]),  1 ] 
       = min[p(θ^*, 𝐦^*𝐫) q(θ^[i]𝐦^[i], θ^*, 𝐦^*) q(𝐦^[i]θ^*, 𝐦^*)/p(θ^[i], 𝐦^[i]𝐫) q(θ^*𝐦^*, θ^[i], 𝐦^[i]) q(𝐦^*θ^[i], 𝐦^[i]),  1 ].

The set of conditional independencies implied by the construction of the proposal distribution in step 2 is represented graphically in Figure <ref>, using a directed acyclic graph (DAG). Using the DAG, the expression for u can be simplified to

    u = min[ p(θ^*, 𝐦^*𝐫) q(θ^[i]𝐦^[i]) q(𝐦^[i]𝐦^*)/p(θ^[i], 𝐦^[i]𝐫) q(θ^*𝐦^*) q(𝐦^*𝐦^[i]),  1 ].

Using the fact that the proposal density for the knot configuration 𝐦^*, given by (<ref>), is symmetric, the acceptance probability is further reduced to

    u = min[ p(θ^*, 𝐦^*𝐫) q(θ^[i]𝐦^[i])/p(θ^[i], 𝐦^[i]𝐫) q(θ^*𝐦^*),  1 ].







§ SIMULATION STUDY


A simulation study is conducted to investigate the comparative performance of the proposed SP-GARCH model and the popular parametric alternatives. The parametric models considered are the first order GARCH, GJR-GARCH (GJR), and Beta-t-GARCH (Beta-t). The aim is to investigate whether the SP-GARCH model is able to offer more accurate estimates of the conditional volatilities {σ_t  t ∈{1  T}} such that σ_t = [(r_tℱ_t-1)]^1/2 for certain family of processes followed by {r_t}. The focus is on comparing structural features of the estimated coefficient functions and measures of loss in both the in-sample and out-of-sample estimates of conditional volatilities. 

The model definitions of the comparing parametric models are given by (<ref>) and Table <ref>. The innovation distribution for each model is F_sdt,ν. The parameters are constrained to ensure the stationarity of {r_t} and the positivity of {σ_t^2}, by truncating the likelihood, similar to (<ref>). For each parametric model, the flat prior is used for most of the parameters with an exception of ν, where the prior is proportional to ν^-2.

Independent sequences of T observations are generated from four data generating processes (DGPs), labelled DGP 1, DGP 2, DGP 3, and DGP 4. The DGPs can be written in the functional coefficient form of (<ref>), with DGP i, for i ∈{1  4}, given by

    r_t   = σ_tϵ_i,t, 
    σ^2_t   = 0.1 + g_i(ϵ_i,t-1) σ^2_t-1,

where {ϵ_i,t} is an i.i.d. sequence from F_sdt,8 for i ∈{1,2}, and from F_sdt,5 for i ∈{3,4}. The coefficient functions g_1(·)  g_4(·) are specified as follows.

    g_1(ϵ)    = 1.1 - 0.48(ϵ + 0.77)^2_+ + 0.58(ϵ + 0.473)^2_+; 
    
    g_2(ϵ)    = 0.85 + 0.1ϵ^2; 
    
    g_3(ϵ)    = 0.82 + 0.15[6ϵ^2/(3 + ϵ^2)]; 
    
    g_4(ϵ)    = 0.8 + [0.1 + 0.15 I_(-∞,0)(ϵ)]ϵ^2.

The function g_1(·) is a quadratic spline with two knots at -0.77 and -0.473. The motivation for g_1(·) is to construct a coefficient function with structural features that cannot be fully captured by most parametric GARCH models, yet are not completely unrealistic from an empirical perspective. For example, g_1(·) implies that a relatively small negative return impacts the next day conditional volatility differently from a small positive return of the same magnitude, while the impact of a negative return stops to grow with its magnitude once the return is beyond a certain threshold. The functions g_2(·), g_3(·), and g_4(·) correspond to the coefficient functions of the first order GARCH, Beta-t, and GJR models, respectively. It can be checked that all DGPs satisfy the covariance stationarity condition in (<ref>), with [g_1(ϵ_1,t)] = 0.977, [g_2(ϵ_2,t)] = 0.95, [g_3(ϵ_3,t)] = 0.97, and [g_4(ϵ_4,t)] = 0.975.

We generate N_sim = 500 independent sequences of length T = 4001 from each of the two DGPs. The posterior mean estimates of the conditional volatilities (as well as other quantities of interest presented in this section) are computed using the first 4000 observations of each generated sequence; the last observation of each sequence is left out for quantifying the out-of-sample performance.

To compute the posterior mean estimates for the SP-GARCH model via (<ref>), the adaptive MCMC sampling algorithm detailed in Section <ref> is employed. The sampling algorithm is run for 5.5 × 10^5 iterations, where the first 0.5 × 10^5 burn-in iterations are ignored when computing the posterior mean estimates. The potential knots k_1 k_9 are chosen as quantiles of the standardised Student-t distribution with 8 degrees-of-freedom, at levels 0.1  0.9, respectively; the step-size parameter π in (<ref>) is set to 0.1; the shrinkage parameter σ^2_β in (<ref>) is set to 2500.

To compute the posterior mean estimates for the parametric models, we use an independence sampler with a multivariate Gaussian mixture proposal configured in the same way as (<ref>). The mean vector and the base covariance matrix of the Gaussian mixture are set to the sample mean and covariance matrix of a pilot chain generated by the same adaptive RWM sampler as the one used in step 2(b) of the sampling algorithm in Section <ref>. The independence sampler is shown to work well for parametric GARCH models in <cit.> and <cit.>.

The Monte Carlo (MC) mean of the 500 posterior mean estimates of the coefficient functions is plotted for each model in Figure <ref>, along with the 95% MC confidence band. Both the MC mean and the MC confidence band are computed pointwise using an equally spaced grid of points over [-4,  4] with a spacing of 0.01. For DGP 1, only the SP-GARCH is able to capture the essential features of the true coefficient function. It is interesting to see that both the GJR and Beta-t models are trying to capture the asymmetry of g_1(·) near zero, where most of the innovations are found, while missing the important features for large innovations. For DGP 2, which corresponds to the the GARCH model, all the models except for the Beta-t are able to approximate the simple quadratic function. The Beta-t model is performing poorly due to the fact that its coefficient function deviates from the quadratic from if the innovation distribution is heavier tailed than the Gaussian. Conversely, when the true model is Beta-t, as for DGP 3, the quadratic models (GARCH and GJR) provide poor approximations as g_3(·) increases faster than quadratic for small innovations, but grows much slower than quadratic for large innovations. The SP-GARCH is able to offer a good approximation to g_3(·) except for the tails of the innovation distribution where the data is sparse. Lastly, when the true model is GJR, SP-GARCH is able to learn the shape of g_4(·) reasonably well. Expectedly, the approximation becomes increasingly biased as ϵ moves deeper into the tails of the innovation distribution.




To compare the quality of the posterior mean estimates of the conditional volatilities between the models for each of the two DGPs, both the in-sample and out-of-sample versions of the L_p loss are considered. For the ith simulated sequence, the in-sample L_p loss is defined, for p ≥ 1, as

    L_p,i^(In) = [1/T-1∑_t=1^T-1|σ̂_i,t - σ_i,t|^p]^1/p,

where σ̂_i,t is the posterior mean estimate of the volatility for tth observation of the ith simulated sequence, and σ_i,t is the corresponding true volatility; for the holdout observations across different sequences, the out-of-sample L_p loss is defined, for p ≥ 1, as

    L_p^(Out) = [1/N_sim∑_i=1^N_sim|σ̂_i,T - σ_i,T|^p]^1/p.

Notice that the out-of-sample L_p loss is calculated cross-sectionally rather than time-series-wise. The cases p ∈{1, 2} are considered in this study: When p = 2, the loss measure L_p^(·) corresponds to the usual root-mean-square-error (RMSE); when p = 1, L_p^(·) corresponds to the usual mean-absolute-error (MAE).

The MC mean of the 500 realisations of the in-sample L_p loss is reported for each p ∈{1, 2}, each model, and each DGP in Table <ref>, where L̅_p^(In) denotes (1/500)∑_i=1^500L_p,i^(In). Table <ref> is the out-of-sample analogue of Table <ref>. We make some observations: (i) Comparing Table <ref> with Figure <ref>, it is evident that conditional volatility estimates are sensitive to mis-specification of the coefficient function. (ii) Compared to the parametric models, SP-GARCH is robust to such mis-specification as it can adapt to a wide range of coefficient functions. (iii) For certain class of coefficient functions, the incurred efficiency loss (by not using the correct parametric model) is very small. For example, when the true DGP is GARCH, SP-GARCH is as efficient as the GJR model, given the fact that GJR has only one more degree-of-freedom than that of GARCH, and that SP-GARCH is a much more flexible model. (iv) The relative accuracy of SP-GARCH is retained out-of-sample. This suggests that the accuracy gain is not caused by overfitting.









































§ EMPIRICAL STUDY


As an empirical study, the SP-GARCH model is applied to the daily returns of five stock market indices and five individual stocks. The aim is to gain insights into the conditional volatility process of each series by analysing the estimated coefficient function. The five stock market indices are: S&P 500 (US), FTSE 100 (UK), DAX (Germany), Nikkei 225 (Japan), and Hang Seng (Hong Kong); the five individual stocks are: Apple, ARM, Intel, Nvidia, and SanDisk. A sample of daily log-returns is obtained for each index or stock. Most of the samples start in January 1996, with exceptions of ARM and Nvidia, for which the samples start in April 1998 and January 1999, respectively. All of the samples end in November 2015. Nvidia has the smallest sample size of 4236 observations, while FTSE has the largest sample size of 5176 observations.

For each sample, the posterior mean estimate of the coefficient function and the 95% credible band are computed using the MCMC sampling algorithm detailed in Section <ref>. The sequence of potential knots, the configuration of the sampling algorithm, and the value of the hyperparameter are the same as those used for the simulation study in Section <ref>.

In Figure <ref>, the posterior mean estimate of the coefficient function is plotted for each of the ten samples, together with the pointwise 95% credible band. It is striking to see that the samples can be easily separated into two groups based on the structural features of the coefficient functions; there appear to be a clear distinction between the coefficient functions for the stock indices and those for the individual stocks. For the stock indices, the functions appear to be quadratic, with asymmetries created by shifting the minimums away from zero; for the individual stocks, the functions appear to be piecewise-linear-like, with asymmetries created by having slopes with varying degrees of steepness for negative and positive innovations. There also appear to be a “threshold effect" for large negative innovations for some individual stocks such as Intel and SanDisk.




The posterior summaries of ν and μ are reported in Table <ref>. For each parameter-sample combination, we report the posterior mean (Mean) and the 95% credible interval (Lower, Upper). One can observe that the innovation distributions for the individual stocks are heavier tailed than those for the stock indices, as indicated by the estimates of ν. The estimates of μ indicate that the long-run (unconditional) means of the daily return processes are all positive but small in magnitude, with Apple generating the largest expected daily return of around 0.1%.




One may be interested in estimating the long-run (unconditional) volatilities of the return processes as quantitative measures of the long-run daily risk; σ = (r_t)^1/2. This can be done using (<ref>) and (<ref>) for the SP-GARCH model. In Figure <ref>, a scatter plot is generated by plotting each sample in the product space of unconditional volatility and unconditional mean; {σ}×{μ}. Such plot is useful for comparing the long-run daily risk-return profiles of the samples, which could correspond to assets in an investment portfolio. It can be seen from the plot that both DAX and Apple show favourable risk-return characteristics compared to the rest of the samples.




As an illustration of the mixing properties of the constructed Markov chain when the sampling algorithm is applied to real data, the realised values of the chain in dimensions b_2 and β_3, after discarding the initial burn-in iterations, are plotted against the iteration counts minus 5.5 × 10^5 in panels (a) and (b) of Figure <ref>, when the posterior distribution is conditional on the Apple data. It appears that the chain is settled down to rapid mixing in both dimensions. In panel (b), it is reassuring to see that the algorithm appears to be able to effectively switch between multiple modes including a point mass at exactly zero. The scatter plot in panel (c) reveals the complex landscape of the joint marginal posterior distribution of b_2 and β_3.




It is useful to also examine the mixing properties of the Markov chain of the binary vectors, i.e., the Markov chain in model space. One way to visualise the sample path (or trace) of such chain is to first map each realisation of the K-dimensional binary vector 𝐦^[i] to a scalar τ^[i], for i ∈{1  N}, such that the mapping is one-to-one, and then plot τ^[i] against i, for i ∈{1  N}. One such mapping is easily found by treating each binary vector as a K-bit binary number. The mapped scalar is then the decimal equivalent of the binary number. To illustrate the idea, this mapping is used to generate the trace plots in Figure <ref> for S&P 500 and Apple. For both samples, the chain appears to be rapidly mixing in model space. 




A simple measure of model complexity for the SP-GARCH model is the number of knots chosen for the regression spline, given by

    K_𝐦 = ∑_i=1^Km_i.

Using the realised binary vectors of the Markov chain, one can easily compute the marginal posterior probability of K_𝐦. Such posterior probabilities are reported in Table <ref> for each sample for K_𝐦≤ 3. For stock indices, the most favoured model is the one without any knots; for individual stocks, more complex models are preferred, as most of the posterior mass is on K_𝐦 = 2.




It is possible that the structural features of the estimated coefficient functions shown in Figure <ref> are not unique, due to multimodality of the joint posterior distributions, and that our MCMC sampler is not effectively exploring all the modes. This possibility can potentially render our conclusions concerning the shapes of the coefficient functions invalid. For example, there could be a parametric GARCH model that fits the data equally well as the SP-GARCH, but giving a structurally different estimate of the coefficient function. One way to guard against this possibility is to compare SP-GARCH against the parametric models using a measure of model adequacy. We need to ensure that SP-GARCH is indeed a more adequate description of the data when the inference about the coefficient function is different under SP-GARCH.

We perform a model comparison study based on the deviance information criterion (DIC) of <cit.> which can be considered as a Bayesian generalisation of the Akaike information criterion (AIC) <cit.>. Successful applications of the DIC include <cit.> and <cit.> for comparing stochastic volatility models, and <cit.> for comparing Markov-switching GARCH models. When computing the DIC, one obtains both a measure of model fit and model complexity. The model with the smallest DIC achieves the best fit-complexity trade-off.

Usually, a flexible semiparametric model like the SP-GARCH would not perform well under the DIC due to its high complexity. However, the proposed BMA approach will cause the posterior mass to concentrate over simpler models when more complex models do not significantly improve the likelihood (i.e., model fit). It would thus be interesting to investigate how SP-GARCH performs against the parametric models when both fit and complexity are taken into consideration.

A few studies have investigated the definition and the computation of the DIC for latent variable models <cit.>, however the DIC has not been examined in the context of BMA. For example, it is not clear how one may compute the DIC using the output of the reversible jump MCMC sampler of <cit.>, where the models are sampled together with their parameters. A natural generalisation of the DIC in the context of BMA may be the averaged DIC defined as follows (using the more general notation in Section <ref>).

    DIC_ave = _τ(DIC_τ𝐫) = ∑_τ=1^MDIC_τ p(τ𝐫),

where τ∈{1  M} is the model index, p(τ𝐫) is the posterior probability of the model τ. The DIC for model τ, denoted by DIC_τ, is the usual single model DIC given by

    DIC_τ = D̅_τ + p_D^τ.

The first term is the conditional posterior mean of the deviance

    D̅_τ = _θ_τ[D(θ_τ, τ) τ, 𝐫],

which can be considered as a Bayesian measure of model fit. The deviance is defined as

    D(θ_τ, τ) = -2 log p(𝐫θ_τ, τ),

where p(𝐫θ_τ, τ) is the likelihood. The second term in (<ref>) is the effective number of parameters – a measure of model complexity, and is defined as

    p_D^τ = D̅_τ - D(θ̅_τ, τ),

where θ̅_τ is the conditional posterior mean of θ_τ; θ̅_τ = _θ_τ(θ_ττ, 𝐫). It can be seen from the definition in (<ref>) that, for any consistent model selection procedure, the averaged DIC in (<ref>) is asymptotically equivalent to the single model DIC of <cit.>.

It is straightforward to compute an estimate of DIC_ave whenever the posterior sample {(θ_τ^[i], τ^[i])  i ∈{1  N}} is available. Let the index set ℐ_τ⊆{1  N} be defined as follows.

    ℐ_τ = {i_1 i_N_τ τ^[i_1] = ⋯ = τ^[i_N_τ] = τ}.

The estimate of DIC_τ can be computed using the posterior draws {(θ_τ^[i], τ^[i])  i ∈ℐ_τ}, where D̅_τ≈ N_τ^-1∑_i ∈ℐ_τ D(θ_τ^[i], τ^[i]) and θ̅_τ≈ N_τ^-1∑_i ∈ℐ_τθ_τ^[i]. The posterior model probability p(τ𝐫) is approximated by N_τ / N. For any τ such that ℐ_τ = ∅, it can be assumed that p(τ𝐫) ≈ 0; thus such τ is excluded during computation. Alternatively, one can compute DIC_ave by separately computing the averaged posterior mean of deviance

    D̅_ave = _τ(D̅_τ𝐫) = _τ{_θ_τ[D(θ_τ, τ) τ, 𝐫] 𝐫} = _θ_τ, τ[D(θ_τ, τ) 𝐫],

and the averaged effective number of parameters

    p_D^ave = _τ(p_D^τ𝐫) = D̅_ave - ∑_τ=1^M D(θ̅_τ, τ) p(τ𝐫).

The averaged DIC is then given by DIC_ave = D̅_ave + p_D^ave. Similarly, both terms can be approximated using a sample of posterior draws. For example, D̅_ave≈ N^-1∑_i=1^N D(θ_τ^[i], τ^[i]).

In Table <ref>, for each parametric model, we report the estimates of the single model DIC, the posterior mean of deviance, and the effective number of parameters; for the SP-GARCH model, we report the estimates of the averaged DIC, the averaged posterior mean of deviance, and the averaged effective number of parameters. For the table and the rest of the section, the notations DIC, D̅, and p_D are overloaded to mean either the single model or the averaged version of the measures, depending on the model in question.

For all the stock indices, the DIC ranks SP-GARCH as the best model, followed by (in the order of increasing DIC) Beta-t, GJR, and GARCH models. Both measures D̅ and DIC agree on the ranking of the models, indicating that the differences in model fit dominate the differences in complexity when ranked by the DIC. For all the models that allow g(·) to be asymmetric, the DIC estimates are exceedingly lower than those for the GARCH model. This supports the well known finding that it is important to model asymmetry for equity returns.

The fact that the SP-GARCH model stands out as the best amongst the asymmetric models calls for further investigation. Recall that, in Figure <ref>, the coefficient function estimates (a) – (e) are minimised at points between 0 and 2, and that, in Table <ref>, most of the marginal posterior mass is on K_𝐦 = 0 for the indices. In other words, the most adequate model for the stock indices is obtained when g(·) is a simple shifted quadratic function. This leads to the conclusion that, when modelling for stock indices, asymmetry is perhaps best modelled by allowing the minimum of g(·) to shift away from zero, as opposed to by allowing each side of zero to scale differently. To contrast the coefficient functions estimated under these two approaches, Figure <ref>(a) shows the posterior mean estimates of g(·) of both the GJR and SP-GARCH models for the S&P 500 returns.

For the individual stocks, both the DIC and D̅ estimates for the Beta-t and SP-GARCH models are substantially lower than those for the GARCH and GJR models. Unlike the index case, the DIC and D̅ estimates for the GJR model are very similar to those for the GARCH model. This suggests either that asymmetry is of less importance for the chosen stocks, or that the type of asymmetry present in these returns can not be adequately captured by the GJR model.

For Apple, ARM, and Intel, the estimates of D̅ for Beta-t are almost identical to those for SP-GARCH, while Beta-t achieves slightly lower DIC estimates. For Nvidia and SanDisk, SP-GARCH is ranked ahead of Beta-t by both DIC and D̅.

By comparing the plots of the estimated coefficient functions of the SP-GARCH model to those of the Beta-t model, we find that the estimates given by the two models are very similar in their shapes. An example of such comparison is shown in Figure <ref>(b) for the Intel returns. Notice that the two functions are almost identical for ϵ_t-1∈ [-2,  2], where most of the observations are found. Moreover, for large negative innovations, the coefficient function of the SP-GARCH shows a similar “damping" behaviour to that of the Beta-t model. This finding provides evidence, at least for the chosen stocks, supporting the empirical validity of the Beta-t model.

The BMA approach allows SP-GARCH to adapt its complexity according to the complexity of the data. This is evident in the estimated effective number of parameters p_D. For the indices, the p_D estimates for SP-GARCH are all close to 6, which are similar to those for the GJR and Beta-t models; for the individual stocks, the p_D estimates for SP-GARCH are increased slightly to approximately 7, indicating that more complexity is needed to accommodate the more complex dynamics of the data.









§ CONCLUSION


In its original form, the volatility function governing the dynamics of the GARCH model has a very restricted form. However, financial markets may be subjected to dramatic changes, calling for flexible models that can adapt to the changes. The main contribution of this paper is to propose the SP-GARCH model, where the sequence of conditional variances {σ^2_t} follows a functional coefficient autoregressive (FAR) model. Flexibility is achieved by allowing the coefficient function to be an arbitrary smooth function of the the lagged innovation. This smooth function is then approximated using regression spline basis functions based on a certain knot configuration. A strong feature of the proposed model is that the GARCH model and many of its popular extensions are special cases of the SP-GARCH model, differentiated by only the shape of the coefficient function g(·); inference on g(·) allows implicit testing of the simpler specifications. To perform inference, while accounting for the uncertainty in knot configuration, an approach based on Bayesian model averaging (BMA) is employed. To implement the BMA based approach, a carefully designed MCMC algorithm is developed to sample from the joint posterior distribution over the knot configurations and the parameter space. As illustrated by the trace plots, the sampler is able to move rapidly through the space on which the joint posterior distribution is defined. A simulation study shows that, when flexibility is required, the SP-GARCH model is, indeed, able to learn the shape of the coefficient function from the data. On the other hand, when the true data generating process is simple (namely, the original GARCH model), the SP-GARCH model, as a conditional volatility estimator, is as efficient as the simpler parametric models. Applications to ten real financial time series reveal that the coefficient functions for the stock indices are different in structural features from those for the individual stocks. A model comparison study based on the DIC confirms that the SP-GARCH model is, indeed, the most adequate description of the data. Furthermore, the SP-GARCH model is able to adapt its complexity according to the complexity of the data; when applied to real return series, the SP-GARCH model is often as parsimonious as the popular parametric models, which makes it very interpretable.

The SP-GARCH model can be naturally extended to incorporate realised measures of volatility constructed using high-frequency intra-daily data, such as realised variance, realised kernel, and realised range (See, for example, <cit.>, <cit.>, <cit.>, and <cit.>). For example, the volatility function of the SP-Realised-GARCH model can be written as σ^2_t = ω + g(x_t-1) σ^2_t-1, where x_t-1 is a lagged realised measure.



§ ACKNOWLEDGEMENT


We thank Gareth W. Peters and Chris J. Oates for their comments on the manuscript. WYC was supported by the Australian Research Council Centre of Excellence for Mathematical and Statistical Frontiers (ACEMS).


chicago


