9.8. O método de Milstein

Vimos que o método de Euler-Maruyama para equações estocásticas multiplicativas é de ordem forte \(1/2\), ao contrário do método de Euler clássico, que é de ordem \(1\). É natural buscarmos de métodos de ordem mais alta também para equações estocásticas. Vamos ver, aqui, o método de Milstein, que é de ordem forte \(1\).

Para simplificar, vamos considerar equações diferenciais estocásticas autônomas, da forma

\[ \mathrm{d}X_t = f(X_t)\;\mathrm{d}t + g(X_t)\;\mathrm{d}W_t. \]

Mas o caso não autônomo pode ser tratado de forma análoga. Vamos considerar, ainda, uma condição inicial

\[ \left. X_t \right|_{t = 0} = X_0. \]

Expansão de Taylor estocástica

A equação diferencial estocástica é interpretada como significando

\[ X_t = X_0 + \int_0^t f(X_s)\;\mathrm{d}s + \int_0^t g(X_s)\;\mathrm{d}W_s. \]

A partir de um determinado ponto \(t_0 \geq 0\), podemos escrever

\[ X_t = X_{t_0} + \int_{t_0}^t f(X_s)\;\mathrm{d}s + \int_{t_0}^t g(X_s)\;\mathrm{d}W_s. \]

A ideia é usar a isometria de Itô para reescrever os integrandos acima e obter uma expansão da solução nos pontos da malha.

Pela isometria de Itô, temos

