Lo intentaré, pero aún no estoy 100% seguro de que sea el camino a seguir.
Ansatz:
Busquemos la distribución de la integral de un movimiento browniano con respecto al tiempo (llámalo $x$ ) y encontrar la expectativa del producto de dos integrales de este tipo $x$ y $y$ . A continuación, calcula la covarianza como $Cov(x,y)=E(xy)-E(x)E(y)$ .
1. Distribución de $I(t,T)\equiv\int\limits_{s=t}^{T}W_s\mathrm{d}s$
Desde esta respuesta sabemos que $$ I(t,T)\equiv \int\limits_{s=t}^{T}W_s\mathrm{d}s=\int\limits_{s=t}^{T}(T-s)dW_s $$
y que $I(t,T)$ se distribuye normalmente con $$I(t,T)\sim \mathrm{N}\left(0,\frac{1}{3}(T-t)^3\right)$$
2. Expectativa de $I(t,T)I(t,U)$
Por la misma vía de razonamiento (y algún débil recuerdo de Iso isometría yo diría que sí:
$$ \begin{align} \mathrm{E}\left(I(t,T)I(t,U)\right)&=\mathrm{E}\left(\int\limits_{s=t}^{T}(T-s)\mathrm{d}W_s\int\limits_{x=t}^{U}(U-x)\mathrm{d}W_x\right)\\ &=\mathrm{E}\left(\int\limits_{s=t}^{T}\int\limits_{x=t}^{U}(T-s)(U-x)\mathrm{d}W_s\mathrm{d}W_x\right)\\ &=\int\limits_{s=t}^{T}\int\limits_{x=t}^{U}(T-s)(U-x)\mathrm{E}\left(\mathrm{d}W_s\mathrm{d}W_x\right)\\ &=\int\limits_{s=t}^{U}(U-s)^2\rho\mathrm{d}t\\ &=\frac{1}{3}(U-t)^3 \end{align} $$ N.B.: suponemos que $U\leq T$ .
3. $I(t,T)$ y $I(t,U)$ se distribuyen normalmente de forma bivariada
Simplifiquemos y dejemos que $x_1\equiv I(t,T)$ , $x_2\equiv I(t,U)$ y $\mathbf{x}=\left(x_1,x_2\right)^T$ También deja que $\sigma_1^2=\frac{1}{3}(T-t)^3$ , $\sigma_2^2=\frac{1}{3}(U-t)^3$ y $\sigma_{1,2}=\frac{1}{3}\rho (U-t)^3$ . Entonces $\mathbf{x}$ se distribuye normalmente de forma bivariada como
$$ \mathbf{x}\equiv\begin{pmatrix}x_1\\x_2\end{pmatrix}\sim \mathrm{N}\left(\mathbf{0},\begin{pmatrix}\sigma_1^2 & \sigma_{1,2}\\ \sigma_{1,2} & \sigma_2^2\end{pmatrix}\right) $$
Ahora vamos a utilizar un pequeño truco que aprendió hace poco . Dado un vector real $\mathbf{t}$ la MGF de la distribución normal multivariante se define como $\varphi_X(t)\equiv\mathrm{E}\left(e^{\mathbf{t}^T\mathbf{x}}\right)=e^{\mathbf{t}^T\mathbf{\mu}+\frac{1}{2}\mathbf{t}^T\mathbf{\Sigma}\mathbf{t}}$ y, en nuestro caso, esto es
$$ \varphi(t_1,t_2)=\mathrm{E}\left(e^{t_1x_1+t_2x_2}\right)=e^{\frac{1}{2}t_1^2\sigma_1^2+\frac{1}{2}t_2^2\sigma_2^2+t_1t_2\sigma_{1,2}} $$
Tenga en cuenta que
$$ \left.\frac{\partial \left(e^{x+ty}\right)}{\partial t}\right|_{t=0}=ye^x $$
así,
$$ \begin{align} \mathrm{E}\left(e^{I(t,T)}I(t,U)\right)&=\mathrm{E}\left(e^{x_1}x_2\right)\\ &=\left.\frac{\partial \varphi(t_1=1,t_2)}{\partial t_2}\right|_{t_2=0}\\ &=e^{\frac{1}{2}\sigma_1^2}\sigma_{1,2}\\ &=\frac{1}{3}\rho(U-t)^3e^{\frac{1}{6}(T-t)^3} \end{align}$$
5. Poner todo junto
Así, para los movimientos brownianos $W^1_t, W_2^t$ con $dW_1dW_2=\rho dt$
$$ \begin{align} \mathrm{Cov}\left(e^{\int\limits_{s=t}^TW^1_s\mathrm{d}s}\int\limits_{x=t}^UW^2_x\mathrm{d}x\right)&=\mathrm{E}\left(e^{x_1}x_2\right)-\mathrm{E}\left(e^{x_1}\right)\mathrm{E}(x_2)\\ &=\mathrm{E}\left(e^{x_1}x_2\right)\\ &=1/3\rho(U-t)^3e^{\frac{1}{6}(T-t)^3} \end{align} $$