Evaluation of extremal properties of GARCH(p,q) processes

TL;DR

This paper develops new methods to analyze the extremal properties of GARCH(p,q) processes, enabling better understanding of tail behavior and clustering in financial return models, including complex cases previously uncharacterized.

Contribution

The paper introduces reliable algorithms for evaluating tail indices, extremal index, and cluster functionals for all GARCH(p,q) processes, extending analysis to broader classes of stochastic recurrence equations.

Findings

Derived the forward tail chain for all GARCH(p,q) processes.

Provided methods to evaluate extremal index and tail behavior.

Extended extremal property analysis to integrated and asymmetric GARCH processes.

Abstract

Generalized autoregressive conditionally heteroskedastic (GARCH) processes are widely used for modelling features commonly found in observed financial returns. The extremal properties of these processes are of considerable interest for market risk management. For the simplest GARCH(p,q) process, with max(p,q) = 1, all extremal features have been fully characterised. Although the marginal features of extreme values of the process have been theoretically characterised when max(p, q) >= 2, much remains to be found about both marginal and dependence structure during extreme excursions. Specifically, a reliable method is required for evaluating the tail index, which regulates the marginal tail behaviour and there is a need for methods and algorithms for determining clustering. In particular, for the latter, the mean number of extreme values in a short-term cluster, i.e., the reciprocal of the extremal index, has only been characterised in special cases which exclude all GARCH(p,q) processes that are used in practice. Although recent research has identified the multivariate regular variation property of stationary GARCH(p,q) processes, currently there are no reliable methods for numerically evaluating key components of these characterisations. We overcome these issues and are able to generate the forward tail chain of the process to derive the extremal index and a range of other cluster functionals for all GARCH(p, q) processes including integrated GARCH processes and processes with unbounded and asymmetric innovations. The new theory and methods we present extend to assessing the strict stationarity and extremal properties for a much broader class of stochastic recurrence equations.

Tables

Table 1. Table 1: Values of key stationarity and extremal properties for models A-E for three innovation distributions: 1 t 3 subscript 𝑡 3 t_{3} ; 2 skew t 3 subscript 𝑡 3 t_{3} with ξ = 1 𝜉 1 \xi=1 ; and 3 Gaussian, with Model A - 3 denoting GARCH model formulation A with innovation distribution 3.

Model	$\gamma$	$\eta$	$\kappa$	${\theta }_{X^2}$	${\theta }_{X^U}$	${\theta }_{X^L}$	$\delta$
A - 1	-0.4186	0.071	1.27	0.64	0.76		0.5
A - 2	-0.4611	0.041	1.23	0.69	0.72	0.89	0.80
A - 3	-0.3358	0.023	2.37	0.59	0.72		0.5
B - 1	-0.0400	0.003	1.26	0.38	0.49		0.5
B - 2	-0.0305	0.016	1.09	0.37	0.41	0.57	0.74
B - 3	-0.0155	0.005	1.92	0.16	0.24		0.5
C - 1	-0.0300	0	1	0.21	0.29		0.5
C - 2	-0.0335	0	1	0.29	0.33	0.45	0.71
C - 3	-0.0082	0	1	0.03	0.05		0.5
D - 1	-0.0208	0.007	1	0.21	0.29		0.5
D - 2	-0.0234	0.007	1	0.27	0.31	0.44	0.72
D - 3	-0.0062	0.002	1	0.03	0.05		0.5
E - 1	-0.7461	0.137	0.65	0.27	0.40		0.5
E - 2	-0.7595	0.129	0.68	0.29	0.40	0.45	0.56
E - 3	-0.2411	0.152	0.25	0.13	0.22		0.5

Equations

$${\chi }_X(\tau )=\underset{x\to \infty }{\lim }\mathrm{Pr}(X_{t+\tau }>x\text{leavevmode}\mid \text{leavevmode}{X}_t>x),$$

$$X_t={\sigma }_t{Z}_t$$

$${\sigma }_t^2={\alpha }_0+\sum _{i=1}^q{\alpha }_i{X}_{t-i}^2+\sum _{j=1}^p{\beta }_j{\sigma }_{t-j}^2,t\in \mathbb{Z},$$

$$\sum _{i=1}^q{\alpha }_i+\sum _{j=1}^p{\beta }_j=1$$

$$\frac{\mathrm{Pr}(X_t>ux)}{\mathrm{Pr}(X_t>u)}\to x^{-2\kappa }\text{mbox}and\frac{\mathrm{Pr}(X_t^2>ux)}{\mathrm{Pr}(X_t^2>u)}\to x^{-\kappa }.$$

$$N_n(B)=\#\{k/n\in B:X_k>u_n\},$$

$$\sum _{i=1}^{\infty }i{\pi }_X(i)=1/{\theta }_X,$$

$$(X_0/u,X_1/u,\dots ,X_t/u)\mid X_0>u,$$

$${\chi }_X(\tau )=\mathrm{Pr}({\hat{X}}_{\tau }>1\mid {\hat{X}}_0>1)\text{mbox}and{\theta }_X=\underset{t\to \infty }{\lim }\mathrm{Pr}({\hat{X}}_1<1,\dots ,{\hat{X}}_t<1\mid {\hat{X}}_0>1).$$

$${\mathbf{Y}}_t=\left(\begin{array}{c}{X}_t^2\\ ⋮\\ X_{t-q+1}^2\\ {\sigma }_t^2\\ ⋮\\ {\sigma }_{t-p+1}^2\end{array}\right){\mathbf{A}}_t=\left(\begin{array}{cccc}{\alpha }^{(q-1)}{Z}_t^2& {\alpha }_q{Z}_t^2& {\beta }^{(p-1)}{Z}_t^2& {\beta }_p{Z}_t^2\\ {\mathbf{I}}_{q-1}& 0_{q-1}& 0_{(q-1)\times (p-1)}& 0_{q-1}\\ {\alpha }^{(q-1)}& {\alpha }_q& {\beta }^{(p-1)}& {\beta }_p\\ 0_{(p-1)\times (q-1)}& 0_{p-1}& {\mathbf{I}}_{p-1}& 0_{p-1}\end{array}\right){\mathbf{B}}_t=\left(\begin{array}{c}{\alpha }_0{Z}_t^2\\ 0_{q-1}\\ {\alpha }_0\\ 0_{q-1}\end{array}\right)$$

$${\mathbf{Y}}_t={\mathbf{A}}_t{\mathbf{Y}}_{t-1}+{\mathbf{B}}_t,t\in \mathbb{Z},$$

$${\mathbf{Y}}_t=\left(\begin{array}{c}{X}_t^2\\ {\sigma }_t^2\end{array}\right){\mathbf{A}}_t=\left(\begin{array}{cc}{\alpha }_1{Z}_t^2& {\beta }_1{Z}_t^2\\ {\alpha }_1& {\beta }_1\end{array}\right){\mathbf{B}}_t=\left(\begin{array}{c}{\alpha }_0{Z}_t^2\\ {\alpha }_0\end{array}\right),$$

$$\gamma =\underset{t\to \infty }{\lim }\frac{1}{t}E\left(\ln \Vert {\mathbf{A}}_t{\mathbf{A}}_{t-1}\cdots {\mathbf{A}}_1\Vert \right),$$

$${\gamma }_t=\frac{1}{t}\ln \Vert {\mathbf{A}}_t{\mathbf{A}}_{t-1}\cdots {\mathbf{A}}_1\Vert \text{and}\gamma =\underset{t\to \infty }{\lim }{\gamma }_t,$$

$$\frac{\mathrm{Pr}(\Vert {\mathbf{Y}}_t\Vert >rx,{\mathbf{Y}}_t/\Vert {\mathbf{Y}}_t\Vert \in \cdot )}{\mathrm{Pr}(\Vert {\mathbf{Y}}_t\Vert >x)}\stackrel{v}{\to }{r}^{-\kappa }\mathrm{Pr}({\hat{\mathbf{\Theta }}}_t\in \cdot ),\text{mbox}asx\to \infty ,$$

$$\mathrm{Pr}(X_t^2>rx\mid X_t^2>x)\to r^{-\kappa },\text{mbox}and\mathrm{Pr}({\sigma }_t^2>rx\mid {\sigma }_t^2>x)\to r^{-\kappa },\text{mbox}asx\to \infty .$$

$$R_t$$

$${\mathbf{\Theta }}_t$$

$${\mathbf{\Theta }}_t^{-}$$

$$H_{{\mathbf{\Theta }}_t}(\mathbf{w})=\mathrm{Pr}({\mathbf{\Theta }}_t^{-}\le \mathbf{w}),$$

$$\mathrm{Pr}(R_t>rx,{\mathbf{\Theta }}_t^{-}\le \mathbf{w}\mid R_t>x)\stackrel{v}{\to }\mathrm{Pr}({\hat{R}}_t>r)\mathrm{Pr}({\hat{\mathbf{\Theta }}}_t^{-}\le \mathbf{w})=r^{-\kappa }{H}_{{\hat{\mathbf{\Theta }}}_t}(\mathbf{w}).$$

$$\mathrm{Pr}({\hat{R}}_t>r)=r^{-\kappa }\text{mbox}forr\ge 1.$$

$$H_{{\hat{\mathbf{\Theta }}}_t}(w)=\frac{2}{E(|Z{|}^{2\kappa })}{\int }_0^{{\left(\frac{w}{1-w}\right)}^{1/2}}(1+s^2{)}^{\kappa }{F}_Z(ds),\text{mbox}for0\le w\le 1,$$

$$E(\Vert \mathbf{A}{\hat{\mathbf{\Theta }}}_t{\Vert }^{\kappa })=1$$

$$\mathrm{Pr}({\hat{\mathbf{\Theta }}}_t\in \cdot )=E(\Vert \mathbf{A}{\hat{\mathbf{\Theta }}}_t{\Vert }^{\kappa };\mathbf{A}{\hat{\mathbf{\Theta }}}_t/\Vert \mathbf{A}{\hat{\mathbf{\Theta }}}_t\Vert \in \cdot )$$

$$H_{{\hat{\mathbf{\Theta }}}_t}(\mathbf{w})=E(\Vert \mathbf{A}{\hat{\mathbf{\Theta }}}_t{\Vert }^{\kappa };(\mathbf{A}{\hat{\mathbf{\Theta }}}_t/\Vert \mathbf{A}{\hat{\mathbf{\Theta }}}_t\Vert {)}^{-}\le \mathbf{w}).$$

$$\underset{t\to \infty }{\lim }\frac{1}{t}\ln E\left(\Vert {\mathbf{A}}_t{\mathbf{A}}_{t-1}\cdots {\mathbf{A}}_1{\Vert }^{\kappa }\right)=0$$

$${\mathbf{A}}_t{\mathbf{A}}_{t-1}\cdots {\mathbf{A}}_1={\mathbf{A}}_t\prod _{i=1}^{t-1}({\alpha }_1{Z}_i^2+{\beta }_1),$$

$$E\left[{\left({\alpha }_1{Z}_t^2+{\beta }_1\right)}^{\kappa }\right]=1.$$

$${\int }_{{\mathbb{S}}^{p+q}}E\left[\Vert \mathbf{Aw}{\Vert }^{k}g\left(\frac{\mathbf{Aw}}{\Vert \mathbf{Aw}\Vert }\right)\right]H_k(d\mathbf{w})={\rho }_k{\int }_{{\mathbb{S}}^{p+q}}g\left(\mathbf{w}\right)H_k(d\mathbf{w}),$$

