Extreme cases of normal random numbers and NaN

Question

While trying to implement my version of Euler's method for simulating a SDE in C++, I came up with a problem. It occurs in some cases that the path generated by my method ends up giving values which are NaN (not-a-number). This happens rarely, but usually around 1 or 2 of these cases appear if I compute 100 paths.

I believe the issue comes up with the normal random numbers $Z_i$ that I'm using for the iteration $$S_{i+1}=S_i+f(S_i)\cdot \sqrt{dt}\cdot Z_i$$ (where the SDE is for example $dS_t=f(S_t)dW_t$). I tried using both std::normal_distribution with generator std::mt199737 and implementing my own Box-Muller method (generating uniform numbers with std::uniform_real_distribution).

With both approaches, there is always the odd value that makes the whole simulation go to infinity and produce NaN. The thing is that I need to do some computations with the final values $S_1$ produced by my paths.

Should I just ignore the NaN ones? There is only a very small number of them compared with the total number of paths.

Do people have these sort of problems when running this type of simulation?

Here is my code:

#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;
}

I simply compile it with

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

This is on a terminal on 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

From time to time I get a path that terminates with a NaN. I also implemented this with other functions and get similar occurrences. I tried for example taking $f(x)=4x$ (so $dS_t=4S_tdW_t$) and sometimes I get unusually huge numbers (even ignoring the NaN paths). Also in this particular case of $f(x)=4x$, I get completely inconsistent results from run to run, which kind of annoys me, but I haven't thought about why yet.

Finally, if I initialize the generator with 0, I get, with this exact same code, 2 NaN's.

I never encountered this in this context. I suggest you post a minimum working example of your code (e.g. the one using std::normal_distribution) that reproduces your problem. — LocalVolatility, Jan 21 '17 at 20:51
I can't reproduce your problem (I get no NaNs). In order to make it fully reproducible please: 1) fix the seed (e.g. at 0) instead of randomizing it, 2) show your full code including main - I don't e.g. know what you used for numberOfIntervals and 3) show your compile command line incl. flags. — LocalVolatility, Jan 21 '17 at 21:44

score 3 · Accepted Answer · edited Apr 13 '17 at 12:46

The problem is not with your code but with the SDE itself. The process

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

is a non-negative local martingale. For $\beta \leq 1$ it is a martingale and for $\beta > 1$ it is a strict local martingale. This manifests in your simulation by the diffusion term exploding for some paths when using the function square. It leads to S = Inf or S = -Inf and on the next step and all following ones you get S = NaN.

When using your function fourx, you shouldn't encounter this. A better discretization of the SDE

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

is by simulating the logs. I.e. you approximate

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

by

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

This way you get rid of the state-dependence in the diffusion coefficient and incur no discretization error for any step size.

Extreme cases of normal random numbers and NaN

1 Answers1