Is a downsampled ARMA process also an ARMA process? If $X_n$ is an ARMA process, defined say by
$$
\varphi(L)X_n=\vartheta(L)\varepsilon_n,
$$
and $Y_m=X_{km}$ for some integer $k>1$, can I always represent $Y_m$ as an ARMA process?
If so, is there an explicit procedure to compute the coefficients of $Y_m$?
 A: The downsampled process $y_t=x_{kt}$ is fully specified by its autocovariance function $\dot{\gamma}_j$ which is clearly equal to the autocovariance function $\gamma_j$ of the original process $x_t$ at multiples of $k$, that is, $\dot{\gamma}_j=\gamma_{jk}$.  All properties of $y_t$ can be inferred from the downsampled autocovariance function $\dot \gamma_j$.
Autoregressive part: When $x_t$ is ARMA$(p,q)$, its autocovariance function satisfies the difference equation $\phi(B)\gamma_j = 0$ for lags higher than $q$.  In almost all cases (for exceptions see below), this difference equation has solution on the form
$$
\gamma_j = b_1 r_1^{j} + \dots + b_p r_p^{j}.
$$
where $r_i$, $i=1,2,\dots,p$ are the reciprocals of the roots of $\phi(z)$.
Hence, for lags higher than some order $\dot q$, the autocovariance function of the downsampled process $y_t$ has the form
$$
\dot\gamma_j = b_1 r_1^{kj} + \dots + b_p r_p^{kj} = b_1 (r_1^k)^j + \dots + b_p (r_p^k)^j. \tag{1}
$$
This form of autocovariance function characterises an ARMA$(p,\dot q)$ process so $y_t$ must be an ARMA process.  Except in the cases considered below, its autoregressive order is thus $\dot p=p$.
The autoregressive polynomial of the downsampled process $y_t$ is thus given by
$$
\dot\phi(z) = (1-r_1^k z)(1-r_2^k z)\cdots(1-r_p^k z).
$$
Expanding this we obtain the new autoregressive coefficients.
Roots of multiplicity higher than 1, say $2$, lead to terms such as $b_i (kj) r_i^{kj}=(b_ik)j(r_i^k)^j$ in (1) so these carry over to the downsampled model in the same way.
If the autoregressive polynomial have pairs of complex conjungate roots, then the autocovariance function exhibit exponentially decaying oscillations and it may happen that the period of those oscillations coinside with $k$.  For each such pair, the order $\dot p$ of the downsampled process will drop by 1 with the two complex conjungate roots replaced by a single real positive root.
Moving average part: The order $\dot q$ of the moving average part of the downsampled process will be at most equal to $q$ but in general lower.  For example, if $x_t$ is ARMA$(0,5)$ such that $\gamma_j=0$ (cuts off) for lags $j>5$ and $y_t=x_{2t}$ then $\dot\gamma_j=\gamma_{2j}$ clearly cuts off for lags $j>\lfloor 5/2\rfloor = 2$ such that $y_t$ is ARMA$(0,2)$.  But if $x_t$ is ARMA$(1,1)$, then $y_t$ is also ARMA$(1,1)$ as both $\gamma_j$ and $\dot\gamma_j$ will decline geometrically only for lags $j\ge 1$.  I think maybe the general formula is $\dot q=\max(\min(\dot p,q),\lfloor q/k\rfloor)$.
The MA coefficients and the new white noise variance can be found be equating the autocovariance function of the downsampled process at lag $0,1,2,\dots,\dot q$ to that of the original process $x_t$ at lag $0,k,2k,\dots,k\dot q$.  This results in a set of $\dot q + 1$ non-linear equations in the new unknown MA-coefficients and the new unknown white noise variance that in general must be solved numerically, for example, in much the same way as in this post.
New common roots: A further complication is that the MA polynomial $\dot\theta(z)$ of the downsampled process may have common roots with $\dot\phi(z)$ in which case there will be a further drop in both $\dot p$ and $\dot q$.
