Explosive Euler Paths for Superlinear Diffusions and Log-SDE Simulation

# Extreme cases of normal random numbers and NaN


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

## Answer by LocalVolatility (score 3, accepted)

https://quant.stackexchange.com/a/32018

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.