$${\int }_{{\mathbb{S}}^{p+q}}E\left[\Vert \mathbf{Aw}{\Vert }^k\right]H_k(d\mathbf{w})={\rho }_k{\int }_{{\mathbb{S}}^{p+q}}{H}_k(d\mathbf{w})={\rho }_k.$$

Full text

March 8, 2024

Keywords: Cluster of extremes; extremal index; fixed point distributions; GARCH process; multivariate regular variation, particle filtering, stochastic recurrence equations

1 Introduction

Risk management in the stock markets, commonly called market risk management, requires the use of statistical tools and models which aim at reducing the potential size of losses, occurring by sudden drops or growth in the value of stock. Losses can be amplified during periods of large volatility. Risk managers routinely use strategies to handle, model and predict the volatility of daily log-returns, defined as $X_{t}=\log{P_{t}}-\log{P_{t-1}}$, where $P_{t}$, $t=1,2,\ldots,$ is the price of a generic asset. The behaviour of the extreme values of daily log-returns is critically important for market risk management. Isolated extreme values of daily log-return can often be managed, but there is major risk when there is a clustering of these extreme values, and so the study of this dependence structure during extreme events is essential.

It is standard to assume that the series $\{X_{t}\}$ is a stationary series. The most widely adopted models for $\{X_{t}\}$ are the generalised autoregressive conditionally heteroskedastic (GARCH) models (Bollerslev, 1986) and stochastic volatility (SV) models (Taylor, 1986). These models are capable of capturing many of the empirical features of daily log-returns. Both processes have heavy tailed marginal distributions with the leading decay rates the same for both models. Where they differ is in terms of their extremal dependence structure. One of the most common ways to measure this is through the lag $\tau$ tail dependence

[TABLE]

proposed by Ledford and Tawn (2003), with Davis and Mikosch (2009b) terming $\{\chi_{\tau}\}_{\tau\geq 0}$ the extremogram. For SV models Breidt and Davis (1998) show that there is no clustering of extreme values, so that $\chi_{\tau}=0$ for all $\tau>0$. Thus extreme values from SV processes occur in temporal isolation. In contrast, for any GARCH($p,q$) process $\chi_{\tau}>0$ for at least one value of $\tau>0$. But, the values of the extremogram, and other extremal dependence features, are only known for a very restricted subclass of GARCH($p,q$) processes. The aim of this paper is to derive these extremal features for all GARCH($p,q$) models used in typical financial applications and to present algorithms for their evaluation.

We consider GARCH($p,q$) models, for $p\in\mathbb{N}$ and $q\in\mathbb{N}_{+}$, of the form

[TABLE]

where, for every fixed $t\in\mathbb{Z}$, the random variables $Z_{t}$ and $\sigma_{t}$ are independent. Furthermore, we assume that $\{Z_{t}\}$ are independent and identically distributed continuous random variables with $E(Z_{t})=0$ and $\mbox{Var}(Z_{t})=1$. The process $\{\sigma_{t}\}$, commonly referred to as the conditional volatility of $\{X_{t}\}$, is given by

[TABLE]

where $\alpha_{0}>0$ and the parameters $\alpha_{i}\geq 0,i=1,\ldots,q-1$, $\alpha_{q}>0$ and $\beta_{j}\geq 0,j=1,\ldots,p-1$, $\beta_{p}>0$ for $p\geq 1$, have to satisfy additional constraints for the process to be strictly stationary; see Section 2. When $p=0$, equations (1.2) and (1.3) correspond to an ARCH($q$) process and when

[TABLE]

they correspond to an integrated GARCH($p,q$) process, denoted IGARCH($p,q$), which is strictly stationary but not second-order stationary. We develop new numerically robust and efficient methods for assessing whether any GARCH($p,q$) process is strictly stationary. Our results cover all of the GARCH($p,q$) processes including these special cases.

Existing theoretical and computational methods for deriving extremal properties are well established for special cases of the GARCH($p,q$) process, namely: for symmetric $Z_{t}$ with $p=0,q=1$, corresponding to the ARCH(1) process (de Haan et al., 1989) and for $p=q=1$, corresponding to a GARCH(1,1) (Laurini and Tawn, 2012); and for asymmetric $Z_{t}$ with $p=q=1$ (Ehlert et al., 2015). Additional results of other tails probabilities are derived in the two-dimensional case by Collamore et al. (2014) but are effective only for ARCH(1) as further complications arise for the GARCH($1,1$).

For general GARCH($p,q$) models, with arbitrary $(p,q)$ many theoretical extremal properties have been derived by Basrak et al. (2002), Davis and Mikosch (2009a) and Basrak and Segers (2009), including the tails of marginal distributions and some results for the extremal clustering properties. At first sight it seems that these results give everything that is needed for numerical evaluation of the required extremal properties. But this is far from the case, as we will show.

Firstly consider the marginal tail behaviour of $X_{t}$ and $X_{t}^{2}$. Basrak et al. (2002) showed that for GARCH($p,q$) processes there is an explicit theoretical expression for $\kappa\geq 0$ such that for fixed $x>1$ as $u\rightarrow\infty$ then

[TABLE]

These papers only give an asymptotic limiting expression for the evaluation of $\kappa$, but they do not illustrate its application. We find that direct computation using their expression gives very poor numerical performance. Janssen (2010) presents an alternative approach to evaluate $\kappa$, however that approach applies only under the assumption that the innovation $Z_{t}$ has bounded support, ruling out many important distributions used by practitioners, e.g., $Z_{t}$ being Gaussian or $t$-distributed. Furthermore, the associated numerical methods are very slow. We propose the first reliable and computational efficient numerical algorithms for the valuation of $\kappa$, which are valid irrespective of whether $Z_{t}$ are unbounded or bounded. We also find a new formulation for $\kappa$ which gives new insights and we show that for all IGARCH($p,q$) processes $\kappa=1$.

Now consider the extremal dependence/clustering features of the process. Basrak and Segers (2009) and Basrak and Segers (2011) propose computational algorithms for their evaluation. However, these methods have major limitations, which they identify, and only apply to some stochastic recurrence equations with bounded innovations, but do not hold for any GARCH($p,q$) processes, see Section 2.2. So currently no extremal clustering features for GARCH($p,q$) processes, when $\max(p,q)\geq 2$, can be evaluated. We propose an entirely new numerical algorithm to evaluate a range of cluster features for any GARCH($p,q$) process, regardless of the values of $p$ and $q$ and without imposing any restrictive assumptions.

There are a range of extremal dependence features of GARCH process that are of interest to practitioners. The most standard features are summarised by the time-normalized point process $N_{n}$ of exceedances of a level $u_{n}$, defined by

[TABLE]

where $u_{n}$ tends to the upper endpoint of the distribution of $X_{t}$ as $n\rightarrow\infty$, such that $n\bar{F}_{X}(u_{n})\rightarrow\psi$, for any finite $\psi>0$, where $F_{X}$ and $\bar{F}_{X}$ are the marginal distribution and survivor functions of $\{X_{t}\}$ respectively. As $n\to\infty$, $N_{n}$ converges to a compound Poisson process $N$, where events occur as in a homogeneous Poisson process with intensity $\psi\theta_{X}$, where $0<\theta_{X}\leq 1$ is termed the extremal index and with multiplicities distribution denoted by $\pi_{X}(k)$ for $k\geq 1$ (Hsing et al., 1988). We use the term clusters to describe the independent extreme events, with the associated multiplicities corresponding to the number/size of extreme values in each cluster. It follows that

[TABLE]

i.e., the extremal index $\theta_{X}$ is the reciprocal of the limiting mean cluster size of extreme values. The smaller the extremal index then the larger the average number of extreme values per cluster. The special case $\theta_{X}=1$ corresponds to there being no clustering of extremes, so extreme values occur in isolation in time. Thus, minimally there is a need to derive $\{\chi_{\tau}\}_{\tau>0}$ and $\theta_{X}$, and to get the latter we need $\pi_{X}(\cdot)$. Other functionals are also of interest to financial institutions for managing the duration of a stress period or predicting the total amount of losses that can be faced in such stress period, such as the aggregate of excesses over a cluster (total loss). These can also be derived using our methods.

All of these cluster functionals can be obtained from the tail chain of the process, which has been widely used for studying extremal clustering (Rootzén (1988), Smith et al. (1997), Segers (2003), Planinić and Soulier (2018)). The tail chain is defined for a heavy tailed process $\{X_{t}\}$ in the following way. When $u\rightarrow\infty$, if for any $t\in\mathbb{N}_{+}$

[TABLE]

converges weakly to $(\hat{X}_{0},\hat{X}_{1},\ldots,\hat{X}_{t})$, with $\hat{X}_{0}$ non-degenerate, then the limit process $\{\hat{X}_{t}\}_{t\geq 0}$ is termed the forward tail chain. For SV models $\hat{X}_{t}=0$ almost surely for all $t\geq 1$, so large values are not followed by large values for SV processes. In contrast, for any GARCH($p,q$) process at least one $\hat{X}_{t}$, for $t\geq 1$, is non-degenerate, and every element of the tail chain is non-degenerate if $\min(\alpha_{1},\ldots,\alpha_{q})>0$ and $\min(\beta_{1},\ldots,\beta_{p})>0$, see Section 6.1 for further details.

Similarly, there is a backward tail chain $\{\hat{X}_{t}\}_{t\leq 0}$ with identical definition for negative $t$. Here $\hat{X}_{0}$ is the same value for both forward and backward chains. The connections between forward and backward chains were determined by Janssen and Segers (2014). However, in many cases only the forward chain is required to derive the functionals of interest. For example, the extremogram $\{\chi_{X}(\tau);\tau\geq 0\}$ and the extremal index of $\{X_{t}\}$ can be expressed, respectively, as

[TABLE]

In this paper we derive the theory for obtaining the forward tail chain, the extremogram, the extremal index and cluster size distribution for any GARCH($p,q$) process with bounded or unbounded support for the innovations. We provide a new fast, yet accurate, Monte Carlo algorithm for the numerical evaluation of these extremal features. We first obtain the forward tail chain for the squared GARCH($p,q$) process and then use this in a filtering argument, similar to de Haan et al. (1989), to derive the features of interest.

The paper is structured in the following way. In Section 2 we give the required background details for the properties of stationary GARCH($p,q$) processes and the theory of multivariate regular variation that is required for our methodology. In Section 3.1 we give new results for testing stationarity and a new formulation for the tail index $\kappa$. In Sections 3.3 and 4 we derive the tail chain for the series squared and original GARCH($p,q$) processes respectively, with Section 3.2 containing the key particle filtering algorithm which enables us to sample from a $(p+q)$-dimensional extreme state of the tail chain. Critical to the development of this algorithm is the theory of fixed point distributions. Section 5 discusses the novel numerical evaluation of all components of the method, including checking for stationarity and evaluating $\kappa$, and it also illustrates the rapid convergence of the particle filter algorithm. Section 6 has a study of a range of extremal dependence features of the GARCH($p,q$) process over a variety of parameter values. In Section 7 we identify that the methods and algorithms that are developed here in the context of GARCH($p,q$) processes immediately extend to a much broader class of stochastic recurrence equations, and so they are likely to have a much wider impact. The proofs of the theorems are given in the Appendix A.

2 Known properties of GARCH($p,q$) processes

2.1 Strict Stationarity

