Ito's formula for $B^2$ doesn't work with Monte Carlo

Question

I am trying to "verify" Ito's formula for $B^2$ with simple Monte Carlo method.
According to my textbook, $ B^2 $ satisfies this SDE: $d(B^2) = dt + 2BdB$.

My simulation (notebook) results in following distribution for $B^2$ with $dt \in \{1/100,1/1000,1/10000\}$:

And distribution for solution of $dy = dt + 2BdB$:

Clearly they are different and it looks like I am making some very simple mistake, am I?

Here's my code in Julia (notebook):

# k instances of ito integral simulations of f from t=0. to t=1.
function ito_sim(f, n, k)
    y   = zeros(k)
    dt  = 1./n
    sdt = sqrt(dt)
    t = 0.
    b = zeros(k)
    for i in 1:n
        db = randn(k) * sdt
        b += db
        t += dt
        y += f(t, dt, b, db) 
    end
    y
end;
plt[:hist](ito_sim((t,dt,b,db)-> db, n, 10000) .^ 2);
plt[:hist](ito_sim((t,dt,b,db)-> dt + 2 * b .* db, n, 10000));

PS I searched through stackoverflow and it seem to me that math.stachexchange is better suited for this question.

Answer 1

Looking at your code , your scheme is biased , and it is because of the way you calculate the Ito integral.

You built a function that calculate samples of an Ito process at time $1$. Let $g(t,B_t)$ an Ito process defined as

$$dg(t,B_t)=a(t,B_t)dt+\sigma(t,B_t)dB_t$$ where $B_t$ is a Brownian motion.

You say that given a subdivision $\{t_i\}$, defined as $t_i=\frac{i}{n}$

$$g(1,B_1)=\sum_{k=1}^{n}{[g(t_k,B_{t_k})-g(t_{k-1},B_{t_{k-1}})]}\approx\sum_{k=1}^{n}{[a(t_k,B_{t_{k}})(t_k-t_{k-1})+\sigma(t_k,B_{t_{k}})(B_{t_k}-B_{t_{k-1}})]}$$

Then we define a function $$f(t,dt,b,db)=[a(t,b)dt+\sigma(t,b)db]$$

$$g(1,B_1)\approx\sum_{k=1}^{n}{f(t_k,t_k-t_{k-1},B_{t_{k}},B_{t_{k}}-B_{t_{k-1}})}$$

The above formula can be a good approximation if $g(0,B_0)=0$, which is the case here. When you apply your code to the process $$dg(t,B_t)=dt+2B_tdB_t$$

Your final result will look like

$$X=\sum_{k=1}^{n}{[\frac{1}{n}+2B_{t_{k}}(B_{t_k}-B_{t_{k-1}})]}$$ and its expectation is $3$, while you expect $1$.

The reason is that you use this approximation in your code $$\int_{t_k}^{t_{k+1}}{B_tdB_t}\approx B_{t_{k+1}}(B_{t_{k+1}}-B_{t_{k}})$$

If you had chosen

$$\int_{t_k}^{t_{k+1}}{B_tdB_t}\approx B_{t_{k}}(B_{t_{k+1}}-B_{t_{k}})$$

the bias is removed, and the expectation would have been $1$.