Quantifying Tail Risk

import matplotlib.pyplot as plt
import warnings

import numpy as np
import pandas as pd
import yfinance as yf
import seaborn as sns
from tabulate import tabulate

plt.style.use('seaborn')
plt.rcParams['figure.figsize'] = [8, 4.5]
plt.rcParams['figure.dpi'] = 300
warnings.simplefilter(action='ignore', category=FutureWarning)

%matplotlib inline
%config InlineBackend.figure_format = 'retina'

# Import cufflinks
import cufflinks as cf
cf.set_config_file(offline=True)

ASSETS = ['EURUSD=X', 'EURJPY=X', 'EURSEK=X']
START_DATE = '2016-11-14'
END_DATE = '2019-11-14'

df = yf.download(ASSETS, start=START_DATE, 
                 end=END_DATE, adjusted=True)
print(f'Downloaded {df.shape[0]} rows of data.')

adj_close = df['Adj Close']
adj_close = adj_close.dropna()
adj_close.head()

adj_close.tail()

plot_title = f'Daily Prices:    {START_DATE} - {END_DATE}'
adj_close.iplot(title=plot_title,subplots=True)
plt.tight_layout()
plt.show()

returns = adj_close.pct_change()
returns = returns.dropna()
plot_title = f'Daily Returns:    {START_DATE} - {END_DATE}'
returns.iplot(title=plot_title,subplots=True)
plt.tight_layout()
plt.savefig('ch6_im3.png')
plt.show()

cum_ret = (adj_close[0:]/adj_close.iloc[0])-1
plot_title = f'Cumulative Returns:    {START_DATE} - {END_DATE}'
cum_ret.iplot(title=plot_title)
plt.tight_layout()
plt.savefig('ch6_im3.png')
plt.show()

# Select first prices
first_prices = adj_close.iloc[0]

# Create normalized_prices
normalized = adj_close.div(first_prices)


# Create tickers
tickers = ['EURJPY=X', 'EURSEK=X']

# Subtract the normalized benchmark from the normalized stock prices, and plot the result
normalized[tickers].sub(normalized['EURUSD=X'], axis=0).iplot(title='Relative performance Vs EURUSD')

# Mean Annual Returns
(returns.mean()*252).iplot(kind='bar', title='Annualized Mean Daily Return')

(returns.std()*(252**0.5)).iplot(kind='bar', title='Annualized Mean Daily Volatility')

sh_ratio = (returns.mean()*252)/(returns.std()*(252**0.5))
sh_ratio.iplot(kind='bar', title='Sharpe Ratio')

# Calculate 4 moments of distribution to evaluate normality of returns distribution 
kurt = returns.kurt()
skew = returns.skew()
sigma = returns.std()
mu =  returns.mean()
Norm_df = pd.DataFrame(data=kurt, columns= ["KURT"])
Norm_df["SKEW"] = skew
Norm_df["SIGMA"] = sigma
Norm_df["MU"] = mu
Norm_df

# Calculate the scaled return (return expresses as Z-Score) for chosen security
# to view normality of returns distribution
Sec = 'EURJPY=X'

EmpNorm = pd.DataFrame(data=adj_close[Sec])
EmpNorm['Return'] = EmpNorm[Sec].pct_change()
EmpNorm  = EmpNorm.dropna()

sigma_sec = EmpNorm['Return'].std()
mu_sec =  EmpNorm['Return'].mean()

EmpNorm['Scaled_Return'] = EmpNorm['Return'].apply(lambda x: (x-mu_sec)/sigma_sec)

# Calculate minimum and maximum bin range
sr_min = np.min(EmpNorm['Scaled_Return'])
sr_max = np.max(EmpNorm['Scaled_Return'])

print(f'Minimum {sr_min:.6f}, Maximum {sr_max:.6f}')

# Define bins
from scipy.stats import norm

x = np.linspace(sr_min, sr_max, 100)

# Calculate normal probability density function
y = norm.pdf(x,0,1)

# Plot histogram of scaled returns
plt.hist(EmpNorm['Scaled_Return'], bins=200, density=True, color = 'blue', label = 'Empirical', alpha=1)

# Plot norm pdf
plt.plot(x, y, color = 'orange', label = 'Normal', alpha=1)

# Set x and y axis limits
plt.xlim(-4,4)
plt.ylim(0,0.8)

# Set title
plt.title('Empirical vs Normal Distribution')

# Set legends
plt.legend();

# To add the Pearson Correlation Coefficient between two variables onto the scatterplot. 
# Write a function that takes in two arrays, calculates the statistic, and then draws it on the graph.

from scipy.stats import pearsonr

