Forward Monte Carlo Algorithm in C++, increasing size, program breaks

Background Information:

Here is an outline of the algorithm known as forward Monte Carlo for pricing American options which is from the paper, "A Forward Monte Carlo Method for American Options Pricing" by Daniel Wei-Chung Miao and Yung-Hsin Lee.

Question:

My program works correctly when the time steps N = 100 or anything less and when M = 100 or anything less. But when I increase N or M to say 1000 then my program breaks and does not run and I am not sure why.

Here is my code:

#include <iostream>
#include <cmath>
#include <math.h>
#include <limits>
#include <algorithm>
#include <Eigen/Dense>
#include <Eigen/Geometry>
#include <random>
#include <vector>
#include <time.h>

using namespace Eigen;
using namespace std;


void FM(double T, double r, double q, double sigma, double S0, double K, int M, int N);
MatrixXd generateGaussianNoise(int M, int N);        // Generates Normally distributed random numbers
double Black_Scholes(double T, double K, double S0, double r, double sigma);
double phi(long double x);
VectorXd time_vector(double min, double max, int N);
MatrixXd call_payoff(MatrixXd S, double K);

int main(){
    double r = 0.03;        // Riskless interest rate
    double q = 0.0;         // Divident yield
    double sigma = 0.15;    // Volatility of stock
    int T = 1;              // Time (expiry)
    int N = 100;            // Number of time steps
    double K = 100;         // Strike price
    double S0 = 102;        // Initial stock price
    int M = 1000;           // Number of paths              // Current issue

    FM(T,r,q,sigma,S0,K,M,N);

    return 0;
}

MatrixXd generateGaussianNoise(int M, int N){
    MatrixXd Z(N,M);
    random_device rd;
    mt19937 e2(time(0));
    normal_distribution<double> dist(0.0, 1.0);
    for(int i = 0; i < M; i++){
        for(int j = 0; j < N; j++){
            Z(i,j) = dist(e2);
        }
    }
    return Z;
}

double phi(double x){
    return 0.5 * erfc(-x * M_SQRT1_2);
}

double Black_Scholes(double T, double K, double S0, double r, double sigma){
    double d_1 = (log(S0/K) + (r+sigma*sigma/2.)*(T))/(sigma*sqrt(T));
    double d_2 = (log(S0/K) + (r-sigma*sigma/2.)*(T))/(sigma*sqrt(T));
    double C = S0*phi(d_1) - phi(d_2)*K*exp(-r*T);
    return C;
}

VectorXd time_vector(double min, double max, int N){
    VectorXd m(N + 1);
     double delta = (max-min)/N;
     for(int i = 0; i <= N; i++){
             m(i) = min + i*delta;
     }
    return m;
}

MatrixXd call_payoff(MatrixXd S, double K){
    MatrixXd result(S.rows(),S.cols());
    for(int i = 0; i < S.rows(); i++){
        for(int j = 0; j < S.cols(); j++){
            if(S(i,j) - K > 0){
                result(i,j) = S(i,j) - K;
            }else{
                result(i,j) = 0.0;
            }
        }
    }
    return result;
}

