# Comptes Rendus
Mathématique

Christopher S. Withers and Saralees Nadarajah

The distribution of the maximum of an ARMA(1, 1) process

Volume 358, issue 8 (2020), p. 909-916

Published online: 3 December 2020

https://doi.org/10.5802/crmath.111

<figure><img src="image_1.png" /></figure>

This article is licensed under the
CREATIVE COMMONS ATTRIBUTION 4.0 INTERNATIONAL LICENSE.
http://creativecommons.org/licenses/by/4.0/

<figure><img src="image_2.png" /></figure>

Les Comptes Rendus. Mathématique sont membres du
Centre Mersenne pour l'édition scientifique ouverte
www.centre-mersenne.org

e-ISSN : 1778-3569
Probabilities / Probabilités

# The distribution of the maximum of an ARMA(1, 1) process

Christopher S. Withers<sup>a</sup> and Saralees Nadarajah<sup>b</sup>

<sup>a</sup> Callaghan Innovation, Lower Hutt, New Zealand

<sup>b</sup> University of Manchester, Manchester M13 9PL, UK

E-mails: kit.withers@gmail.com, mbbsssn2@manchester.ac.uk

**Abstract.** We give the cumulative distribution function of $M_n = \max(X_1, \dots, X_n)$, the maximum of a sequence of $n$ observations from an ARMA(1, 1) process. Solutions are first given in terms of repeated integrals and then for the case, where the underlying random variables are absolutely continuous. The distribution of $M_n$ is then given as a weighted sum of the $n$th powers of the eigenvalues of a non-symmetric Fredholm kernel. The weights are given in terms of the left and right eigenfunctions of the kernel.

These results are large deviations expansions for estimates, since the maximum need not be standardized to have a limit. In fact, such a limit need not exist.

Manuscript received 16th August 2020, accepted 3rd September 2020.

## 1. Introduction

There has been little work with respect to extreme value theory for ARMA (autoregressive-moving average) processes. The authors are aware only of the work of Rootzén [2, 3]. Both these papers provide the limiting extreme value distributions or assume that the errors come from a specific class. We are aware of no work giving the *exact* distribution of the maximum of ARMA processes.

Knowing the exact distribution of the maximum is important. It can lead to improved point estimation, improved interval estimation and improved prediction. In the extreme value literature, the exact distribution is approximated by a limiting distribution with the latter assumed to exist as the sample size approaches infinity. The derivation of the exact distribution for ARMA processes has not been attempted in the literature.

This note continues the application of a powerful new method for obtaining the *exact* distribution of extremes of $n$ correlated observations as weighted sums of $n$th powers of certain eigenvalues. The method was first illustrated for an autoregressive process of order one in Withers and Nadarajah [4] and a moving average of order one in Withers and Nadarajah [5].
Let {$e_i$} be independent and identically distributed random variables from some cumulative distribution function (cdf) $F$ on $R$. Let $f(x)$ denote the probability density function (pdf) with respect to Lebesque measure. We consider the ARMA process of order (1, 1),

$$X_i - rX_{i-1} = e_i + se_{i-1}.$$

We restrict ourselves to the most important case $r > 0$. (When this condition does not hold the method can be adapted as done in Withers and Nadarajah [4].) In Section 2, we give expressions for the cdf of the maximum

$$M_n = \max(X_1, \dots, X_n), \quad n \ge 1,$$

in terms of repeated integrals. This is obtained via the recurrence relationship

$$G_n(\mathbf{y}) = \mathcal{K} G_{n-1}(\mathbf{y}), \quad \mathbf{y} = (y_0, y_1), \quad n \ge 1, \qquad (1)$$

where

$$G_n(\mathbf{y}) = P(M_n \le x, X_n \le y_0, e_n \le y_1), \qquad (2)$$

$$\mathcal{K} h(\mathbf{y}) = r \int_{a_{\mathbf{y},z}}^{\infty} dz_0 f(y_{0,x} - r z_0 - sz_1) h(z_0, dz_1), \qquad (3)$$