def corrfunc(x,y, ax=None, **kws):
    """Plot the correlation coefficient in the top left hand corner of a plot."""
    r, _ = pearsonr(x, y)
    ax = ax or plt.gca()
    # Unicode for lowercase rho (ρ)
    rho = '\u03C1'
    # Text, position and coordinate system for annotation
    ax.annotate(f'{rho} = {r:.3f}', xy=(.05, .95), xycoords=ax.transAxes)

corr = returns.corr()
ax = sns.heatmap(
    corr,
    annot = True,
    fmt='.2g',
    linewidths=.5,
    vmin=-1, vmax=1, center=0,
    cmap='RdYlGn_r',
    square=True
    )
ax.set_title("Correlation Matrix", fontsize=20)
ax.set_xticklabels(
    ax.get_xticklabels(),
    rotation=45,
    horizontalalignment='right'
);

# To add the Pearson Correlation Coefficient between two variables onto the scatterplot. 
# Write a function that takes in two arrays, calculates the statistic, and then draws it on the graph.

from scipy.stats import pearsonr


def corrfunc(x,y, ax=None, **kws):
    """Plot the correlation coefficient in the top left hand corner of a plot."""
    r, _ = pearsonr(x, y)
    ax = ax or plt.gca()
    # Unicode for lowercase rho (ρ)
    rho = '\u03C1'
    # Text, position and coordinate system for annotation
    ax.annotate(f'{rho} = {r:.2f}', xy=(.05, .95), xycoords=ax.transAxes)

g = sns.pairplot(returns, kind='scatter', diag_kind="hist")
g.map_lower(corrfunc)
g.map_upper(corrfunc)
plt.show()

from scipy import stats

def regfunc(x,y, ax=None, **kws):
    """Plot the reg coefficient in the top left hand corner of a plot."""
    slope, intercept, rvalue, pvalue, stderr = stats.linregress(x,y)
    ax = ax or plt.gca()
    # R Squared
    r_sq = 'Reg Coeff'
    # Text, position and coordinate system for annotation
    ax.annotate(f'{r_sq} = {rvalue**2:.2f}', xy=(.05, .95), xycoords=ax.transAxes)

h = sns.pairplot(returns, kind='reg', diag_kind="hist")
h.map_upper(regfunc)
h.map_lower(regfunc)
plt.show()

# Covariance matrix (Sigma)
cov_mat = returns.cov()
cov_mat

# Create the weights vector
# Convert 1D vector to 2D arrays (numofassetx1) using np.newaxis and than transpose it
numofasset = len(ASSETS)
wts = numofasset * [1./numofasset]
wts = np.array(wts)[:,np.newaxis]
wts

# Overide Equal Weight?
#wts = np.array([0.33333333,0.33333333,0.33333333])[:,np.newaxis]
#wts

# Calculate Simple Daily Portfolio Returns
port_ret = np.dot(returns,wts)
port_ret_df = pd.DataFrame(data=port_ret, columns=["Simple Port Ret"])
port_ret_df= port_ret_df.dropna()

# Calculate Daily Log Returns for portfolio
port_ret_df["Log Port Ret"] = np.log(1+port_ret_df["Simple Port Ret"])
port_ret_df.head()

# Convert log returns dataframe to numpy array for later use.
port_logret_np = port_ret_df["Log Port Ret"].to_numpy()

# Portfolio Volatility
from numpy.linalg import multi_dot
port_stdev = np.sqrt(multi_dot([wts.T, cov_mat, wts]))
port_stdev = port_stdev.flatten()[0]

print('1D Portfolio Volatility: {0:.4f}%'.format(port_stdev))

# Portfolio returns for securities with given weights
port_ret = np.dot(returns,wts)
port_ret.flatten()
# Expected (Mean) Portfolio Return
port_mean = port_ret.mean()
print("1D Portfolio Expected Return: {0:.4f}%".format(port_mean))

# Calculate 1D Portfolio VaR at different confidence levels
from scipy.stats import norm

pVaR_90 = norm.ppf(1-0.90,port_mean,port_stdev).flatten()[0]
pVaR_95 = norm.ppf(1-0.95,port_mean,port_stdev).flatten()[0]
pVaR_99 = norm.ppf(1-0.99,port_mean,port_stdev).flatten()[0]

# Calculate 10D Portfolio VaR at different confidence levels
pVaR_90_10D = pVaR_90*(10**(1/2))
pVaR_95_10D = pVaR_95*(10**(1/2))
pVaR_99_10D = pVaR_99*(10**(1/2))

