European Option Pricing with Python, Java and C++

Please make sure to read the disclaimer at the bottom of the page before continuing reading.

Plain vanilla call and put european options are one of the simplest financial derivatives existing. A european call (put) option is essentially a contract which by paying a fee gives you the right to buy (sell) a stock at a predetermined price, the strike price, at a future date. This kind of contracts was originally intended as an insurance tool for companies to fix the selling price of their goods and hedge the price risk.

Options have some interesting features, for starters, if you are the buyer of an option the payoff is potentially unlimited and the loss is limited to the option price as you can see by the payoff diagram below

of course the reverse is true for the seller (huge downside and limited upside).

An interesting question is: how to price options? What is the ‘fair’ price to be paid for an option? Two answers are possible. The first one is using the Black and Scholes formula and the second one is using the Monte Carlo approach. I am going to attempt to price a european call option using the Monte Carlo approach with Python, Java, and C++.

Assuming the stock can be simulated as I have explained in this article, we can calculate a huge number of payoffs and then take the average value as the expected payoff. The present value of the expected payoff will be the price we are looking for. Then again, we face another tough choice: returns should be assumed normally distributed or should returns distribution be estimated using some statistical method such as KDE? Better to do both right? I am going to use both in Python while in Java and C++ I am going to use the first approach only.

1) Normally distributed returns, Python.
Let’s first start assuming returns are normally distributed and use the Apple stock. The drift and the volatility are estimated using historical prices from Google for a call option which expires on the 7th of August 2015 (price is quoted on the 5th of August).

Assuming normally distributed returns the Python script is as follows

	"""
	MONTE CARLO PLAIN VANILLA OPTION PRICING
	This script is used to estimate the price of a plain vanilla
	option using the Monte Carlo method and assuming normally
	distributed returns using a parametric normal distribution
	approach.
	Call option quotations are available at:
	http://www.google.com/finance/option_chain?q=NASDAQ%3AAAPL&ei=fNHBVaicDsbtsAHa7K-QDQ
	In this script the following assumptions are made:
	- Returns are normally distributed and can be modeled using
	a normal distribution. The price of the stock at time t+1
	can be assumed to be:
	s1 = s0*drift + s0*stdv*Z
	where: t=1 and Z is a normally distributed random variable (0,1)
	The drift term is estimated by averaging historical returns.
	The stdv term is estimated by calculating the standard deviation
	of historical returns.
	"""
	import numpy as np
	# Optional seed
	np.random.seed(12345678)
	# Stocastic walk
	# This function calculates the stochastic integral after periods
	# and returns the final price.
	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
	# Parameters
	s0 = 114.64 # Actual price
	drift = 0.0016273 # Drift term (daily)
	volatility = 0.088864 # Volatility (daily)
	t_ = 365 # Total periods in a year
	r = 0.033 # Risk free rate (yearly)
	days = 2 # Days until option expiration
	N = 100000 # Number of Monte Carlo trials
	zero_trials = 0 # Number of trials where the option payoff = 0
	k = 100 # Strike price
	avg = 0 # Temporary variable to be assigned to the sum
	# of the simulated payoffs
	# Simulation loop
	for i in range(N):
	temp = stoc_walk(s0,drift,volatility,days)
	if temp > k:
	payoff = temp-k
	payoff = payoff*np.exp(-r/t_*days)
	avg += payoff
	else:
	zero_trials += 1
	# Averaging the payoffs
	price = avg/float(N)
	# Priting the results
	print("MONTE CARLO PLAIN VANILLA CALL OPTION PRICING")
	print("Option price: ",price)
	print("Initial price: ",s0)
	print("Strike price: ",k)
	print("Daily expected drift: ",drift*100,"%")
	print("Daily expected volatility: ",volatility*100,"%")
	print("Total trials: ",N)
	print("Zero trials: ",zero_trials)
	print("Percentage of total trials: ",zero_trials/N*100,"%")
	# Output
	# MONTE CARLO PLAIN VANILLA CALL OPTION PRICING
	# Option price: 16.1412296587
	# Initial price: 114.64
	# Strike price: 100
	# Daily expected drift: 0.16272999999999999 %
	# Daily expected volatility: 8.8864 %
	# Total trials: 100000
	# Zero trials: 14664
	# Percentage of total trials: 14.664 %
	# >>>

view raw option_pricing_normal.py hosted with ❤ by GitHub

So far the results seems ok although Google option price for the same option was about 14.65 but there might be some differences related to the risk free rate I assumed and other factors that I did not account for.

2) Estimated returns distribution using KDE, Python.
In this old post I explained how to use Scipy to roughly estimate the probability distribution of a random variable. By using the same approach and some more code, we obtain the following

We can clearly see that the estimated distribution is very different from the standard Normal, and you should bear in mind that the estimated distribution has some probability mass in the tails (not shown in the picture) which is not at all represented by the Normal distribution. Using the auxiliary scripts I used (you can download them here) you can check the high kurtosis and some other interesting statistics.

