Let $$ X(k) = \sum_{j = 0}^\infty a_j\, \varepsilon(k - j) , \quad k \in \mathbb{N}, $$ where $\sum_{j = 0}^\infty a_j^2 < \infty$ and $\big(\varepsilon(j) \big)_{j \in \mathbb{Z}}$ is a sequence of uncorrelated random variables with common mean zero and common finite variance $\sigma > 0$.
If $(X(1), \ldots, X(n))$ are jointly Gaussian for all $n \in \mathbb{N}$, is $\big(\varepsilon(j) \big)_{j \in \mathbb{Z}}$ a sequence of Gaussian random variables?
Related question
This question is a bit related to this math SE question but there the proof relied on the invertability of $X$. It is not clear, if this is possible here.
Hi: If you watch the video below, you will see what your $a_j$ are. From that, you can infer that you have the AR(1) below:
$X(k) = \theta \times X(k-1) + \epsilon_k$
where $X(k)$ is Gaussian.
Then, since you have an AR(1) where the response is Gaussian ( jointly Gaussian $X(k)$ implies marginally Gaussian), this implies that the $\epsilon_k$ term has to be Gaussian. Note that, When watching, disregard the constant term since you don't have one. Also, assume that the person in the video goes all the way out to $t = \infty$ rather than stopping at a finite $t$.
https://www.youtube.com/watch?v=mBYgDY7BkBc