Proof of Mean-Square Riemann Integrals of a Jointly Gaussian Family are Jointly Gaussian

Fix versions $I_t$ of the mean-square Riemann integrals $\int_a^t H_u\,du$ for $t\in[a,b]$, with $I_a=0$. By Jointly Gaussian Families of Random Variables and Gaussian Processes it suffices to show: for every choice of finitely many distinct indices of the enlarged family, the tuple of corresponding members is a Gaussian random vector. Since the index set is the disjoint union of $J$ and two copies of $[a,b]$, such a tuple has the form $(X_{j_1},\dots,X_{j_p},H_{t_1},\dots,H_{t_q},I_{s_1},\dots,I_{s_r})$ with distinct $j_l\in J$, distinct $t_w\in[a,b]$, and distinct $s_v\in[a,b]$, where any of $p,q,r$ may be $0$ but not all three. (The same random variable may occur under different indices; this is permitted, as Gaussian random vectors may have repeated coordinates.)

For each natural number $n\ge1$ and each $v\in\{1,\dots,r\}$ define a random variable $S^{(v)}_n$ as follows: if $s_v=a$, set $S^{(v)}_n=0$; otherwise let $S^{(v)}_n$ be the mean-square Riemann sum of $(H_u)_{u\in[a,s_v]}$ over the partition of $[a,s_v]$ into $n$ equal intervals with left endpoints as tags. Since $I_{s_v}$ is a mean-square Riemann integral over $[a,s_v]$ and the mesh $(s_v-a)/n$ tends to $0$, Mean-Square Riemann Integral of a Family of Random Variables gives $\lVert S^{(v)}_n-I_{s_v}\rVert_2\to0$ as $n\to\infty$; when $s_v=a$ this holds trivially because $I_a=0=S^{(v)}_n$.

Fix $n$. The random variables appearing in the tuple $(X_{j_1},\dots,X_{j_p},H_{t_1},\dots,H_{t_q},S^{(1)}_n,\dots,S^{(r)}_n)$ are built from finitely many members of the combined family: the $X_{j_l}$, the $H_{t_w}$, and the values $H_u$ at the finitely many tags of the $r$ partitions. Choose finitely many distinct indices of the combined family whose members include all of these; by the joint Gaussianity hypothesis and Jointly Gaussian Families of Random Variables and Gaussian Processes, the corresponding tuple is a Gaussian random vector $V$. Every coordinate of the tuple above is a finite linear combination of coordinates of $V$ (a single coordinate for the $X_{j_l}$ and $H_{t_w}$; the combination $\sum_i H_{\tau_i}(x_i-x_{i-1})$ for $S^{(v)}_n$ with $s_v>a$; the empty combination, that is $0$, for $S^{(v)}_n$ with $s_v=a$). Hence the tuple is the image of $V$ under a linear map and is a Gaussian random vector by Affine Transformations of Gaussian Random Vectors are Gaussian.

Letting $n\to\infty$, the coordinates of these Gaussian random vectors converge in mean square to the coordinates of $(X_{j_1},\dots,X_{j_p},H_{t_1},\dots,H_{t_q},I_{s_1},\dots,I_{s_r})$: the first $p+q$ coordinates are constant sequences, and $\lVert S^{(v)}_n-I_{s_v}\rVert_2\to0$ for each $v$. The limit coordinates are square-integrable, as required by Mean-Square Limits of Gaussian Random Vectors are Gaussian: the $X_{j_l}$ and $H_{t_w}$ are Gaussian, hence square-integrable by Square-Integrability, Moments, and Covariance Matrix of a Gaussian Random Vector, and the $I_{s_v}$ are square-integrable by Mean-Square Riemann Integral of a Family of Random Variables. By Mean-Square Limits of Gaussian Random Vectors are Gaussian, the limit tuple is a Gaussian random vector. Since the choice of finitely many distinct indices was arbitrary, the enlarged family is jointly Gaussian. $\blacksquare$