By using this approach we obtain a price which is closer to the one in Google Finance

	"""
	MONTE CARLO PLAIN VANILLA OPTION PRICING
	This script is used to estimate the price of a plain vanilla
	option using the Monte Carlo method and assuming that returns
	can be simulated using an estimated probability density (KDE estimate)
	Call option quotations are available at:
	http://www.google.com/finance/option_chain?q=NASDAQ%3AAAPL&ei=fNHBVaicDsbtsAHa7K-QDQ
	Risk free* rates can be found here:
	http://www.treasury.gov/resource-center/data-chart-center/interest-rates/Pages/TextView.aspx?data=yield
	In this script the following assumptions are made:
	- Returns are not exactly normally distributed, but they following
	a historical distribution that can be estimated from the historical prices.
	Therefore, the price of the stock at time t+1
	can be assumed to be:
	s1 = s0 + s0*K
	where: t=1 and K is a random variable whose density is described by the
	empirical estimated density of the actual stock returns.
	*I choose these as risk free (even though they are not), the concept of
	risk free may be subjective.
	"""
	import pickle
	import numpy as np
	# Load the estimated probability density
	pickle_in = open("KDE_","rb")
	kde_dist = pickle.load(pickle_in)
	pickle_in.close()
	# Check that the total area under the probability density sums up to 1
	print("Area below the curve from -10 to 10:",kde_dist.integrate_box(-10,10),'\n')
	# Optional seed
	#np.random.seed(12345678)
	# Stocastic walk
	# This function calculates the simulated price after periods
	# and returns the final price.
	def stoc_walk(p,periods):
	w = kde_dist.resample(size=periods)[0] #compare to parametric: w = np.random.normal(0,1,size=periods)
	for i in range(periods):
	p += p*w[i] #compare to parametric: p += dr*p + w[i]*vol*p
	return p
	# Parameters
	s0 = 114.64 # Actual price
	t_ = 365 # Total periods in a year
	r = 0.033 # Risk free rate (yearly)
	days = 2 # Days until option expiration
	N = 100000 # Number of Monte Carlo trials
	zero_trials = 0 # Number of trials where the option payoff = 0
	k=100 # Strike price
	avg = 0 # Temporary variable to be assigned to the sum
	# of the simulated payoffs
	# Simulation loop
	for i in range(N):
	temp = stoc_walk(s0,days)
	if temp > k:
	payoff = temp-k
	payoff = payoff*np.exp(-r/t_*days)
	avg += payoff
	else:
	zero_trials += 1
	# Averaging the payoffs
	price = avg/float(N)
	# Priting the results
	print("MONTE CARLO PLAIN VANILLA CALL OPTION PRICING")
	print("Option price: ",price)
	print("Initial price: ",s0)
	print("Strike price: ",k)
	print("Total trials: ",N)
	print("Zero trials: ",zero_trials)
	print("Percentage of total trials: ",zero_trials/N*100,"%")
	# Output
	# Area below the curve from -10 to 10: 1.0
	# MONTE CARLO PLAIN VANILLA CALL OPTION PRICING
	# Option price: 15.1243063405
	# Initial price: 114.64
	# Strike price: 100
	# Total trials: 100000
	# Zero trials: 563
	# Percentage of total trials: 0.563 %
	# >>>
	# Indeed 15.12 is closer to the real price of 14.65, we improved our model.

view raw Option_pricing_KDE_estimate.py hosted with ❤ by GitHub

A plot of the payoffs could be visually interesting, for instance note that at the bottom of the plot you can see that many payoffs that were 0.

 

I am pretty satisfied with the result, however now it is time to use C++ which is way faster than Python at executing low level tasks such as loops. The Java implementation can be downloaded here: it includes also put options and straddle (call + put) evaluation although the Black and Scholes formula is implemented for the Call option only. Some amendments should be done to the Java class, my original project was to make a simple library for pricing derivatives however since it is a time consuming process It may take some time before it is completed.

This is my C++ implementation:

	// 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)
	*/
	std::default_random_engine generator;
	double stoc_walk(double p,double dr,double vol,int periods)
	{
	double mean = 0.0;
	double stdv = 1.0;
	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.0; //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 = 100000; //Number of Monte Carlo trials
	double zero_trials = 0.0;
	double k = 12.0; //Strike price
	//Initialize random number generator
	//Calculate N payoffs
	double avg = 0.0;
	for(int j=0; j < N; j++)
	{
	double temp = stoc_walk(s0,drift,volatility,days);
	if(temp > k)
	{
	double payoff = temp - k;
	payoff = payoff * exp(-r);
	avg += payoff;
	}
	else
	{
	zero_trials += 1;
	}
	}
	//Average the results
	double price_ = avg/(double)N;
	//Print results
	cout << "MONTE CARLO PLAIN VANILLA CALL OPTION PRICING" << endl;
	cout << "Option price: " << price_ << 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;
	}

view raw option_pricing.cpp hosted with ❤ by GitHub

Disclaimer
This article is for educational purpose only. The author is not responsible for any consequence or loss due to inappropriate use. The article may well contain mistakes and errors. The numbers used might not be accurate. You should never use this article for purposes different from the educational one.