Quantitative Analysis, Risk Management, Modelling, Algo Trading, and Big Data Analysis

Student t Distributed Linear Value-at-Risk

One of the most underestimated feature of the financial asset distributions is their kurtosis. A rough approximation of the asset return distribution by the Normal distribution becomes often an evident exaggeration or misinterpretations of the facts. And we know that. The problem arises if we investigate a Value-at-Risk (VaR) measure. Within a standard approach, it is computed based on the analytical formula:
$$
\mbox{VaR}_{h, \alpha} = \Phi ^{-1} (1-\alpha)\sqrt{h}\sigma – h\mu
$$ what expresses $(1-\alpha)$ $h$-day VaR, $\Phi ^{-1} (1-\alpha)$ denotes the standardised Normal distribution $(1-\alpha)$ quantile, and $h$ is a time horizon. The majority of stocks display significant excess kurtosis and negative skewness. In plain English that means that the distribution is peaked and likely to have a heavy left tail. In this scenario, the best fit of the Normal probability density function (pdf) to the asset return distribution underestimates the risk accumulated in a far negative territory.

In this post we will introduce the concept of Student t Distributed Linear VaR, i.e. an alternative method to measure VaR by taking into account “a correction” for kurtosis of the asset return distribution. We will use the market stock data of IBM as an exemplary case study and investigate the difference in a standard and non-standard VaR calculation based on the parametric models. It will also be an excellent opportunity to learn how to do it in Python, quickly and effectively.

1. Leptokurtic Distributions and Student t VaR

A leptokurtic distribution is one whose density function has a higher peak and greater mass in the tails than the normal density function of the same variance. In a symmetric unimodal distribution, i.e. one whose density function has only one peak, leptokurtosis is indicated by a positive excess kurtosis (Alexander 2008). The impact of leptokurtosis on VaR is non-negligible. Firstly, assuming high significance levels, e.g. $\alpha \le 0.01$, the estimation of VaR for Normal and leptokurtic distribution would return:
$$
\mbox{VaR}_{lepto,\ \alpha\le 0.01} \gt \mbox{VaR}_{Norm,\ \alpha\le 0.01}
$$
which is fine, however for low significance levels the situation may change:
$$
\mbox{VaR}_{lepto,\ \alpha\ge 0.05} \lt \mbox{VaR}_{Norm,\ \alpha\ge 0.05} \ .
$$
A good example of the leptokurtic distribution is Student t distribution. If it describes the data better and displays significant excess kurtosis, there is a good chance that VaR for high significance levels will be a much better measure of the tail risk. The Student t distribution is formally given by its pdf:
$$
f_\nu(x) = \sqrt{\nu\pi} \Gamma\left(\frac{\nu}{2}\right)^{-1} \Gamma\left(\frac{\nu+1}{2}\right)
(1+\nu^{-1} x^2)^{-(\nu+1)/2}
$$ where $\Gamma$ denotes the gamma function and $\nu$ the degrees of freedom. For $\nu\gt 2$ its variance is $\nu(\nu-2)^{-1}$ and the distribution has a finite excess kurtosis of $\kappa = 6(\nu-4)^{-1}$ for $\nu \gt 4$. It is obvious that for $\nu\to \infty$ the Student t pdf approaches the Normal pdf. Large values of $\nu$ for the fitted Student t pdf in case of the financial asset daily return distribution are unlikely, unless the Normal distribution is the best pdf model indeed (possible to obtain for very large samples).

We find $\alpha$ quantile of the Student t distribution ($\mu=0$ and $\sigma=1$) by the integration of the pdf and usually denote it as $\sqrt{\nu^{-1}(\nu-2)} t_{\nu}^{-1}(\alpha)$ where $t_{\nu}^{-1}$ is the percent point function of the t distribution. Accounting for the fact that in practice we deal with random variables of the form of $X = \mu + \sigma T$ where $X$ would describe the asset return given as a transformation of the standardised Student t random variable, we reach to the formal definition of the $(1-\alpha)$ $h$-day Student t VaR given by:
$$
\mbox{Student}\ t\ \mbox{VaR}_{\nu, \alpha, h} = \sqrt{\nu^{-1}(\nu-2)h}\ t_{\nu}^{-1}(1-\alpha)\sigma – h\mu
$$ Having that, we are ready to use this knowledge in practice and to understand, based on data analysis, the best regimes where one should (and should not) apply the Student t VaR as a more “correct” measure of risk.

2. Normal and Student t VaR for IBM Daily Returns

Within the following case study we will make of use of Yahoo! Finance data provider in order to fetch the price-series of the IBM stock (adjusted close) and transform it into return-series. As the first task we have to verify before heading further down this road will be the sample skewness and kurtosis. We begin, as usual:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
# Student t Distributed Linear Value-at-Risk
#   The case study of VaR at the high significance levels
# (c) 2015 Pawel Lachowicz, QuantAtRisk.com
 
import numpy as np
import math
from scipy.stats import skew, kurtosis, kurtosistest
import matplotlib.pyplot as plt
from scipy.stats import norm, t
import pandas_datareader.data as web
 
 
# Fetching Yahoo! Finance for IBM stock data
data = web.DataReader("IBM", data_source='yahoo',
                  start='2010-12-01', end='2015-12-01')['Adj Close']
cp = np.array(data.values)  # daily adj-close prices
ret = cp[1:]/cp[:-1] - 1    # compute daily returns
 
# Plotting IBM price- and return-series
plt.figure(num=2, figsize=(9, 6))
plt.subplot(2, 1, 1)
plt.plot(cp)
plt.axis("tight")
plt.ylabel("IBM Adj Close [USD]")
plt.subplot(2, 1, 2)
plt.plot(ret, color=(.6, .6, .6))
plt.axis("tight")
plt.ylabel("Daily Returns")
plt.xlabel("Time Period 2010-12-01 to 2015-12-01 [days]")

where we cover the most recent past 5 years of trading of IBM at NASDAQ. Both series visualised look like:
IBMseries
where, in particular, the bottom one reveals five days where the stock experienced 5% or higher daily loss.

In Python, a correct computation of the sample skewness and kurtosis is possible thanks to the scipy.stats module. As a test you may verify that for a huge sample (e.g. of 10 million) of random variables $X\sim N(0, 1)$, the expected skewness and kurtosis,

from scipy.stats import skew, kurtosis
X = np.random.randn(10000000)
print(skew(X))
print(kurtosis(X, fisher=False))

is indeed equal 0 and 3,

0.0003890933008049
3.0010747070253507

respectively.

As we will convince ourselves in a second, this is not the case for our IBM return sample. Continuing,

31
32
33
34
35
36
37
38
39
40
41
print("Skewness  = %.2f" % skew(ret))
print("Kurtosis  = %.2f" % kurtosis(ret, fisher=False))
# H_0: the null hypothesis that the kurtosis of the population from which the
# sample was drawn is that of the normal distribution kurtosis = 3(n-1)/(n+1)
_, pvalue = kurtosistest(ret)
beta = 0.05
print("p-value   = %.2f" % pvalue)
if(pvalue < beta):
    print("Reject H_0 in favour of H_1 at %.5f level\n" % beta)
else:
    print("Accept H_0 at %.5f level\n" % beta)

we compute the moments and run the (excess) kurtosis significance test. As already explained in the comment, we are interested whether the derived value of kurtosis is significant at $\beta=0.05$ level or not? We find that for IBM, the results reveal:

Skewness  = -0.70
Kurtosis  = 8.39
p-value   = 0.00
Reject H_0 in favour of H_1 at 0.05000 level

id est, we reject the null hypothesis that the kurtosis of the population from which the sample was drawn is that of the normal distribution kurtosis in favour of the alternative hypothesis ($H_1$) that it is not. And if not then we have work to do!

In the next step we fit the IBM daily returns data sample with both Normal pdf and Student t pdf in the following way:

43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
# N(x; mu, sig) best fit (finding: mu, stdev)
mu_norm, sig_norm = norm.fit(ret)
dx = 0.0001  # resolution
x = np.arange(-0.1, 0.1, dx)
pdf = norm.pdf(x, mu_norm, sig_norm)
print("Integral norm.pdf(x; mu_norm, sig_norm) dx = %.2f" % (np.sum(pdf*dx)))
print("Sample mean  = %.5f" % mu_norm)
print("Sample stdev = %.5f" % sig_norm)
print()
 
# Student t best fit (finding: nu)
parm = t.fit(ret)
nu, mu_t, sig_t = parm
pdf2 = t.pdf(x, nu, mu_t, sig_t)
print("Integral t.pdf(x; mu, sig) dx = %.2f" % (np.sum(pdf2*dx)))
print("nu = %.2f" % nu)
print()

where we ensure that the complete integral over both pdfs return 1, as expected:

Integral norm.pdf(x; mu_norm, sig_norm) dx = 1.00
Sample mean  = 0.00014
Sample stdev = 0.01205
 
Integral t.pdf(x; mu, sig) dx = 1.00
nu = 3.66

The best fit of the Student t pdf returns the estimate of 3.66 degrees of freedom.

From this point, the last step is the derivation of two distinct VaR measure, the first one corresponding to the best fit of the Normal pdf and the second one to the best fit of the Student t pdf, respectively. The trick is

61
62
63
64
65
66
67
68
69
# Compute VaR
h = 1  # days
alpha = 0.01  # significance level
StudenthVaR = (h*(nu-2)/nu)**0.5 * t.ppf(1-alpha, nu)*sig_norm - h*mu_norm
NormalhVaR = norm.ppf(1-alpha)*sig_norm - mu_norm
 
lev = 100*(1-alpha)
print("%g%% %g-day Student t VaR = %.2f%%" % (lev, h, StudenthVaR*100))
print("%g%% %g-day Normal VaR    = %.2f%%" % (lev, h, NormalhVaR*100))

that in the calculation of the $(1-\alpha)$ $1$-day Student t VaR ($h=1$) according to:
$$
\mbox{Student}\ t\ \mbox{VaR}_{\nu, \alpha=0.01, h=1} = \sqrt{\nu^{-1}(\nu-2)}\ t_{\nu}^{-1}(0.99)\sigma – \mu
$$ we have to take $\sigma$ and $\mu$ given by sig_norm and mu_norm variable (in the code) and not the ones returned by the t.fit(ret) function in line #54. We find out that:

99% 1-day Student t VaR = 3.19%
99% 1-day Normal VaR    = 2.79%

i.e. VaR derived based on the Normal pdf best fit may underestimate the tail risk for IBM, as suspected in the presence of highly significant excess kurtosis. Since the beauty of the code is aesthetic, the picture is romantic. We visualise the results of our investigation, extending our code with:

71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
plt.figure(num=1, figsize=(11, 6))
grey = .77, .77, .77
# main figure
plt.hist(ret, bins=50, normed=True, color=grey, edgecolor='none')
plt.hold(True)
plt.axis("tight")
plt.plot(x, pdf, 'b', label="Normal PDF fit")
plt.hold(True)
plt.axis("tight")
plt.plot(x, pdf2, 'g', label="Student t PDF fit")
plt.xlim([-0.2, 0.1])
plt.ylim([0, 50])
plt.legend(loc="best")
plt.xlabel("Daily Returns of IBM")
plt.ylabel("Normalised Return Distribution")
# inset
a = plt.axes([.22, .35, .3, .4])
plt.hist(ret, bins=50, normed=True, color=grey, edgecolor='none')
plt.hold(True)
plt.plot(x, pdf, 'b')
plt.hold(True)
plt.plot(x, pdf2, 'g')
plt.hold(True)
# Student VaR line
plt.plot([-StudenthVaR, -StudenthVaR], [0, 3], c='g')
# Normal VaR line
plt.plot([-NormalhVaR, -NormalhVaR], [0, 4], c='b')
plt.text(-NormalhVaR-0.01, 4.1, "Norm VaR", color='b')
plt.text(-StudenthVaR-0.0171, 3.1, "Student t VaR", color='g')
plt.xlim([-0.07, -0.02])
plt.ylim([0, 5])
plt.show()