Let us start by defining strict stationarity for GARCH($p,q$) processes. We focus on the squared GARCH process, $X_{t}^{2}$, and rewrite the process as a stochastic recurrence equation (SRE) as this enables the exploitation of a range of established results (Kesten, 1973) for such processes, e.g., the existence of results for the marginal distribution.

Let the $(p+q)$ vector $\mathbf{Y}_{t}$, the $(p+q)\times(p+q)$ matrix $\mathbf{A}_{t}$ and the $(p+q)$ vector $\mathbf{B}_{t}$ be

[TABLE]

where $\alpha^{(s)}=(\alpha_{1},\ldots,\alpha_{s})\in\mathbb{R}^{s}$, $\beta^{(s)}=(\beta_{1},\ldots,\beta_{s})\in\mathbb{R}^{s}$, $\boldsymbol{I}_{s}$ is the identity matrix of size $s$, $\boldsymbol{0}_{(r\times s)}$ is a matrix of zeros with $r$ rows and $s$ columns and $0_{s}$ is a column vector of zeros having length $s$. In each case here if $s<0$ then these terms are to be interpreted as being dimensionless. Then it follows that the squared GARCH($p,q$) processes satisfies the SRE

[TABLE]

where $\{\mathbf{A}_{t}\}$ and $\{\mathbf{B}_{t}\}$ are each sequences of independent and identically distributed stochastic matrices and vectors.

The formulation of the SRE via (2.1) is due to Francq and Zakoïan (2010, Section 2.2.2, page 29). This SRE formulation is less parsimonious than that of Bougerol and Picard (1992), but has the benefit of covering all GARCH($p,q$) processes, unlike that of Bougerol and Picard (1992) which does not include the case $p=q=1$. In contrast here when $p=q=1$ we have that the terms in expression (2.1) simplify to

[TABLE]

where $\alpha^{(s)}$ and $\beta^{(s)}$ are scalar and none of the $\boldsymbol{I}_{s},\boldsymbol{0}_{r\times s}$ or $0_{s}$ are included.

For general SRE of the form (2.2), but without the specific specification of $\mathbf{A}_{t}$ and $\mathbf{B}_{t}$ corresponding to a GARCH($p,q$) process, Francq and Zakoïan (2010, Theorem 2.4 - page 30/32) show that it is necessary and sufficient that there is a negative top Lyapunov exponent of $\mathbf{A}_{t}$ for the existence of a unique, strictly stationary solution. Before defining the Lyapunov exponent of a stochastic matrix first consider the spectral radius of a deterministic square matrix $\mathbf{A}$, denoted $\rho(\mathbf{A})$. Here $\rho(\mathbf{A})$ is the greatest modulus of its eigenvalues, and an important algebraic result is $\lim_{t\to\infty}t^{-1}\log\lVert\mathbf{A}^{t}\rVert=\log\rho(\mathbf{A})$, where $\lVert\cdot\rVert$ is any norm on the space of $\mathbf{A}$. The extension to a sequence of strictly stationary and ergodic random matrices $\{\mathbf{A}_{t},t\in\mathbb{Z}\}$, for which $E\ln^{+}\lVert\mathbf{A}_{t}\rVert<\infty$ (here $\ln^{+}x=\ln x,$ if $x\geq 1$ and 0 otherwise), is such that the top Lyapunov exponent is

[TABLE]

and $\exp(\gamma)$ is the spectral radius of the sequence $\{\mathbf{A}_{t},t\in\mathbb{Z}\}$. Hence if $E\ln^{+}\lVert\mathbf{A}_{t}\rVert<\infty$ for all $t$, it is necessary and sufficient that $\gamma<0$ for a strictly stationary process $\mathbf{Y}_{t}$.

So for strict stationarity of the squared and original GARCH($p,q$) processes $X_{t}^{2}$ and $X_{t}$ we need $E\ln^{+}\lVert\mathbf{A}_{t}\rVert<\infty$ for $\mathbf{A}_{t}$ given by expression (2.1) and $\gamma<0$. The finite moment condition holds for GARCH($p,q$) processes as $\lVert\mathbf{A}_{t}\rVert<CZ^{2}_{t}+D$, for constants $C>0$ and $D>0$ and so $E\ln^{+}\lVert\mathbf{A}_{t}\rVert<E\ln^{+}(CZ^{2}_{t}+D)=E(\ln(CZ^{2}_{t}+D)\mid CZ^{2}_{t}+D>1)\Pr(CZ^{2}_{t}+D>1)<KE(\ln(Z_{t})\mid CZ^{2}_{t}+D>1)<\infty$, for a suitable constant $K$ and where the last inequality holds as $E(Z_{t})$ is finite by the model definition (1.2). Consequently we only need to check if $\gamma<0$.

Unfortunately, expression (2.4) is not an ideal starting point for evaluating $\gamma$. Instead we also have

[TABLE]

see Francq and Zakoïan (2010, Theorem 2.3 - page 30). So, via expression (2.5), it would appear a relatively simple simulation can be performed to obtain a reliable Monte Carlo estimate of $\gamma$. However, as we will show in Section 5.2, this is far from the case and a mix of careful asymptotic approximation analysis and numerical evaluation is required to evaluate $\gamma$.

In some special cases we do not need to evaluate $\gamma$ to find if the process is strictly stationary, e.g., GARCH($p,q$) processes are always strictly stationary when $\sum_{i=1}^{q}\alpha_{i}+\sum_{j=1}^{p}\beta_{j}\leq 1$; this includes all IGARCH($p,q$) processes. It is also known that $\sum_{j=1}^{p}\beta_{j}<1$ is necessary but not sufficient for strict stationarity. Therefore, numerical evaluation of $\gamma$ is required whenever $\sum_{i=1}^{q}\alpha_{i}+\sum_{j=1}^{p}\beta_{j}>1$ when $\sum_{j=1}^{p}\beta_{j}<1$.

2.2 Existing Results

Basrak et al. (2002) show that there exists a unique stationary solution to the SRE (2.2) and that this solution exhibits a multivariate regular variation property, i.e., for any $t$, any norm $\lVert\cdot\rVert$ and all $r>0$,

[TABLE]

where $\stackrel{{\scriptstyle v}}{{\to}}$ denotes vague convergence (Kallemberg, 1983), $\kappa\geq 0$, and $\hat{\mathbf{\Theta}}_{t}$ is a $(p+q)$-dimensional random vector in the unit sphere (with respect to a norm $\lVert\cdot\rVert$) defined by $\mathbb{S}^{p+q}\subset\mathbb{R}^{p+q}$, and their $(p+q)$ elements will be denoted by $\hat{\mathbf{\Theta}}_{t}=(\hat{\vartheta}_{t}^{(1)},\ldots,\hat{\vartheta}_{t}^{(p+q)})$. If condition (2.6) holds then $\mathbf{Y}_{t}$ is said to exhibit multivariate regular variation with index $\kappa$ and the distribution of $\hat{\mathbf{\Theta}}_{t}$ is termed the spectral measure of the vector $\mathbf{Y}_{t}$. See Resnick (1987) for further details on multivariate regular variation. A consequence of the multivariate regular variation property (2.6) for GARCH($p,q$) processes is that all the marginal variables of $\mathbf{Y}_{t}$ have regularly varying tails with index $\kappa$, so in particular for $r\geq 1$ and all $t$

[TABLE]

So both the squared GARCH($p,q$) process and its variance have regularly varying tails of index $\kappa$.

It is insightful to consider a slightly rearranged version of limit (2.6) and to be specific about which norm we will use. We take the $L_{1}$ norm, and define radial, $R_{t}$, and angular (two variants $\mathbf{\Theta}_{t}$ and $\mathbf{\Theta}^{-}_{t}$) random variables by

[TABLE]

with $\mathbb{S}^{p+q}$ the $(p+q)$ dimensional unit simplex. We have two angular variables as the $p+q$ dimensional variable $\mathbf{\Theta}_{t}$ has redundancy in its final dimension as $\lVert\mathbf{\Theta}_{t}\rVert=1$, and so for studying the distribution of angular variables it is simpler to work with the $p+q-1$ dimensional variable $\mathbf{\Theta}^{-}_{t}$, which is related to $\mathbf{\Theta}_{t}$ by $\mathbf{\Theta}_{t}=(\mathbf{\Theta}^{-}_{t},1-\lVert\mathbf{\Theta}^{-}_{t}\rVert)$ and $\mathbf{\Theta}^{-}_{t}$ being $\mathbf{\Theta}_{t}$ without its last component. We use this $\mathbf{W}^{-}$ notation to create a $(p+q-1)$ dimensional vector from any $(p+q)$ dimensional vector $\mathbf{W}$ on the simplex $\mathbb{S}^{p+q}$ throughout. Furthermore, for $\mathbf{w}\in\mathbb{R}^{p+q-1}$, we use the notation

[TABLE]

with vector algebra, here and elsewhere, interpreted as being componentwise.

We will denote the limit random variables, that arise in limit (2.6), for $(R_{t},\mathbf{\Theta}_{t},\mathbf{\Theta}^{-}_{t})$ by $(\hat{R}_{t},\hat{\mathbf{\Theta}}_{t},\hat{\mathbf{\Theta}}^{-}_{t})$, with the distribution function of $\hat{\mathbf{\Theta}}_{t}$, denoted by $H_{\hat{\mathbf{\Theta}}_{t}}$, defined similarly to distribution (2.9). Subsequently $H_{\hat{\mathbf{\Theta}}_{t}}$ is referred to as the spectral measure, the term coming from multivariate regular variation terminology (Resnick, 1987). Then, for $r\geq 1$, as $x\to\infty$, the limit (2.6) becomes,

[TABLE]

From the first expression for the asymptotic form in limit (2.10) we see that the radial variable $R_{t}$ and the angular variables $\mathbf{\Theta}_{t}$ become asymptotically independent, as the radial variable $R_{t}$ grows due to $x\to\infty$, i.e., the variables $\hat{R}_{t}$ and $\hat{\mathbf{\Theta}}_{t}$ are independent. The second term in this limit shows that $\hat{R}_{t}$ is a Pareto random variable with tail index $\kappa$, i.e.,

[TABLE]

There is additional structure imposed on both $H_{\hat{\mathbf{\Theta}}_{t}}$ and $\kappa$ by the GARCH($p,q$) process, which is identified by Basrak and Segers (2009) and Janssen (2010) respectively. We discuss this structure in each case below.

For the GARCH($1,1$) process, Laurini and Tawn (2012) provided an expression for the spectral measure, for a different description of the angular variable to that used here. For the choice of the angular variable (2.2), for all $t$, their result translates to

[TABLE]

where $F_{Z}$ is the distribution function of the innovations $Z_{t}$. When $\max(p,q)\geq 2$, through highly skilled use of the multivariate regular variation structure, Basrak and Segers (2009, Propositions 3.3, 5.1) show, that when $\mathbf{A}$ is independent and identically distributed to $\mathbf{A}_{t}$, that

[TABLE]

and uniquely

[TABLE]

where the notation $E(X;Y):=E(X\mathbf{1}_{Y})$ where $\mathbf{1}_{Y}$ is the indicator of the event $Y$. Thus

[TABLE]

