The finite linear-Gaussian state-space path is jointly Gaussian

Each Gaussian marginal law in the model exists on its finite-dimensional Euclidean Borel space. Repeatedly applying the finite-product measure construction creates a finite product probability space with these marginals. Its coordinate maps $x_0,q_1,r_1,\dots ,q_T,r_T$ have the prescribed marginals, and the rectangle formula gives product probabilities on all measurable coordinate rectangles. Therefore the generated sigma-algebras are mutually independent, and the finite family has the product joint law. Thus this space supplies the required independent Gaussian noise vectors.

For a Borel set $A\subseteq {\mathbb{R}}^{d_t}$, define $K_t(x,A)=\int 1_A(F_{t}x+z)N_{d_t}(0,Q_t)(dz)$; define $G_t$ analogously using $H_t$ and $N_{p_t}(0,R_t)$. The maps $(x,z)\mapsto F_{t}x+z$ and $(x,z)\mapsto H_{t}x+z$ are continuous, so the integrands are product-measurable. Each fixed-state section is a pushforward probability measure. Regard each fixed Gaussian noise law as a constant probability kernel. Kernel integration then makes the evaluations measurable in the state argument, so both families are probability kernels. For fixed $x$ and any $u$, the scalar affine characteristic-function identity applied to $u^T{q}_t$, together with the multivariate-normal characteristic-function formula, gives $\mathbb{E}[e^{i{u}^T(F_{t}x+q_t)}]=\exp (i{u}^T{F}_{t}x-\frac{1}{2}{u}^T{Q}_{t}u)$. This is the multivariate-normal characteristic function with mean $F_{t}x$ and covariance $Q_t$, so its uniqueness identifies the translated law with $N_{d_t}(F_{t}x,Q_t)$. The same calculation with $H_t,r_t,R_t$ identifies $G_t$. No covariance inverse is used.

Let $\varphi$ be any nonnegative Borel function on the finite trajectory space. The recursion is a finite composition of continuous affine maps, so $\varphi (Z_T)$ is measurable. By Tonelli, integrate the noise product law from step 1.1 backward, starting with the last pair $(q_T,r_T)$. For each fixed preceding history, the integral over $q_T$ is the pushforward through $x_T=F_T{x}_{T-1}+q_T$, hence equals integration against $K_T$; the following integral over $r_T$ is the pushforward through $y_T=H_T{x}_T+r_T$, hence equals integration against $G_T$. Repeating this argument for $t=T-1,\dots ,1$ and finally integrating $x_0$ against ${\mu }_0$ gives $\mathbb{E}[\varphi (Z_T)]=\int {\mu }_0(d{x}_0){\prod }_{t=1}^T{K}_t(x_{t-1},d{x}_t)G_t(x_t,d{y}_t)\varphi (x_0,x_1,y_1,\dots ,x_T,y_T)$. For $T=0$, it is simply integration against ${\mu }_0$. This is exactly the unique general model law, and the trajectory space is standard Borel.

The stacked-noise lemma gives the vector $W_T$ the multivariate normal law $N(a_T,{\mathrm{\Sigma }}_T)$. Successive substitution in the two recursions shows by finite induction that each coordinate of $Z_T$ is a linear function of $W_T$; collecting those coefficients gives a deterministic matrix $M_T$. For any vector $u$, $u^T{Z}_T=(M_T^{T}u{)}^T{W}_T$, whose law is normal with mean $u^T{M}_T{a}_T$ and variance $u^T{M}_T{\mathrm{\Sigma }}_T{M}_T^{T}u$. The multivariate-normal projection definition gives $Z_T\sim N(M_T{a}_T,M_T{\mathrm{\Sigma }}_T{M}_T^T)$. The covariance is positive semidefinite because $v^T{\mathrm{\Sigma }}_{T}v\ge 0$ implies $u^T{M}_T{\mathrm{\Sigma }}_T{M}_T^{T}u=(M_T^{T}u{)}^T{\mathrm{\Sigma }}_T(M_T^{T}u)\ge 0$. Thus singular covariances are included and no inverse is needed.