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 $$
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_g,\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$.

$\begingroup$ I understand it has been a year since the answer, but I am not convinced this answers the original question. The calculated variance and covariance are not the conditional variances that are requested, but they are averaged over the realizations of both $I_t(g_1)$ and $I_t(g_2)$. However, what is asked here is to find the conditional variance (and mean) of $I_t(g_2)$ given one specific value of $I_t(g_1)$. I have asked a similar question as I haven't found this one: mathoverflow.net/questions/315781/… $\endgroup$ – Dr_Zaszuś Nov 20 '18 at 16:17

1$\begingroup$ @Dr_Zaszuś : If the joint distribution is Gaussian (as shown in my answer) and known, then all the conditional distributions are known (and Gaussian)  this is a textbook fact. $\endgroup$ – Iosif Pinelis Nov 20 '18 at 17:13

$\begingroup$ Yes, right. I have realized that later. What I mean is that your answer proves that the joint distribution is normal. It however formally didn't answer the OP question about how to sample the two functions together or to calculate the conditional distribution. The answer can be deduced from your answer, but is not there... It took me a while to make the connection. You might want to update your answer and spare further readers the thinking. $\endgroup$ – Dr_Zaszuś Nov 21 '18 at 15:39

1$\begingroup$ @Dr_Zaszuś : I have added the explicit expression for the conditional distribution. $\endgroup$ – Iosif Pinelis Nov 22 '18 at 21:26