Let

  be a vector of all asset weights
and

 be a vector of all asset volatility (square root of the diagnal of V, the covariance matrix)

Then define diversification ratio as

It is basically the weighted average of volatility devided by the portfolio volatility

An interesting hypothesis/observation: if the assest volatility is positively correlated with asset expected excess return, then maximizing diversification ratio is related to maximizing ex ante Sharpe ratio of the portfolio.

For example, Using CAPM

and beta is

Obviously  is not a constant because correlation of stock return to market return varies across assets/stocks. But Beta is indeed positively correlated with asset volatility, thus the expected excess return is positively correlated with asset volatility under the CAPM.

Python implementation using scipy optimize

def calc_diversification_ratio(w, V):
    # average weighted vol
    w_vol = np.dot(np.sqrt(np.diag(V)), w.T)
    # portfolio vol
    port_vol = np.sqrt(calculate_portfolio_var(w, V))
    diversification_ratio = w_vol/port_vol
    # return negative for minimization problem (maximize = minimize -)
    return -diversification_ratio


#####################################################################
#               PORTFOLIO Optimization functions                    #
#####################################################################

def max_div_port(w0, V, bnd=None, long_only=True):
    # w0: initial weight
    # V: covariance matrix
    # bnd: individual position limit
    # long only: long only constraint
    cons = ({'type': 'eq', 'fun': total_weight_constraint},)
    if long_only: # add in long only constraint
        cons = cons + ({'type': 'ineq', 'fun':  long_only_constraint},)
    res = minimize(calc_diversification_ratio, w0, bounds=bnd, args=V, method='SLSQP', constraints=cons)
    return res

In this example. We going to apply principal component analysis on equity return covariance matrix to construct principal component portfolios because they have some interesting characteristics. The weights are based on the eigenvectors.

For example, the first PC is a factor that captures maximal amount of covariation = linear combination of assets that has highest possible variance. Second PC factor is second most variable portfolio that is orthogonal to first PC.

The principal components are just statistical factors that contain no real meaning. However, many findings suggest that the first principal component portfolio is related the market portfolio (at least close to) or the market risk premium from the CAPM model

Step by Step Analysis:

1.Download Dow Jones Constituents daily stock prices from Yahoo Finance or Bloomberg

2.Calculate returns and the covariance matrix of the returns

3. Apply PCA on the covariance matrix and the explained variance % looks like following:

Notice the explained variance ratio quickly decayed, so let’s just look at the first 3 PC

4. Re-scale the PC’s eigenvectors to sum up to 1 so they can be used as portfolio weights

5. Construct a portfolio using the weights derived from PC

6. Then we calculate the cumulative return of the portfolio and compare the PC portfolio return to the market return

The first PC portfolio return looks identical to the market portfolio

The weights of the equities in our portfolio are:

The second and third PC portfolios return look like this ( I have no idea what they are: value & momentum maybe?):

To verify whether PC1 portfolio still looks like the market for other time horizon, I repeated the same process using weekly and monthly data, the results are pretty robust

Weekly

Monthly

import numpy as np
import matplotlib.pyplot as plt
import pandas as pd
from sklearn.decomposition import PCA as sklearnPCA

# read in data
stocks_ = pd.read_csv('C:\indu_dly.csv',index_col =0)
indu_index = pd.read_csv('C:\indu_index_dly.csv',index_col =0) # DJOW total return index

# get cumulative total return index normalized at 1
indu_index_ret = indu_index.pct_change()
indu_ret_idx = indu_index_ret[1:]+1
indu_ret_idx= indu_ret_idx.cumprod()
indu_ret_idx.columns =['dow']

# calculate stock covariance
stocks_ret = stocks_.pct_change()
stocks_ret = stocks_ret[1:] # skip the first row (NaN)
stocks_cov = stocks_ret.cov()

# USING SKLEARN
sklearn_pca = sklearnPCA(n_components=20) # let's look at the first 20 components
pc = sklearn_pca.fit_transform(stocks_ret)

# plot the variance explained by pcs
plt.bar(range(20),sklearn_pca.explained_variance_ratio_)
plt.title('variance explained by pc')
plt.show()