$$g_{\mathbf{y}}(t) = (y_{0,x} - t)/r, \quad y_{0,x} = \min(y_0, x), \quad a_{\mathbf{y},z} = (y_{0,x} - sz_1 - y_1)/r,$$

and dependency on $x$ is suppressed. So, $\mathcal{K}$ is a linear integral operator depending on $x$. For (1) to work at $n = 1$ we define $M_0 = -\infty$ so that

$$G_0(\mathbf{y}) = P(X_0 \le y_0, e_0 \le y_1). \qquad (4)$$

Similarly,

$$G_1(\mathbf{y}) = P(X_1 \le y_{0,x}, e_1 \le y_1) = G_0(y_{0,x}, y_1).$$

In Section 3, we consider the case when $F$ is absolutely continuous. In this case we show that corresponding to $\mathcal{K}$ is a Fredholm kernel $K(\mathbf{y}, \mathbf{z})$. We give a solution in terms of its eigenvalues and eigenfunctions. This leads easily to the asymptotic results stated in the abstract. However, there are two problems: the kernel is a generalized function and numerical solution by direct Gaussian quadrature fails. In Section 4, we show that these problems can be avoided by using the iterated Fredholm kernel $K_2(\mathbf{y}, \mathbf{z})$.

Our expansions for $P(M_n \le x)$ for fixed $x$ are large deviation results. If $x$ is replaced by $x_n$ such that $P(M_n \le x_n)$ tends to the generalized extreme value cdf, then the expansion still holds, but not the asymptotic expansion in terms of a single eigenvalue, since this may approach 1 as $n \to \infty$.

For $a$ and $b$ functions on $R^2$, we set

$$\int a = \int a(\mathbf{y})d\mathbf{y} = \int_{R^2} a(\mathbf{y})d\mathbf{y}$$

and similarly for $\int ab$. Let $I\{A\}$ denotes the indicator function; that is, $I\{A\} = 1$ if $A$ is true and $I\{A\} = 0$ if $A$ is false.

## 2. Solutions using repeated integrals

**Theorem 1.** We have $G_n$ of (2) satisfying the recurrence relation (1) in terms of the integral operator $\mathcal{K}$ of (3).

**Proof.** For $n \ge 1$, $G_n$ of (2) satisfies

$$\begin{align*}
G_n(\mathbf{y}) &= P(M_{n-1} \le x, X_n \le y_0, rX_{n-1} + e_n + se_{n-1} \le y_{0,x}, e_n \le y_1) \\
&= E P(M_{n-1} \le x, X_{n-1} \le g_{\mathbf{y}}(se_{n-1} + e_n) | e_n) I\{e_n \le y_1\} \\
&= \mathcal{K} G_{n-1}(\mathbf{y}).
\end{align*}$$

This ends the proof. □
Our goal is to determine $u_n = P(M_n \le x) = G_n(\infty)$, where $\infty = (\infty, \infty)$.

**Theorem 2.** Set

$$a_n = [\mathcal{H}^n G_0(\mathbf{y})]_{\mathbf{y}=\infty}, \quad n \ge 0,$$ 
$$u_n = a_n, \quad n \ge 0. \tag{5}$$

**Proof.** By Theorem 1, for $n \ge 0$,

$$G_n(\mathbf{y}) = \mathcal{H}^n G_0(\mathbf{y}).$$

Putting $\mathbf{y} = \infty$ gives (5).

For example,

$$u_0 = a_0 = 1, \quad P(X_1 \le x) = u_1 = a_1 = E I\{e_0 \le y_1\} \int G_0(g_y \mathbf{z}_1 + e_0) \, dz_1.$$

**3. The case of F absolutely continuous**

Our solution Theorem 2 does not tell us how $u_n$ behaves for large $n$. Also calculating $a_n$ requires repeated integration. Here, we give another solution that overcomes these problems, using Fredholm integral theory given in Appendix A of Withers and Nadarajah [5], referred to below as the appendix.