Basrak and Segers (2009, p. 1075) propose an approach to simulate from $H_{\hat{\mathbf{\Theta}}_{t}}(\mathbf{w})$ for an SRE of the form (2.2), with the required distribution $H_{\hat{\mathbf{\Theta}}_{t}}$ being the invariant distribution of a Markov chain and MCMC methods used for its evaluation. However, this method cannot be used for general GARCH($p,q$) processes for the following reasons. Firstly, they make an assumption that $\mathbf{A}_{t}$ is bounded, which excludes the possibility of $Z_{t}$ being, for example, Gaussian or $t_{\nu}$ distributed. Much more critically though, Basrak and Segers (2011) note that the proof that $H_{\hat{\mathbf{\Theta}}_{t}}(\mathbf{w})$ is the stationary distribution of the Markov Chain that they proposed was flawed, and their claimed results only hold under one of two very specific conditions on the matrix $\mathbf{A}_{t}$ in the SRE framework (2.2). Neither of these conditions are satisfied by the form of $\mathbf{A}_{t}$ for GARCH$(p,q)$ processes, when $\max(p,q)\geq 2$, even with a bounded $Z_{t}$. Thus the algorithm proposed by Basrak and Segers (2009) cannot be used for obtaining $H_{\hat{\mathbf{\Theta}}_{t}}(\mathbf{w})$ for a GARCH$(p,q)$ process. Our approach in Section 3.2 overcomes both of these restrictions.

Next we focus on how $\kappa$ is determined. In particular, Basrak et al. (2002) showed that there exist a $\kappa>0$ which is the unique positive solution of the equation

[TABLE]

where here and throughout the matrix norm we use will use is $\lVert\mathbf{A}\rVert=\sum\lvert a_{ij}\rvert$, where $a_{ij}$ is $(i,j)$th element of matrix $\mathbf{A}$. For the GARCH($1,1$) process Mikosch and Stărică (2000) show that $\kappa$ is simple to evaluate using expression (2.15). Specifically, taking $\mathbf{A}_{t}$ as in expression (2.3) we have that

[TABLE]

from which it simply follows that expression (2.15) holds and $\kappa$ satisfies

[TABLE]

Setting $\beta_{1}=0$ for the GARCH($1,1$)  process gives the same result for $\kappa$ derived by de Haan et al. (1989) for the ARCH(1) process. For general GARCH($p,q$) processes no such existing simplification of equation (2.15) gives an easier expression for $\kappa$. So it is natural to try to find $\kappa$ by a numerical solution of the limit equation (2.15). However, direct numerical solution is non-trivial due to numerical instabilities.

The only existing feasible method to evaluate $\kappa$ was proposed by Janssen (2010, Proposition 4.3.1), which exploits Kesten (1973, Proof of Theorem 3). With $\mathbf{A}$ specified in (2.2) the conditions required for the results of Kesten (1973) apply and the equality

[TABLE]

holds for all continuous functions $g$, all unit measures $H_{k}$ on the space $\mathbb{S}^{p+q}$, and where $\rho_{k}$ is a constant. The special case of $g\equiv 1$ in (2.17) gives the simplification

[TABLE]

For any given $k\in(0,\infty)$, whatever the chosen unit measure $H_{k}$, if the pair $(\rho_{k},H_{k})$ satisfies equality (2.17), then $\rho_{k}$ is determined solely by $k$, i.e., not by the choice of $H_{k}$.

As we know from condition (2.13) that the $\kappa$ moment of $\lVert\mathbf{A}\hat{\mathbf{\Theta}}_{t}\rVert$ is equal to $1$, where $\hat{\mathbf{\Theta}}_{t}\sim H_{\hat{\mathbf{\Theta}}_{t}}$ on the space $\mathbb{S}^{p+q}$ then it follows from property (2.18) that $\rho_{k}=1$ when $k=\kappa$. Kesten (1973) and Janssen (2010) shows that there is only one solution to equation (2.13), so $\kappa$ is the unique solution of

[TABLE]

and that the unit measure $H_{\kappa}$ must correspond to the distribution function $H_{\hat{\mathbf{\Theta}}_{t}}$. Thus if we can find, or simulate from, $H_{\hat{\mathbf{\Theta}}_{t}}$ we can find $\kappa$. To evaluate $\kappa$ all that is required is to define a class of unit measures $H_{k}$, over $k\in(0,\infty)$, which contains within it as an interior point $H_{\hat{\mathbf{\Theta}}_{t}}$, and then vary $k$ until property (2.19) is found.

Janssen (2010) proposes an algorithm to simulate from a class of functions $H_{k}$ which adapts the invalid algorithm of Basrak and Segers (2009). This gives a valid method for finding $\kappa$ and $H_{\hat{\mathbf{\Theta}}_{t}}$ but critically it only applies when the innovations $Z_{t}$ have bounded support. Therefore there is a need for an algorithm to cover cases where $Z_{t}$ have unbounded support. In Section 3.2 we describe such an algorithm that applies whatever the support of $Z_{t}$. Furthermore, our methods for calculating $\kappa$, have substantial computational performance and efficiency gains compared to the algorithm of Janssen (2010); see Section 5.6.

3 Extremal properties of squared GARCH processes

3.1 New formulations for $\gamma$ and $\kappa$

Expression (2.5) suggests using Monte Carlo for the evaluation of $\gamma$, by taking $t$ to be very large. This approach suffers from serious numerical instabilities. As the norm of the product is tending to zero it seems sensible to first normalise the size of the individual terms in the product. We do this be scaling each matrix $\mathbf{A}_{i}$ by its largest (in magnitude) eigenvalue, which we denote by $\lambda_{i}$, with $\lambda_{i}>0$ for all $i$. From Kesten and Spitzer (1984) for an arbitrary $\mathbf{A}_{i}$, satisfying conditions of Section 2.1, it is guaranteed that $\lambda_{i}$ is simple and exceeds all other eigenvalues in absolute value. Let

[TABLE]

so that the product (2.5) can be re-written as

[TABLE]

and let

[TABLE]

We study the limit behaviour of both components in equation (3.1) in Theorem 3.1, whose proof is postponed to Appendix A.

Theorem 3.1**.**

If $\mathbf{A}_{t},t\in\mathbb{Z}$, is a sequence of independent and identically distributed random matrices, with non-negative entries, and

[TABLE]

where $\lambda_{i}$ is the magnitude of the largest eigenvalue of $\mathbf{A}_{i}$ and $\eta\in\mathbb{R}$ is such that

[TABLE]

then $\gamma=\lim_{t\to\infty}\gamma_{t}$ if and only if

[TABLE]

To assess this result in terms of what is already known we first compare with the GARCH($1,1$) process. In that case Mikosch and Stărică (2000) show that $\gamma=E[\ln(\alpha_{1}Z_{t}^{2}+\beta_{1})]$, but the only, and hence largest in magnitude, eigenvalue of $\mathbf{A}_{t}$ is $\alpha_{1}Z_{t}^{2}+\beta_{1}$, thus Theorem 3.1 is identical to their result when $\eta=0$. To show that $\eta=0$ for all GARCH($1,1$) processes note that

[TABLE]

so

[TABLE]

Hence $\eta=0$ from Theorem 3.1.

The practical evaluation of the tail index $\kappa$ is not discussed by Basrak and Segers (2009) or subsequent authors. Using representation (2.15) as the basis for numerical evaluation of $\kappa$ for a GARCH($p,q$) process turns out to be trivial only for very ARCH(1) and GARCH($1,1$) processes. Monte Carlo is essential, but solving the limiting equation (2.15) is non-trivial when $\max(p,q)\geq 2$ due to major numerical instabilities, In Theorem 3.2 we present a new representation which provides both insight into which factors determine $\kappa$ as well as a basis for a method of evaluation with greater numerical stability.

Theorem 3.2**.**

Under the same notation and conditions of Theorem 3.1 and that $\gamma<0$, the unique solution $\kappa>0$, of the limiting equation (2.15), satisfies

[TABLE]

Furthermore, it follows that

[TABLE]

For all strictly stationary GARCH(1,1) processes, Mikosch and Stărică (2000) show that $\kappa$ must satisfy $E[(\alpha_{1}Z_{t}^{2}+\beta_{1})^{\kappa}]=1$, but this is simply $E[\lambda^{\kappa}]=1$, so Theorem 3.2 gives that $\eta=0$ for all GARCH($1,1$) processes, as shown directly above.

The results of Theorem 3.2 are particularly powerful as they allow the simple evaluation of $\kappa$, if $\eta$ is known, or vice-versa. Theorem 3.1 gives a limiting expression from which to approximate $\eta$ and hence $\kappa$ can be found approximately, but better still, Section 3.2 gives a reliable numerical method to calculate $\kappa$ and then Theorem 3.2 provides the ideal way to find $\eta$. These approaches are illustrated in Section 5. Finally as an immediate consequence of Theorems 3.1 and 3.2 if we know the process is stationary, and we know $\kappa$, we can directly calculate $\gamma$ using the following result.

Theorem 3.3**.**

Under the same notation and conditions of Theorem 3.1, if $\gamma<0$ and limiting equation (2.13) gives $\kappa>0$, then

[TABLE]

Unlike Theorem 3.1, we cannot use Theorem 3.3 to test for stationarity of the process as this result only provides an expression for $\gamma$ given that $\gamma<0$, i.e., stationarity needs to be confirmed prior to its use. The only comparable existing result to Theorem 3.3 is by Kesten and Spitzer (1984) where it is shown that $\gamma\leq\ln(E(\lambda))$ for general random matrix Markov processes. Finally, Theorem 3.4 shows that for all IGARCH($p,q$) processes $\kappa=1$, in which case expression (3.4) gives that $\eta=-\ln(E(\lambda))$ and Theorem 3.3 gives that $\gamma=E(\ln\lambda)-\ln(E\left(\lambda\right))$.

Theorem 3.4**.**

For all IGARCH($p,q$) with $\gamma<0$, we have that $\kappa=1$. Furthermore, if we have a stationary GARCH($p,q$) process with $\kappa=1$ then it must be an IGARCH($p,q$) process.

This seems to be the first time that it has been claimed that any IGARCH($p,q$) process with $\max(p,q)\geq 2$ has $\kappa=1$, although it was proved for $p=q=1$ by Mikosch and Stărică (2000). The finite mean and infinite variance of the IGARCH($p,q$) process implies that $0.5<\kappa\leq 1$. So our result gives much more, e.g., all IGARCH($p,q$) processes have $E(|X_{t}|^{2-\epsilon})<\infty$ for any $\epsilon\in(0,2]$. This finding for $\kappa$ is not too surprising though as the variance of $X_{t}$ is infinite when $\sum_{i=1}^{q}\alpha_{i}+\sum_{j=1}^{p}\beta_{j}=1$ but was finite when this sum is less than $1$, suggesting $\kappa$ for the IGARCH($p,q$) was a critical boundary point for having a finite variance.

3.2 Evaluating the Spectral Measure and the Tail Index

This section gives the details of our algorithm for sampling from the limit distribution $H_{\hat{\mathbf{\Theta}}_{t}}$ and then uses this algorithm repeatedly to find $\kappa$. The algorithm requires no assumptions on the support for $Z_{t}$. Throughout we take $t=0$, both to help simplify notation here and as it will be from time $t=0$ that we start the tail chains in Section 3.3. We will first assume that $\kappa$ is known and present Algorithm 1 for generating from $H_{\hat{\mathbf{\Theta}}_{0}}$ and then discuss the case when $\kappa$ is unknown.