# check the explained variance reatio
sklearn_pca.explained_variance_ratio_

# get the Principal components
pcs =sklearn_pca.components_

# first component
pc1 = pcs[0,:]
# normalized to 1 
pc_w = np.asmatrix(pc1/sum(pc1)).T

# apply our first componenet as weight of the stocks
pc1_ret = stocks_ret.values*pc_w

# plot the total return index of the first PC portfolio
pc_ret = pd.DataFrame(data =pc1_ret, index= stocks_ret.index)
pc_ret_idx = pc_ret+1
pc_ret_idx= pc_ret_idx.cumprod()
pc_ret_idx.columns =['pc1']

pc_ret_idx['indu'] = indu_index[1:]
pc_ret_idx.plot(subplots=True,title ='PC portfolio vs Market',layout =[1,2])

# plot the weights in the PC
weights_df = pd.DataFrame(data = pc_w*100,index = stocks_.columns)
weights_df.columns=['weights']
weights_df.plot.bar(title='PCA portfolio weights',rot =45,fontsize =8)

This post presents a script implementation of CreditMetrics VaR calculation in python. The code follows the calculations and standards in R ‘CreditMetrics’ Package from CRAN.

CreditMetrics was developed by J.P Morgan in 1997 and is used as a tool for accessing portfolio risk due to changes in debt value caused by changes in credit quality.

It takes a Migration Matrix (usually provided by a credit rating agency)

• The Migration Matrix uses historical transition probabilities between rating classes

To compute the loss distribution, the tool needs the following information

• Probability of Default (PD): risk of losing money
• Exposure At Default (EAD): amount of risk
• Recovery Rate (RR): the fraction of EAD that may be recovered when default happens
• Credit Spread: used to calculate expected market value for securities with different ratings

With simulation we can also get a distribution of loss

• Idea: company return V are assumed to be standard normal distributed V ~ N(0,1) and simulated
• Company value V at time t determines the rating at time t
• Migration thresholds can be computed from the transition probabilities
  • Taking a BB rated issuer from the migration Matrix, we can build a standard normal density function and the thresholds (cut-off) of each rating using normal inverse function

• Migration Thresholds have to be computed for each output ratings
  • e.g. Default threshold is the quintile for the probability of default
  • Z_D = N^-1(P_D) , N^-1 is the inverse function of the standard normal distribution and P_D is the probability of default (given by the Migration Matrix)
  • Threshold for migration to CCC is
  • Z_CCC = N^-1(P_D + P_ccc)

• Credit Spread is the risk premium demanded by the market
• 
• Under risk neutral probabilities
• 
• And if Default happens, the expected value of V is
• 
• Rearrange the above equation we get
• 
• t = 1 for one year migration matrix

To simulate correlated random variables, use Cholesky Decomposition

(Note, this is a simplified approach of CreditMetrics and it’s assuming a constant default correlation, in practice this can be estimated using asset correlations; the VaR calculated from this framework does not capture market risk – interest rate movement)

Example:
Assume, there are 3 assets with exposure of [1590899,1007000,188000] in the portfolio
Total exposure:2,785,899
Their ratings are [BBB, AA, B]; In code we use their indices in the transition matrix.
Recover rate = 55%

1% CreditVaR = 140,388.88
1% Shortfall  = 310,713.24

 from __future__ import division
import pandas as pd
import numpy as np
from scipy.stats import norm

from numpy.linalg import cholesky

# User inputs
TransMat = np.matrix("""
90.81, 8.33, 0.68, 0.06, 0.08, 0.02, 0.01, 0.01;
0.70, 90.65, 7.79, 0.64, 0.06, 0.13, 0.02, 0.01;
0.09, 2.27, 91.05, 5.52, 0.74, 0.26, 0.01, 0.06;
0.02, 0.33, 5.95, 85.93, 5.30, 1.17, 1.12, 0.18;
0.03, 0.14, 0.67, 7.73, 80.53, 8.84, 1.00, 1.06;
0.01, 0.11, 0.24, 0.43, 6.48, 83.46, 4.07, 5.20;
0.21, 0, 0.22, 1.30, 2.38, 11.24, 64.86, 19.79""")/100

