Analyze Big Financial Data with Python¶
The Python Quant¶
- Founder and MD of The Python Quants GmbH
- Lecturer Mathematical Finance at Saarland University
- Co-organizer of "For Python Quants" conference in New York
- Organizer of "Python for Quant Finance" meetup in London
- Organizer of "Python for Big Data Analytics" meetup in Berlin
- Book (2013) "Derivatives Analytics with Python"
- Book (2014) "Python for Finance", O'Reilly
- Dr.rer.pol in Mathematical Finance
- Graduate in Business Administration
- Martial Arts Practitioner and Fan
See www.hilpisch.com.
IPython¶
IPython Notebook is an interactive development environment (IDE). Among others it is characterized by the following features:
- interactive
Pythoncoding magiccommands- system shell
- markdown capabilities
Latexsupport- export capabilities (HTML, slides, Latex/PDF)
You can integrate different kinds of content into an IPython Notebook.
from IPython.display import HTML, Image, YouTubeVideo, display
Image('http://hilpisch.com/tpq_logo.png', width=400)
HTML('<iframe src=http://quant-platform.com width=800 height=350></iframe>')
vid = YouTubeVideo("99Xcs5nvxYw")
display(vid)
As of today, IPython is the most popular and most powerful (?!) development environment for interactive Python development. It comes in three technological flavours:
IPython ShellIPython QTConsoleIPython Notebook
The only thing that you need in addition is a good text editor (my choice Sublime Text, http://www.sublimetext.com).
And a consistent, easy to use and maintain Python distribution (my choice Anaconda, http://continuum.io/downloads).
IPython as an Interactive Environment¶
You can do calculations, define Python objects and do file operations.
3 + 4 ** 0.5
l = [1, _, 'Some Text']
l
import pickle
with open(path + 'save.pkl', 'w') as f:
pickle.dump(l, f)
ll /flash/data
You can also easily do simulations and plot data.
import matplotlib.pyplot as plt
%matplotlib inline
from random import gauss
data = [gauss(0, 1) for _ in xrange(10000)]
plt.hist(data, bins=35)
plt.grid(True)
plt.xlabel('value')
plt.ylabel('frequency')
There are many useful magic commands that IPython adds to the mix.
%loadpy http://matplotlib.org/mpl_examples/pie_and_polar_charts/polar_scatter_demo.py
"""
Demo of scatter plot on a polar axis.
Size increases radially in this example and color increases with angle (just to
verify the symbols are being scattered correctly).
"""
import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline
N = 150
r = 2 * np.random.rand(N)
theta = 2 * np.pi * np.random.rand(N)
area = 200 * r**2 * np.random.rand(N)
colors = theta
ax = plt.subplot(111, polar=True)
c = plt.scatter(theta, r, c=colors, s=area, cmap=plt.cm.hsv)
c.set_alpha(0.75)
plt.show()
NumPy – The Numerical Workhorse¶
NumPy represents the basic building of the whole scientific stack in Python. The major class is the ndarray class which models a multi-dimensional data container for array-like data (vectors, matrices, cubes, tables).
import numpy as np
data = np.random.standard_normal((50, 20))
data[:, 1]
# indexing, slicing
NumPy provides powerful methods and so-called universal functions (all implemented in C for performance).
data.mean()
# powerful methods ...
np.std(data)
# ... and universal functions (ufuncs)
All other major libraries, like matplotlib interact seamlessly with NumPy arrays.
plt.plot(data.cumsum(axis=0))
plt.grid()
pandas – For Time Series Data and More¶
One major strength of the pandas library is the management and analysis of time series data.
import pandas as pd
index = pd.date_range(start='2014/7/1', periods=len(data), freq='W')
# define DatetimeIndex object
index
df = pd.DataFrame(data, index=index)
# define DataFrame object
pandas brings lots of convenience to data analysis.
df.iloc[:, :3].describe()
Typical ufuncs of NumPy can in general be applied to DataFrame objects. pandas also provides a convenience wrapper for the matplotlib library.
df.cumsum().plot(legend=False)
PyTables for Hardware Bound IO¶
PyTables is a Python wrapper for the high performing HDF5 file/database format (http://www.hdfgroup.org)
import tables as tb
We work with a NumPy array containing 1 GB of random data.
n = 12500
m = 10000
%time data = np.random.standard_normal((n, m))
data.nbytes
# a giga byte worth of data
PyTables reaches in general reading and writing speeds that are only bound by the performance of the hardware.
h5 = tb.open_file(path + 'data.h5', 'w')
%time h5.create_array('/', 'data', data)
h5.get_filesize()
h5.close()
Reading the data is generally pretty fast.
h5 = tb.open_file(path + 'data.h5', 'r')
%time data = h5.root.data.read()
data.flags.owndata
h5.close()
Application – Optimal Mean-Variance Portfolios¶
Let us implement a Markowitz portfolio optimization procedure with the help of NumPy, pandas & PyTables.
symbols = ['AAPL', 'GOOG', 'MSFT', 'DB', 'GLD']
# the symbols for the portfolio optimization
noa = len(symbols)
pandas provides functions to retrieve stock price information from eg Yahoo! or Google.
import pandas.io.data as web
data = pd.DataFrame()
for sym in symbols:
data[sym] = web.DataReader(sym, data_source='google')['Close']
data.columns = symbols
pandas is also well integrated with PyTables for performant DataFrame storage and retrieval.
h5 = pd.HDFStore(path + 'stocks.h5', 'w')
%time h5['stocks'] = data
h5.close()
A graphical inspection of the normalized stock price data.
(data / data.ix[0] * 100).plot(figsize=(8, 5))
For the portfolio analysis, return data is needed.
rets = np.log(data / data.shift(1))
rets.mean() * 252
With pandas, the covariance matrix is only one method call away.
rets.cov() * 252
The return of a portfolio is given as
\[ \begin{eqnarray*} \mu_p &=& \mathbf{E}\left(\sum_I w_i r_i \right) \\ &=& \sum_I w_i \mathbf{E}\left( r_i \right) \\ &=& \sum_I w_i \mu_i \\ &=& w^T \mu \end{eqnarray*} \]
Here, \(w_i\) and \(r_i\) are the weight and the (expected, historical) return of asset \(i\), respectively.
The covariance matrix is defined as follows.
\[ \begin{eqnarray*} \Sigma = \begin{bmatrix} \sigma_{1}^2 \ \sigma_{12} \ \dots \ \sigma_{1I} \\ \sigma_{21} \ \sigma_{2}^2 \ \dots \ \sigma_{2I} \\ \vdots \ \vdots \ \ddots \ \vdots \\ \sigma_{I1} \ \sigma_{I2} \ \dots \ \sigma_{I}^2 \end{bmatrix} \end{eqnarray*} \]
Here, it holds \(\sigma_{ij} = \sigma_{ji} = E\left(\left(r_i-\mu_i\right)\left(r_j-\mu_j\right)\right)\).
Given the covariance matrix, the portfolio variance is then given below.
\[ \begin{eqnarray*} \sigma_p^2 &=& \mathbf{E}\left( (r - \mu )^2 \right) \\ &=& \sum_{i \in I}\sum_{j \in I} w_i w_j \sigma_{ij} \\ &=& w^T \Sigma w \end{eqnarray*} \]
The Python implementations are even a bit more compact than the mathematical definitions ...
weights = np.random.random(noa)
weights /= np.sum(weights)
# getting random weights and normalization
weights
np.sum(weights * rets.mean() ) * 252
# expected portfolio return
np.dot(weights.T, np.dot(rets.cov() * 252, weights))
# expected portfolio variance
np.sqrt(np.dot(weights.T, np.dot(rets.cov() * 252, weights)))
# expected portfolio standard deviation/volatility
We can easily simulate possible mean-variance/volatility combinations.
%%time
prets = []
pvols = []
for p in range (2500):
weights = np.random.random(noa)
weights /= np.sum(weights)
prets.append(np.sum(rets.mean() * weights) * 252)
pvols.append(np.sqrt(np.dot(weights.T,
np.dot(rets.cov() * 252, weights))))
prets = np.array(prets)
pvols = np.array(pvols)
Let us inspect the results graphically.
plt.figure(figsize=(8, 4))
plt.scatter(pvols, prets, c=prets / pvols, marker='o')
plt.grid(True)
plt.xlabel('expected volatility')
plt.ylabel('expected return')
plt.colorbar(label='Sharpe ratio')
Let us derive the efficient frontier. First, some helper functions.
def statistics(weights):
pret = np.sum(weights * rets.mean()) * 252
pvol = np.sqrt(np.dot(weights.T, np.dot(rets.cov() * 252, weights)))
return np.array([pret, pvol])
def min_portfolio_vola(weights):
return statistics(weights)[1]
Second, the optimization procedure.
import scipy.optimize as sco
bnds = tuple((0, 1) for _ in xrange(noa))
# bounds for the asset weights
%%time
trets = np.linspace(0.0, 0.21, 50)
tvols = []
for tret in trets:
cons = ({'type': 'eq', 'fun': lambda x: statistics(x)[0] - tret},
{'type': 'eq', 'fun': lambda x: np.sum(x) - 1})
res = sco.minimize(min_portfolio_vola, noa * [1. / noa,], method='SLSQP',
bounds=bnds, constraints=cons)
tvols.append(res['fun'])
tvols = np.array(tvols)
And the efficient frontier visualized.
plt.figure(figsize=(8, 4))
plt.scatter(pvols, prets, c=prets / pvols, marker='o')
plt.scatter(tvols, trets, c=trets / tvols, marker='x')
plt.grid(True)
plt.xlabel('expected volatility')
plt.ylabel('expected return')
plt.colorbar(label='Sharpe ratio')
Final file inspections.
ll /flash/data/
Application – Principal Component Analysis¶
Principal Component Analysis (PCA) is another important statistical tool in finance. Using the same libraries as before and a "bit more", we are going to implement a PCA based on the names in the German DAX index.
from sklearn.decomposition import KernelPCA
First, we need the relevant stock price and index data.
symbols = ['ADS.DE', 'ALV.DE', 'BAS.DE', 'BAYN.DE', 'BEI.DE',
'BMW.DE', 'CBK.DE', 'CON.DE', 'DAI.DE', 'DB1.DE',
'DBK.DE', 'DPW.DE', 'DTE.DE', 'EOAN.DE', 'FME.DE',
'FRE.DE', 'HEI.DE', 'HEN3.DE', 'IFX.DE', 'LHA.DE',
'LIN.DE', 'LXS.DE', 'MRK.DE', 'MUV2.DE', 'RWE.DE',
'SAP.DE', 'SDF.DE', 'SIE.DE', 'TKA.DE', 'VOW3.DE',
'^GDAXI']
%%time
data = pd.DataFrame()
for sym in symbols:
data[sym] = web.DataReader(sym, data_source='yahoo')['Close']
data = np.log(data / data.shift(1))
# calculating the log returns
data = data.dropna()
# getting rid of rows with missing data
We now have all the return data in a single DataFrame object.
data[data.columns[:5]].head()
For convenience, we separate the DAX index data.
dax = pd.DataFrame(data.pop('^GDAXI'))
dax.hist(bins=35)
Let us implement the PCA.
scale_function = lambda x: (x - x.mean()) / x.std()
pca = KernelPCA().fit(data.apply(scale_function))
Let's check the explanatory power of the first components.
get_we = lambda x: x / x.sum()
# weighting function
get_we(pca.lambdas_)[:10]
# explanatory power of first ten components
get_we(pca.lambdas_)[:10].sum()
# power of first five components (>75%)
A PCA index with 1 component replicating the DAX index.
pca = KernelPCA(n_components=1).fit(data.apply(scale_function))
dax['PCA_1'] = pca.transform(-data)
dax.cumsum().apply(scale_function).plot(figsize=(8, 4))
A PCA index with 10 components replicating the DAX index.
pca = KernelPCA(n_components=10).fit(data.apply(scale_function))
pca_components = pca.transform(data)
weights = get_we(pca.lambdas_)
dax['PCA_10'] = np.dot(pca_components, weights)
dax[['^GDAXI', 'PCA_10']].cumsum().apply(scale_function).plot(figsize=(8, 4))
Parallel Process Simulation¶
Monte Carlo simulation is a computationally demanding task that nowadays is implemented generally on a large scale (eg for Value-at-Risk or Credit-Value-Adjustment calculations).
import math
NumPy per se is well suited to implement efficient numerical algorithms, like Monte Carlo simulation.
def simulate_geometric_brownian_motion(p):
M, I = p
# time steps, paths
S0 = 100; r = 0.05; sigma = 0.2; T = 1.0
# model parameters
dt = T / M
paths = np.zeros((M + 1, I))
paths[0] = S0
for t in range(1, M + 1):
paths[t] = paths[t - 1] * np.exp((r - 0.5 * sigma ** 2) * dt +
sigma * math.sqrt(dt) * np.random.standard_normal(I))
return paths
Vectorization leads to compact code, C implementation under the hood to good performance.
A sample simulation with the simulation function.
%time paths = simulate_geometric_brownian_motion((50, 100000))
# example simulation
paths[:, 0]
paths[-1].mean()
# 5% drift p.a.
And a graphical look at selected data.
plt.plot(paths[:, :10])
plt.grid(True)
# some paths visualized
In many financial application areas, not a single process but a larger number of processes has to be simulated. Such a task lends itself pretty well for parallelization.
from time import time
import multiprocessing as mp
I = 10000 # number of paths
M = 50 # number of time steps
t = 100 # number of tasks/simulations
# running with 8 cores
times = []
for w in range(1, 9):
t0 = time()
pool = mp.Pool(processes=w)
result = pool.map(simulate_geometric_brownian_motion, t * [(M, I), ])
times.append(time() - t0)
A visual comparison of the results.
plt.plot(range(1, 9), times)
plt.plot(range(1, 9), times, 'ro')
plt.grid(True)
plt.xlabel('number of threads')
plt.ylabel('time in seconds')
plt.title('%d Monte Carlo simulations' % t)
Parallel Monte Carlo Valuation¶
A major application area for Monte Simulation is option pricing. The multiprocessing approach translates pretty well to this application area.
# European option class for simulation and valuation
class european_option(object):
import math
import numpy as np
def __init__(self, S0, T, r, sigma, h):
self.S0 = S0
self.T = T
self.r = r
self.sigma = sigma
self.h = h
def generate_paths(self, p):
M, I = p
dt = self.T / M
paths = np.zeros((M + 1, I))
paths[0] = self.S0
for t in range(1, M + 1):
paths[t] = paths[t - 1] * np.exp((self.r - 0.5 * self.sigma ** 2) * dt +
self.sigma * math.sqrt(dt) * np.random.standard_normal(I))
return paths
def value_option(self, p, K):
paths = self.generate_paths(p)
payoff = eval(self.h)
return math.exp(-self.r * self.T) * np.sum(payoff) / p[1]
Let us define a standard European call option.
h = 'np.maximum(paths[-1] - K, 0)'
# payoff for standard European call option
o = european_option(100, 1.0, 0.05, 0.2, h)
# the option object
p = 50, 10000
# number of time steps, number of paths
The class allows, as before, the simulation of the geometric Brownian motion as well as the valuation of European options with arbitrary payoff.
o.generate_paths(p)[:, 0]
# sample simulation
o.value_option(p, 100)
# example valuation
Working with an instance of a class, we need to wrap methods in a function.
# wrap method in function
def worker(K):
return o.value_option((50, 25000), K)
times = []
for w in range(1, 9):
t0 = time()
pool = mp.Pool(processes=w)
result = pool.map(worker, range(80, 120))
# option values for strikes K from 80 to 120
times.append(time() - t0)
And the time results again compared graphically.
plt.plot(range(1, 9), times)
plt.plot(range(1, 9), times, 'ro')
plt.grid(True)
plt.xlabel('number of threads')
plt.ylabel('time in seconds')
plt.title('%d Monte Carlo valuations' % len(result))
Processing Large CSV Files¶
We generate a CSV file on disk that is relatively large. We process this file with the help of pandas and PyTables.
First, some imports.
import os
import numpy as np
import pandas as pd
import datetime as dt
Generating an Example CSV File¶
Number of rows to be generated for random data set.
N = int(1e7)
N
Using both random integers as well as floats.
ran_int = np.random.randint(0, 10000, size=(2, N))
ran_flo = np.random.standard_normal((2, N))
Filename for csv file.
csv_name = path + 'data.csv'
Writing the data row by row.
%%time
try:
open(csv_name, 'r')
except:
with open(csv_name, 'wb') as csv_file:
header = 'date,int1,int2,flo1,flo2\n'
csv_file.write(header)
for _ in xrange(10):
# 10 times the original data set
for i in xrange(N):
row = '%s, %i, %i, %f, %f\n' % \
(dt.datetime.now(), ran_int[0, i], ran_int[1, i],
ran_flo[0, i], ran_flo[1, i])
csv_file.write(row)
print os.path.getsize(csv_name)
Delete the original NumPy ndarray objects.
del ran_int
del ran_flo
Reading some rows to check the content.
with open(csv_name, 'rb') as csv_file:
for _ in xrange(5):
print csv_file.readline(),
Reading and Writing with pandas¶
The filename for the pandas HDFStore.
pd_name = path + 'data.h5p'
h5 = pd.HDFStore(pd_name, 'w')
pandas allows to read data from (large) files chunk-wise via a file-iterator.
it = pd.read_csv(csv_name, iterator=True, chunksize=N / 20)
Reading and storing the data chunk-wise (this processes roughly 12 GB of data).
%%time
for i, chunk in enumerate(it):
h5.append('data', chunk)
if i % 20 == 0:
print os.path.getsize(pd_name)
The resulting HDF5 file – with 100,000,000 rows and five columns.
h5
Disk-Based Analytics with pandas¶
The disk-based pandas DataFrame mainly behaves like an in-memory object – but these operations are not memory efficient.
%time h5['data'].describe()
Data selection and plotting works as with regular pandas DataFrame objects – again not really memory efficient.
%%time
%matplotlib inline
h5['data']['flo2'][0:1000].cumsum().plot()
h5.close()
The major reason is that the DataFrame data structure is broken up (e.g. columns) during storage. For analytics it has to be put together in-memory again.
import tables as tb
h5 = tb.open_file(path + 'data.h5p', 'r')
h5
h5.close()
Reading with pandas and Writing with PyTables¶
The PyTables database file.
import tables as tb
tb_name = path + 'data.h5t'
h5 = tb.open_file(tb_name, 'w')
Using a rec array object of NumPy to provide the row description for the PyTables table. To this end, a custom dtype object is needed.
dty = np.dtype([('date', 'S26'), ('int1', '<i8'), ('int2', '<i8'),
('flo1', '<f8'), ('flo2', '<f8')])
# change dtype for date from object to string
Adding compression to the mix (less storage, better backups, better data transfer, etc.).
filters = tb.Filters(complevel=2, complib='blosc')
Again reading and writing chunk-wise, this time appending to a PyTables table object (this processes again roughly 12GB of uncompressed data).
it = pd.read_csv(csv_name, iterator=True, chunksize=N / 20)
%%time
tab = h5.create_table('/', 'data', np.array(it.read().to_records(index=False), dty),
filters=filters)
# initialize table object by using first chunk and adjusted dtype
for chunk in it:
tab.append(chunk.to_records(index=False))
tab.flush()
The resulting table object – 100,000,000 rows and five columns.
h5.get_filesize()
tab
Out-of-Memory Analytics with PyTables¶
Data on disk can be used as if it would be both in-memory and uncompressed. De-compression is done at run-time.
tab[N:N + 3]
# slicing row-wise
tab[N:N + 3]['date']
# access selected data points
Counting of rows is easily accomplished (although here not really needed).
%time len(tab[:]['flo1'])
# length of column (object)
Aggregation operations, like summing up or calculating the mean value, are another application area.
%time tab[:]['flo1'].sum()
# sum over column
%time tab[:]['flo1'].mean()
# mean over column
Typical, SQL-like, conditions and queries can be added.
%time sum([row['flo2'] for row in tab.where('(flo1 > 3) & (int2 < 1000)')])
# sum combined with condition
h5.close()
All operations have been on quite large data sets that might not fit (if uncompressed) into the memory of a single machine.
ll /flash/data/*
Using compression of course reduces the size of the PyTables table object relative to the csv and the pandas HDFStore files.
!rm /flash/data/*.h5*
Conclusions¶
Python and libraries such as NumPy, pandas, PyTables provide useful means and approaches to work with (large, big) financial data. Using these libraries leads to a number of advantages:
- high level abstraction:
Pythonprovides "easy" access to advanced technologies; even non-expert programmers can become productive quite fast; the power of multi-core architectures is easily harnessed; similarly, analyst and developer productivity is relatively high in general - numerical performance: libraries such as
NumPy, pandascircumvent the performance burden of an interpreted language likePythonby delegating performance-critical operations to routines implemented inC - I/O performance: using
PyTableswhich in turn relies on theHDF5data storage standard allows to do IO at a speed only limited by the performance of the available hardware (eg 2/1 GB/s reading/writing speed with new SSDs from Intel)
The Python Quants GmbH – the company Web site
Dr. Yves J. Hilpisch – my personal Web site
Python for Finance – my NEW book (out as Early Release)
Derivatives Analytics with Python – my derivatives analytics book
www.derivatives-analytics-with-python.com
Contact Us