To simulate from the spectral measure $H_{\hat{\mathbf{\Theta}}_{0}}$, defined via (2.14), our approach is to introduce a stochastic process whose invariant distribution is $H_{\hat{\mathbf{\Theta}}_{0}}$. We will then use sequential importance sampling Doucet et al. (2000) to generate approximate samples from the state of this stochastic process at consecutive time-steps. We denote the stochastic process by $\widetilde{\mathbf{\Theta}}_{s}$, for $s=0,1,\ldots,$ and denote its joint distribution function at iteration $s$ by

[TABLE]

This stochastic process is constructed such that as $s\to\infty$, $\widetilde{\mathbf{\Theta}}_{s}\stackrel{{\scriptstyle d}}{{\to}}\hat{\mathbf{\Theta}}_{0}\sim H_{\hat{\mathbf{\Theta}}_{0}}$. We perform sequential updates until it appears that the distribution of the state of the stochastic process has converged to the invariant distribution of the process, $H_{\hat{\mathbf{\Theta}}_{0}}(\mathbf{w})$. The samples at this final iteration are then taken as samples from $H_{\hat{\mathbf{\Theta}}_{0}}$. This algorithm is similar to the use of sequential Monte Carlo for sampling from Feynman-Kac distribution (e.g. Del Moral and Miclo, 2000) and fixed point distributions (e.g. Del Moral and Miclo, 2003).

Let $\widetilde{\mathbf{\Theta}}_{s}$, for $s=0,1,\ldots,$ be a Markov process, with initial state an arbitrary distribution on $\mathbb{S}^{p+q}$ and whose transitions for $s\geq 1$ are given by

[TABLE]

where, as above, $E(X;Y):=E(X\mathbf{1}_{Y})$. By construction, the invariant distribution of this process is $H_{\hat{\mathbf{\Theta}}_{0}}$. To see this notice that if $\widetilde{\mathbf{\Theta}}_{s-1}$ is drawn from $H_{\hat{\mathbf{\Theta}}_{0}}$, then the right-hand side of (3.6) is equal to

[TABLE]

As $E(\lVert\mathbf{A}\hat{\mathbf{\Theta}}_{0}\rVert^{\kappa})=1$, this is equal to the definition of $H_{\hat{\mathbf{\Theta}}_{0}}$ given by expression (2.14).

Furthermore, if we have a sample from $\Pr(\widetilde{\mathbf{\Theta}}_{s-1}\in\cdot)$, we can use importance sampling to generate a sample from $\Pr(\widetilde{\mathbf{\Theta}}_{s}\in\cdot)$. This would involve first simulating a value for $\widetilde{\mathbf{\Theta}}_{s}$ via

[TABLE]

and assigning this value an importance sampling weight proportional to $\lVert\mathbf{A}\widetilde{\mathbf{\Theta}}_{s-1}\rVert^{\kappa}$. Thus we can use sequential importance sampling to generate samples of $\widetilde{\mathbf{\Theta}}_{s}$ values for $s\geq 1$. Specifically, we implement Algorithm 1, with our choice of initial distribution, in step 1, being chosen to be close to $H_{\hat{\mathbf{\Theta}}_{0}}(\mathbf{w})$, so as to speed up convergence, see Section 5.3 for details. For details of how we determine convergence in step 6 of Algorithm 1 see Section 5.4.

Our approach is closely related to approaches for sampling from quasi-stationary distributions (see Griffin et al., 2017, and references therein). This can be most clearly seen in situations where $\lVert\mathbf{A}\rVert$ is bounded. In this case we can define a stochastic process with transitions given by (3.6) and with killing at each iteration, with the probability of survival being proportional to $\lVert\mathbf{A}\widetilde{\mathbf{\Theta}}_{s}\rVert^{\kappa}$. In this case the spectral measure, $H_{\hat{\mathbf{\Theta}}_{0}}$, is the quasi-stationary distribution of the process, i.e., the stationary distribution of the process conditional on survival.

Now consider the situation when $\kappa$ is not known. For a trial value of $k$ (for $\kappa$), apply Algorithm 1 until convergence, giving a sample of weighted particles $\{\widetilde{\mathbf{\Theta}}^{(j)}(k),m^{(j)}(k)\}_{j=1}^{J}$ after the chain is deemed to have converged. Using these particles approximate the expectation $E(\lVert\mathbf{A}\widetilde{\mathbf{\Theta}}_{0}\rVert^{k})$ using the Monte Carlo sample by

[TABLE]

where $\tilde{\rho}_{k}$ is the Monte Carlo approximation of $\rho_{k}$ and $F_{Z}$ is the distribution function of the GARCH($p,q$) process. We repeat this evaluation over $k>0$ until we find the unique value of $k$ which gives this weighted mean to be equal to 1. This value is $k=\kappa$.

3.3 Generation of the tail chain of the squared process

The tail chain $\{\hat{X}^{2}_{t}\}_{t\geq 0}$ can be evaluated using Algorithm 2. There are two stages to the algorithm, initialisation and propagation of the chain. Key to getting the tail chain is finding the joint behaviour of $\hat{\mathbf{\Theta}}_{t}=(\hat{\vartheta}_{t}^{(1)},\ldots,\hat{\vartheta}_{t}^{(p+q)})$ over time.

For initialisation we first need to consider the behaviour of the process conditional on it being in an extreme state, and we take the time of this, for convenience, to be $t=0$. Focusing on limit (2.10) when $t=0$, we have that the limit variables $(\hat{R}_{0},\hat{\mathbf{\Theta}}_{0})$, are independent with distributions (2.11) and $H_{\hat{\mathbf{\Theta}}_{0}}$, given by expression (2.14), respectively. We initialise the chain when $\hat{X}^{2}_{0}>1$, which is equivalent to $\hat{R}_{0}\hat{\vartheta}_{0}^{(1)}>1$. For propagation of the chain we can use the established results of Basrak et al. (2002), which give that for $t\geq 1$ that $\hat{\mathbf{\Theta}}_{t}=\mathbf{A}_{t}\hat{\mathbf{\Theta}}_{t-1}$ and hence

[TABLE]

We extract the tail chain $\{\hat{X}^{2}_{t}\}_{t\geq 0}$ from the product of $\hat{R}_{0}$ and $\hat{\vartheta}_{t}^{(1)}$ for all $t$.

If interest is in the tail chain of the $\hat{\sigma}_{t}^{2}$ process instead, i.e., $\{\hat{\sigma}^{2}_{t}\}_{t\geq 0}$ given that $\hat{\sigma}^{2}_{0}>1$, then exactly the same approach can be taken as in Algorithm 2 but with the condition in step 2 changed to $\hat{R}_{0}\hat{\vartheta}_{0}^{(q+1)}>1$.

Here $T$ is selected so that $\hat{X}^{2}_{t}=\hat{R}_{0}\hat{\vartheta}_{t}^{(1)}<1$ for all $t>T$ with probability as close to $1$ as possible. This is achievable as all components of $\hat{\mathbf{\Theta}}_{t}$ have negative drift and converge to 0 almost surely. In practice, $T$ is taken as large as possible subject to limits of storage and computational time, we took $T=1000$ to save running chains unnecessarily long, but for processes with weak extremal dependence $T=50$ is more than sufficient.

3.4 The evaluation of cluster functionals

From repeated realisations of the tail chain $\{\hat{X}_{t}^{2}\}_{t\geq 0}$ for the squared GARCH process we can derive the properties of key cluster functionals for the $X_{t}^{2}$ process, e.g., the extremogram, the extremal index and the cluster size distribution.

First note that the extremogram for the squared GARCH$(p,q)$ process is

[TABLE]

Thus we can numerically determine $\chi_{X^{2}}(\tau)$ as the proportion of tail chains starting above 1 at $t=0$ with an exceedance at $t=\tau$ over different replicate tail chains. Any required precision of this value can be achieved by a suitable selection of the number of Monte Carlo replicate tail chains. To derive both the extremal index and the cluster size distribution we first define the measure introduced by Rootzén (1988), namely

[TABLE]

i.e., the probability that there are at least $i$ values in a cluster given that we look at a cluster only forwards in time from an arbitrary exceedance. O’Brien (1987) showed that the extremal index $0<\theta_{X^{2}}\leq 1$, is given by $\theta_{X^{2}}=\theta_{X^{2}}^{(1)}$ and Rootzén (1988) showed that the cluster size distribution is given by

[TABLE]

with the reciprocal of the mean of this distribution being $\theta_{X^{2}}$.

4 Tail chain properties for GARCH processes

First note that if $X_{t}^{2}$ is regularly varying with index $\kappa>0$ and if

[TABLE]

where $0<\delta<1$ then it follows that $X_{t}$ is a regularly varying random variable, with index $2\kappa$, in both its upper and lower tails. Details of the evaluation of $\delta$ are given in Section 5.7.

To translate results about the tail chain $\hat{X}^{2}_{t}$ of the squared GARCH process, into properties for the tail chain $\hat{X}_{t}$ of the GARCH process we adopt a similar strategy to de Haan et al. (1989) and Ehlert et al. (2015). It is key to recognise that there are two tails chains for $X_{t}$, an upper and a lower tail chain $\hat{X}^{L}_{t}$ and $\hat{X}^{U}_{t}$ respectively, with

[TABLE]

where $I_{t}$ is a sequence of independent and identically distributed Bernoulli$(\delta)$ variables, with $I_{t}=\{0,1\}$ with respective probabilities $\{1-\delta,\delta\}$ and where $\delta$ is given by limit (4.1), The sequence $I_{t}$ is also independent of $\{\hat{Y}_{t}\}$.

Many properties of the cluster functions for $\hat{X}^{U}_{t}$ and $\hat{X}^{L}_{t}$ chains can simply be derived using Monte Carlo methods from the $\hat{X}^{2}_{t}$ tail chains by using Bernoulli thinning, implied by expression (4.2), but some functionals can be explicitly determined, we study a few of these below.

First, we focus on the $\hat{X}^{U}_{t}$ tail chain, corresponding to positive events in the GARCH process. The extremogram for $\hat{X}^{U}_{t}$ is given by

[TABLE]

where $\chi_{X^{2}}(\tau)$ is the lag $\tau$ extremogram for the $X^{2}_{t}$ process. An event in the tail chain for $\{\hat{X}^{2}_{t}\}$ with $i$ exceedances of the level $1$ does not occur in the tail chain of $\{\hat{X}^{U}_{t}\}$ with probability $\delta^{i}$. Therefore summing over all event lengths, the probability of no exceedances of level $1$ from an event in the tail chain of $\{\hat{X}^{2}_{t}\}$ is given by

[TABLE]

where $\pi_{X^{2}}(i)$ is the probability that a cluster of length $i$ in the $\{X_{t}^{2}\}$ series, see Section 3.4. The probability of a cluster of length $j$ in the $\hat{X}^{U}_{t}$ tail chain is

[TABLE]

where the denominator corresponds to conditioning on the cluster for $\{\hat{X}^{2}_{t}\}$ process being retained for the $\{\hat{X}^{U}_{t}\}$ series. Then the mean cluster size, $1/\theta_{X^{U}}$ for $\{\hat{X}^{U}_{t}\}$ is given by

[TABLE]

where $\theta_{X^{2}}$ is the extremal index of the squared GARCH process, see Section 3.4. So the extremal index of $\{X^{U}_{t}\}$ is

[TABLE]

Similarly, for the lower tail behaviour of $\{X_{t}\}$ it follows that the extremogram is

[TABLE]

the probability of no values below the level $-1$ by $X_{t}^{L}$ from an event in the tail chain of $\{\hat{X}^{2}_{t}\}$ is given by