**Theorem 3.** Let $\mathbf{z} = (z_0, z_1)$ in $\mathbb{R}^2$. Suppose that F is absolutely continuous with pdf f and that $h(\mathbf{z}) \to 0$ as $z_1 \to \infty$. Let $\delta(\cdot)$ denote the Dirac delta function on $\mathbb{R}$. Set

$$b_{\mathbf{y}}(\mathbf{z}) = f(y_{0,x} - r z_0 - s z_1), \quad \gamma_{\mathbf{y}}(z_0) = (y_{0,x} - y_1 - r z_0) / s, \\ A_{\mathbf{y},\mathbf{z}} = \{r z_0 + s z_1 > y_{0,x} - y_1\}, \\ C_1(\mathbf{y}, \mathbf{z}) = \delta(z_1 - \gamma_{\mathbf{y}}(z_0)) b_{\mathbf{y}}(\mathbf{z}), \quad C_2(\mathbf{y}, \mathbf{z}) = I\{A_{\mathbf{y},\mathbf{z}}\} \frac{d}{dz_1} b_{\mathbf{y}}(\mathbf{z}), \\ C(\mathbf{y}, \mathbf{z}) = \sum_{j=1}^2 C_j(\mathbf{y}, \mathbf{z}), \quad K(\mathbf{y}, \mathbf{z}) = -rC(\mathbf{y}, \mathbf{z}).$$

Then we can write (3) in the form

$$\mathcal{H} r(\mathbf{y}) = \int K(\mathbf{y}, \mathbf{z}) r(\mathbf{z}) \, d\mathbf{z}. \tag{6}$$

**Proof.** Set $c_{\mathbf{y}}(\mathbf{z}) = I\{A_{\mathbf{y},\mathbf{z}}\} b_{\mathbf{y}}(\mathbf{z})$. Then $\frac{d}{dz_1} c_{\mathbf{y}}(\mathbf{z}) = C(\mathbf{y}, \mathbf{z})$ and

$$\mathcal{H} h(\mathbf{y}) = r \iint dz_0 c_{\mathbf{y}}(\mathbf{z}) h(z_0, dz_1) = -r \int h(\mathbf{z}) \frac{d}{dz_1} c_{\mathbf{y}}(\mathbf{z}) dz,$$

integrating by parts. This ends the proof. □

Although $K(\mathbf{y}, \mathbf{z})$ is a generalized function, it satisfies

**Theorem 4.** For $r, |s| \ne 1$,

$$0 < \iint K(\mathbf{y}, \mathbf{z}) K(\mathbf{z}, \mathbf{y}) dy dz < \infty.$$

**Proof.** Note that

$$\iint C(\mathbf{y}, \mathbf{z}) C(\mathbf{z}, \mathbf{y}) dy dz = \sum_{i,j=1}^{2} \alpha_{i,j},$$

where

$$\alpha_{i,j} = \iint C_i(\mathbf{y}, \mathbf{z}) C_j(\mathbf{z}, \mathbf{y}) dy dz.$$
Note $\alpha_{1,1}$ involves two delta functions, so the four integrations over $y_0, y_1, z_0, z_1$ reduce to two
over $y_0, z_0$ at $z_1 = \gamma_y(z_0), y_1 = \gamma_z(z_0)$, that is at $y_1 = y_1^*, z_1 = z_1^*$, where

$$
y_1^* = [y_{0,x} - r z_0 - s (z_{0,x} - r y_0)] / (1 - s^2)
$$

and

$$
z_1^* = [z_{0,x} - r y_0 - s (y_{0,x} - r z_0)] / (1 - s^2).
$$

So,

$$
\alpha_{1,1} = 2 \iint_{y_0 < z_0} [K(\mathbf{y}, \mathbf{z}) K(\mathbf{z}, \mathbf{y})]_{y_1=y_1^*, z_1=z_1^*} d\mathbf{y} d\mathbf{z} = 2 (J_1 + J_2 + J_3),
$$

where $J_1, J_2$ and $J_3$ integrate over $A_1 = \{x < y_0 < z_0\}, A_2 = \{y_0 < x < z_0\},$ and $A_3 = \{y_0 < z_0 < x\}$,
respectively. A transformation of variables gives