pVaR_90_10D = pVaR_90_10D.round(5)
pVaR_95_10D = pVaR_95_10D.round(5)
pVaR_99_10D = pVaR_99_10D.round(5)

# Ouput results in tabular format
ptable = [['90%', pVaR_90_10D],['95%', pVaR_95_10D],['99%', pVaR_99_10D]]
header = ['\033[1mConfidence Level\033[0m', '\033[1mDelta-Normal 10D Value-at-Risk\033[0m']
print(tabulate(ptable,headers=header))

# Visualize Returns Distribution and VaR threshold

mu= port_mean*10
sigma= port_stdev*(10**(1/2))
s = np.random.normal(mu, sigma, 10000000)

plt.figure(figsize=(12,6))
count, bins, ignored = plt.hist(s, 1000, density=True)
plt.plot(bins, 1/(sigma * np.sqrt(2 * np.pi)) *
               np.exp( - (bins - mu)**2 / (2 * sigma**2) ),
         linewidth=3, color='orange')
plt.xlabel('Return rate')
plt.title('Parametric Distribution of Returns with Delta-Normal VaR @ 95 Confidence Interval', fontsize=20)
plt.axvline(x = pVaR_95_10D, color = 'r', linestyle = '--', label = "95% CI of VaR: {0:.4f}%".format(pVaR_95_10D))
plt.legend();
plt.show()

print ("     At a 95% confidence level, losses will not exceed {0:.3f} over a 10 day period.".format(-pVaR_95_10D))

stock_data = yf.download(ASSETS, start=START_DATE, 
                 end=END_DATE, adjusted=True)
print(f'Downloaded {stock_data.shape[0]} rows of data.')

stock_data = stock_data.dropna()
stock_data.head()

# Collect Stock Price Data in an array
stock_price = stock_data.iloc[:,0:3]
stock_price = stock_price.values

array_sum = np.sum(stock_price)
array_has_nan = np.isnan(array_sum)
print(array_has_nan)

# Calculate Daily Returns
nrows = len(stock_price)
# From In all columns, from Stock Price t+1 onwards, divide price by previous price and subtract 1
stock_returns = stock_price[1:nrows,:] / stock_price[0:nrows-1,:] - 1

# Specify number of simulations
mc_rep = 100000
# Specify number of time steps
time_steps = 10

# Specify portfolio weights
portf_WT = np.array([1/3, 1/3, 1/3])
# Call numpy.array.flatten() to convert numpy.array to a 1D array.
portf_WT = wts.flatten()


# Compute Covariance Matrix
cov = np.cov(np.transpose(stock_returns))

# Compute Mean 1D Return for each security
mu = np.mean(stock_returns, axis=0)

# Form a 10 x 3 Matrix, based on 10 time steps for rows and 3 Securities for columns
# Populate with 1D mean return (drift) of each security
Mu = np.full((time_steps,3),mu)

# Transpose to obtain a 3 row x 10 column array representing securities and time steps respectively
Mu = np.transpose(Mu)

# Visualize as Dataframe
Mu_df = pd.DataFrame(data=Mu, index=ASSETS)
Mu_df

# Create initial matrix for price path of portfolio
# With dimensions of 10 rows (for time steps) and 1000 columns (for simulated price paths)
# Populate with zeros
portf_returns_30_m = np.full((time_steps,mc_rep),0.)

# Seed for reproducibility
np.random.seed(100)

for i in range(0,mc_rep):
    # Draw random numbers from Standard Normal distribution for each time step and for each asset
    # forming a 30 x 3 matrix
    Z = np.random.normal(size=3*time_steps)
    # Reshape into a 3 x 30 matrix
    Z = Z.reshape((3,time_steps))
    # Perform the Cholesky decomposition of the covariance matrix and obtain lower matrix
    L = np.linalg.cholesky(cov)
    # Fill the Mu Matrix by shocking the drift with correlated random values
    daily_returns = Mu + np.inner(L,np.transpose(Z))
    # Compute the path of portfolio values by calculating the 
    # cumulative product of the transpose of correlated daily returns by portfolio weights
    portf_Returns_30 = np.cumprod(np.inner(portf_WT,np.transpose(daily_returns)) + 1)
    # Each row of the matrix  represents a simulated price process. 
    # Obtain each simulated Portfolio Return at each time step
    portf_returns_30_m[:,i] = portf_Returns_30

# Correct x-axis for days
x_ticks = np.arange(1,time_steps+1)
tick_list = list(x_ticks)
str_tick_list =[]
for i in tick_list:
    str_tick_list.append(str(i))

