(Non-limit) distribution of maxima from different univariate, discrete and stationary time series

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^{ntheta} (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$.