A monte-carlo approach to the VaR (Value at Risk)

A monte-carlo approach to the VaR (Value at Risk)

Today I decided to take a break from Python technicalities and bring you, instead, a tutorial on a risk analysis technique: The Value at Risk or the VaR.

What is it? The VaR is a measure of risk. For a given financial product or a portfolio of financial products, the VaR is the worst loss you can undergo 99% of the time.

Let me explain further. Say you bought an Apple share with $100\$$.

Now you are exposed to the fluctuations of the stock price. If tomorrow the value of one share becomes $120\$$ you made a profit of $20\%$, on the other hand, if the value drops to $70\$$, you've lost $30\%$.

Being a carefull person, you want to quantify the risk you are exposed to. You can ask yourself: what's the worst possible scenario?. Well, that's silly, because the worst case scenario is clearly loosing all your investement. Another question you may ask yourself is: not considering the $1\%$ worst possibilities, in which case, god be with you, how bad my loss can be? The answer is the Value at Risk.

In the following, will see how you can calculate the VaR for one stock, then for a portfolio of two stocks, which you can easily generalize to the $n$ stocks case.

What you will learn from this tutorial

  1. Getting and processing financial data
  2. Finding the appropriate distribution for your data
  3. Generating correlated* random data

Note: * Generating independent (pseudo-)random data is easy, but generating correlated data, that's some cool sh** right there!

Confession:

I am sorry to say this, but I didn't invent the VaR technique I am going to expose. In fact, I've known the technique for a while, but I've never understood it quit well until I found this video on youtube:

In [1]:
from IPython.display import YouTubeVideo
YouTubeVideo('92WaNz9mPeY')
Out[1]:

