Aplicando correctamente GARCH en Python

Question

Problema: Uso correcto de GARCH(1,1)
Objetivo de la investigación: Pronóstico de volatilidad/varianza.
Herramientas utilizadas: Python
Instrumento: SPX (precios de cierre ajustados específicamente)
Material de referencia: Sobre la Estimación de Modelos GARCH con una Aplicación a los Precios de Acciones de Nordea (Chao Li, 2007)

Nota: He revisado casi todas las publicaciones de Quant.SE que discuten GARCH, pero no he visto ninguna de ellas con la naturaleza aproximada de lo que estoy preguntando.

Descripción:

Estoy realizando una investigación personal sobre la volatilidad histórica de los precios de las acciones cuando me encontré con el concepto de pronóstico de volatilidad. Naturalmente, esto llamó mi interés. Aterricé en el material de referencia mencionado anteriormente. Después de un par de lecturas, decidí ver si mi comprensión de la aplicación de GARCH(1,1) es correcta y si puedo obtener alguna idea sobre cómo usarlo al investigar la volatilidad histórica.

Para esto, utilizaré los precios de SPX y las bibliotecas bt, pandas y arch en Python.

Mis pasos de programación iniciales son los siguientes.

In [1]: import pandas as pd
   ...: import bt
   ...: import arch
   ...: 

In [2]: df = bt.get('SPX', start='1990-01-01')

In [3]: df.head()
Out[3]: 
               spx
Date              
1990-01-02  359.69
1990-01-03  358.76
1990-01-04  355.67
1990-01-05  352.20
1990-01-08  353.79

In [4]: df['pct_change'] = df['spx'].pct_change().dropna()
   ...: df['stdev21'] = pd.rolling_std(df['pct_change'], 21)
   ...: df['hvol21'] = df['stdev21']*(252**0.5) # Anualizar.
   ...: df['variance'] = df['hvol21']**2
   ...: df = df.dropna() # Eliminar filas con celdas en blanco.
   ...: df.head()

Out[4]: 
               spx  pct_change   stdev21    hvol21  variance
Date                                                        
1990-03-02  335.54    0.008415  0.007304  0.115946  0.013443
1990-03-05  333.74   -0.005364  0.007425  0.117863  0.013892
1990-03-06  337.93    0.012555  0.007770  0.123344  0.015214
1990-03-07  336.95   -0.002900  0.007804  0.123889  0.015348
1990-03-08  340.27    0.009853  0.007855  0.124689  0.015547

Lo anterior es bastante simple. Ahora, usaré la función GARCH proporcionada por el módulo de Python arch para obtener omega, beta y alpha.

In [5]: returns = df['pct_change'] * 100
   ...: am = arch.arch_model(returns)
   ...: res = am.fit(iter=5)
   ...: res.params
Iteración:   5,   Conteo Func.:  39,   LLF Neg.: 8447.41751792
Iteración:  10,   Conteo Func.:  74,   LLF Neg.: 8443.32521758
Optimización terminada con éxito.    (Salida modo 0)
            Valor actual función: 8443.31731767
            Iteraciones: 12
            Evaluaciones función: 86
            Evaluaciones de gradiente: 12
Out[5]: 
mu          0.058224
omega       0.011511
alpha[1]    0.079411
beta[1]     0.911240
Nombre: parámetros, dtype: float64

De acuerdo con la fórmula $\sigma_t^2 = \omega + \alpha_1{a^2}_{t-1} + \beta_1{\sigma^2}_{t-1}$, ejecuto el siguiente código.

In [6]: df['C'] = res.params['omega']
   ...: df['B'] = df['variance'] * res.params['beta[1]']
   ...: df['A'] = (df['pct_change']**2) * res.params['alpha[1]']
   ...: df['forecast_var'] = df.loc[:,'C':'A'].sum(axis=1)
   ...: df['forecast_vol'] = df['forecast_var']**0.5
   ...: df.head()
   ...: 
