Distinguishing EWMA Volatility Forecasts from Static VaR Fits

# EWMA VaR, code from Quant Risk


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

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

## Answer by Tim Wilding (score 1)

https://quant.stackexchange.com/a/37773

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.