[TABLE]

the probability of a cluster of length $j$ in the $\hat{X}^{L}_{t}$ tail chain is

[TABLE]

and extremal index of $\{X^{L}_{t}\}$ is $\theta_{X^{L}}=\theta_{X^{2}}(1-\mathit{\Pi}^{L})/(1-\delta)$.

5 Numerical Methods

5.1 Introduction

Throughout this section a range of GARCH($p,q$) models will be illustrated. The details of these models are given here and will be referenced subsequently as GARCH models A-E, where

A

: $p=q=2$ with $(\alpha_{1},\alpha_{2},\beta_{1},\beta_{2})=(0.3,0.15,0.2,0.1)$

B

: $p=q=2$ with $(\alpha_{1},\alpha_{2},\beta_{1},\beta_{2})=(0.07,0.04,0.8,0.08)$

C

: $p=q=1$ with $(\alpha_{1},\beta_{1})=(0.1,0.9)$

D

: $p=q=2$ with $(\alpha_{1},\alpha_{2},\beta_{1},\beta_{2})=(0.07,0.03,0.8,0.1)$

E

: $p=2$, $q=0$ with $(\alpha_{1},\alpha_{2})=(1.2,0.5)$.

We selected these models to give a variety of stationarity and extremal behaviours. Models A and B are second order stationary. Models C and D, which are IGARCH models, and model E are not second order stationary (Francq and Zakoïan, 2010, p. 35).

Let $\phi=\sum_{i=1}^{q}\alpha_{i}+\sum_{i=1}^{p}\beta_{i}$ be the sum of the meaningful GARCH parameters. The parameter $\phi$ is increasing from model A to E with $\phi=1$ for models C and D. We will show that all these models are strictly stationary and that the marginal tail index $\kappa$ decreases with increasing $\phi$ for these models. From Theorem 3.4 we have that when $\phi=1$ then $\kappa=1$, we will also illustrate this numerically for models C and D. Model B corresponds to the model studied by Mikosch and Stărică (2000). Case C, though being an IGARCH(1,1) process, is not covered by previous results of Laurini and Tawn (2012) given its IGARCH form, but is of interest here as it helps to illustrate the new methods in a case where analytical solutions are possible.

In Sections 5 and 6 we take the the distribution of the innovation process $Z_{t}$, to be standard Gaussian, a scaled Student-$t_{\nu}$ distribution, and the skew Student-$t_{\nu}$ distribution introduced in Azzalini and Capitanio (2003). In each case the innovation distribution has zero mean and unit variance. First consider the univariate skew-$t$ distribution, denoted by $\textrm{St}(\mu,\omega,\xi,\nu)$, where $(\mu,\omega,\xi,\nu)\in\mathbb{R}\times\mathbb{R}_{+}\times\mathbb{R}\times(2,\infty)$ are location, scale, skewness and degree of freedom parameters respectively. For the existence of the variance of $Z$ we require that $\nu>2$. The distribution of $Z$ has density

[TABLE]

where $z_{S}=(z-\mu)/\omega$, and $f_{T}$ and $F_{T}$ denote, respectively, the density and distribution function of a Student-$t$ random variable with location and scale parameters 0 and 1 respectively and with $\nu$ degrees of freedom. With this notation we have

[TABLE]

where

[TABLE]

with $\Gamma$ being the gamma function. To impose the moment conditions on $Z$ of model (1.2) we also require that

[TABLE]

where $\nu$ and $\xi$ are constrained further to ensure that $\omega$ is a positive real number. The parameter $\xi\in\mathbb{R}$ controls the skewness: $\xi>0$ and $\xi<0$ correspond to right and left skew respectively, and $\xi=0$ to a symmetric distribution. Important special cases arise when: $\xi=0$, $Z$ is the (scaled) Student-$t_{\nu}$ distribution; $\nu\to\infty$, $Z$ is a skew-Normal distribution; $\xi=0$ and $\nu\to\infty$, a standard normal distribution. If not stated below we assume that $Z_{t}$ follows a standard normal distribution.

Table 1 presents values for the key stationarity and extremal properties, $\gamma,\eta,\kappa,\theta_{X^{2}},\theta_{X^{U}},\theta_{X^{L}}$ and $\delta$ for each of models A-E and for Student-$t$, asymmetric Student-$t$ and Gaussian innovations. These values are derived using the numerical methods in the rest of Section 5, and their values are discussed in Sections 5 and 6.

5.2 Evaluation of $\gamma$ and $\eta$

Expression (2.5) suggests using Monte Carlo for the evaluation of $\gamma$, by taking $t$ to be very large. We have also introduced, through Theorems 3.1 and 3.3, two new ways to evaluate $\gamma$, with the latter only applicable once it is known that the process is stationary, and hence it is known that $\gamma<0$. Here we compare these approaches to illustrate the superior computational stability and reliability of our proposed approaches. First we illustrate the methods in Figure 1 for model A which is known, from the discussions in Section 2.1, to be stationary without the requirement of the evaluation of $\gamma$.

Figure 1 (left) shows that $\gamma_{t}$, evaluated using expression (2.5), has serious numerical instabilities for large $t$. All ten independent realisations of $\gamma_{t}$ appear to be converging to roughly the same negative value as $t$ increases. This finding suggests that the limit $\gamma$ is negative and so the process is correctly found to be strictly stationary. But at different random values of large $t$ each replicate stops at a time when the norm for that replicate is calculated as [math] to machine precision. In these cases $\gamma_{t}=-\infty$ for all subsequent $t$ even though $\gamma$ is known to be finite. Hence wrong conclusions about strict stationarity can be reached for model A using this method. This numerical instability for evaluating $\gamma$ does not appear to have been reported. For example, estimates of $\gamma$ using this approach are presented for ARCH(2) processes in Francq and Zakoïan (2010, p. 34/35), however they stop evaluating $\gamma_{t}$, when $t=1000$, which is before we see is the critical failure of numerical evaluation in Figure 1. By increasing $t$ we find similar numerical problems to those experienced for model A in Figure 1.

Theorem 3.1 shows that actually $\gamma=E(\ln\lambda)+\lim_{t\rightarrow\infty}\eta_{t}$ with $\eta_{t}$ is defined by expression (3.2), if condition (3.3) holds. To use this result first we check that condition (3.3) holds. For model A, Figure 1, shows $\frac{1}{t}\ln\lVert C_{t}\rVert$ for the same replicates as in the left panel. This figure indicates that this quantity appears to be converging to [math] as $t\rightarrow 0$, i.e., condition (3.3) appears to hold.

For model A using the formulation of Theorem 3.1 we find that $\gamma=-0.359+0.019=-0.34<0$, which clearly satisfies the strict stationarity condition. We have evaluated $E(\ln\lambda)$ by both numerical integration and using the Monte Carlo approximation $\sum_{i=1}^{t}\ln(\lambda_{i})/t$ for large $t$ where $\lambda_{i}$ is the largest magnitude of the eigenvalue of $\mathbf{A}_{i}$, and obtained identical and stable results once Monte Carlo errors are accounted for. Evaluations of the $\eta_{t}$ term are shown in Figure 1 right panel, for the same 10 replicates as before. The numerical problems are now resolved with convergence similar in all cases, with greater variation between realisations than between iterations over large ranges of $t$. The value we present for $\eta=\lim_{t\rightarrow\infty}\eta_{t}$ is evaluated as the mean over the ten different realisations using when $t=30000$. Although Monte Carlo results are subject to noise, it is possible to obtain any desired level of accuracy by running sufficient replicates when assessing the variability using central limit results (Goldsheid, 1991). Using the results of Theorem 3.2, we can also evaluate $\eta$, but this requires knowledge of the value of $\kappa$. Using the methods of Section 5.6, based on Algorithm 1, for model A we have that $\kappa=2.37$ and so it follows from Theorem 3.2 that $\eta=0.0233$. Thus this shows that trying to evaluate $\eta$ using the Monte Carlo limit, as in Figure 1, is liable to a small numerical error.

As stationarity of model A has been derived we can also evaluate $\gamma$ using the even more numerically reliable method given by Theorem 3.3, which comes from the expectations of two functions of the largest eigenvalue of the matrix $\mathbf{A}$. This evaluation requires the knowledge of $\kappa$. As for the evaluation of $\eta$ above, for model A we use $\kappa=2.37$, and we then obtain that $\gamma=-0.336$.

All of the values of both $\eta$ and $\gamma$ reported in Table 1 are evaluated using the methods based on Theorems 3.2 and 3.3 respectively, but of course they can only be used once stationarity has been determined, or at least there is strong evidence that $\gamma$ may be negative based on using the method based on Theorem 3.1. Table 1 shows that although $\eta=0$ for all GARCH$(1,1)$ processes (confirming Theorem 3.2) this does not hold for any of our GARCH$(p,q)$ processes with $\max(p,q)\geq 2$. Furthermore, in all models, we have $\gamma$ closest to 0 with the Gaussian innovation, then the symmetric $t_{3}$ distribution. There is no obvious pattern in the behaviour of $\eta$ over the factors we explore in Table 1.

5.3 Initialising Algorithm 1

To be able generate realisations from the tail chain (3.10), through Algorithm 2, we first need to generate samples for $\hat{\mathbf{\Theta}}_{0}$ using Algorithm 1. However, to use Algorithm 1 we need to be able to sample from a suitable random variable $\widetilde{\mathbf{\Theta}}_{0}$ on $\mathbb{S}^{p+q}$ with a distribution which is as close as possible to the target limit distribution function $H_{\hat{\mathbf{\Theta}}_{0}}(\mathbf{w})$, so that the rate of convergence of $\widetilde{\mathbf{\Theta}}_{s}\rightarrow^{d}\hat{\mathbf{\Theta}}_{0}$ as $s\rightarrow\infty$ is maximised. From limit (2.10) we have that

[TABLE]

So for large enough $x$, i.e., $x\geq u$ for some high threshold $u$, if we treat this limiting representation as an equality this gives us an initial estimate $H^{(0)}_{\widetilde{\mathbf{\Theta}}}$ of $H_{\hat{\mathbf{\Theta}}_{0}}$. We select $u$ as a high threshold of $R_{t}$ such that limit property (2.10) appears to be well represented, i.e., radial values appear to have a Pareto tail and radial and angular values appear independent.

In practice to obtain $H^{(0)}_{\widetilde{\mathbf{\Theta}}}$ we generate a sample of length $n$ from the required GARCH($p,q$) process and take the empirical distribution of simulated values of $\mathbf{\Theta}_{t}$ given that $R_{t}>u$ after a burn in period of $n_{b}$ i.e.,

[TABLE]

where $\mathbf{1}(F)$ is the indicator function of event $F$ and with $n_{u}=\sum_{j=n_{b}+1}^{n}\mathbf{1}(R_{j}>u)$. As initial particles for Algorithm 1 we use all the realisations of $\mathbf{\Theta}_{t}$ given that $R_{t}>u$, for $t=1,\ldots,n$. We used $n=1.1\times 10^{7}$, and $u$ to be the $99.99\%$ quantile of $R_{t}$ giving $J=10^{3}$ particles, each with equal weight $J^{-1}$.

5.4 Investigation into convergence of Algorithm 1