Out[6]: 
               spx  pct_change   stdev21    hvol21  variance         C  \
Date                                                                     
1990-03-02  335.54    0.008415  0.007304  0.115946  0.013443  0.011511   
1990-03-05  333.74   -0.005364  0.007425  0.117863  0.013892  0.011511   
1990-03-06  337.93    0.012555  0.007770  0.123344  0.015214  0.011511   
1990-03-07  336.95   -0.002900  0.007804  0.123889  0.015348  0.011511   
1990-03-08  340.27    0.009853  0.007855  0.124689  0.015547  0.011511   

                   B         A  forecast_var  forecast_vol  
Date                                                        
1990-03-02  0.012250  0.000006      0.023767      0.154165  
1990-03-05  0.012659  0.000002      0.024172      0.155474  
1990-03-06  0.013863  0.000013      0.025387      0.159333  
1990-03-07  0.013986  0.000001      0.025498      0.159681  
1990-03-08  0.014167  0.000008      0.025686      0.160269

Ahora llego a un valor de volatilidad pronosticado. Mis preguntas son las siguientes:

1. En el bloque de código [5], se utilizó esto: returns = df['pct_change'] * 100. Lo tomé arbitrariamente tal como está, ya que es la forma en que he visto que se utilizan los retornos en los cálculos de GARCH. Sin embargo, en los cálculos iniciales de la varianza, nunca necesité multiplicar la columna pct_change por 100. _¿Debo dejarlo así o debo proporcionar la columna pct_change tal como está a la función?_

2. Si proporciono la columna pct_change a la función tal como está, los valores de omega, beta y alpha se vuelven considerablemente más pequeños, lo que reduce los valores de forecast_var y forecast_vol en un orden de magnitud. Obviamente, el cambio en la precisión causa una comparación problemática entre mis valores históricos de volatilidad y mis valores pronosticados de volatilidad. ¿Cómo puedo resolver esto?

3. ¿Debería estar usando un enfoque completamente diferente?

Por favor, señale también cualquier error evidente que haya pasado por alto o cualquier lógica faltante

Answer 1

No importa si usas *100 o simplemente pct_change, siempre y cuando seas consistente. Sin embargo, en la práctica, debido a las inestabilidades numéricas de punto flotante subyacentes en los algoritmos de optimización subyacentes/tolerancias predeterminadas utilizadas en scipy/arch, tener los rendimientos expresados en %, es decir, multiplicados por 100, tendrá una mejor probabilidad de converger durante el ajuste del modelo.

Los errores comienzan en In[6]. Una vez que el modelo está ajustado, puedes obtener las volatilidades condicionales previstas en res.conditional_volatility, las cuales necesitas anualizar, es decir, multiplicar por sqrt(252).

Ten en cuenta que en la fórmula de GARCH, a(t-1) es el residuo del modelo, el cual puedes encontrar en res.residual. No es el pct_change**2. Para llegar al residuo a partir de pct_change, tienes que trabajar hacia atrás en las ecuaciones. Pero esto ya se ha hecho por ti en res.resid

De todas formas, la fórmula correcta en In[6] para forecast_vol debería ser:

0.01 * np.sqrt(res.params['omega'] + res.params['alpha[1]'] * res.resid**2 + res.conditional_volatility**2 * res.params['beta[1]'])

^^ Aquí estás pronosticando la volatilidad instantánea para el PRÓXIMO día. Multiplicando por 0.01 para compensar la multiplicación anterior del rendimiento diario (pct_change) por 100

Espero que esto aclare un poco las cosas.

Answer 2

Creo que en la mejor práctica deberías multiplicar por 100 (o incluso por un número mayor, siempre que los valores estén entre 0 - 1000). La razón es para la convergencia del optimizador como mencionaron anteriormente. La volatilidad por supuesto se verá afectada (escalada aproximadamente 100 veces), así que después de calcular las volatilidades, debes reescalarlas a 1/100.

Saludos cordiales,