Using Euler-Maruyama method to solve Ornstein-Uhlenbeck equation (SDE)

Question

I am trying to implement the Euler–Maruyama method and use it to solve the Ornstein–Uhlenbeck process. I am basing my code on the wikipedia page where a python implementation is shown. More generally I'm trying to get to grips with stochastic differential equations in c++.

This is my code

#include <iostream>
#include <fstream>
#include <random>
#include <vector>
#include <cmath>

using namespace std;

int num_sims = 5;
double t_init = 3;
double t_end  = 7;
int N         = 1000;
float dt      = (float)(t_end - t_init) / N;
double y_init = 0;

double c_theta = 0.7;
double c_mu    = 1.5;
double c_sigma = 0.06;

double dW ( double dt ){
    normal_distribution<double> distribution( dt , sqrt(dt) );
    random_device rd;
    default_random_engine generator(rd());

    return distribution(generator);
}

double mu( double y, double t ){
    return c_theta * (c_mu - y);
}

int main()
{
    vector<double> ts;
    for (double t = t_init; t < t_end + dt; t += dt){
        ts.push_back(t);
    }

    vector<double> ys(N+1, 0);
    ys[0] = y_init;

    for (int h = 0; h < num_sims; h+=1){
        for (int i = 1; i < ts.size(); i+=1 ){
            double t = (i - 1) * dt;
            double y = ys[i - 1];
            ys[i] = y + mu(y, t) * dt + c_sigma * dW(dt);
        
            cout << h << '\t' << ys[i] << endl;
        }
    }
    return 0;
}

This code seems to work well when I plot the results in python. My question is if this code can be written in a more efficient way? I feel like I should be using the boost packages on some level, and also maybe structures. I'm quite new to c++, so maybe my translation of the python code is a bit too word for word.

Answer 1

Make constants constexpr

You can tell C++ that some of your variables are constast by marking them const. If the constants can be computed at compile time, you can even mark them constexpr. Doing this helps the compiler generated more optimal code and will result in an error if you accidentily do try to modify a constant. So:

constexpr int num_sims = 5;
constexpr double t_init = 3.0;
...

Mark variables and functions static where appropriate

If a variable or function is not used by any other source files, you should mark them static. This again might help the compiler produce more efficient code, and can prevent conflicts with other variables and functions with identical names in other source files in your project. It looks like everything except main() can be made static.

Optimize random number generation

You are using the right functions to create normal distributed random numbers, but the way you have written it, it will instantiate a new distribution, random device and generator for each call to dW(). You only need to create them once. To do this, move everything except the return statement out of the function (and make them static):

static normal_distribution<double> distribution(dt, sqrt(dt));
static random_device rd;
static default_random_engine generator(rd());

static double dW(double dt) {
    return distribution(generator);
}

Use N consistently

You created a constant for the number of datapoints you are generating. Once you have done that, use it consistently everywhere. Don't start mixing N, ts.size() and other ways to calculate the size of the array, since this might result in slight differences that will cause problems. For example, floating point values have a limited precision, and because of that the following for-loop might actually run more or less than N iterations:

for (double t = t_init; t < t_end + dt; t += dt) {

Use N consistently:

for (int i = 0; i < N; ++i) {
    ts.push_back(t_init + i * dt);
}

...

for (int h = 0; h < num_sims; ++h) {
    for (int i = 1; i < N; ++i) {
        ...

Alternatively, you could also use ts.size() everywhere, as long as you are consistent. In that case, the line to change is:

vector<double> ys(ts.size());

Note that ys doesn't need to be one larger than ts, and you don't need to initialize it with zeroes, as you will overwrite all values anyway.

Use ++ to increment counters

It's customary in C and C++ to use the ++ and -- operators to increment and decrement counters. And since it can give a slight performance benefit in some cases, prefer to use the prefix version to the postfix version if there is no difference to the result. So:

for (int i = ...; ...; ++i)

Avoid using std::endl

You should avoid using std::endl and use '\n' instead. The former is equivalent to the latter, but it will also force a flush of the output, which is often unnecessary and might slow down your application.

Do you need std::vectors at all?

You are storing the time and y-coordinates of the datapoints in the vectors ts and ys, but you never read them back. The only thing you use is the value of y of the previous timestep, but this could easily be held in a single variable:

int main() {
    for (int h = 0; h < num_sims; ++h) {
        double y = y_init;

        for (int i = 0; i < N; ++i) {
            double t = i * dt;
            y += mu(y, t) * dt + c_sigma * dW(dt);
            cout << h << '\t' << y << '\n';
        }
    }
}

Avoid unnecessary complexity

I feel like I should be using the boost packages on some level, and also maybe structures.

This problem absolutely doesn't need Boost packages nor structs. Follow the KISS principle.