The Black-Scholes-Merton model is a mathematical model applicable in financial markets containing derivative investment instruments for the theoretical estimation of European-style options' price. According to such model, the index level at maturity of a call option can be modeled as it follows with z normally distributed
$$S_T(z)=S_0e^{(r-0.5\sigma^2)T+\sigma\sqrt{T}z}$$where
Let's assume the following values
In order to calculate the present value of the call option, it is possible to apply a Monte Carlo approach according to the following procedure:
In [2]:
import numpy as np
def bsm(S0,r,sigma,T,K,R = 100000 , seed=500):
np.random.seed(seed)
z = np.random.standard_normal(R)
ST = S0 * np.exp(( r - 0.5 * sigma ** 2) * T + sigma * np.sqrt(T) * z)
hT = np.maximum(ST - K, 0)
C0 = np.exp(-r * T) * np.sum(hT) / R
return C0
In [3]:
import time
tm = time.time()
C0 = bsm(S0=105,r=0.06,sigma=0.22,T=1.0,K=109,R = 100000 , seed=500)
pm = time.time() - tm
print("Value of European Call Option: {0:.4g}".format(C0)+" - time[{0:.4g} secs]".format(pm))
Let's see how much time is necessary for 70,000,000 iterations intead of 100,000 iterations.
In [4]:
tm = time.time()
C0 = bsm(S0=105,r=0.06,sigma=0.22,T=1.0,K=109,R = 70000000 , seed=500)
pm = time.time() - tm
print("Value of European Call Option: {0:.4g}".format(C0)+" - time[{0:.4g} secs]".format(pm))
Let's see how we can speed up the computation with the numexpr package.
In [5]:
import numexpr as ne
def bsm_ne(S0,r,sigma,T,K,R = 70000000 , seed=500):
np.random.seed(seed)
z = np.random.standard_normal(R)
ST = ne.evaluate('S0 * exp(( r - 0.5 * sigma ** 2) * T + sigma * sqrt(T) * z)')
hT = np.maximum(ST - K, 0)
C0 = np.exp(-r * T) * np.sum(hT) / R
return C0
In [6]:
tm = time.time()
C0 = bsm_ne(S0=105,r=0.06,sigma=0.22,T=1.0,K=109,R = 70000000 , seed=500)
pm = time.time() - tm
print("Value of the European Call Option: {0:.4g}".format(C0)+" - time[{0:.4g} secs]".format(pm))
The daily return of a stock is easily computable as it follows $$dr(t)=\frac{P(t)}{P(t-1)}-1$$ Similarly, the cumulative return of a stock is easily computable as it follows $$cr(t)=\frac{P(t)}{P(0)}-1$$ What is it P(t)? There are basically 2 options for this value, i.e.
We take the adjusted close price (see What Hedge Funds Really Do to understand why).
Typically for evaluating the performance of a portfolio key factors to focus on are
We will see how to compute and plot these factors.
Functions from pandas.io.data and pandas.io.ga extract data from various Internet sources into a DataFrame The following sources are supported:
For further info see pandas documentation
In [9]:
import numpy as np
import pandas as pd
import pandas.io.data as web
df_final = web.DataReader(['GOOG','SPY'], data_source='yahoo',
start='1/21/2010', end='4/15/2016')
print(df_final)
print(df_final.shape)
In [10]:
df_final.ix[:,:,'SPY'].head()
Out[10]:
In [11]:
print(type(df_final.ix[:,:,'SPY']))
print("\n>>> null values:"+str(pd.isnull(df_final.ix[:,:,'GOOG']).sum().sum()))
In [12]:
df_final = web.DataReader(['GOOG','SPY'], data_source='yahoo',
start='1/21/1999', end='4/15/2016')
df_final.ix[:,:,'GOOG'].head()
Out[12]:
In [13]:
print(type(df_final.ix[:,:,'GOOG']))
print("\n>>> null values:"+str(pd.isnull(df_final.ix[:,:,'GOOG']).sum().sum()))
There is a couple of observations to be done:
Hence, we can define the following functions.
In [14]:
import matplotlib.pyplot as plt
def get_data(symbols,
add_ref=True,
data_source='yahoo',
price='Adj Close',
start='1/21/2010',
end='4/15/2016'):
"""Read stock data (adjusted close) for given symbols from."""
if add_ref and 'SPY' not in symbols: # add SPY for reference, if absent
symbols.insert(0, 'SPY')
df = web.DataReader(symbols,
data_source=data_source,
start=start,
end=end)
return df[price,:,:]
get_data(symbols=['GOOG','SPY']).tail()
Out[14]:
Also, notice that it is not necessary to perform an initial join with the data range of interest filtering out non trading days, as padas does it for us , i.e.
In [15]:
df_stock = get_data(symbols=['GOOG','SPY'],start='1/21/1999',end='4/15/2016')
print(">> Trading days from pandas:"+str(df_stock.shape[0]))
dates = pd.date_range('1/21/1999', '4/15/2016')
df = pd.DataFrame(index=dates)
print(">> Calendar days:"+str(df.shape[0]))
df = df.join(df_stock)
print(">> After join:"+str(df.shape[0]))
df = df.dropna(subset=["SPY"])
print(">> After removing non trading days:"+str(df.shape[0]))
In [16]:
ax = get_data(symbols=['GOOG','SPY','IBM','GLD'],start='1/21/1999', end='4/15/2016').plot(title="Stock Data", fontsize=9)
ax.set_xlabel("Date")
ax.set_ylabel("Price")
plt.show()
In [17]:
def fill_missing_values(df_data):
"""Fill missing values in data frame, in place."""
df_data.fillna(method='ffill',inplace=True)
df_data.fillna(method='backfill',inplace=True)
return df_data
ax = fill_missing_values(get_data(symbols=['GOOG','SPY','IBM','GLD'],
start='1/21/1999',
end='4/15/2016')).plot(title="Stock Data", fontsize=9)
ax.set_xlabel("Date")
ax.set_ylabel("Price")
plt.show()
In [18]:
def normalize_data(df):
return df/df.ix[0,:]
ax = normalize_data(
fill_missing_values(
get_data(
symbols=['GOOG','SPY','IBM','GLD'],
start='1/21/1999',
end='4/15/2016'))).plot(title="Stock Data", fontsize=9)
ax.set_xlabel("Date")
ax.set_ylabel("Normalized price")
plt.show()
In [19]:
df = fill_missing_values(
get_data(
symbols=['GOOG','SPY','IBM','GLD'],
start='4/21/2015',
end='7/15/2016'))
# 1. Computing rolling mean using a 20-day window
rm_df = pd.DataFrame.rolling(df, window=20).mean()
ax = rm_df.plot(title="Rolling Mean")
ax.set_xlabel("Date")
ax.set_ylabel("Price")
plt.show()
# 2. Computing rolling standard deviation using a 20-day window
rstd_df = pd.DataFrame.rolling(df, window=20).std()
ax = rstd_df.plot(title="Rolling Standard Deviation")
ax.set_xlabel("Date")
ax.set_ylabel("Price")
plt.show()
# 3. Compute upper and lower bands
def get_bollinger_bands(rm, rstd):
"""Return upper and lower Bollinger Bands."""
upper_band, lower_band = rm + 2 * rstd, rm - 2 * rstd
return upper_band, lower_band
upper_band, lower_band = get_bollinger_bands(rm_df, rstd_df)
# Plot raw SPY values, rolling mean and Bollinger Bands
ax = df['SPY'].plot(title="Bollinger Bands",label='SPY')
rm_df['SPY'].plot(label='Rolling mean', ax=ax)
upper_band['SPY'].plot(label='upper band', ax=ax)
lower_band['SPY'].plot(label='lower band', ax=ax)
# Add axis labels and legend
ax.set_xlabel("Date")
ax.set_ylabel("Price")
ax.legend(loc='lower left')
plt.show()
In [20]:
def compute_daily_returns_2(df):
"""Compute and return the daily return values."""
# Note: Returned DataFrame must have the same number of rows
daily_returns = df.copy()
daily_returns[1:] = (df[1:]/df[:-1].values) - 1
daily_returns.ix[0,:] = 0
return daily_returns
def compute_daily_returns(df):
"""Compute and return the daily return values."""
# Note: Returned DataFrame must have the same number of rows
daily_returns = (df / df.shift(1)) - 1
daily_returns.ix[0,:] = 0
return daily_returns
pd.util.testing.assert_frame_equal(
compute_daily_returns(
fill_missing_values(
get_data(
symbols=['GOOG','SPY','IBM','GLD'],
start='4/21/2016',
end='7/15/2016'))),
compute_daily_returns_2(
fill_missing_values(
get_data(
symbols=['GOOG','SPY','IBM','GLD'],
start='4/21/2016',
end='7/15/2016'))))
ax = compute_daily_returns(fill_missing_values(
get_data(
symbols=['GOOG','SPY','IBM','GLD'],
start='4/21/2016',
end='7/15/2016'))).plot(title="Daily returns")
ax.set_xlabel("Date")
ax.set_ylabel("Daily return")
plt.show()
In [106]:
df = compute_daily_returns(fill_missing_values(get_data(
symbols=['SPY'],
start='4/21/2000',
end='7/15/2016')))
plt.hist(df['SPY'],bins=30,color='c',label=['Daily return'])
plt.axvline(df['SPY'].mean(), color='b', linestyle='dashed', linewidth=2 , label='Mean')
plt.axvline(-df['SPY'].std(), color='r', linestyle='dashed', linewidth=2 , label='Std')
plt.axvline(df['SPY'].std(), color='r', linestyle='dashed', linewidth=2 )
plt.title('SPY daily return distribution')
plt.xlabel('Daily return')
plt.grid(True)
plt.legend()
plt.show()
In [107]:
def cumulative_returns(df):
return df/df.ix[0,:] - 1
ax = cumulative_returns(fill_missing_values(
get_data(
symbols=['GOOG','SPY','IBM','GLD'],
start='4/21/2016',
end='7/15/2016'))).plot(title="Cumulative returns")
ax.set_xlabel("Date")
ax.set_ylabel("Cumulative return")
plt.show()
Sharpe ratio is a way to examine the performance of an investment by adjusting for its risk. The ratio measures the excess return (or risk premium) per unit of deviation in an investment asset or a trading strategy, typically referred to as risk (and is a deviation risk measure), named after William F. Sharpe.
The ex-ante Sharpe ratio is defined as: $$S=\frac{E\{R_p-R_f\}}{Std\{R_p-R_f\}}$$
where $R_p$ is the asset return, $R_f$ is the return on a benchmark asset, such as the risk free rate of return. $E\{R_p-R_f\}$ is the expected value of the excess of the asset return over the benchmark return, and $Std\{R_p-R_f\}$ is the standard deviation of the asset return.
The ex-post Sharpe ratio uses the same equation as the one above but with realized returns of the asset and benchmark rather than expected returns.
Examples of risk free rates are
Using daily returns as $R_p$ the ex-post Sharpe ratio is computable as: $$S=\frac{mean\{R_p^{daily}-R_f^{daily}\}}{Std\{R_p^{daily}\}}$$
as $R_f^{daily}$ is typically costant for several months it is possible to remove from the denominator. Also, $R_f^{daily}$ is typically approximable with 0% but it is possible to compute given $R_f^{yearly}$ as it follows remembering there are 252 trading days in a year
$$R_f^{daily}=\sqrt[252]{1+R_f^{yearly}}$$Sharpe ratio is typically an annual measure. This means if we are using different sample frequencies need to apply the following formula: $$S^{annual}=\sqrt[2]{SPY} \times S$$
where $SPR$ is the number of samples per year considered in computing $S$, e.g.
In [135]:
def sharpe_ratio(df,sample_freq='d',risk_free_rate=0.0):
sr = (df - risk_free_rate).mean() / df.std()
if sample_freq == 'd':
sr = sr * np.sqrt(252)
elif sample_freq == 'w':
sr = sr * np.sqrt(52)
elif sample_freq == 'm':
sr = sr * np.sqrt(12)
else:
raise Exception('unkown sample frequency :'+str(sample_freq))
return sr
In [143]:
# Sharpe ratio
sharpe_ratio(
compute_daily_returns(
fill_missing_values(
get_data(
symbols=['GOOG','SPY','IBM','GLD'],
start='4/21/2015',
end='7/15/2016'))))
Out[143]:
In [162]:
df = fill_missing_values(get_data(symbols=['GOOG','SPY','IBM','GLD'],
start='4/21/2015',
end='7/15/2016'))
In [163]:
# 1. Cumulative return
cumulative_returns(df).ix[-1,:]
Out[163]:
In [164]:
# 2. Average daily return
compute_daily_returns(df).mean()
Out[164]:
In [165]:
# 3. Rsk (Standard deviation of daily return)
compute_daily_returns(df).std()
Out[165]:
In [166]:
# 4. Sharpe ratio
sharpe_ratio(compute_daily_returns(df))
Out[166]: