First see what definitions say
$\{X_t\}$ is strictly stationary if for any $ t_1,t_2,...,t_n \in T$ and any $k \in T$
$$ P(X_{t_1},...,X_{t_n}) = P(X_{t_1+k},...,X_{t_n+k})$$
that is, we have statistical equilibrium.
$\{X_t\}$ is second order/weakly stationary if
- $\mathbb{E}[X_t] = \mu < \infty$, independent of $t$
- $Var(X_t) = \sigma^2 <\infty$ , independent of $t$
- $cov(X_t,X_{t+\tau})$, is a function of $\tau$ only.
Think of what weakly stationarity means. It means that the expected value of the process is finite and constant. It also means that the autocovariance does not depend on where two random variables are positioned but just on their distance! Autocovariance between today and yesterday same as autocovariance between 100 days and 101 days ago.
Now take the MA(1) process
$$\boxed{ X_t = \varepsilon_t + \theta \varepsilon_{t-1} }$$
where $\varepsilon_t \sim WN(0,\sigma^2)$.
This process is weakly stationary as
- $\mathbb{E}[X_t] = \mathbb{E}[\varepsilon_{t}] + \theta \mathbb{E} [\varepsilon_{t-1}] = 0$
- $Var(X_t) = \sigma^2(1+\theta^2)$, i.e. independent of $t$
- $cov(X_t,X_{t-\tau}) = \theta \sigma^2, |\tau|=1$, $cov(X_t,X_{t-\tau}) = 0, \forall |\tau|>1$ hence function of $\tau$ only.
(if you have difficulties deriviting that, let me know. It is very straightforward and you shouldn't have any difficulties)
MA(1) is also strictly stationary as both $P(X_{t_1},...,X_{t_n})$ and $P(X_{t_1+k},...,X_{t_n+k})$ multivariate (1-dependent) Normal distributions with identical parameters as it is a combination of WN random variables.
Edit to provide further explanation:
Each $X_{t_1},..X_{t_n}$ is a linear combination of independent Gaussian random variables,. In particular all $X_{t_i} \sim N(0, \sum_{j=1}^{q} \theta_j^2 \sigma^2) \forall j $. The joint distribution of normal random variables is a multivariate normal. As they are all driven by an MA(q) process their covariance matrix can be easily derived. But as seen above, this does not depend on location, just on distance. Explicitly for MA(1):
$ f_{\mathbf {X} }(X_{t_1} \dots X_{t_n})={\frac {\exp \left(-{\frac {1}{2}}{\mathbf {X} }^{\mathrm {T} }{\boldsymbol {\Sigma }}^{-1}{\mathbf {X} }\right)}{\sqrt {(2\pi )^{k}|{\boldsymbol {\Sigma }}|}}}$
where
$$\boldsymbol {\Sigma } = \begin{bmatrix}
\sigma^2(1+\theta^2) & \theta \sigma^2 & 0 & 0 & \dots & 0 \\
\theta \sigma^2 & \sigma^2(1+\theta^2) & \theta \sigma^2 & 0 & \dots & 0 \\
0 & \theta \sigma^2 & \sigma^2(1+\theta^2) & \theta \sigma^2 & \ddots & \vdots \\
0 & 0 & \ddots & \ddots & \ddots & 0 \\
\vdots & \vdots & \ddots & \ddots & \ddots & \theta \sigma^2 \\
0 & 0 & \dots & 0 & \theta \sigma^2 & \sigma^2(1+\theta^2)
\end{bmatrix}$$
Hence MA(q) processes driven by Gaussian noise are always strictly stationary (only condition is that the sum of MA coefficients should be finite)
In general, all weakly stationary Gaussian processes are strictly stationary too.