It seems that the discussions on ARMA are always focused on weak (second-order) stationary, but what about strong stationary? What are the conditions on the coefficients for it to be strictly stationary and what is the corresponding unconditional distribution?
Furthermore, is there a general procedure for deriving its unconditional distribution? I have tried deriving it in terms of the characteristic function, but not sure if it is correct. For example, for $AR(1)$ given by
$$x_t = ax_{t-1} + z_t$$
we get
$$\psi_X(s) := \mathbb{E}[\exp(i s x_t)] = \mathbb{E}[\exp(i a s x_{t-1})]\mathbb{E}[\exp(i s z_t)] = \psi_X(as)\exp(-0.5\sigma^2s^2)$$
assuming $z_t \sim N(0,\sigma^2)$, but it is not clear to me how to proceed from here. One can guess a solution and check it by substitution, so, using the knowledge of the variance of a weak stationary $AR(1)$, we may guess that
$$\psi_X(s) = \exp\left(-\frac{1}{2}\frac{\sigma^2}{1-a^2}s^2\right)$$
which gives
$$\psi_X(as) = \exp\left(-\frac{1}{2}\frac{\sigma^2 a^2}{1-a^2}s^2\right)$$
and we find that
$$\frac{\psi_X(s)}{\psi_X(as)} = \exp\left(-\frac{1}{2}\frac{\sigma^2}{1-a^2}s^2\right)\exp\left(\frac{1}{2}\frac{\sigma^2 a^2}{1-a^2}s^2\right) = \exp\left(-\frac{1}{2}\sigma^2s^2\right)$$
which satisfies the equation. So, it would seem that weak stationary is a sufficient condition for $AR(1)$ to be also strongly stationary (probably this only holds for this special case of Gaussian noise, and I guess there already must be a theorem that proves this). But how would one proceed in a general case, where noise is not restricted to Gaussian and without guessing? And what about more complex stationary processes, say, like $GARCH$? For example, $ARCH(1)$ is
$$x_t = \sigma_t z_t$$
$$\sigma_t^2 = a_0 + a_1 x_{t-1}^2 $$
which we can re-write as
$$x_t^2 = (a_0 + a_1 x_{t-1}^2) z_t^2 .$$
which then gives
$$\psi(s) := \mathbb{E}[\exp(i s x_t^2)] = \mathbb{E}[\exp(i s a_0 z_t^2)] \mathbb{E}[\exp(i s a_1 x_{t-1}^2 z_t^2)] .$$
If we again assume that $z_t \sim N(0,1)$, then $z_t^2 \sim \chi^2(1)$ and $c z_t^2 \sim \Gamma(0.5,2c)$, then we have
$$\mathbb{E}[\exp(i s a_0 z_t^2)] = (1 – 2 a_0 i s)^{-0.5}$$
and, conditional on $x_{t-1}$,
$$\mathbb{E}_{t-1}[\exp(i s a_1 x_{t-1}^2 z_t^2) | x_{t-1}] = (1 – 2 a_1 x_{t-1}^2 i s)^{-0.5}$$
so that
$$\psi(s) = (1 – 2 a_0 i s)^{-0.5} \mathbb{E}[(1 – 2 a_1 x_{t-1}^2 i s)^{-0.5}]$$
and guessing the solution is not so easy now.
Add 1
It might be possible to use fractional integral (using Cauchy formula) to write the above as
$$\psi(s) = (1 – 2 a_0 i s)^{-0.5} (2 a_1 i s)^{-0.5}
\mathbb{E}\left[\left(\frac{1}{2 a_1 i s} – x^2\right)^{-0.5}\right]
= (1 – 2 a_0 i s)^{-0.5} (2 a_1 i s)^{-0.5}\Gamma(0.5)(J^{0.5}f)\left(\frac{1}{2 a_1 i s}\right)$$
where $f$ is the pdf of $x^2$. So we get a fractional integral equation (although I am not sure about the accuracy of the above).
Best Answer
You are making the problem much harder than needed. A necessary condition that the process ${X_t}$ be stationary is that $X_t$ must have the same distribution as $X_{t-1}$. In particular, $X_t$ must have the same mean and variance as $X_{t-1}$. Now, in order for $E[X_t] = E[X_{t-1}]$ to hold, it must be that $$E[X_t] = aE[X_{t-1}] + E[Z_t] = aE[X_{t-1}] \implies E[X_t] = 0 ~\forall t.$$ Similarly, in order for $\operatorname{var}(X_t) = \operatorname{var}(X_{t-1})$ to hold, it must be that $$\operatorname{var}(X_t) = a^2 \operatorname{var}(X_{t-1}) + \operatorname{var}(Z_t) = a^2 \operatorname{var}(X_t) + \sigma^2 \implies \operatorname{var}(X_t) = \frac{\sigma^2}{1-a^2}.$$ Now, the distribution of $X_t$ is the convolution of the distributions of $aX_{t-1}$ and $Z_t$, and if $X_{t-1}$ and $Z_t$ are (independent) normal random variables, so is $X_t$ a normal random variable.
With regard to strict stationarity when the $X$'s and $Z$'s are normal random variables, see the last two paragraphs of this answer of mine over on dsp.SE.