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.
91,780 characters · 33 sections · 65 citation commands
Bayesian modeling of time-varying conditional heteroscedasticity
\def\spacingset#1{ {#1}} \spacingset{1}
\if00 \fi
\if10 {
} \fi
Keywords: Autoregressive model, B-splines, Hamiltonian Monte Carlo (HMC), Non-stationary, Posterior contraction, Volatility
For datasets observed over a long period, stationarity turns out to be an oversimplified assumption that ignores systematic deviations of parameters from constancy. Thus time-varying parameter models have been studied extensively in the literature of statistics, economics, and related fields. For example, the financial datasets, due to external factors such as war, terrorist attacks, economic crisis, political events, etc. exhibit deviation from time-constant stationary models. Accounting for such changes is crucial as otherwise time-constant models can lead to incorrect policy implications as pointed out by bai97. Thus functional estimation of unknown parameter curves using time-varying models has become an important research topic today. In this paper, we analyze popular conditional heteroscedastic models such as AutoRegressive Conditional Heteroscedasticity (ARCH) and Generalized ARCH (GARCH) from a Bayesian perspective in a time-varying setup. Before discussing our new contributions in this paper, we provide a brief overview of some previous works in this areas.
In the regression regime, time-varying models have garnered a lot of recent attention; see, for instance, fan99, fan2000, hoover98, huang04, lin01, ramsay05, zhang02 among others. The models show time-heterogeneous relationship between response and predictors. Consider the following two regression models $$\text{Model I: } y_i = x_i^{\mkern-1.5mu\mathsf{T}} \theta_i +e_i, \quad \text{Model II: } y_i = x_i^{\mkern-1.5mu\mathsf{T}}\theta_0 + e_i, \quad\quad i = 1, \ldots ,n,$$
where $x_i \in \mathbb{R}^d$ ($i = 1,\ldots,n$) are the covariates, $^{\mkern-1.5mu\mathsf{T}}$ is the transpose, $\theta_0$ and $\theta_i = \theta(i/n)$ are the regression coefficients. Here, $\theta_0$ is a constant parameter and $\theta : [0, 1] \to \mathbb{R}^d$ is a smooth function. Estimation of $\theta(\cdot)$ has been considered by hoover98, cai07) and zhouwu10 among others. One popular way to decide if there is an evidence to favor time-varying models over the time-constant analogue is to perform hypothesis testing. See, for instance, regression2012, regression2015, chow60, brown75, nabeya88, leybourne89, nyblom89, ploberger89, andrews93 and lin99. zhouwu10 discussed obtaining simultaneous confidence bands (SCB) in model I, i.e. with additive errors. However their treatment is heavily based on the closed-form solution and it does not extend to processes defined by a more general recursion.
For time-varying AR, MA, or ARMA processes, the results from time-varying linear regression can be naturally extended. However, such an extension is not obvious for conditional heteroscedastic (CH hereafter ) models. These are, by the simple definition of evolution is difficult to estimate even in the time-constant case. However, one cannot possibly ignore its usefulness in analyzing and predicting financial datasets. These models (even the simple time-constant ones) have remained primary tools for analyzing and forecasting certain trends for stock market datasets since engle82 introduced the classical ARCH model and bollerslev extended it to a more general GARCH model. However, with the rapid dynamics of market vulnerability, the simple classical time-constant models fail in terms of estimation or prediction due to over-compensating the past data. Several references point out the necessity of extending these classical models to a set-up where the parameters can change across time, for example stuaricua2005, engle2005spline and fry08. Consider the simple tvARCH(1) model $$X_i= \sigma_i\zeta_i, \zeta_i \sim N(0,1), \sigma_i^2= \mu_0(i/n)+a_1(i/n)X_{i-1}^2.$$ Similar models can be defined for tvGARCH(1,1) as well where $\sigma_i^2$ has an additional recursive term involving $\sigma_{i-1}^2$ $$X_i= \sigma_i\zeta_i, \zeta_i \sim N(0,1), \sigma_i^2= \mu_0(i/n)+a_1(i/n)X_{i-1}^2+b_1(i/n)\sigma_{i-1}^2.$$ When the two recursive parameters in a GARCH model sum up to 1, i.e. $a_1+b_1=1$ it is usually called an integrated GARCH (iGARCH; or bubble garch/explosive garch by some authors) process which employing the above display can also be extended towards a time-varying analog i.e. $b_1(\cdot)=1-a_1(\cdot)$. A wide range of financial datasets exhibits iGARCH phenomena.
In the parlance of time-varying parameter models in the CH setting, numerous works discussed the CUSUM-type procedure, for instance, kim00 for testing change in parameters of GARCH(1,1). kulperger05 studied the high moment partial sum process based on residuals and applied it to residual CUSUM tests in GARCH models. Interested readers can find some more change–point detection results in the context of CH models in chu95, chen97, lin1999testing, kokoszka00 or andreou06.
A time-varying framework and a pointwise curve estimation using M-estimators for locally stationary ARCH models were provided by dahlhaussubbarao2006. Since then, while several pointwise approaches were discussed in the tvARMA and tvARCH case (cf. dahlhaus2009, dahlhaussubbarao2006, fry08), pointwise theoretical results for estimation in tvGARCH processes were discussed in tvgarch2013 and rohan13 for GARCH(1,1) and GARCH($p$,$q$) models respectively. In a series of recent works sayar2020, sayarthesis such models were discussed in wide generality. However, the focus remained frequentist, and the main goal accomplished there was to build simultaneous inference. One strong criticism for the CH type models remained that one needs a relatively large sample size ($n \sim 2000$) to achieve nominal coverage levels. The recursive definition of the models and a subsequent kernel-based method of estimating make it difficult to achieve satisfying results for relatively smaller sample sizes. This motivated us to explore a Bayesian way of building and estimating these models and use the posteriors to construct posterior estimates of the coefficient curves $\theta(\cdot)$.
In this paper, we develop a Bayesian estimation method for time-varying analogs of ARCH, GARCH, and iGARCH models. We model the time-varying functional parameters using cubic B-splines. In the context of general varying-coefficient modeling, spline bases are a popular choice for its convenience and flexibility hastie1993varying, gu1993smoothing, cai2000functional,biller2001Bayesian, huang2002varying, huang2004functional,amorim2008regression,fan2008statistical,yue2014Bayesian,franco2019unified. Specific to the literature of time-varying volatility modeling, B-spline-based models have also been explored engle2008spline,audrino2009splines,liu2016spline.
Our contributions in this paper are two-fold. Towards the methodological development, note that the tvARCH, tvGARCH, and tviGARCH models require complex shape constraints on the coefficient functions. We achieve those by imposing different hierarchical structures on B-spline coefficients. The constraints are designed to be able to develop an efficient sampling algorithm based on gradient-based Hamiltonian Monte Carlo (HMC) neal2011mcmc, betancourt2015hamiltonian, betancourt2017conceptual,livingstone2019geometric. Strong motivation towards implementing such a Bayesian methodology was to circumvent the requirement of a huge sample size which is almost essential for effective estimation using the frequentist and kernel-based methods. This requirement on sample size has been frequently pointed out in the literature of ARCH/GARCH models and thus this was one of our main motivations to see if a reasonable estimation scheme can be designed in a Bayesian way.
Secondly, the existing literature on obtaining posterior concentration rates for dependent data is thin, even for an extremely simple model. To the best of our knowledge, ours is the first such attempt towards a theoretical development for these models under Gaussian-link. Posterior contraction rates for these models with respect to the average Hellinger metric are established. The main challenge therein is to construct exponentially consistent tests for these classes of models. Using some recently developed tools from jeong2019frequentist, ning2020Bayesian we have developed such tests. We first establish posterior contraction rates with respect to average log-affinity and then the same rate is transferred to the average Hellinger metric. The frequentist literature on inference about time-varying needs very stringent moment assumption and local stationarity assumptions which are often difficult to verify. Moreover, for econometric datasets, the existence of even the fourth moment is often questionable. Thus this paper offers some alternative way to estimate coefficients under lesser assumptions.
The rest of the paper is organized as follows. Section (ref) describes the proposed Bayesian model in detail. Section (ref) discusses an efficient computational scheme for the proposed method. We calculate posterior contraction rate in Section (ref). In Section (ref) we study the performance of our proposed method in the light of. Section (ref) deals with real data application of the proposed methods for the three separate models and concludes with a brief interpretation of the results. We wrap the paper up with discussions, some concluding remarks, and possible future directions in Section (ref). The supplementary materials contain theoretical proofs and some additional results.
We elaborate on the models and our Bayesian framework for time-varying analogs of three specific cases that are popularly used to analyze econometric datasets.
Let $\{X_i\}$ satisfy the following time-varying ARCH($p$) model for $X_i$ given $\mathcal{F}_{i-1}=\{X_j: j\leq (i-1)\}$,
where the parameter functions $\mu(\cdot), a_i(\cdot)$ satisfy
In a Bayesian regime we put priors on $\mu(\cdot)$ and $a_i(\cdot)$. To respect the shape-constraints as imposed by $\mathcal{P}$ we reformulate the problem. With $B_j$ as the B-spline basis functions, let
The prior induced by above construction is $\mathcal{P}$-supported. The verification is very straightforward. In above construction, $\sum_{j=0}^PM_j=1$. Thus $\sum_{j=1}^PM_j\leq 1$. Since $0\leq\theta_{kj}\leq 1$, $\sup_{x}a_i(x)\leq M_i$. Thus $\sup_x\sum_{i=1}^Pa_{i}(x)\leq \sum_{i=1}^PM_i\leq 1$. We have $\sum_{j=1}^PM_j\leq 1$ if and only if $\delta_0=-\infty$, which has probability zero. On the other hand, we also have $\mu(\cdot)\geq 0$ as we have $\exp(\beta_j)\geq 0$. Thus, the induced priors, described above are well supported in $\mathcal{P}$.
Let $\{X_i\}$ satisfy the following time-varying GARCH($p$,$q)$ model for $X_i$ given $\mathcal{F}_{i-1}=\{X_j: j\leq (i-1)\}$,
Additionally we impose the following constraints on parameter space for the time-varying parameters,
The condition on the AR parameters imposed by ((ref)) is somewhat popular in time-varying AR literature. See dahlhaussubbarao2006, subbarao2008, sayar2020 for example. Different from these references, we additionally do not assume existence of any unobserved local-stationary process that are close to the observed process.
To proceed with Bayesian computation, we again put priors on the unknown functions $\mu(\cdot), a_i(\cdot)$ and $b_j(\cdot)$'s such that they are supported in $\mathcal{P}_1$. Again the restrictions imposed by ((ref)) are respected. The complete description of prior is
Here $B_{j}$'s are the B-spline basis functions. The parameters $\delta_{j}$'s are unbounded. The verification of support condition (ref) for the proposed prior is similar.
Although the GARCH(1,1) remains one of the most popular models to analyze econometric datasets, empirical evidence shows that these datasets regularly raise suspicion to the parameter space restriction $\sum_i a_i+\sum_j b_j<1$. Note that we used a time-varying analog of this restriction for the tvGARCH modeling in Section (ref). This often creates a very stringent condition as the validity of $\sum_i a_i(t)+\sum_j b_j(t)<1$ is questionable. The special case for a time-constant GARCH model where this restriction fails is called an iGARCH model in the literature. We consider the following time-varying analog of iGARCH.
We impose the following constraints on parameter space for the time-varying parameters,
The prior functions that allow us to reformulate the problem keeping it consistent with ((ref)) is described below:
The complete likelihood $L$ of the proposed Bayesian method is given by
where $\mu(x) =\sum_{j=1}^{K_1}\exp(\beta_j)B_j(x), a_{k}(x)=\sum_{j=1}^{K_2}\theta_{kj}M_{k}B_j(x)$ and $M_j=\frac{\exp(\delta_j)}{\sum_{k=0}^p{\exp(\delta_k)}}$. We develop efficient Markov Chain Monte Carlo (MCMC) algorithm to sample the parameter $\beta,\theta$ and $\delta$ from the above likelihood. The computation of derivatives allows us to develop an efficient gradient-based MCMC algorithm to sample these parameters. We calculate the gradients of negative log-likelihood $(-\log L)$ with respect to the parameters $\beta$, $\theta$ and $\delta$. The gradients are given below,
where ${\mathbf 1}_{\{j=k\}}$ stands for the indicator function which takes the value 1 when $j=k$.
The complete likelihood $L_2$ of the proposed Bayesian method of (ref) is given by
We calculate the gradients of negative log-likelihood $(-\log L_2)$ with respect to the parameters $\beta$, $\theta$, $\eta$ and $\delta$. The gradients are given below,
While fitting tvGARCH$(p,q)$, we assume for any $t<0$, $X^2_t=0,\sigma^2_t=0$. Thus, we need to additionally estimate the parameter $\sigma^2_0$. The derivative of the likelihood concerning $\sigma^2_0$ is calculated numerically using the {\tt jacobian} function from R package {\tt pracma}. For the tviGARCH, the derivatives are similar so we avoid computing them for the sake of brevity.
Based on these gradient functions, we develop gradient-based Hamiltonian Monte Carlo (HMC) sampling. Note that, parameter spaces of $\theta_{kj}$'s have bounded support. We circumvent this by mapping any Metropolis candidate falling outside the parameter space back to the nearest boundary. HMC has two parameters, required to be specified. These are the leap-frog step and the step-size parameter. It is difficult to tune both of them simultaneously. We choose to tune the step size parameter to maintain an acceptance range between 0.6 to 0.8. After every 100 iterations, the step-length is adjusted (increased or reduced) accordingly if it falls outside. neal2011mcmc showed that a higher leapfrog step is better for estimation accuracy at the expense of greater computation. To maintain a balance between accuracy and computational complexity, we keep it fixed at 30 and obtain good results.
We now focus on obtaining posterior contraction rates for our proposed Bayesian models. The posterior consistency is studied in the asymptotic regime of increasing number of time points $n$. We study the posterior consistency with respect to the average Hellinger distance on the coefficient functions which is $$d_{1,n}^2=\frac{1}{n}d^2_{H}(\kappa_1,\kappa_2)=\frac{1}{n}\int(\sqrt{f_1}-\sqrt{f_2})^2,$$ where $f_1=\prod_{i=1}^nP_{\kappa_1}(X_{i}|X_{i-1})$ and $f_2$ denotes the corresponding likelihoods.
Definition: For a sequence $\epsilon_n$ if $\Pi_n(d(f,f_0)|X^{(n)}\geq M_n\epsilon_n|X^{(n)})\rightarrow 0$ in $F_{\kappa_0}^{(n)}$-probability for every sequence $M_n\rightarrow\infty$, then the sequence $\epsilon_n$ is called the posterior contraction rate.
All the proofs are postponed to the supplementary materials. The proof is based on the general contraction rate result for independent non-i.i.d. observations ghosal2017fundamentals and some results on B-splines based finite random series. The exponentially consistent tests are constructed leveraging on the famous Neyman-Pearson Lemma as in ning2020Bayesian. Thus the first step is to calculate posterior contraction rate with respect to average log-affinity $r_n^2(f_1,f_2)=-\frac{1}{n}\log\int f_1^{1/2}f_2^{1/2}$. Then we show that $r_n^2(f_1,f_2)\lesssim\epsilon_n^2$ implies $\frac{1}{n}d_H^2(f_1,f_2) \lesssim \epsilon_n^2$. We also consider following simplified priors for $\alpha_j$ and $\tau_i$ to get better control over tail probabilities,
Let $\kappa=(\mu, a_1)$ stands for the complete set of parameters. For sake of generality of the method, we put a prior on $K_1$ and $K_2$ with probability mass function given by,
for $i=1,2$. These priors have not been considered while fitting the model as it would require computationally expensive reversible jump MCMC strategy. The contraction rate will depend on the smoothness of true coefficient functions $\mu$ and $a$ and the parameters $b_{13}$ and $b_{23}$ from the prior distributions of $K_1$ and $K_2$. Let $\kappa_0=(\mu_0, a_{01})$ be the truth of $\kappa$.
Assumptions(A): There exists constants $M_X>1, 0<M_{\mu}<M_X$ such that,
Assumptions (A.1)-(A.3) ensure $$\mathbb{E}_{\kappa_0}(X_i^2)=\mathbb{E}_{\kappa_0}(\mathbb{E}_{\kappa_0}(X_i^2|X_{i-1}))<M_{\mu}+\left(1-\frac{M_{\mu}}{M_X}\right)M_X<M_X$$ by recursion.
where $b_{ij}$ are specified in ((ref)).
Let $\kappa=(\mu, a_1, b_1)$ stands for the complete set of parameters. For sake of generality of the method, we put a prior on $K_1$, $K_2$ and $K_3$ with probability mass function given by,
for $i=1,2$. These priors have not been considered while fitting the model as it would require computationally expensive reversible jump MCMC strategy. The contraction rate will depend on the smoothness of true coefficient functions $\mu$ and $a$ and the parameters $b_{13}$ and $b_{23}$ from the prior distributions of $K_1$ and $K_2$. Let $\kappa_0=(\mu_0, a_{01})$ be the truth of $\kappa$.
Assumptions(B): There exists constants $M_X>1, 0<M_{\mu}<M_X$ such that,
Assumptions (B.1) and (B.3) ensure $$\mathbb{E}_{\kappa_0}(X_i^2)=\mathbb{E}_{\kappa_0}(\mathbb{E}_{\kappa_0}(X_i^2|X_{i-1}))<M_{\mu}+\left(1-\frac{M_{\mu}}{M_X}\right)M_X<M_X$$ by recursion. Similarly we have $\mathbb{E}(\sigma_i^2)=\mathbb{E}_{\kappa_0}(\mathbb{E}_{\kappa_0}(X_i^2|\mathcal{F}_i))=\mathbb{E}_{\kappa_0}(X_i^2)<M_X$.
where $b_{ij}$ are specified in ((ref)).
Let $\kappa=(\mu, a_1)$ stands for the complete set of parameters. For sake of generality of the method, we put a prior on $K_1$ and $K_2$ with probability mass function given by,
for $i=1,2$. These priors have not been considered while fitting the model as it would require computationally expensive reversible jump MCMC strategy. The contraction rate will depend on the smoothness of true coefficient functions $\mu$ and $a$ and the parameters $b_{13}$ and $b_{23}$ from the prior distributions of $K_1$ and $K_2$. Let $\kappa_0=(\mu_0, a_{01})$ be the truth of $\kappa$.
where $b_{ij}$ are specified in ((ref)).
We run simulations to study the performance of our proposed Bayesian method in capturing the true coefficient functions under different true models. The hyperparameters $c_1$ and $c_2$ of the normal prior are all set 100, which makes the prior weakly informative. We consider 4, 5 and 6 equidistant knots for the B-splines when $n=200, 500$ and $1000$ respectively. We collect 10000 MCMC samples and consider the last 5000 as post burn-in samples for inferences. We shall compare the estimated functions with the true functions in terms of the posterior estimates of functions along with its 95% pointwise credible bands. The credible bands are calculated from the MCMC samples at each point $t=1/T,2/T,\ldots, 1$. We take the posterior mean as the posterior estimate of the unknown functions.
Since, to the best of our knowledge, there is no other Bayesian model for these time-varying conditional heteroscedastic models, we compare our Bayesian estimates with corresponding frequentist time-varying estimates. For computing the time-varying estimates of these models, we use the kernel-based method from sayar2020. The M-estimator of the parameter vector $\theta(t)$ are obtained using the conditional quasi log-likelihood. For instance, in the tvARCH(1) case, say $\theta(t)=(\mu(t),a_1(t))$
where $\ell(\cdot)$ denotes the Gaussian log-likelihood. Note that these methods are fast but usually need a cross-validated choice of bandwidth $b_n$. We use $K(x)=3/4(1-x^2)\mathbf{I}(|x| \leq 1)$ and an appropriately chosen bandwidth as suggested by the authors therein. Since our discussion also involves iGARCH formulation, we wrote a separate kernel-based frequentist estimation for iGARCH models analogously. Apart from these two time-varying estimates, we also obtain a time-constant fit on the same data to help initiate a discussion on whether there was a necessity of introducing coefficients varying with time. For this tseries and rugarch R packages are used respectively for ARCH/GARCH and iGARCH fits.
To compare these estimates, we evaluate the average mean square errors (AMSE) for the three estimates. Note that in an usual linear regression of response $y$ on predictor $X$ scenario, the fitted MSE is often defined as $\frac{1}{n}\sum (y_i - \hat{y}_i)^2$. Since here, $X_i|\mathcal{F}_{i-1} \sim N(0,\sigma_i^2)$, we use the following definition of AMSE $$\text{AMSE}=\frac{1}{n}\sum_i(X_i^2-\hat{\sigma}_i^2)^2$$ where the $\hat{\sigma}_i^2$ is computed with the fitted parameter values as per the model under consideration. For example, for a tvGARCH(1,1) model we have $$\hat{\sigma}_i^2= \hat{\mu}(i/n)+\hat{a}(i/n)X_{i-1}^2+\hat{b}(i/n)\hat{\sigma}_{i-1}^2$$ where $\hat{\mu}(\cdot), \hat{a}(\cdot)$ and $\hat{b}(\cdot)$ are the estimated curves from the posterior. Replacing the response $y_i$ by $X_i^2$ is natural as often autocorrelations of $X_i^2$ are checked to gauge presence of CH effect. Moreover, one of the early methods to deal with CH models was to view $X_i^2$ approximated by an TVAR(1) process. See bose2003estimating and references therein. Similar estimators as our proposed AMSE to evaluate the fitting accuracy has been used in the literature previously. See starica2003garch, fry08, tvgarch2013, sayar2020 for example.
In the next three subsections, we provide the results for the three models, namely, tvARCH, tvGARCH, and tviGARCH. Our conclusions from these results are illustrated at the end of the section.
We start by considering the following tvARCH(1) model from (ref). Three different choices for $n$ are considered, $n=200,500$ and $1000$. The true functions are,
We compare the estimated functions with the truth for sample size 1000 in Figures (ref). Table (ref) illustrates the performance of our method with respect to other competing methods.
Next we explore the following GARCH(1,1) model (cf. (ref))for different choices of $n$. The true functions are,
Note that, estimation of GARCH, due to the additional $b_i(\cdot)$ parameter curves is a significantly more challenging problem and often requires a much higher sample size to have a reasonable estimation. We show by the means of the following pictures in Figure (ref) that the estimation looks reasonable even for smaller sample sizes. The AMSE score comparisons are shown in Tables (ref). The performance of our method is also contrasted with other competing methods.
Finally we consider the tviGARCH(1,1) model (cf. (ref)) a special case of GARCH. Note that due to the constraint $a_1(\cdot)+b_1(\cdot)=1$ we only consider the mean function and AR(1) function for plotting. For this case, our true functions are as follows
The frequentist computation for tviGARCH method is carried out based on a kernel-based estimation scheme along the same line as sayar2020. The estimated plots along with the 95% credible intervals are shown in Figure (ref) for three sample sizes $n=200,500,1000$ and the AMSE scores in Table (ref).
To summarize, our estimated functions are close to true functions for all the cases. We also find that the credible bands are getting tighter with increasing sample size. Thus estimates are improving in precision as sample size increases as shown in Figures (ref) to (ref). AMSEs of our Bayesian estimates are at least better for all the cases as in Tables (ref) to (ref). For tviGARCH, AMSE* is considered due to the huge and somewhat incomparable values of AMSE due to non-existent variance.
Towards applying our methods on real-life datasets we stick to econometric data for varying time horizons. These datasets show considerable time-variation justifying our models to be suitable for understanding how the parameter functions have evolved. Typically we model the log-return data of the daily closing price of these data to avoid the unit-root scenario. The log-return is defined as follows and is close to the relative return
$$Y_i=\log P_i- \log P_{i-1}=\log \left( 1+ \frac{P_i-P_{i-1}}{P_{i-1}}\right) \approx \frac{P_i-P_{i-1}}{P_{i-1}}, $$ where $P_i$ is the closing price on the $i^{th}$ day. Conditional heteroscedastic models are popularly used for model building, analysis and forecasting. Here we extend such models to a more sophisticated and general scenario by allowing the coefficient functions to vary.
In this section, we analyze two datasets: USD to JPY conversion and NASDAQ, a popular US stock market data. We analyze the NASDAQ data through tvGARCH(1,1) and tviGARCH(1,1) models and USDJPY conversion rate data through tvARCH(1) models. We just fit one lag for these models as multiple lag fits are similar and larger lags seem to be insignificant. This result is consistent with the findings in sayar2020, subbarao2008 etc. Moreover, as subbarao2008 claims, stock indices and Forex rates are more suited to GARCH and ARCH type of models respectively for their superior predictive performance. Each of these datasets was collected up to 31 July 2020. We exhibit our results for the last 200,500 and 1000 days which capture the last 6 months, around 1.5 years, and around 3 years of data respectively. All these datasets were collected from \url{www.investing.com}. Note that these datasets are usually available for weekdays barring holidays and typically there are about 260 data points every single year.
We obtain the following Figure (ref) that shows our estimation for fitting a tvARCH(1) model on the USD to JPY conversion data for the last 200,500 and 1000 days ending on 31 July 2020. The AMSE is also computed and contrasted with other competing methods in Table (ref). Figure (ref) depicts the estimated functions with 95% credible bands for different sample sizes. One can see the bands become much shorter for larger sample sizes. The mean coefficient function $\mu(\cdot)$ is generally time-varying for all three cases as one cannot fit a horizontal line through the 95% posterior bands. There seems to be a rise in the mean value around 100 days ago from July 31, 2020, which is right around the time the COVID-19 pandemic hit the world. With the analysis of $n=1000$ days, we see that the volatility is quite high around October 2016 which coincides with the presidential election time of 2016. The AR(1) coefficient does not show the huge time-varying property. We also tabulate the AMSE for the three sample sizes in Table (ref) and one can see for smaller sample sizes such as $n=200$, the proposed Bayesian tvARCH achieves a significantly lower score but when the sample size grows then the performance becomes similar.
As has become standard in analyzing stock market datasets using GARCH models, we use time-varying GARCH for small orders. We obtain the following Figure (ref) for fitting a tvGARCH(1,1) model on the NASDAQ data for the last 200,500, 1000 days ending on 31 July 2020. One can see the $a_1(\cdot)$ values are generally low and the $b_1(\cdot)$ values are higher which is consistent with how these outcomes turn out for time-constant estimates for econometric datasets. One can also see the role sample size plays in curating these time-varying estimates. For $n=200$, the $b_1(\cdot)$ achieves high value of 0.6 around mid-March 2020 but for higher sample sizes it shows values as high as 0.8. One can also note the striking similarity for the analysis of the last 500 and 1000 days which is fairly consistent with the idea that estimation is more stable for such CH type models with a higher sample size. Nonetheless, the estimates for $n=200$ seem quite smooth as well which can be seen as a benefit of our methodology. Table (ref) provides a comparison of AMSE scores across the three methods for three sample sizes. The Bayesian tvGARCH(1,1) performs relatively better than other methods and estimated curves have smaller credible bands with a growing sample size. The behavior of the mean function also shows higher volatility around the pandemic.
In Figure (ref) the sum of estimated coefficient functions $a(\cdot)+b(\cdot)$ is close to 1 for a significant time-horizon. This motivates us to also fit tviGARCH(1,1) to analyze the same NASDAQ data. The estimated functions are presented in Figure (ref) for the last $n=200, 500$ and $1000$ days. Table (ref) compares the AMSE scores for the same three methods as before with varying sample sizes. The estimated mean and AR(1) functions of Figure (ref) change a little from the estimated functions of tvGARCH(1,1) fit in Figure (ref). Moreover, the effect of the three sample sizes is clear here with $n=1000$ showing very precise bands and can reveal an interesting time-varying pattern.
In terms of AMSE, one can see in Table (ref) that the frequentist methods did worse than even the time-constant versions. The time-constant estimates were computed using the rugarch package in R. The Bayesian tviGARCH method provides significantly better AMSE uniform overall sample sizes. Here the mean function also shows higher volatility around the time when the pandemic struck us. Volatility due to the presidential election in 2016 can also be observed here.
For the analysis of NASDAQ data, we have used two different models and thus it is pertinent to answer how should one choose between a competing class of models. We provide some measures in this subsection to decide between these two competing models. We start by comparing the performance of tvGARCH and tviGARCH models in terms of Bayes factor kass1995bayes. Our calculation of the Bayes factor is based on the posterior samples using the harmonic mean identity of neton1994approximate. Let us denote $B_{200}, B_{500}$ and $B_{1000}$ as the Bayes factors for three sample sizes where $$B_{i}=\frac{P(D^{(i)}|\textrm{tvGARCH})}{P(D^{(i)})|\textrm{tviGARCH})}$$ for sample size $i$ and the corresponding dataset $D^{(i)}$. The values we obtain are $2\log(B_{200})=8.16$, $2\log(B_{500})=19.08$ and $2\log(B_{1000})=24.14$. According to guidelines from section 3.2 of kass1995bayes, there is `positive' evidence in favor of tvGARCH for sample sizes 200 and 500. However, the same evidence becomes `strong' for sample size 1000.
We also try to address out-of-sample predictive performance comparison here. Note that for a time-varying GARCH or time-varying iGARCH model this is generally a difficult task due to the assumed non-stationarity of the model. Thus, we take following approach to calculate out of sample joint predictive log-likelihoods for model comparison. Let us assume we have the data $D^{(n)}$ with $n$ data points. To evaluate the joint predictive log-likelihood for the last $m (< n)$ most recent data points, we fit the models in (ref) and (ref) with the first $n-m$ data. Note that the assumed time horizons for these two models are $n$. Based on the estimated B-spline coefficients, and other parameters from each model, we can compute the joint predictive log-likelihood of the last $m$ data points as $$L_m^{(n)}=\frac{1}{m}\sum_{i=n-m+1}^{n}\frac{1}{2}\left\{-X_{i}^2/\hat{\sigma}_{i}^2-\log(\hat{\sigma}_{i})-\log(2\pi)\right\},$$ where $\hat{\sigma}_{i}^2=\hat{\mu}(i/n)+\hat{a}_1(i/n)X_{i-1}^2+\hat{b}_1(i/n)\hat{\sigma}_{i-1}^2$. For the tviGARCH model, we have $\hat{b}_1(\cdot)=1-\hat{a}_{1}(\cdot)$.
Using this predictive log-likelihood we decide to evaluate the two fits from tvGARCH and tviGARCH in the following manner. For each of the sample sizes, we run it on three separate regimes of the data, the full data, and two halves of the data. In all these 8 settings, ( three sample sizes, three possible regions of the data, but the latest half of the 1000-sized data is the same as the full data for sample size 500) we compute 10, 20, and 50 steps ahead forecast. We tabulate these results in Table (ref). One can see that generally speaking, there is somewhat conclusive evidence towards the iGARCH model for a smaller sample size. This supports our motivation why we additionally provide a tviGARCH(1,1) modeling on the same dataset.
Based on our model comparison exercises, we have an interesting phenomenon where for in-sample model fit, tvGARCH is better but in terms of out-of-sample prediction, tviGARCH outperforms tvGARCH in almost all the cases. Note that, tvGARCH has one additional free parameter and thus is expected to fit the data better but since the estimated $a_1(\cdot)$ and $b_1(\cdot)$ coefficients are close to one satisfying the iGARCH formulation, the out-of-sample performance for tviGARCH may have exceeded that for tvGARCH.
As per the suggestion from a reviewer, we also add a one-step-ahead point forecasting exercise between these models. Here the computation method remains the same as outlined in the predictive log-likelihood computation however we only restrict ourselves to $m=1$ to make the discussion concise. For this part of the exercise we choose to compute posterior mean of $(X_n^2-\hat{\sigma}^2)^2$ where to ensure out-of-sample prediction $\hat{\sigma}^2$ is estimated solely based on $X_1,\ldots,X_{n-1}$. As one-step-ahead forecasts can be prohibitively misleading given it depends so much on one single location, we decide to take an average of over 15 random starting points over the entire time spectrum of 10 years resulting in 2518 points). For each of the sample sizes, we tabulate the performance in the following Table (ref). Note that, here we are only comparing the two Bayesian time-varying models to see which one fits our data better. The advantage of predicting the future coefficients using B-spline is not available in the kernel-based frequentist method and thus is not included here in the discussion.
One can see, we again observe the same advantage of tviGARCH modeling over the tvGARCH one for smaller sample sizes. This is an interesting find of this paper in the context of the Bayesian model fitting to these datasets.
In this paper, we consider a Bayesian estimation framework for time-varying conditional heteroscedastic models. Our prior specifications are amenable to Hamiltonian Monte Carlo for efficient computation. One of the key motivations towards going to a Bayesian regime was to achieve reasonable estimation even for a small sample size. Our simulation coverage shows good performance for all three models tvARCH, tvGARCH, tviGARCH for both small and large sample sizes. Importantly, in all three of the cases, we were able to establish posterior contraction rates. These calculations are, to the best of our knowledge, the first such work in even simple dependent models let alone the complicated recursions that these conditional heteroscedastic models demand. Moreover, the assumptions on the true functions and the number of moments needed were very minimal. An interesting future theoretical work would be to calculate posterior contraction rate with respect to empirical $\ell_2$-distance which is a more desirable metric for function estimation. While analyzing real data, we see that the parameter curves vary significantly for the intercept terms, but not that much for AR or CH parameters. The associated codes to fit the three models are available at \url{https://github.com/royarkaprava/tvARCH-tvGARCH-tviGARCH}.
As future work, it will be interesting to explore multivariate time-varying volatility modeling multigarch,multigarch2 through a Bayesian framework similar to ours. Another interesting time-heterogeneity that we plan to explore through the glass of a Bayesian framework is regime-switching CH models where instead of the smooth time-varying functions the parameters change abruptly. We have a brief discussion in section (ref) about how to choose between competing models. Those discussions can easily be extended to choose a proper number of lags or to choose between different types of ARCH/GARCH models. We believe this would provide an interesting parallel to the usual penalized likelihood-based methods for model selection in time-series. Finally note that, even though we do provide some insights onto future prediction for these datasets for real data applications, that was not the main focus in this paper. Forecasting for the time-varying model is extremely tricky since it requires `estimation' of the future parameter values. Although in-filled asymptotics can help in this regard, still the literature so far is very sparse in this direction for both Bayesian and frequentist regimes. We plan to explore this extensively in near future.
We would like to thank the editor, the associate editor, and anonymous referees for their constructive suggestions that improved the quality of the manuscript.
We study the frequentist property of the posterior distribution in increasing $n$ regime assuming that the observations are coming from a true density $f_0$ characterized by the parameter $\kappa_0$. We follow the general theory of ghosal2000convergence to study posterior contraction rate for our problem. In Bayesian framework, the density $f$ is itself a random measure and has distribution $\Pi$ which is the prior distribution induced by the assumed prior distribution on $\kappa$. The posterior distribution of a neighborhood $U_n=\{f:d(f,f_0)<\epsilon_n\}$ around $f_0$ given the observation $X^{(n)}=\{X_0,X_1,\ldots,X_n\}$ is $$ \Pi_n(U_n^c|X^{(n)})=\frac{\int_{U_n^c}f(X^{(n)})d \Pi(\kappa)}{\int f(X^{(n)})d\Pi(\kappa)} $$
The posterior consistency would hold if above posterior probability almost surely goes to zero in $F_{\kappa_0}^{(n)}$ probability as $n$ goes to $\infty$, where $F_{\kappa_0}^{(n)}$ is the true distribution of $X^{(n)}$. Recall the definition of posterior contraction rate; for a sequence $\epsilon_n$ if $\Pi_n(d(f,f_0)|X^{(n)}\geq M_n\epsilon_n|X^{(n)})\rightarrow 0$ in $F_{\kappa_0}^{(n)}$-probability for every sequence $M_n\rightarrow\infty$, then the sequence $\epsilon_n$ is called the posterior contraction rate. If the assertion is true for a constant $M_n=M$, then the corresponding contraction rate becomes slightly stronger.
Note that for two densities $f_{0}, f$ characterized by $\kappa_0$ and $\kappa$ respectively, the Kullback-Leibler divergences are given by
Assume that there exists a sieve in parameter space such that $\Pi(W_n^c)\leq \exp(-(C_n+2)n\epsilon_n^2)$ and we have tests $\chi_n$ such that $$\mathbb{E}_{\kappa_0}(\chi_n)\leq e^{-L_nn\epsilon_n^2/2}\quad \sup_{\kappa\in W_n: d^2(f,f_0)>L_{n}\epsilon_n^2}\mathbb{E}_{\kappa}(1-\chi_n)\lesssim e^{-L_{n}n\epsilon_n^2}$$ for some $L_{n}>C_n+2$. Say $U_n=\{f:d^2(f,f_0)\leq L_n\epsilon_n^2\}$ and $S_n=\{\int \frac{f(X^n)}{f_0(X^n)}d\Pi(\kappa)\geq\Pi_{n}(\frac{1}{n}KL(\kappa_0,\kappa)<\epsilon_n)\exp(-C_{n}n\epsilon_n^2)\}$. We can bound the posterior probability from above by,
Taking expectation with respect to $\kappa_0$, first term goes to zero by construction of $\chi_n$. The second term $\mathbb{E}_{\kappa_0}\mathbb{1}\{S_n^c\}$ goes to zero due to Lemma 8.21 of ghosal2017fundamentals for any sequence $C_n\rightarrow \infty$. We would require that $\Pi_{n}\{\frac{1}{n}KL(\kappa_0,\kappa)< \epsilon_n\} \geq \exp(-n\epsilon_n^2)$. Then for the third term,
Thus we adhere to the following plan
We first study the contraction properties with respect to $d^2(f,f_0)=r_n^2(f, f_0)=-\frac{1}{n}\log\int\sqrt{ff_0}$ and then show that the same rate holds for average Hellinger $\frac{1}{n}d_H^2(f,f_0)$. Note that $L_n$ can be taken as $L_n=M_n^2$.
For the sake of technical convenience we show our proof by fixing lag order $p$ at $p=1$, however the results are easily generalizable for any fixed order $p$. All the proofs go through for higher lags with the same technical tools.
The likelihood based on the parameter space $\kappa$ is given by $ \mathbb{P}_{\kappa}(X_0)\prod_{i=1}^n\mathbb{P}_{\kappa}(X_i|X_{i-1}). $ Let $Q_{\kappa, t}(X_i)$ be the distribution of $X_i$ with parameter space $\kappa$.
We have
where $\sigma_i^2=\mu(i/n)+a_1(i/n)X_{i-1}^2$ and $\sigma_{0i}^2=\mu_0(i/n)+a_{01}(i/n)X_{i-1}^2$. Then
Our first goal is to estimate these two quantities in terms of the distances between $\mu(\cdot), \mu_0(\cdot)$ and $a_1(\cdot), a_{01}(\cdot)$. In the light of $\mathbb{E}_{\kappa_0}(I+II)=\mathbb{E}_{\kappa_0}(I)+\mathbb{E}_{\kappa_0}(II)$, we bound the expectation of $I$ and $II$ separately.
Bounding the first term I: Note that, by the means of mean value theorem, for a random variable $\sigma_{*t}^2$ between $\sigma_i^2$ and $\sigma_{0i}^2$,
where $\sigma_{*i}^2>\rho$ for the first term and $\sigma_{*i}^2>\rho X_{i-1}^2$ is used for the second term due to the assumption (A.2). Thus the first term $I$ satisfies
$$\mathbb{E}_{\kappa_0}(I) \leq \mathbb{E}_{\kappa_0}(|I|) \lesssim T(\|\mu-\mu_0\|_{\infty} +\|a_1-a_{01}\|_{\infty}),$$ in the light of assumption (A.3).
Bounding the second term II: For the second term we proceed as follows: $$X_i^2(1/\sigma_{0i}^2- 1/\sigma_i^2)\} \leq \frac{1}{\sigma_i^2\sigma_{0i}^2}|X_i^2\sigma_{0i}^2-X_i^2\sigma_i^2|\leq \frac{X_i^2}{\sigma_{0i}^2}\left(\frac{|\mu(i/n)-\mu_0(i/n)|}{\mu(i/n)}+\frac{|a_1(i/n)- a_{01}(i/n)|X_{i-1}^2}{a(i/n)X_{i-1}^2}\right)$$
where we use the fact that $\mu(i/n)>\rho$ and $a(i/n)>\rho$ due to closeness of $(\mu,a)$ and $(\mu_0,a_0)$. Consequently, we have the deterministic inequality
$$X_i^2(1/\sigma_{0i}^2- 1/\sigma_i^2)\}\leq \frac{X_i^2}{\sigma_{0i}^2}\frac{\|\mu-\mu_0\|_{\infty}+\|a_1-a_{01}\|_{\infty}}{\rho}.$$
Taking expectation under truth we arrive at similar upper bound for $\mathbb{E}_{\kappa_0}(II)$ as $\mathbb{E}_{\kappa_0}(I)$ since $$\mathbb{E}(X_i^2/\sigma_{0i}^2)=\mathbb{E}(\mathbb{E}(X_i^2/\sigma_{0i}^2|X_{i-1}))=1.$$
Thus we conclude that,
In this section, we focus on the requirements to calculate posterior contraction rate as in Section (ref). We first show posterior consistency in terms of average negative log-affinity which is defined as $r_n^2(f_1,f_2)=-\frac{1}{n}\log\int f_1^{1/2}f_2^{1/2}$ between $f_1$ and $f_2$. Here, we have $f_1=\prod_{i=1}^n P_{\kappa_1}(X_{i}|X_{i-1})$. Then we show that, having $r_n^2(f_1,f_0)\lesssim\epsilon_n^2$ implies that our distance metric $d_{2,n}^2(f_1,f_0)\lesssim \epsilon_n^2$.
Proceeding with the rest of the proof of Theorem 1, we use the results of B-Splines, $\|\mu-\mu_0\|_{\infty}\leq \|\alpha-\alpha_0\|_{\infty},$ where $\alpha = \{\alpha_j\}$ and $\|a_1-a_{01}\|_{\infty}\leq \|\gamma-\gamma_{0}\|_{\infty}$, where $\gamma_{j}=\theta_{1j}M_1,$ such that $\gamma_{j}<1$. The H\"older smooth functions with regularity $\iota$ can be approximated uniformly up to order $K^{-\iota}$ with $K$ many B-splines. Thus we have $\epsilon_n\gtrsim\max\{K_{1n}^{-\iota_1},K_{2n}^{-\iota_2}\}$.
We start by providing a lower bound of the prior probability as required by (i) in Plan (ref). We also have the result (ref) and the prior probabilities $\Pi(\|\alpha-\alpha_0\|_{\infty} \lesssim \epsilon_n, \|\gamma-\gamma_{0}\|_{\infty} \lesssim \epsilon_n) \gtrsim \epsilon_n^{K_{1n}+K_{2n}}$ based on the discussion of A2 from shen2015adaptive. The rate of contraction cannot be better than the parametric rate $n^{-1/2}$, and so $\log (1/\epsilon_n)\lesssim \log T$. Thus (i) requires that in terms of pre-rate $\bar{\epsilon}_n$ as $(K_{1n}+K_{2n})\log n\lesssim n\bar{\epsilon}_n^2$. Now we construct a sequence of test functions having exponentially decaying Type I and Type II error. In our problem, we consider following sieve
where $A_n, B_n$ are at least polynomial in $n$ and $\sigma_0^2$ is the standard deviation of $X_0^2$ and $K_n=\max\{K_{1n},K_{2n}\}$. We take $\rho_n=O(n^{-a})$ with $a<1$, $A_n=O(n^{a_1}), B_n=O(n^{a_2})$ with $a_2>a_1$ for technical need.
We need to choose these bounds carefully so that we have $\Pi(W_n^c)\leq \exp(-(1+C_1)n\epsilon_n^2)$, which depend on tail properties of the prior. We have,
where $\alpha_{K_{1n}}$ is the vector of full set of coefficients of length $K_{1n}$ and $\gamma_{K_{2n}}$ is the vector of coefficients of length $K_{2n}$. The quantity $\Pi[\alpha_{K_{1n}}\notin[\rho_n), A_n]^{K_{1n}}$ can be further upper bounded by $K_{1n}\Pi(\alpha_{1}\notin[\rho_n, A_n)])\leq K_{1n}\exp\{-R_1n^{a_3}\}$, for some constant $R_1, a_3>0$ which can be verified from the discussion of the assumption A.2 of shen2015adaptive for our choice of exponential prior. On the other hand, $$\Pi\{\gamma_{K_{2n}}\notin\left[\rho_n, 1-\frac{A_n}{B_n}\right]^{K_{2n}}\}\leq K_{2n}\Pi(\gamma_1\notin[\rho_n, 1-\frac{A_n}{B_n}])\leq K_{2n}\exp\{-R_2n^{a_4}\},$$ for some constant $R_2, a_4>0$ which can be verified from the proof of roy2018high. Hence, $\Pi(W_n^c)\lesssim F_1(K_{1n})+F_2(K_{2n})+(K_{1n}+K_{2n})\exp\{-Rn^{a_5}\}$. The two functions $F_1$ and $F_2$ in the last expression stand for the tail probabilities of the prior of $K_1$ and $K_2$. We can calculate their asymptotic order as, $F_1(x)=\Pi(K_{1}>x)\asymp\exp\{-x(\log x)^{b_{13}}\}$ and $F_2(x)=\Pi(K_{2}>x)\asymp\exp\{-x(\log x)^{b_{23}}\}$. We need $\Pi(W_n^c)\lesssim\exp\{-(1+C_n)n\epsilon_n^2\}$. Hence, we calculate pre-rate from the following equation for some sequence $H_n\rightarrow\infty$,
Now, we construct test $\chi_n$ such that $$\mathbb{E}_{\kappa_0}(\chi_n)\leq e^{-L_nn\epsilon_n^2/2}\quad \sup_{\kappa\in W_n: r_n^2(\kappa,\kappa_0)>L_n\epsilon_n^2}\mathbb{E}_{\kappa}(1-\chi_n)\lesssim e^{-L_nn\epsilon_n^2}$$ for some $L_n>C_n+2$. To construct the test, we first construct the for point alternative $H_0:\kappa=\kappa_0$ vs $H_1:\kappa=\kappa_1$. The most powerful test for such problem is Neyman-Pearson test $\phi_{1n}=\mathbb{1}\{f_1/f_0\geq 1\}$. For $r_n^2> L_n\epsilon_n^2$, we have
$$\mathbb{E}_{\kappa_0}\phi_{1n}=\mathbb{E}_{\kappa_0}(\sqrt{f_1/f_0}\geq 1)\leq \int \sqrt{f_1f_0}\leq \exp(-L_nn\epsilon_n^2),$$
$$\mathbb{E}_{\kappa_1}(1-\phi_{1n})=\mathbb{E}_{\kappa_1}(\sqrt{f_0/f_1}\geq 1)\leq \int \sqrt{f_0f_1}\leq \exp(-L_nn\epsilon_n^2).$$
It is natural to have a neighborhood around $\kappa_1$ such the Type II error remains exponentially small for all the alternatives in that neighborhood under the test function $\phi_{1n}$. By Cauchy-Schwarz inequality, we can write that $$\mathbb{E}_{\kappa}(1-\phi_{1n})\leq \{\mathbb{E}_{\kappa_1}(1-\phi_{1n})\}^{1/2}\{\mathbb{E}_{\kappa_1}(f/f_1)^2\}^{1/2}.$$In the above expression, the first factor is already exponentially decaying. The second factor can be allowed to grow at most of order $e^{cn\epsilon_n^2}$ for some positive small constant $c$. We, in fact establish that the second factor is bounded by our choice of the sieve through the following claim.
The test function $\chi_n$ satisfying exponentially decaying Type I and Type II probabilities is then obtained by taking maximum over all tests $\phi_{jn}$'s for each ball, having above radius. Take $\chi_n=\max_j\phi_{jn}$. Type I and Type II probabilities are given by $\mathbb{P}_0(\chi_n)\leq \sum_{j}\mathbb{P}_0\phi_{jn}\leq D_{n}\mathbb{P}_0\phi_{jn}$ and $\sup_{\kappa\in W_n: r_n(\kappa,\kappa_0)>Mn\epsilon_n^2}\mathbb{P}(1-\chi_n)\leq \exp(-cn\epsilon_n^2)$. Hence, we need to show that $\log D_{n}\lesssim n\epsilon_n^2$, where $D_n$ is the required number of balls of above radius needed to cover our sieve $W_n$. We have
where $r_1,r_2$ are defined in ((ref)). Given our choices of $A_n, B_n$ and $\rho_n$, the two radii $r_1$ and $r_2$ are some fractional polynomials in $n$. Thus $\log D_n\lesssim (K_{1n}+K_{2n})\log T$ which is required to be $\lesssim n\epsilon_n^2$ as in the prior mass condition. Based on (ref), we have $$K_{1n}\asymp n^{1/(2\iota_1+1)}(\log n)^{\iota/(2\iota+1)+(1-b_{13})/2}, K_{2n}\asymp n^{1/(2\iota_2+1)}(\log n)^{\iota/(2\iota_2+1)+1-b_{23}}$$ and a pre-rate $$\bar{\epsilon}_n=\max\bigg\{n^{-\iota_1/(2\iota_1+1)} (\log n)^{\iota_1/(2\iota_1+1)},n^{-\iota_2/(2\iota_2+1)} (\log n)^{\iota_2/(2\iota_2+1)}\bigg\}.$$ The actual rate will be slower that pre-rate. Now, the covering number condition, prior mass conditions and basis approximation give us $(K_{1n}+K_{2n})\log n\lesssim n\epsilon_n^2$ and $\epsilon_n\gtrsim\max\{K_{1n}^{-\iota_1},K_{2n}^{-\iota_2}\}$. Combining all these conditions, we calculate the posterior contraction rate $\epsilon_n$ equal to
Write Reyni divergence as $$r_n^2=-\frac{1}{n}\log\int\sqrt{f_0f_1}=-\frac{1}{n}\log \mathbb{E}_{\kappa_0}\sqrt{\frac{f_1}{f_0}}.$$We need to show $r_n^2\lesssim\epsilon_n^2$ implies that $d_{2,n}^2(\kappa_0,\kappa)\lesssim \epsilon_n^2$ as $\epsilon_n$ goes to zero. If $r_n^2\leq\epsilon_n^2$, we have $\left(\mathbb{E}_{\kappa_0}\sqrt{\frac{f_1}{f_0}}\right)^{-1/n}\leq\exp(\epsilon_n^2)$ which implies for small $\epsilon_n^2$, we have $\left(\mathbb{E}_{\kappa_0}\sqrt{\frac{f_1}{f_0}}\right)^{1/n}\geq 1-\epsilon_n^2$. By Cauchy-Squarz inequality $\left(\int\sqrt{f_0f_1}\right)^2\leq \int f_0\int f = 1$. Thus we have,
$$ 1-\epsilon_n^2\leq \left(\mathbb{E}_{\kappa_0}\sqrt{\frac{f_1}{f_0}}\right)^{1/n}\leq 1, $$
Since $d_{H}^2(f_1,f_0)=2(1-\mathbb{E}_{\kappa_0}\sqrt{\frac{f_1}{f_0}})$ $$\left(\mathbb{E}_{\kappa_0}\sqrt{\frac{f_1}{f_0}}\right)^{1/n}=\left\{1-\left(1-\mathbb{E}_{\kappa_0}\sqrt{\frac{f_1}{f_0}}\right)\right\}^{1/n}\approx1-\frac{1}{2n}d_{H}^2(f_1,f_0).$$ Thus $\frac{1}{n}d_{H}^2(f_1,f_0)\lesssim\epsilon_n^2$. Thus it is consistent under average Hellinger distance.
Note that, the proof of Theorem (ref) follows via exactly same route. We just jot down the important differences from the proof of Theorem (ref). For Theorem (ref) also, we restrict ourselves to tvGARCH(1,1) situation with the assurance that the proof easily extends to a general tvGARCH(p,q) case.
Note that the KL support step from subsubsection (ref) follows almost similarly with the modified $$\sigma_i^2=\mu(i/n)+a_1(i/n)X_{i-1}^2+b_1(i/n)\sigma_{i-1}^2 \text{ and }\sigma_{0i}^2=\mu_0(i/n)+a_{01}(i/n)X_{i-1}^2+b_{01}(i/n)\sigma_{0,i-1}^2$$ To our advantage, we also have the lower bound for the additional $b_1(\cdot)$ function from assumption B.2. First we point out the proof is not exactly straightforward by adding an additional radii for the $b_1(\cdot)$ coefficient since the third term in the above expression of $\sigma^2_i$ also involves $\sigma_{i-1}^2$ which has evolved differently for $\kappa=(\mu(\cdot), a_1(\cdot), b_1(\cdot))$ and $\kappa=(\mu_0(\cdot), a_{10}(\cdot), b_{10}(\cdot))$. We begin by estimating the difference of the third term
This leads to the recursion
We have
As $M_{\mu}<M_X$,
which will be used for bounding the terms $I$ and $II$ in ((ref)). Towards bounding the first term $I$, we have, along the lines of ((ref)),
For bounding the second term $II$, however, we first look at the summand as follows
Taking sum over $t$, we have
Combining the bounds for $I$ and $II$ we arrive at
In the next part we modify the construction of the sieve to suit the tvGARCH structure.
For the part where we construct exponentially consistent test, note that the only challenge that remains for GARCH processes is to obtain a claim as Claim (ref). Towards that, we use a new definition of sieve in the light of ((ref)) as following:
where $A_n, B_n$ are at least polynomial in $n$ and $\sigma_0^2$ is the standard deviation of $X_0^2$ and $K_n=\max\{K_{1n},K_{2n}\}$. We take $\rho_n \asymp n^{-a}$ with $a<1$, $A_n=B_n(1-n^{1/n}\rho_n), B_n \asymp n^{a_2}$ for sufficiently large $n$ such that $n^{1/n}\rho_n<1$.
The equivalence of Reyni divergence and Hellinger is exactly same as shown in Subsubsection (ref).
We prove Theorem (ref) along very similar lines as outlined for Theorem (ref). The only difference pops up in the KL difference step where existence of second moment is necessary but unfortunately for the tviGARCH case, the variance is infinite. Thus we handle the KL bound in a different way.
Our goal is to obtain a deterministic bound on $\sum_{i=1}^n \frac{|\sigma_i^2-\sigma_{i0}^2|}{\sigma_i^2}$. Denoting $a_1(i/n)$ by $a_{1i}$ we have
where $x_i= (1-a_{1i0})/(1-a_{1i})$. In the first inequality of the above derivation, we have used $\sigma_i^2$ satisfies the following
and $\inf_x(\mu(x))>\rho_{\mu},\sup_xa_1(x)<1-\rho$ and $\inf_x a_1(x)>\rho$ thanks to the closeness of $(\mu,a)$ with $(\mu_0,a_{10})$. For the third term, we have bounds on $x_i$ as follows $x_i\leq 1+ (\|a_1-a_{10}\|_{\infty})/\rho$. Consequently,
by a Taylor series approximation of the binomial term in the third line since $\|a_1-a_{10}\|_{\infty}$ is small. This along with ((ref)) leads to the following bound $$\sum_{i=1}^n\frac{|\sigma_i^2-\sigma_{i0}^2|}{\sigma_i^2}\leq \frac{n}{\rho^2}\|\mu-\mu_{0}\|_{\infty}+\frac{1}{\rho^2}|\sigma_0^2-\sigma_{00}^2|+\frac{n\|a_1-a_{10}\|_{\infty}}{\rho(1-\rho)} \lesssim n\|\mu-\mu_0\|+|\sigma_0^2-\sigma_{00}^2|+n\| a_1-a_{10}\|_{\infty}.$$
Note that for bounding $I$ in the KL term from ((ref)), $\sigma_{i}^{*2}$ satisfy the exact same bounds as ((ref)). For term $II$ as well the bounds follow along the lines of ((ref)).
We consider following sieve for iGARCH:
Within the sieve, $$\sum_{i=1}^n\frac{|\sigma_i^2-\sigma_{i1}^2|}{\sigma_i^2}\leq \frac{n}{\rho_n^2}\|\mu-\mu_{1}\|_{\infty}+\frac{1}{\rho_n^2}|\sigma_0^2-\sigma_{01}^2|+\frac{n\|a_1-a_{11}\|_{\infty}}{\rho_n(1-\rho_n)}.$$ Along the line of Claim (ref) and Claim (ref), the rest of the proof goes through if the radii satisfy $\|\mu-\mu_{1}\|_{\infty}=r_1\leq \frac{\rho_n^2}{8n}$, $|\sigma_0^2-\sigma_{01}^2|=r_2 \leq \rho_n^2/8$ and $\|a_1-a_{11}\|_{\infty}=r_3 \leq \frac{\rho_n(1-\rho_n)}{8n}$.