N = 3 # number of asset
Nsim = 5000 # num sim for CVaR
r = 0 # risk free rate
exposure = np.matrix([1590899,1007000,188000]).T
recover = 0.55
LGD = 1- recover
idx = [3, 1, 5] # rating index start at 0
t= 1

# correlation matrix
rho = 0.4
sigma = rho*np.ones((N,N))
sigma = sigma -np.diag(np.diag(sigma)) + np.eye(N)

# compute the cut off for each credit rating
Z=np.cumsum(np.flipud(TransMat.T),0)
Z[Z>=1] = 1-1/1e12;
Z[Z<=0] = 0+1/1e12;

CutOffs=norm.ppf(Z,0,1) # compute cut offes by inverting normal distribution

# credit spread implied by transmat
PD_t = TransMat[:,-1] # default probability at t
credit_spread = -np.log(1-LGD*PD_t)/1

# simulate jointly normals with sigma as vcov matrix
# use cholesky decomposition

c = cholesky(sigma)
# cut off matrix for each bond based on their ratings
cut = np.matrix(CutOffs[:,idx]).T
# reference value 
EV = np.multiply(exposure, np.exp(-(r+credit_spread[idx])*t))
PORT_V = sum(EV)

# bond state variable for security Value
cp = np.tile(credit_spread.T,[N,1])
state = np.multiply(exposure,np.exp(-(r+cp)*t))
state = np.append(state,np.multiply(exposure,recover),axis=1) #last column is default case
states = np.fliplr(state) # keep in same order as credit cutoff

Loss=np.zeros((N,Nsim)) # initialization of value array for MC

# Monte Carlo Simulation Nsim times
for i in range(0,Nsim):
    YY = np.matrix(np.random.normal(size=3))
    rr = c*YY.T
    rating = rr<cut
    rate_idx = rating.shape[1]-np.sum(rating,1) # index of the rating
    row_idx = xrange(0,N)
    col_idx = np.squeeze(np.asarray(rate_idx))
    V_t = states[row_idx,col_idx] # retrieve the corresponding state value of the exposure
    Loss_t = V_t-EV.T
    Loss[:,i] = Loss_t

Portfolio_MC_Loss = np.sum(Loss,0)
Port_Var = -1*np.percentile(Portfolio_MC_Loss,1)
ES = -1*np.mean(Portfolio_MC_Loss[Portfolio_MC_Loss<-1*Port_Var])

I will go over a very basic example of what risk parity is and how to construct a simple risk parity (equal risk) portfolio  and extend it to a risk budgeting portfolio (target risk allocation)

First define the marginal risk contribution as:

Then, the risk contribution of asset j to the total portfolio is:

Risk Parity portfolio is a portfolio which RC are equal across all individual assets

For example, going back to the previous example and assume we have an equal weight portfolio of four assets with weight = 25% each, it’s RC (in total portfolio risk %) is:

To compute the weight of a risk parity portfolio, we could use the optimise function from python

Let the sum of squared error of a portfolio assets RC  be:ThenIn python

from __future__ import division
import numpy as np
from matplotlib import pyplot as plt
from numpy.linalg import inv,pinv
from scipy.optimize import minimize

 # risk budgeting optimization
def calculate_portfolio_var(w,V):
    # function that calculates portfolio risk
    w = np.matrix(w)
    return (w*V*w.T)[0,0]

def calculate_risk_contribution(w,V):
    # function that calculates asset contribution to total risk
    w = np.matrix(w)
    sigma = np.sqrt(calculate_portfolio_var(w,V))
    # Marginal Risk Contribution
    MRC = V*w.T
    # Risk Contribution
    RC = np.multiply(MRC,w.T)/sigma
    return RC

