Showing posts with label flexible probabilities. Show all posts
Showing posts with label flexible probabilities. Show all posts

Tuesday, 14 June 2016

mini-Meucci : Applying The Checklist - Steps 6-7

"There are some days when I think I'm going to die from an overdose of satisfaction."
Salvador Dali, Artist (1904-1989)


Today we'll be visiting 2 sites along Via Meucci, Evaluation and Attribution.

Evaluation

We need some way to measure the goodness of the ex-ante portfolio distribution from the scenarios in the Aggregation step, and for this Meucci introduces the concept of a Satisfaction index.

Given the distribution of the ex-ante performance we can compute the satisfaction index (or its opposite, the risk index) associated with the portfolio in a number of ways:

  • expected utility/certainty equivalent: mean-variance, higher moments, prospect theory
  • spectral/distortion: VaR (economic capital), CVaR, Wang
  • non-dimensional ratios: Sharpe, Sortino, Omega and Kappa ratios

Given that our toy example has the goal of creating a low volatility portfolio, we'll choose the (mathematical) opposite of volatility as our measure of satisfaction e.g. if the ex-ante portfolio volatility is \$40,000 then the satisfaction index is -\$40,000.

In this simple example, the investor is only concerned with risk and not return i.e. less volatility is more satisfying than more volatility.

In the Construction step, this satisfaction index will be used to select the optimal portfolio to invest in.

See slide #11 Ex-ante Evaluation (2-stock example) and slide #48 (general case).

For more details, read chapter 5, "Evaluating Allocations", of A. Meucci, 2005, Risk and Asset Allocation, (Springer).

Example Python Code

In [3]:
%matplotlib inline
from pandas_datareader import data
import numpy as np
import pandas as pd
import datetime
import math
import matplotlib.pyplot as plt
import seaborn

# Get Yahoo data on 30 DJIA stocks and a few ETFs
tickers = ['MMM','AXP','AAPL','BA','CAT','CVX','CSCO','KO','DD','XOM','GE','GS',
           'HD','INTC','IBM','JNJ','JPM','MCD','MRK','MSFT','NKE','PFE','PG',
           'TRV','UNH','UTX','VZ','V','WMT','DIS','SPY','DIA','TLT','SHY']
start = datetime.datetime(2008, 4, 1)
end = datetime.datetime(2016, 5, 31)
rawdata = data.DataReader(tickers, 'yahoo', start, end) 
prices = rawdata.to_frame().unstack(level=1)['Adj Close']

# Quest for Invariance (random walk model) and Estimation (historical approach)
risk_drivers = np.log(prices)

# Set estimation interval = investment horizon (tau)
tau = 21 # investment horizon in days
invariants = risk_drivers.diff(tau).drop(risk_drivers.index[0:tau])

# Projection to the Investment Horizon
# Using the historical simulation approach and setting estimation interval = 
# investment horizon, means that projected invariants = invariants
# Recover the projected scenarios for the risk drivers at the tau-day horizon
risk_drivers_prjn = risk_drivers.loc[end,:] + invariants

# Pricing at the Investment Horizon
# Compute the projected $ P&L per unit of each stock for all scenarios
prices_prjn = np.exp(risk_drivers_prjn)
pnl = prices_prjn - prices.loc[end,:]

# Aggregation at the Investment Horizon
# Aggregate the individual stock P&Ls to projected portfolio P&L for all
# scenarios. Assume equally weighted protfolio at beginning of investment period
capital = 1e6
n_asset = 30
asset_tickers = tickers[0:30]
asset_weights = np.ones(n_asset) / n_asset
# initial holdings ie number of shares
h0 = capital * asset_weights / prices.loc[end, asset_tickers]
pnl_portfolio = np.dot(pnl.loc[:, asset_tickers], h0)

# Apply flexible probabilities to portfolio P&L scenarios
n_scenarios = len(pnl_portfolio)