what brings us to:
IBMrisk
where in the inset, we zoomed in the left tail of the distribution marking both VaR results. Now, it is much clearly visible why at $\alpha=0.01$ the Student t VaR is greater: the effect of leptokurtic distribution. Moreover, the fit of the Student t pdf appears to be a far better parametric model describing the central mass (density) of IBM daily returns.

There is an ongoing debate which one out of these two risk measures should be formally applied when it comes to reporting of 1-day VaR. I do hope that just by presenting the abovementioned results I made you more aware that: (1) it is better to have a closer look at the excess kurtosis for your favourite asset (portfolio) return distribution, and (2) normal is boring.

Someone wise once said: “Don’t be afraid of being different. Be afraid of being the same as everyone else.” The same means normal. Dare to fit non-Normal distributions!

DOWNLOAD
     tStudentVaR.py

REFERENCES
     Alexander, C., 2008, Market Risk Analysis. Volume IV., John Wiley & Sons Ltd

Predicting Heavy and Extreme Losses in Real-Time for Portfolio Holders (1)

The probability of improbable events. The simplicity amongst complexity. The purity in its best form. The ultimate cure for those who trade, for those who invest. Does it exist? Can we compute it? Is it really something impossible? In this post we challenge ourselves to the frontiers of accessible statistics and data analysis in order to find most optimal computable solutions to this enigmatic but mind-draining problem. Since there is no certainty, everything remains shrouded in the veil of probability. The probability we have to face and hope that something unlikely, determined to take place anyhow, eventually will not happen.


Our goal is to calculate the probability of a very rare event (e.g. a heavy and/or extreme loss) in the trading market (e.g. of a stock plummeting 5% or much more) in a specified time-horizon (e.g. on the next day, in one week, in one month, etc.). The probability. Not the certainty of that event.

In this Part 1, first, we look at the tail of an asset return distribution and compress our knowledge on Value-at-Risk (VaR) to extract the essence required to understand why VaR-stuff is not the best card in our deck. Next, we move to a classical Bayes’ theorem which helps us to derive a conditional probability of a rare event given… yep, another event that (hypothetically) will take place. Eventually, in Part 2, we will hit the bull between its eyes with an advanced concept taken from the Bayesian approach to statistics and map, in real-time, for any return-series its loss probabilities. Again, the probabilities, not certainties.

1. VaR (not) for Rare Events

In the framework of VaR we take into consideration $T$ days of trading history of an asset. Next, we drive a number (VaR) that describes a loss that is likely to take place with the probability of approximately $\alpha$. “To take place” does not mean here that it will take place. In this approach we try to provide some likelihood (a quantitative measure) of the rare event in a specified time-horizon (e.g. on the next day if daily return-series are under our VaR investigation; a scenario considered in this post).

If by $L$ we denote a loss (in percent) an asset can experience on the next day, then:
$$
\mbox{Pr}(L \le -\mbox{VaR}_{1-\alpha}) = \alpha
$$ would be the probability of a loss of $-L\times D$ dollars where $-L\times D\ge -\mbox{VaR}\times D$, equal, for instance, $\alpha=0.05$ (also referred to as $(1-\alpha)$% VaR measure; $\mbox{VaR}_{95}$, etc.) and $D$ is the position size (money invested in the asset in terms of physical currency, e.g. in dollars). In other words, the historical data can help us to find $\mbox{VaR}_{95}$ given $\alpha$ assuming 5% of chances that $\mbox{VaR}_{1-\alpha}$ will be exceeded on the next day.

In order to illustrate that case and its shortcomings when it comes to the analysis of rare events, let’s look at the 10 year trading history of two stocks in the NASDAQ market: highly volatile CAAS (China Automotive Systems, Inc.; of the market capital of 247M) and highly liquid AAPL (Apple Inc.; of the market capital of 750B). First, we fetch their adjusted close price-series from Yahoo! Finance and derive the corresponding daily return-series utilising Python:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
# Predicting Heavy and Extreme Losses in Real-Time for Portfolio Holders
# (c) 2015 QuantAtRisk.com, by Pawel Lachowicz
#
# heavy1.py
 
import pandas.io.data as web
import matplotlib.pyplot as plt
import numpy as np
from pyvar import findvar, findalpha
 
 
# ---1. Data Processing
 
# fetch and download daily adjusted-close price series for CAAS
#  and AAPL stocks using Yahoo! Finance public data provider
caas = web.DataReader("CAAS", data_source='yahoo',
                      start='2005-05-13', end='2015-05-13')['Adj Close']
aapl = web.DataReader("AAPL", data_source='yahoo',
                      start='2005-05-13', end='2015-05-13')['Adj Close']
 
CAAScp = np.array(caas.values)
AAPLcp = np.array(aapl.values)
 
f = file("data1.dat","wb")
np.save(f, CAAScp)
np.save(f, AAPLcp)
f.close()
 
# read in the data from a file
f = file("data1.dat","rb")
CAAScp = np.load(f)
AAPLcp = np.load(f)
f.close()
 
# compute return-series
retCAAS = CAAScp[1:]/CAAScp[:-1]-1
retAAPL = AAPLcp[1:]/AAPLcp[:-1]-1

The best way to understand the data is by plotting them:

39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
# plotting (figure #1)
#  adjusted-close price-series
fig, ax1 = plt.subplots(figsize=(10, 6))
plt.xlabel("Trading days 13/05/2005-13/05/2015")
plt.plot(CAAScp, '-r', label="CAAS")
plt.axis("tight")
plt.legend(loc=(0.02, 0.8))
plt.ylabel("CAAS Adj Close Price (US$)")
ax2 = ax1.twinx()
plt.plot(AAPLcp, '-', label="AAPL")
plt.legend(loc=(0.02, 0.9))
plt.axis("tight")
plt.ylabel("AAPL Adj Close Price (US$)")
 
# plotting (figure #2)
#  daily return-series
plt.figure(num=2, figsize=(10, 6))
plt.subplot(211)
plt.grid(True)
plt.plot(retCAAS, '-r', label="CAAS")
plt.axis("tight")
plt.ylim([-0.25,0.5])
plt.legend(loc="upper right")
plt.ylabel("CAAS daily returns")
plt.subplot(212)
plt.grid(True)
plt.plot(retAAPL, '-', label="AAPL")
plt.legend(loc="upper right")
plt.axis("tight")
plt.ylim([-0.25,0.5])
plt.ylabel("AAPL daily returns")
plt.xlabel("Trading days 13/05/2005-13/05/2015")

We obtain the price-series
figure_1
and return-series
figure_2

respectively. For the latter plot, by fixing the scaling of both $y$-axes we immediately gain an chance to inspect the number of daily trades closing with heavy losses. Well, at least at first glance and for both directions of trading. In this post we will be considering the long positions only.

Having our data pre-processed we may implement two different strategies to make use of the VaR framework in order to work out the probabilities for tail events. The first one is based on setting $\alpha$ level and finding $-VaR_{1-\alpha}$. Let’s assume $\alpha=0.01$, then the following piece of code

73
74
75
76
77
78
79
80
# ---2. Compuation of VaR given alpha
 
alpha = 0.01
 
VaR_CAAS, bardata1 = findvar(retCAAS, alpha=alpha, nbins=200)
VaR_AAPL, bardata2 = findvar(retAAPL, alpha=alpha, nbins=100)
 
cl = 100.*(1-alpha)

aims at computation of the corresponding numbers making use of the function:

def findvar(ret, alpha=0.05, nbins=100):
    # Function computes the empirical Value-at-Risk (VaR) for return-series
    #   (ret) defined as NumPy 1D array, given alpha
    # (c) 2015 QuantAtRisk.com, by Pawel Lachowicz
    #
    # compute a normalised histogram (\int H(x)dx = 1)
    #  nbins: number of bins used (recommended nbins>50)
    hist, bins = np.histogram(ret, bins=nbins, density=True)
    wd = np.diff(bins)
    # cumulative sum from -inf to +inf
    cumsum = np.cumsum(hist * wd)
    # find an area of H(x) for computing VaR
    crit = cumsum[cumsum <= alpha]
    n = len(crit)
    # (1-alpha)VaR
    VaR = bins[n]
    # supplementary data of the bar plot
    bardata = hist, n, wd
    return VaR, bardata

Here, we create the histogram with a specified number of bins and $-VaR_{1-\alpha}$ is found in an empirical manner. For many reasons this approach is much better than fitting the Normal Distribution to the data and finding VaR based on the integration of that continuous function. It is a well-know fact that such function would underestimate the probabilities in far tail of the return distribution. And the game is all about capturing what is going out there, right?

The results of computation we display as follows:

82
83
print("%g%% VaR (CAAS) = %.2f%%" % (cl, VaR_CAAS*100.))
print("%g%% VaR (AAPL) = %.2f%%\n" % (cl, VaR_AAPL*100.))

i.e.

99% VaR (CAAS) = -9.19%
99% VaR (AAPL) = -5.83%

In order to gain a good feeling of those numbers we display the left-tails of both return distributions

85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
# plotting (figure #3)
#  histograms of daily returns; H(x)
#
plt.figure(num=3, figsize=(10, 6))
c = (.7,.7,.7)  # grey color (RGB)
#
# CAAS
ax = plt.subplot(211)
hist1, bins1 = np.histogram(retCAAS, bins=200, density=False)
widths = np.diff(bins1)
b = plt.bar(bins1[:-1], hist1, widths, color=c, edgecolor="k", label="CAAS")
plt.legend(loc=(0.02, 0.8))
#
# mark in red all histogram values where int_{-infty}^{VaR} H(x)dx = alpha
hn, nb, _ = bardata1
for i in range(nb):
    b[i].set_color('r')
    b[i].set_edgecolor('k')
plt.text(-0.225, 30, "VaR$_{%.0f}$ (CAAS) = %.2f%%" % (cl, VaR_CAAS*100.))
plt.xlim([-0.25, 0])
plt.ylim([0, 50])
#
# AAPL
ax2 = plt.subplot(212)
hist2, bins2 = np.histogram(retAAPL, bins=100, density=False)
widths = np.diff(bins2)
b = plt.bar(bins2[:-1], hist2, widths, color=c, edgecolor="k", label="AAPL")
plt.legend(loc=(0.02, 0.8))
#
# mark in red all histogram bars where int_{-infty}^{VaR} H(x)dx = alpha
hn, nb, wd = bardata2
for i in range(nb):
    b[i].set_color('r')
    b[i].set_edgecolor('k')
plt.text(-0.225, 30, "VaR$_{%.0f}$ (AAPL) = %.2f%%" % (cl, VaR_AAPL*100.))
plt.xlim([-0.25, 0])
plt.ylim([0, 50])
plt.xlabel("Stock Returns (left tail only)")
plt.show()

