<jupyter_start><jupyter_text>1.0 Financial returns Investment risk- Risk is a measure of uncertainty. Dispersion / variance of financial returns Measurements1. Standard deviation / variance of daily returns1. Kurtosis of the daily returns distribution1. Skewness of the daily returns distribution1. Historical drawdown A tale of 2 returns- Discrete returns (simple returns) represent periodic (daily, weekly, monthly, etc) price movement- Log returns are often used in financial modeling, they assume continuous compounding    - Log returns are always smaller than discrete returnsDiscrete returns: $ R_{t_2} = \frac{P_{t_2}-P_{t_1}}{P_{t_1}} $Log returns: $ Rl_{t_2} = \frac{ln(P_{t_1})}{ln(P_{t_2})} = ln(P_{t_1}) - ln(P_{t_2}) $<jupyter_code>import pandas as pd
msftPrices = pd.read_csv('data/MSFTPrices.csv', parse_dates=['Date'])
msftPrices = msftPrices.sort_values(by='Date')
msftPrices.set_index('Date', inplace=True)
msftPrices.head()
msftPrices['Returns'] = msftPrices['Adjusted'].pct_change()
msftPrices.head()
import matplotlib.pyplot as plt
plt.hist(msftPrices['Returns'].dropna(), bins=75, density=False)
plt.show()<jupyter_output><empty_output><jupyter_text>1.1 Mean, variance, and normal distributionsProbability distributions have the following moments:1. Mean1. Variance1. Skewness1. KurtosisComparing against a normal distribution- Normal distributions have a skewness near 0 and a kurtosis near 3- Financial returns tend not to be normally distributed- Financial returns can have high kurtosisStandard deviation (volatility)- Investment with higher SD is higher risk investment<jupyter_code>import numpy as np
# Average daily return
np.mean(msftPrices['Returns'])
# Average annualized return assuming 252 trading days
((1+np.mean(msftPrices['Returns']))**252)-1
# Standard deviation
np.std(msftPrices['Returns'])
# Variance
# np.std(msftPrices['Returns'])**2
np.var(msftPrices['Returns'])<jupyter_output><empty_output><jupyter_text>Scaling volatility- Volatility scales with the square root of time- Assume 252 trading days in a year and 21 trading days in a month$ \sigma_{annual} = \sigma_{daily} * \sqrt{252} $$ \sigma_{monthly} = \sigma_{daily} * \sqrt{21} $<jupyter_code># Annualized volatility of returns
np.std(msftPrices['Returns']) * np.sqrt(252)<jupyter_output><empty_output><jupyter_text>1.3 Skewness and kurtosis Skewness- It's the 3rd moment of a distribution- **Negative skew**: Mass of distribution concentrated on the right side, a right-leaning curve- **Positive skew**: Opposite of negative skew<jupyter_code>from scipy.stats import skew
# If skewness > 0, it suggests non-normality
skew(msftPrices['Returns'].dropna())<jupyter_output><empty_output><jupyter_text>Kurtosis- Measure of the thickness of the tails of a distribution- Most financial returns are leptokurtic- **Leptokurtic**: distribution has positive excess kurtosis (> 3)- **Excess kurtosis**: subtract 3 from sample kurtosis<jupyter_code>from scipy.stats import kurtosis
# Excess kurtosis > 0, suggests non-normality
kurtosis(msftPrices['Returns'].dropna())<jupyter_output><empty_output><jupyter_text>Shapiro-Wilk test- Testing for normality- H0: data are normally distributed<jupyter_code>from scipy import stats
pValue = stats.shapiro(msftPrices['Returns'].dropna())[1]
if pValue <= 0.05:
    print('H0 of normality is rejected.')
else:
    print('H0 of normality is accepted.')<jupyter_output>H0 of normality is rejected.<jupyter_text>2.0 Portfolio composition Portfolio return formula$ R_{p} = R_{a_1}w_{a_1} + R_{a_2}w_{a_2} + ... + R_{a_n}w_{a_n} $- $R_{p}$: portfolio return- $R_{a_n}$: return for asset n- $W_{a_n}$: weight for asset n<jupyter_code>import numpy as np
stockReturns = pd.read_csv('data/Big9Returns2017.csv', index_col=['Date'], parse_dates=True).sort_index()
print(stockReturns.info())
stockReturns.head()
stockWeights = pd.read_csv('data/EfficientFrontierPortfoliosSlim.csv')
print(stockWeights.info())
stockWeights.head()
portfolioReturns = pd.DataFrame(stockReturns.iloc[:39].values * stockWeights.iloc[:,:9].values, 
                        index=stockReturns.iloc[:39].index, columns=stockReturns.columns)

portfolioReturns['portfolio_total_return'] = portfolioReturns.sum(axis=1)

portfolioReturns.head()<jupyter_output><empty_output><jupyter_text>Equally weighted portfolioAssuming stockReturns is a Pandas Dataframe of stock returns, calculate the portfolio return for an equally weighted portfolio<jupyter_code>numStocks = len(stockReturns.columns)
equalStockWeights = np.repeat(1/numStocks, numStocks)
equalPortfolioReturns = stockReturns.mul(equalStockWeights, axis=1).sum(axis=1)
equalPortfolioReturns.head()
msftReturns = msftPrices['Returns']
msftReturns.plot(grid=True)
plt.show()
msftCumulativeReturns = ((1+msftReturns).cumprod()-1)
msftCumulativeReturns.plot(grid=True)
plt.show()<jupyter_output><empty_output><jupyter_text>Market capitalization**Market Cap**: The value of a company's publically traded shares Market-cap weighted portfolios$ w_{mcap_n} = \frac{mcap_{n}}{\sum\limits^n_{i=1}mcap_i} $ Market-cap weights```pythonmCap = np.array([100, 200, 100, 100])mCapWeights = mCap / sum(mCap)mCapWeights```  2.1 Correlation and covariance<jupyter_code># Correlation matrix
correlationMatrix = stockReturns.corr()
correlationMatrix
# Correlation matrix heatmap
import seaborn as sns
ax = sns.heatmap(correlationMatrix, cmap='PuBu')<jupyter_output><empty_output><jupyter_text>Portfolio standard deviation$ \sigma_p = \sqrt{w_1^2\sigma_1^2 + w_2^2\sigma_2^2 + 2w_1w_2\rho_{1,2}\sigma_1\sigma_2} $- $\sigma_p$: portfolio standard deviation- $w$: asset weight- $\sigma$: asset volatility- $\rho_{1,2}$: correlation between assets 1 and 2To calculate the portfolio volatility:```pythonportfolioVolatility = np.sqrt(np.dot(weights.T, np.dot(covarianceMatrix, weights)))```<jupyter_code>covarianceMatrix = stockReturns.cov()
covarianceMatrix
# Annualizing the covariance matrix
covarianceMatrixAnnual = covarianceMatrix*252
covarianceMatrixAnnual<jupyter_output><empty_output><jupyter_text>Portfolio standard deviation using covariance Portfolio volatility:$ \sigma_{portfolio} = \sqrt{w_T \bullet \Sigma \bullet w} $- $\sigma_{portfolio}$: portfolio volatility- $\Sigma$: covariance matrix of returns- $w$: portfolio weights ($w_T$ is transposed portfolio weights)- $\bullet$: dot-multiplication operator Matrix transpose$ \begin{bmatrix} 1 & 2 \end{bmatrix}^T = \begin{bmatrix} 1 \\ 2 \end{bmatrix} $ Dot product$ a \bullet b = \sum\limits^{n}_{i=1}a_ib_i = a_1b_1 + a_2b_2 + \ldots + a_nb_n $  2.2 Markowitz portfolios Sharpe ratioThis is a measure of risk-adjusted return$ S = \frac{R_a-r_f}{\sigma_a} $- $S$: sharpe ratio- $R_a$: asset return- $r_f$: risk-free rate of return- $\sigma_a$: asset volatility![Markowitz portfolios](the-markowitz-portfolios.PNG)- **MSR**: max sharpe ratio portfolio- **GMV**: global minimum volatility portfolioTo choose the best portfolio:- Portfolio on the bounding edge of the efficient frontier- Higher return is available if you can stomach higher risk<jupyter_code>numStocks = 9
riskFree = 0
stockWeights['Sharpe'] = (stockWeights['Returns']-riskFree) / stockWeights['Volatility']
stockWeights.head()
MSR = stockWeights.sort_values(by=['Sharpe'], ascending=False)
MSRWeights = MSR.iloc[0, 0:numStocks]
np.array(MSRWeights)<jupyter_output><empty_output><jupyter_text>Past performance is not a guarantee of future returnsMax Sharpe ratio portfolio sound nice, but in reality, returns are difficult to predict<jupyter_code>GMV = stockWeights.sort_values(by=['Volatility'], ascending=True)
GMVWeights = GMV.iloc[0,0:numStocks]
np.array(GMVWeights)<jupyter_output><empty_output><jupyter_text>3.0 The capital asset pricing asset CAPM (Capital asset pricing model)CAPM is the fundamental building block for other models in finance$ E(R_P) - RF = \beta_P(E(R_M) - RF) $- $E(R_p) - RF$: excess expected return of a stock or portfolio P- $E(R_M) - RF$: excess expected return of the broad market portfolio B- $RF$: regional risk-free rate- $\beta_p$: portfolio beta / exposure to the braod market portfolio B Calculating historical beta using covariance$ \beta_P = \frac{cov(R_P,R_B)}{var(R_B)} $- $cov(R_P,R_B)$: covariance between portfolio P and the benchmark market index B- $var(R_B)$: variance of the benchmark market index Excess returns`Excess return = Portfolio return - Risk free return`<jupyter_code>famaFrenchFactors = pd.read_csv('data/FamaFrenchFactors.csv', index_col='Date', parse_dates=['Date'])
print(famaFrenchFactors.info())
famaFrenchFactors.head()
convarianceMatrix = famaFrenchFactors[['Portfolio','Market_Excess']].cov()
convarianceMatrix
convarianceCoefficient = convarianceMatrix.iloc[0,1]
convarianceCoefficient
benchmarkVariance = famaFrenchFactors['Market_Excess'].var()
benchmarkVariance
portfolioBeta = convarianceCoefficient / benchmarkVariance
portfolioBeta<jupyter_output><empty_output><jupyter_text>Calculating beta using linear regression<jupyter_code>import statsmodels.formula.api as smf
model = smf.ols(formula='Portfolio ~ Market_Excess', data=famaFrenchFactors)
fit = model.fit()
beta = fit.params['Market_Excess']
beta<jupyter_output><empty_output><jupyter_text>R-squared vs Adjusted R-squared<jupyter_code>rSquared = fit.rsquared
rSquared
adjustedRSquared = fit.rsquared_adj
adjustedRSquared<jupyter_output><empty_output><jupyter_text>3.1 Alpha and multi-factor models The Fama-French 3 factor model$ R_P = RF + \beta_M(R_M - RF) + b_{SMB} \bullet SMB + b_{HML} \bullet HML + \alpha $- $SMB$: small minus big factor- $b_{SMB}$: exposure to SMB factor- $HML$: high minus low factor- $b_{HML}$: exposure to HML factor- $\alpha$: performance which is unexplained by any other factors- $\beta_M$: beta to the broad market portfolio B<jupyter_code>famaFrenchFactors.head()
model2 = smf.ols(formula='Portfolio ~ Market_Excess + SMB + HML', data=famaFrenchFactors)
fit2 = model2.fit()
adjustedRSquared2 = fit2.rsquared_adj
adjustedRSquared2
fit2.pvalues['HML']
fit2.pvalues['HML'] < 0.05
fit2.params['HML']
fit2.params['SMB']
portfolioAlpha = fit2.params['Intercept']
portfolioAlpha
portfolioAlphaAnnualized = ((1+portfolioAlpha)**252)-1
portfolioAlphaAnnualized<jupyter_output><empty_output><jupyter_text>3.2 Expanding the 3-factor model The Fama-French 5 factor model2 additional factors:- **RMW**: profitability    - returns of companies with high vs low operating profitability- **CMA**: investment    - returns of companies with aggressive vs conservative investments<jupyter_code>model3 = smf.ols(formula='Portfolio ~ Market_Excess + SMB + HML + RMW + CMA', data=famaFrenchFactors)
fit3 = model3.fit()
adjustedRSquared3 = fit3.rsquared_adj
adjustedRSquared3<jupyter_output><empty_output><jupyter_text>4.0 Estimating tail risk**Tail risk** is the risk of extreme investment outcomes- Historical drawdown- Value at risk- Conditional value at risk- Monte-Carlo simulation Historical drawdown![Historical drawdown](historical-drawdown.PNG)**Drawdown** is the percentage loss from the highest cumulative historical point$ Drawdown = \frac{r_t}{RM} - 1 $- $r_t$: cumulative return at time t- **RM**: running maximum<jupyter_code>uso = pd.read_csv('data/USO.csv', index_col='Date', parse_dates=['Date'])
uso.head()
uso.plot(grid=True)
plt.show()
runningMax = np.maximum.accumulate(uso['USO'])
runningMax[runningMax < 1] = 1
drawdown = uso['USO'] / runningMax - 1
drawdown.plot(grid=True)
plt.show()<jupyter_output><empty_output><jupyter_text>Historical value at risk![Historical value](historical-shortfall-vs-value.PNG)**Value at risk (VaR)** is a threshold with a given confidence level that losses will not historically  exceed a certain levelE.g.`VaR(95) = -2.3%`95% certain that losses will not exceed `-2.3%` in a given day based on historical value Historical expected shortfall**Conditional value at risk (cVaR)**: estimate of expected losses sustained in the worst 1-x% of scenarios`cVaR(95) = -2.5%`In worst 5% of cases, losses were on average exceed -2.5% historically.<jupyter_code>varLevel = 95
var95 = np.percentile(stockReturns, 100 - varLevel)
var95
cVar95 = stockReturns[stockReturns <= var95].mean()
cVar95<jupyter_output><empty_output><jupyter_text>4.1 VaR extensions![VaR quantiles](var-quantiles.PNG)**Empirical historical values** are those that have actually occurred. For simulation, use sample from probability distribution<jupyter_code>from scipy.stats import norm

# parametric VaR(95)
mu = np.mean(stockReturns)
std = np.std(stockReturns)
confLevel = 0.05
VaR = norm.ppf(confLevel, mu, std)
VaR<jupyter_output><empty_output><jupyter_text>Scaling risk![Scaling risk](scaling-risk.PNG)Assuming you have a one-day estimate of VaR(95). Estimate 5-day VaR(95)<jupyter_code>forecastDays = 5
VaR5Days = VaR*np.sqrt(forecastDays)
VaR5Days<jupyter_output><empty_output><jupyter_text>4.2 Random walks<jupyter_code>T = 252
S0 = 10
randomReturns = np.random.normal(mu[0], std[0], T) + 1
forecastValues = S0 * randomReturns.cumprod()
forecastValues
# count, bins, ignored = plt.hist(forecastValues, 30, density=True)

# plt.plot(bins, 
#          1/(std * np.sqrt(2 * np.pi)) * np.exp( - (bins - mu[0])**2 / (2 * std[0]**2)), 
#          linewidth=2, color='r')

# plt.show()
'''Calculate the VaR(95) of 100 Monte Carlo simulations'''

mu = 0.0005
vol = 0.001
T = 252
simulationReturns = []

for i in range(100):
    randomReturns = np.random.normal(mu, vol, T)
    simulationReturns.append(randomReturns)

VaR95 = np.percentile(simulationReturns, 5)
VaR95<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.1209271327014218, 'fraction_numerical': 0.020364336492890996, 'mean_word_length': 6.668938103350369, 'pattern_counts': {'":': 0, '<': 92, '<?xml version=': 0, '>': 91, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '22', 'dir': 'jupyter-structured-clean-dedup', 'id': '294421', 'license': 'BSD-2-Clause', 'n_tokens_mistral': 4718, 'n_tokens_neox': 4370, 'n_words': 1150, 'path': 'intro_to_portfolio_risk_management_in_python/intro-to-portfolio-risk-management-in-python.ipynb', 'repo_name': 'shawnngtq/applied-finance', 'size': '13504'}