ml-finance-python

edhec_risk_kit_205.py

(23410B)

 import pandas as pd
 import numpy as np
 
 def get_ffme_returns():
     """
     Load the Fama-French Dataset for the returns of the Top and Bottom Deciles by MarketCap
     """
     me_m = pd.read_csv("data/Portfolios_Formed_on_ME_monthly_EW.csv",
                        header=0, index_col=0, na_values=-99.99)
     rets = me_m[['Lo 10', 'Hi 10']]
     rets.columns = ['SmallCap', 'LargeCap']
     rets = rets/100
     rets.index = pd.to_datetime(rets.index, format="%Y%m").to_period('M')
     return rets
 
 def get_fff_returns():
     """
     Load the Fama-French Research Factor Monthly Dataset
     """
     rets = pd.read_csv("data/F-F_Research_Data_Factors_m.csv",
                        header=0, index_col=0, na_values=-99.99)/100
     rets.index = pd.to_datetime(rets.index, format="%Y%m").to_period('M')
     return rets
 
 
 def get_hfi_returns():
     """
     Load and format the EDHEC Hedge Fund Index Returns
     """
     hfi = pd.read_csv("data/edhec-hedgefundindices.csv",
                       header=0, index_col=0, parse_dates=True)
     hfi = hfi/100
     hfi.index = hfi.index.to_period('M')
     return hfi
 
 def get_ind_file(filetype, weighting="vw", n_inds=30):
     """
     Load and format the Ken French Industry Portfolios files
     Variant is a tuple of (weighting, size) where:
         weighting is one of "ew", "vw"
         number of inds is 30 or 49
     """    
     if filetype is "returns":
         name = f"{weighting}_rets" 
         divisor = 100
     elif filetype is "nfirms":
         name = "nfirms"
         divisor = 1
     elif filetype is "size":
         name = "size"
         divisor = 1
     else:
         raise ValueError(f"filetype must be one of: returns, nfirms, size")
     
     ind = pd.read_csv(f"data/ind{n_inds}_m_{name}.csv", header=0, index_col=0, na_values=-99.99)/divisor
     ind.index = pd.to_datetime(ind.index, format="%Y%m").to_period('M')
     ind.columns = ind.columns.str.strip()
     return ind
 
 def get_ind_returns(weighting="vw", n_inds=30):
     """
     Load and format the Ken French Industry Portfolios Monthly Returns
     """
     return get_ind_file("returns", weighting=weighting, n_inds=n_inds)
 
 def get_ind_nfirms(n_inds=30):
     """
     Load and format the Ken French 30 Industry Portfolios Average number of Firms
     """
     return get_ind_file("nfirms", n_inds=n_inds)
 
 def get_ind_size(n_inds=30):
     """
     Load and format the Ken French 30 Industry Portfolios Average size (market cap)
     """
     return get_ind_file("size", n_inds=n_inds)
 
 
 def get_ind_market_caps(n_inds=30, weights=False):
     """
     Load the industry portfolio data and derive the market caps
     """
     ind_nfirms = get_ind_nfirms(n_inds=n_inds)
     ind_size = get_ind_size(n_inds=n_inds)
     ind_mktcap = ind_nfirms * ind_size
     if weights:
         total_mktcap = ind_mktcap.sum(axis=1)
         ind_capweight = ind_mktcap.divide(total_mktcap, axis="rows")
         return ind_capweight
     #else
     return ind_mktcap
 
 def get_total_market_index_returns(n_inds=30):
     """
     Load the 30 industry portfolio data and derive the returns of a capweighted total market index
     """
     ind_capweight = get_ind_market_caps(n_inds=n_inds)
     ind_return = get_ind_returns(weighting="vw", n_inds=n_inds)
     total_market_return = (ind_capweight * ind_return).sum(axis="columns")
     return total_market_return
                          
 def skewness(r):
     """
     Alternative to scipy.stats.skew()
     Computes the skewness of the supplied Series or DataFrame
     Returns a float or a Series
     """
     demeaned_r = r - r.mean()
     # use the population standard deviation, so set dof=0
     sigma_r = r.std(ddof=0)
     exp = (demeaned_r**3).mean()
     return exp/sigma_r**3
 
 
 def kurtosis(r):
     """
     Alternative to scipy.stats.kurtosis()
     Computes the kurtosis of the supplied Series or DataFrame
     Returns a float or a Series
     """
     demeaned_r = r - r.mean()
     # use the population standard deviation, so set dof=0
     sigma_r = r.std(ddof=0)
     exp = (demeaned_r**4).mean()
     return exp/sigma_r**4
 
 
 def compound(r):
     """
     returns the result of compounding the set of returns in r
     """
     return np.expm1(np.log1p(r).sum())
 
                          
 def annualize_rets(r, periods_per_year):
     """
     Annualizes a set of returns
     We should infer the periods per year
     but that is currently left as an exercise
     to the reader :-)
     """
     compounded_growth = (1+r).prod()
     n_periods = r.shape[0]
     return compounded_growth**(periods_per_year/n_periods)-1
 
 
 def annualize_vol(r, periods_per_year):
     """
     Annualizes the vol of a set of returns
     We should infer the periods per year
     but that is currently left as an exercise
     to the reader :-)
     """
     return r.std()*(periods_per_year**0.5)
 
 
 def sharpe_ratio(r, riskfree_rate, periods_per_year):
     """
     Computes the annualized sharpe ratio of a set of returns
     """
     # convert the annual riskfree rate to per period
     rf_per_period = (1+riskfree_rate)**(1/periods_per_year)-1
     excess_ret = r - rf_per_period
     ann_ex_ret = annualize_rets(excess_ret, periods_per_year)
     ann_vol = annualize_vol(r, periods_per_year)
     return ann_ex_ret/ann_vol
 
 
 import scipy.stats
 def is_normal(r, level=0.01):
     """
     Applies the Jarque-Bera test to determine if a Series is normal or not
     Test is applied at the 1% level by default
     Returns True if the hypothesis of normality is accepted, False otherwise
     """
     if isinstance(r, pd.DataFrame):
         return r.aggregate(is_normal)
     else:
         statistic, p_value = scipy.stats.jarque_bera(r)
         return p_value > level
 
 
 def drawdown(return_series: pd.Series):
     """Takes a time series of asset returns.
        returns a DataFrame with columns for
        the wealth index, 
        the previous peaks, and 
        the percentage drawdown
     """
     wealth_index = 1000*(1+return_series).cumprod()
     previous_peaks = wealth_index.cummax()
     drawdowns = (wealth_index - previous_peaks)/previous_peaks
     return pd.DataFrame({"Wealth": wealth_index, 
                          "Previous Peak": previous_peaks, 
                          "Drawdown": drawdowns})
 
 
 def semideviation(r):
     """
     Returns the semideviation aka negative semideviation of r
     r must be a Series or a DataFrame, else raises a TypeError
     """
     if isinstance(r, pd.Series):
         is_negative = r < 0
         return r[is_negative].std(ddof=0)
     elif isinstance(r, pd.DataFrame):
         return r.aggregate(semideviation)
     else:
         raise TypeError("Expected r to be a Series or DataFrame")
 
 
 def var_historic(r, level=5):
     """
     Returns the historic Value at Risk at a specified level
     i.e. returns the number such that "level" percent of the returns
     fall below that number, and the (100-level) percent are above
     """
     if isinstance(r, pd.DataFrame):
         return r.aggregate(var_historic, level=level)
     elif isinstance(r, pd.Series):
         return -np.percentile(r, level)
     else:
         raise TypeError("Expected r to be a Series or DataFrame")
 
 
 def cvar_historic(r, level=5):
     """
     Computes the Conditional VaR of Series or DataFrame
     """
     if isinstance(r, pd.Series):
         is_beyond = r <= -var_historic(r, level=level)
         return -r[is_beyond].mean()
     elif isinstance(r, pd.DataFrame):
         return r.aggregate(cvar_historic, level=level)
     else:
         raise TypeError("Expected r to be a Series or DataFrame")
 
 
 from scipy.stats import norm
 def var_gaussian(r, level=5, modified=False):
     """
     Returns the Parametric Gauusian VaR of a Series or DataFrame
     If "modified" is True, then the modified VaR is returned,
     using the Cornish-Fisher modification
     """
     # compute the Z score assuming it was Gaussian
     z = norm.ppf(level/100)
     if modified:
         # modify the Z score based on observed skewness and kurtosis
         s = skewness(r)
         k = kurtosis(r)
         z = (z +
                 (z**2 - 1)*s/6 +
                 (z**3 -3*z)*(k-3)/24 -
                 (2*z**3 - 5*z)*(s**2)/36
             )
     return -(r.mean() + z*r.std(ddof=0))
 
 
 def portfolio_return(weights, returns):
     """
     Computes the return on a portfolio from constituent returns and weights
     weights are a numpy array or Nx1 matrix and returns are a numpy array or Nx1 matrix
     """
     return weights.T @ returns
 
 
 def portfolio_vol(weights, covmat):
     """
     Computes the vol of a portfolio from a covariance matrix and constituent weights
     weights are a numpy array or N x 1 maxtrix and covmat is an N x N matrix
     """
     vol = (weights.T @ covmat @ weights)**0.5
     return vol 
 
 
 def plot_ef2(n_points, er, cov):
     """
     Plots the 2-asset efficient frontier
     """
     if er.shape[0] != 2 or er.shape[0] != 2:
         raise ValueError("plot_ef2 can only plot 2-asset frontiers")
     weights = [np.array([w, 1-w]) for w in np.linspace(0, 1, n_points)]
     rets = [portfolio_return(w, er) for w in weights]
     vols = [portfolio_vol(w, cov) for w in weights]
     ef = pd.DataFrame({
         "Returns": rets, 
         "Volatility": vols
     })
     return ef.plot.line(x="Volatility", y="Returns", style=".-")
 
 
 from scipy.optimize import minimize
 
 def minimize_vol(target_return, er, cov):
     """
     Returns the optimal weights that achieve the target return
     given a set of expected returns and a covariance matrix
     """
     n = er.shape[0]
     init_guess = np.repeat(1/n, n)
     bounds = ((0.0, 1.0),) * n # an N-tuple of 2-tuples!
     # construct the constraints
     weights_sum_to_1 = {'type': 'eq',
                         'fun': lambda weights: np.sum(weights) - 1
     }
     return_is_target = {'type': 'eq',
                         'args': (er,),
                         'fun': lambda weights, er: target_return - portfolio_return(weights,er)
     }
     weights = minimize(portfolio_vol, init_guess,
                        args=(cov,), method='SLSQP',
                        options={'disp': False},
                        constraints=(weights_sum_to_1,return_is_target),
                        bounds=bounds)
     return weights.x
 
 
 def tracking_error(r_a, r_b):
     """
     Returns the Tracking Error between the two return series
     """
     return np.sqrt(((r_a - r_b)**2).sum())
 
                          
 def msr(riskfree_rate, er, cov):
     """
     Returns the weights of the portfolio that gives you the maximum sharpe ratio
     given the riskfree rate and expected returns and a covariance matrix
     """
     n = er.shape[0]
     init_guess = np.repeat(1/n, n)
     bounds = ((0.0, 1.0),) * n # an N-tuple of 2-tuples!
     # construct the constraints
     weights_sum_to_1 = {'type': 'eq',
                         'fun': lambda weights: np.sum(weights) - 1
     }
     def neg_sharpe(weights, riskfree_rate, er, cov):
         """
         Returns the negative of the sharpe ratio
         of the given portfolio
         """
         r = portfolio_return(weights, er)
         vol = portfolio_vol(weights, cov)
         return -(r - riskfree_rate)/vol
     
     weights = minimize(neg_sharpe, init_guess,
                        args=(riskfree_rate, er, cov), method='SLSQP',
                        options={'disp': False},
                        constraints=(weights_sum_to_1,),
                        bounds=bounds)
     return weights.x
 
 
 def gmv(cov):
     """
     Returns the weights of the Global Minimum Volatility portfolio
     given a covariance matrix
     """
     n = cov.shape[0]
     return msr(0, np.repeat(1, n), cov)
 
 
 def optimal_weights(n_points, er, cov):
     """
     Returns a list of weights that represent a grid of n_points on the efficient frontier
     """
     target_rs = np.linspace(er.min(), er.max(), n_points)
     weights = [minimize_vol(target_return, er, cov) for target_return in target_rs]
     return weights
 
 
 def plot_ef(n_points, er, cov, style='.-', legend=False, show_cml=False, riskfree_rate=0, show_ew=False, show_gmv=False):
     """
     Plots the multi-asset efficient frontier
     """
     weights = optimal_weights(n_points, er, cov)
     rets = [portfolio_return(w, er) for w in weights]
     vols = [portfolio_vol(w, cov) for w in weights]
     ef = pd.DataFrame({
         "Returns": rets, 
         "Volatility": vols
     })
     ax = ef.plot.line(x="Volatility", y="Returns", style=style, legend=legend)
     if show_cml:
         ax.set_xlim(left = 0)
         # get MSR
         w_msr = msr(riskfree_rate, er, cov)
         r_msr = portfolio_return(w_msr, er)
         vol_msr = portfolio_vol(w_msr, cov)
         # add CML
         cml_x = [0, vol_msr]
         cml_y = [riskfree_rate, r_msr]
         ax.plot(cml_x, cml_y, color='green', marker='o', linestyle='dashed', linewidth=2, markersize=10)
     if show_ew:
         n = er.shape[0]
         w_ew = np.repeat(1/n, n)
         r_ew = portfolio_return(w_ew, er)
         vol_ew = portfolio_vol(w_ew, cov)
         # add EW
         ax.plot([vol_ew], [r_ew], color='goldenrod', marker='o', markersize=10)
     if show_gmv:
         w_gmv = gmv(cov)
         r_gmv = portfolio_return(w_gmv, er)
         vol_gmv = portfolio_vol(w_gmv, cov)
         # add EW
         ax.plot([vol_gmv], [r_gmv], color='midnightblue', marker='o', markersize=10)
         
         return ax
 
                          
 def run_cppi(risky_r, safe_r=None, m=3, start=1000, floor=0.8, riskfree_rate=0.03, drawdown=None):
     """
     Run a backtest of the CPPI strategy, given a set of returns for the risky asset
     Returns a dictionary containing: Asset Value History, Risk Budget History, Risky Weight History
     """
     # set up the CPPI parameters
     dates = risky_r.index
     n_steps = len(dates)
     account_value = start
     floor_value = start*floor
     peak = account_value
     if isinstance(risky_r, pd.Series): 
         risky_r = pd.DataFrame(risky_r, columns=["R"])
 
     if safe_r is None:
         safe_r = pd.DataFrame().reindex_like(risky_r)
         safe_r.values[:] = riskfree_rate/12 # fast way to set all values to a number
     # set up some DataFrames for saving intermediate values
     account_history = pd.DataFrame().reindex_like(risky_r)
     risky_w_history = pd.DataFrame().reindex_like(risky_r)
     cushion_history = pd.DataFrame().reindex_like(risky_r)
     floorval_history = pd.DataFrame().reindex_like(risky_r)
     peak_history = pd.DataFrame().reindex_like(risky_r)
 
     for step in range(n_steps):
         if drawdown is not None:
             peak = np.maximum(peak, account_value)
             floor_value = peak*(1-drawdown)
         cushion = (account_value - floor_value)/account_value
         risky_w = m*cushion
         risky_w = np.minimum(risky_w, 1)
         risky_w = np.maximum(risky_w, 0)
         safe_w = 1-risky_w
         risky_alloc = account_value*risky_w
         safe_alloc = account_value*safe_w
         # recompute the new account value at the end of this step
         account_value = risky_alloc*(1+risky_r.iloc[step]) + safe_alloc*(1+safe_r.iloc[step])
         # save the histories for analysis and plotting
         cushion_history.iloc[step] = cushion
         risky_w_history.iloc[step] = risky_w
         account_history.iloc[step] = account_value
         floorval_history.iloc[step] = floor_value
         peak_history.iloc[step] = peak
     risky_wealth = start*(1+risky_r).cumprod()
     backtest_result = {
         "Wealth": account_history,
         "Risky Wealth": risky_wealth, 
         "Risk Budget": cushion_history,
         "Risky Allocation": risky_w_history,
         "m": m,
         "start": start,
         "floor": floor,
         "risky_r":risky_r,
         "safe_r": safe_r,
         "drawdown": drawdown,
         "peak": peak_history,
         "floor": floorval_history
     }
     return backtest_result
 
 
 def summary_stats(r, riskfree_rate=0.03):
     """
     Return a DataFrame that contains aggregated summary stats for the returns in the columns of r
     """
     ann_r = r.aggregate(annualize_rets, periods_per_year=12)
     ann_vol = r.aggregate(annualize_vol, periods_per_year=12)
     ann_sr = r.aggregate(sharpe_ratio, riskfree_rate=riskfree_rate, periods_per_year=12)
     dd = r.aggregate(lambda r: drawdown(r).Drawdown.min())
     skew = r.aggregate(skewness)
     kurt = r.aggregate(kurtosis)
     cf_var5 = r.aggregate(var_gaussian, modified=True)
     hist_cvar5 = r.aggregate(cvar_historic)
     return pd.DataFrame({
         "Annualized Return": ann_r,
         "Annualized Vol": ann_vol,
         "Skewness": skew,
         "Kurtosis": kurt,
         "Cornish-Fisher VaR (5%)": cf_var5,
         "Historic CVaR (5%)": hist_cvar5,
         "Sharpe Ratio": ann_sr,
         "Max Drawdown": dd
     })
 
                          
 def gbm(n_years = 10, n_scenarios=1000, mu=0.07, sigma=0.15, steps_per_year=12, s_0=100.0, prices=True):
     """
     Evolution of Geometric Brownian Motion trajectories, such as for Stock Prices through Monte Carlo
     :param n_years:  The number of years to generate data for
     :param n_paths: The number of scenarios/trajectories
     :param mu: Annualized Drift, e.g. Market Return
     :param sigma: Annualized Volatility
     :param steps_per_year: granularity of the simulation
     :param s_0: initial value
     :return: a numpy array of n_paths columns and n_years*steps_per_year rows
     """
     # Derive per-step Model Parameters from User Specifications
     dt = 1/steps_per_year
     n_steps = int(n_years*steps_per_year) + 1
     # the standard way ...
     # rets_plus_1 = np.random.normal(loc=mu*dt+1, scale=sigma*np.sqrt(dt), size=(n_steps, n_scenarios))
     # without discretization error ...
     rets_plus_1 = np.random.normal(loc=(1+mu)**dt, scale=(sigma*np.sqrt(dt)), size=(n_steps, n_scenarios))
     rets_plus_1[0] = 1
     ret_val = s_0*pd.DataFrame(rets_plus_1).cumprod() if prices else rets_plus_1-1
     return ret_val
 
                          
 import statsmodels.api as sm
 def regress(dependent_variable, explanatory_variables, alpha=True):
     """
     Runs a linear regression to decompose the dependent variable into the explanatory variables
     returns an object of type statsmodel's RegressionResults on which you can call
        .summary() to print a full summary
        .params for the coefficients
        .tvalues and .pvalues for the significance levels
        .rsquared_adj and .rsquared for quality of fit
     """
     if alpha:
         explanatory_variables = explanatory_variables.copy()
         explanatory_variables["Alpha"] = 1
     
     lm = sm.OLS(dependent_variable, explanatory_variables).fit()
     return lm
 
 def portfolio_tracking_error(weights, ref_r, bb_r):
     """
     returns the tracking error between the reference returns
     and a portfolio of building block returns held with given weights
     """
     return tracking_error(ref_r, (weights*bb_r).sum(axis=1))
                          
 def style_analysis(dependent_variable, explanatory_variables):
     """
     Returns the optimal weights that minimizes the Tracking error between
     a portfolio of the explanatory variables and the dependent variable
     """
     n = explanatory_variables.shape[1]
     init_guess = np.repeat(1/n, n)
     bounds = ((0.0, 1.0),) * n # an N-tuple of 2-tuples!
     # construct the constraints
     weights_sum_to_1 = {'type': 'eq',
                         'fun': lambda weights: np.sum(weights) - 1
     }
     solution = minimize(portfolio_tracking_error, init_guess,
                        args=(dependent_variable, explanatory_variables,), method='SLSQP',
                        options={'disp': False},
                        constraints=(weights_sum_to_1,),
                        bounds=bounds)
     weights = pd.Series(solution.x, index=explanatory_variables.columns)
     return weights
 
 
 def ff_analysis(r, factors):
     """
     Returns the loadings  of r on the Fama French Factors
     which can be read in using get_fff_returns()
     the index of r must be a (not necessarily proper) subset of the index of factors
     r is either a Series or a DataFrame
     """
     if isinstance(r, pd.Series):
         dependent_variable = r
         explanatory_variables = factors.loc[r.index]
         tilts = regress(dependent_variable, explanatory_variables).params
     elif isinstance(r, pd.DataFrame):
         tilts = pd.DataFrame({col: ff_analysis(r[col], factors) for col in r.columns})
     else:
         raise TypeError("r must be a Series or a DataFrame")
     return tilts
 
 def weight_ew(r, cap_weights=None, max_cw_mult=None, microcap_threshold=None, **kwargs):
     """
     Returns the weights of the EW portfolio based on the asset returns "r" as a DataFrame
     If supplied a set of capweights and a capweight tether, it is applied and reweighted 
     """
     n = len(r.columns)
     ew = pd.Series(1/n, index=r.columns)
     if cap_weights is not None:
         cw = cap_weights.loc[r.index[0]] # starting cap weight
         ## exclude microcaps
         if microcap_threshold is not None and microcap_threshold > 0:
             microcap = cw < microcap_threshold
             ew[microcap] = 0
             ew = ew/ew.sum()
         #limit weight to a multiple of capweight
         if max_cw_mult is not None and max_cw_mult > 0:
             ew = np.minimum(ew, cw*max_cw_mult)
             ew = ew/ew.sum() #reweight
     return ew
 
 def weight_cw(r, cap_weights, **kwargs):
     """
     Returns the weights of the CW portfolio based on the time series of capweights
     """
     w = cap_weights.loc[r.index[1]]
     return w/w.sum()
 
 def backtest_ws(r, estimation_window=60, weighting=weight_ew, verbose=False, **kwargs):
     """
     Backtests a given weighting scheme, given some parameters:
     r : asset returns to use to build the portfolio
     estimation_window: the window to use to estimate parameters
     weighting: the weighting scheme to use, must be a function that takes "r", and a variable number of keyword-value arguments
     """
     n_periods = r.shape[0]
     # return windows
     windows = [(start, start+estimation_window) for start in range(n_periods-estimation_window)]
     weights = [weighting(r.iloc[win[0]:win[1]], **kwargs) for win in windows]
     # convert List of weights to DataFrame
     weights = pd.DataFrame(weights, index=r.iloc[estimation_window:].index, columns=r.columns)
     returns = (weights * r).sum(axis="columns",  min_count=1) #mincount is to generate NAs if all inputs are NAs
     return returns
 
 def sample_cov(r, **kwargs):
     """
     Returns the sample covariance of the supplied returns
     """
     return r.cov()
 
 def weight_gmv(r, cov_estimator=sample_cov, **kwargs):
     """
     Produces the weights of the GMV portfolio given a covariance matrix of the returns 
     """
     est_cov = cov_estimator(r, **kwargs)
     return gmv(est_cov)
 
 def cc_cov(r, **kwargs):
     """
     Estimates a covariance matrix by using the Elton/Gruber Constant Correlation model
     """
     rhos = r.corr()
     n = rhos.shape[0]
     # this is a symmetric matrix with diagonals all 1 - so the mean correlation is ...
     rho_bar = (rhos.values.sum()-n)/(n*(n-1))
     ccor = np.full_like(rhos, rho_bar)
     np.fill_diagonal(ccor, 1.)
     sd = r.std()
     return pd.DataFrame(ccor * np.outer(sd, sd), index=r.columns, columns=r.columns)
 
 def shrinkage_cov(r, delta=0.5, **kwargs):
     """
     Covariance estimator that shrinks between the Sample Covariance and the Constant Correlation Estimators
     """
     prior = cc_cov(r, **kwargs)
     sample = sample_cov(r, **kwargs)
     return delta*prior + (1-delta)*sample