$$
J_i = b_i \int \int_{A_i} f(\mathbf{u}) f(\mathbf{v}) \, \mathrm{d}u \mathrm{d}v,
$$

where $b_1 = 1/b_2 = (1-s^2)/r^2$ and $b_3 = (1-s^2)/(1-r^2)$. The other $\alpha_{i,j}$ can be dealt with similarly. $\square$

This theorem implies that $K(\mathbf{y},\mathbf{z})$ is a (non-symmetric) Fredholm kernel with respect to Lebesgue measure, allowing the Fredholm theory of the appendix to be applied, in particular the functional forms of the Jordan form and singular value decomposition.

Let {$\lambda_j, r_j, l_j : j \ge 1$} be the eigenvalues and associated right and left eigenfunctions of $\mathcal{K}$ ordered so that $|\lambda_j| \ge |\lambda_{j+1}|$. If {$\lambda_j$} are real then {$r_j, l_j$} can be taken as real. By the appendix referred to, these satisfy

$$
\mathcal{K} r_j(\mathbf{y}) = \lambda_j r_j(\mathbf{y}), \quad \bar{l}_j(\mathbf{z}) \mathcal{K} = \lambda_j \bar{l}_j(\mathbf{z}), \quad \int r_j \bar{l}_k = I \{j=k\}, \qquad (7)
$$

where $\bar{\zeta}$ is the complex conjugate of $\zeta$, and

$$
\bar{l}(\mathbf{z}) \mathcal{K} = \int \bar{l}(\mathbf{y}) K(\mathbf{y}, \mathbf{z}) dy.
$$

So, {$r_j(\mathbf{y}), l_k(\mathbf{y})$} are biorthogonal functions with respect to Lebesgue measure.

We now assume that $K(\mathbf{y}, \mathbf{z})$ has diagonal Jordan form. (This holds, for example, when the eigenvalues are distinct. This will generally be the case in applications.) The functional equivalent of the Jordan form is, by (3.6) of Withers and Nadarajah [4],

$$
K(\mathbf{y}, \mathbf{z}) = \sum_{j=1}^{\infty} \lambda_j r_j(\mathbf{y}) \bar{l}_j(\mathbf{z}).
$$

This implies that

$$
K_n(\mathbf{y}, \mathbf{z}) = \mathcal{K}^{n-1} K(\mathbf{y}, \mathbf{z}) = \sum_{j=1}^{\infty} \lambda_j^n r_j(\mathbf{y}) \bar{l}_j(\mathbf{z}), \quad (8)
$$

where $\mathcal{K}^n$ is the operator corresponding to the iterated kernel $K_n(\mathbf{y}, \mathbf{z})$. By (A.8) of Withers and Nadarajah [5] with $\mu$ Lebesgue measure on $R^2$, if $\mathcal{K}G$ is in $L_2(R^2)$ then

$$
\mathcal{K}^n G(\mathbf{y}) = \sum_{j=1}^{\infty} B_j(G) r_j(\mathbf{y}) \lambda_j^n, \quad n \ge 1, \qquad (9)
$$

where

$$
B_j(G) = \int_{R^2} G \bar{l}_j.
$$

Putting $\mathbf{y} = \infty$ and $G = G_0$ in (9) gives

**Theorem 5.** For $B_j$ of (9) and $n \ge 1$,

$$
a_n = \sum_{j=1}^{\infty} r_j(\infty) B_j(G_0) \lambda_j^n.
$$
**Corollary 6.** Suppose that the eigenvalue $\lambda_1$ of largest magnitude has multiplicity M. For $n \ge 1$,

$$a_n = B(G_0) \lambda_1^n (1 + \epsilon_n),$$

where $\epsilon_n \to 0$ exponentially as $n \to \infty$ and

$$B(G_0) = \sum_{j=1}^{M} r_j(\infty) B_j(G_0).$$

So, for $n \ge 1$, by (5), $u_n = B(G_0) \lambda_1^n (1 + \epsilon_n)$.

