Unconditional distribution of ARMA process with t-student errors 
In the $Y_t\sim ARMA(p,q)$ model, when the errors have Normal distribution, the unconditional distribution of $Y_t$ is Normal. When the errors have a
  t-student distribution with $\nu$ degrees of freedom. What is the
  unconditional distribution of $Y_t$?

So
$$Y_t=\phi_1Y_{t-1}+\dots+\phi_pY_{t-p}+e_t-\theta_1e_{t-1}-\dots-\theta_q e_{t-q}$$
where $e_t\sim t_\nu$.
I have no idea how to find the distribution of it and the books that I have mostly cover only the case with Gaussian errors.
Some reference would also be interesting.
 A: The ARMA model is within the general class of linear models where your observable vector is a linear function of an underlying vector of IID error terms.  Consider the general linear model form with IID errors following a T-distribution:
$$Y_t = \sum_{k=0}^\infty A_k \varepsilon_{t-k} \quad \quad \quad \varepsilon_k \sim \text{IID Student T}(\text{df} = \nu). \quad \quad \quad \quad$$
One of the useful properties of the Student's T-distribution is that it can be written as a mixture of normal with a gamma-distributed precision parameter.  With this representation, the above model form can be written equivalently as:
$$Y_t = \sum_{k=0}^\infty \frac{A_k \epsilon_{t-k}}{\sqrt{\lambda_k}} \quad \quad \quad \epsilon_k \sim \text{IID N}(0, 1) \quad \lambda_k \sim \text{Gamma}(\tfrac{\nu}{2}, \tfrac{\nu}{2}).$$
You can see from this form that the value $Y_t$ is a sum of independent terms which are each ratios of normal random variables and gamma random variables (giving scaled T random variables).  The difference between the present model and the standard Gaussian linear model is the presence of the denominator terms in the sum.  (In the standard Gaussian case we have fixed $\lambda = \lambda_k$.)
The distribution for this quantity is a complicated convolution, but the CLT ensures that it converges to normality under mild conditions.  It is possible to simulate the distribution by applying the random root-gamma denominators to the summation terms, which gives a bit more variability to the quantity than occurs in the standard Gaussian linear model.
