This notebook is based on matlab code provided through the course "Monte Carlo Methods in Finance".

https://iversity.org/my/courses/monte-carlo-methods-in-finance/

I attempt to convert all matlab code to a Python version, and try to stay as close as possible as the original implementation. Of course, in many cases I make use of the functionality of Python libraries, such as NumPy, SciPy and pandas.

To ensure you have all these packages (correctly) installed I recommend installing Python + all relevant packages through a package manager, such as Anaconda, Python(X,Y) or Enthought Canopy.

In [2]:

from __future__ import division
import numpy as np
import pandas as pd

Week 1¶

europeanOption -- demo 1¶

In this demonstration the price of a European asset is determined through an expected payoff function. This model assumes stocks are described through a log normal distribution.

The payoff function is passed as an argument to the pricing function 'priceEuropeanOption'.

In [3]:

# quad is a numerical integrator. You pass it a function and integration limits.
from scipy.integrate import quad

# this object 'norm' is be used to calculate the PDF or CDF or other properties of the normal distribution
from scipy.stats import norm 

def priceEuropeanOption(S0,r,T,sigma,payoff):
    def ST(x):
        return S0 * np.exp( (r-.5 * sigma**2) * T + sigma * np.sqrt(T) * x )
    R = 10
    def integrand(x):
        return norm.pdf(x) * payoff(ST(x))
    discountFactor = np.exp(-r*T)
    return discountFactor * quad( integrand , -R, R)[0] 

demo_priceEuropeanOption_asset: Pricing the asset¶

Sanity check: the pricing of a stock using priceEuropeanOption should equal its initial price S0.

In [4]:

def demo_priceEuropeanOption_asset():
    S0 = 100
    sigma = .4
    T = 2
    def payoff(ST):
        return ST
    r = .05
    return priceEuropeanOption(S0,r,T,sigma,payoff)

print demo_priceEuropeanOption_asset()

100.0

demo_priceEuropeanOption_call: Price of a European call option¶

Two demonstrations: the price of a European Call Option (the right but not the obligation to buy the stock at strike K at maturity T) and the pricing of a normal stock.

In [5]:

def demo_priceEuropeanOption_call():
    S0 = 100
    sigma = .4
    K = 90
    T = 2
    def payoff(ST):
        return np.max([ST - K, 0])
    r = .05
    return priceEuropeanOption(S0,r,T,sigma,payoff)


print demo_priceEuropeanOption_call()

30.761890184

investmentPortfolio -- demo 2¶

This demonstration constructs a portfolio of three stocks.

This requires stock data, which we can load from the yahoo website using the pandas library.

In [6]:

from pandas.io.data import DataReader
from datetime import datetime
import matplotlib.pyplot as plt
S = DataReader(["IBM", "GOOG", "SI"],  "yahoo", datetime(2007,7,1), datetime(2013,6,30))['Adj Close']

In [24]:

S.head()

Out[24]:

	GOOG	IBM	SI
Date
2007-07-02	530.38	93.52	117.05
2007-07-03	534.34	94.92	120.98
2007-07-05	541.63	96.23	120.87
2007-07-06	539.40	97.10	120.13
2007-07-09	542.56	97.05	122.16

In [7]:

S.plot();

In [8]:

initialValue = 100
normalizedS = initialValue * S / S.ix[0]
normalizedS.plot();

demo_investmentPortfolio: Time series of portfolio values from the series of asset prices.¶

In [9]:

# define composition of portfolio 
c = [500, 200, 200]
# construct time series of portfolio; we use the list to ensure the order of the stock matches
P = (S[["IBM", "GOOG", "SI"]] * c).sum(axis = 1)
# and normalize
normalizedP = initialValue * P / P.ix[0]

normalizedP.name = "Portfolio"

plt.figure();
normalizedS.plot();
normalizedP.plot();
plt.show();

<matplotlib.figure.Figure at 0xbe26ed0>

In [10]:

w = (S[["IBM", "GOOG", "SI"]] * c) / P.ix[0]
w.plot();

Returns -- demo 3¶

demo_Bachelier_BlackScholes: Compare Bachelier and Black-Scholes models¶

In [11]:

import pandas as pd
from pandas.io.data import DataReader
from datetime import datetime
import matplotlib.pyplot as plt
S = DataReader(["IBM", "GOOG", "SI"],  "yahoo", datetime(2007,7,1), datetime(2013,6,30))['Adj Close']

Defining the changes dS and the log changes r. We combine these into a single dataframe

In [12]:

Stocks = S.join(S.diff(), rsuffix='_diff' )
Stocks = Stocks.join(np.log(S).diff(), rsuffix='_log_diff' )
Stocks.head()

Out[12]:

	GOOG	IBM	SI	GOOG_diff	IBM_diff	SI_diff	GOOG_log_diff	IBM_log_diff	SI_log_diff
