I would like to simulate three Gaussian processes $X$, $Y$ and $Z$ defined as
$$dX_t=-(a_xX_t+b_x)dt+\sigma_XdW_t^x$$
$$dY_t=-(a_yY_t+b_y)dt+\sigma_YdW_t^y$$
$$dZ_t=(\int_{0}^tX_udu-\int_{0}^tY_udu-\frac{1}{2}\sigma_Z^2)dt+\sigma_ZdW_t^z$$
with $(a_x,b_x,a_y,b_y)\in\mathbb{R_{+}^4}$, $(\sigma_X,\sigma_Y,\sigma_Z) \in \mathbb{R_{+}^3}$ and $$d<W_t^x,W_t^y>=\rho_{xy}dt$$
$$d<W_t^x,W_t^z>=\rho_{xz}dt$$
$$d<W_t^y,W_t^z>=\rho_{yz}dt$$
$(\rho_{xy},\rho_{xz},\rho_{yz}) \in [-1,1]^3$ (chosen in a way to generate a correlation matrix)
Ultimately , I want to estimate several expectations $$E(e^{a_1^l(t)X_t+a_2^l(t)Y_t+a_3^l(t)Z_t})$$ at any time $t$, with $(a_1^l(t),a_2^l(t),a_3^l(t)) \in\mathbb{R_{+}^3}$,$l >0$ and make sure that the $g$ functions defined as follows are smooth functions $$g_{l}(\sigma_x,\sigma_y,\sigma_z)=E(e^{a_1^l(t)X_t+a_2^l(t)Y_t+a_3^l(t)Z_t})$$
For the sake of simplicity, we can choose $a_x=a_y=0.01$, $\sigma_X=\sigma_Y=0.01$, $\sigma_Z=0.1$ $b_x=0$ , $b_y=0.0002$, $(\rho_{xy}=0.5,\rho_{xz}=-0.15,\rho_{yz}=0.2)$
By construction, $U_t=(X_t,Y_t,Z_t)$ is a Gaussian vector, hence we know its exact distribution at any time $t$. One can work out its mean $m_t$ and covariance matrix $C_t=L_tL_t^T$ (fastidious). Moreover, we can define its conditional( to $U_s$) mean $m_{t,s}$ and conditional covariance matrix $C_{t,s}=L_{t,s}L_{t,s}^T$, with $t>s$.
The error that one would expect when simulating the processes must come from the statistical error due to the Monte-Carlo scheme.
I build a time grid $\{t_i\}_{k=0}^{N}$ such as $0=t_o<t_1<...<t_N$, and we want a sample of size $N_{sample}$.
At time $t_i$, the $k^{th}$ scenario out of $N_{sample}$ is generated using the formula
$$U_{t_i}^{k}=m_{t_i,t_{i-1}}^{k}+L_{t_i,t_{i-1}}R_{i,k}$$
where $U_{t_i}^{k}$ is the $k^{th}$ sample of $U$ at time $t_i$, $m_{t_i,t_{i-1}}^{k}$ the mean of $U_{t_i}$ conditional to $U_{t_{i-1}}$ on the $k^{th}$ scenario. $R_{i,k}$ is standard normal random vector.
To achieve my goal, I made sure that the first and second moment of $U_t$ are replicated for the different time $t$. I used the generated $R_{i,k}$ straightforwardly , I realized that the mean of $U_t$ was a bit off...
I then assumed that $$R_{i,k}=-R_{i-1,k}$$ The results was much better for $X$ and $Y$, but I noticed a bias that flaws the dynamics of $Z_t$.
To illustrate the problem , we can take two correlated Brownian motions $W_1,W_2$,
We define two processes :
$X_1(t)=X_1(s)+W_1(t)-W_1(s) , X_1(0)=m $ $X_2(t)=X_2(s)+X_1(s)f(s,t)+W_2(t)-W_2(s) , s<t$
$f$ is deterministic, or whatever suits me, including $f(t,t)=0$ . We consider two periods $[0,t_1]$ and $[t_1,t_2]$, and two scenarios ( one scenario $X_1^N(t)$ and the antithetic one $X_1^A(t)$).
On $[0,t_1]$, we have $\frac{X_1^N(t_1)+X_1^A(t_1)}{2}=m=E(X_1(t_1))$ and $\frac{X_2^N(t_1)+X_2^A(t_1)}{2}=X_2(0)+mf(0,t_1)=E(X_2(t_1))$.
We move to the next period:
On $[t_1,t_2]$, we have $\frac{X_1^N(t_2)+X_1^A(t_2)}{2}=m=E(X_1(t_2))$ and $\frac{X_2^N(t_2)+X_2^A(t_2)}{2}=X_2(0)+m(f(0,t_1)+f(t_1,t_2))$ but $E(X_2(t_2))=mf(0,t_2)$ and $f(0,t_1)+f(t_1,t_2)$ is not equal to$f(0,t_2)$ (I constructed it that way). We have the same issue in my original problem.
I know that having a bias on $U_t$ is not an issue as long as we estimate $E(e^{a_1(t)X_t+a_2(t)Y_t+a_3(t)Z_t})$ correctly and there is no noise when stressing the "sigmas": How to build the right bias?
I am aware of the possibilities to readjust $U_t$ (importance sampling) in a way that I match the target expectation, but as defined above, I want to replicate several expectations.
Thanks.