확률 미분 방정식 시뮬레이션

1 변수의 확률 미분 방정식은 일반적으로 다음과 같이 표현됩니다.

$$
dX (t) = f (X (t)) dt + g (X (t)) dW (t),\\X (0) = X_0, 0\leq t\leq T\\.
$$

여기서 $f,g$는 스칼라 함수, $W(t)$는 위너 과정입니다. 이 방정식의 수치 계산을 해 보겠습니다.

수치 계산을 하기 위해서는, 이산화를 할 필요가 있습니다. 확률 미분 방정식의 이산화에는 크게 두 가지 방법이 있으며, 각각 Euler-Maruyama 방법과 Milstein 방법이라고합니다. 이번에는 정확도가 떨어지지만 간단한 Euler-Maruyama 방법을 사용하여 계산합니다.

먼저 구간 $[0,T]$를 이산화합니다. 꽤 큰 양의 정수 $N$를 사용하여, $\Delta t := T/N$, $\tau_j := j\Delta t$로 합니다. 또한 $X_j:=X(\tau_j)$로 표현하기로 합니다.

Euler-Maruyama 방법은 다음 형식으로 방정식을 이산화합니다.

$$
X_j = X_{j-1}+f(X_{j-1})\Delta t + g(X_{j-1}) (W(\tau_j)-W(\tau_{j-1})), j = 1,2,\cdots L\\.
$$

여기서 $W(t)$는 위너 과정이므로 $W(\tau_j)-W(\tau_{j-1})$ 는 평균 0, 분산 $\tau_j-\tau_{j-1 }$ 의 정규 분포 $N(0,\tau_j -\tau_{j-1})$ 에 따른 확률 변수입니다.

Euler-Maruyama 방법을 사용하여 다음 방정식으로 표시되는 기하학적 브라운 운동을 계산해보십시오.

dX(t) = \sigma X(t) dt +  \mu X(t) dW(t), \ \ X(0)=X_0 \  \ .

여기서 $\sigma,\mu$는 상수입니다. 기하학적 브라운 운동은 해석 해를 구할 수 있으며 해는 다음과 같습니다.

X(t) = X(0) \exp ( (\sigma-\frac{1}{2}\mu^2)t+\mu W(t) ) \ \ .

Euler-Maruyama법의 결과와 해석해의 결과를 비교하는 코드는 다음과 같습니다.

import numpy as np
import matplotlib.pyplot as plt
import seaborn as sns

sigma = 2
mu = 1
Xzero = 1

T = 1.0
N = 2**8
dt = T / N

t = np.arange(0.0, T, dt)
dW = np.sqrt(dt) *  np.random.randn(N)
dW[0] = 0
W = np.cumsum(dW)
Xtrue = Xzero * np.exp(  (sigma - 0.5 * mu**2) * t +(mu*W)  )

Xem = np.zeros(N)
Xtemp = Xzero
Xem[0] = Xtemp
for j in range(1, N):
    Xtemp = Xtemp + dt * sigma * Xtemp + mu * Xtemp * dW[j]
    Xem[j] = Xtemp

plt.plot(t, Xtrue, label='Analytic')
plt.plot(t, Xem, label='Euler-Maruyama')
plt.legend()

참고문헌

Learning SDEs in Python