Casos extremos de números aleatorios normales y NaN

Question

Mientras intentaba implementar mi versión del método de Euler para simular una EDE en C++, me encontré con un problema. Ocurre en algunos casos que la trayectoria generada por mi método acaba dando valores que son NaN (no es un número). Esto ocurre en raras ocasiones, pero normalmente aparecen alrededor de 1 o 2 de estos casos si calculo 100 trayectorias.

Creo que el problema surge con los números aleatorios normales $Z_i$ que estoy usando para la iteración $$S_{i+1}=S_i+f(S_i)\cdot \sqrt{dt}\cdot Z_i$$ (donde el SDE es por ejemplo $dS_t=f(S_t)dW_t$ ). He probado a utilizar ambos std::normal_distribution con generador std::mt199737 e implementando mi propio método Box-Muller (generando números uniformes con std::uniform_real_distribution ).

Con ambos enfoques, siempre hay un valor extraño que hace que toda la simulación vaya al infinito y produzca NaN . El caso es que necesito hacer algunos cálculos con los valores finales $S_1$ producido por mis caminos.

¿Debo ignorar el NaN ¿? El número de ellos es muy reducido en comparación con el número total de caminos.

¿Hay gente que tenga este tipo de problemas al realizar este tipo de simulación?

Aquí está mi código:

#include <iostream>
#include <vector>
#include <random>

typedef double ( *funct )( double );

double square( double x ) { return x * x; }
double fourx( double x ) { return 4 * x; }

double generateGaussian();
std::vector< double > createPath( unsigned, double, funct );

int main()
{
  double initial = 1.0;
  unsigned intervals = 3000;
  unsigned paths = 500;

  std::vector< double > path;
  int n = 0;
  for ( auto j = 0; j < paths; ++j )
  {
    path = createPath( intervals, initial, square );
    if ( isnan( path.back() ) )
      ++n;
  }

  std::cout << n << std::endl;
}

double generateGaussian()
{
  static std::mt19937 generator( std::random_device{}() );
  std::normal_distribution< double > dist( 0.0, 1.0 );

  return dist( generator );
}

std::vector< double > createPath( unsigned numberOfIntervals, double initialCondition, funct f )
{
  double sqdt = sqrt( 1.0 / numberOfIntervals );
  std::vector< double > path;
  path.push_back( initialCondition );
  for ( auto i = 0; i < numberOfIntervals; ++i )
  {
    double Z = generateGaussian();
    double Si = path.back();
    double S = Si + f( Si ) * sqdt * Z;

    path.push_back( S );
  }
  return path;
}

Simplemente lo compilo con

g++ -std=c++11 test.cpp

Esto es en un terminal en Mac OS X 10.10.5:

g++ -v
Configured with: --prefix=/Library/Developer/CommandLineTools/usr --with-gxx-include-dir=/usr/include/c++/4.2.1
Apple LLVM version 7.0.2 (clang-700.1.81)
Target: x86_64-apple-darwin14.5.0
Thread model: posix

De vez en cuando me sale una ruta que termina con un NaN . También implementé esto con otras funciones y obtuve ocurrencias similares. Probé, por ejemplo, tomando $f(x)=4x$ (así $dS_t=4S_tdW_t$ ) y a veces obtengo números inusualmente grandes (incluso ignorando el NaN caminos). También en este caso particular de $f(x)=4x$ En cuanto a los resultados, los resultados son totalmente incoherentes de una carrera a otra, lo que me molesta, pero aún no he pensado en el motivo.

Finalmente, si inicializo el generador con 0, obtengo, con este mismo código, 2 NaN 's.

Answer 1

El problema no está en tu código sino en el propio SDE. El proceso

\begin{equation} \mathrm{d}S_t = \sigma S^\beta \mathrm{d}W_t \end{equation}

es una martingala local no negativa. Para $\beta \leq 1$ es una martingala y para $\beta > 1$ es una martingala local estricta. Esto se manifiesta en su simulación por el término de difusión que explota para algunas trayectorias cuando se utiliza la función square . Lleva a S = Inf o S = -Inf y en el siguiente paso y todos los siguientes se obtiene S = NaN .

Al utilizar su función fourx no deberías encontrarte con esto. Una mejor discretización de la SDE

\begin{equation} \mathrm{d}S_t = \sigma S_t \mathrm{d}W_t \end{equation}

es mediante la simulación de los registros. Es decir, se aproxima

\begin{equation} \mathrm{d} \ln S_t = -\frac{1}{2} \sigma^2 \mathrm{d}t + \sigma \mathrm{d}W_t \end{equation}

por

\begin{equation} \ln S_i = \ln S_{i - 1} - \frac{1}{2} \sigma^2 \Delta t + \sigma \sqrt{\Delta t} Z. \end{equation}

De esta forma se elimina la dependencia del estado en el coeficiente de difusión y no se produce ningún error de discretización para cualquier tamaño de paso.

Ver también este pregunta relacionada.