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;