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.
76,404 characters · 23 sections · 67 citation commands
Simultaneous inference for time-varying models
Time-varying dynamical systems have been studied extensively in the literature of statistics, economics and related fields. For stochastic processes observed over a long time horizon, stationarity is often an over-simplified assumption that ignores systematic deviations of parameters from constancy. For example, in the context of financial datasets, empirical evidence shows that external factors such as war, terrorist attacks, economic crisis, some political event etc. introduce such parameter inconstancy. As bai97 points out, `failure to take into account parameter changes, given their presence, may lead to incorrect policy implications and predictions'. Thus functional estimation of unknown parameter curves using time-varying models has become an important research topic recently. In this paper, we propose a general setting for simultaneous inference of local linear M-estimators in semi-parametric time-varying models. Our formulation is general enough to allow unifying time-varying models from the usual linear regression, generalized regression and several auto-regression type models together. Before discussing our new contributions in this paper, we provide a brief overview of some previous works in these areas.
In the regression context, time-varying models are discussed over the past two decades to describe non-constant relationships between the response and the predictors; see, for instance, fan99, fan2000, hoover98, huang04, lin01, ramsay05, zhang02 among others. 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\in\mathbb{R}^d$ 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. Hypothesis testing is widely used to choose between model I and model II, 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. Little has been known for time-varying models in this direction previously.
The results from time-varying linear regression can be naturally extended to time-varying AR, MA or ARMA processes. However, such an extension is not obvious for conditional heteroscedastic (CH) models. These are difficult to estimate but also often more useful in analyzing and predicting financial datasets. Since engle82 introduced the classical ARCH model and bollerslev extended it to a more general GARCH model, these have remained primary tools for analyzing and forecasting certain trends for stock market datasets. As the market is vulnerable to frequent changes, non-uniformity across time is a natural phenomenon. The necessity of extending these classical models to a set-up where the parameters can change across time has been pointed out in several references; for example stuaricua2005, engle05 and fry08. Towards time-varying parameter models in the CH setting, numerous works discussed the CUSUM-type procedure, for instance, kim00 for testing for changes in the parameters of a GARCH(1,1) time series. 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, linbook99, kokoszka00 or andreou06.
Historically in the analysis of financial datasets, the common practice to account for the time-varying nature of the parameter curves was to transfer a stationary tool/method in some ad hoc way. For example, in mikosch2004, the authors analyzed S&P500 data from 1953-1990 and suggested that time-varying parameters are more suitable due to such a long time-horizon. They re-estimated the parameters for every block of 100 sample points and to account for the abrupt fluctuation of the coefficients, they generated re-estimates of parameters for samples of size $100,200, \ldots.$ This treatment suffers from different degree of reliability of the estimators at different parts of the time horizon. There are examples outside the analysis of economic datasets, where similar approach of splitting the time-horizon has been adapted to fit CH type models. For example, in giaco12, the authors analyzed Italian mortality rates from 1960-2003 using an AR(1)-ARCH(1) model and observed abrupt behavior of yearwise coefficients. Our framework can simultaneously capture these models and provide significant improvements over such heuristic treatments.
A time-varying framework and a pointwise curve estimation using M-estimators for locally stationary ARCH models was 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. Even though the conditional heteroscedastic model remained widely popular in analyzing many different types of econometric data, the topic of simultaneous inference in this field remains relatively untouched. Consider the simple tvARCH(1) model $$X_i= \sigma_i\zeta_i,\quad \zeta_i \sim N(0,1),\quad \sigma_i^2= \alpha_0(i/n)+\alpha_1(i/n)X_{i-1}^2.$$ Typically, it is considered that for large number of realizations, the corresponding parameters $\alpha_0, \alpha_1$ vary smoothly over time and can be modeled as smooth functions $\alpha_0,\alpha_1:[0,1] \to \mathbb{R}$. Pointwise confidence bands for these functions do not help to infer about their overall pattern (like testing for constancy or some specific parametric form). While one remedy could be to subjectively assume a certain class of functions for $\alpha_0,\alpha_1$ such as linear or polynomial and perform a hypothesis test, this can be problematic for many real life datasets. See for example the intercept function for the USGBP analysis in Section (ref). We rather take an objective approach where we do not assume any parametric form as such and wish to establish valid simultaneous inference. In this paper, we therefore derive simultaneous confidence bands which cover $\alpha_0, \alpha_1$ over the whole time interval $t \in (0,1)$ with a given confidence. After construction, one can perform many hypothesis tests such as time-constancy, linearity etc. in one go. To the best of our knowledge, no theoretical results for simultaneous confidence intervals for nonstationary time series were derived before this work.
We next summarize our contributions in this paper. We use Bahadur representations, a Gaussian approximation theorem from zhouwu09 and extreme value theory for Gaussian processes to obtain simultaneous confidence bands for contrasts of parameter curves in very general time-varying models. These intervals provide a generalization from testing parameter constancy to testing any particular parametric form such as linear, quadratic, exponential etc. To deal with bias expansions, we use a theory for locally stationary processes which was recently formalized in dahlhaus2017.
Moving on to some practical applicability of our results, we show how our result applies to time-varying ARCH and GARCH models. For tv(G)ARCH models, we improve the existing conditions in subbarao2008 (we only need that the innovation process has $4+a$ moments for some $a > 0$ compared to 8 moments needed therein) for constructing confidence intervals and provide simultaneous instead of pointwise confidence intervals. We provide an empirical justification of how the coverage can be significantly improved by a wild bootstrap technique and use Gaussian approximation theory to theoretically establish it. Finally we also provide some data analysis and volatility forecasting. First we show for numerous real-life datasets that the time-varying fit does better than the time-constant ones in short-range forecasts. This underlines the importance to decide whether a constant or a time-varying model should be used. One interesting find from our analysis is that simultaneous inference can lead to models where a subset is time-varying and these semi-time varying model can sometimes achieve both statistical confidence and better forecasting ability.
The rest of the article is organized as follows. In Section (ref), we state two specific classes of time series models and the related assumptions. For the sake of better focus and readability, we decided to narrow down the scope of the paper to these specific models. However our theoretical results of M-estimation and the SCBs allow to treat much more general models. The more general assumptions are given in the Appendix (cf. Assumption (ref) therein). In Section (ref) we provide our main results, namely a Bahadur representation of the estimators of the parameter functions and a SCB result for the related contrasts. Section (ref) is dedicated to practical issues which arise when using the SCBs, like estimation of the dispersion matrix of the estimator, bandwidth selection and a wild Bootstrap procedure to overcome the slow logarithmic convergence from the theoretical SCB. Some summarized simulation studies and real data applications can be found in Section (ref). The proofs of the main results are deferred to Appendix, while the proof of several more elementary lemmata and a more general assumption set for tvGARCH processes can be found in the Supplementary material.
Suppose that $\zeta_i$, $i\in\mathbb{Z}$ is a sequence of i.i.d. random variables. We consider the following two time series models. In both cases, $\Theta$ denotes a parameter space specified below in Section (ref).
Case 1 does not directly cover the tvGARCH model, we therefore operate with it separately as Case 2 throughout the paper. In either case, our goal is to estimate $\theta(\cdot)$ from the observations $Z_i^c = (Y_i, X_i^c)$, $i = 1,\ldots,n$.
In this paper, we focus on local M-estimation: Let $K(\cdot)\in \mathcal{K}$, where $\mathcal{K}$ is the family of non-negative symmetric kernels with support $[-1,1]$ which are continuously differentiable on $[-1,1]$ such that $\int_{-1}^{1} |K'(u)|^2 d u > 0$. We consider as objective function $\ell(z,\theta)$ the negative conditional Gaussian likelihood. This reads
For some bandwidth $b_n > 0$, define the local linear likelihood function
where $K_{b_n}(\cdot) := K(\cdot/b_n)$. Let $\Theta' := [-R,R]^k$ with some $R > 0$. A local linear estimator of $\theta(t)$, $\theta'(t)$ is given by
For our main results, we need the following assumptions on our time series models.
For the tvAR($k$) model (cf. richterdahlhaus2017, Example 4.1), one may choose $p = k$, $m_1(x) = x_1$, ..., $m_k(x) = x_k$, $l = 0$, $\nu_0(x) = 1$, leading to the rather strong condition $\sum_{i=1}^{k}|\alpha_i| < 1$ in ((ref)). However, as it can be seen in the proof of Proposition (ref) in the appendix, the condition ((ref)) is only needed to guarantee the existence of the process and corresponding moments. By using techniques which are more specific to the model, one can obtain much less strict assumptions such as $\Theta$ being a compact subset of \[ \{\theta = (\alpha_1,...,\alpha_k,\beta_0) \in \mathbb{R}^{k} \times (0,\infty): \alpha(z) = 1 + \sum_{i=1}^{k}\alpha_i z^i \text{ has only zeros outside the unit circle}\}, \] cf. richterdahlhaus2017, Example 4.1. In the tvARCH case, the above Assumption (ref) asks for $\mathbb{E} |\zeta_1|^{6+a} < \infty$ with some $a > 0$.
In the following, we consider Case 2, the tvGARCH model. In this specific model, the moment conditions can be relaxed to $\mathbb{E}|\zeta_1|^{4+a} < \infty$. The tvGARCH model was for instance studied in the stationary case in garch2004. More recently, pointwise asymptotic results were obtained in tvgarch2013. For a matrix $A$, we define $\|A\|_q := (\|A_{ij}\|_q)_{ij}$ as a component-wise application of $\|\cdot\|_q$. For matrices $A,B$, let $A \otimes B$ denote the Kronecker product and
denote the $k$-fold Kronecker product. Let $\rho(A)$ denote the spectral norm of $A$.
In the important GARCH(1,1) case, a straightforward calculation shows that the condition ((ref)) can be translated to
If $\zeta_0 \sim N(0,1)$, it holds that $\|\zeta_0\|_4^2 = \sqrt{3} \approx 1.73$. bollerslev proved that stationary GARCH(1,1) processes have 4th moments under the exact same condition ((ref)). In Section (ref), Remark (ref) therein, we further talk about the applicability of ((ref)).
We conjecture that also for general GARCH($l,m$) models, ((ref)) is equivalent to the condition \[ \text{ for all }\theta \in \tilde \Theta: \quad\quad \rho(\mathbb{E}[M_0(\theta)^{\otimes 2}]) < 1, \] which would then exactly meet the condition from bollerslev. Note that estimation and the true curve $\theta(\cdot)$ lie in $\Theta$ which is has to be a compact subset of $\tilde \Theta$. Therefore, we automatically ask that all parameters of the GARCH process are nonzero. Again, this condition could in principle be relaxed which would add a significant amount of technicalities.
We discuss the theoretical confidence band result in this section. We directly start with a weak Bahadur representation which plays a key role for introducing simultaneity. For $l \ge 0$, define \[ \mu_{K,l} := \int K(x) x^l dx, \quad\quad\sigma_{K,l}^2 := \int K(x)^2 x^l dx. \] We now have to define some quantities $V(t), I(t), \Lambda(t)$ which are needed to provide the theoretical results. They correspond to the so-called (miss-specified) Fisher information matrices which occur naturally as variance of the M-estimators. These quantities need not to be known in practice because they are estimated. They depend on the so-called stationary approximation $\tilde Y_i(t)$ of the considered time-varying process $Y_i$. In case 1 and case 2, this is given as follows: For $t \in [0,1]$,
For $t\in [0,1]$, let $\tilde Z_j(t) := (\tilde Y_j(t),\tilde Y_{j-1}(t),...)$ denote the infinite vector containing the stationary approximations. We now define
In our theoretical models, these quantities can be related to each other. The following lemma (a direct implication of Propositions (ref) and (ref) in the appendix) summarizes these forms.
In the following, we obtain a weak Bahadur representation of $\hat \theta_{b_n}$ which will be used to construct simultaneous confidence bands. The first part of Theorem (ref) shows that $\hat \theta_{b_n}(t) - \theta(t)$ can be approximated by the expression $V(t)^{-1}\nabla_{\theta}L_{n,b_n}^c(t,\theta(t),\theta'(t))$ as expected due to a standard Taylor argument. The second part of Theorem (ref) deals with approximating this term by a weighted sum of $t$-free terms, namely \[ (nb_n)^{-1}\sum_{i=1}^{n}K_{b_n}(i/n-t)h_i, \quad\quad h_i := \nabla_{\theta}\ell(\tilde Z_i(i/n),\theta(i/n)), \] which is necessary to apply some earlier results from zhouwu10. Let $\mathcal{T}_n := [b_n, 1-b_n]$. For some vector or matrix $x$, let $|x|:= |x|_2$ denote its Euclidean or Frobenius norm, respectively.
Based on the weak Bahadur result, we use results from MR2827528 to obtain a Gaussian analogue of \[ \frac{1}{n b_n}\sum_{i=1}^{n}K_{b_n}(t-i/n) C^{\mkern-1.5mu\mathsf{T}} V(t)^{-1} \nabla_{\theta}\ell(\tilde Z_i(i/n), \theta(i/n)) =:\frac{1}{n b_n}\sum_{i=1}^{n}K_{b_n}(t-i/n)\tilde h_i(i/n) \] for some $C \in \mathbb{R}^{s \times k}$. For a positive semidefinite matrix $A$ with eigendecomposition $A = QDQ^{\mkern-1.5mu\mathsf{T}}$, where $Q$ is orthonormal and $D$ is a diagonal matrix, define $A^{1/2} = Q D^{1/2}Q^{\mkern-1.5mu\mathsf{T}}$, where $D^{1/2}$ is the elementwise root of $D$. Then the following asymptotic statement for simultaneous confidence bands for $\theta(\cdot)$ holds.
Note that for practical use of the SCB in ((ref)), one needs to estimate the bias term, choose a proper bandwidth $b_n$ and estimate $\Sigma_C(t)$. Furthermore, the theoretical SCB only has slow logarithmic convergence, thus one requires huge $n$ to achieve the desired coverage probability. To tackle these type of problems, we discuss practical issues in the next Section (ref).
In this section, we discuss some issues which arise by implementing the procedure from Theorem (ref). We focus on estimation of $\hat \theta_{b_n}$ and optimization of the corresponding SCBs.
There are several possible ways to eliminate the bias term in Theorem (ref). A natural way is to estimate $\theta''(t)$ by using a local quadratic estimation routine with some bandwidth $b_n' \ge b_n$. However the estimation of $\theta''(t)$ may be unstable due to the convergence condition $nb_n^5 \to \infty$ which may be hard to realize together with $n b_n^7 \log(n) \to 0$ from Theorem (ref) in practice. Here instead we propose a bias correction via the a jack-knife method inspired from hardle1986. We define
Since the weak Bahadur representation from Theorem (ref) holds both for $\hat \theta_{b_n/\sqrt{2}}$ and $\hat \theta_{b_n}(t)$, we obtain \[ \sup_{t\in \mathcal{T}_n}\big|V(t)\cdot \{\tilde \theta_{b_n}(t) - \theta(t)\} - (nb_n)^{-1}\sum_{i=1}^{n}\tilde K_{b_n}(i/n-t)h_i\big| = O_{\mathbb{P}}(\tau_n^{(2)} + \beta_n b_n^2 + b_n^3 + (nb_n)^{-1}), \] where $\tilde K(x) := 2\sqrt{2}K(\sqrt{2}x) - K(x)$. Note that the bias term of order $b_n^2$ is eliminated by construction. This shows that Theorem (ref) still holds true for $\tilde \theta_{b_n}(\cdot)$ with kernel $K$ replaced by the fourth-order kernel $\tilde K$ and with no bias term of order $b_n^2$.
In this subsection, we discuss the estimation of $\Sigma_C^2(t)$ since this term is generally unknown but arises in the SCB in Theorem (ref). By Lemma (ref), one has in both cases that $\Lambda(t) = I(t)$, which shows that
As pointed out by a referee, $\Sigma_C^2(t)$ is of the well-known 'sandwich'-form (cf. bollerslevwooldridge1992). Even if the distribution of the innovations $\zeta_i$ is misspecified by the likelihood, one can typically simplify the representation ((ref)). If the distribution of $\zeta_0$ is correctly specified, one has $V(t) = I(t)$ and thus
If the distribution of $\zeta_0$ is misspecified by the likelihood, Lemma (ref) shows that $V(t) = c_0\cdot I(t)$ with some constant matrix $c_0$ which only depends on the fourth moment $\mathbb{E}[\zeta_0^4]$ of $\zeta_0$. Then one has
This also means that the representation ((ref)) is stable under misspecification of the innovation distribution as long as one corrects the expression with the factor $c_0$. To do so, one needs a possibility to estimate $\mathbb{E}[\zeta_0^4]$ from the data. A possibility how to do this for linear processes was discussed in bootstrap4cumulant.
In summary, the representation ((ref)) holds always true, the simpler representations ((ref)) and ((ref)) can be used under additional assumptions or if stable estimators of $c_0$ are available.
To cover all possible situations above, we discuss both estimation of $V(t)$ and $I(t)$. We propose the (boundary-corrected) estimators
where $\hat \mu_{K,0,b_n}(t) := \int_{-t/b_n}^{(1-t)/b_n} K(x) dx$. The convergence of these estimators is given in the next Proposition. Note that the following Proposition also holds if $\widehat \theta_{b_n}'$ in ((ref)) and ((ref)) is replaced by $0$.
This shows uniform consistency of $\hat V_{b_n}(\cdot)$, $\hat I_{b_n}(\cdot)$ if $(\beta_n + b_n)\log(n)^2 \to 0$. Note that in (ii), we need more moments to discuss $\nabla_{\theta}\ell \cdot \nabla_{\theta}\ell^{\mkern-1.5mu\mathsf{T}} \in \mathcal{H}(2M_y,2M_x,\chi,\bar{\bar C})$ ($\bar{\bar C} > 0$). In many special cases, this may be relaxed.
In either case ((ref)) or ((ref)), we define $\hat \Sigma_{C}(t)$ by replacing $V(t),I(t)$ by the corresponding estimators $\hat V_{b_n}(t)$, $\hat I_{b_n}(t)$.
Based on the asymptotic squared error decomposition \[ \big|\hat \theta_{b_n,C}(t) - \theta_{C}(t)\big| \approx \Big|\frac{b_n^2}{2}\mu_{K,2}\theta''_C(t)\Big|^2 + \frac{\sigma_{K,0}^2}{nb_n}\mathrm{tr}(\Sigma_C(t)), \] which can be read off the weak Bahadur representation ((ref)), the squared error global optimal bandwidth choice reads
In practice, $\hat b_n$ is not available due to the unknown quantities on the right hand side, in particular $\theta''(t)$. We therefore adapt a model-based cross validation method from richterdahlhaus2017, which was shown to work even if the underlying parameter curve is only H{\" o}lder continuous and $\nabla_{\theta}\ell(\tilde Z_i(t),\theta(t))$ is uncorrelated. Here, we reformulate this selection procedure for the local linear setting. For $j = 1,\ldots,n$, define the leave-one-out local linear likelihood
and the corresponding leave-one-out estimator \[ (\hat \theta_{b_n,-j}(t), \hat \theta_{b_n,-j}'(t)) = \mathop{\rm argmin}_{\theta \in \Theta, \theta' \in \Theta'} L_{n,b_n,-j}^c(t,\theta,\theta'). \] The bandwidth $\hat b_n^{CV}$ is chosen via minimizing
where $w(\cdot)$ is some weight function to exclude boundary effects. A possible choice is $w(\cdot) := \mathbf{1}_{[\gamma_0,1-\gamma_0]}$ with some fixed $\gamma_0 > 0$. Note that it is important to use the modified local linear approach due to the different bias terms. In richterdahlhaus2017, it was shown that the local constant version of this procedure selects asymptotically optimal bandwidths and works even if a model misspecification is present, i.e. if the function $\ell$ leads to estimators $\hat \theta_{b_n}$ which are not consistent. This motivates that a similar behavior should hold for the local constant version.
The SCB for $\theta_C(t)$ obtained in Theorem (ref) provides a slow logarithmic rate of convergence to the Gumbel distribution. Thus, even for moderately large values of sample size $n$, it is practically infeasible to use such a theoretical SCB as the coverage will possibly be lower than the specified nominal level. First we show an empirical coverage comparison of how far the theoretical confidence intervals lag behind in achieving their nominal coverage. We use the same simulation setting (cf. Section (ref)) for the tvGARCH case:
$$X_i = \sigma_i \zeta_i, \sigma_i^2 = \alpha_0(i/n) + \alpha_1(i/n) X_{i-1}^2 + \beta_1(i/n) \sigma_{i-1}^2,$$ where $\alpha_0(t) = 1.0 + 0.2 \sin(2\pi t)$, $\alpha_1(t) = 0.45 + 0.1 \sin(\pi t)$ and $\beta_1(t) = 0.1 + 0.1\sin(\pi t)$, $\zeta_i$ is i.i.d. standard normal distributed. For estimation, we choose $K(x) = \frac{3}{4}(1 - x^2)\mathbf{1}_{[-1,1]}(x)$ to be the Epanechnikov kernel, $n = 2000, 5000 $ for several different $b_n$. From Table (ref) one can see that the simultaneous coverage is never even positive for the SCB specified in Theorem (ref). The individual coverages are very low for small bandwidth and with higher bandwidth they over-compensate. The performance for $n= 5000$ observations is slightly better, hinting at the logarithmic rate of convergence in Theorem (ref).
We circumvent this convergence issue in this subsection by proposing a wild bootstrap algorithm. Recall the jackknife-based bias corrected estimator $\tilde \theta_{b_n}$ from ((ref)). Let $\tilde{\theta}_C(t)=C^T \tilde{\theta}_{b_n}(t)$. We have the following proposition as the key idea behind the bootstrap method.
The proof of Proposition (ref) is immediate from the approximation rates ((ref)), ((ref)), ((ref)) and ((ref)) in the appendix which, ignoring the $\log(n)$ terms, are of the form $c_n \cdot (nb_n)^{-1/2}\log(n)^{-1/2}$ with \[ c_n \in \{\big(b_n n^{(2\gamma+\varsigma \gamma-\varsigma)/(\varsigma+4\gamma+2\gamma \varsigma)}\big)^{-1/2}, b_n^{1/2}, b_n, (n b_n^{7})^{1/2}, (n b_n^2)^{-1/2}\}, \] where $\gamma > 1$ is arbitrarily large and $\varsigma > 0$.
One can interpret ((ref)) in the sense that $\Sigma_C(t)Q_{b_n}^{(0)}(t)$ approximates the stochastic variation in $\hat \theta_{b_n,C}(t)-\theta_C(t)$ uniformly over $t \in \mathcal{T}_n$. Thus it can be used as a margin for the noise to construct confidence bands, provided one can consistently estimate $\Sigma_C(t)$.
The results shown above only hold for $t \in \mathcal{T}_n$. For inference of some time series models like ARCH or GARCH, large bandwidths are needed to get sufficiently smooth and stable estimators even for a large number of observations. It seems hard to generalize the SCB result Theorem (ref) to the whole interval $t \in (0,1)$. However it is possible to generalize the bootstrap procedure which may be more important in practice:
Note that the additional term in ((ref)) reduces to $Q_{b_n}^{(0)}(t)$ for $t \in \mathcal{T}_n$.
To eliminate the bias inside $t \in \mathcal{T}_n$, it is still recommended to use the jack-knife estimator $\tilde \theta_C(t)$. From Proposition (ref) we obtain
where
The additional factor $N_{b_n}^{(0)}(t) N_{b_n/\sqrt{2}}^{(0)}(t)$ in ((ref)) serves as an indicator how near $t$ is to the boundary. For $t \in \mathcal{T}_n$, this factor is 1 while for $t \in (0,1)\backslash \mathcal{T}_n$, $N_{b_n}^{(0)}(t) N_{b_n/\sqrt{2}}^{(0)}(t)$ may be very small, inducing large diameters of the band near the boundary. Note that the bias correction of the jack-knife estimator $\tilde \theta_C(t)$ may be useless in $t \in (0,1) \backslash \mathcal{T}_n$ since $N_{b_n/\sqrt{2}}^{(1)}(t)N_{b_n}^{(0)}(t)\not= N_{b_n}^{(1)}(t)N_{b_n/\sqrt{2}}^{(0)}(t)$. However it is necessary from a theoretical point of view to use the same estimator for the whole region $(0,1)$ to get a uniform band based on the approximation ((ref)).
In practice, the result ((ref)) can be used as follows: We can create a large number of i.i.d. copies $W_{b_n}^{(boot,debias)}(t)$ of $W_{b_n}^{(debias)}(t)$ by creating i.i.d. copies
where $V_1^*, V_2^*,\ldots ,$ are i.i.d. $N(0, I_{s \times s})$-distributed random variables, and computing $W_{b_n}^{(boot,debias)}(t)$ according to ((ref)). Quantiles of $W_{b_n}^{(debias)}(t)$ then can be determined by using the corresponding empirical quantile of the copies $W_{b_n}^{(boot,debias)}(t)$. Then one can use ((ref)) to construct the confidence band for $\theta_C(t)$. For convenience of the readers, we provide a summarized algorithm of the above discussion.
Algorithm for constructing SCBs of $\theta_C(t)$:
This section consists of some summarized simulations and some real data applications related to our theoretical results. Because of the generality of our theoretical framework, it is impossible to report simulation performance even for the most prominent examples in these different classes. Therefore we restrict ourselves to conditional heteroscedasticity (CH) models for simulations and real data applications. For the time-varying simultaneous band, to the best of our knowledge, there is no or little simulation results reported. For the tvAR, tvMA and tvARMA processes we obtained quite satisfactory results, but they are omitted here to keep this discussion concise.
In this section, we study the finite sample coverage probabilities of our SCBs for theoretical coverage $\alpha = 0.9$ and $\alpha = 0.95$ in the following tvARCH(1) and tvGARCH(1,1) models:
where $\zeta_i$ is i.i.d. standard normal distributed. For estimation, we choose $K(x) = \frac{3}{4}(1 - x^2)\mathbf{1}_{[-1,1]}(x)$ to be the Epanechnikov kernel, $n = 500, 1000, 2000, 5000 $ for several different $b_n$ (the optimal bandwidths ((ref)) are also reported for model (a) and model (b)). For each situation, $N = 2000$ replications are performed and it is checked if the obtained SCB based on ((ref)) contains the true curves in $t \in (0,1)$. In both models we have $\Lambda(t) = I(t) = V(t)$ and therefore estimate $\Sigma_C^2(t) = C^{\mkern-1.5mu\mathsf{T}} I(t)^{-1} C$ via replacing $I(t)$ by $\hat I_{b_n}(t)$ from ((ref)). We obtained the results given in Tables (ref) and (ref). The estimation, for smaller sample sizes $n$, sometimes may lead to difficulties since the optimization routine (optim in programming language R) may not converge. We decided to discard these pathological cases for simplicity. It can be seen that the empirical coverage probabilities are reasonably close to the nominal level for bandwidths close to the optimal ones and they do not differ too much for other bandwidths as well.
In this section, we consider a few real-data applications of our procedure. As mentioned in Section (ref), there are abundant results in the literature about time-varying regression but the results for time-varying autoregressive conditional heteroscedastic models are scarce. Thus it is important to evaluate the performance of our constructed SCBs for these type of models in both theoretical and real data scenarios. Among the popular heteroscedastic models, usually GARCH type models are most difficult to estimate due to the recursion of the variance term.
We consider two examples from the class of conditional heteroscedastic models with two types of financial datasets: one foreign exchange and one stock market daily pricing dataset. As fry08 found out, ARCH models have good forecasting ability for currency exchange type data whereas for data coming from the stock market, GARCH models are preferred. Typically, these daily closing price datasets show unit root behavior and thus instead of using the daily price data, we model the log-return data. 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. Because of the apparent time-varying nature of volatility these log-return data typically show, conditional heteroscedastic models are used for analysis and forecasting.
For the first application, we consider a tvARCH($p$) model with $p=1,2$. It has the following form \[ Y_{i}^2 = \sigma_i^2 \zeta_i^2, \quad\quad \sigma_i^2 = \alpha_0(i/n) + \alpha_1(i/n) Y_{i-1}^2 + \ldots + \alpha_p(i/n)Y_{i-p}^2. \]
Many different exchange rates from 1990-1999 for USD with other currencies were analyzed in fry08 using tvARCH($p$) models with $p=0,1,2$. We collect the data for USD-GBP exchange rates from \url{www.federalreserve.gov/releases/h10/Hist/default1999.htm}. The authors suggested choosing $p=1$ for USD-GBP exchange rates and we also decided to restrict ourselves to fitting a tvARCH(1) model only. Note that in principle, our simultaneous bands can be used to decide whether the additional parameter in a tvARCH(2) model is needed or not. The dataset has sample size 2514 and we use a cross-validated bandwidth of $b_n=0.26$. We also provide the plots for the log-returns and an ACF plot of the squared time series that shows the evidence of conditional heteroscedasticity.
Based on Figure (ref), time-constancy for the parameter curve $\alpha_0(\cdot)$ is rejected at 5% level of significance. For $\alpha_1(\cdot)$, the estimate generally stays below the stationary fit. This can be explained by the geometry of the parameter space of the GARCH model, cf. hillebrand_garch. Also, one can see from the plot of actual log-returns that there are large shocks from 1990 to 1993 compared to those seen in 1993-1999. This can be explained through the high (low) values shown for the estimated curve $\alpha_0(\cdot)$ for the time-period 1990-1993 (1993-1999).
In the empirical analysis of log-returns for stock market data, palm96 and others have found that lower order GARCH models account sufficiently for conditional heteroscedasticity. Moreover, GARCH(1,1) and in a very few cases GARCH(1,2) and GARCH(2,1) models are used and higher order GARCH models are typically not necessary. Another advantage of using GARCH(1,1) over ARCH($p$) models is that one does not need to worry about choosing a proper lag $p$ as GARCH(1,1) can be thought as an ARCH model with $p=\infty$. In this subsection, we implement a time-varying version of GARCH(1,1) and obtain the bootstrapped SCB. A tvGARCH(1,1) model has the following form: \[ Y_{i}^2 = \sigma_i^2 \zeta_i^2, \quad\quad \sigma_i^2 = \alpha_0(i/n) + \alpha_1(i/n) Y_{i-1}^2 + \beta_1(i/n) \sigma_{i-1}^2. \]
As our second example, we choose to analyze the log returns of NASDAQ from January 2011 to December 2018. This is an important index in the US stock market. We collect this one and all other stock index datasets in later analysis from \url{www.investing.com}. Our cross-validated bandwidth is $b_n=0.405$ for this dataset of size $n=1751$. Since our simulations show excellent performance for sample sizes around $n=2000$ and the estimated parameter functions satisfy the parameter restriction $\sup_{0 \leq t \leq 1}(\hat{\beta}_1(t)^2+2 \hat{\alpha}_1(t)\hat{\beta}_1(t)+3\hat{\alpha}_1(t)^2)<1$, it is reasonable to say our simultaneous confidence bands would also be valid here. As one can see from Figure (ref), the time series shows significant lags in its ACF plot after squaring; indicating conditional heteroscedasticity.
One can see that the estimates for $\alpha_0(\cdot)$ is mostly above the corresponding time-constant fit. As mentioned in the caption $\alpha_1(\cdot)$ is time-varying since the SCB does not contain a horizontal line. The fit fluctuates around the time-constant fit. For $\beta_1(t)$, the time-varying fit is below the corresponding time-constant fit. Overall, since $\alpha_1(t)$ is deemed time-varying through this analysis the time-constant hypothesis can be rejected at 5% level of significance.
It is a legitimate question whether time-varying models in forecasting econometric time-series are more useful compared to their time-constant analogue. Note that the main goal of this paper is not to build better forecasting models. The extension to predictive intervals from confidence intervals for conditional heteroscedastic models is not very straight-forward. Moreover, it is unclear how in-fill asymptotics discussed in this paper would extend to forecasting future trends or estimate time-varying functions with time arguments $t > 1$ in the future. Any asymptotic theory would need to consider the rescaling mechanism rigorously, keeping in mind the data observed up to a certain point. In this subsection we show empirically that time-varying models can indeed lead to better forecasts compared to time-constant analogues.
Following starica2003garch and subbarao2008, we show for a wide range of econometric datasets that time-varying models can provide better short range forecasts. We also allow multiple windows of forecasting and multiple start points to highlight why most of these datasets call for a time-varying fit. In the following Tables (ref) and (ref), we use ARCH(1) models for the forex datsets and GARCH(1,1) for the stock market indices. Our POOS (pseudo-out-of-sample) evaluation of forecasting is chalked out as follows:
Define, for a $h-$step ahead forecasting scheme, $\bar{\sigma}^2_{t,t+h}= \frac{1}{h}\sum_{i=t+1}^{t+h}\sum\sigma_{i|t}^2$ where $\sigma_{i|t}^2$ are the $(i-t)$ step ahead forecasts of time-constant or time-varying fit at time $t$. We compare this with the `realized' volatility $\bar{X}^2_{t,t+h}=\frac{1}{h}\sum_{i=t+1}^{t+h}X_i^2$. We then compute the aggregated measure for a start point $s$ as following: $$ AMSE = \frac{1}{n-h-s}\sum_{s+1}^{n-h} (\bar{\sigma}^2_{t,t+h} - \bar{X}^2_{t,t+h}).$$
For the forecasting horizon $h$ values, we choose $h \in \{25,50,75,100,150,200\}$ and start points $s \in \{500,1000\}$. Our forecasting method is the same as the one implemented in the fGARCH R-package. For the time-constant fit at time $t$ we use the data from 1 to $t$ to predict the $h$-step ahead forecast. For the time-varying fit however, it is unclear what the time-varying projection will be. Following subbarao2008, we assume the last $m$ points to be stationary for a small $m$ and use that to obtain the future forecasts. Since $m$ is a tuning parameter, we choose the value that produces the minimum $AMSE$ over $m\in \{100,200, \ldots, 500\}$.
One can see from Table (ref) and (ref) how for a wide range of datasets, starting points and forecasting horizons the AMSE of the time-varying forecast is considerably smaller than the one for the time-constant version. We believe, even if prediction and forecasting is more important from an economist's perspective, this POOS analysis provides a strong motivation to choose a time-varying model over a time-constant one.
We use the following one-step ahead $AMSE_1$, inspired from tvgarch2013 to validate our models:
$$AMSE_1= \frac{1}{n}\sum_{i=1}^n (X_i^2 - \hat{\sigma}^2(i/n))^2$$.
Here, $X_t$ are the log-returns and $\hat{\sigma}^2(\cdot)$ refers to the fitted model using ARCH(1) for foreign exchange datasets, and GARCH(1,1) for stock market indices. In each row we exhibit the best model in bold.
Note that the examples exhibited here show that often the time-constant model has poor one-step ahead forecasting quality compared to their time-varying analogue. However, from Table (ref), one can see that in some of these time-varying models, we have a subset of parameters not rejecting the hypothesis of time-constancy. We suspect that setting a subset of parameters to be time-constant and allowing the rest to vary over time can improve forecasting over both models. One finding of this analysis is that for some of the datasets such as USGBP, the semi-time-varying model may outperform the time-constant model in terms of forecasting. Note that it is easy to tailor and find time-constant fits that allow for even better forecasts, but those models do not have proper confidence (in terms of closeness to the true model) for the already observed data. For the numbers in the above table on the semi-time-varying column, we kept the time-varying coefficients as they are and searched for the best (in terms of $AMSE_1$) constant for the time-constant coefficients among the horizontal lines that fit within the bands entirely. Here we would like also to put a word of caution: Note that the semi-time-varying analysis is somewhat adhoc. In principle, one can also re-run the optimization by fitting only a proper subset as time-varying and the rest as time-constant. We have checked this with multiple of the above datasets and the AMSE were not too different from that reported above. The major takeaway from this analysis remains that our time-varying fit, albeit not meant for prediction and constructed only for building simultaneous confidence intervals, can achieve better forecasts than the corresponding time-constant fits. Additionally, our theory can also lead to new models which have only a subset of coefficients time-varying.
We are grateful to the editor, associate editor and two anonymous referees for their valuable comments and feedback in different rounds which has helped in significantly improving this paper. This research was partially supported by NSF/DMS 1405410.