```import numpy as np
from scipy.stats import norm, t

def fit_normal(ret):
mu_norm, sig_norm = norm.fit(ret)
return mu_norm, sig_norm

def fit_t(ret):
parm = t.fit(ret)
nu, mu_t, sig_t = parm
nu = np.round(nu)
return mu_t, sig_t, nu

def polar_returns(ret, k):
"""
Calculate polar return
:param obs: pandas DataFrame
"""
ret= np.mat(ret)
# Find the radius and the angle decomposition on price relative vectors

# Select the 'window' greater values on the observation
index = np.argpartition(radius, -(int(ret.shape[0] * k) + 1))[-(int(ret.shape[0] * k) + 1):]

# Return the radius and the angle for extreme found values

# Pareto Extreme Risk Index
def ERI(R, Z, w):
# alpha
alpha = safe_div((R.shape[0] - 1), np.log(safe_div(R[:-1], R[-1])).sum())
# gamma
gamma = (1 / (Z.shape[0] - 1)) * np.power(np.clip(w * Z[:-1].T, 0.0, np.inf), alpha).sum()
return gamma

# Normal CVaR
def CVaR(mu, sig, alpha=0.01):
return alpha ** -1 * norm.pdf(norm.ppf(alpha)) * sig - mu

# Student T CVaR
def TCVaR(mu, sig, nu, h=1, alpha=0.01):
xanu = t.ppf(alpha, nu)
return -1 / alpha * (1 - nu) ** (-1) * (nu - 2 + xanu ** 2) * t.pdf(xanu, nu) * sig - h * mu

"""
The following risk and adjusted returns metrics belongs to:
http://www.turingfinance.com/computational-investing-with-python-week-one/
"""

import math
import numpy
import numpy.random as nrand

"""
Note - for some of the metrics the absolute value is returns. This is because if the risk (loss) is higher we want to
discount the expected excess return from the portfolio by a higher amount. Therefore risk should be positive.
"""

def vol(returns):
# Return the standard deviation of returns
return numpy.std(returns)

def beta(returns, market):
# Create a matrix of [returns, market]
m = numpy.matrix([returns, market])
# Return the covariance of m divided by the standard deviation of the market returns
return numpy.cov(m)[0][1] / numpy.std(market)

def lpm(returns, threshold, order):
# This method returns a lower partial moment of the returns
# Create an array he same length as returns containing the minimum return threshold
threshold_array = numpy.empty(len(returns))
threshold_array.fill(threshold)
# Calculate the difference between the threshold and the returns
diff = threshold_array - returns
# Set the minimum of each to 0
diff = diff.clip(min=0)
# Return the sum of the different to the power of order
return numpy.sum(diff ** order) / len(returns)

def hpm(returns, threshold, order):
# This method returns a higher partial moment of the returns
# Create an array he same length as returns containing the minimum return threshold
threshold_array = numpy.empty(len(returns))
threshold_array.fill(threshold)
# Calculate the difference between the returns and the threshold
diff = returns - threshold_array
# Set the minimum of each to 0
diff = diff.clip(min=0)
# Return the sum of the different to the power of order
return numpy.sum(diff ** order) / len(returns)

def var(returns, alpha):
# This method calculates the historical simulation var of the returns
sorted_returns = numpy.sort(returns)
# Calculate the index associated with alpha
index = int(alpha * len(sorted_returns))
# VaR should be positive
return abs(sorted_returns[index])

def cvar(returns, alpha):
# This method calculates the condition VaR of the returns
sorted_returns = numpy.sort(returns)
# Calculate the index associated with alpha
index = int(alpha * len(sorted_returns))
# Calculate the total VaR beyond alpha
sum_var = sorted_returns[0]
for i in range(1, index):
sum_var += sorted_returns[i]
# Return the average VaR
# CVaR should be positive
return abs(sum_var / index)

def prices(returns, base):
# Converts returns into prices
s = [base]
for i in range(len(returns)):
s.append(base * (1 + returns[i]))
return numpy.array(s)

def dd(returns, tau):
# Returns the draw-down given time period tau
values = prices(returns, 100)
pos = len(values) - 1
pre = pos - tau
drawdown = float('+inf')
# Find the maximum drawdown given tau
while pre >= 0:
dd_i = (values[pos] / values[pre]) - 1
if dd_i < drawdown:
drawdown = dd_i
pos, pre = pos - 1, pre - 1
# Drawdown should be positive
return abs(drawdown)

def max_dd(returns):
# Returns the maximum draw-down for any tau in (0, T) where T is the length of the return series
max_drawdown = float('-inf')
for i in range(0, len(returns)):
drawdown_i = dd(returns, i)
if drawdown_i > max_drawdown:
max_drawdown = drawdown_i
# Max draw-down should be positive
return abs(max_drawdown)

def average_dd(returns, periods):
# Returns the average maximum drawdown over n periods
drawdowns = []
for i in range(0, len(returns)):
drawdown_i = dd(returns, i)
drawdowns.append(drawdown_i)
drawdowns = sorted(drawdowns)
total_dd = abs(drawdowns[0])
for i in range(1, periods):
total_dd += abs(drawdowns[i])

def average_dd_squared(returns, periods):
# Returns the average maximum drawdown squared over n periods
drawdowns = []
for i in range(0, len(returns)):
drawdown_i = math.pow(dd(returns, i), 2.0)
drawdowns.append(drawdown_i)
drawdowns = sorted(drawdowns)
total_dd = abs(drawdowns[0])
for i in range(1, periods):
total_dd += abs(drawdowns[i])

def treynor_ratio(er, returns, market, rf):
return (er - rf) / beta(returns, market)

def sharpe_ratio(er, returns, rf):
return (er - rf) / vol(returns)

def information_ratio(returns, benchmark):
diff = returns - benchmark
return numpy.mean(diff) / vol(diff)

def modigliani_ratio(er, returns, benchmark, rf):
np_rf = numpy.empty(len(returns))
np_rf.fill(rf)
rdiff = returns - np_rf
bdiff = benchmark - np_rf
return (er - rf) * (vol(rdiff) / vol(bdiff)) + rf

def excess_var(er, returns, rf, alpha):
return (er - rf) / var(returns, alpha)

def conditional_sharpe_ratio(er, returns, rf, alpha):
return (er - rf) / cvar(returns, alpha)

def omega_ratio(er, returns, rf, target=0):
return (er - rf) / lpm(returns, target, 1)

def sortino_ratio(er, returns, rf, target=0):
return (er - rf) / math.sqrt(lpm(returns, target, 2))

def kappa_three_ratio(er, returns, rf, target=0):
return (er - rf) / math.pow(lpm(returns, target, 3), float(1/3))

def gain_loss_ratio(returns, target=0):
return hpm(returns, target, 1) / lpm(returns, target, 1)

def upside_potential_ratio(returns, target=0):
return hpm(returns, target, 1) / math.sqrt(lpm(returns, target, 2))

def calmar_ratio(er, returns, rf):
return (er - rf) / max_dd(returns)

def sterling_ration(er, returns, rf, periods):
return (er - rf) / average_dd(returns, periods)

def burke_ratio(er, returns, rf, periods):
return (er - rf) / math.sqrt(average_dd_squared(returns, periods))

# tests
def test_risk_metrics():
# This is just a testing method
r = nrand.uniform(-1, 1, 50)
m = nrand.uniform(-1, 1, 50)
print("vol =", vol(r))
print("beta =", beta(r, m))
print("hpm(0.0)_1 =", hpm(r, 0.0, 1))
print("lpm(0.0)_1 =", lpm(r, 0.0, 1))
print("VaR(0.05) =", var(r, 0.05))
print("CVaR(0.05) =", cvar(r, 0.05))
print("Drawdown(5) =", dd(r, 5))
print("Max Drawdown =", max_dd(r))

# Returns from the portfolio (r) and market (m)
r = nrand.uniform(-1, 1, 50)
m = nrand.uniform(-1, 1, 50)
# Expected return
e = numpy.mean(r)
# Risk free rate
f = 0.06
# Risk-adjusted return based on Volatility
print("Treynor Ratio =", treynor_ratio(e, r, m, f))
print("Sharpe Ratio =", sharpe_ratio(e, r, f))
print("Information Ratio =", information_ratio(r, m))
# Risk-adjusted return based on Value at Risk
print("Excess VaR =", excess_var(e, r, f, 0.05))
print("Conditional Sharpe Ratio =", conditional_sharpe_ratio(e, r, f, 0.05))
# Risk-adjusted return based on Lower Partial Moments
print("Omega Ratio =", omega_ratio(e, r, f))
print("Sortino Ratio =", sortino_ratio(e, r, f))
print("Kappa 3 Ratio =", kappa_three_ratio(e, r, f))
print("Gain Loss Ratio =", gain_loss_ratio(r))
print("Upside Potential Ratio =", upside_potential_ratio(r))
# Risk-adjusted return based on Drawdown risk
print("Calmar Ratio =", calmar_ratio(e, r, f))
print("Sterling Ratio =", sterling_ration(e, r, f, 5))
print("Burke Ratio =", burke_ratio(e, r, f, 5))

if __name__ == "__main__":
test_risk_metrics()