Conditional stochastic integration Let's say we have two functions $h(s)$ and $g(s)$. We can easily simulate a stochastic integral, e.g.
$$t \mapsto \int_0^t h(s) dB(s) \sim \mathcal{N}\bigg(0, \int_0^t h(s)^2 ds \bigg). $$
What is the conditional distribution of stochastic integral of $g(s)$ with respect to $B(s)$ then?
$$ t \mapsto \int_0^t g(s) dB(s) \bigg | \int_0^t h(s) dB(s) \sim ?$$
I have an intuition that these two integrals might have a joint normal distribution with a covariance equal to $\int_0^t h(s)g(s)ds$, but I need a rigorous derivation. Generally, I want to find a way to simulate a bunch of $n$ integrals with respect to the same Brownian motion:
$$ t \mapsto \bigg( \int_0^t h_1(s) dB(s),  \int_0^t h_2(s) dB(s) , \dots ,  \int_0^t h_n(s) dB(s) \bigg)^T  $$
 A: (From the context of your post, it is apparent that $B$ is assumed to be a (say standard) Brownian motion.)
The joint normal distribution follows from the way the stochastic integral is defined. However, if you are already convinced that
\begin{equation}
    I_t(h):=\int_0^t h(s) dB(s) \sim N\bigg(0, \int_0^t h(s)^2 ds \bigg)
\end{equation}
for all $h\in L^2([0,t])$, then you have
\begin{equation}
    \sum_1^n r_i I_t(g_i)=I_t\Big(\sum_1^n r_i g_i\Big) \sim N\bigg(0, \int_0^t \Big(\sum_1^n r_i g_i(s)\Big)^2 ds \bigg) \tag{*}
\end{equation}
for all real $r_i$'s and all functions $g_i\in L^2([0,t])$. So, all linear combinations of the $I_t(g_i)$'s are Gaussian and hence the $I_t(g_i)$'s are jointly Gaussian (think, e.g., of the joint characteristic function). Formula (*) for $n=1,2$ also yields the covariances:
\begin{align*}
    Var(I_t(g_1)+I_t(g_2))&=\int_0^t (g_1(s)+g_2(s))^2 ds \\ 
    &=\int_0^t  g_1(s)^2 ds+\int_0^t  g_2(s)^2 ds
    +2\int_0^t g_1(s)g_2(s)\,ds \\ 
    &=Var(I_t(g_1))+Var(I_t(g_2))+2\int_0^t g_1(s)g_2(s)\,ds,   
\end{align*}
so that
$Cov(I_t(g_1),I_t(g_2))=\int_0^t g_1(s)g_2(s)\,ds$.
Added in response to a comment by Dr_Zaszuś:
Thus, for any $g,h$ in $L^2([0,t])$, the pair $(I_t(g),I_t(h))$ has the bivariate normal distribution $N(\mu_g,\mu_h,\sigma^2_g,\sigma^2_h,\rho)$, where $\mu_g=\mu_h=0$, $\sigma_g=\|g\|$, $\sigma_h=\|h\|$, and $\rho=\frac{g\cdot h}{\|g\|\|h\|}$, where $g\cdot h:=\int_0^t g(s)h(s)\,ds$ and $\|f\|:=\sqrt{f\cdot f}$. So, the conditional distribution of $I_t(g)$ given $I_t(h)=y$ is the normal distribution $N(\mu_{g;y},(1-\rho^2)\sigma_g^2)$ for any real $y$, where $\mu_{g;y}:=\rho\sigma_g y/\sigma_h$.
