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

1. Generate M paths of stock prices, where each path i = 1,...,M evolves in discrete time with index j = 1,...,N (time interval) $\Delta t = T/N$ as follows:

$$S = S_{i,j} = S_{i,j-1}\exp((r - q - \sigma^2/2)*\Delta t + \sigma \sqrt{\Delta t}Z_{i,j}, \ Z_{i,j}\sim N(0,1)$$

2. If a given path i is alive (option not yet exercised) at time index j - 1 < N, generate the price for the time index j, denoted as $S = S_{ij}$.

1. If it's a call option:

• If j = N, the option is expired with value $V_i = \exp(-r T)(S-K)^{+}$ and path i is finished.

• If j < N, calculate $\hat{S_c} = f_{C}(S)$.

• If S > $\hat{S}_C$, the option is exercised with value $V_i = \exp(-r t)(S - K)^{+}$ and path i is stopped. Otherwise, the option is held and path continues to live to the next step j+1.

2. If it's a put option:

• [omitted; not needed in this program]
3. When all the simulation paths are completed, the American option is valued by averaging the discounted payoff as $V = \frac{1}{M}\sum_{i=1}^{M}V_i$.

I am using C++ and the Eigen library to carry out this algorithm.

# Questions:

The code contains a number of functions needed to price the American option. Is it better to make each function used into a separate class?

I want to make my code run faster and take into consideration storage complexity. I wouldn't mind losing the dependency on the Eigen library in favor of Standard Library types.

# The 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;

double 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 = 100;            // 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;
}

double 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;

}

• Does it crash due to low memory? If not, I'm afraid the post is off topic. My wild guess is that there is something happening with numeric limits, NaNs, +/-infs and other very funny stuff about floating point arithmetic. – Incomputable Mar 23 '17 at 14:29
• @Incomputable That is what I am guessing the problem is. Let me edit that, but I am not compltely sure but the amount of memory I have on this computer sadly is only 3GB – Wolfy Mar 23 '17 at 14:32
• you can run the program and look at the task manager. It is impossible to fill up all the memory faster than human reaction. – Incomputable Mar 23 '17 at 14:37
• The part about removing the Eigen library is basically asking for code to be written, which is off-topic. – 200_success Oct 12 '17 at 6:29
• @Incomputable Yes the code works now – Wolfy Oct 12 '17 at 16:15