def risk_budget_objective(x,pars):
    # calculate portfolio risk
    V = pars[0]# covariance table
    x_t = pars[1] # risk target in percent of portfolio risk
    sig_p =  np.sqrt(calculate_portfolio_var(x,V)) # portfolio sigma
    risk_target = np.asmatrix(np.multiply(sig_p,x_t))
    asset_RC = calculate_risk_contribution(x,V)
    J = sum(np.square(asset_RC-risk_target.T))[0,0] # sum of squared error
    return J

def total_weight_constraint(x):
    return np.sum(x)-1.0

def long_only_constraint(x):
    return x

x_t = [0.25, 0.25, 0.25, 0.25] # your risk budget percent of total portfolio risk (equal risk)
cons = ({'type': 'eq', 'fun': total_weight_constraint},
{'type': 'ineq', 'fun': long_only_constraint})
res= minimize(risk_budget_objective, w0, args=[V,x_t], method='SLSQP',constraints=cons, options={'disp': True})
w_rb = np.asmatrix(res.x)

Then the weights in the risk parity portfolio are

[ 0.19543974,  0.21521557,  0.16260951,  0.42673519]

and the risk contributions to the total portfolio are

The problem is actually set up to  calculate any risk budget allocation since it actually uses you can set up your target risk allocation for each asset then minimise the objective function of squared error

e.g.

x_t = [0.3, 0.3, 0.1, 0.3] # your risk budget percent of total portfolio risk

The weights are
[ 0.22837243, 0.25116466, 0.08875776, 0.43170515]

individual asset RC percentages are:

Suppose stock price S satisfies the following SDE:
 we define

The following is part of the Ito’s Lemma (3)
Let y be a stochastic variable that follows the process
 Then
 Then we plug the following variables into the Ito’s Lemma
 After rearrangements, and substituting the PDE of Y=log(S) with respect to S and t, we get
 and Since

This is the analytic solution to the SDE, We can simulate asset price using the above equation

from __future__ import division
from random import gauss
from math import exp, sqrt
from matplotlib import pyplot as plt

def generate_asset_price(S,v,r,T):
    return S * exp((r - 0.5 * v**2) * T + v * sqrt(T) * gauss(0,1.0))

# USER INPUT
S0 = 28.65 # underlying price
v = 0.2076 # vol of 20.76%
mu = 0.02 # mu
dt = 1/252 # 1 day
T = 2 # period end
n = int(T/dt) # number of steps

S_path=[]
S=S0 # starting price
for i in xrange(1,n+1):
    S_t = generate_asset_price(S,v,mu,dt)
    S= S_t
    S_path.append(S_t)

plt.plot(S_path)

Run the same code but add bounds for the optimiser

# positive weight portfolio
bnd=[(0, 1),(0, 1),(0, 1),(0, 1)] # only positive weights
res2= minimize(calculate_portfolio_var, w0, args=V, bounds = bnd, method='SLSQP',constraints=cons)
w_g_long_only = res2.x

The weights of unconstrained MVP

The weights of long only MVP, Asset 1 has a slightly positive weight

And it’s easy to limit weights on individual holdings
E.g. an asset cannot be more than  50% of the portfolio

bnd=[(0, .5),(0, .5),(0, .5),(0, .5)]

from __future__ import division
import numpy as np
from matplotlib import pyplot as plt
from numpy.linalg import inv,pinv
from scipy.optimize import minimize

# USER INPUT
V = np.matrix('123 37.5 70 30; 37.5 122 72 13.5; 70 72 321 -32; 30 13.5 -32 52')/100  # covariance
R = np.matrix('14; 12; 15; 7')/100
rf = 3/100

w0= [0.25,0.25,0.25,0.25]
# min var optimization
def calculate_portfolio_var(w,V):
    w = np.matrix(w)
    return (w*V*w.T)[0,0]

# unconstrained portfolio (only sum(w) = 1 )
cons = ({'type': 'eq', 'fun': lambda x:  np.sum(x)-1.0})
res= minimize(calculate_portfolio_var, w0, args=V, method='SLSQP',constraints=cons)
w_g = res.x
mu_g = w_g*R
var_g = np.dot(w_g*V,w_g)

The variables and calculation are from APPENDIX OF “A CRITIQUE OF THE ASSET PRICING THEORY’S TESTS” ROLL (1977)

