2
\$\begingroup\$

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;
}
mdfst13
22.4k6 gold badges34 silver badges70 bronze badges
asked Mar 23, 2017 at 14:21
\$\endgroup\$
11
  • \$\begingroup\$ 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. \$\endgroup\$ Commented Mar 23, 2017 at 14:29
  • \$\begingroup\$ @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 \$\endgroup\$ Commented Mar 23, 2017 at 14:32
  • \$\begingroup\$ you can run the program and look at the task manager. It is impossible to fill up all the memory faster than human reaction. \$\endgroup\$ Commented Mar 23, 2017 at 14:37
  • 2
    \$\begingroup\$ The part about removing the Eigen library is basically asking for code to be written, which is off-topic. \$\endgroup\$ Commented Oct 12, 2017 at 6:29
  • 2
    \$\begingroup\$ @Incomputable Yes the code works now \$\endgroup\$ Commented Oct 12, 2017 at 16:15

0

Know someone who can answer? Share a link to this question via email, Twitter, or Facebook.

Your Answer

Draft saved
Draft discarded

Sign up or log in

Sign up using Google
Sign up using Email and Password

Post as a guest

Required, but never shown

Post as a guest

Required, but never shown

By clicking "Post Your Answer", you agree to our terms of service and acknowledge you have read our privacy policy.

Start asking to get answers

Find the answer to your question by asking.

Ask question

Explore related questions

See similar questions with these tags.