9.5. Simulações ilustrando ordem de convergência do método de Euler-Maruyama

Podemos ilustrar e conferir a order de convergência do método de Euler-Maruyama considerando equações diferenciais estocásticas para as quais conhemos explicitamente a solução exata.

Vamos considerar duas equações, uma com ruído multiplicativo, que teoricamente nos dá uma convergência forte de ordem \(1/2\), e outra com ruído aditivo, que teoricamente nos dá uma convergência forte de ordem \(1\). Em ambas temos, também, convergência fraca de ordem \(1\).

O caso com ruído multiplicativo é ilustrado com o movimento Browniano geométrico, enquanto que o caso com ruído aditivo é ilustrado com o processo de Orstein-Uhlenbeck.

Convergência forte de ordem 1/2 no movimento browniano geométrico

Nesse caso, a equação do movimento browniano geométrico é

\[ \mathrm{d}X_t = \mu X_t \;\mathrm{d}t + \sigma X_t \;\mathrm{d}W_t, \]

cuja solução é

\[ X_t = X_0 e^{(\mu + \sigma^2/2)t + \sigma W_t}. \]

Aproximamos a solução pelo método de Euler-Maruyama, em um intervalo de tempo \([0, T]\), dado por

\[ X_j^n = X_{j-1}^n + \mu X_{j-1}^n\Delta t + \sigma X_{j-1}^n \Delta W_j, \quad j = 1, \ldots, n, \]

onde

\[ n\in \mathbb{N}, \quad \Delta j = \frac{T}{n}, \quad \Delta W_j = Z_j \sqrt{\Delta t}, \quad Z_j \sim \mathcal{N}(0, 1). \]

Os incrementos \(\Delta W_j\) nos dão um processo discreto \(\{W_j^n\}_{j = 0, \ldots, n}\) via

\[ W_0^n = 0, \quad W_j^n = W_{j-1}^n + \Delta W_j. \]

Para cada \(n\) fixo, podemos interpretar \(\{W_j^n\}_j\) como, de fato, uma amostra exata de um processo de Wiener \(\{W_t\}_{t \geq 0}\), com \(W_j^n = W_{t_j}\).

De posse de \(\{W_j^n\}_j\), podemos considerar uma amostra da solução exata

\[ X_{t_j}^n = X_0 e^{(\mu + \sigma^2/2)t_j + \sigma W_j^n}. \]

Usamos o método de Monte-Carlo para resolver a equação um certo número \(M\in \mathbb{N}\) de vezes, obtendo caminhos amostrais \(\{X_j^n(\omega_m)\}_{j = 1}^n\) do método de Euler-Maruyama, assim como amostras de \(M\) caminhos amostrais \(\{X_{t_j}^n\}_{j = 1}^n\) da solução exata.

Ao final, podemos estimar o erro forte via

\[ e_n^{\mathrm{forte}} = \max_{j = 0, \ldots, n} \mathbb{E}_m\left[ \left| X_j^n - X_{t_j}^n\right|\right] = \max_{j=0, \ldots, n}\frac{1}{M}\sum_{m=1}^M \left| X_j^n(\omega_m) - X_{t_j}^n(\omega_m)\right|, \]

onde o valor esperado é tomado em relação às amostras \(\omega_1, \ldots, \omega_M\).

Observe que devemos considerar o erro em relação ao passo \(\Delta t\), mas também temos o erro em relação à amostragem, que pode ser significativo se \(M\) não for suficientemente grande.

Nas simulações abaixo, tomamos alguns valores de \(\mu\), \(\sigma\), \(T\), \(M\) e \(n\) para ilustrar os efeitos desse parâmetros. Fazemos uma regressão linear de mínimos quadrados nos \((\Delta t^n, e_n^{\mathrm{forte}})\) para encontrar \(\ln(C)\) e a ordem de convergência \(p\) tais que

\[ \ln(e_n^{\mathrm{forte}}) \approx \ln(C) + p \ln(\Delta t_n), \]

de modo que

\[ e_n^{\mathrm{forte}} \approx C\Delta t_n^p, \]

correspondendo a uma taxa de convergência da ordem de \(\Delta t^p\).

As amostras do processo de Wiener são calculadas na malha mais fina. Os passos correspondentes são dados a partir dessas amostras, seja na malha mais fina ou nas mais grossas.

Se o ruído for relativamente pequeno, o erro pode não sofrer tanto e acabar sendo de ordem um pouco melhor:

Nesse caso, não é questão de falta de amostras. A ordem permanece mais alta, se aumentarmos o número de amostras:

Mas se dermos mais tempo ao processo, para o ruído ter um efeito mais significativo, obtemos a ordem forte \(1/2\):

Ou, podemos manter um tempo mais curto mas aumentarmos a difusão, observando novamente a ordem forte \(1/2\):

Se por outro lado, diminuirmos a amostragem, enxergamos uma ordem forte diferente, mas isso é por conta do erro amostral, mascarando a taxa real de convergência:

Convergência fraca de ordem 1 no movimento Browniano geométrico

A ordem da convergência fraca também sofre influências do número de amostras e da intensidade de ruído. Mas não vamos ilustrar isso aqui. Exibimos apenas uma simulação, com uma ordem de convergência próxima de \(1\). Observe, também, que, além da ordem mais alta de convergência, o erro em si tem uma magnitude bem menor.

Convergência forte de ordem 1 no processo de Orstein-Uhlenbeck

A equação, nesse caso, é

\[ \mathrm{d}O_t = -\nu O_t \;\mathrm{d}t + \sigma \;\mathrm{d}W_t, \]

cuja solução é

\[ O_t = e^{-\nu t} O_0 + \int_0^t e^{-\nu (t - s)} \;\mathrm{d}W_s. \]

Como o ruído é aditivo, esperamos recuperar a convergência forte de ordem \(1\).