For the case of input process $\{X(t)\}$ being white Gaussian noise with two-sided power spectral density $\frac{N_0}{2}$, the output process $\{Y(t)\}$ is a strictly stationary zero-mean Gaussian process in which all the random variables have the same variance $\frac{N_0}{2}\int_{-\infty}^\infty |H(f)|^2 \,\mathrm df$ almost as you say. But the key point here is that the $Y$'s are zero-mean Gaussian random variables with known variance, and it is known that if $Y \sim N(0,\sigma_Y^2)$ then $$E[Y^n] = \begin{cases}\sigma_Y^n \cdot (n-1) \times (n-3) \times (n-5)\times \cdots \times 3 \times 1, & n ~ \text{even},\\ 0, & n ~ \text{odd}.\end{cases}$$
Now, if the input process is just white noise (not necessarily Gaussian white noise), then the output process is not necessarily Gaussian, and while the mean and variance are as specified above, we cannot say anything about the distribution of the $Y(t)$'s and cannot infer anything about the higher moments of the $Y(t)$'s the way we could for Gaussian output processes.
If the input process is a zero-mean WSS Gaussian process (which is also a strictly stationary process) but not necessarily a _white noise process, than $\{Y(t)\}$ is also a zero-mean WSS Gaussian process (and so strictly stationary too) with variance $$\sigma_Y^2 = \int_{-\infty}^\infty S_X(f)|H)f)|^2 \,\mathrm df$$ where $S_X(f)$ is the power spectral density of the input process and so once again we are in business and can use $$E[Y^n] = \begin{cases}\sigma_Y^n \cdot (n-1) \times (n-3) \times (n-5)\times \cdots \times 3 \times 1, & n ~ \text{even},\\ 0, & n ~ \text{odd}.\end{cases}$$