We illustrate the convergence of Algorithm 1 for models A and C. First consider model C where the true distribution of $\hat{\mathbf{\Theta}}_{0}$, here a scalar, is given by expression (2.12). Hence we can compare the $s$ iteration estimate $\hat{H}^{(s)}_{\widetilde{\mathbf{\Theta}}}(w)$ against the truth $\hat{H}_{\hat{\mathbf{\Theta}}_{0}}(w)$. Figure 2 illustrates this distributional convergence as well as that of the distribution of the $J$ particle weights, $\mathbf{m}^{(s)}$ of expression (3.9) on iteration $s$. First note that the 95% pointwise confidence intervals of the initial estimate $\hat{H}^{(0)}_{\widetilde{\Theta}}$ given in Section 5.3 does not the contain the true target distribution $\hat{H}_{\hat{\mathbf{\Theta}}_{0}}(w)$. Despite our efforts to obtain a good initial guess for Algorithm 1, the statistically significant difference between them is due to the slow convergence of the distribution of $\Pr(\mathbf{\Theta}_{0}<w|R_{0}>u)$ to $H_{\hat{\mathbf{\Theta}}_{0}}(w)$ as $u\rightarrow\infty$.

After one step of Algorithm 1 we have that $H^{(1)}_{\widetilde{\mathbf{\Theta}}}(w)$ is equal to $H_{\hat{\mathbf{\Theta}}_{0}}(w)$ to within visible detection. Further iterations of the algorithm lead to no visible changes in $\hat{H}^{(s)}_{\widetilde{\Theta}}(w)$ for $s>1$. In fact, in this example, we found essentially a perfect convergence after one iteration whatever the initial distribution estimate indicating a unique solution with the algorithm being robust and highly efficient in converging to it. Now focus on the particle weights obtained in the algorithm. Initially, i.e., for $s=0$, all weights are equal $J^{-1}$, but, as Figure 2 shows, within an iteration they have quite a different distribution of weights and that this distribution essentially has converged at $s=2$ to its limit form. Thus Algorithm 1 works exceptionally well in this case where we know the answer. Similar tests over other GARCH(1,1) processes gave identical convergence performances.

Next we assess the convergence of $H^{(s)}_{\widetilde{\mathbf{\Theta}}}(\mathbf{w})$, over $s$, for model A. Here $\hat{\mathbf{\Theta}}_{0}$ is four dimensional and its distribution is not known, so we cannot easily show graphically the convergence of the full joint distribution convergence and even for lower dimensional summaries we can only show the algorithm converges to some limit. Figure 3 illustrates convergence for each of the marginal distributions of $H^{(s)}_{\widetilde{\mathbf{\Theta}}}(\mathbf{w})$ over $s$. Other than for the second component the marginals appear to stabilise to their limit after just one iteration, for that margin it occurs in two iterations. The final panel of Figure 3 similarly shows that the distribution of the particle weights of the particles also converges after two iterations. We also assessed (not shown) the convergence of the dependence structure of $H^{(s)}_{\widetilde{\mathbf{\Theta}}}(\mathbf{w})$ through monitoring how $\mbox{corr}(\widetilde{\vartheta}^{(i)}_{s},\widetilde{\vartheta}^{(j)}_{s})$ converges to $\mbox{corr}(\hat{\vartheta}^{(i)}_{0},\hat{\vartheta}^{(j)}_{0})$, where $\widetilde{\mathbf{\Theta}}=(\widetilde{\vartheta}^{1},\ldots,\widetilde{\vartheta}^{p+q})$. In all cases we found rapid convergence.

We studied a number of other GARCH($p,q$) processes and found excellent convergence of the algorithm, with convergence appearing to occur in $p+q-1$ iterations in all cases with our initialisation method, and convergence to the same value over a range of other initialisation distributions. The theory behind Algorithm 1 indicates that there is a unique solution and that the algorithm will find this, hence our numerical studies support this and show that it works with very high efficiency.

5.5 Assessment of the convergence for the Extremal Index

We show that our solution for $H_{\hat{\mathbf{\Theta}}_{0}}(w)$ gives values of functionals, such as the extremal index, which are consistent with estimates obtained from long-run simulations from the GARCH($p,q$) process. The convergence of Algorithm 1 can be assessed for the cluster functionals by combining outputs from Algorithms 1 and 2. We illustrate this by finding $\theta_{X^{2}}$, the extremal index of the squared process for model A. Using the tail chain $\{\hat{X}_{t};t=0,1,\ldots\}$ we have

[TABLE]

We can estimate the extremal index using the runs estimator $\tilde{\theta}_{X^{2}}(u,m)$, proposed by Smith and Weissman (1994), based on a sample from the GARCH($p,q$) process of length $n$, where

[TABLE]

Figure 4 shows $\theta_{X^{2}}$ together with the runs estimate, based on $n=10^{7}$, for a range of values of $u$ and $m$. The limit value is within the estimated 95% confidence intervals for the runs estimate for all $u$ and both values of $m$ and the values of the runs estimate approach the true value $\theta_{X^{2}}$ as $u$ increases for both $m$ values, although the uncertainty in the estimators increases. This plot also goes to show why it is not possible to derive these extremal features of the GARCH($p,q$) process simply from very long runs as the competing needs of large $u$, for convergence, and large numbers of exceedances of $u$, for numerical stability, makes getting numerically reliable values essentially impossible without taking $n$ to be multiple orders of magnitude larger than here.

5.6 Evaluation of $\kappa$

Basrak and Segers (2009), and subsequent authors, imply that the way to evaluate $\kappa$ is to numerical solve the limiting equation (2.15), although they do not illustrate this. In Figure 1 (left panel) we showed that there are major numerical instabilities in evaluating $\lVert\mathbf{A}_{t}\cdots\mathbf{A}_{1}\rVert$ for large $t$; so in practice it is impossible to solve equation (2.15) directly. In this paper we have discussed three alternative approaches for determining $\kappa$: the algorithm of Janssen (2010); using the formulation for $\kappa$ given by Theorem 3.2; and exploiting Algorithm 1. Here we describe, and illustrate, the relative merits of these methods.

The algorithm of Janssen (2010) is only for bounded innovation variables. From a numerical efficiency perspective it suffers from the critical problem that as it is based on rejection sampling, meaning it can get seriously stuck. Finally, the routine was written in pure R and has a naive initialisation, so it is very slow (taking 2/3 days) to evaluate $\kappa$ to an accuracy of three significant figures even when applied for a GARCH(2,1) model. The speed slows at a cubic rate as the number of the GARCH parameters grows. So this algorithm cannot be used for arbitrary GARCH($p,q$) processes, even with bounded innovations.

Both the new approaches that we present for evaluating $\kappa$ are not restricted by the choice of the GARCH dimensions $p$ and $q$, they apply whether the innovations are bounded or unbounded, and they are relatively much faster as they are coded in C wrapped by R and run in parallel.

The first of our methods is based on the equivalent representation to limiting equation (2.15), i.e., that $\kappa>0$ satisfies $E[(\lambda\exp(\eta))^{\kappa}]=1$ as given by Theorem 3.2. As $\lambda$ can be derived analytically (or, less efficiently, numerically) from $\mathbf{A}$ it remains to find $\eta$ and then $\kappa$ can be found when using either numerical integration or Monte Carlo methods to evaluate the required expectation. We derive an estimate of $\eta$ using the methods presented in Section 5.2. Unfortunately, this approach is not ideal as is illustrated in Figure 1, which shows that the 10 replicates of $\eta_{t}$ are not sufficiently stable and self-consistent in their values at large iterations to accurately deduce the precise value of the limit $\eta$.

We find that this approach only works well for calculating $\kappa$ for models where $\lvert\phi-1\rvert>0.05$ as there is too much sensitivity to the uncertainty of $\eta$ otherwise. Although this is not the ideal way to evaluate $\kappa$, its form gives helpful intuition into what influences $\kappa$. As seen in Section 5.2 we can actually evaluate $\eta$ much more accurately, but that needs $\kappa$ to be found, so that would lead to a circular argument.

Our preferred approach to calculating $\kappa$ is to use Algorithm 1, iterating over $k$ to give $\kappa$, as this provides no numerical problems whatever the dimension of the GARCH ($p,q$) process. Key to the solution is the evaluation of the Monte Carlo estimate $\tilde{\rho}_{k}$ of $\rho_{k}$ in expression (3.7). Figure 5 shows $\tilde{\rho}_{k}$ against $k$ for each of the models A-E. There is clearly a unique solution for $k>0$ to the equation $\tilde{\rho}_{k}=1$, with the values of $k=\kappa$ that solve this equation given in Table 1. However to achieve this we need to reduce the noise in the Monte Carlo estimates $\tilde{\rho}_{k}$ of $\rho_{k}$. For each value of $k$ shown in Figure 5 we used $J=10^{6}$ and evaluated the Monte Carlo integral (3.7) with $10^{4}$ replicates on $Z$ to get $\kappa$ to the required precision. To find $\kappa$, from the curve of $\tilde{\rho}_{k}$, we used an initial grid search coupled with a bisection method.

Figure 5 and Table 1 illustrate that $\phi$, the sum of the meaningful GARCH parameters, has a substantial impact on the value of $\kappa$ with for $\phi>1$, we find $\kappa<1$; when $\phi<1$, then $\kappa>1$; and for $\phi=1$, $\kappa=1$, with the latter consistent with Theorem 3.4. Unfortunately, when $\phi\not=1$ no explicit relationship appears to hold between $\phi$ and $\kappa$, as $\kappa$ changes markedly with the innovation distribution. From Table 1 is can be seen that when $\phi<1$ we have that the shorter the tail of the innovation distribution gives the larger $\kappa$ and hence shorter tails of the GARCH$(p,q)$ marginal distribution; whereas the reverse holds when $\phi>1$; and when $\phi=1$ then $\kappa$ is invariant to the innovation distribution. The case when $\phi>1$ is somewhat surprising as at first thought you would expect that having a heavier tail innovation would result in a heavier tailed GARCH$(p,q)$ process, whereas in fact the opposite occurs.

We finish with empirical diagnostic checks to illustrate that the derived value of $\kappa$ is consistent with the observable tail of the GARCH($p,q$) process. The observable tail can be derived from long run simulations. Specifically we compare the probabilities limiting $\Pr(\hat{X}_{t}^{2}>r\mid\hat{X}_{t}^{2}>1)=r^{-\kappa}$ with the empirical estimate of the probabilities $\Pr(X_{t}^{2}>rx\mid X_{t}^{2}>x)$ for very large $x$, over a range of $r>1$. Figure 6 shows this comparison on a log scale. With the choice of such scaling the true relationship between them has a gradient $\kappa$. The results show that at this far into the distributional tail, and subject to Monte Carlo noise, the empirical distribution is consistent with the limit formulation, and hence also is consistent with the value for $\kappa$ that we have derived.

5.7 Evaluation of $\delta$

In Section 4 we introduced a tail-skewness parameter $\delta$ as a limiting conditional probability (4.1). We have not found any previous discussion on the evaluation of $\delta$, which is an important parameter in the calculation of extremal features of GARCH($p,q$) processes when the innovations are asymmetrically distributed, despite Ehlert et al. (2015) covering this class for GARCH(1,1) processes. A natural starting point to evaluate $\delta$ is to take a long-run simulation from the GARCH($p,q$) process and simply estimate the probability (4.1) empirically for a large enough value of $x$. However, this is likely to be unreliable in practice.

Key to the method we propose is the following expression

[TABLE]

