# Quantifying Uncertainty: Evaluating Trading Algorithms using Probabilistic Programming

## Quantopian Inc.

• PhD from Brown University -- research on computational neuroscience and machine learning using Bayesian modeling.
• Quantitative Researcher at Quantopian Inc: Building the world's first algorithmic trading platform in the web browser.

# Outline

1. Motivation / Problem: Identify stellar quants with limited data
2. The naive approach: point estimates
3. Quick primer on Random Variables and Bayesian statistics
4. Quantifying uncertainty
5. Including prior knowledge
6. Hierarchical modeling

# The naive approach

Use one of the many risk metrics that people have devised, like the Sharpe ratio:

$$\text{Sharpe} = \frac{\text{returns}}{\text{volatility}}$$

In [5]:
sns.distplot(hier_df.groupby('live_algo_id').last().sharpe)

Out[5]:
<matplotlib.axes._subplots.AxesSubplot at 0x7fa45347a5d0>

In [3]:
for algo_id, df_algo in df.groupby(level=0):
df_algo.reset_index().sharpe.plot(lw=1)
plt.xlabel('days running')
plt.ylabel('Sharpe ratio')

Out[3]:
<matplotlib.text.Text at 0x7f409313eb10>


# 2. Short primer on random variables

• Represents our beliefs about an unknown state.
• Probability distribution assigns a probability to each possible state.
• Not a single number (e.g. most likely state).

## You already know what a variable is...

In [8]:
coin = 0 # 0 for tails
coin = 1 # 1 for heads


## A random variable assigns all possible values a certain probability

In []:
coin = {0: 50%,
1: 50%}


### Alternatively:

coin ~ Bernoulli(p=0.5)

• coin is a random variable
• Bernoulli is a probability distribution
• ~ reads as "is distributed as"

## This was discrete (binary), what about the continuous case?

returns ~ Normal($\mu$, $\sigma^2$)

In [30]:
from scipy import stats
sns.distplot(df.benchmark, kde=False, fit=stats.norm)
plt.xlabel('returns')


# How to estimate $\mu$ and $\sigma$?

• Naive: point estimate
• Set mu = np.mean(data) and sigma = np.std(data)
• Maximum Likelihood Estimate
• Correct answer as $n \rightarrow \infty$

# Bayesian analysis

• Most of the time $n \neq \infty$...
• Uncertainty about $\mu$ and $\sigma$
• Turn $\mu$ and $\sigma$ into random variables
• How to estimate?

# Bayes Formula!

Use prior knowledge and data to update our beliefs.

In [6]:
interactive(gen_plot, n=(0, 600), bayes=True)


# Probabilistic Programming

• Model unknown causes (e.g. $\mu$) of a phenomenon as random variables.
• Write a programmatic story of how unknown causes result in observable data.
• Use Bayes formula to invert generative model to infer unknown causes.

## Approximating the posterior with MCMC sampling

In [8]:
plot_want_get()


# PyMC3

• Probabilistic Programming framework written in Python.
• Allows for construction of probabilistic models using intuitive syntax.
• Fast: Just-in-time compiled by Theano.
• Extensible: easily incorporates custom MCMC algorithms and unusual probability distributions.
• Authors: John Salvatier, Chris Fonnesbeck, Thomas Wiecki
• Still in alpha (lack of documentation)!

## Lets select a single algorithm to work with.

In [9]:
algo_df.sharpe.plot()

Out[9]:
<matplotlib.axes._subplots.AxesSubplot at 0x7fa452f51290>


## Statistical model

$\mu \sim normal(0, 1)$

$\sigma \sim gamma(0.01, 2)$

$returns \sim normal(\mu, \sigma^2)$

$Sharpe = \frac{\mu}{\sigma}$

## This is what the data looks like

In [39]:
print returns_data

[ 0.00021801  0.00065388  0.00272272 ..., -0.00572433  0.00782789
-0.0161355 ]


In [14]:
import pymc as pm

with pm.Model() as model:
# Priors on Random Variables
mean_return = pm.Normal('mean return', mu=0, sd=1)
volatility = pm.Gamma('volatility', mu=.01, sd=2)

# Model returns as Normal
obs = pm.Normal('returns',
mu=mean_return,
sd=volatility,
observed=returns_data)

sharpe = pm.Deterministic('sharpe ratio',
(mean_return / volatility) * np.sqrt(252) )

In [15]:
with model:
# Instantiate MCMC sampler
step = pm.NUTS()
# Draw 500 samples from the posterior
trace = pm.sample(500, step)

 [-----------------100%-----------------] 500 of 500 complete in 0.4 sec

In [13]:
sns.tsplot(b_sharpe, err_style='unit_points', interpolate=True)
plt.ylabel('Bayesian Sharpe ratio');

In [30]:
print 'P(Sharpe Ratio > 0) = %.2f%%' % \
(np.mean(trace['sharpe ratio'] > 0) * 100)

P(Sharpe Ratio > 0) = 95.40%



## Prior knowledge

In [17]:
bmark_mean = df.bmark_rf.mean()
bmark_sd = df.bmark_rf.std()

with pm.Model() as model:
# Priors on Random Variables
mean_return = pm.Normal('mean return', mu=bmark_mean, sd=.001)
volatility = pm.Gamma('volatility', mu=bmark_sd, sd=.001)

# Model returns as Normal
returns = pm.Normal('returns',
mu=mean_return,
sd=volatility,
observed=returns_data)

sharpe = pm.Deterministic('sharpe ratio',
(mean_return / volatility) * np.sqrt(252) )

In [19]:
sns.tsplot(b_sharpe_info, err_style='unit_points', interpolate=True)

pm.forestplot(trace_unpooled, vars=['sharpe ratio']);