# EWMA VaR, code from Quant Risk

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',
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]")

# 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)
else:
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)
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()

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


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]
EWMAstdev[-i-1]=math.sqrt(sum(Weights*sqrdData))


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

• Isn't this quite well treated here: stats.stackexchange.com/questions/111851/… – Ric Feb 16 '18 at 14:07
• There is some Python code there as well and psuedo code in a mentioned link – Ric Feb 16 '18 at 14:08

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.