(link: https://www.youtube.com/watch?v=92WaNz9mPeY)

The video is part of a free MIT course on finance, very promissing material, recommended to anyone interested in finance and having a basic background in mathematics.

To work: Getting the data

In my opinion, the best way to getting financial data is finding a web API that you can easily query within your code. To our luck there is a totaly free API that does just that, it serves historical financial data for up to the past 20 years.

Go here: http://www.alphavantage.co and click on "GET YOUR FREE API KEY TODAY". Fill the form validate.

Under the validation button a line will appear saying "Welcome to Alpha Vantage! Your API key is: XXXX. Please record this API key for lifetime access to Alpha Vantage.". The XXXX part is your API. Copy it and past it in place of XXXX in the code below:

In [8]:
#Importing the stack we'll use throughout the tutorial
import requests
import datetime
import matplotlib.pyplot as plt
from scipy.stats import norm, cauchy, levy_stable, gennorm
from scipy.stats.mstats import mquantiles
import numpy as np
import statsmodels.api as sm
import math
import pandas as pd

from pandas.plotting import register_matplotlib_converters
register_matplotlib_converters()

api_key = 'XXXX'

The request library allows you to send http requests and fetch http responses. You can play with it and the documentation of the API here.

As to not waste your time, here is a function that fetches the stock daily data for a given company.

In [3]:
def get_data(symbol, api_key):
    params = {
        'function'  : 'TIME_SERIES_DAILY',
        'symbol'    : symbol,
        'outputsize': 'full',
        'apikey'    : api_key
    }
    r = requests.get( "http://www.alphavantage.co/query", params=params  )

    json_data = r.json()

    return json_data[u'Time Series (Daily)']

The first argument of the function above is a symbol for the company. For Apple it is 'AAPL'

In [200]:
#Getting apple data
data = get_data('AAPL', api_key)

You can check that the data is in a dict of dicts. The function bellow takes that dict of dicts and returns a 2-d Numpy array of the dates and closing prices for the past year. (That's all we'll need for today)

In [204]:
def last_year_data(data):
    #We generate the dates of one year in the past (_d)
    today           = datetime.date(2019, 12, 31)
    a_year_before   = today - datetime.timedelta(days=365)

    _d = today
    all_dates = []

    while _d > a_year_before:
        all_dates += [_d]
        _d -= datetime.timedelta(days=1)
    
    #We extract from the data the observations of the year
    one_year_data = []
    
    for _d in all_dates:
        try:
            _date    = _d
            _text_date = datetime.date.strftime(_d , "%Y-%m-%d")
            _close = float(data[_text_date]['4. close'])
            one_year_data += [(_date, _close)]
        except KeyError:
            pass
    
    
    #We return an array of the observations sorted from earliest to older
    return one_year_data

This is the data for the past year:

In [205]:
one_year_data = last_year_data(data)

As a little reward for your efforts, here is a beautifull graph of the APPL shares closing prices for the past year.

In [206]:
dates, prices = (np.array(x) for x in zip(*one_year_data))
In [207]:
plt.plot(dates, prices)
plt.xlabel("Date")
plt.ylabel("Price in $")
plt.title("Apples stock closing prices for the past year")
plt.show()

Statistical analysis

Now that you have the data tied up in the basement, it's time to make it speak!

Naively, we'll start by having a look at the distribution of the closing prices. Maybe the shape of their histogram will look familiar?

In [208]:
plt.hist(prices, bins=40, alpha=0.7, rwidth=0.85, color='#0504aa')
plt.xlabel("Price")
plt.ylabel("Frequency")
plt.title("The distribution of Apples stock closing prices")
plt.show()

No it doesn't look familiar. It looks like some trimodal distribution as have pointed out Kenneth Abbott in the video lecture I metioned earlier.

Less naively, we'll have a look at the distribution of relative (percent) changes of the stock closing price throughout the year.

In [209]:
#Calculating the percent changes
y_diff = -100*np.diff(prices)/prices[1:]
y_diff = y_diff.astype(float)
In [210]:
plt.hist(y_diff, bins=50, rwidth=0.8)
plt.xlabel("Relative changes (%)")
plt.ylabel("Frequency")
plt.title("The distribution of Apples stock closing prices relative changes")
plt.show()

Yeay! this looks familiar! bell shaped and symetric(ish), could it be the Normal distribution?

The most rigourous way to check if the data has a normal distribution would be to conduct a statistical hypothesis test. Which: 1- Is boring 2- Takes time to explain 3- Not pretty

Given these three scientific, well founded arguments, we'll skip hypothesis testing and do my favorite way to judge the goodness of fit of data to a distribution: The qq-plot.

Quick definition: the qq-plot is the graph of the empirical quantiles against the theoritical ones.

With a little background in statistics you'll know that, the better the fit, the closest the qq-plot is to the 45° line.

Here is how you can do it:

In [211]:
# Fit the data to the normal distribution
loc, scale = norm.fit(y_diff)


# Dist plot
plt.figure(figsize=(7, 10))
plt.subplot(2, 1, 1)

m, M = min(y_diff), max(y_diff)
_x = np.linspace(m, M, 500)
_y = norm.pdf(_x, loc=loc, scale=scale)

plt.xlabel("Relative changes (%)")
plt.ylabel("Frequency")
plt.title("Empirical VS Theoritical distribution")

plt.hist(y_diff, bins=50, rwidth=0.8, density='normal')
plt.plot(_x, _y)


# Calculate the empirical quantiles
y_diff_qs = mquantiles(y_diff, np.arange(0.01, 0.99, 0.01))

# And the theoritical ones, thanks LARGE NUMBERS THEOREM
norm_qs   = mquantiles(norm.rvs(loc=loc, scale=scale, size=10000), np.arange(0.01, 0.99, 0.01))

plt.subplot(2, 1, 2)
# qq-plot
plt.scatter(norm_qs, y_diff_qs)
x_lim = min(y_diff_qs), max(y_diff_qs)
plt.plot(x_lim, x_lim, ls="--", c=".5")
plt.xlabel("Theoritical quantiles")
plt.ylabel("Empirical quantiles")
plt.title("QQ-plot against normal distribution")
plt.show()

The fit seems ok, but I am not liking the S-shape of the qq-plot. Thats an indication of a heavy-tailed ditribution (wikipedia is your friend).

The Cauchy distribution is a heavy tailed one. Let's check it out.

In [212]:
# Fit the data to the Cauchy distribution
loc, scale = cauchy.fit(y_diff)

plt.figure(figsize=(7, 10))
plt.subplot(2, 1, 1)

# Dist plot
m, M = min(y_diff), max(y_diff)
_x = np.linspace(m, M, 500)
_y = cauchy.pdf(_x, loc=loc, scale=scale)

plt.xlabel("Relative changes (%)")
plt.ylabel("Frequency")
plt.title("Empirical VS Theoritical distribution")

plt.hist(y_diff, bins=50, rwidth=0.8, density='normal')
plt.plot(_x, _y)


# Calculate the empirical quantiles
y_diff_qs = mquantiles(y_diff, np.arange(0.01, 0.99, 0.01))

# And the theoritical ones, thanks LARGE NUMBERS THEOREM
cauchy_qs   = mquantiles(cauchy.rvs(loc=loc, scale=scale, size=10000), np.arange(0.01, 0.99, 0.01))

plt.subplot(2, 1, 2)
# qq-plot
plt.scatter(cauchy_qs, y_diff_qs)
x_lim = min(y_diff_qs), max(y_diff_qs)
plt.plot(x_lim, x_lim, ls="--", c=".5")
plt.xlabel("Theoritical quantiles")
plt.ylabel("Empirical quantiles")
plt.title("QQ-plot against normal distribution")
plt.show()

Although the pdf of the Cauchy distributions looks like it fits the distribution better, that qq-plot is not quite a straight line.

Another heavy-tailed distribution is the exponential power distribution (also called the generalized normal distribution).

Let's check it out!

In [213]:
# This dustribution has three parameters
beta, loc, scale = gennorm.fit(y_diff)

plt.figure(figsize=(7, 10))
plt.subplot(2, 1, 1)

# Dist plot
m, M = min(y_diff), max(y_diff)
_x = np.linspace(m, M, 500)
_y = gennorm.pdf(_x, beta=beta, loc=loc, scale=scale)

plt.xlabel("Relative changes (%)")
plt.ylabel("Frequency")
plt.title("Empirical VS Theoritical distribution")

plt.hist(y_diff, bins=50, rwidth=0.8, density='normal')
plt.plot(_x, _y)



# Calculate the empirical quantiles
y_diff_qs = mquantiles(y_diff, np.arange(0.01, 0.99, 0.01))

# And the theoritical ones, thanks LARGE NUMBERS THEOREM
gennorm_qs = mquantiles(gennorm.rvs(loc=loc, beta=beta, scale=scale, size=10000), np.arange(0.01, 0.99, 0.01))

plt.subplot(2, 1, 2)
# qq-plot
plt.scatter(gennorm_qs, y_diff_qs)
x_lim = min(y_diff_qs), max(y_diff_qs)
plt.plot(x_lim, x_lim, ls="--", c=".5")
plt.xlabel("Theoritical quantiles")
plt.ylabel("Empirical quantiles")
plt.title("QQ-plot against normal distribution")
plt.show()

Now this is a good fit! The theoritical pdf is close tho the empirical distribution and the qq-plot is pretty much a straight 45° line

The VaR

Now that we have setteled for a distribution for the Apple stock percent changes, we can find the VaR by finding the $1%$ quantile of the distribution.

We can do that like this:

In [214]:
var_distribution = mquantiles(gennorm.rvs(loc=loc, beta=beta, scale=scale, size=1000000),
                              [0.01])
var_distribution
Out[214]:
array([-4.84942823])

The VaR of the Appl Stock is then: $4.99\%$

That means that, $99\%$ of the times, our losses will be no more than $4.99\%$ if we invest solely in Apple stock.

Enter Correlation!

What we saw so far wasn't so hard! The trickiest part was finding the right the distribution.

Finding the VaR was done by generating (pseudo-)random data according to the theoritical distribution and finding the first percentile.

Considering a portfolio of two different companies shares, if we generate the random data for each company we'll end-up with two independent sets of data. Which is generaly not the case.

So the question becomes: How to generate correlated data?

Let's see!

Again I urge you to complete the lecture on VaR I shared with you earlier.

Fetching Microsoft's data

We'll switch from Numpy 2-D arrays to Pandas dataframes which has a great support merging datasets.

In [223]:
msft_data = last_year_data(get_data('MSFT', api_key))

msft_table = pd.DataFrame(msft_data, columns=['date', 'MSFT'])
aapl_table = pd.DataFrame(one_year_data, columns=['date', 'AAPL'])

data = pd.merge(aapl_table, msft_table, on='date', how='inner')
data.sort_values('date', inplace=True)

data[['AAPL', 'MSFT']] = data[['AAPL', 'MSFT']].apply(pd.to_numeric)

Are Apple's closing prices and Microsoft's independent?

The simplest way to investigate such question is to get the correlation matrix of the two time series.

In [224]:
corr_mat = np.corrcoef(data['AAPL'], data['MSFT'])
print(corr_mat)

[[1.         0.95528353]
 [0.95528353 1.        ]]

We can see that the correlation coeficient between the two closing prices is $0.3969$ or $39.69\%$ which indicates that that the two prices are not independent.

You can have a feel for that dependence just by looking at the plots of the prices.

In [225]:
fig1, ax1 = plt.subplots()

ax1.plot(data['date'], data['AAPL'])
ax2 = ax1.twinx()
ax2.plot(data['date'], data['MSFT'], color='r')

ax1.set_ylabel("Apple's price")
ax2.set_ylabel("Microsoft's price")

plt.xlabel("Dates")
plt.title("Apple's an Microsoft's closing prices")
plt.show()

Finding the right distribution for Microsoft's data

We try to fit Microsoft's data to the normal, Cauchy and the exponential power distribution, the later proves out to be the best.

I'll skip the trials and show you just the plots for the exponential power distribution.

In [226]:
# Calculating the percentage changes for Microsoft's closing price
msft_diff = data['MSFT'].pct_change(periods=1).dropna()
In [227]:
# This dustribution has three parameters
beta, loc, scale = gennorm.fit(msft_diff)

plt.figure(figsize=(7, 10))
plt.subplot(2, 1, 1)

# Dist plot
m, M = min(msft_diff), max(msft_diff)
_x = np.linspace(m, M, 500)
_y = gennorm.pdf(_x, beta=beta, loc=loc, scale=scale)

plt.xlabel("Relative changes (%)")
plt.ylabel("Frequency")
plt.title("Empirical VS Theoritical distribution")

plt.hist(msft_diff, bins=50, rwidth=0.8, density='normal')
plt.plot(_x, _y)



# Calculate the empirical quantiles
msft_diff_qs = mquantiles(msft_diff, np.arange(0.01, 0.99, 0.01))

# And the theoritical ones, thanks LARGE NUMBERS THEOREM
gennorm_qs = mquantiles(gennorm.rvs(loc=loc, beta=beta, scale=scale, size=10000), np.arange(0.01, 0.99, 0.01))

plt.subplot(2, 1, 2)
# qq-plot
plt.scatter(gennorm_qs, msft_diff_qs)
x_lim = min(msft_diff_qs), max(msft_diff_qs)
plt.plot(x_lim, x_lim, ls="--", c=".5")
plt.xlabel("Theoritical quantiles")
plt.ylabel("Empirical quantiles")
plt.title("QQ-plot against normal distribution")
plt.show()

Generating correlated data

Here is the tricky part of this tutorial.

Let's fit Apple's an Microsoft's data to their theoritical distribution, generate data according to those distributions and estimate the correlation matrix for the generated data

In [228]:
# Fitting to distribution for Microsoft
beta, loc, scale = gennorm.fit(msft_diff)
fake_msft = gennorm.rvs(loc=loc, beta=beta, scale=scale, size=1000000)

# Fitting to distribution for Apple
aapl_diff = data['AAPL'].pct_change(periods=1).dropna()
beta, loc, scale = gennorm.fit(aapl_diff)
fake_aapl = gennorm.rvs(loc=loc, beta=beta, scale=scale, size=1000000)

np.corrcoef(fake_aapl, fake_msft)
Out[228]:
array([[ 1.0000000e+00, -3.6715976e-04],
       [-3.6715976e-04,  1.0000000e+00]])

See that the generated date are not correlated, which does not reflect the reality that the closing prices for Microsoft and Apple are correlated.

To see how we can generate correlated data from uncorrelated data, let's delve into the definition of the correlation matrix.

Matimatiks

Let $X=[X_1, X_2...X_n]$ be a random vector of dimension $n$.

(In our case, $n=2$, and $X_0$ is Apple's stock random variable and $X_1$ is Microsoft's)

Definition: The correlation matrix, $P_X$, of $X$ is defined as:

$$P_X = \mathrm{E}[XX^{\intercal}]$$

Where $\mathrm{E}$ is the expectation and $X^{\intercal}$ is the transpose of $X$.

If the random variables $X_1, X_2...X_n$ are uncorrelated two-by-two then their correlation matrix is the identity matrix $I$.

The trick

Step 1: Get the eigenvalues $\Lambda=[\lambda_1, \lambda_2...\lambda_3]$ of a correlation matrix $P$ and their eigenvectors $V=[v_1, v_2...v_3]$.

Step 2: Write the correlation matrix in the canonical form (Eigendecomposition) as follows:

$$P = V diag(\Lambda) V^{\intercal} = \left( V \sqrt{diag(\Lambda)} \right) \left(\sqrt{diag(\Lambda)} V \right)^{\intercal} $$

Where $diag(\Lambda)$ is the diagonal matrix made of the elements of $\Lambda$

Step 3: Generate an uncorrelated random vector $Z$. That it's correlation matrix is equal to the identity matrix.

$$P_Z = \mathrm{I}$$

Step 4: Take $X=\left( V \sqrt{diag(\Lambda)} \right)Z$. The correlation matrix of $X$ is $P$. Indeed:

$P_X = \mathrm{E}[XX^{\intercal}] $

$P_X$ $=\mathrm{E}\left[\left( V \sqrt{diag(\Lambda)} \right)Z Z^{\intercal}\left( V \sqrt{diag(\Lambda)} \right)^{\intercal}\right] $

$P_X$ $= \left( V \sqrt{diag(\Lambda)} \right) \mathrm{E}[Z Z^{\intercal}] \left( V \sqrt{diag(\Lambda)} \right)^{\intercal} $

$P_X$ $= \left( V \sqrt{diag(\Lambda)} \right) \mathrm{I} \left( V \sqrt{diag(\Lambda)} \right)^{\intercal} $

$P_X$ $= V diag(\Lambda) V^{\intercal} $

$P_X$ $= P $

Application

Here is the code for the steps above:

In [229]:
# Calculating the correlation matrix for the two stock prices
P = np.corrcoef(data['AAPL'], data['MSFT'])

# Eigendecomposition
w, m = np.linalg.eig(P)
V_sqrtL = np.dot(m, np.diag(np.sqrt(w)))

# Z is the uncorrelated generated data
Z = np.array([fake_aapl, fake_msft])

# X should be correlated
X = np.dot(V_sqrtL, Z)

# Calculating the correlation matrix for X
P_X = np.corrcoef(X[0], X[1])
P_X
Out[229]:
array([[1.        , 0.97354997],
       [0.97354997, 1.        ]])

As you may see, now the data in the vector $X$ are correlated, with a correlation coefficient close to the one estimated for the original data.

In [230]:
P
Out[230]:
array([[1.        , 0.95528353],
       [0.95528353, 1.        ]])

That difference is (I think) due to the sampling error. It is quite hard to generate totally uncorelated data.

The VaR

Now given this precious correlated data that we generated, we can estinate the var for a portfolio made of shares from Apple's stock and other from Microsofts.

Let's sat I invest $40\%$ of my money in Apple and $60\%$ in Microsoft. My position on the market can be sumarized with the followint array:

In [231]:
positions_vector = np.matrix([0.4, 0.6])

#We group the generated data in one matrix
fake_data = np.matrix([fake_aapl, fake_msft])

# Calculating the returns generated by my position
returns = np.array(np.dot(positions_vector, fake_data))

# And now thw VaR
VaR = mquantiles(returns[0], [0.01])
VaR
Out[231]:
array([-0.02635678])

Our port folio has a var of $2.63\%$. Or $99%$ percent of the time, out losses will be no more than $2.22%$. Which is better than owning only one stock! This is an illustration of the diversification principle that states, roughly, that the more your portfolio is diversified the less risk you are exposed to

Our tutorial has come to an end, again, scroll to the top and go watch the video lecture.

Thanks for reading!

back to top