Significant price difference while calculating call option price using Monte Carlo approach in C++ and Python?

Question

I am trying to calculate the price of a european call option using the Monte Carlo approach. I coded the algorithm in C++ and Python. As far as I know the implementation is correct and as N (the number of trials) gets bigger, the price should converge to a similar value in both the programs.

My problem is that as N gets bigger, say just from 1000 to 10000 trials the prices converge to two different values. In C++ the price converges towards the value of 3.30 while with Python it converges towards 3.70.

I think that gap of 0.40 is too wide, I should get more similar resutls. Why is this gap so big? What did I do wrong? I cannot seem to find my mistake.

Here is the code I used:

Python

import numpy as np
import matplotlib.pyplot as plt


def stoc_walk(p,dr,vol,periods):
    w = np.random.normal(0,1,size=periods)
    for i in range(periods):
        p += dr*p + w[i]*vol*p
    return p

s0 = 10;
drift = 0.001502
volatility = 0.026
r = 0.02
days = 255
N = 10000
zero_trials = 0

k=12
payoffs = []

for i in range(N):
    temp = stoc_walk(s0,drift,volatility,days)
    if temp > k:
        payoff = temp-k
        payoffs.append(payoff*np.exp(-r))
    else:
        payoffs.append(0)
        zero_trials += 1

payoffs = np.array(payoffs)
avg = payoffs.mean()

print("MONTE CARLO PLAIN VANILLA CALL OPTION PRICING")
print("Option price: ",avg)
print("Initial price: ",s0)
print("Strike price: ",k)
print("Daily expected drift: ",drift)
print("Daily expected volatility: ",volatility)
print("Total trials: ",N)
print("Zero trials: ",zero_trials)
print("Percentage of total trials: ",zero_trials/N)

C++

//Call option Monte Carlo evaluation;

#include <iostream>
#include <random>
#include <math.h>
#include <chrono>

using namespace std;

/*  double stoc_walk: returns simulated price after periods

    p = price at t=t0
    dr = drift
    vol = volatility
    periods (days)
*/
double stoc_walk(double p,double dr,double vol,int periods)
{
    double mean = 0.0;
    double stdv = 1.0;

    /* initialize random seed: */
    int seed = rand() %1000 + 1;
    //unsigned seed = std::chrono::system_clock::now().time_since_epoch().count();
    std::default_random_engine generator(seed);
    std::normal_distribution<double> distribution(mean,stdv);

    for(int i=0; i < periods; i++)
    {
        double w = distribution(generator);
        p += dr*p + w*vol*p;
    }
    return p;
}

int main()
{
    //Initialize variables
    double s0 = 10;             //Initial price
    double drift = 0.001502;    //daily drift
    double volatility = 0.026;  //volatility (daily)
    double r = 0.02;            //Risk free yearly rate
    int days = 255;             //Days
    int N = 10000;              //Number of Monte Carlo trials
    double zero_trials = 0;

    double k = 12;               //Strike price
    int temp = 0;                //Temporary variable
    double payoffs[N];           //Payoff vector
    double payoff = 0;

    srand (time(NULL));         //Initialize random number generator

    //Calculate N payoffs
    for(int j=0; j < N; j++)
    {
        temp = stoc_walk(s0,drift,volatility,days);
        if(temp > k)
        {
            payoff = temp - k;
            payoffs[j] = payoff * exp(-r);
        }
        else
        {
            payoffs[j] = 0;
            zero_trials += 1;
        }
    }

    //Average the results
    double sum_ = 0;
    double avg_ = 0;
    for(int i=0; i<N; i++)
    {
        sum_ += payoffs[i];
    }
    avg_ = sum_/N;

    //Print results
    cout << "MONTE CARLO PLAIN VANILLA CALL OPTION PRICING" << endl;
    cout << "Option price: " << avg_ << endl;
    cout << "Initial price: " << s0 << endl;
    cout << "Strike price: " << k << endl;
    cout << "Daily expected drift: " << drift*100 << "%" << endl;
    cout << "Daily volatility: " << volatility*100 << "%" << endl;
    cout << "Total trials: " << N << endl;
    cout << "Zero trials: " << zero_trials << endl;
    cout << "Percentage of total trials: " << zero_trials/N*100 << "%";

    return 0;
}

Answer 1

There is a bug in your C++ implementation.

int seed = rand() %1000 + 1;

Seeds will be in the range [1...1000], which means your gaussian sampling is HIGHLY correlated

And second bug is that temp variable is declared int

OK, after simple modification, see code below, avg in C++ version is 3.67

C++: http://codepad.org/D5UZql2P

Python: http://codepad.org/TeAYSwkV

Answer 2

Because this is Monte Carlo, and because the implementation of the distribution algorithms is outside your own code, there is every possibility that your problem has to do with random number generation and not specifically with your code as shown.

This means you are very unlikely to get code help on this here.

What you must do is be sure the two implementations work the same for the same data. To do that, you should generate a set of random input variables and save them separately-- then feed that same random stream into both algorithms. It's the only way to control for the random number generation part of Monte Carlo, and prove to yourself that the implementations are indeed logically identical.