void FM(double T, double r, double q, double sigma, double S0, double K, int M, int N){
    MatrixXd Z = generateGaussianNoise(M,N);
    double dt = T/N;
    VectorXd t = time_vector(0.0,T,N);

    // Generate M paths of stock prices
    MatrixXd S(M,N+1);
    for(int i = 0; i < M; i++){
        S(i,0) = S0;
        for(int j = 1; j <= N; j++){
            S(i,j) = S(i,j-1)*exp((double) (r - q - pow(sigma,2.0))*dt + sigma*sqrt(dt)*(double)Z(i,j-1));
        }
    }

    //
    // If path i is "alive" at time index j - 1 < N, generate the price for time index j, denoted as S = S_ij
    // Case for call option:
    // If j = N, the option is expired with value V = exp(-rT)(S-K)^+ and path i is finished
    // If j < N, calculate S_c = f_C(S)
    // If S > S_c, the option is exercised with value V_i = exp(-rT)(S-K)^+ and path i is stopped. Otherwise,
    // the option is held and path continues to live to the next step j+1
    //
    // Case for put option:
    // If j = N, the option is expired with value V = exp(-rT)(K-S)^+ and path i is finished
    // If j < N, calculate S_p = f_p(S)
    // if S < S_p, the option is exercised with value V_i and path i is stopped. Otherwise,
    // the option is held and path continues to live to the next step j+1.

    // Compute S_c parameters and S_p
    double m = 2*r/(pow(sigma,2.0));
    double n = 2*(r-q)/(pow(sigma,2.0));
    VectorXd k(t.size());
    for(int i = 0; i < k.size(); i++){
        k(i) = 1.0 - exp((double) -r*(double)(T - t(i)));                       // Note the t vector (not sure if this is correct)
    }
    VectorXd Q_2(k.size());
    VectorXd Q_1(k.size());
    for(int i = 0; i < Q_2.size(); i++){
        Q_1(i) = (-1*(n-1) + sqrt((double)(n-1)*(n-1) + (double)4*m/(double)(k(i))))/2.0;                       // Q_1 < 0
        Q_2(i) = (-1*(n-1) + sqrt((double)(n-1)*(n-1) + (double)4*m/(double)(k(i))))/2.0;                       // Q_2 > 0
    }
    double d_1 = (log(S0/K) + (r+sigma*sigma/2.)*(T))/(sigma*sqrt(T));
    double C_e = Black_Scholes(T, K, S0, r, sigma);                 // C_e(S) is the European call option price calculated by Black-Scholes
    double Delta = exp(-q*T)*phi(d_1);
    MatrixXd V(M,N+1);
    VectorXd S_c(Q_2.size());
    MatrixXd call_fun = call_payoff(S,K);
    for(int j = 0; j < N + 1; j++){
        for(int i = 0; i < M; i++){
            if(j == N){
                V(i,j) = exp(-r*T)*call_fun(i,j); //////////////
                //cout << "The option is expired with value " << V(i) << " and path " << i << " is finished" << endl;
            }
            else if(j < N){
                S_c(j) = Q_2(j)*(C_e + K)/(Q_2(j) - (1 - Delta));
            }
            else if (S(i,j) > S_c(j)){
                V(i,j) = exp(-r*T)*call_fun(i,j); ///////////////
                //cout << "The option is expired with value " << V(i) << " and path " << i << " is finished" << endl;
            }
        }
    }
    double Value = 0.0;
    for(int i = 0; i < V.rows(); i++){
        for(int j = 0; j < V.cols(); j++){
            Value += V(i,j);
        }
    }
    Value = 1.0/M * Value;
    cout << C_e << endl;
    cout << endl;
    cout << Value << endl;









}

I am pretty new with C++ so I am not sure how to debug my program when this sort of problem arises. This has happened to me to another algorithm I wrote but when I restarted my computer then it worked fine. Any suggestions are greatly appreciated.

From what the user Incomputable asked I believe this is crashing because of low memory, here are the specifications of my computer:

Update:

Taking the advice from user Daniel Jour, I changed the FM function to void. Following the same sample where I set M = 1000 and leave N = 100 then I get this crash message:

    This application has requested the Runtime to terminate it in an unusual way.
Please contact the application's support team for more information.
Assertion failed!

Program: C:\Users\Morgan Weiss\workspace\Forward_Monte_Carlo\Debug\Forward_Monte_Carlo.exe
File: c:\mingw\include\c++\6.2.0\eigen\src/Core/DenseCoeffsBase.h, Line 365

Expression: row >= 0 && row < rows() && col >= 0 && col < cols()

Update 2:

I set N = 1000 and M = 1000 and it ran just fine with no issue, so I am not sure why if I set N not equal to M the program will crash... Any ideas?

Solution

Looking onto your code the following seems to be quite strange:

MatrixXd call_payoff(MatrixXd S, double K){
    MatrixXd result(S.rows(),S.cols()); <-- result size is exact the same as input size
    ....
}

Then:

VectorXd S_c(Q_2.size()); <---- Vector (one of the dimensions is 1)
MatrixXd call_fun = call_payoff(S,K); <--- Matrix 1xN (or Nx1, I am not sure)

And then:

for(int j = 0; j < N + 1; j++){
    for(int i = 0; i < M; i++){
       ...
       V(i,j) = exp(-r*T)*call_fun(i,j); <---- i and j may be significantly bigger than 1
       ...
    }
 }