EWMA VaR, código de Quant Risk

Question

Actualmente estoy tratando de mejorar mi python y la comprensión de VaR. He estado usando el sitio web Quant Risk: http://www.quantatrisk.com/2015/12/02/student-t-distributed-linear-value-at-risk/ . Concretamente el gran trozo de código:

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

Me gustaría ser capaz de implementar una superposición que calcula un EWMA de la varianza, pero no estoy seguro de cómo ir sobre it.i sabe que probablemente necesitamos algo a lo largo de las líneas de:

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

Pero no estoy seguro de cómo introducirlo. ¿Alguien puede ayudar?

Answer 1

Creo que la volatilidad ponderada exponencialmente es un modelo de volatilidad ligeramente diferente al de su enlace web y sería difícil hacer simples cambios en el código de ese enlace. La volatilidad ponderada exponencialmente es una previsión dinámica de la volatilidad que cambia cada período ( $\sigma_t$ = volatilidad en el período $t$ ), mientras que el código python de tu enlace parece ocuparse de lo bien que un único modelo estático se ajusta a todos los períodos. El código python se ajusta a $\sigma$ para una distribución Normal que es una estimación de la volatilidad para todo el periodo.

Si utiliza una vida media bastante corta, podría comparar la volatilidad prevista para los últimos períodos y determinar si $r_t/\sigma_t$ es más coherente con un $\mathcal{N}(0,1)$ distribución que $r_t/\sigma$ . Así, por ejemplo, si su vida media es de aproximadamente 13 días, entonces podría utilizar los primeros 100 días de rendimientos para calcular una volatilidad para el período 100. La vida media de 13 días significa que el efecto de los 100 $^{th}$ en la volatilidad de EW es bastante pequeña.