where $F_{\hat{\sigma}_{t}^{2}}$ is the distribution of $\hat{\sigma}^{2}_{t}$, the $(q+1)$the component of the limit vector $\hat{\mathbf{\Theta}}_{t}$ which is defined by limit (2.6).

For the $\textrm{St}(\mu,\omega,\xi,\nu)$ innovation distribution then expression (5.1) becomes

[TABLE]

With a sample of $s_{1},\ldots,s_{m}$ from $\hat{\sigma}_{t}^{2}$ derived using Algorithm 1 we can get a Monte Carlo approximation, to any desired accuracy though the choice of $m$, as follows

[TABLE]

Intuitively it seems as though $\delta$ should be equal to $\delta_{Z}=\lim_{x\rightarrow z_{F}}\Pr(Z_{t}>x\mid|Z_{t}|>x)$, where $z_{F}$ is the upper end point of $|Z_{t}|$. This is not the case though, as knowing that $|X_{t}|$ is large is quite different from knowing that $|Z_{t}|$ is large, as the former also can be achieved with large volatility and small innovations. For $\textrm{St}(\mu,\omega,\xi,\nu)$ we have that $\delta_{Z}$ is given by

[TABLE]

Figure 7 shows how $\delta=\delta(\xi)$ varies with the skew-$t_{3}$ distribution parameter $\xi$ for $\xi\geq 0$; the values of $\delta$ for $\xi<0$ follow due to $\delta$ being symmetrical about $0.5$, i.e., for $\xi<0$ then $\delta(\xi)$ is equal to $1-\delta(|\xi|)$. The figure shows that for a given level of $\xi$, i.e., skewness in the innovation distribution, as $\phi$ increases there is a diminishing level of skewness in the tails of the GARCH$(p,q)$ process as measured by $\delta$. When $\phi=0$ then $\delta=\delta_{Z}$, and this value is seen to be an upper bound for $\delta$ in Figure 7. These results are consistent with intuition as the larger the value of $\phi$ the more the process is driven by the past values of the process and the less the new innovations (and their skewness) matter. Table 1 illustrates how $\delta$ changes over models when $\xi=1$, with $\delta$ typically decreasing as $\phi$ increases, as the persistence of volatility is more important than the innovation structure as $\phi$ increases.

6 Results for the GARCH($p,q$) process

6.1 Extremogram

Figure 8 gives the extremogram $\chi_{X^{2}}(\tau)$ for the squared of GARCH process for models A-D with normal and $t$ innovations with $\nu=3$ degrees of freedom. Firstly notice the impact of the innovation distribution on $\chi_{X^{2}}(\tau)$. In all cases the heavier tailed innovation distribution leads to weaker extremal dependence at all lags. Models C and D, both IGARCH processes (with $\phi=1$), exhibit much slower decay rates in extremal dependence as lag $\tau$ increases than for models A and B with $\phi<1$, with the level of extremal dependence appearing to be strongly related to $\phi$. Furthermore, we see for models B and D that $\chi_{X^{2}}(2)>\chi_{X^{2}}(1)$, with $\chi_{X^{2}}(\tau)$ decaying monotonically for $\tau\geq 2$. The reason for this seems to be that $\beta_{2}>\max(\alpha_{1},\alpha_{2})$ here. From Section 4 we have that the evaluation of $\chi_{X^{U}}(\tau)$ and $\chi_{X^{L}}(\tau)$ is simple from $\chi_{X^{2}}(\tau)$ once $\delta$ is known, which we have from Section 5.7 and Table 1.

An empirical estimate $\tilde{\chi}_{X^{2}}(\tau,u)$ of the extremogram of $\chi_{X^{2}}(\tau)$ based on a sample of length $n$ from a GARCH($p,q$) process is given by

[TABLE]

where $u$ is a threshold. Figure 8 shows $\tilde{\chi}_{X^{2}}(\tau,u)$ for large $n$ and for three threshold choices $u$ corresponding to $0.99,0.999$ and $0.9999$ quantiles of $X_{t}^{2}$. The agreement with limit values $\chi_{X^{2}}(\tau)$ that we have evaluated is very good generally, with the empirical estimates suffering from bias and variance trade-off, as with all threshold based estimates. Model B has the slowest convergence of the empirical estimator, but even here at the highest threshold there is almost perfect overlap between empirical estimates and the true values for all lags $\tau$. In contrast for model C the highest threshold produces the least good estimate, presumably due to its high variance. This gives strong evidence that our evaluation of $\chi_{X^{2}}(\tau)$ is accurate, but it also shows how difficult it is to get accurate values from direct simulations due to different convergence rates from apparently rather similar models.

6.2 Extremal Index

We finish by looking at how the extremal index $\theta_{X^{2}}$, $\theta_{X^{U}}$ and $\theta_{X^{L}}$, of the upper tail of the series $\{X^{2}_{t}\},\{X_{t}\}$ and $\{-X_{t}\}$, change over GARCH($p,q$) processes. Results for each of these characteristics are given for models A-E and three innovation distributions are given in Table 1. Firstly consider the effect of the innovation distribution for a given model on the extremal index of the $\{X^{2}_{t}\}$ process, $\theta_{X^{2}}$ through the associated average size of clusters extremes values, i.e., $1/\theta_{X^{2}}$. With shorter tailed innovations clusters last longer on average and introducing skewness further reduces the mean cluster size. For increasing $\phi$, for $0<\phi\leq 1$, we have increasing average cluster sizes, but that pattern does not follow when $\phi>1$. In all cases, $\min(\theta_{X^{U}},\theta_{X^{U}})\geq\theta_{X^{2}}$, indicating the extremes of the processes $\{X_{t}\}$ and $\{-X_{t}\}$ exhibit less clustering on average than the $\{X^{2}_{t}\}$ process. We have equality, in this inequality, only when $\delta=0$ or $1$, and find that as $\delta$ tends to these limits one or other of the processes $\{X_{t}\}$ and $\{-X_{t}\}$ has similar cluster of extremes events as the $\{X^{2}_{t}\}$ process. For $\xi>0$ more clustering occurs in the upper tail than in the lower tail of the GARCH$(p,q)$ process, with the reverse happening with $\xi<0$.

To give a better idea of how the GARCH($p,q$) parameters affect the extremal index Figure 9 presents a contour plot of $\theta_{X^{U}}$ for GARCH(2,2) process over different key parameters $(\alpha_{1},\alpha_{2},\beta_{1},\beta_{2})$. In particular, in each panel we hold fixed two parameters and contour over the other two key parameters of a GARCH(2,2) process. Figure 9, shows that broadly $\theta_{X^{U}}$ decreases, i.e., average cluster sizes increase, with increasing $\phi$, up to the boundary case of an IGARCH(2,2) model. Consequently contours are near linear in the parameters. Unless $\max(\beta_{1},\beta_{2})$ is large then small values of $\alpha_{1}$ and $\alpha_{2}$ tend to lead to very limited clustering. This seems logical given the GARCH formulation, as small $\alpha_{1}$ and $\alpha_{2}$ mean that the effect of the large $X_{t}^{2}$ value can have limited impact of the subsequent volatilities, so without $\beta_{1}$ and $\beta_{2}$ being large, to pick up the momentum of the evolution of the event, the large event is very likely to die out rapidly. We also have that we obtain stronger extremal dependence with larger values of the pair $(\alpha_{2},\beta_{2})$ than for $(\alpha_{1},\beta_{1})$ for equal values of $\phi$, as seen by values of $\theta_{X^{U}}$ being smaller in Figure 9 right panel by comparison to the left panel. The reason for this is the stronger effect of $\beta$ coefficients compared to that of the $\alpha$ as the persistency in the volatility is translated into a stronger extremal dependence.

7 Discussion

The new theory and methods we present extend to assessing the strict stationarity and extremal properties for a much broader class of stochastic recurrence equations. Specifically, for any process

[TABLE]

where $\mathbf{A}_{t}$ and $\mathbf{B}_{t}$ are stochastic, independent and identically distributed, sequences of matrices and vectors respectively, satisfying the conditions of Kesten (1973), then our paper gives methods which: help assess stationarity of the process; determine its tail index of regular variation for the marginals distributions of $\mathbf{Y}_{t}$; determine ways to simulate from the spectral measure from the limiting joint distribution of $\mathbf{Y}_{t}$ given some norm, $\lVert\mathbf{Y}_{t}\rVert$, tends to infinity; and are able to generate a wide range of properties, such as the extremal index of any marginal of the process $\mathbf{Y}_{t}$ from its tail chain. Previously none of these properties could be derived due to a combination of a lack of appropriate numerically stable algorithms. The methods presented here overcome these limitations and provide a broad toolbox of numerically robust approaches to derive the extremal analysis of a wide class of stochastic recurrence equations including all GARCH($p,q$) processes with bounded and unbounded innovation variables.

Appendix A Proof of Theorems

Proof of Theorem 3.1.

First rewrite the product of independent matrices as

[TABLE]

Thus

[TABLE]

so as $t\rightarrow\infty$

[TABLE]

hence $\gamma=E(\ln\lambda)+\eta$. ∎

Proof of Theorem 3.2.

With the same notation as for Theorem 3.1 we have that

[TABLE]

Under condition (3.3) we have that there exists a sequence $s_{t}>0$ such that $s_{t}\rightarrow 0$ and $s_{t}t\rightarrow\infty$ as $t\rightarrow\infty$

[TABLE]

for all $t$. Combining these inequalities with expression (A.2) we obtain that

[TABLE]

Focusing first on the upper bound we have that

[TABLE]

where the second step comes from $\lambda_{i}$, $i=1,\ldots,t$ being independent and identically distributed. Hence, as $t\rightarrow\infty$

[TABLE]

By an identical argument, the lower bound is found to be asymptotically equal to the upper bound. Hence

[TABLE]

so this limit is equal to 0, as required by condition (2.15), only when $\kappa$ satisfies $E[(\lambda\exp(\eta))^{\kappa}]=1$.

∎

Proof of Theorem 3.4.

First suppose that $\kappa=1$, then from property (2.13) we have that $E(\lVert\mathbf{A}\hat{\mathbf{\Theta}}_{t}\rVert)=1$, where $\mathbf{A}$ is independent of $\hat{\mathbf{\Theta}}_{t}$. Then, by considering the vector $\mathbf{Y}_{t}$, we see that marginally the first $q$ and the last $p$ components of $\hat{\mathbf{\Theta}}_{t}$ each have identical marginal distributions. Furthermore, as $E(Z_{t}^{2})=1$, it follows that all $E(\hat{\vartheta}_{t}^{(i)})=1/(p+q)$ for all $i=1,\ldots,p+q$, where recall $\hat{\mathbf{\Theta}}_{t}=(\hat{\vartheta}_{t}^{(1)},\ldots,\hat{\vartheta}_{t}^{(p+q)})$. We also have that

[TABLE]

Let $a_{i,j}$ denote the $(i,j)$th element of $\mathbf{A}$, using the independence of $a_{i,j}$ from $\hat{\vartheta}_{t}^{(j)}$ then

[TABLE]

Thus $E(\lVert\mathbf{A}\hat{\mathbf{\Theta}}_{t}\rVert)=1$ only when $\sum_{i=1}^{q}\alpha_{i}+\sum_{i=1}^{p}\beta_{i}=1$, i.e., when the process is IGARCH($p,q$). The argument is simply reversed giving that $\kappa=1$ for any IGARCH($p,q$) process. ∎

Acknowledgements

We would like to thank Feridun Turkman for encouraging and helpful discussions.