確率積分(Ito integral)を算出してみた

■ 確率積分として以下のような方式があります

Riemann sum
式)
$ \sum_{j=0}^{N-1} h(t_j) (t_{j+1}-t_j) $

Ito integral.
伊藤清（Kiyoshi Ito)によって拡張された確率過程の方式
式)
$ \sum_{j=0}^{N-1}h(t_j) (W(t_{j+1})-W(t_j)) $.
<=> $ \int_0^{T} h(t)dW(t) $.

Stratonovich integrals
式)
$\sum_{j=0}^{N-1}h\big(\frac{t_j+t_{j+1}}{2}\big)(W(t_{j+1})-W(t_j))$

　Ito integralの計算

$h(t)\equiv W(t)$の時、Ito integralは、

\begin{align*}
\sum_{j=0}^{N-1} W(t_j) (W(t_{j+1})-W(t_j)) \\
&=\sum_{j=0}^{N-1} ({W(t_{j+1})}^2 - {W(t_{j+1})}^2 + 2 W(t_j) W(t_{j+1}) - {W(t_j)}^2  - {W(t_j)}^2) \\
&= \frac{1}{2} \sum_{j=0}^{N-1} ({W(t_{j+1})}^2 - {W(t_j)}^2 - (W(t_{j+1}) - W(t_j))^2 )\\
&= \frac{1}{2} (W(T)^2 - W(0)^2)  -  \frac{1}{2} \sum_{j=0}^{N-1} (W(t_{j+1}) - W(t_j))^2
\end{align*}

$\sum_{j=0}^{N-1} (W(t_{j+1}) - W(t_j))^2$はWiener Processの分散に等しく、$T$とおけるため、
$$\sum_{j=0}^{N-1} (W(t_{j+1}) - W(t_j))^2=\frac{1}{2}(W(T)^2-T)$$と表現できる。

Ito integralを算出してみる。

N=10000;
M=1;
T = 1.0
dt = T / N;

t = np.arange(0.0,1.0, dt);
dW = np.sqrt(dt)*randn(N,M); # (N, M)行列
W = np.cumsum(dW,axis=0);
e = np.array([[0]])
W_=np.concatenate((e,W[:-1]), axis=0)  # 初期値に0をおく　 W(0)として。

ito = np.dot(W_.T,dW) # (10000, 1).T * (10000, 1)  => (1, 1)
np.abs(ito - 0.5*(W[-1]**2 - T) )[0][0]

## output
##  0.00414652405737

参考記事

• Desmond J. Higham "An Algorithmic Introduction to Numerical Simulation of Stochastic Differential Equations"

• Wikipedia「伊藤清」 https://ja.wikipedia.org/wiki/%E4%BC%8A%E8%97%A4%E6%B8%85
• Learning SDEs in Python
  https://www.quantopian.com/posts/learning-sdes-in-python