# Extreme cases of normal random numbers and NaN

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

## 1 Answer

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.

See also this related question.

• 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? – dbluesk Jan 21 '17 at 23:54
• 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. – LocalVolatility Jan 22 '17 at 0:01