\[ \mathrm{d}f(X_t) = f'(X_t)\;\mathrm{d}X_t + \frac{1}{2}f''(X_t)\;\mathrm{d}t = f'(X_t)f(X_t)\;\mathrm{d}t + f'(X_t)g(X_t)\;\mathrm{d}W_t + \frac{1}{2}f''(X_t)g(X_t)^2\;\mathrm{d}t. \]

e

\[ \mathrm{d}g(X_t) = g'(X_t)\mathrm{d}X_t + \frac{1}{2}g''(X_t)\;\mathrm{d}t = g'(X_t)f(X_t)\;\mathrm{d}t + g'(X_t)g(X_t)\;\mathrm{d}W_t + \frac{1}{2}g''(X_t)g(X_t)^2\;\mathrm{d}t. \]

A partir do ponto \(t_0\), obtemos

\[ f(X_s) = f(X_{t_0}) + \int_{t_0}^s f'(X_s)f(X_s)\;\mathrm{d}s + \int_{t_0}^s f'(X_s)g(X_s)\;\mathrm{d}W_s + \int_{t_0}^s \frac{1}{2}f''(X_s)g(X_s)^2\;\mathrm{d}s \]

e

\[ g(X_s) = g(X_{t_0}) + \int_{t_0}^s g'(X_s)f(X_s)\;\mathrm{d}s + \int_{t_0}^s g'(X_s)g(X_s)\;\mathrm{d}W_s + \int_{t_0}^s \frac{1}{2}g''(X_s)g(X_s)^2\;\mathrm{d}s. \]

Substituindo essas expressões para \(f(X_s)\) e \(g(X_s)\) na equação integral para \(\{X_t\}_t\), obtemos

\[ X_t = X_{t_0} + \int_{t_0}^t f(X_{t_0})\;\mathrm{d}s + \int_{t_0}^t g(X_{t_0})\;\mathrm{d}W_s + R_{t_0, t}, \]

onde o resto \(R_{t_0, t}\) é dado por

\[ R_{t_0, t} = R_{t_0, t}^f + R_{t_0, t}^g, \]

com

\[ R_{t_0, t}^f = \int_{t_0}^{t} \left(\int_{t_0}^s \left(f'(X_\tau)f(X_\tau) + \frac{1}{2}f''(X_\tau)g(X_\tau)^2\right)\;\mathrm{d}\tau + \int_{t_0}^\tau f'(X_\tau)g(X_\tau)\;\mathrm{d}W_\tau\right)\;\mathrm{d}s \]

e

\[ R_{t_0, t}^g = \int_{t_0}^{t} \left(\int_{t_0}^s \left(g'(X_\tau)f(X_\tau) + \frac{1}{2}g''(X_\tau)g(X_\tau)^2\right)\;\mathrm{d}\tau + \int_{t_0}^s g'(X_\tau)g(X_\tau)\;\mathrm{d}W_\tau\right)\;\mathrm{d}W_s. \]

Cada termo tem uma determinada ordem de grandeza em termos de \(\mathrm{d}t\) e \(\mathrm{d}W_t\). Como \(\mathrm{d}W_t \sim \sqrt{\mathrm{d}t}\), as diferentes integrais presentes na expansão não têm a mesma ordem. Podemos dizer que, em uma aproximação com passo \(\Delta t = t - t_0\), as integrais simples tem ordem

\[ \int_{t_0}^t f(X_s)\;\mathrm{d}s \sim \Delta t \]

e

\[ \int_{t_0}^t g(X_s)\;\mathrm{d}W_s \sim \Delta t^{1/2}. \]

Por sua vez, as integrais duplas presentes na expansão desses dois termos têm ordem

\[ \begin{align*} & \int_{t_0}^t \int_{t_0}^s F(X_\tau) \;\mathrm{d}\tau \;\mathrm{d}s \sim \Delta t^2. \\ & \int_{t_0}^t \int_{t_0}^s F(X_\tau) \;\mathrm{d}W_\tau \;\mathrm{d}s \sim \Delta t^{3/2} \\ & \int_{t_0}^t \int_{t_0}^s F(X_\tau) \;\mathrm{d}W_\tau \;\mathrm{d}W_s \sim \Delta t. \end{align*} \]

Assim, podemos, de acordo com o objetivo, expandir \(F(X_\tau)\) apenas em alguns termos, de ordem mais baixa, começando pela integral estocástica dupla.

Expansão nos pontos da malha

Em pontos \(t_j = jT/n\), \(j = 0, 1, \ldots, n\), podemos escrever

\[ X_{t_j} = X_{t_{j-1}} + \int_{t_{j-1}}^{t_j} f(X_s)\;\mathrm{d}s + \int_{t_{j-1}}^{t_j} g(X_s)\;\mathrm{d}W_s. \]

Substituindo \(f(X_s)\) e \(g(X_s)\) de acordo com a fórmula de Itô a partir de um ponto \(t_{j-1}\) da malha, obtemos

\[ X_{t_j} = X_{t_{j-1}} + \int_{t_{j-1}}^{t_j} f(X_{t_{j-1}})\;\mathrm{d}s + \int_{t_{j-1}}^{t_j} g(X_{t_{j-1}})\;\mathrm{d}W_s + R_j, \]

onde o resto \(R_j\) é dado por

\[ \begin{align*} R_j & = \int_{t_{j-1}}^{t_j} \left(\int_{t_{j-1}}^s \left(f'(X_\tau)f(X_\tau) + \frac{1}{2}f''(X_\tau)g(X_\tau)^2\right)\;\mathrm{d}\tau + \int_{t_{j-1}}^s f'(X_\tau)g(X_\tau)\;\mathrm{d}W_\tau\right)\;\mathrm{d}s \\ & \qquad + \int_{t_{j-1}}^{t_j} \left(\int_{t_{j-1}}^s \left(g'(X_\tau)f(X_\tau) + \frac{1}{2}g''(X_\tau)g(X_\tau)^2\right)\;\mathrm{d}\tau + \int_{t_{j-1}}^s g'(X_\tau)g(X_\tau)\;\mathrm{d}W_\tau\right)\;\mathrm{d}W_s. \end{align*} \]

Revisitando o método de Euler-Maruyama

Observe que a aproximação para o método de Euler-Maruyama é obtida ao descartarmos completamente o resto \(R_j\), ficando, apenas

\[ X_{t_j} \approx X_{t_{j-1}} + \int_{t_{j-1}}^{t_j} f(X_{t_{j-1}})\;\mathrm{d}s + \int_{t_{j-1}}^{t_j} g(X_{t_{j-1}})\;\mathrm{d}W_s = X_{t_{j-1}} + f(X_{t_{j-1}})\Delta t + g(X_{t_{j-1}}) \Delta W_{j-1}, \]

onde \(\Delta W_{j-1} = W_{t_j} - W_{t_{j-1}}\).

Isso nos leva ao método de Euler-Maruyama

\[ X_j = X_{j-1} + f(X_{j-1}) \Delta t + g(X_{j-1})\Delta W_{j-1}, \qquad j = 1, \ldots, n, \]

com \(X_0\) dado.

Em termos da aproximação, podemos escrever

\[ \int_{t_0}^t f(X_s)\;\mathrm{d}s = f(X_{t_0})\Delta t + \mathcal{O}(\Delta t^{3/2}) \]

e

\[ \int_{t_0}^t g(X_s)\;\mathrm{d}W_s = g(X_{t_0})\Delta t^{1/2} + \mathcal{O}(\Delta t) \]

A integral em \(f\) tem uma ordem mais alta do que a em \(g\), mas não podemos descartá-la por completo, pois aproximação em \(f\) cairia para ordem zero, na verdade (Pense em aproximar a função \(f(t) = t = \int_0^t \;\mathrm{d}s\) por zero, \(\tilde f(t) = 0\), em que cada passo \(\int_t^{\Delta t} \;\mathrm{d}s = \Delta t\) mas, após a integração em um intervalo \(0\leq t \leq 1,\) o erro é \(f(1) - \tilde f(1) = 1 - 0 = 1\)). Temos que aproximar os dois termos até uma ordem mínima desejada.

O método de Milstein

Para o método de Milstein, vamos reter o termo de ordem mais baixa do resto \(R_j\), que é o termo com a integral estocástica dupla:

\[ \int_{t_{j-1}}^{t_j} \int_{t_{j-1}}^s g'(X_\tau)g(X_\tau)\;\mathrm{d}W_\tau\;\mathrm{d}W_s. \]

Considerando \(u(x) = g'(x)g(x)\) e aplicando a isometria de Itô a \(\{u(X_t)\}_{t \geq 0}\), obtemos

\[ g'(X_\tau)g(X_\tau) = g'(X_{t_{j-1}})g(X_{t_{j-1}}) + \int_{t_{j-1}}^\tau \mathrm{d}u(X_\eta), \]

com

\[ \begin{align*} \mathrm{d}u(X_t) & = u'(X_t)\;\mathrm{d}X_t + \frac{1}{2}u''(X_t)g(X_t)^2\;\mathrm{d}t \\ & = (g''(X_t)g(X_t) + g'(X_t)^2)f(X_t)\;\mathrm{d}t \\ & \qquad + (g''(X_t)g(X_t) + g'(X_t)^2)g(X_t)\;\mathrm{d}W_t \\ & \qquad + \frac{1}{2}(g'''(X_t)g(X_t) + 3g''(X_t)g'(X_t) + g'(X_t)^3)g(X_t)^2\;\mathrm{d}t. \end{align*} \]

Isso implica em

\[ \int_{t_{j-1}}^{t_j} \int_{t_{j-1}}^s g'(X_\tau)g(X_\tau)\;\mathrm{d}W_\tau\;\mathrm{d}W_s = \int_{t_{j-1}}^{t_j} \int_{t_{j-1}}^s g'(X_{t_{j-1}})g(X_{t_{j-1}})\;\mathrm{d}W_\tau\;\mathrm{d}W_s + \int_{t_{j-1}}^{t_j} \int_{t_{j-1}}^s \int_{t_{j-1}}^\tau \;\mathrm{d}u(X_\eta)\;\mathrm{d}W_\tau\;\mathrm{d}W_s, \]

que nos leva à expansão

\[ X_{t_j} = X_{t_{j-1}} + \int_{t_{j-1}}^{t_j} f(X_{t_{j-1}})\;\mathrm{d}s + \int_{t_{j-1}}^{t_j} g(X_{t_{j-1}})\;\mathrm{d}W_s + g'(X_{t_{j-1}})g(X_{t_{j-1}})\int_{t_{j-1}}^{t_j} \int_{t_{j-1}}^s \;\mathrm{d}W_\tau\;\mathrm{d}W_s + \tilde R_j, \]

com um resto \(\tilde R_j\) contendo a integral tripla no lugar da integral estocástica dupla. Descartando o resto \(\tilde R_j\), obtemos o aproximação utilizada no método de Milstein

\[ \begin{align*} X_{t_j} & \approx X_{t_{j-1}} + \int_{t_{j-1}}^{t_j} f(X_{j-1})\;\mathrm{d}s + \int_{t_{j-1}}^{t_j} g(X_{j-1})\;\mathrm{d}W_s + \int_{t_{j-1}}^{t_j} \int_{t_{j-1}}^s g'(X_{t_{j-1}})g(X_{t_{j-1}})\;\mathrm{d}W_\tau\;\mathrm{d}W_s \\ & = X_{t_{j-1}} + f(X_{j-1})\Delta t + g(X_{j-1}) (W_{t_j} - W_{t_{j-1}}) + g'(X_{t_{j-1}})g(X_{t_{j-1}})\int_{t_{j-1}}^{t_j} \int_{t_{j-1}}^s \;\mathrm{d}W_\tau\;\mathrm{d}W_s. \end{align*} \]

Nesse caso em particular, a integral estocástica dupla pode ser calculada mais explicitamente. De fato, a integral interior é simplesmente

\[ \int_{t_{j-1}}^s \;\mathrm{d}W_\tau = W_s - W_{t_{j-1}}. \]

Assim,

\[ \int_{t_{j-1}}^{t_j} \int_{t_{j-1}}^s \;\mathrm{d}W_\tau\;\mathrm{d}W_s = \int_{t_{j-1}}^{t_j} (W_s - W_{t_{j-1}})\;\mathrm{d}W_s. \]

Podemos usar a propriedade de invariância por translação dos processos de Wiener, que diz que \(\{\tilde W_\tau\}_{\tau \geq 0}\) dado por \(\tilde W_\tau = W_{t_{j-1} + \tau} - W_{t_{j-1}}\) também é um processo de Wiener, para escrever

\[ \int_{t_{j-1}}^{t_j} (W_s - W_{t_{j-1}})\;\mathrm{d}W_s = \int_{0}^{\Delta t} \tilde W_\tau\;\mathrm{d}\tilde W_\tau = \frac{\tilde W_{\Delta t}^2}{2} - \frac{\Delta t}{2} = \frac{(W_{t_j} - W_{t_{j-1}})^2}{2} - \frac{\Delta t}{2}. \]

Ou, mais diretamente,

\[ \begin{align*} \int_{t_{j-1}}^{t_j} \int_{t_{j-1}}^s \;\mathrm{d}W_\tau\;\mathrm{d}W_s & = \int_{t_{j-1}}^{t_j} (W_s - W_{t_{j-1}})\;\mathrm{d}W_s \\ & = \int_{t_{j-1}}^{t_j} W_s \;\mathrm{d}W_s - \int_{t_{j-1}}^{t_j} W_{t_{j-1}}\;\mathrm{d}W_s \\ & = \frac{1}{2}\left(W_{t_j}^2 - W_{t_{j-1}}^2 - t_j + t_{j-1} \right) - W_{t_{j-1}}\left(W_{t_j} - W_{t_{j-1}}\right) \\ & = \frac{1}{2}\left((W_{t_j} + W_{t_{j-1}})(W_{t_j} - W_{t_{j-1}}) - t_j + t_{j-1} \right) - W_{t_{j-1}}\left(W_{t_j} - W_{t_{j-1}}\right) \\ & = \frac{1}{2}\left((W_{t_j} + W_{t_{j-1}})(W_{t_j} - W_{t_{j-1}}) - 2W_{t_{j-1}}\left(W_{t_j} - W_{t_{j-1}}\right) \right) - \frac{\Delta t}{2} \\ & = \frac{1}{2}\left(W_{t_j} - W_{t_{j-1}}\right)^2 - \frac{\Delta t}{2} \end{align*} \]

De um jeito ou de outro, podemos escrever

\[ \int_{t_{j-1}}^{t_j} (W_s - W_{t_{j-1}})\;\mathrm{d}W_s = \frac{\Delta W_{j-1}^2}{2} - \frac{\Delta t}{2}, \]

onde \(\Delta W_{j-1} = W_{t_j} - W_{t_{j-1}}\). Dessa forma, obtemos aproximação

\[ X_{t_j} \approx X_{t_{j-1}} + f(X_{j-1})\Delta t + g(X_{j-1}) \Delta W_{j-1} + \frac{1}{2}g'(X_{t_{j-1}})g(X_{t_{j-1}})\left(\Delta W_{j-1}^2 - \Delta t\right), \]

Isso nos leva ao método de Milstein

\[ X_j = X_{j-1} + f(X_{j-1}) \Delta t + g(X_{j-1})\Delta W_{j-1} + \frac{1}{2}g'(X_{j-1})g(X_{j-1})\left(\Delta W_{j-1}^2 - \Delta t\right), \qquad j = 1, \ldots, n, \]

com \(X_0\) dado. Mais explicitamente, como \(\Delta W_{j-1} \sim \sqrt{\Delta t} \mathcal{N}(0, 1)\), escrevemos

\[ X_j = X_{j-1} + f(X_{j-1}) \Delta t + g(X_{j-1})Z_{j-1}\sqrt{\Delta t} + \frac{1}{2}g'(X_{j-1})g(X_{j-1})\left(Z_{j-1}^2 - 1 \right)\Delta t, \qquad j = 1, \ldots, n, \]

onde \(Z_i \sim \mathcal{N}(0, 1)\) são independentes.

Ordem de convergência do método de Milstein

Não faremos as estimativas aqui, mas, como dito acima, pode-se mostrar que o método de Milstein tem ordem forte \(1\) e ordem fraca \(1\), também.