(Non-limit) distribution of maxima from different univariate, discrete and stationary time series Motivation: I'm currently studying the convergence of maxima from simulated time series to max-stable distributions, and in order to do so, I want to better understand the penultimate distribution of such maxima in non-limit cases. I tried to read Leadbetter et al. (1983) "Extremes and related properties of random sequences and processes", but did not get too far.
Starting point: I understand that maxima $M_n$ from a sequence of n iid r.v. $X_1, X_2, ..., X_n$ are distributed according to $P(M_n\leq x)=P(X_1\leq x, X_2\leq x,...,X_n\leq x))=F^n(x)$. Thus I manage to derive the penultimate distribution of maxima from a white-noise process with Gaussian innovations.
Question 0 (New): Is $F^n(x)$ even the penultimate distribution of the Gaussian, or is this a first misunderstanding? According to Cohen (1982), the penultimate distribution of maxima from a sequence of Gaussians is the Type III Extreme Value distribution / GEV with shape parameter $\xi<0$.
Question 1: Is there are closed-form distribution of maxima from an ARMA(2,2) process with Gaussian innovations? I found something on arxiv on maxima of ARMA(1,1) processes. But if there is not such distribution, ...
Question 2: ... what is the unconditional (i.e. $P(X_i<x)$, not $P(X_i<x | x_{i-1})$) distribution of ARMA(2,2)-generated r.v., such that one could employ $F^n(x)$ (either ignoring that fact that these r.v. are not independent, or adjusting $n$ such that it accounts for the dependence, e.g. such that $n^*=n/k$, where $k$ is the lag where the ACF/PACF are close enough to zero).
I'm thankful for both every hint you can give me, but also for pointing out if my questions are nonsense.
 A: Question 2. You want the stationary distribution of the
Gaussian AR process $X_t$
$$
    (1 - \phi_1 B - \dots - \phi_p B^p) \, X_t
    = (1 + \theta_1 B + \dots + \theta_q B^q) \,\varepsilon_t
  $$
for the special case $p=q=2$.  This distribution a.k.a. the
invariant distribution is a Gaussian distribution: its mean
$\mu_X$ and sd $\sigma_X$ can be found. In the case where
$\varepsilon_t$ has mean zero we have $\mu_X = 0$ and $$ \sigma_X^2
  = \sigma_\zeta^2 \sum_{k \geq 0} \psi_k^2 $$ where the coefficients $\psi_k$
are the "psi weights" of the $\text{MA}(\infty)$ representation $X_t
  = \sum_{k \geq 0} \psi_k \zeta_{t-k}$ where $\zeta_t$ is a Gaussian
white noise. The "psi weights" are computed by many R packages.
An alternative derivation uses the ARMA model in state-space form:
the state equation defines a vector AR(1) process with $r:= \max\{p,
  \, q + 1\}$. We can assume that the observed series is the first
component of the state $\boldsymbol{\alpha}_t$ in the model
\begin{align*}
     \boldsymbol{\alpha}_t &= \mathbf{T} \boldsymbol{\alpha}_{t-1} +
     \boldsymbol{\eta}_t\\
     X_t &= \alpha_{1,t}
  \end{align*}
where both the $r \times r$ transition matrix $\mathbf{T}$ and the
covariance of the Gaussian white noise $\boldsymbol{\eta}_t$ depend
on the ARMA coefficients $\phi_i$, $\theta_j$. The stationary
covariance of the state $\boldsymbol{\alpha}_t$ can be computed by
solving a linear system. See e.g. Chap. 4 of Harvey A.C. Time Series
Models.
For the special case $p = q= 2$ you can find a closed form
for the variance if needed.
Question 0. No, $F_X^n(x)$ is not the cited penultimate
distribution, which is a a Generalized Extreme Value (GEV) with
negative shape $\xi_n < 0$ depending on $n$.  The penultimate
approximation depending on $n$ improves the convergence rate
compared to the ultimate distribution (here Gumbel). See p. 151 in
Embrechts P., Klüppelberg C. and Mikosch
T. for a
discussion. In the article by Cohen (1982) cited in OP, a
penultimate approximation is found for a sequence of i.i.d.  normal
and is shown to be such that an approximation with rate $O\{(\log
  n)^{-2}\}$ results instead of the $O\{(\log n)^{-1}\}$ rate known to
hold for the Gumbel approximation.  In Theorem 3, the case of a
Gaussian stationary times series $X_t$ is considered; It is shown
that under mild conditions on the autocorrelation sequence, the
distribution of the maximum differs from that of the maximum $n$ i.i.d. rv.s with
the same margin by $O\{(\log n)^{-2}\}$. So by triangle inequality
the penultimate approximation still leads to the better rate of
convergence when applied to the maximum of stationary Gaussian
sequences.
Question 1. I doubt that a closed-form expression would be of a
great practical interest. I think that a good approximation can be
obtained as
$$
  F_{M_n}(x)\approx  F_X^{n\theta} (x)
  $$
where $\theta \in (0,\,1)$ depends on $n$ and on the ARMA
coefficients. For a given size $n$ and given parameters we can find
a $\theta$ that leads to a good approximation for $x$ large enough,
say for $x > 0.95$. Indeed the Gaussian $\text{ARMA}(p,\,q)$ process
with given coefficients is easy to simulate from and thus it is easy
to simulate a sample of maxima $M_n$ and then find a good value for
$\theta$ by censoring the small values of $M_n$.
