Python for Finance: Part II: 7 Monte Carlo Simulations
Monte Carlo simulations
We are interested in observing the different possible realizations of a future event.
– Scenario 1
– Scenario 2
– Scenario 3
– Scenario 4
– Scenario 5
– Scenario 6
Historical data => A larger data set with “fictional” data
Current Revenues = Last Year Revenues * (1 + y-o-y growth rate)
– Revenue growth rate – Historical Data or User Intuition
– Revenue volatility – Historical Data or User Intuition
Cogs (Cost of Goods Sold): Modeled as a percentage of revenues
Opex: Modeled as a percentage of revenues
Revenues – Cogs = Gross Profit
Revenues – Opex = Operating Profit
|
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 |
import numpy as np import matplotlib.pyplot as plt rev_m = 170 rev_stdev = 20 iterations = 1000 rev = np.random.normal(rev_m, rev_stdev, iterations) rev plt.figure(figsize=(15,6)) plt.plot(rev) plt.show() COGS = - (rev * np.random.normal(0.6,0.1)) plt.figure(figsize=(15,6)) plt.plot(COGS) plt.show() COGS.mean() // -89.588838273282661 COGS.std() // 10.581370507673284 |
|
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 |
Gross_Profit = rev + COGS Gross_Profit plt.figure(figsize=(15,6)) plt.plot(Gross_Profit) plt.show() max(Gross_Profit) // 111.04813446627591 min(Gross_Profit) // 52.352121850588489 Gross_Profit.mean() // 80.688398780785718 Gross_Profit.std() // 9.5301363387028868 plt.figure(figsize=(10,6)) plt.hist(Gross_Profit, bins =[40,50,60,70,80,90,100,110,120]) plt.show() plt.figure(figsize=(10,6)) plt.hist(Gross_Profit, bins =20) plt.show() |
Asset pricing with Monte Carlo
Price Today = Price Yesterday * er
r: log return of share price between yesterday and today.
In(price today / price yesterday)
e.g. eIn(x) = x
Price Today = Price Yesterday * e
Logarithm Basics
log2(16) = x
2x = 16
x = 4
log100(1) = 0
1000 = 1
log2(2) = 1/3
81/3 = 2
log2(1/8) = -3
2-3 = 1/8
log8(1/2) = -1/3
8-1/3 = 1/81/3 = 1/2
Brownian motion
We can use Brownian motion in order to model r:
– Drift:
=> The direction rates of return have been headed in the past.
In(Current Price / Previous Price)
=> Calculate average, standard deviation and variance of daily returns in the historical period.
Drift = (μ – 1/2σ2)
Drift is the expected daily return of the stock.
– Volatility:
Random variable.
Random variable = σ * Z(Rand(0;1))
Price Today = Price Yesterday * eDrift + Random variable
Repeat the calculation 1,000 times.
|
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 |
import numpy as np import pandas as pd from pandas_datareader import data as wb import matplotlib.pyplot as plt from scipy.stats import norm %matplotlib inline ticker = 'AMZN' data = pd.DataFrame() data[ticker] = wb.DataReader(ticker, data_source = 'google', start='2007-1-1')['Close'] log_returns = np.log(1+data.pct_change()) log_returns.tail() data.plot(figsize=(10,6)) log_returns.plot(figsize=(10,6)) u=log_returns.mean() u // AMZN 0.001228 // dtype: float64 var = log_returns.var() var // AMZN 0.000635 // dtype: float64 drift = u - (0.5*var) drift // AMZN 0.00091 // dtype: float64 stdev = log_returns.std() stdev // AMZN 0.025207 // dtype: float64 |
|
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 |
type(drift) // pandas.core.series.Series type(stdev) // pandas.core.series.Series np.array(drift) // array([ 0.0009102]) // object.values: transfers the object into a numpy array drift.values // array([ 0.0009102]) stdev.values // array([ 0.02520659]) norm.ppf(0.95) // 1.6448536269514722 x = np.random.rand(10,2) x norm.ppf(x) z = norm.ppf(np.random.rand(10,2)) z t_intervals = 1000 iterations = 10 daily_returns = np.exp(drift.values + stdev.values * norm.ppf(np.random.rand(t_intervals, iterations))) daily_returns |
Euler’s Method
Differential equations introduction
yII + 2yI = 3y
fII(x) + 2fI(x) = 3f(x)
Leibniz’s notation
d2y / dx2 + 2 (dy / dx) = 3y
|
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 |
S0 = data.iloc[-1] S0 // AMZN 995.17 // Name: 2017-06-19 00:00:00, dtype: float64 price_list = np.zeros_like(daily_returns) price_list price_list[0] = 50 price_list for t in range(1, t_intervals): price_list[t] = price_list[t - 1]*daily_returns[t] price_list plt.figure(figsize=(10,6)) plt.plot(price_list) |
An Introduction to Derivative Contracts
A derivative is a financial instrument, whose price is derived based on the development of one or more underlying assets.
Originally, derivatives served as a hedging instrument.
– Hedging
– Speculating
– Aribtrageurs
Four Types of derivatives
– Forwards
Two parties agree that one party will sell to the other an underlying asset at a future point of time.
– Futures
Highly standardized forward contracts.
– Swaps
Two parties agree to exchange cash flows based on an underlying asset.
e.g. Interest rate, Stock Price, Bond Price, Commodity
– Options
An option contract enables its owner to buy or sell an underlying asset at a given price.
The Black Scholes formula
– A tool for derivatives pricing.
– calculates the value of an option.
– The holder of the option may decide he wants to buy the stock, but he may also decide he is better off without doing it. This freedom is valuable to every investor. Hence, it has a price.
A Call Option’s Payoff
– Strike Price vs. Share Price
– Share Price > Strike Price –> Exercise
– Strike Price > Share Price –> Don’t Exercise
The Black Scholes Formula
The Black Scholes formula calculates the value of a call by taking the difference between the amount you get if you exercise the option minus the amount you have to pay if you exercise the option.
|
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 |
import numpy as np import pandas as pd from pandas_datareader import data as wb from scipy.stats import norm def d1(S, K, r, stdev, T): return (np.log(S / K) + (r - stdev ** 2/2) * T) / (stdev * np.sqrt(T)) def d2(S, K, r, stdev, T): return (np.log(S / K) + (r - stdev ** 2/2) * T) / (stdev * np.sqrt(T)) norm.cdf(0) // 0.5 norm.cdf(0.25) // 0.5987063256829237 norm.cdf(0.75) // 0.77337264762313174 norm.cdf(9) // 1.0 def BSM(S, K, r, stdev, T): return (S*norm.cdf(d1(S, K, r, stdev, T))) - (K * np.exp(-r * T) * norm.cdf(d2(S, K, r, stdev, T))) ticker = 'AMZN' data = pd.DataFrame() data[ticker] = wb.DataReader(ticker, data_source='google', start='2007-1-1', end='2017-3-21')['Close'] S = data.iloc[-1] S // AMZN 843.2 // Name: 2017-03-21 00:00:00, dtype: float64 log_returns = np.log(1+data.pct_change()) stdev = log_returns.std()*250**0.5 stdev // AMZN 0.402353 // dtype: float64 r = 0.025 K = 110.0 T = 1 d1(S, K, r, stdev, T) // AMZN 4.922989 // dtype: float64 d2(S, K, r, stdev, T) // AMZN 4.922989 // dtype: float64 BSM(S, K, r, stdev, T) // AMZN 735.915596 // Name: 2017-03-21 00:00:00, dtype: float64 |
|
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 |
import numpy as np import pandas as pd from pandas_datareader import data as wb from scipy.stats import norm import matplotlib.pyplot as plt %matplotlib inline ticker = 'AMZN' data = pd.DataFrame() data[ticker] = wb.DataReader(ticker, data_source='google', start='2007-1-1', end='2007-3-21')['Close'] log_returns = np.log(1 + data.pct_change()) r = 0.025 stdev = log_returns.std() * 250 ** 0.5 stdev // AMZN 0.28593 // dtype: float64 type(stdev) // pandas.core.series.Series stdev = stdev.values stdev // array([ 0.28592956]) T=1.0 t_intervals = 250 delta_t = T / t_intervals iterations = 10000 Z = np.random.standard_normal((t_intervals + 1, iterations)) S = np.zeros_like(Z) S0 = data.iloc[-1] S[0] = S0 for t in range(1, t_intervals + 1): S[t]=S[t-1] * np.exp((r - 0.5 * stdev ** 2) * delta_t + stdev * delta_t ** 0.5 * Z[t]) S S.shape // (251, 10000) plt.figure(figsize=(10,6)) plt.plot(S[:,:10]) |
|
1 2 3 4 5 6 7 8 9 10 |
p = np.maximum(S[-1] - 110, 0) p // array([ 0., 0., 0., ..., 0., 0., 0.]) p.shape // (10000,) C = np.exp(-r * T) * np.sum(p) / iterations C // 0.0 |