# Time-conditioned flexible probs with exponential decay
half_life = 252 * 2 # half life of 2 years
es_lambda = math.log(2) / half_life
exp_probs = np.exp(-es_lambda * (np.arange(0, n_scenarios)[::-1]))
exp_probs = exp_probs / sum(exp_probs)
# effective number of scenarios
ens_exp_probs = np.exp(sum(-exp_probs * np.log(exp_probs))) 

# Projected Distribution of Portfolio P&L at Horizon  with flexible probs
import rnr_meucci_functions as rnr
mu_port_e, sigma2_port_e = rnr.fp_mean_cov(pnl_portfolio.T, exp_probs)  

# Ex-ante Evaluation
# Evaluate the ex-ante portfolio by some satisfaction index
# For example, assume the investor evaluates allocations based only on
# volatility, as measured by standard deviation, and does not take into 
# account expected returns. In this case, satisfaction is the opposite of
# the projected volatility of the portfolio
satisfaction = - np.sqrt(sigma2_port_e)
print('Ex-ante satisfaction index : {:,.0f} (in $ terms)'\
      .format(-np.sqrt(sigma2_port_e)))
print('Ex-ante satisfaction index : {:,.2%} (in % terms)'\
      .format(-np.sqrt(sigma2_port_e)/capital))
Ex-ante satisfaction index : -39,747 (in $ terms)
Ex-ante satisfaction index : -3.97% (in % terms)

Attribution

There are 2 goals in this step:

  1. Linearly attribute the portfolio ex-ante P&L to risk factors and a residual
  2. Additively attribute the risk/satisfaction index to the risk factors and a residual

In other words, the objective is to analyse the sources of return and risk that drive the portfolio P&L.

Applications of this analysis include identifying hotspots and unintended bets in the portfolio, monitoring exposure to a given risk factor and hedging a given risk.

Attribution of Portfolio Ex-ante P&L

In practice, this involves a linear regression of the projected portfolio returns against some set of risk factor projected returns . The regression coefficients are called the exposures (loadings / beta) to the risk factors, and this is the output of this sub-step.

Attribution of Risk/Satisfaction Index

