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() ) )

  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.

  • 2
    $\begingroup$ 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. $\endgroup$ Jan 21, 2017 at 20:51
  • $\begingroup$ 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. $\endgroup$ Jan 21, 2017 at 21:44

1 Answer 1


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}


\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.

See also this related question.

  • $\begingroup$ Regarding the element $\beta$ you introduced, I tried using $\beta=1/2$, i.e., $dS_t=\sqrt{S_t}dW_t$, and still get lots of nan. What other simulation methods should I look into? $\endgroup$
    – dbluesk
    Jan 21, 2017 at 23:54
  • 1
    $\begingroup$ Now you run into a different problem due to your simulation scheme. You try to take the square root of a number that became negative in your last step. Search for the keywords "Monte Carlo" together with "CEV model" or "SABR model" - you'll find lots of papers on this. $\endgroup$ Jan 22, 2017 at 0:01

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.