Let return and covariance matrix be
 
and let I be a 4 by 1 vector of all ones

From definition A.9 in ROLL(1977)

The Efficient Frontier:

The Global Minimum Variance Portfolio properties are defined as follows:

The Tangent Portfolio
  then we can use Corollary 1 and
 with , and the weight of the tangent portfolio  is

from __future__ import division
from matplotlib import pyplot as plt
from numpy.linalg import inv,pinv

""" USING THE MATRIX NOTATION DEFINED IN
APPENDIX OF "A CRITIQUE OF THE ASSET PRICING THEORY'S TESTS" ROLL (1977) 
"""

# USER INPUT
V = np.matrix('123 37.5 70 30; 37.5 122 72 13.5; 70 72 321 -32; 30 13.5 -32 52')/100  # covariance
R = np.matrix('14; 12; 15; 7')/100 # return
rf = 3/100 # risk free

# define the 
I = np.ones((len(R),1))
SD = np.sqrt(np.diag(V))

# Variable for Efficient frontier (SEE (A.9) ROLL 1977 page 160)
C = I.T*pinv(V)*I 
B = R.T*pinv(V)*I
A = R.T*pinv(V)*R
D = A*C-B**2

########################################################
#EFFICIENT FRONTIER
########################################################
# Efficient Frontier Range for Return
mu = np.arange(-max(R),max(R)*5,max(R)/100);
 
# Plot Efficient Frontier
minvar = (A-2*B*mu+(C*mu**2))/D;
minstd = np.sqrt(minvar)[0];
minstd = np.squeeze(np.asarray(minstd))
plt.plot(minstd,mu,SD,R,'*')
plt.title('Efficient Frontier with Individual Securities',fontsize=20)
plt.ylabel('Expected Return (%)', fontsize=14)
plt.xlabel('Standard Deviation (%)', fontsize=14)
plt.show()


########################################################
#MVP
########################################################
# Mean and Variance of Global Minimum Variance Portfolio
mu_g = B/C
var_g = 1/C
std_g = np.sqrt(var_g)

# Minimum Variance Portfolio Weights
w_g = (pinv(V)*I)/C

# Plot Efficient Frontier
minvar = (A-2*B*mu+(C*mu**2))/D;
minstd = np.sqrt(minvar)[0];
minstd = np.squeeze(np.asarray(minstd))
plt.plot(minstd,mu,SD,R,'*',std_g,mu_g,'ro')
plt.title('Efficient Frontier with Individual Securities',fontsize=20)
plt.ylabel('Expected Return (%)', fontsize=14)
plt.xlabel('Standard Deviation (%)', fontsize=14)
plt.show()


########################################################
#TANGENT PORTFOLIO
########################################################
# Expected Return of Tangency Portfolio
mu_tan = (B*rf-A)/(C*rf-B);
 
# Variance and Standard Deviation of Tangency Portfolio
vartan = (A-2*rf*B + (rf**2*C))/((B-C*rf)**2);
stdtan = np.sqrt(vartan);
 
# Weights for Tangency Portfolio
w_tan = (pinv(V)*(R - rf*I))/(B-C*rf)

# Tangency Line
m_tan = mu[mu >= rf];
minvar_rf = (m_tan-rf)**2/(A-2*rf*B+C*rf**2);
minstd_rf = np.sqrt(minvar_rf);
minstd_rf = np.squeeze(np.asarray(minstd_rf))

# Plot with tangency portfolio
plt.plot(minstd,mu,SD,R,'*',minstd_rf,m_tan,'r',std_g,mu_g,'bo',stdtan,mu_tan,'go')
plt.title('Efficient Frontier with Individual Securities',fontsize=18)
plt.ylabel('Expected Return (%)', fontsize=12)
plt.xlabel('Standard Deviation (%)', fontsize=12)
plt.text(0.5,rf,'rf',fontsize=12);
plt.text(0.5+std_g,mu_g,'Global Minimum Variance Portfolio',fontsize=12);
plt.text(0.5+stdtan,mu_tan,'Tangency Portfolio',fontsize=12);
plt.show()