If the chosen risk/satisfaction index has the nice property of positive homogeneity, then it can be decomposed and allocated to the risk factors. For example, we can decompose portfolio volatility, Value-at-Risk (VaR), and conditional VaR in this manner (but can't do it for the Sharpe ratio).

See slide #12 Ex-ante Attribution (2-stock example) and slide #49 (general case), and The Prayer (former Checklist).

For a more advanced and innovative discussion on Attribution, check out Meucci's Factors on Demand.

Also do note that it is possible to do the above attributions to the individual assets in the portfolio.

Example Python Code

In [5]:
# Ex-ante Attribution
# Linearly attribute the portfolio ex-ante PnL to the S&P500, long bond ETF
# + a residual
# Additively attribute the volatility of the portfolio's PnL to S&P500, 
# long bond ETF + a residual

# Set factors 
factor_tickers = ['SPY', 'TLT', 'Residual']
n_factor = 2
# Calculate linear returns - historical simulation
asset_rets = np.array(prices.pct_change(tau).ix[tau:, asset_tickers]) 
factor_rets = np.array(prices.pct_change(tau).ix[tau:, ['SPY', 'TLT']]) 

# Calculate portfolio standard deviation (in percentage terms)
port_std = np.sqrt(sigma2_port_e) / capital

# Factor attribution exposures and risk contributions (using flexible probs)
beta, vol_contr_Z = rnr.factor_attribution(asset_rets, factor_rets,
                                           asset_weights, exp_probs, n_factor)

print('Ex-ante factor exposure (beta):')
for i, factor in enumerate(factor_tickers):
    print('\t\t{}:\t{:.2f}'.format(factor, beta[i]))
print('')
print('Ex-ante portfolio volatility = {:.2%}'.format(port_std))
print('\tFactor risk contribution:')
for j, factor in enumerate(factor_tickers):
    print('\t\t{}:\t{:.2%}'.format(factor, vol_contr_Z[j]))

# Plot factor risk contribution chart
rnr.plot_waterfall_chart(pd.Series(vol_contr_Z, index=factor_tickers),
                       'Factor Contribution To Portfolio Volatility')
Ex-ante factor exposure (beta):
  SPY: 0.99
  TLT: 0.03
  Residual: 1.00

Ex-ante portfolio volatility = 3.97%
 Factor risk contribution:
  SPY: 3.84%
  TLT: -0.03%
  Residual: 0.17%

Next stop on the tour is Construction...

P.S. If you do sign-up for the Bootcamp, please let the good folks at ARPM know you learnt about it at returnandrisk.com!

References to Attilio Meucci's Work

  1. The Checklist slides
  2. The Prayer (former Checklist)
  3. A. Meucci, 2005, Risk and Asset Allocation, (Springer)
  4. Factors on Demand http://symmys.com/node/164

Download Python Code

Click here for the GitHub repo

Friday, 10 June 2016

mini-Meucci : Applying The Checklist - Steps 3-5

"In the future, instead of striving to be right at a high cost, it will be more appropriate to be flexible and plural at a lower cost. If you cannot accurately predict the future then you must flexibly be prepared to deal with various possible futures."
Edward de Bono, author and thinker extraordinaire (born 1933)


In this third leg of The Checklist tour, we will take 3 more steps, Projection, Pricing and Aggregation.

Quick Recap

Goal: Apply Meucci's The Checklist to create a toy example low volatility equity portfolio using DJIA stocks, with an investment horizon of 21 days.

Risk drivers: log of stock prices (values): $X_t = \ln V_t$, where $X$ is risk driver, $V$ is stock price (value).

Invariants: modelled using a simple random walk: $X_{t+1} = X_t + \epsilon_{t+1}$, where $\epsilon$ is invariant.

Estimation: Historical (simulation) with Flexible Probabilities approach (HFP)

Projection to the Horizon

The goal of this step is to create various possible scenarios (i.e. de Bono's "various possible futures") of the risk drivers at the investment horizon.

"Ultimately we are interested in the value of our positions at the investment horizon. In order to determine the distribution of our positions, we must first determine the distribution of the risk drivers at the investment horizon. This distribution, in turn, is obtained by projecting to the horizon the invariants distribution, obtained in the Estimation Step 2."
The Prayer (former Checklist)

So 2 sub-steps are involved - first project the invariants, then recover the risk drivers.

In our toy example, we can take a 'short-cut' and set the estimation interval equal to the investment horizon i.e. 21 days (in the previous 2 posts the estimation interval was 1 day). In this case, the estimation step also performs projection of the invariants.

To recover the future risk driver, simply apply the right-hand side of our random walk model i.e. risk driver at horizon = risk driver today plus invariant at horizon.

In general though, if you have an estimation interval shorter than your investment horizon (e.g. 5-day estimation interval and 20-day horizon, you'll need to apply the 2-step recipe of 'project then recover' recursively (e.g. 4 times one after the other).

Read The Prayer (former Checklist) and see slide #33 Projection to the Horizon (general case) and slide #8 (2-stock example).

Like to know more? Go to the ARPM Bootcamp!

Python Code Example

In [1]:
%matplotlib inline
from pandas_datareader import data
import numpy as np
import datetime
import math
import matplotlib.pyplot as plt
import seaborn

# Get Yahoo data on 30 DJIA stocks and a few ETFs
tickers = ['MMM','AXP','AAPL','BA','CAT','CVX','CSCO','KO','DD','XOM','GE','GS',
           'HD','INTC','IBM','JNJ','JPM','MCD','MRK','MSFT','NKE','PFE','PG',
           'TRV','UNH','UTX','VZ','V','WMT','DIS','SPY','DIA','TLT','SHY']
start = datetime.datetime(2008, 4, 1)
end = datetime.datetime(2016, 5, 31)
rawdata = data.DataReader(tickers, 'yahoo', start, end) 
prices = rawdata.to_frame().unstack(level=1)['Adj Close']

# Quest for Invariance (random walk model) and Estimation (historical approach)
risk_drivers = np.log(prices)
# Set estimation interval = investment horizon (tau)
tau = 21 # investment horizon in days
invariants = risk_drivers.diff(tau).drop(risk_drivers.index[0:tau])

# Projection to the Investment Horizon
# Using the historical simulation approach and setting estimation interval = 
# investment horizon, means that projected invariants = invariants.
# Recover the projected scenarios for the risk drivers at the tau-day horizon
risk_drivers_prjn = risk_drivers.loc[end,:] + invariants

Pricing at the Horizon

We now have the projected distribution of the risk drivers and the goal of this step is calculate the projected profit and loss (P&L in currency terms not percentage returns) per unit of each asset.

First sub-step is to derive the projected prices from the projected risk drivers. Recall that in our case the risk drivers were log prices, so we can just exponentiate them to get back the prices.

Second sub-step of calculating the $P&L per unit of each asset is simply subtract today's price from the projected price.

Do this for all tickers/scenarios and you get the joint distribution of the ex-ante $P&L's.

Here is the code to do it - a bit different from the formulas in the slides but equivalent...

Python Code Example

In [2]:
# Compute the projected $ P&L per unit of each stock for all scenarios
prices_prjn = np.exp(risk_drivers_prjn)
pnl = prices_prjn - prices.loc[end,:]

Aggregation

In this step we take the ex-ante $P&L's per unit of each stock and multiply by the quantity of each stock held in the portfolio, and add them up to get the portfolio (aggregated) P&L. Do this for all scenarios and we get the distribution of the projected portfolio P&L at the horizon.

Two examples are shown below for an equally weighted portfolio of the 30 DJIA stocks. The first shows the distribution of the portfolio P&L using equal probabilities for the scenarios, while the second uses time-conditioned flexible probabilities...

Python Code Example

In [8]:
# Aggregation at the Investment Horizon
# Aggregate the individual stock P&Ls into projected portfolio P&L for all scenarios
# Assume equally weighted protfolio at beginning of investment period
capital = 1e6
n_asset = 30
asset_tickers = tickers[0:30]
asset_weights = np.ones(n_asset) / n_asset
# initial holdings ie number of shares
h0 = capital * asset_weights / prices.loc[end, asset_tickers]
pnl_portfolio = np.dot(pnl.loc[:, asset_tickers], h0)

# Apply flexible probabilities to portfolio P&L scenarios
n_scenarios = len(pnl_portfolio)

# Equal probs
equal_probs = np.ones(n_scenarios) / n_scenarios

# Time-conditioned flexible probs with exponential decay
half_life = 252 * 2 # half life of 2 years
es_lambda = math.log(2) / half_life
exp_probs = np.exp(-es_lambda * (np.arange(0, n_scenarios)[::-1]))
exp_probs = exp_probs / sum(exp_probs)
# effective number of scenarios
ens_exp_probs = np.exp(sum(-exp_probs * np.log(exp_probs))) 

# Projected Distribution of Portfolio P&L at Horizon  with flexible probabilities
import rnr_meucci_functions as rnr
mu_port, sigma2_port = rnr.fp_mean_cov(pnl_portfolio.T, equal_probs)  
mu_port_e, sigma2_port_e = rnr.fp_mean_cov(pnl_portfolio.T, exp_probs)  

print('Ex-ante portfolio $P&L mean over horizon (equal probs) : {:,.0f}'\
      .format(mu_port))
print('Ex-ante portfolio $P&L volatility over horizon (equal probs) : {:,.0f}'\
      .format(np.sqrt(sigma2_port)))
print('')
print('Ex-ante portfolio $P&L mean over horizon (flex probs) : {:,.0f}'\
      .format(mu_port_e))
print('Ex-ante portfolio $P&L volatility over horizon (flex probs) : {:,.0f}'\
      .format(np.sqrt(sigma2_port_e)))

fig = plt.figure(figsize=(9, 8))
ax = fig.add_subplot(111)
ax.hist(pnl_portfolio, 50, weights=exp_probs) 
ax.set_title('Ex-ante Distribution of Portfolio P&L \
             (flexbile probabilities with exponential decay)') 
plt.show()
Ex-ante portfolio $P&L mean over horizon (equal probs) : 10,431
Ex-ante portfolio $P&L volatility over horizon (equal probs) : 49,829

Ex-ante portfolio $P&L mean over horizon (flex probs) : 10,395
Ex-ante portfolio $P&L volatility over horizon (flex probs) : 39,750

Monday, 6 June 2016

mini-Meucci : Applying The Checklist - Step 2

"Guessing before proving! Need I remind you that it is so that all important discoveries have been made?”
Henri Poincaré, French mathematician (1854-1912)


In this second leg of The Checklist tour, Estimation, we are going to make some educated guesses about the true unknown distribution of the invariants. But first...

Recap and a bit more on Quest for Invariance

"The quest for invariance is necessary for the practitioners to learn about the future by observing the past in a stochastic environment."
The Prayer (former Checklist)

The essence is that the pattern (distribution) of the invariants will repeat.

"Being able to identify the invariants that steer the dynamics of the risk drivers is of crucial importance because it allows us to project the market randomness to the desired investment horizon."
The Prayer (former Checklist)

Also, recall that we chose a simple random walk to model the invariants (same model for all the stocks and ETFs). But if it isn't sufficient and your time series shows for example mean reversion or volatility clustering, then you should use more sophisticated models like ARMA or GARCH (see reference #3 at end of post).

Estimation

Estimation is a yuge! topic - see chapter 4 of A. Meucci, 2005, Risk and Asset Allocation, (Springer) for an introduction.

In a nutshell, the goal is to find numbers (estimators) that describe the distribution of the invariants e.g. location and dispersion can be estimated by the sample mean and sample covariance matrix.

In general there are 3 approaches to Estimation: historical (non-parametric), analytical (parametric) and copula-marginal (mixed).

Historical Approach

If you have enough data relative to your investment horizon (e.g. 10 yrs of daily data and a 21-day horizon), then it's recommended to go the historical non-parametric way, and this is the road we'll take throughout this tour - but with a twist!

"The simplest of all estimators for the invariants distribution is the nonparametric empirical distribution, justified by the law of large numbers, i.e. "i.i.d. history repeats itself". The empirical distribution assigns an equal probability to each of the past observations in the series of the historical scenarios."
The Prayer (former Checklist)

For an example of the analytical approach, see The Checklist slides (slide #7) where the 2-stock example assumes a bivariate normal distribution. For more on copula-marginal, go to the ARPM Bootcamp!

Flexible Probabilities - the twist explained...

I put Flexible Probabilities into the must-see category of attractions in Meucci-land. Using them is a quick and powerful method of enhancing estimation techniques.

The twist is that instead of applying equal weight to all the historical observations/scenarios (as is common practice in calculating the sample mean and covariance), you can apply different weighting schemes. But make sure the weights are normalised to be like probabilities (i.e. range between 0 and 1 and sum to 1). In other words, use weighted estimators e.g. weighted mean and weighted covariance.

Python Code Examples of Flexible Probabilities

The Python code below shows 2 examples of flexible probabilities:

  1. Time-conditioned - apply higher weight to more recent observations, with the weights decaying exponentially
  2. State-conditioned - using the VIX as an example, give weight to past scenarios where the VIX is greater than 20 and no weight otherwise
In [1]:
%matplotlib inline
from pandas_datareader import data
import numpy as np
import datetime
import math
import matplotlib.pyplot as plt
import matplotlib.gridspec as gridspec
import seaborn

# Get Yahoo data on 30 DJIA stocks and a few ETFs
tickers = ['MMM','AXP','AAPL','BA','CAT','CVX','CSCO','KO','DD','XOM','GE','GS',
           'HD','INTC','IBM','JNJ','JPM','MCD','MRK','MSFT','NKE','PFE','PG',
           'TRV','UNH','UTX','VZ','V','WMT','DIS','SPY','DIA','TLT','SHY']
start = datetime.datetime(2005, 12, 31)
end = datetime.datetime(2016, 5, 30)
rawdata = data.DataReader(tickers, 'yahoo', start, end) 
prices = rawdata.to_frame().unstack(level=1)['Adj Close']
risk_drivers = np.log(prices)
invariants = risk_drivers.diff().drop(risk_drivers.index[0])
T = len(invariants)

# Get VIX data
vix = data.DataReader('^VIX', 'yahoo', start, end)['Close']
vix.drop(vix.index[0], inplace=True)

# Equal probs
equal_probs = np.ones(len(vix)) / len(vix)

# Time-conditioned flexible probs with exponential decay
half_life = 252 * 2 # half life of 2 years
es_lambda = math.log(2) / half_life
exp_probs = np.exp(-es_lambda * (np.arange(0, len(vix))[::-1]))
exp_probs = exp_probs / sum(exp_probs)
# effective number of scenarios
ens_exp_probs = math.exp(sum(-exp_probs * np.log(exp_probs))) 

# State-conditioned flexible probs based on VIX > 20
state_probs = np.zeros(len(vix)) / len(vix)
state_cond = np.array(vix > 20)
state_probs[state_cond] = 1 / state_cond.sum()

# Plot charts
fig = plt.figure(figsize=(9, 8))
gs = gridspec.GridSpec(2, 2)
ax = fig.add_subplot(gs[0, 0])
ax2 = fig.add_subplot(gs[0, 1])
ax3 = fig.add_subplot(gs[1, 0])
ax4 = fig.add_subplot(gs[1, 1])
ax.plot(vix.index, equal_probs)
ax.set_title("Equal Probabilities (weights)")
ax2.plot(vix.index, vix)
ax2.set_title("Implied Volatility Index (VIX)")
ax2.axhline(20, color='r')
ax3.plot(vix.index, exp_probs)
ax3.set_title("Time-conditioned Probabilities with Exponential Decay")
ax4.plot(vix.index, state_probs, marker='o', markersize=3, linestyle='None', alpha=0.7)
ax4.set_title("State-conditioned Probabilities (VIX > 20)")
plt.tight_layout()
plt.show()

There are many other ways to create sets of flexible probabilities (e.g. rolling window), and indeed these can also be combined together using Meucci's Entropy-based technique (see reference #4 at end of post).

Example Application of Historical with Flexible Probabilities Approach

Let's now apply the Historical with Flexible Probabilities (HFP) approach to perform a simple 'stress' analysis, and look at the difference between correlations using equal probabilities and state-conditioned probabilities when the VIX > 20 (stressed market condition)

To keep it simple, we'll choose a few tickers only.

In [2]:
# Stress analysis
import rnr_meucci_functions as rnr

tmp_tickers = ['AAPL', 'JPM', 'WMT', 'SPY', 'TLT']

# HFP distribution of invariants using equal probs
mu, sigma2 = rnr.fp_mean_cov(invariants.ix[:,tmp_tickers].T, equal_probs)

# HFP distribution of invariants using state-conditioned probs (VIX > 20)
mu_s, sigma2_s = rnr.fp_mean_cov(invariants.ix[:,tmp_tickers].T, state_probs)

# Calculate correlations
from statsmodels.stats.moment_helpers import cov2corr
corr = cov2corr(sigma2)
corr_s = cov2corr(sigma2_s)

# Plot correlation heatmaps
rnr.plot_2_corr_heatmaps(corr, corr_s, tmp_tickers, 
                     "HFP Correlation Heatmap - equal probs",
                     "HFP Correlation Heatmap - state probs (VIX > 20)")

The main result is that when volatility is elevated (i.e. stressed market condition), correlations within equities rise, but they fall between bonds (TLT) and equities. This quantifies what we expect to see happen, and such scenario driven simulations are easily accomplished with the HFP approach. Of course, there are many more applications of HFP, and in the Aggregation step we'll look at a portfolio risk example.

Next stop on the tour is Projection...


P.S. If you do sign-up for the Bootcamp, please let the good folks at ARPM know you learnt about it at returnandrisk.com!

References to Attilio Meucci's Work

  1. The Checklist slides
  2. The Prayer (former Checklist)
  3. Review of Discrete and Continuous Processes in Finance: Theory and Applications http://symmys.com/node/131
  4. Historical Scenarios with Fully Flexible Probabilities http://symmys.com/node/150

Download Python Code

Click here for the GitHub repo