Date
2007-07-02	530.38	93.52	117.05	NaN	NaN	NaN	NaN	NaN	NaN
2007-07-03	534.34	94.92	120.98	3.96	1.40	3.93	0.007439	0.014859	0.033024
2007-07-05	541.63	96.23	120.87	7.29	1.31	-0.11	0.013551	0.013707	-0.000910
2007-07-06	539.40	97.10	120.13	-2.23	0.87	-0.74	-0.004126	0.009000	-0.006141
2007-07-09	542.56	97.05	122.16	3.16	-0.05	2.03	0.005841	-0.000515	0.016757

Make a list of the columns which correspond to returns.

In [13]:

returns = [colnames for colnames in Stocks.columns if colnames.endswith('_diff') ]

In [15]:

nBins = 100
R = 10
Bins = np.linspace(-R, R, nBins)

# This function computes and plots the histogram, and returns axis objects for each figure
# The axis objects can be used to set things like plotranges and titles.
axes = (Stocks[returns] / (Stocks[returns].quantile(.75) - Stocks[returns].quantile(.25) )).hist(\
                bins = Bins, figsize=(10, S.columns.size * 5), facecolor='w', edgecolor='k', layout = (S.columns.size, 2));
for ax in axes.flatten(): 
    ax.set_xlim(left = -R, right = R) # enforces axis limits
    ax.set_title(ax.get_title() + '_normalized')

demo_logReturns: Properties of log returns in a time series¶

Here we determine and plot the log returns of the stock data.

In [16]:

# As the name suggests, we use the shift operator to shift the index of S. This computes S[1:] / S[:-1]. 
Stocks = S.join( np.log(S / S.shift(1)), rsuffix='_logreturns' )
Stocks.head()

Out[16]:

	GOOG	IBM	SI	GOOG_logreturns	IBM_logreturns	SI_logreturns
Date
2007-07-02	530.38	93.52	117.05	NaN	NaN	NaN
2007-07-03	534.34	94.92	120.98	0.007439	0.014859	0.033024
2007-07-05	541.63	96.23	120.87	0.013551	0.013707	-0.000910
2007-07-06	539.40	97.10	120.13	-0.004126	0.009000	-0.006141
2007-07-09	542.56	97.05	122.16	0.005841	-0.000515	0.016757

In [17]:

# To fix the order of plotting we construct a list
column_order = []
for i in S.columns: column_order.extend([i, i + '_logreturns'])

In [18]:

# The list column_order selects the columns of Stocks in a particular order. No data is copied in this process.
axes = Stocks[column_order].plot(subplots=True, figsize=(10, S.columns.size * 5))

demo_simpleReturns: Properties of returns in a time series¶

Here we also include the ordinary returns of the stocks, and construct a portfolio.

In [19]:

Stocks = Stocks.join( S / S.shift(1), rsuffix='_returns' )
Stocks.head()

Out[19]:

	GOOG	IBM	SI	GOOG_logreturns	IBM_logreturns	SI_logreturns	GOOG_returns	IBM_returns	SI_returns
Date
2007-07-02	530.38	93.52	117.05	NaN	NaN	NaN	NaN	NaN	NaN
2007-07-03	534.34	94.92	120.98	0.007439	0.014859	0.033024	1.007466	1.014970	1.033575
2007-07-05	541.63	96.23	120.87	0.013551	0.013707	-0.000910	1.013643	1.013801	0.999091
2007-07-06	539.40	97.10	120.13	-0.004126	0.009000	-0.006141	0.995883	1.009041	0.993878
2007-07-09	542.56	97.05	122.16	0.005841	-0.000515	0.016757	1.005858	0.999485	1.016898

In [20]:

# To fix the order of plotting we construct a list
column_order = []
for i in S.columns: column_order.extend([i, i + '_returns'])

In [21]:

Stocks[column_order].plot(subplots=True, figsize=(10, S.columns.size * 5));

Next we look at a portfolio P, and plot its return

In [22]:

# portfolio weight
c = [500, 200, 200] 

# add portfolio to Stocks
Stocks['P'] = (S[["IBM", "GOOG", "SI"]] * c).sum(axis = 1)

# and compute the return of the portfolio
Stocks['P_returns'] = 100 * ( Stocks['P'] / Stocks['P'].shift(1) - 1)
Stocks.head()

Out[22]:

	GOOG	IBM	SI	GOOG_logreturns	IBM_logreturns	SI_logreturns	GOOG_returns	IBM_returns	SI_returns	P	P_returns
Date
2007-07-02	530.38	93.52	117.05	NaN	NaN	NaN	NaN	NaN	NaN	176246	NaN
2007-07-03	534.34	94.92	120.98	0.007439	0.014859	0.033024	1.007466	1.014970	1.033575	178524	1.292512
2007-07-05	541.63	96.23	120.87	0.013551	0.013707	-0.000910	1.013643	1.013801	0.999091	180615	1.171271
2007-07-06	539.40	97.10	120.13	-0.004126	0.009000	-0.006141	0.995883	1.009041	0.993878	180456	-0.088033
2007-07-09	542.56	97.05	122.16	0.005841	-0.000515	0.016757	1.005858	0.999485	1.016898	181469	0.561356

In [23]:

# A plot of how our portfolio is performing
Stocks[['P','P_returns']].plot(subplots=True, figsize=(10,  5));