Unfortunately, we cannot use the method of Withers and Nadarajah [6] for the numerical solution of the equations for the eigenvalues and eigenfunctions. For example, the first equation in (7) for $r(\mathbf{y}) = r_j(\mathbf{y})$ at $\lambda = \lambda_j$ can be written $-\lambda r(\mathbf{y})/r = c_1(\mathbf{y}) + c_2(\mathbf{y})$, where

$$c_i(\mathbf{y}) = \int C_i(\mathbf{y}, \mathbf{z}) r(\mathbf{z}) d\mathbf{z}.$$

Suppose that we use Gaussian quadrature

$$\int_{R^2} a(\mathbf{z}) d\mathbf{z} \approx \sum_{j=1}^{q} w_j a(\mathbf{z}_j), \qquad (10)$$

where $\{\mathbf{z}_1, \dots, \mathbf{z}_q\}$ are given points in $R^2$ and $\{w_1, \dots, w_q\}$ are given weights. Then

$$c_2(\mathbf{y}) \approx \sum_{j=1}^{q} w_j C_2(\mathbf{y}, \mathbf{z}_j) r(\mathbf{z}_j).$$

However,

$$c_1(\mathbf{y}) = \int r(z_0, \gamma_y(z_0)) dz_0$$

only has one single integral, so would need a different approximation, say