where we mark our $\alpha=$1% regions in red:
heavy3

As you may notice, so far, we haven’t done much new. A classical textbook example coded in Python. However, the last figure reveals the main players of the game. For instance, there is only 1 event of a daily loss larger than 15% for AAPL while CAAS experienced 4 heavy losses. Much higher 99% VaR for CAAS takes into account those 4 historical events and put more weight on 1-day VaR as estimated in our calculation.

Imagine now that we monitor all those tail extreme/rare events day by day. It’s not too difficult to notice (based on the inspection of Figure #2) that in case of AAPL, the stock recorded its first serious loss of $L \lt -10$% approximately 650 days since the beginning of our “monitoring” which commenced on May 13, 2005. In contrast, CAAS was much more volatile and you needed to wait only ca. 100 days to record the loss of the same magnitude.

If something did not happen, e.g. $L \lt -10$%, the VaR-like measure is highly inadequate measure of probabilities for rare events. Once the event took place, the estimation of VaR changes (is updated) but decreases in time until a new rare event occurs. Let’s illustrate it with Python. This is our second VaR strategy: finding $\alpha$ given the threshold for rare events. We write a simple function that does the job for us:

def findalpha(ret, thr=1, nbins=100):
    # Function computes the probablity P(X<thr)=alpha given threshold
    #   level (thr) and return-series (NumPy 1D array). X denotes the
    #   returns as a rv and nbins is number of bins used for histogram
    # (c) 2015 QuantAtRisk.com, by Pawel Lachowicz
    #
    # compute normalised histogram (\int H(x)dx=1)
    hist, bins = np.histogram(ret, bins=nbins, density=True)
    # compute a default histogram
    hist1, bins1 = np.histogram(ret, bins=nbins, density=False)
    wd = np.diff(bins1)
    x = np.where(bins1 < thr)
    y = np.where(hist1 != 0)
    z = list(set(x[0]).intersection(set(y[0])))
    crit = np.cumsum(hist[z]*wd[z])
    # find alpha
    try:
        alpha = crit[-1]
    except Exception as e:
        alpha = 0
    # count number of events falling into (-inft, thr] intervals
    nevents = np.sum(hist1[z])
    return alpha, nevents

We call it in our main program:

126
127
128
129
130
131
132
133
134
135
136
137
138
# ---3. Computation of alpha, given the threshold for rare events
 
thr = -0.10
 
alpha1, ne1 = findalpha(retCAAS, thr=thr, nbins=200)
alpha2, ne2 = findalpha(retAAPL, thr=thr, nbins=200)
 
print("CAAS:")
print("  Pr( L < %.2f%% ) = %.2f%%" % (thr*100., alpha1*100.))
print("  %g historical event(s)" % ne1)
print("AAPL:")
print("  Pr( L < %.2f%% ) = %.2f%%" % (thr*100., alpha2*100.))
print("  %g historical event(s)" % ne2)

what returns the following results:

CAAS:
  Pr( L < -10.00% ) = 0.76%
  19 historical event(s)
AAPL:
  Pr( L < -10.00% ) = 0.12%
  3 historical event(s)

And that’s great however all these numbers are given as a summary of 10 years of analysed data (May 13/2005 to May/13 2015 in our case). The final picture could be better understood if we could dynamically track in time the changes of both probabilities and the number of rare events. We achieve it by the following not-state-of-the-art code:

141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
# ---4. Mapping alphas for Rare Events
 
alphas = []
nevents = []
for t in range(1,len(retCAAS)-1):
    data = retCAAS[0:t]
    alpha, ne = findalpha(data, thr=thr, nbins=500)
    alphas.append(alpha)
    nevents.append(ne)
 
alphas2 = []
nevents2 = []
for t in range(1,len(retAAPL)-1):
    data = retAAPL[0:t]
    alpha, ne = findalpha(data, thr=thr, nbins=500)
    alphas2.append(alpha)
    nevents2.append(ne)
 
# plotting (figure #4)
#   running probability for rare events
#
plt.figure(num=4, figsize=(10, 6))
ax1 = plt.subplot(211)
plt.plot(np.array(alphas)*100., 'r')
plt.plot(np.array(alphas2)*100.)
plt.ylabel("Pr( L < %.2f%% ) [%%]" % (thr*100.))
plt.axis('tight')
ax2 = plt.subplot(212)
plt.plot(np.array(nevents), 'r', label="CAAS")
plt.plot(np.array(nevents2), label="AAPL")
plt.ylabel("# of Events with L < %.2f%%" % (thr*100.))
plt.axis('tight')
plt.legend(loc="upper left")
plt.xlabel("Trading days 13/05/2005-13/05/2015")
plt.show()

revealing the following picture:
heavy4

It’s a great way of looking at and understanding the far left-tail volatility of the asset under your current investigation. The probability between two rare/extreme events decreases for the obvious reason: along the time axis we include more and more data therefore the return distribution evolves and shifts its mass to the right leaving left tail events less and less probable.

The question remains: if an asset has never experienced an extreme loss of given magnitute, the probability of such event, within our VaR framework, simply remains zero! For example, for losses $L \lt -20$%,
heavy4a
CAAS displays 2 extreme events while AAPL none! Therefore, our estimation of superbly rare daily loss of -20% (or more) for AAPL stock based on 10 years of data is zero or undefined or… completely unsound. Can we do better than this?

2. Classical Conditional Prediction for Rare Events

Let’s consider a case where we want to predict the probability of the asset/stock (CAAS; traded at NASDAQ exchange) falling down more than $-L$%. Previously we have achieved such estimation through the integration of its probability density function. An alternative way is via derivation of the conditional probability, i.e. that CAAS will lose more than $-L$% given that NASDAQ index drops more than $-L$% on the next day.

Formula? Well, Bayes’ formula. That all what we need:
$$
\mbox{Pr}(B|R) = \frac{ \mbox{Pr}(R|B)\ \mbox{Pr}(B) } { \mbox{Pr}(R) } .
$$
Great, now, the meaning of each term. Let $A$ denotes the event of the stock daily return to be between -L% and 0%. Let B is the event of the stock return to be $\lt -L$%. Let $R$ is the event of NASDAQ daily return to be $\lt -L$%. Given that, based on both CAAS and NASDAQ historical time-series, we are able to compute $\mbox{Pr}(A)$, $\mbox{Pr}(B)$, $\mbox{Pr}(R|A)$, $\mbox{Pr}(R|B)$, therefore
$$
\mbox{Pr}(R) = \mbox{Pr}(R|A)\mbox{Pr}(A) + \mbox{Pr}(R|B)\mbox{Pr}(B)
$$ as well. Here, $\mbox{Pr}(R|A)$ would stand for the probability of NASDAQ falling down more than $-L$% given the observation of the CAAS return to be in $(-L;0)$% interval and $\mbox{Pr}(R|B)$ would denote the probability of NASDAQ falling down more than $-L$% given the observation of the CAAS return also dropping more than $-L$% on the same day.

With Bayes’ formula we inverse the engineering, aiming at providing the answer on $\mbox{Pr}(B|R)$ given all available data. Ready to code it? Awesome. Here we go:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
# Predicting Heavy and Extreme Losses in Real-Time for Portfolio Holders
# (c) 2015 QuantAtRisk.com, by Pawel Lachowicz
#
# heavy2.py
 
import pandas.io.data as web
import matplotlib.pyplot as plt
import numpy as np
 
# ---1. Data Processing
 
# fetch and download daily adjusted-close price series for CAAS stock
#  and NASDAQ index using Yahoo! Finance public data provider
'''
caas = web.DataReader("CAAS", data_source='yahoo',
                   start='2005-05-13', end='2015-05-13')['Adj Close']
nasdaq = web.DataReader("^IXIC", data_source='yahoo',
                   start='2005-05-13', end='2015-05-13')['Adj Close']
 
CAAScp = np.array(caas.values)
NASDAQcp = np.array(nasdaq.values)
 
f = file("data2.dat","wb")
np.save(f,CAAScp)
np.save(f,NASDAQcp)
f.close()
'''
 
f = file("data2.dat","rb")
CAAScp = np.load(f)
NASDAQcp = np.load(f)
f.close()
 
# compute the return-series
retCAAS = CAAScp[1:]/CAAScp[:-1]-1
retNASDAQ = NASDAQcp[1:]/NASDAQcp[:-1]-1

The same code as used in Section 1 but instead of AAPL data, we fetch NASDAQ index daily close prices. Let’s plot the return-series:

37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
# plotting (figure #1)
#  return-series for CAAS and NASDAQ index
#
plt.figure(num=2, figsize=(10, 6))
plt.subplot(211)
plt.grid(True)
plt.plot(retCAAS, '-r', label="CAAS")
plt.axis("tight")
plt.ylim([-0.25,0.5])
plt.legend(loc="upper right")
plt.ylabel("CAAS daily returns")
plt.subplot(212)
plt.grid(True)
plt.plot(retNASDAQ, '-', label="NASDAQ")
plt.legend(loc="upper right")
plt.axis("tight")
plt.ylim([-0.10,0.15])
plt.ylabel("NASDAQ daily returns")
plt.xlabel("Trading days 13/05/2005-13/05/2015")
plt.show()

i.e.,
heavy5
where we observe a different trading dynamics for NASDAQ index between 750th and 1000th day (as counted from May/13, 2005). Interestingly, the heaviest losses of NASDAQ are not ideally correlated with those of CAAS. That make this case study more exciting! Please also note that CAAS is not the part of the NASDAQ index (i.e., its component). Bayes’ formula is designed around independent events and, within a fair approximation, we may think of CAAS as an asset ticking off that box here.

What remains and takes a lot of caution is the code that “looks at” the data in a desired way. Here is its final form:

58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
# ---2. Computations of Conditional Probabilities for Rare Events
 
# isolate return-series displaying negative returns solely
#  set 1 for time stamps corresponding to positive returns
nretCAAS = np.where(retCAAS < 0, retCAAS, 1)
nretNASDAQ = np.where(retNASDAQ < 0, retNASDAQ, 1)
 
# set threshold for rare events
thr = -0.065
 
# compute the sets of events
A = np.where(nretCAAS < 0, nretCAAS, 1)
A = np.where(A >= thr, A, 1)
B = np.where(nretCAAS < thr, retCAAS, 1)
R = np.where(nretNASDAQ < thr, retNASDAQ, 1)
nA = float(len(A[A != 1]))
nB = float(len(B[B != 1]))
n = float(len(nretCAAS[nretCAAS != 1]))  # n must equal to nA + nB
# (optional)
print(nA, nB, n == (nA + nB))  # check, if True then proceed further
print(len(A), len(B), len(R))
print
 
# compute the probabilities
pA = nA/n
pB = nB/n
 
# compute the conditional probabilities
pRA = np.sum(np.where(R+A < 0, 1, 0))/n
pRB = np.sum(np.where(R+B < 0, 1, 0))/n
 
pR = pRA*pA + pRB*pB
 
# display results
print("Pr(A)\t = %5.5f%%" % (pA*100.))
print("Pr(B)\t = %5.5f%%" % (pB*100.))
print("Pr(R|A)\t = %5.5f%%" % (pRA*100.))
print("Pr(R|B)\t = %5.5f%%" % (pRB*100.))
print("Pr(R)\t = %5.5f%%" % (pR*100.))
 
if(pR>0):
    pBR = pRB*pB/pR
    print("\nPr(B|R)\t = %5.5f%%" % (pBR*100.))
else:
    print("\nPr(B|R) impossible to be determined. Pr(R)=0.")

Python’s NumPy library helps us in a tremendous way by its smartly designed function of where which we employ in lines #69-72 and #86-87. First we test a logical condition, if it’s evaluated to True we grab the right data (actual returns), else we return 1. That opens for us a couple of shortcuts in finding the number of specific events and making sure we are still on the right side of the force (lines #73-79).

As you can see, in line #66 we specified our threshold level of $L = -6.5$%. Given that, we derive the probabilities:

(1216.0, 88.0, True)
(2516, 2516, 2516)
 
Pr(A)	 = 93.25153%
Pr(B)	 = 6.74847%
Pr(R|A)	 = 0.07669%
Pr(R|B)	 = 0.15337%
Pr(R)	 = 0.08186%
 
Pr(B|R)	 = 12.64368%

The level of -6.5% is pretty random but delivers an interesting founding. Namely, based on 10 years of data there is 12.6% of chances that on May/14 2015 CAAS will lose more than -6.5% if NASDAQ drops by the same amount. How much exactly? We don’t know. It’s not certain. Only probable in 12.6%.

The outcome of $\mbox{Pr}(R|B)$ which is close to zero may support our assumption that CAAS has a negligible “influence” on NASDAQ itself. On the other side of the rainbow, it’s much more difficult to interpret $\mbox{Pr}(R)$, the probability of an “isolated” rare event (i.e., $L<-6.5$%) since its estimation is solely based on two players in the game: CAAS and NASDAQ. However, when computed for $N\ge 100$ stocks outside the NASDAQ index but traded at NASDAQ, such distribution of $\mbox{Pr}(R)$'s would be of great value as an another alternative estimator for rare events (I should write about it a distinct post). Now, back the business. There is one problem with our conditional estimation of rare event probabilities as outlined within this Section. A quick check for $L = -10$% reveals:

Pr(A)	 = 98.69632%
Pr(B)	 = 1.30368%
Pr(R|A)	 = 0.00000%
Pr(R|B)	 = 0.00000%
Pr(R)	 = 0.00000%
 
Pr(B|R) impossible to be determined. Pr(R)=0.

Let’s remind that $R$ is the event of NASDAQ daily return to be $\lt -L$%. A quick look at NASDAQ return-series tells us the whole story in order to address the following questions that come to your mind: why, why, why zero? Well, simply because there was no event in past 10 year history of the index that it slid more than -10%. And we are cooked in the water.

A visual representation of that problem can be obtained by executing the following piece of code:

104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
from pyvar import cpr  # non-at-all PEP8 style ;)
 
prob = []
for t in range(2,len(retCAAS)-1):
    ret1 = retCAAS[0:t]
    ret2 = retNASDAQ[0:t]
    pBR, _ = cpr(ret1, ret2, thr=thr)
    prob.append(pBR)
 
# plotting (figure #2)
#   the conditional probability for rare events given threshold
#
plt.figure(num=2, figsize=(10, 6))
plt.plot(np.array(prob)*100., 'r')
plt.ylabel("Pr(B|R) [%%]")
plt.axis('tight')
plt.title("L$_{thr}$ = %.2f%%" % (thr*100.))
plt.xlabel("Trading days 13/05/2005-13/05/2015")
plt.show()

where we moved all Bayes’-based calculations (lines #58-102) into a function cpr (a part of pyvar.py local library; see Download section below) and repeated our previous experiment, however, this time, $\mbox{Pr}(R|B)$ changing in time given the threshold of $L=-6.5$%. The resulting plot would be:
heavy7
Before ca. 800th day the NASDAQ index was lacking any event of $L\lt -6.5$% therefore the $\mbox{Pr}(R|B)$ could not be determined. After that period we got some heavy hits and the conditional probability could be derived. The end value (as for May 13, 2015) is 12.64%.

Hope amongst Hopelessness?

There is no better title to summarise the outcomes we derived. Our goal was to come up with some (running) probability for very rare event(s) that, in our case, would be expressed as a 1-day loss of $-L$% in trading for a specified financial asset. In the first attempt we engaged VaR-related framework and the estimation of occurrence of the rare event based on the PDF integration. In our second attempt we sought for an alternative solution making use of Bayes’ conditional prediction of the asset’s loss of $-L$% (or more) given the observation of the same event somewhere else (NASDAQ index). For the former we ended up the a running probability of $\mbox{Pr}(L \le L_{thr})$ while for the latter with $\mbox{Pr}(R|B)$ given $L_{thr}$ for both B and R events.

Assuming $L_{thr}$ to be equal to -6.5%, we can present the corresponding computations in the following chart:
heavy8

Well, the situation looks hopeless. It ain’t get better even if, as discussed earlier, the probability of $\mbox{Pr}(R)$ has been added. Both methodologies seem to deliver the answer to the same question however, somehow, we are left confused, disappointed, misled, in tasteless despair…

Stay tuned as in Part 2 we will see the light at the end of the tunnel. I promise.

DOWNLOAD
    heavy1.py, heavy2.py, pyvar.py

RELATED POSTS
    VaR and Expected Shortfall vs. Black Swan
    Extreme VaR for Portfolio Managers

Applied Portfolio VaR Decomposition. (2) Impact vs Moving Elements.


Calculations of daily Value-at-Risk (VaR) for any $N$-asset portfolio, as we have studied it already in Part 1, heavily depend on the covariance matrix we need to estimate. This estimation requires historical return time-series. Often negligible but superbly important question one should ask here is: How long?! How long those return-series ought to be in order to provide us with fair and sufficient information on the daily return distribution of each asset held in the portfolio?

It is intuitive that taking data covering solely past 250 trading days and, for example, past 2500 days (moving element) would have a substantial impact on the calculation of current portfolio VaR. In this short post, we will examine this relation.

Case Study (Supplementary Research)

Using the example of 7-asset portfolio from Part 1, we modify our Matlab code that allows us to download the history of stock (adjusted) close prices covering past 10 calendar years (ca. 3600 days),

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
% Applied Portfolio VaR Decomposition. (2) Impact vs Moving Elements.
% (c) 2015 QuantAtRisk.com, by Pawel Lachowicz
 
 
clear all; close all; clc;
format short
 
% Read the list of the portfolio components
fileID = fopen('portfolio.lst');
tmp = textscan(fileID, '%s');
fclose(fileID);
pc=tmp{1};  % a list as a cell array
 
% fetch stock data for last 10 years since:
t0=735976; % 12-Jan-2015
date2=datestr(t0,'yyyy-mm-dd');        % from
date1=datestr(t0-10*365,'yyyy-mm-dd');  % to
 
%{
% create an empty array for storing stock data
stockd={};
% scan the tickers and fetch the data from Quandl.com
for i=1:length(pc)
    quandlc=['WIKI/',pc{i}];
    fprintf('%4.0f %s\n',i,quandlc);
    % fetch the data of the stock from Quandl
    % using recommanded Quandl's command and
    % saving them directly into Matlab's FTS object (fts)
    fts=0;
    [fts,headers]=Quandl.get(quandlc,'type','fints', ...
                  'authcode','YourQuandlCode',...
                  'start_date',date1,'end_date',date2);
    stockd{i}=fts; % entire FTS object in an array's cell
end
save data2
%}
load data2

Assuming the same initial holdings (position sizes) for each asset:

39
40
41
42
43
44
45
46
% an initial dollar expose per position
d=[   55621.00; ...
     101017.00; ...
      23409.00; ...
    1320814.00; ...
     131145.00; ...
     321124.00; ...
    1046867.00];

we aim at calculation of the portfolio VaR,
$$
\mbox{VaR}_P = 1.65 \sigma_P C = 1.65\left(\boldsymbol{d}’ \boldsymbol{M}_2 \boldsymbol{d} \right)^{1/2}
$$ in a function of time or, alternatively speaking, in a function of varying length of return time-series:

48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
% examine VaR_P as a function of historical data
% taken into account
e=[];
for td=(10*250):-1:60   % loop over trading days
    rs=[];
    for i=1:length(pc)
        cp=fts2mat(stockd{i}.Adj_Close,0);
        rv=cp(2:end,1)./cp(1:end-1,1)-1;
        rs=[rs rv(end-td:end)];
    end
    rs(isnan(rs))=0.0;
    % covariance matrix
    M2=cov(rs);
    % the portfolio volatility
    vol_P=sqrt(d'*M2*d);
    % diversified portfolio VaR
    VaR_P=1.65*vol_P;
    e=[e; td VaR_P];
end

The meaning of the loop is straightforward: How does $\mbox{VaR}_P$ change if we take only last 60 trading days into account, and what is its value when we include more and more data going back in time as far as 10 years?

The best way to illustrate the results is to compare running portfolio VaR with the stock market itself, here reflected by ETF of SPY,

68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
% download close price time-series for SPY (ETF)
% tracking S&P500 Index
[fts,headers]=Quandl.get('GOOG/NYSE_SPY','type','fints', ...
                  'authcode','YourQuandlCode',...
                  'start_date',date1,'end_date',date2);
SPY_cp=fts2mat(fts.Close,0);
 
% plot results
subplot(2,1,1)
xtmp=(2501:-1:1)./250;
plot(xtmp,SPY_cp(end-10*250:end,1));
xlim([min(xtmp) max(xtmp)]);
set(gca,'xdir','reverse');
ylabel('SPY Close Price [$]')
 
subplot(2,1,2)
plot(e(:,1)./250,e(:,2),'r');
set(gca,'xdir','reverse');
ylabel('VaR_P [$]')
xlabel('Time since 12-Jan-2015 [trading years]')

uncovering the big picture as follows:
ex
And now all becomes clear! Firstly, $\mbox{VaR}_P$ of our 7-asset portfolio is a function of data. Secondly, $\mbox{VaR}_P$ changes in the way that we might expect as contrasted with the market behaviour. In this point, note an evident increase in $\mbox{VaR}_P$ value if we include not 6 but 7 years of data. Why? The market went down in 2008 and the impact of it had been imprinted in the daily losses of single stocks. Therefore this contribution enriched the left tail of their return distributions. An inclusion of this information in our covariance matrix in a natural way increases the level of $\mbox{VaR}_P$.

So, how long? How many data should we include to get fair estimation of $\mbox{VaR}_P$? It’s not so easy to address. I would go with a lower bound of at least 1 year. As for the upper limit, it depends. If you are optimistic that in 2015 stocks will keep their upward momentum, look back up to 5 years. If you are a cautious investor, include 7 to 10 years of data.

Footnote

Some of you may still wonder why Impact vs Moving Elements?! Well, I believe the inspiration for a good title may come from the most unexpected direction. It was 2013 when I heard that song for the first time and I loved it! Through ages, beautiful women have always inspired men, therefore let me share my inspiration with you. You can become inspired too ;-)

YouTube Preview Image

Applied Portfolio VaR Decomposition. (1) Marginal and Component VaR.

Risk. The only ingredient of life that makes us growing and pushing outside our comfort zones. In finance, taking the risk is a risky business. Once your money have been invested, you need to keep your eyes on the ball that is rolling. Controlling the risk is the art and is the science: a quantitative approach that allows us to ease the pressure. With some proper tools we can manage our positions in the markets. No matter whether you invest in a long-term horizon or trade many times per day, your portfolio is a subject to the constantly changing levels of risk and asset volatility.


VaR. Value-at-Risk. The most hated and most adorable quantitative measure of risk. Your running security amongst the unpredictability of the floor. A drilling pain in the ass; a few sleepless nights you’ve been so familiarised with. The omnipresent and rarely confirmed assumptions of the assets returns to be normally distributed. Well, we all know the truth but, somehow, we love to settle for the analytical solutions allowing us to “see a less attractive girl” — more attractive.

In this post we will stay in the mystic circle of those assumptions as a starting point. Don’t blame me. This is what they teach us in the top-class textbooks on financial risk management and test our knowledge within FRM or CFA exams. You may pass the exam but do you really, I mean, really understand the effects of what you have just studied in real portfolio trading?

We kick off from a theory on VaR decomposition in the framework of an active $N$-asset portfolio when our initial capital of $C$ has been fully invested in the stock market (as an example). Next, we investigate the portfolio VaR and by computation of its marginal VaR we derive new positions in our portfolio that minimise the portfolio VaR.

The Analytical VaR for N-Asset Portfolio. The Beginning.

Let’s keep our theoretical considerations to the absolute but essential minimum. At the beginning you have $C$ dollars in hand. You decide to invest them all into $N$ assets (stocks, currencies, commodities, etc.). Great. Say, you pick up US equities and buy shares of $N=7$ stocks. The number of shares depends on stock price and each position size (in dollars) varies depending on your risk-return appetite. Therefore, the position size can be understood as a weighting factor, $w_i$, for $i=1,…,N$, such:
$$
\sum_{i=1}^N w_i = 1 \ \ \mbox{or} \ \ 100\mbox{%}
$$ or $d_i=w_iC$ in dollars and $\sum_i d_i=C$. Once the orders have been fulfilled, you are a proud holder of $N$-asset portfolio of stocks. Congratulations!

From now on, every day, when the markets close, you check a daily rate of return for individual assets, $r_{t,i} = P_{i,t}/P_{i,{t-1}}-1$ where by $P$ we denote the stock close price time-series. Your portfolio rate of return, day-to-day i.e. from $t$-1 to $t$, will be:
$$
R_{p,t} = \sum_{i=1}^N w_ir_{i,t} = \boldsymbol{w}’\boldsymbol{r}
$$ or alternatively:
$$
R_{p,t} = \sum_{i=1}^N d_ir_{i,t} = \boldsymbol{d}’\boldsymbol{r}
$$ what provides us with the updated portfolio value. Now, that differs from the portfolio expected return defined as:
$$
\mu_P = E\left(\sum_{i=1}^{N} w_ir_i \right) = \sum_{i=1}^{N} w_i\mu_i = \boldsymbol{w}’\boldsymbol{\mu}
$$ where the expected daily return of the $i$-th asset, $\mu_i$, needs to be estimated based on its historical returns (a return-series). This is the most fragile spot in the practical application of the Modern Portfolio Theory. If we estimate,
$$
\mu_i = \sum_{t’={(t-1)}-T}^{{t-1}} r_{t’}/T
$$ looking back at last $T$ trading days, it is obvious that we will derive different $\mu_i$’s for $T$=252 (1 year) and for $T$=1260 (5 years). That affects directly the computation of the portfolio covariance matrix thus, as we will see below, the (expected) portfolio volatility. Keep that crucial remark in mind!

The expected portfolio variance is:
$$
\sigma_P^2 = \boldsymbol{w}’ \boldsymbol{M}_2 \boldsymbol{w}
$$
where $\boldsymbol{M}_2$ is the covariance matrix $N\times N$ with the individual covariances of $c_{ij}=\mbox{cov}(\{r_i\},\{r_j\})$ between security $i$ and $j$ computed based on the return-series. Therefore,
$$
\sigma_P^2 = [w_1,…,w_N]
\left[
\begin{array}{cccc}
c_{11} & c_{12} & … & c_{1N} \\
… & … & … & … \\
c_{N1} & c_{N2} & … & c_{NN}
\end{array}
\right]
\left[
\begin{array}{cccc}
w_{1} \\
… \\
w_{N}
\end{array}
\right]
$$ what can be expressed in the terms of the dollar exposure as:
$$
\sigma_P^2C^2 = \boldsymbol{d}’ \boldsymbol{M}_2 \boldsymbol{d} \ .
$$ In order to move from the portfolio variance to the portfolio VaR we need to know the distribution of the portfolio returns. And we assume it to be normal (delta-normal model). Given that, the portfolio VaR at the 95% confidence level ($\alpha=1.65$) is:
$$
\mbox{VaR}_P = 1.65 \sigma_P C = 1.65\left(\boldsymbol{d}’ \boldsymbol{M}_2 \boldsymbol{d} \right)^{1/2}
$$ what can be understood as diversified VaR since we use our portfolio to reduce the overall risk of the investment. On the other hand, the 95% individual VaR describing individual risk of each asset in the portfolio is:
$$
\mbox{VaR}_i = 1.65 |\sigma_i d_i|
$$ where $\sigma_i$ represents the asset volatility over past period of $T$.

It is extremely difficult to obtain a portfolio with a perfect correlation (unity) among all its components. In such a case, we would talk about undiversified VaR possible to be computed as a sum of $\mbox{VaR}_i$ when, in addition as a required condition, there is no short position held in our portfolio. Reporting the risk for a portfolio currently alive is often conducted by the risk managers by providing a quotation of both diversified and undiversified VaR.

Now, imagine that we add one unit of an $i$-th asset to our portfolio and we aim at understanding the impact of this action on the portfolio VaR. In this case we talk about a marginal contribution to risk, possible to be derived through the differentiation of the portfolio variance with respect to $w_i$:
$$
\frac{\mbox{d}\sigma_P}{\mbox{d}w_i} = \frac { \left[ \sum_{i=1}^{N} w_i^2\sigma_i^2 + \sum_{i=1}^N \sum_{j=1, j\ne i}^{N} w_iw_jc_{ij} \right]^{1/2} }
{ \mbox{d}w_i }
$$
$$
= \frac{1}{2} \left[ 2w_i\sigma_i^2 + 2 \sum_{j=1, j\ne 1}^{N} w_jc_{ij} \right] \times \\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \
\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \left[ \sum_{i=1}^{N} w_i^2\sigma_i^2 + \sum_{i=1}^N \sum_{j=1, j\ne i}^{N} w_iw_jc_{ij} \right]^{-1/2}
$$
$$
= \frac{w_i^2\sigma_i^2 + \sum_{j=1, j\ne 1}^{N} w_j c_{ij} } {\sigma_P}
= \frac{c_{iP}}{\sigma_P}
$$ where, again, $c_{ij}$ denotes the covariance between any $i$-th and $j$-th asset and $c_{iP}$ between $i$-th component and portfolio. That allows us to define the marginal VaR for $i$-th asset in the portfolio in the following way:
$$
\Delta\mbox{VaR}_i = \frac{\mbox{dVaR}_P}{\mbox{d}d_i} = \frac{\mbox{dVaR}_P}{\mbox{d}w_iC} = \frac{\mbox{d}(\alpha\sigma_P C)}{\mbox{d}w_iC} = \alpha\frac{C}{C}
\frac{\mbox{d}\sigma_P}{\mbox{d}w_i} = \alpha \frac{c_{iP}}{\sigma_P}
$$ that has a close relation to the systematic risk of $i$-th asset as confronted with the portfolio itself, described by:
$$
\beta_i = \frac{c_{iP}}{\sigma_P^2}
$$ therefore
$$
\Delta\mbox{VaR}_i = \alpha \sigma_P \beta_i = \beta_i \frac{ \mbox{VaR}_P } {C}
$$ The best and most practical interpretation of the marginal VaR calculated for all positions in the portfolio would be: the higher $\Delta\mbox{VaR}_i$ the corresponding exposure of the $i$-th component should be reduced to lower the overall portfolio VaR. Simple as that. Hold on, we will see that at work in the calculations a little bit later.

Since the portfolio volatility is a highly nonlinear function of its components, a simplistic computation of individual VaRs and adding the up all together injects the error into portfolio risk estimation. However, with a help of the marginal VaR we gain a tool to capture the fireflies amongst the darkness of the night.

We define the component VaR as:
$$
\mbox{cVaR}_i = \Delta\mbox{VaR}_i d_i = \Delta\mbox{VaR}_i w_i C = \beta_i \frac{ \mbox{VaR}_P w_i C } {C} = \beta_i w_i \mbox{VaR}_P
$$ utilising the observation that:
$$
\sigma_P = \sum_{i=1}^{N} w_i c_{iP} = \sum_{i=1}^{N} w_i \beta_i \sigma_P = \sigma_P \sum_{i=1}^{N} w_i \beta_i = \sigma_P
$$ what very nicely leads us to:
$$
\sum_{i=1}^{N} \mbox{cVaR}_i = \mbox{VaR}_P \sum_{i=1}^{N} w_i \beta_i = \mbox{VaR}_P \ .
$$ Having that, we get the percent contribution of the $i$-th asset to the portfolio VaR as:
$$
i = \frac{ \mbox{cVaR}_i } { \mbox{VaR}_P } = w_i \beta_i
$$ It is possible to calculate the required changes in the positions sizes based on our VaR analysis. We obtain it by re-evaluation of the portfolio marginal VaRs meeting now the following condition:
$$
\Delta\mbox{VaR}_i = \beta_i \frac{ \mbox{VaR}_P } { C } = \ \mbox{constant .}
$$ a subject of finding a solution for a “Risk-Minimising Position” problem. All right. The end of a lecture. Time for coding!

 

Case Study

Do You like sex? I do! So, let’s do it in Matlab this time! First, we pick up randomly seven stocks among US equites (for the simplicity of this example), and their tickers are given in the portfolio.lst text file as follows:

AAPL
DIS
IBM
JNJ
KO
NKE
TXN

and next we download their 3 years of historical data from Quandl.com as the data provider:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
% Applied Portfolio VaR Decomposition. (1) Marginal and Component VaR.
%
% (c) 2015 QuantAtRisk.com, by Pawel Lachowicz
 
clear all; close all; clc;
format short
 
% Read the list of the portfolio components
fileID = fopen('portfolio.lst');
tmp = textscan(fileID, '%s');
fclose(fileID);
pc=tmp{1};  % a list as a cell array
 
% fetch stock data for last 3 years since:
t0=735976; % 12-Jan-2015
date2=datestr(t0,'yyyy-mm-dd');        % from
date1=datestr(t0-3*365,'yyyy-mm-dd');  % to
 
% create an empty array for storing stock data
stockd={};
for i=1:length(pc) % scan the tickers and fetch the data from Quandl.com
    quandlc=['WIKI/',pc{i}];
    fprintf('%4.0f %s\n',i,quandlc);
    % fetch the data of the stock from Quandl
    % using recommanded Quandl's command and
    % saving them directly into Matlab's FTS object (fts)
    fts=0;
    [fts,headers]=Quandl.get(quandlc,'type','fints', ...
                  'authcode','YourQuandlCode',...
                  'start_date',date1,'end_date',date2);
    stockd{i}=fts; % entire FTS object in an array's cell
end
 
% limit data to 3 years of trading days, select the adjusted
% close price time-series and calculate return-series
rs=[];
for i=1:length(pc)
    cp=fts2mat(stockd{i}.Adj_Close,0);
    rv=cp(2:end,1)./cp(1:end-1,1)-1;
    rs=[rs rv(end-(3*5*4*12):end)];
end
rs(isnan(rs))=0.0;

The covariance matrix, $\boldsymbol{M}_2$, computed based on 7 return-series is:

44
45
% covariance matrix
M2=cov(rs);

Let’s now assume some random dollar exposure, i.e. position size for each asset,

47
48
49
50
51
52
53
54
55
56
57
% exposure per position
d=[   55621.00; ...
     101017.00; ...
      23409.00; ...
    1320814.00; ...
     131145.00; ...
     321124.00; ...
    1046867.00]
 
% invested capital of C
C=sum(d)

returning our total exposure of $C=\$$2,999,997 in the stock market which is based on the number of shares purchased and/or our risk-return appetite (in Part II of this series, we will witness how that works in a simulated trading).

We turn a long theory on the VaR decomposition into numbers just within a few lines of code, namely:

59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
% the portfolio volatility
vol_p=sqrt(d'*M2*d);
 
% diversified portfolio VaR
VaR_P=1.65*vol_p;
 
% volatility of assets (based on daily returns)
v=sqrt(diag(cov(rs)));
 
% individual and undiversified portfolio VaR
VaRi=1.65*abs(v.*d);
uVaR_P=sum(VaRi);
 
% the portfolio betas
beta=C*(M2*d)/vol_p^2;
 
% the portfolio marginal VaR
DVaR=1.65*(M2*d)/vol_p;
 
% the component VaR [percent]
cVaR=100*DVaR.*d/VaR_P;
 
% initial positions of assets in portfolio [percent]
orgpos=100*d/C

The last column vector provides us with an information on the original (currently held) rate of exposure of our invested capital in a function of the asset number in portfolio:

orgpos =
          1.85
          3.37
          0.78
         44.03
          4.37
         10.70
         34.90

Before we generate our very first risk report, let’s develop a short code that would allow us to find a solution for “Risk-Minimising Position” problem. Again, we want to find new position sizes (different from orgpos) at the lowest possible level of the portfolio risk. We construct intelligently a dedicated function that does the job for us:

% Finding a solution for the "Risk-Minimasing Position" problem
%   based on portfolio VaR
% (c) 2015 QuantAtRisk, by Pawel Lachowicz
 
function [best_x,best_dVaR,best_VaRp,best_volp]=rmp(M2,d,VaR_p)
    c=sum(d);
    N=length(d);
    best_VaRp=VaR_p;
    for i=1:100                     % an arbitrary number (recommended >50)
        k=0;
        while(k~=1)
            d=random('uniform',0,1,N,1);
            d=d/sum(d)*c;
            vol_p=sqrt(d'*M2*d);
            % diversified VaR (portfolio)
            VaRp=1.65*vol_p;
            dVaR=1.65*(M2*d)/vol_p;
            test=fix(dVaR*100);     % set precision here
            m=fix(mean(test));
            t=(test==m);
            k=sum(t)/N;
            if(VaRp<best_VaRp)
                best_x=d;
                best_dVaR=dVaR;
                best_VaRp=VaRp;
                best_volp=vol_p;
            end
        end
    end
end

We feed it with three input parameters, namely, the covariance matrix of $\boldsymbol{M}_2$, the current dollar exposure per asset of $\boldsymbol{d}$, and the current portfolio VaR. Inside we generate randomly a new dollar exposure column vector, next re-calculate the portfolio VaR and marginal VaR until,
$$
\Delta\mbox{VaR}_i = \beta_i \frac{ \mbox{VaR}_P } { C } = \ \mbox{constant}
$$ condition is met. We need to be sure that a new portfolio VaR is lower than the original one and somehow lowest among a large number of random possibilities. That is why we repeat that step 100 times. In this point, as a technical remark, I would suggest to use more than 50 iterations and to experiment with a precision that is required to find marginal VaRs (or $\beta$’s) to be (nearly) the same. The greater both golden numbers are the longer time of waiting for the final solution.

It is also that apex of the day in daily bank’s portfolio evaluation process that assesses the bank’s market risk and global exposure. A risk manager knows how long it takes to find new risk-minimaised positions and how fast C/C++/GPU solutions need to be implemented to cut the time of those computations. Here, with our Matlab function of rmp, we can feel the pain of this game.

Eventually, we finalise our computations:

84
85
86
[nd,nDVaR,nVaR_P,nvol_P]=rmp(M2,d,VaR_P);
redVaR=100*(nVaR_P-VaR_P)/VaR_P;
newpos=100*nd/C;

followed by a nicely formatted report:

88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
fprintf('\n============================================\n');
fprintf('  Portfolio Risk Report as of %s\n',datestr(t0));
fprintf('============================================\n\n');
fprintf('1-day VaR estimation at %1.0f%% confidence level\n\n',95.0)
fprintf('Number of assets:\t\t %1.0f \n',length(d))
fprintf('Current exposure:\t\t $%1s\n',bankformat(C));
fprintf('Portfolio VaR (undiversified):\t $%1s\n',bankformat(uVaR_P));
fprintf('Portfolio VaR (diversified):\t $%1s\n',bankformat(VaR_P));
fprintf('Component VaR [individual VaR]\n');
for i=1:length(d)
    fprintf('   %s\t%6.2f%%   $%-10s\t[$%-10s]\n',pc{i},cVaR(i), ...
            bankformat(cVaRd(i)),bankformat(VaRi(i)));
end
fprintf('\nRisk-Minimising Position scenario\n');
fprintf('--------------------------------------------------------------');
fprintf('\n\t\t   Original   Marginal     New        Marginal\n');
fprintf('\t\t   position   VaR          position   VaR\n');
for i=1:length(d)
    fprintf('\t  %-4s\t    %6.2f%%   %1.5f\t    %6.2f%%   %1.2f\n', ...
            pc{i},orgpos(i), ...
            DVaR(i),newpos(i),nDVaR(i))
end
fprintf('div VaR\t\t $%1s\t\t $%1s\n',bankformat(VaR_P), ...
        bankformat(nVaR_P));
fprintf('\t\t\t\t\t%10.2f%%\n',redVaR);
fprintf('ann Vol\t\t%10.2f%%\t\t%10.2f%%\n',  ...
        sqrt(252)*vol_P/1e6*1e2,sqrt(252)*nvol_P/1e6*1e2);
fprintf('--------------------------------------------------------------');
fprintf('\n\n');

where bankformat function you can find here.

As a portfolio risk manager you are interested in numbers, numbers, and numbers that, in fact, reveal a lot. It is of paramount importance to analyse the portfolio risk report with attention before taking any further steps or decisions:

 
============================================
  Portfolio Risk Report as of 12-Jan-2015
============================================
 
1-day VaR estimation at 95% confidence level
 
Number of assets:		 7 
Current exposure:		 $2,999,997.00
Portfolio VaR (undiversified):	 $54,271.26
Portfolio VaR (diversified):	 $33,027.94
Component VaR [individual VaR]
   AAPL	  0.61%   $202.99    	[$1,564.29  ]
   DIS	  1.71%   $564.13    	[$1,885.84  ]
   IBM	  0.26%   $84.43     	[$429.01    ]
   JNJ	 30.99%   $10,235.65 	[$17,647.63 ]
   KO	  1.19%   $392.02    	[$2,044.05  ]
   NKE	  9.10%   $3,006.55  	[$7,402.58  ]
   TXN	 56.14%   $18,542.15 	[$23,297.83 ]
 
Risk-Minimising Position scenario
--------------------------------------------------------------
		   Original   Marginal     New        Marginal
		   position   VaR          position   VaR
	  AAPL	      1.85%   0.00365	      6.02%   0.01
	  DIS 	      3.37%   0.00558	      1.69%   0.01
	  IBM 	      0.78%   0.00361	     10.79%   0.01
	  JNJ 	     44.03%   0.00775	     37.47%   0.01
	  KO  	      4.37%   0.00299	     27.36%   0.01
	  NKE 	     10.70%   0.00936	      8.13%   0.01
	  TXN 	     34.90%   0.01771	      8.53%   0.01
div VaR		 $33,027.94		 $26,198.39
					    -20.68%
ann Vol		     31.78%		     25.21%
--------------------------------------------------------------

From the risk report we can spot the difference of $\$$21,243.32 between the undiversified and diversified portfolio VaR and the largest contribution of Texas Instruments Incorporated stock (TXN) to the overall portfolio VaR though it is Johnson & Johnson stock (JNJ) that occupies the largest dollar position in our portfolio.

The original marginal VaR is the largest for TXN, NKE, and JNJ. Therefore, in order to minimise portfolio VaR, we should cut them and/or reduce their positions. On the other hand, the exposure should be increased in case of KO, IBM, AAPL, and DIS which display the lowest marginal VaR. Both suggestions find the solution in a form of derived new holdings that would reduce the portfolio VaR by $-$20.68% decreasing concurrently the annualised portfolio volatility by 6.57%.

Indeed. Risk is like sex. The more the merrier! Unless, you are a risk manager and you calculate your “positions”.

Extreme VaR for Portfolio Managers


Regardless of who you are, either an algo trader hidden in the wilderness of Alaska or an active portfolio manager in a large financial institution in São Paulo, your investments are the subject to sudden and undesired movements in the markets. There is one fear. We know it very well. It’s a low-probability, high-impact event. An extreme loss. An avalanche hitting your account balance hard. Really hard.

The financial risk management has developed a set of quantitative tools helping to rise an awareness of the worst case scenarios. The most popular two are Value-at-Risk (VaR) and Expected Shortfall (ES). Both numbers provide you with the amount of stress you put on your shoulders when the worst takes place indeed.

Quoting $VaR_{0.05}^{\rm 1-day}$ value we assume that in 95% of cases we should feel confident that the maximal loss for our portfolio, $C\times VaR_{0.05}^{\rm 1-day}$, will not be exceed in trading on the next day (1-day VaR) where $C$ is a current portfolio value (e.g. $\$$35,292.13). However if unforeseen happens, our best guess is that we will lose $C\times ES_{0.05}^{\rm 1-day}$ of dollars, where $ES_{0.05}^{\rm 1-day}$ stands for the 1-day expected shortfall estimate. Usually $ES_{0.05}^{\rm 1-day} > VaR_{0.05}^{\rm 1-day}$. We illustrated this case previously across VaR and Expected Shortfall vs. Black Swan article using IBM stock trading data. We also discussed a superbly rare case of extreme-loss in the markets on Oct 19, 1987 (-22%), far beyond any measure.

In this post we will underline an additional quantitative risk measure, quite often overlooked by portfolio managers, an extreme value-at-risk (EVaR). Using Matlab and an exemplary portfolio of stocks, we gonna confront VaR, ES, and EVaR for a better understanding of the amount of risk involved in the game. Hitting hard is inevitable under drastic market swings so you need to know the highest pain threshold: EVaR.

Generalised Extreme-Value Theory and Extreme VaR

The extreme-value theory (EVT) plays an important role in today’s financial world. It provides a tailor-made solution for those who wish to deal with highly unexpected events. In a classical VaR-approach, the risk manager is more concentrated and concerned around central tendency statistics. We are surrounded by them every day. The mean values, standard deviations, $\chi^2$ distributions, etc. They are blood under the skin. The point is they are governed by central limit theorems, which, on a further note, cannot be applied to extremes. That is where the EVT enters the game.

From the statistical point of view, the problem is simple. If we assume that $X$ is random variable of the loss and $X$ is independent and identically distributed (iid) by that we also assume a general framework suggesting us that $X$ comes from an unknown distribution of $F(x)=Pr(X\le x)$. EVT makes use of the fundamental Fisher-Trippett theorem and delivers practical solution for a large sample of $n$ extreme random loss variables knows also as a generalised extreme-value (GEV) distribution defined as:
$$
H(z;a,b) = \left\{
\begin{array}{lr}
\exp\left[-\left( 1+ z\frac{x-b}{a} \right)^{-1/z} \right] & : z \ne 0\\
\exp\left[-\exp\left(\frac{x-b}{a} \right) \right] & : z = 0 \\
\end{array}
\right.
$$ Under the condition of $x$ meeting $1+z(x-b)/a > 0$ both parameters $a$ and $b$ in GEV distribution ought to be understood as location and scale parameters of limiting distribution $H$. Their meaning is close but distinct from mean and standard deviation. The last parameter of $z$ corresponds to a tail index. Tht tail is important in $H$. Why? Simply because it points our attention to the heaviness of extreme losses in the data sample.

In our previous article on Black Swan and Extreme Loss Modeling we applied $H$ with $z=0$ to fit the data. This case of GEV is known as the Gumbel distribution referring to $F(x)$ to have exponential (light) tails. If the underlying data fitted with $H$ model reveal $z<0$, we talk about Fréchet distribution. From a practical point of view, the Fréchet case is more useful as it takes into account heavy tails more carefully. By plotting standardised Gumbel and Fréchet distributions ($a=1$, $b=0$), you can notice the gross mass of the $H$ contained between $x$ values of -2 and 6 what translates into $b-2a$ to $b+6a$ interval. With $z>0$, the tail contributes more significantly in production of very large $X$-values. The remaining case for $z<0$, Weibull distribution, we will leave out of our interest here (the case of tails lighter than normal).

If we set the left-hade side of GEV distribution to $p$, take logs of both sides of above equations and rearrange them accordingly, we will end up with a setup of:
$$
\ln(p) = \left\{
\begin{array}{lr}
-\left(1+zk \right)^{-1/z} & : z \ne 0\\
-\exp(-k) & : z = 0 \\
\end{array}
\right.
$$ where $k=(x-b)/a$ and $p$ denotes any chosen cumulative probability. From this point we can derive $x$ values in order to obtain quantiles associated with $p$, namely:
$$
x = \left\{
\begin{array}{ll}
b – \frac{a}{z} \left[1-(-\ln(p))^{-z}\right] & : ({\rm Fréchet}; z>0) \\
b – a\ln\left[-\ln(p)\right] & : ({\rm Gumbel}; z=0) \\
\end{array}
\right.
$$ Note, when $z\rightarrow 0$ then Fréchet quantiles tend to their Gumbel equivalents. Also, please notice that for instance for the standardised Gumbel distribution the 95% quantile is $-\ln[-\ln(0.95)] = 2.9702$ and points at its location in the right-tail, and not in the left-tail as usually associated with heaviest losses. Don’t worry. Gandhi says relax. In practice we fit GEV distribution with reversed sign or we reverse sign for 0.95th quantile for random loss variables of $x$. That point will become more clear with some examples that follow.

We allow ourselves to define the Extreme Value-at-Risk (EVaR) associated directly with fitted GEV distribution, $H_n$, to the data ($n$ data points), passing
$$
Pr[H_n < H'] = p = {Pr[X < H']}^n = [\alpha]^n $$ where $\alpha$ is the VaR confidence level associated with the threshold $H'$, in the following way: $$ \mbox{EVaR} = \left\{ \begin{array}{ll} b_n - \frac{a_n}{z_n} \left[1-(-n\ln(\alpha))^{-nz_n}\right] & : ({\rm Fréchet}; z>0) \\
b_n – a_n\ln\left[-n\ln(\alpha)\right] & : ({\rm Gumbel}; z=0) \\
\end{array}
\right.
$$ These formulas provide a practical way of estimating VaR based the analysis of the distribution and frequency of extreme losses in the data samples and are meant to work best for very high confidence levels. Additionally, as pointed by Dowd and Hauksson, the EVaR does not undergo the rule of square root of a holding period what is applicable for a standard VaR derived for (log-)normal distributions of return series.

1-day EVaR for Assets in your Portfolio

All right. Time to connect the dots. Imagine, you hold a portfolio of 30 bluechips being a part of Dow Jones Industrial index. Before market opens, regardless of quantity and direction of positions you hold in each of those 30 stocks, you, as a cautious portfolio manager, want to know the estimation of extreme loss that can take place on the next trading day.

You look at last 252 trading days of every stock, therefore you are able to extract a sample of $N=30$ largest losses (one per stock over last 365 calendar days). Given this sample, we can fit it with GEV distribution and find the best estimates for $z, a$ and $b$ parameters. However, in practice, $N=30$ is a small sample. There is nothing wrong to extract for each stock 5 worst losses that took place within last 252 trading days. That allows us to consider $n=150$ sample (see the code below).

Let’s use Matlab and the approach described in Create a Portfolio of Stocks based on Google Finance Data fed by Quandl post, here re-quoted in a shorter form:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
% Extreme Value-at-Risk
%
% (c) 2014 QuantAtRisk.com, by Pawel Lachowicz
 
 
close all; clear all; clc;
 
% Read the list of Dow Jones components
fileID = fopen('dowjones.lst');
tmp = textscan(fileID, '%s');
fclose(fileID);
dowjc=tmp{1};  % a list as a cell array
 
% Read in the list of tickers and internal code from Quandl.com
[ndata, text, alldata] = xlsread('QuandlStockCodeListUS.xlsx');
quandlc=text(:,1); % again, as a list in a cell array
quandlcode=text(:,3) % corresponding Quandl's Price Code
 
% fetch stock data for last 365 days
date2=datestr(today,'yyyy-mm-dd')     % from
date1=datestr(today-365,'yyyy-mm-dd') % to
 
stockd={};
% scan the Dow Jones tickers and fetch the data from Quandl.com
for i=1:length(dowjc)
    for j=1:length(quandlc)
        if(strcmp(dowjc{i},quandlc{j}))
            fprintf('%4.0f %s\n',i,quandlc{j});
            fts=0;
            [fts,headers]=Quandl.get(quandlcode{j},'type','fints', ...
                          'authcode','YourQUANDLtocken',...
                          'start_date',date1,'end_date',date2);
            cp=fts2mat(fts.Close,0); % close price
            R{i}=cp(2:end)./cp(1:end-1)-1; % return-series
        end
    end
end
N=length(dowjc);

where make use of external files: dowjones.lst (tickers of stocks in our portfolio) and QuandlStockCodeListUS.xlsx (see here for more details on utilising Quandl.com‘s API for fetching Google Finance stock data).

The Matlab’s cell array of $R$ holds 30 return-series (each 251 point long). Having that, we find a set of top 5 maximal daily losses for each stock ending up with a data sample of $n=150$ points as mentioned eariler.

40
41
42
43
44
45
Rmin=[];
for i=1:N
    ret=sort(R{i});
    Mn=ret(1:5); % extract 5 worst losses per stock
    Rmin=[Rmin; Mn];
end

Now, we fit the GEV distribution,
$$
H(z;a,b) = \exp\left[-\left( 1+ z\frac{x-b}{a} \right)^{-1/z} \right] \ \ ,
$$ employing a ready-to-use function, gevfit, from Matlab’s Statistics Toolbox:

47
48
49
50
51
% fit GEV distribution
[parmhat,parmci] = gevfit(Rmin)
z=parmhat(1); % tail index
a=parmhat(2); % scale parameter
b=parmhat(3); % location parameter

which returns for date2 = 2014-03-01,

parmhat =
   -0.8378    0.0130   -0.0326
parmci =
   -0.9357    0.0112   -0.0348
   -0.7398    0.0151   -0.0304

id est, the best estimates of the model’s parameters:
$$
z_{150} = -0.8378^{+0.0980}_{-0.0979}, \ \ a_{150} = 0.0130^{+0.0021}_{-0.0018}, \ \ b_{150} = -0.0326^{+0.0022}_{-0.0022} \ \ .
$$ At this point, it is essential to note that negative $z$ indicates Fréchet distribution since we fitted data with negative signs.

Plotting data with model,

53
54
55
56
57
58
59
60
61
62
63
64
65
66
% plot data with the GEV distribution for found z,a,b
x=-1:0.001:1;
hist(Rmin,-1:0.01:0);
h=findobj(gca,'Type','patch');
set(h,'FaceColor',[0.7 0.7 0.7],'EdgeColor',[0.6 0.6 0.6]);
h=findobj(gca,'Type','box');
set(h,'Color','k');
% GEV pdf 
pdf1=gevpdf(x,z,a,b);
y=0.01*length(Rmin)*pdf1;
line(x,y,'color','r'); box on;
xlabel('R_{min}');
ylabel(['f(R_{min}|z,a,b)']); 
xlim([-0.15 0]);

illustrates our efforts in a clear manner:

gev01

Eventually, performing derivation,

68
69
z=-z;
EVaR=b-(a/z)*(1-(-n*log(0.95))^(-n*z))

we calculate the best estimation of 1-day EVaR as discussed:
$$
\mbox{EVaR} =
b_{150} – \frac{a_{150}}{-z_{150}} \left[1-(-n\ln(0.95))^{nz_{150}}\right] = -0.0481 \ \ .
$$ This number tells us that among 30 bluechips in our portfolio we should expect extreme loss of 4.81% on the following day (estimated using last 252 days of trading history).

Now, is 4.81% a lot or not? That reminds me an immortal question asked by one of my colleagues many years ago, a guy who loved studying philosophy: Is 4 kilograms of potatoes a lot or not? The best way to address it, is by comparison. Let us again quickly look at return series of 30 stocks in our portfolio (R cell array in the code). For each individual series we derive its normalised distribution (histogram), fit it with a normal probability density function, and calculate corresponding 95% VaR (var local variable) in the following way:

71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
%--figure;
VaR=[];
for i=1:N
    rs=R{i};
    % create the histogram for return series
    [H,h]=hist(rs,100);
    sum=0;
    for i=1:length(H)-1
        sum=sum+(H(i)*(h(i+1)-h(i)));
    end
    H=H/sum;
    %--bar(h,H,'FaceColor',[0.7 0.7 0.7],'EdgeColor',[0.7 0.7 0.7])
    % fit the histogram with a normal distribution
    [muhat,sigmahat] = normfit([rs]');
    Y=normpdf(h,muhat(1),sigmahat(1));
    %--hold on; plot(h,Y,'b-'); % and plot the fit using a blue line
    % find VaR for stock 'i'
    sum=0;
    i=1;
    while (i<length(h))
        sum=sum+(Y(i)*(h(i+1)-h(i)));
        if(sum>=0.05) 
            break;
        end
        i=i+1;
    end
    var=h(i)
    %--hold on; plot(h(i),0,'ro','markerfacecolor','r'); % mark VaR_fit
    VaR=[VaR; var];
end

The vector of VaR then holds all thirty 95% VaR. By running the last piece of code,

102
103
104
105
106
107
108
figure;
[H,h]=hist(VaR*100,20);
line=[EVaR*100 0; EVaR*100 6];
bar(h,H,'FaceColor',[0.7 0.7 0.7],'EdgeColor',[0.7 0.7 0.7])
hold on; plot(line(:,1),line(:,2),'r:');
xlim([-6 0]); ylim([0 6]);
xlabel('95% VaR [percent]'); ylabel('Histogram');

we plot their distribution and compare them to EVaR of -4.81%, here marked by red dotted line in the chart:

gev02

The calculated 1-day 95% EVaR of -4.81% places itself much further away from all classical 95% VaR derived for all 30 stocks. It’s an extreme case indeed.

You can use 1-day EVaR to calculate the expected loss in your portfolio just by assuming that every single stock will drop by -4.81% on the next day. This number (in terms of currency of your portfolio) for sure will make you more alert.

Richard Branson once said: If you gonna fail, fail well. Don’t listen to him.

Accelerated Matlab for Quants. (4) Values, Vectors, and Value-at-Risk

YouTube Preview Image

 

VaR and Expected Shortfall vs. Black Swan


It is one of the most fundamental approaches in measuring the risk, but truly worth revising its calculation. Value-at-Risk (VaR). A magical number quoted at the end of the day by the banks’ financial risk managers, portfolios’ risk managers, and risk managers simply concerned about the expected risk threshold not to be, hopefully, exceeded on the next day. What is so magical about it? According to definition, given a specific time horizon and taking into account last $N$ days of our trading activity, we can always calculate a number which would provide us with a simple answer to the question: if something really unexpected happened on the next day, what would be the loss margin we could for sure suffer? Well, welcome to the world of VaR!

More formally, given a daily (or, weekly, monthly, etc.) distribution of returns (of, for example, a single stock, a portfolio of assets, etc.), we are interested in finding the value of $VaR_\alpha$ of, say, $\alpha=0.05$ (five percent) which would say to as that there is 5% of chances that the value of $VaR_{0.05}$ would be exceeded in trading on the next day. This value is located in the left tail of the distribution and, by convention, it is given as positive number. Therefore, we define $VaR_\alpha$ as:
$$
P[X\le VaR_\alpha] = 1 – \alpha
$$ where $X$ is a random variable and can, as in our assumed case, represent a 1-day return on a traded asset or a portfolio of assets.

Beautiful! We have the risk threshold provided as a number. But in practice, how should we calculate $VaR_\alpha$? Here come into spotlight some of the shortcomings of VaR in general. VaR depends on the time-frame of the daily returns we examine: therefore $VaR_\alpha$ as a number will be different if you include last 252 trading days in estimation comparing to last 504 trading days (2 years). Intuitively, the more data (more days, years, etc.) we have the better estimation of $VaR_\alpha$, right? Yes and no. Keep in mind that the dynamic of the trading market changes continuously. In early 80s there was no high-frequency trading (HFT). It dominated trading traffic after the year of 2000 and now HFT plays a key role (an influential factor) in making impact on prices. So, it would be a bit unfair to compare $VaR_\alpha$ estimated in 80s with what is going on right now. $VaR_\alpha$ assumes that the distribution of returns is normal. Not most beautiful assumption but the most straightforward to compute $VaR_\alpha$. Should we use $VaR_\alpha$? Yes, but with caution.

Okay then, having $VaR_\alpha$ calculated we know how far the loss could reach at the $1-\alpha$ confidence level. The next super important question in risk management every risk manager should ask or at least be interested in reporting is: if $VaR_\alpha$ event occurs on the next day, what is the expected loss we might expect to suffer (say, in dollars)? $VaR_\alpha$ is the threshold. Now we are interested in the expected value of loss given the exceedance. It is the definition of expected shortfall and is based on the concept of conditional probability as follows:
$$
E[X\ |\ X>VaR_\alpha] = ES \ .
$$ In general, if our daily distribution of returns can be described by a function $f(x)$ which would represent a power density function (pdf), then:
$$
ES = \frac{1}{\alpha} \int_{-\infty}^{VaR_\alpha} xf(x)dx \ .
$$

Given any data set of returns, $R_t\ (t=1,…,N)$, calculated from our price history,
$$
R_t = \frac{P_{t+1}}{P_t} -1
$$ both numbers, $VaR_\alpha$ and $ES$ can be, in fact, calculated in two ways. The first method is based on the empirical distribution, i.e. using data as given:
$$
VaR_\alpha = h_i^{VaR} \ \ \mbox{for}\ \ \sum_{i=1}^{M-1} H_i(h_{i+1}-h_i) \le \alpha
$$ where $H$ represents the normalized histogram of $R_t$ (i.e., its integral is equal 1) and $M$ is the number of histograms bins. Similarly for $ES$, we get:
$$
ES = \sum_{i=1}^{h_i^{VaR}} h_iH_i(h_{i+1}-h_i) \ .
$$ The second method would be based on integrations given the best fit to the histogram of $R_t$ using $f(x)$ being a normal distribution. As we will see in the practical example below, both approaches returns different values, and an extra caution should be undertaken while reporting the final risk measures.

Case Study: IBM Daily Returns in 1987

The year of 1987 came into history as the most splendid time in the history of stock trading. Why? Simply because of so-called Black Monday on Oct 19, where markets denoted over 20% losses in a single trading day. Let’s analyse the case of daily returns and their distribution for IBM stock, as analysed by a risk manager after the close of the market on Thu, Dec 31, 1978. The risk manager is interested in the estimation of 1-day 95% VaR threshold, i.e. $VaR_0.05$, and the expected shortfall if on the next trading day (Jan 4, 1988) the exceedance event would occur. Therefore, let’s assume that our portfolio is composed solely of the IBM stock and we have the allocation of the capital of $C$ dollars in it as for Dec 31, 1987.

Using a record of historical trading, IBM.dat, of the form:

IBM
DATES 	      OPEN 	HIGH 	  LOW 	    CLOSE 	VOLUME
03-Jan-1984   0.000000 	0.000000  0.000000  30.437500 	0.000000
04-Jan-1984   0.000000 	0.000000  0.000000  30.968750 	0.000000
05-Jan-1984   0.000000 	0.000000  0.000000  31.062500 	0.000000
06-Jan-1984   0.000000 	0.000000  0.000000  30.875000 	0.000000
...
04-Mar-2011   0.000000 	0.000000  0.000000  161.830000 	0.000000
07-Mar-2011   0.000000 	0.000000  0.000000  159.930000 	0.000000

containing the close price of IBM stock, in Matlab, we extract the data for the year of 1987, then we construct a vector of daily returns $R_t$ on the stock,

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
% Reading in the data
fn=['IBM.dat'];
ftsD=ascii2fts(fn,1,2);
range=[datenum('02-Jan-1987'),datenum('31-Dec-1987')];
drr=[datestr(range(1)),'::',datestr(range(2))];
fts=ftsD(drr);
data=fts2mat(fts,1);
 
% R_t vector
rdata=data(2:end,5)./data(1:end-1,5)-1;
 
% Plotting
figure(1);
subplot(2,1,1);
plot(data(:,1)-datenum('2-Jan-1987'),data(:,5));
xlim([0 365]); ylabel('IBM Close Price ($)');
subplot(2,1,2);
plot(data(2:end,1)-datenum('2-Jan-1987'),rdata);
xlim([0 365]); ylabel('R_t')
xlabel('t (days)')

and present the data graphically:

varesR

Having our data, let’s find $VaR_{0.05}$ and corresponding $ES$,

21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
% Construct the normalized histogram
[H,h]=hist(rdata,100);
hist(rdata,100);
sum=0;
for i=1:length(H)-1
    sum=sum+(H(i)*(h(i+1)-h(i)));
end
figure(2)
H=H/sum;
bar(h,H,'FaceColor',[0.7 0.7 0.7],'EdgeColor',[0.7 0.7 0.7])
 
% Calculate VaR 95 based on empirical distribution
sum=0;
i=1;
while (i<length(H))
    sum=sum+(H(i)*(h(i+1)-h(i)));
    if(sum>=0.05) 
        break;
    end
    i=i+1;
end
VaR=h(i)
hold on; plot(h(i),0,'ro'); % mark VaR in the plot
 
% Fit the normal distribution to R_t data
[muhat,sigmahat] = normfit([rdata]');
Y = normpdf(h,muhat(1),sigmahat(1));
hold on; plot(h,Y,'b-'); % and plot the fit using blue line
 
% Find for the fitted N(muhat,sigmahat) VaR 95%
sum=0;
i=1;
while (i<length(h))
    sum=sum+(Y(i)*(h(i+1)-h(i)));
    if(sum>=0.05) 
        break;
    end
    i=i+1;
end
VaR_fit=h(i)
hold on; plot(h(i),0,'ro','markerfacecolor','r'); % mark VaR_fit
ivar=i;
 
% Calculate the Expected Shortfall 
%  based on fitted normal distribution
sum=0;
i=1;
while (i<=ivar)
    sum=sum+(h(i)*Y(i)*(h(i+1)-h(i)));
    i=i+1;
end
ES_fit=sum/0.05
hold on; plot(ES_fit,0,'ko','markerfacecolor','k');
 
% Add numbers to the plot
text(-0.23,25,['1-Day VaR and Expected Shortfall']);
text(-0.23,23.5,['as for 31/12/1987']);
text(-0.23,20,['VaR = ',num2str(VaR*-1,3),' (',num2str(VaR*-100,2),'%)']);
text(-0.23,18,['ES  = ',num2str(ES_fit*-1,3),' (',num2str(ES_fit*-100,2),'%)']);
text(-0.23,13,['\mu = ',num2str(muhat,2)]);
text(-0.23,11,['\sigma = ',num2str(sigmahat,2)]);
xlabel('R_t');

marking all numbers in resulting common plot:

VarES

In the plot, filled red and black markers correspond to $VaR_{0.05}$ and $ES$, respectively, and are calculated using the fitted normal distribution. Both values have been displayed in the plot as well. In addition to that, the plot contains a red open circle marker denoting the value of $VaR_{0.05}$ obtained from the empirical histogram discrete integration. Please note how the choice of the method influences the value of VaR.

The derived numbers tell the risk manager that in 95% of cases one should feel confident that the value of $VaR_{0.05}$ will not be exceed in trading of IBM on Jan 4, 1988 (next trading day). $C\times VaR_{0.05} = 0.028C$ would provide us with the loss amount (in dollars) if such undesired event had occurred. And if it was the case indeed, the expected loss (in dollars) we potentially would suffer, would be $C\times ES = 0.058C$.

Now it easy to understand why the Black Swan event of Oct 19, 1987 falls into the category of super least likely moves in the markets. It places itself in the far left tail of both empirical and fitted distributions denoting an impressive loss of 23.52% for IBM on that day. If we assume that the daily price changes are normally distributed then probability of a drop of at least 22 standard deviations is
$$
P(Z\le -22) = \Phi(-22) = 1.4 \times 10^{-107} \ .
$$ This is astronomically rare event!

They say that a lightening does not strike two times, but if a 1-day drop in the market of the same magnitude had occurred again on Jan 4, 1988, our long position in IBM would cause us to suffer not the loss of $0.058C$ dollars but approximately $0.23C$. Ahhhuc!

Contact Form Powered By : XYZScripts.com