Have trouble understanding Matlab code of stochastic integrals

Question

Hell, I have trouble understanding the matlab code of approximating stochastic integrals.

rng(100)
T = 1; N = 500; dt = T/N;

dW = sqrt(dt) * randn(1,N);
W = cumsum(dW);

ito = sum([0, W(1:end-1)].*dW)
strat = sum((0.5*([0, W(1:end-1)] + W) + 0.5*sqrt(dt)*randn(1,N)).*dW)

The ito refers Ito integrals, whose formula is $\sum_{j=-0}^{N-1} W(t_j)(W(t_{j+1})-W(t_j))$

But I just cannot understanding it from the coding perspective. I think W refers to $W(t_{j+1})+W(t_j)$ which is not part of the formula.

Similarly, strat refers to Stratonovich integrals, whose formula is $ \sum_{j=-0}^{N-1} W(t_j)(\frac{t_j + t_{j+1}}{2})(W(t_j) - W(t_j))$. But I cannot find the $W(t_j) - W(t_j)$ part in the coding.

Please help me by point out which part of the coding refers to which part of the formula.

Thank you so much

Answer 1

dW are increments, which are iid normals with mean $0$ and variance $\Delta t$. In math notation we write the increments as $W(t_{j+1})-W(t_j)$ but numerically the information flows the other way: you generate increments and then reconstruct the values of W from the increments.

You get W at a particular time by summing up the increments up to that time, since the process starts up at zero. This is what "cumsum" is doing.

Note that in the Ito form of the integral, you multiply the next increment by the present value of $W$. The vector W in the code only contains $W(t_j)$ for $j \geq 1$, so the factor that should hit the first increment is actually $W(t_0)=0$. This is why there is that zero in the first entry. (I actually would say this is bad notation at least for the Ito case and that it would be better to have this $0$ inside the $W$ vector i.e. W=[0 cumsum(dW(1:end-1))] but it's not a big deal.)

That strat line actually involves much more math to see why it is correct. The idea is that if you know $W(s)=a$ and $W(t)=b$ with $s<t$ then $W(s+u(t-s))$ for $u \in [0,1]$ is a Brownian bridge. A single time sample from a Brownian bridge is normally distributed with mean $a+u(b-a)$ and variance $u(1-u)(t-s)$, which is ultimately what they implemented in that line, though it is somewhat obfuscated that this is what is happening.

I think a more obviously correct implementation of what they achieved in the code snippet above for the Stratonovich integral would be:

N=1000; dt=T/N;
dW_fine=sqrt(dt)*randn(1,N); 
dW_coarse=dW_fine(1:2:end)+dW_fine(2:2:end); %N must be even
W=cumsum(dW_fine);
strat=dW_coarse' * W(1:2:end)

This costs more memory but doesn't really cost significantly more CPU time.