Simulating a Stochastic Integral of OU process

Question

The stochastic integral I want to simulate is $$\int_{0}^{1}J_c(s)dJ_c(s)$$ where $J_c(s) = \int_{0}^{s}e^{-c(s-r)}dB(r)$, is an OU process.

I simulate the data using Matlab and the sample codes are as follows. n is the sample size.

dB_x = mvnrnd(zeros(K, 1), Omega_xx/n, n); % n-by-m matrix
B_x = cumsum(dB_x,1); % n-by-m matrix
B_x = [zeros(1,K); B_x];
J_c = exp(-C*[1:n]/n)'.*cumsum(exp(C*[1:n]/n).*[zeros(K,1) dB_x(1:n-1,:)'])';
dJ_c=[J_c(1,:); diff(J_c)]; % m-by-n matrix
int_J =  1/n*J_c'*dJ_c;

I am not sure whether my calculation is correct or not, can anybody give me a correct way of simulating this integral?

Answer 1

Finally, I found it wrong. The correct way is

dB_x = mvnrnd(zeros(K, 1), Omega_xx/n, n); % n-by-m matrix
B_x = cumsum(dB_x,1); % n-by-m matrix
B_x = [zeros(1,K); B_x];
J_c = exp(C*[1:n]/n)'.*cumsum(exp(-C*[1:n]/n).*[zeros(K,1) dB_x(1:n-1,:)'])';
dJ_c=[J_c(1,:); diff(J_c)]; % m-by-n matrix
int_J = J_c'*dJ_c;

Answer 2

Note that \begin{align*} dJ_c(s) = dB_s -cJ_c(s) ds. \end{align*} Moreover, $d\langle J_c, J_c\rangle_t = dt$. Then, by Ito's lemma, \begin{align*} \int_0^1J_c(s) dJ_c(s) &= \int_0^1 d\left(J_c^2(s) -s \right)\\ &= J_c^2(1) - 1. \end{align*} Since \begin{align*} J_c(1) &= \int_0^1 e^{-c(1-s)}dB_s \end{align*} is normal with zero mean and the variance given by \begin{align*} E\left(J_c^2(1) \right) &= \int_0^1e^{-2c(1-s)} ds\\ &=\frac{1}{2c}(1-e^{-2c}). \end{align*} Then \begin{align*} J_c(1) = \sqrt{\frac{1}{2c}(1-e^{-2c})}\, Z, \end{align*} and \begin{align*} \int_0^1J_c(s) dJ_c(s) &= \frac{1}{2c}(1-e^{-2c})\, Z^2 -1. \end{align*} where $Z$ is a standard normal random variable. The simulation becomes a draw of a sample from a standard normal random variable.