$$\int_R a(z_0) dz_0 \approx \sum_{j=1}^{q'} w'_j a(z'_j),$$

where $\{z'_1, \dots, z'_{q'}\}$ are given points in $R$ and $\{w'_1, \dots, w'_{q'}\}$ are given weights. This gives

$$c_1(\mathbf{y}) \approx \sum_{j=1}^{q'} w'_j r(z'_j, \gamma_y(z'_j)).$$

Putting $\mathbf{y} = \mathbf{z}_j$ now gives the system of equations

$$-\lambda r(\mathbf{z}_i) / r \approx \sum_{j=1}^{q'} w'_j r(z'_j, \gamma_{\mathbf{z}_i}(z'_j)) + \sum_{j=1}^{q} w_j C_2(\mathbf{z}_i, \mathbf{z}_j) r(\mathbf{z}_j),$$

that is

$$-\lambda \mathbf{r} \approx \boldsymbol{\theta} + \mathbf{C}\mathbf{r},$$

where $\mathbf{r}$ has $i$th element $r(\mathbf{z}_i) / r$, $\mathbf{C}$ has $(i, j)$th element $w_j C_2(\mathbf{z}_i, \mathbf{z}_j)$, but $\boldsymbol{\theta}$ has $i$th element

$$\sum_{j=1}^{q'} w'_j r(z'_j, \gamma_{\mathbf{z}_i}(z'_j)),$$

which is not a multiple of $\mathbf{r}$.

Differentiating the first equation in (7) for $r(\mathbf{y}) = r_j(\mathbf{y})$ at $\lambda = \lambda_j$ gives the differential-integral equation for $r_{.1}(\mathbf{y}) := \frac{d}{dy_1} r(\mathbf{y})$,

$$\lambda r_{.1}(\mathbf{y}) / f(y_1) = \int r_{.1}((y_{0,x} - sz_1 - y_1) / r, z_1) dz_1.$$
## 4. A numerical solution

In the last section we saw that the kernel $K(y, z)$ is a generalized function and that numerical solution by direct Gaussian quadrature fails. Here, we show how to get around these problems by using the iterated Fredholm kernel $K_2(y, z)$ given by (8).

By (6),

$$K_2(\mathbf{y}, \mathbf{z}) = r^2 \sum_{j,k=1}^{2} C_{j,k}(\mathbf{y}, \mathbf{z}),$$

where

$$C_{j,k}(\mathbf{y}, \mathbf{z}) = \int C_j(\mathbf{y}, \mathbf{t}) C_k(\mathbf{t}, \mathbf{z}) \mathrm{d}\mathbf{t}.$$

We first show that these are ordinary functions, not generalized functions. This is clearly true for $j = k = 2$. Also $C_{1,1}(\mathbf{y}, \mathbf{z}) = b_{\mathbf{y}}(\mathbf{T}) b_{\mathbf{T}}(\mathbf{z})$ at $\mathbf{T} = \mathbf{T}(\mathbf{y}, \mathbf{z})$ given by $T_1 = \gamma_{\mathbf{y}}(T_0)$ and $z_1 = \gamma_{\mathbf{T}}(z_0)$. Eliminate $T_1$ using

$$y_{0,x} - y_1 - rT_0 = sT_1 = s(T_{0,x} - rz_0 - sz_1). \quad (11)$$

So, $rT_0 + sT_{0,x} = y_{0,x} - y_1 + s(rz_0 + sz_1) = a$ say. When $r > 0$, $r + s > 0$ this has a unique solution

$$T_0 = \min[a/(r+s), (a-sx)/r].$$

Note that $T_1$ is then given by (11). Also

$$C_{1,2}(\mathbf{y}, \mathbf{z}) = \int [b_{\mathbf{y}}(\mathbf{t}) C_2(\mathbf{t}, \mathbf{z})]_{t_1 = \gamma_{\mathbf{y}}(t_0)} \mathrm{d}t_0,$$

$$C_{2,1}(\mathbf{y}, \mathbf{z}) = \int [C_2(\mathbf{y}, \mathbf{t}) b_{\mathbf{t}}(\mathbf{z})]_{z_1 = \gamma_{\mathbf{t}}(z_0)} \mathrm{d}t_0 = \int [C_2(\mathbf{y}, \mathbf{t}) b_{\mathbf{t}}(\mathbf{z})]_{t_1 = t_{0,x} - r z_0 - s z_1} \mathrm{d}t_0.$$

This gives $K_2(\mathbf{y}, \mathbf{z})$ as an ordinary function. It follows that convergence in (8) holds and Theorem 5 holds for $n \ge 2$ although perhaps convergence does not hold at $n = 1$. Iterations can be done using:

$$a_{2n} = [\mathcal{K}^{2n} G_0(\mathbf{y})]_{\mathbf{y}=\infty}, \quad a_{2n+1} = [\mathcal{K}^{2n} G_1(\mathbf{y})]_{\mathbf{y}=\infty}, \quad n \ge 0.$$

Also $G_1(\mathbf{y}) = \mathcal{K} G_0(\mathbf{y})$ is an ordinary function since by (6)

$$\mathcal{K} G(\mathbf{y}) = -r \sum_{j=1}^{2} \mathcal{C}_j G(\mathbf{y}), \quad \mathcal{C}_j G(\mathbf{y}) = \int C_j(\mathbf{y}, \mathbf{z}) G(\mathbf{z}) \mathrm{d}\mathbf{z},$$

and for $G(\mathbf{y})$ an ordinary function, $\mathcal{C}_j G(\mathbf{y})$ is also an ordinary function for $j = 1, 2$. This is obvious for $j = 2$. It is true for $j = 1$ since

$$\mathcal{C}_1 G(\mathbf{y}) = \int [b_{\mathbf{y}}(\mathbf{z}) G(\mathbf{z})]_{z_1 = \gamma_{\mathbf{y}}(z_0)} \mathrm{d}z_0.$$

Since

$$0 < \iint K_2(\mathbf{y}, \mathbf{z}) K_2(\mathbf{z}, \mathbf{y}) dy dz < \infty,$$

(8) holds for $n = 2, 4, \dots$. So, the left and right eigenfunctions $r(\mathbf{y}) = r_j(\mathbf{y})$ and $\bar{l}(\mathbf{y}) = \bar{l}_j(\mathbf{y})$ for the eigenvalue $\theta = \lambda_j^2$ of $K_2(\mathbf{y}, \mathbf{z})$ satisfy

$$\mathcal{K}^2 r(\mathbf{y}) = \theta r(\mathbf{y}), \quad \bar{l}(\mathbf{z}) \mathcal{K}^2 = \theta \bar{l}(\mathbf{z}), \quad \int r \bar{l} = 1.$$

By the Gaussian quadrature approximation, (10), we can approximate this as

$$\mathbf{K}_2 \mathbf{r} \approx \theta \mathbf{r}, \quad \mathbf{K}_2' \mathbf{l} \approx \bar{\theta} \mathbf{l},$$

where $\mathbf{K}_2$ is the $q \times q$ matrix with $(i, j)$ element $w_j \mathbf{K}_2(\mathbf{z}_i, \mathbf{z}_j)$, and $\mathbf{r}$ and $\mathbf{l}$ are the $q$-vectors with $j$th elements $w_j r(\mathbf{z}_j)$ and $w_j l(\mathbf{z}_j)$. So, the first $q$ eigenvalues and right and left eigenfunctions of $K_2(\mathbf{y}, \mathbf{z})$ can be approximated by the eigenvalues and right and left eigenvectors of $\mathbf{K}_2$ standardized so that they are biorthonormal.
So, finally we obtain the distribution of $M_n$ to be

$$u_n = a_n \approx \hat{a}_n, \quad n \ge 2,$$

where

$$\hat{a}_{2n+i} = \sum_{j=1}^{q} \hat{\theta}_{j}^{n} \hat{r}_{j,q} \hat{B}_{j}(G_{i}), \quad n \ge 1, \quad i = 0,1$$

and

$$\hat{B}_{j}(G_{i}) = \sum_{k=1}^{q} w_{k} \bar{l}_{j,k} G_{i}(z_{k}),$$

where $\hat{\theta}_j$ is the jth eigenvalue of $K_2$, and $\hat{r}_{j,k}$ and $\bar{l}_{j,k}$ are the kth components of its left and right eigenvectors.

For a more precise result one can let $q$ increase to $\infty$ with $n$ and use known expressions for the remainder in the Gaussian approximation. Compare with equation (25.4.9) of Abramowitz and Stegun [1].

<figure><img src="image_3.png" /><figcaption>**Figure 1.** Plot of $u_n = P(M_n < x)$ versus $x$ for $n = 1000$, $s = 1$ and $r = 0.1, 0.2, \dots, 0.8, 0.9$ when $G_0$ is a product of two independent standard normal cdfs. The curves from the left to right correspond to increasing values of $r$.</figcaption></figure>

Using Theorem 2, Theorem 3, Theorem 5 and the numerical tools developed above one can calculate $u_n = P(M_n < x)$ for absolutely continuous cdfs $F$ and $H$. Figures 1 and 2 show plots of $u_n$ for $n = 1000$, $s = 1, 5$ and $r = 0.1, 0.2, \dots, 0.8, 0.9$ when $G_0$ is a product of two independent standard normal cdfs. In each figure, the distribution of $u_n$ becomes less dominant as $r$ increases from 0 to 1.
<figure><img src="image_4.png" /><figcaption>**Figure 2.** Plot of $u_n = P(M_n < x)$ versus $x$ for $n = 1000$, $s = 5$ and $r = 0.1, 0.2, \dots, 0.8, 0.9$ when $G_0$ is a product of two independent standard normal cdfs. The curves from the left to right correspond to increasing values of $r$.</figcaption></figure>

References

[1] M. Abramowitz, I. A. Stegun (eds.), *Handbook of Mathematical Functions*, Applied Mathematics Series, vol. 55, U.S. Department of Commerce, National Bureau of Standards, 1964.

[2] H. Rootzén, “The rate of convergence of extremes of stationary normal sequences”, *Adv. Appl. Probab.* **15** (1983), p. 54-80.

[3] ——, “Extreme value theory for moving average processes”, *Ann. Probab.* **14** (1986), p. 612-652.

[4] C. S. Withers, S. Nadarajah, “The distribution of the maximum of a first order autoregressive process: The continuous case”, *Metrika* **74** (2011), no. 2, p. 247-266.

[5] ——, “The distribution of the maximum of a first order moving average: The continuous case”, *Extremes* **17** (2014), no. 1, p. 1-24.

[6] ——, “The distribution of the maximum of a second order autoregressive process: The continuous case”, 2020, *Technical Report*, Department of Mathematics, University of Manchester, UK.