Local timeの数値計算
2021/10/13
#math #stochastic analysis #python
JST, UTCとかではなく確率解析のLocal time

TL;DR

  • 田中の公式を使う
  • Lta(X)=∣Xt−a∣−∣X0−a∣−∫0tsgn(Xs−a)dXsL_t^a(X) = |X_t - a| - |X_0 - a| - \int_0^t \mathrm{sgn}(X_s - a)dX_s

Local time

Local timeは以下のように定義される.

  • WtW_t: ブラウン運動
  • LtaL_t^a: Level aaでのLocal time
Lta=lim⁡ϵ→012ϵ∫0t1{x−ϵ≤Ws≤x+ϵ}dsL_t^a = \lim_{\epsilon \to 0}\frac{1}{2\epsilon}\int_0^t \mathbb{1}_{\{x - \epsilon \le W_s \le x + \epsilon \}}ds

Informalにはデルタ関数を使って次のようにも書ける.

Lta=∫0tδa(Ws)dsL_t^a = \int_0^t \delta_a(W_s)ds

つまりLtaL_t^aはブラウン運動WWが時点ttまでにaaに滞在した合計時間を表している.

上の定義はブラウン運動に対するものだが, Semimartingaleに対してもLocal timeを定義できる.

数値計算

上の定義を直接数値計算するのは現実的ではない. そこで以下の田中の公式を用いる.

田中の公式

  • XtX_t: Semimartingale
  • a∈Ra \in \mathbb{R}
  • LtaL_t^a: Level aaでのLocal time
  • sgn(x)\mathrm{sgn}(x): xxが正のとき1, 負のとき-1を返す符号関数
∣Xt−a∣=∣X0−a∣+∫0tsgn(Xs−a)dXs+Lta(X)|X_t - a| = |X_0 - a| + \int_0^t \mathrm{sgn}(X_s - a)dX_s + L_t^a(X)

これを式変形すれば

Lta(X)=∣Xt−a∣−∣X0−a∣−∫0tsgn(Xs−a)dXsL_t^a(X) = |X_t - a| - |X_0 - a| - \int_0^t \mathrm{sgn}(X_s - a)dX_s

これで数値計算ができる. 右辺の伊藤積分はオイラー丸山法とかで数値計算すれば良い.

Examples

以下はPythonでの数値計算例.

ブラウン運動の場合

import numpy as np
import matplotlib.pyplot as plt


# 標準ブラウン運動のパスを生成して返す
def standard_brownian(dt: float, T: float) -> tuple[np.ndarray, np.ndarray]:
    N = int(T / dt)
    t = np.arange(0, T, step=dt)
    t = np.insert(t, len(t), T)

    dB = np.sqrt(dt) * np.random.randn(N)
    return t, np.insert(dB.cumsum(), 0, 0)


# 幾何ブラウン運動のパスを生成して返す
def geometric_brownian(dt: float, T: float, S0: float = 1, mu: float = 0, sigma: float = 1) -> tuple[np.ndarray, np.ndarray]:
    t, B = standard_brownian(dt, T)
    return t, S0 * np.exp((mu - 0.5*sigma**2)*t + sigma*B)


# 伊藤積分を計算する
# f: 被積分関数
# X: 確率過程
# I: Integrator
def ito_integral(f, X: np.ndarray, I: np.ndarray) -> np.ndarray:
    ret = np.zeros_like(X)
    for i in range(1, len(X)):
        dI = I[i] - I[i - 1]
        ret[i] = ret[i - 1] + f(X[i - 1]) * dI
    return ret


if __name__ == "__main__":
    dt = 0.001 # 時間分割幅
    T = 1
    a = 0

    t, X = standard_brownian(dt, T)
    L_t = np.abs(X - a) - np.abs(X[0] - a) - ito_integral(lambda x: np.sign(x - a), X, X)

    plt.plot(t, X, label="$W_t$")
    plt.plot(t, L_t, label="$L_t^0$")
    plt.xlabel("$t$")
    plt.grid()
    plt.legend()
    plt.show()
ブラウン運動での数値計算結果
ブラウン運動での数値計算結果

WtW_tが0付近を通る度にLt0L_t^0が増加していることが分かる. ちゃんと計算出来ていそう.

幾何ブラウン運動の場合

if __name__ == "__main__":
    dt = 0.001
    T = 1
    a = 1

    t, X = geometric_brownian(dt, T)
    L_t = np.abs(X - a) - np.abs(X[0] - a) - ito_integral(lambda x: np.sign(x - a), X, X)

    plt.plot(t, X, label="$W_t$")
    plt.plot(t, L_t, label="$L_t^1$")
    plt.xlabel("$t$")
    plt.grid()
    plt.legend()
    plt.show()
幾何ブラウン運動での数値計算結果
幾何ブラウン運動での数値計算結果

幾何ブラウン運動の場合でも1付近を通る度にLt1L_t^1が増加していることが分かる. 式では分かっていても実際に数値計算してちゃんと結果が出てくるのが面白い.

References