[proofplan]
We pass from the original Gaussian coordinates $X_1,\dots,X_n$ to the innovation coordinates $E_1,\dots,E_n$ by the Gram-Schmidt construction in the Hilbert space $L^2(\Omega,\mathcal F,\mathbb P)$. Orthogonality of the innovations makes their covariance matrix diagonal, and Gaussianity turns this diagonal covariance into independence. The coordinate transformation from observations $x$ to innovations $e$ is lower triangular with unit diagonal, so its Jacobian determinant is $1$. The density therefore becomes the product of the one-dimensional Gaussian innovation densities evaluated at $e_t=x_t-\widehat x_t$.
[/proofplan]
[step:Construct the innovation coordinates by orthogonal projection]
For each $t \in \{1,\dots,n\}$, the space $H_{t-1}$ is finite-dimensional and hence closed in $L^2(\Omega,\mathcal F,\mathbb P)$, so the orthogonal projection $P_{t-1}$ is well-defined. Since $P_{t-1}X_t \in H_{t-1}$, there are real coefficients $a_{tj} \in \mathbb R$, indexed by $1 \le j < t$, such that
\begin{align*}
\widehat X_t = P_{t-1}X_t = \sum_{j=1}^{t-1} a_{tj}X_j.
\end{align*}
For $t=1$, this sum is empty and $\widehat X_1=0$.
Define the innovation vector
\begin{align*}
E: \Omega &\to \mathbb R^n \\
\omega &\mapsto (E_1(\omega),\dots,E_n(\omega)).
\end{align*}
Each $E_t$ is a linear combination of $X_1,\dots,X_t$, namely
\begin{align*}
E_t = X_t - \sum_{j=1}^{t-1} a_{tj}X_j.
\end{align*}
Thus $E$ is obtained from $X$ by a lower triangular linear map $B:\mathbb R^n \to \mathbb R^n$ whose diagonal entries are all $1$.
[/step]
[step:Prove that the innovations have positive variances]
Fix $t \in \{1,\dots,n\}$. Since $E_t=X_t-P_{t-1}X_t$, the defining property of orthogonal projection gives
\begin{align*}
\mathbb E[E_t Y]=0
\end{align*}
for every $Y \in H_{t-1}$.
Suppose, toward a contradiction, that $v_t=0$. Since $v_t=\mathbb E[E_t^2]$, this implies $E_t=0$ in $L^2(\Omega,\mathcal F,\mathbb P)$, hence
\begin{align*}
X_t = \sum_{j=1}^{t-1} a_{tj}X_j
\end{align*}
in $L^2(\Omega,\mathcal F,\mathbb P)$. Therefore the nonzero vector $c \in \mathbb R^n$ defined by
\begin{align*}
c_t = 1, \qquad c_j=-a_{tj}\ \text{for }1 \le j<t, \qquad c_j=0\ \text{for }j>t
\end{align*}
satisfies
\begin{align*}
c^\top X = 0
\end{align*}
in $L^2(\Omega,\mathcal F,\mathbb P)$. Taking second moments gives
\begin{align*}
0 = \mathbb E[(c^\top X)^2] = c^\top \Gamma c.
\end{align*}
This contradicts the nonsingularity of the covariance matrix of a Gaussian vector, which makes $\Gamma$ positive definite. Hence $v_t>0$.
[/step]
[step:Use Gaussianity and orthogonality to factor the innovation density]
If $s<t$, then $E_s \in H_s \subset H_{t-1}$. Since $E_t$ is orthogonal to $H_{t-1}$, we have
\begin{align*}
\mathbb E[E_tE_s]=0.
\end{align*}
Thus the covariance matrix of $E$ is diagonal:
\begin{align*}
\operatorname{Cov}(E)=\operatorname{diag}(v_1,\dots,v_n).
\end{align*}
Because $E=BX$ is a linear image of the centred Gaussian vector $X$, the random vector $E$ is centred Gaussian. A centred Gaussian vector with diagonal covariance matrix has independent components, since its characteristic function factors:
\begin{align*}
\mathbb E[e^{i u\cdot E}]
=
\exp\left(-\frac12\sum_{t=1}^n v_tu_t^2\right)
=
\prod_{t=1}^n \exp\left(-\frac12 v_tu_t^2\right)
\end{align*}
for every $u=(u_1,\dots,u_n)\in \mathbb R^n$. Therefore $E_1,\dots,E_n$ are independent centred Gaussian random variables with variances $v_1,\dots,v_n$. Their joint density with respect to $\mathcal L^n$ is
\begin{align*}
f_E(e_1,\dots,e_n)
=
\prod_{t=1}^{n}
\frac{1}{(2\pi v_t)^{1/2}}
\exp\left(-\frac{e_t^2}{2v_t}\right).
\end{align*}
[guided]
The projection construction gives orthogonality, and for Gaussian random variables orthogonality is exactly the same as independence. We verify this directly.
Take indices $s<t$. The earlier innovation $E_s$ belongs to $H_s$, because it is a linear combination of $X_1,\dots,X_s$. Since $H_s \subset H_{t-1}$ and $E_t=X_t-P_{t-1}X_t$ is orthogonal to all of $H_{t-1}$, we obtain
\begin{align*}
\mathbb E[E_tE_s]=0.
\end{align*}
Therefore distinct innovation coordinates have zero covariance, while the diagonal covariance entries are
\begin{align*}
\mathbb E[E_t^2]=v_t.
\end{align*}
Hence
\begin{align*}
\operatorname{Cov}(E)=\operatorname{diag}(v_1,\dots,v_n).
\end{align*}
Now use the Gaussian hypothesis. Since each $E_t$ is a linear combination of $X_1,\dots,X_t$, the vector $E$ is a linear image of the Gaussian vector $X$, so $E$ is again centred Gaussian. For a centred Gaussian vector with covariance matrix $\operatorname{diag}(v_1,\dots,v_n)$, the characteristic function is
\begin{align*}
\mathbb E[e^{i u\cdot E}]
=
\exp\left(-\frac12 u^\top \operatorname{diag}(v_1,\dots,v_n)u\right)
=
\exp\left(-\frac12\sum_{t=1}^n v_tu_t^2\right)
\end{align*}
for every $u=(u_1,\dots,u_n)\in \mathbb R^n$. The final expression factors as
\begin{align*}
\exp\left(-\frac12\sum_{t=1}^n v_tu_t^2\right)
=
\prod_{t=1}^n \exp\left(-\frac12 v_tu_t^2\right),
\end{align*}
which is the product of the one-dimensional characteristic functions of centred Gaussian random variables with variances $v_t$. Therefore $E_1,\dots,E_n$ are independent and have joint density
\begin{align*}
f_E(e_1,\dots,e_n)
=
\prod_{t=1}^{n}
\frac{1}{(2\pi v_t)^{1/2}}
\exp\left(-\frac{e_t^2}{2v_t}\right).
\end{align*}
[/guided]
[/step]
[step:Transform from innovation coordinates back to observed coordinates]
Define the coordinate map
\begin{align*}
T:\mathbb R^n &\to \mathbb R^n \\
x &\mapsto e
\end{align*}
by
\begin{align*}
e_t = x_t-\sum_{j=1}^{t-1}a_{tj}x_j
\end{align*}
for $t \in \{1,\dots,n\}$. This is the deterministic version of the random transformation $E=BX$. The matrix of $T$ is lower triangular with diagonal entries equal to $1$, so
\begin{align*}
|\det T|=1.
\end{align*}
The change-of-variables formula for Lebesgue density under the invertible linear map $T$ gives
\begin{align*}
f_X(x)=f_E(Tx)|\det T|.
\end{align*}
Since $|\det T|=1$ and $(Tx)_t=x_t-\widehat x_t$, we obtain
\begin{align*}
f_X(x_1,\dots,x_n)
=
\prod_{t=1}^{n}
\frac{1}{(2\pi v_t)^{1/2}}
\exp\left(-\frac{(x_t-\widehat x_t)^2}{2v_t}\right).
\end{align*}
This is the claimed likelihood factorization.
[/step]