plt.figure(figsize=(14,8))
positions = (np.arange(0,time_steps))
labels = str_tick_list
plt.xticks(positions, labels)
plt.plot(portf_returns_30_m)
plt.ylabel('Portfolio Returns')
plt.xlabel('Days')
plt.title('1,000,000 Simulated Portfolio Returns over 10 day horizon', fontsize=20)
plt.show()

# Terminal Values of Portfolio
P_T = portf_returns_30_m[time_steps-1]

P_T_sorted = np.sort(P_T)
percentiles = [10., 5., 1.]
var = np.percentile(P_T, percentiles)

var_list = []
for x, y in zip(percentiles, var):
    var_list.append(y-1)

# Calculate 10D Portfolio VaR at different confidence levels
mVaR_90_10D = var_list[0]
mVaR_95_10D = var_list[1]
mVaR_99_10D = var_list[2]

mVaR_90_10D = mVaR_90_10D.round(5)
mVaR_95_10D = mVaR_95_10D.round(5)
mVaR_99_10D = mVaR_99_10D.round(5)

# Ouput results in tabular format
mtable = [['90%', mVaR_90_10D],['95%', mVaR_95_10D],['99%', mVaR_99_10D]]
header = ['\033[1mConfidence Level\033[0m', '\033[1mMonte Carlo 10D Value-at-Risk\033[0m']
print(tabulate(mtable,headers=header))

# Visualize Returns Distribution
P_T_df = pd.DataFrame(1-P_T)
P_T_df.hist(figsize = (12, 6), bins = 1000, alpha = 1, grid = True)
plt.xlabel('Return rate ')
plt.ylabel('Frequency')
plt.title('Distribution of Simulated Returns with MC VaR @ 95 Confidence Interval', fontsize=20)
plt.axvline(x = mVaR_95_10D, color = 'r', linestyle = '--', label = "95% CI of VaR: {0:.4f}".format(mVaR_95_10D))
plt.legend(framealpha=1, frameon=True);
plt.show()
print ("     Max loss with 95% probability will not exceed than {0:.3f} in single day.".format(-mVaR_95_10D))

# Create function to randonly select subset of 10 integers out of total set integers assocoated with daily returns
def generate_random_index(n=10):
    return np.random.choice(np.arange(len(port_logret_np)), size=n, replace=True)  

# Create list to store simulated terminal values 
VaR_simulated_10 = []
# Specify number of simulations
n_bootstrap = 100000
# For each iteration calculate terminal value after 10 days and store in list
for _ in range(n_bootstrap):
    VaR = port_logret_np[generate_random_index(n=10)].sum()
    VaR_simulated_10.append(VaR)

# Compute historic VaR at desired quantiles
VaR_90 = np.percentile(VaR_simulated_10, q=10).round(5)
VaR_95 = np.percentile(VaR_simulated_10, q=5).round(5)
VaR_99 = np.percentile(VaR_simulated_10, q=1).round(5)

# Ouput results in tabular format
from tabulate import tabulate
htable = [['90%', VaR_90 ],['95%', VaR_95],['99%', VaR_99]]
header = ['\033[1mConfidence Level\033[0m', '\033[1mHistoric 10D Value-at-Risk\033[0m']
print(tabulate(htable,headers=header))

# Visualize Returns Distribution
VaR_10_df = pd.DataFrame(VaR_simulated_10)
VaR_10_df.hist(figsize = (12, 6), bins = 10000, alpha = 1, grid = True)
plt.xlabel('Return rate ')
plt.ylabel('Frequency')
plt.title('Distribution of Bootstrapped Returns with Hist. VaR @ 95 Confidence Interval', fontsize=20)
plt.axvline(x = VaR_95, color = 'r', linestyle = '--', label = "95% CI of VaR: {0:.4f}%".format(VaR_95))
plt.legend(framealpha=1, frameon=True);
plt.show()
print ("     At a 95% confidence level, losses will not exceed {0:.3f} over a 10 day period.".format(-VaR_95))

# Compare Results
my_dict = {'90% Conf. Level' : [VaR_90 ,pVaR_90_10D,  mVaR_90_10D], 
           '95% Conf. Level' : [VaR_95 ,pVaR_95_10D,  mVaR_95_10D], 
           '99% Conf. Level' : [VaR_99 ,pVaR_99_10D,  mVaR_99_10D]}
df = pd.DataFrame.from_dict(my_dict, orient='index')
df.columns = ['Historic', 'Delta-Normal', 'Monte Carlo']
df.round(4)

df.style.set_caption("Comparison of VaR Results").set_table_styles([{
    'selector': 'caption',
    'props': [('color', 'red'),
              ("text-align", "center"),
              ('font-size', '16px')]}])