I am currently attempting to improve my python and understanding of VaR. I have been using the website Quant Risk: http://www.quantatrisk.com/2015/12/02/student-t-distributed-linear-value-at-risk/. Specifically the great bit of code:

# 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.ylabel("IBM Adj Close [USD]")
plt.subplot(2, 1, 2)
plt.plot(ret, color=(.6, .6, .6))
plt.ylabel("Daily Returns")
plt.xlabel("Time Period 2010-12-01 to 2015-12-01 [days]")

# if ret ~ N(0,1), the expected skewness = 0, kurtosis = 3
# ret = np.random.randn(10000000)

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)
 print("Accept H_0 at %.5f level\n" % beta)

# 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)

 # 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)

 # 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))

  # plt.close("all")
   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.plot(x, pdf, 'b', label="Normal PDF fit")
    plt.plot(x, pdf2, 'g', label="Student t PDF fit")
    plt.xlim([-0.2, 0.1])
    plt.ylim([0, 50])
    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.plot(x, pdf, 'b')
    plt.plot(x, pdf2, 'g')
    # 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])

I would really like to be able to implement an overlay which calculates an EWMA of the variance, but I am not sure how to go about it.i know that we probably need something along the lines of:

  EWMAstdev = np.empty([len(ret)-Period_Interval,])
    stndrData = pd.Series(index=ret.index)

    # For efficiency here we should square returns first so the loop does not do it repeadetly 
    sqrdReturns = ret**2

    # Computations here happen in different times, because we first need all the EWMAstdev
    # First get the stdev according to the EWMA
    for i in range(0,len(Returns)-Period_Interval):
        if i == 0: sqrdData = sqrdReturns[-(Period_Interval):]
        else: sqrdData = sqrdReturns[-(Period_Interval+i):-i]                 

But am unsure how to slide it in? Can anyone help?


I think the exponentially-weighted volatility is a slightly different volatility model to your website link and it would be tough to make simple changes to the code at that link. The exponentially-weighted volatility is a dynamic forecast of the volatility that changes each period ($\sigma_t$ = volatility at period $t$), whereas the python code at your link appears to deal with how well a single static model fits all of the periods. The python code fits $\sigma$ for a Normal distribution which is an estimate of the volatility for the entire period.

If you are using a fairly short half-life, you could compare the forecast volatility for the latter periods and determine whether $r_t/\sigma_t$ is more consistent with a $\mathcal{N}(0,1)$ distribution than $r_t/\sigma$. So, for example, if your half-life is approx. 13 days, then you could use the first 100 days of returns to calculate a volatility for the 100th period.The 13-day half-life means that the effect of the 100$^{th}$ period on the EW volatility is quite small.


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.