# Approximation of an Itô integral with python

Exercise 3.11 (Approximation of an Itô Integral). In this example, the stochastic integral $$\int^t_0tW(t)dW(t)$$ is considered. The expected value of the integral and the expected value of the square of the integral are estimated using $$M=1,000,000$$ sample paths. The number of intervals used on $$[0,1]$$ is $$N$$ where $$N=4,8,16,\dots,256$$. The preceding approximate method is used, which for this problem, has the explicit form \begin{align*} E((I(f))^\alpha)&=E\left(\int^1_0tW(t)dW(t)\right)^\alpha\\ &\approx\frac{1}{M}\sum^M_{j=1}\left(\sum^{N-1}_{i=0}t_i^{(N)}W^{(j)}(t_i^{(N)})(W^{(j)}(t^{(N)}_{t+1})-W^{(j)}(t^{(N)}_i))\right)^\alpha\quad(2)\\ \end{align*} for $$\alpha=1$$ or $$2$$ and $$t_i^{(N)}=i/N$$ for $$i=0,1,2,\dots,N$$. Notice that $$(W^{(j)}(t^{(N)}_{i+1})-W^{(j)}(t_i^{(N)}))=\eta_i^{(j)}/\sqrt{N}$$ where $$\eta_i^{(j)}\sim N(0,1)$$ and also that $$(W^{(j)}(t^{(N)}_{i})-W^{(j)}(t_{i-1}^{(N)}))=\eta_{i-1}^{(j)}/\sqrt{N}$$. A computer program that performs this calculation is listed at the end of this chapter. The calculational results are given in Table 3.2 for $$M=1,000,000$$ sample paths. Notice that the results improve as $$N$$ increases. The exact values are $$E(I(f))=0$$ and $$E((I(f))^2)=0.25.$$ $$\begin{array}{ccc} \hline \text{Value of N}&|E(I(f))|&|E((I(f))^2)|\\ \hline \color{red}{2^2}&\color{red}{0.00107}&\color{red}{0.14071}\\ 2^3&0.00012&0.19151\\ 2^4&0.00069&0.21906\\ 2^5&0.00038&0.23508\\ 2^6&0.00007&0.24209\\ 2^7&0.00001&0.24486\\ 2^8&0.00002&0.24821\\ \hline \end{array}$$ Source: Modeling with Itô Stochastic Differential Equations by Edward Allen, pg. 74

## My Question

Per Abezhiko's wonderful reply in my other question, one has $$\mathbb{E}[X_t] = 0$$, since $$W_s$$ and $$\mathrm{d}W_s$$ are independent and have zero mean, while Itô's isometry gives $$\mathrm{Var}[X_t] = \mathbb{E}[X_t^2] = \mathbb{E}\left[\int_0^t s^2W_s^2 \,\mathrm{d}s\right] = \int_0^t s^2\,\mathbb{E}[W_s^2] \,\mathrm{d}s = \int_0^t s^3 \,\mathrm{d}s = \frac{1}{4}t^4$$ which is equal to $$1/4$$ when $$t=1$$.

I've provided the Fortran code here from the textbook which was provided by the author, but I'd like to run this stochastic integral in python. There is a code which I used from here, to get the bottom values in the table, but the output is nonsense:

In [5]:
T = 1.0 # end time
N = 1000000 # number of steps, here it's
dt = T / N

dW = np.sqrt(dt) * np.random.randn(N)
W = np.cumsum(dW)

shiftedW = np.zeros(N)
shiftedW[1:] = W[:-1]

# The two different integrals.
# The Ito integral is roughly the Riemann integral evaluated at the left edge of the subinterval
ito = np.sum(shiftedW*dW)
# The Stratonovich integral is roughly the Riemann integral evaluated at the centre of the subinterval
stratonovich = np.sum((0.5*(shiftedW+W) + 0.5*np.sqrt(dt)*np.random.randn(N))*dW)

# Note that the exact solutions are different - markedly so!
print("The Ito integral is {} with error {}".format(ito, np.abs(ito - 0.5*(W[-1]**2-T))))
print("The Stratonovich integral is {} with error {}".format(stratonovich, np.abs(stratonovich - 0.5*W[-1]**2)))

## The Ito integral is 2.264544487686194 with error 0.0008807299528790224
## The Stratonovich integral is 2.765668595686916 with error 0.00024337804784302364

plt.plot(t, W)
plt.xlabel(r'$$t$$')
plt.ylabel(r'$$W(t)$$')
plt.savefig('itointegral', dpi=300, bbox_inches='tight')
plt.show()


Would anybody be able to provide a hint?

it seems that your coding is not calculating the expected value of the Ito integral but instead, it is only calculating the Ito integral for one trajectory. You need to add a loop to your coding for looping through the number of simulations. So far, you just have coded the logic for simulating one trajectory. I have done a quick modification to your coding and I retrieve the theoretical variance of the Ito integral:


import numpy as np
import scipy as sp
import matplotlib.pyplot as plt

T = 1.0 # end time
N = 100000 # number of steps, here it's
dt = T / N
Ito = 0.
nbSim = 10000
t = np.linspace(0,1,N)

for j in range(nbSim):
dW = np.sqrt(dt) * np.random.randn(N)
W = np.cumsum(dW)

shiftedW = np.zeros(N)
shiftedW[1:] = W[:-1]

# The Ito integral is roughly the Riemann integral evaluated at the left edge of the subinterval
ito_trajectory = np.sum(t*shiftedW*dW)

Ito = Ito + ito_trajectory **2

Ito = Ito / nbSim

# Note that the exact solutions are different - markedly so!
print("The Ito integral is {} ".format(Ito))
'''

• Many thanks! Would there be a way to code all of the trajectories, or would this be not particularly useful? Also, if I wanted to plot the variance via python, would this code be correct? plt.plot(t, shiftedW) plt.xlabel(r'$t$') plt.ylabel(r'$W(t)$') plt.savefig('itointegral', dpi=300, bbox_inches='tight') plt.show() I noticed the trajectory remains the same for this code. Dec 30, 2023 at 15:12