Procesos estocásticos

Como hemos mencionado anteriormente, nosotros no nos centraremos en el estucio clásico de series temporales, sino que consideraremos nuestros datos como fenómenos que evolucionan en el tiempo. Para ello, los modelos adecuados son los llamados procesos estocásticos.

Definición Un proceso estocástico es una colección de variables aleatorias indexadas por un conjunto \(\mathbb{T}\), que habitualmente representa el tiempo: \({X_t, t \in \mathbb{T}}\). El proceso puede ser a tiempo discreto, es decir \(\mathbb{T}\) es \(\mathbb{N}\) ó \(\mathbb{Z}\), o por otra parte a tiempo continuo, dónde \(\mathbb{T}\) es \([0,\infty)\) ó \(\mathbb{R}\).

Nosotros a partir de ahora veremos las series temporales como una realización de un proceso estocástico.

Procesos de segundo orden

Definición: Un proceso \({X_t, t\, \in \mathbb{N} \text{ ó } \mathbb{Z}}\) es de segundo orden si se cumple que

\[E(X^2_t)<\infty\, \forall t\]

Lo cual garantiza que var\((X_t)\) está bien definida y es finita, ya que se cumple que:

\[|E(X_t)|\leq E(|X_t|)\leq (E(x^2_t))^{1/2}<\infty\]

y por tanto

\[\text{Var}(X_t):=E(x^2_t)-E(X_t)^2<\infty\]

A partir de ahora, cualquier proceso estocástico será para nosotros de segundo orden.

Procesos estacionarios

Cada serie temporal tiene su particular forma, cosa que se puede observar al hacer su respectiva gráfica. Necesitamos que los procesos estocásticos sean la representación matemática de algún tipo de regularidad, para que puedan ser utilizados como modelos de nuestras series.

Definición: Un proceso \({X_t, t \in \mathbb{Z}}\) es estrictamente estacionario si para cualquier conjunto de variables incluido en este, éstas presentan la misma ley que las que obtenemos desplazándolas una determinada longitud, es decir, que para cualquier \(t_1, . . . , t_n\) y \(l\), los vectores

\[x_{t_1},...,x_{t_n} \]

y

\[x_{t_{1+l}},...,x_{t_{n+l}} \]

presenten la misma ley.

Nota: Para \(n = 1\) implica que todas las \(X_t\) tienen la misma ley.

Definición: Un conjunto es estacionario (o estacionario en sentido débil) si todas sus variables aleatorias tienen la misma esperanza y la covarianza que presentan solamente depende de la distancia que las separa, es decir:

Nota: Para \(l = 0\) tenemos que \(\gamma(0) = \text{Var}(X_t), \forall t \geq 1\) y también se cumple que \(\gamma(l) =\gamma(-l)\).

A partir de la función \(\gamma(l)\), llamada función de autocovarianza, obtenemos la función

\[\rho(l):=\frac{\gamma(l)}{\gamma(0)}\]

llamada función de autocorrelación.

Por lo tanto, para poder analizar nuestras series, necesitaremos hacer un estudio previo y aplicarles las transformaciones que sean necesarias para que presenten las propiedades que acabamos de definir.

Ejemplos de procesos estocásticos

Ruido IID

Un proceso \({X_t , t \geq 1}\) es un ruido i.i.d si las variables son independientes e idénticamente distribuidas, y con media \(\mu\) y desviación típica \(\sigma\). En este caso, la función de autocorrelación toma los valores \(\rho(l) = 1\) para \(l = 0\) y \(\rho(l) = 0, ∀l > 0\).

Ruido blanco

Un proceso \({X_t, t ≥ 1}\) es un ruido blanco si las variables tienen media \(\mu\), varianza \(\sigma^2\) y están incorreladas, es decir, no tienen relación lineal \((Cov(X_t, X_k) = 0, ∀t, ∀k)\), pero podrían presentar algún tipo de relación no lineal.

Ruido blanco

Nota: Un ruido i.i.d es un caso particular de ruido blanco.

Ruido blanco Gaussiano

Un proceso \({X_t, t ≥ 1}\) es un ruido blanco gaussiano si a parte de ser blanco, todo vector \((X_{t_1}, . . . , X_{t_n})\) tiene ley normal.

Un ruido blanco gaussiano es un caso particular de ruido i.i.d, ya que en este caso la incorrelación implica independencia.

Modelos ARCH

La heterocedasticidad condicional auto-regresiva (en inglés, Autoregressive conditional heteroscedasticity; ARCH) es la condición de que hay uno o más puntos de datos en una serie para los cuales la varianza del término de error actual o innovación es una función de los tamaños reales de los términos de error de los períodos de tiempo anteriores: se relaciona con los cuadrados de las innovaciones anteriores. En la econometría, los modelos ARCH se utilizan para caracterizar y modelar series temporales. Una variedad de otras siglas se aplican a las estructuras particulares que tienen una base similar.

Los modelos de ARCH se emplean comúnmente en el modelado de series de tiempo financieras que presentan agrupaciones de volatilidad variables en el tiempo, es decir, períodos de oscilaciones entremezclados con períodos de relativa calma. A veces se considera que los modelos de tipo ARCH pertenecen a la familia de modelos de volatilidad estocástica , aunque esto es estrictamente incorrecto ya que en el tiempo t la volatilidad es completamente predeterminada (determinista) dados los valores anteriores.

Para ilustrar el significado de los momentos condicionales y no condicionales de un proceso estocástico consideraremos el caso más simple de un proceso AR(1) estacionario.

\[y_t=\phi_1 y_{t-1}+u_t\]

donde \(0<\phi_1<1\), siendo \(u_t\sim iid(0,\sigma^2)\)

Los momentos condicionales de ese proceso dependen de los valores pasados de \(y_t\), dado que el proceso es autorregresivo de primer orden la información del pasado se limita a la existente en el período inmediatamente anterior \(y_{t-1}\).

Por lo tanto, al tomar la esperanza matemática condicional del proceso vamos a obtener:

\[E_{t-1}[y_t]=E[y_t|y_{t-1},y_{t-2},\cdots]=E[y_{t}|y_{t-1}]=E[\phi_1y_{t-1}+u_t]=\phi_1y_{t-1}\]

Por su parte la varianza condicional sería la siguiente expectativa:

\[\text{var}_{t-1}[y_{t}]=E[(y_y-E(y_t))^2|y_{t-1},y_{t-2},\cdots]=E[(y_1-\phi_1y_{t-1})^2|y_{t-1}]=E[u^2_t]=\sigma^2\]

Los momentos no condicionales ya no dependen de las realizaciones en el tiempo de \(y_t\), y toda la información procede del término de perturbación aleatoria.

Para la media no condicional tenemos:

\[\begin{align*} y_t-\phi_1y_{t-1}&=u_t\\ y_t(1-\phi_1 L)&=u_t\\ y_t&=(1-\phi_1 L)^{-1}u_t\\ y_t&=(1+\phi_1 L+\phi_1^2L^2+\cdots)u_t\\ y_t&=u_t+\phi_1u_{t-1}+\phi_1^2u_{t-2}+\cdots\\ E[y_t]&=0 \end{align*}\]

Para la varianza no condicional:

\[\begin{align*} \text{var}(y_t)&=E(u_t+\phi_1u_{t-1}+\phi_1^2u_{t-2}+\cdots)^2\\ \text{var}(y_t)&=E(u^2_t+\phi_1^2u_{t-1}^2+\phi_1^4u_{t-2}^2+\cdots)\\ \text{var}(y_t)&=\sigma^2+\phi_1^2\sigma^2+\phi_1^4\sigma^2+\cdots\\ \text{var}(y_t)&=\sigma^2(1+\phi_1^2+\phi_1^4+\cdots)\\ \text{var}(y_t)&=\frac{\sigma^2}{1-\phi_1^2} \end{align*}\]

Los productos cruzados en la operación previa son nulos en virtud de que los términos de perturbación aleatoria son independientes en el tiempo.

En su trabajo fundamental sobre los procesos ARCH, Engle (1982) plantea que los intervalos de predicción de un proceso estocástico podrían mejorar si pudiéramos utilizar más información para la predicción de la varianza. El modelo propuesto por Engle es el siguiente:

\[\begin{align*} y_t&=u_t\\ u_t&=v_t h_t^{1/2}\\ h_t&=\alpha_0+\alpha_1u^2_{t-1} \end{align*}\]

siendo \(v_t\sim N(0,1)\) y \(y_t|\Theta_{t-1}\sim N(0,h_t)\) siendo \(\Theta_t\) el conjunto de información disponible en \(t\).

Donde la primera ecuación es la de la media, la segunda es el término de perturbación y la tercera es \(h_t\) la varianza de \(u_t\) condicional a la información disponible en el período \(t\).

En su trabajo original Engle supone que la ecuación de la media \(y_t\) podría tener procesos más complejos involucrando variables explicatorias de la forma \(y_t=x_t\beta\).

De las tres ecuaciones previas podemos obtener la varianza condicional de la siguiente manera.

\[u_t=v_t\sqrt{\alpha_0+\alpha_1u^2_{t-1}}\]

Los momentos no condicionales de esta varianza son los siguientes.

La esperanza matemática:

\[E[u_t]=E\left[v_t\sqrt{\alpha_0+\alpha_1u^2_{t-1}}\right]=0\]

La varianza

\[\text{var}[u_t]=E[v_t^2(\alpha_0+\alpha_1u^2_{t-1})]=\alpha_0+\alpha_1\text{var}[u_t]=\frac{\alpha_0}{1-\alpha_1}\]

Claramente tenemos un proceso estacionario en virtud de que \(0<\alpha_1<1\).

Si ahora obtenemos los momentos condicionales:

Media condicional:

\[E_{t-1}[u_t]=E_{t-1}\left[v_t\sqrt{\alpha_0+\alpha_1u^2_{t-1}}\right]=0\]

La varianza condicional:

\[\text{var}_{t-1}[u_t]=E_{t-1}[v_t^2(\alpha_0+\alpha_1u^2_{t-1})]=\alpha_0+\alpha_1u^2_{t-1}\]

Destaca en este resultado que la varianza condicional sigue un proceso en el cual la media condicional depende del valor previo del término de perturbación al cuadrado. La varianza condicional será un valor positivo y será estable si se cumple la condición de estacionariedad \(0<\alpha_1<1\) y la de su deriva positiva \(\alpha_0>0\).

El modelo puede generalizarse a procesos AR de mayor orden, si lo generalizamos tendremos un ARCH\((p)\):

\[u=v_t\sqrt{\alpha_0+\alpha_1u^2_{t-1}+\cdots+\alpha_pu^2_{t-p}}\]

La condición de estabilidad es que \(\displaystyle \sum_{i=1}^p\alpha_i<1\) y las de no egatividad \(\alpha_0>0\), \(\alpha_i\geq 0\).

Una representación equivalente para la ecuación de la varianza condicional y que es la que se emplea comúnmente en los libros de texto es la siguiente:

\[h_t^2=\alpha_0+\alpha_1u^2_{t-1}+\cdots+\alpha_pu^2_{t-p}\]

En la ecuación ARCH resulta claro que la varianza es heterocedastica y depende del cuadrado de los choques aleatorios pasados y, en consecuencia, tendrán el mismo impacto en ella tanto los choques positivos como negativos.

library(tidyquant)
library(TTR)
library(ggplot2)
library(knitr)
library(nortest)
library(MASS)
library(stats)
library(moments)
library(fGarch)
library(fBasics)


AAPL <- tq_get("AAPL", get="stock.prices", from= "2021-07-01", to= "2021-07-31") 
AAPL <- data.frame(AAPL) 
precio<-AAPL$close
n<-length(precio)
rendimiento<-c()
for (i in 1:n) {
  rendimiento[i]<-(precio[i+1]-precio[i])/(precio[i])
}
rendimiento<-na.omit(rendimiento)



modelo<-garchFit(~garch(1,0),data=rendimiento)
## 
## Series Initialization:
##  ARMA Model:                arma
##  Formula Mean:              ~ arma(0, 0)
##  GARCH Model:               garch
##  Formula Variance:          ~ garch(1, 0)
##  ARMA Order:                0 0
##  Max ARMA Order:            0
##  GARCH Order:               1 0
##  Max GARCH Order:           1
##  Maximum Order:             1
##  Conditional Dist:          norm
##  h.start:                   2
##  llh.start:                 1
##  Length of Series:          20
##  Recursion Init:            mci
##  Series Scale:              0.01437149
## 
## Parameter Initialization:
##  Initial Parameters:          $params
##  Limits of Transformations:   $U, $V
##  Which Parameters are Fixed?  $includes
##  Parameter Matrix:
##                      U          V    params includes
##     mu     -2.18314636   2.183146 0.2183146     TRUE
##     omega   0.00000100 100.000000 0.1000000     TRUE
##     alpha1  0.00000001   1.000000 0.1000000     TRUE
##     gamma1 -0.99999999   1.000000 0.1000000    FALSE
##     delta   0.00000000   2.000000 2.0000000    FALSE
##     skew    0.10000000  10.000000 1.0000000    FALSE
##     shape   1.00000000  10.000000 4.0000000    FALSE
##  Index List of Parameters to be Optimized:
##     mu  omega alpha1 
##      1      2      3 
##  Persistence:                  0.1 
## 
## 
## --- START OF TRACE ---
## Selected Algorithm: nlminb 
## 
## R coded nlminb Solver: 
## 
##   0:     53.634740: 0.218315 0.100000 0.100000
##   1:     28.443939: 0.219229  1.00633 0.522557
##   2:     28.152616: 0.246720  1.20664 0.0820831
##   3:     28.136499: 0.247282  1.16020 0.164411
##   4:     28.091898: 0.248239  1.13900 0.168434
##   5:     27.697665: 0.261373 0.917179 0.229167
##   6:     27.555195: 0.261564 0.836282 0.183063
##   7:     27.482053: 0.271755 0.665759 0.241589
##   8:     27.475333: 0.279710 0.685739 0.278002
##   9:     27.473502: 0.278457 0.686049 0.261599
##  10:     27.473478: 0.278749 0.682866 0.262526
##  11:     27.473474: 0.278660 0.683789 0.261994
##  12:     27.473474: 0.278663 0.683774 0.262026
##  13:     27.473474: 0.278663 0.683774 0.262024
## 
## Final Estimate of the Negative LLH:
##  LLH:  -57.3767    norm LLH:  -2.868835 
##           mu        omega       alpha1 
## 0.0040048001 0.0001412266 0.2620243186 
## 
## R-optimhess Difference Approximated Hessian Matrix:
##                  mu         omega       alpha1
## mu     -104897.2540    -314913.25    158.79432
## omega  -314913.2454 -316341548.83 -37657.93758
## alpha1     158.7943     -37657.94    -13.16173
## attr(,"time")
## Time difference of 0.004029989 secs
## 
## --- END OF TRACE ---
## 
## 
## Time to Estimate Parameters:
##  Time difference of 0.07790112 secs

ahora observemos los coeficientes del modelo

coef(modelo)
##           mu        omega       alpha1 
## 0.0040048001 0.0001412266 0.2620243186

luego el modelo tiene la siguiente forma

\[\begin{align*} h_t^2&=\alpha_0+\alpha_1u^2_{t-1}\\ &=0.0001412266+ 0.2620243186 u^2_{t-1} \end{align*}\]

con dicho modelo podemos hacer una predicción de como puede variar el rendimiento de la acción en 5 momentos después

predict(modelo,n.ahead=5)

veamoslo en un gráfico

predict(modelo,n.ahead=5,plot=TRUE)

podemos revisar otros modelos ARCH\((p)\), tomemos como ejemplo un modelo ARCH(2) y un ARCH(3)

ARCH(2)

modelo2<-garchFit(~garch(2,0),data=rendimiento)
## 
## Series Initialization:
##  ARMA Model:                arma
##  Formula Mean:              ~ arma(0, 0)
##  GARCH Model:               garch
##  Formula Variance:          ~ garch(2, 0)
##  ARMA Order:                0 0
##  Max ARMA Order:            0
##  GARCH Order:               2 0
##  Max GARCH Order:           2
##  Maximum Order:             2
##  Conditional Dist:          norm
##  h.start:                   3
##  llh.start:                 1
##  Length of Series:          20
##  Recursion Init:            mci
##  Series Scale:              0.01437149
## 
## Parameter Initialization:
##  Initial Parameters:          $params
##  Limits of Transformations:   $U, $V
##  Which Parameters are Fixed?  $includes
##  Parameter Matrix:
##                      U          V    params includes
##     mu     -2.18314636   2.183146 0.2183146     TRUE
##     omega   0.00000100 100.000000 0.1000000     TRUE
##     alpha1  0.00000001   1.000000 0.0500000     TRUE
##     alpha2  0.00000001   1.000000 0.0500000     TRUE
##     gamma1 -0.99999999   1.000000 0.1000000    FALSE
##     gamma2 -0.99999999   1.000000 0.1000000    FALSE
##     delta   0.00000000   2.000000 2.0000000    FALSE
##     skew    0.10000000  10.000000 1.0000000    FALSE
##     shape   1.00000000  10.000000 4.0000000    FALSE
##  Index List of Parameters to be Optimized:
##     mu  omega alpha1 alpha2 
##      1      2      3      4 
##  Persistence:                  0.1 
## 
## 
## --- START OF TRACE ---
## Selected Algorithm: nlminb 
## 
## R coded nlminb Solver: 
## 
##   0:     53.913292: 0.218315 0.100000 0.0500000 0.0500000
##   1:     29.549964: 0.219251 0.851396 0.595053 0.421891
##   2:     28.657955: 0.238361  1.02076 0.627672 1.00000e-08
##   3:     28.583242: 0.250071  1.13956 0.446674 1.00000e-08
##   4:     28.541189: 0.267030  1.38944 0.0324276 1.00000e-08
##   5:     28.446055: 0.262291  1.20618 0.288741 1.00000e-08
##   6:     28.341345: 0.263408  1.15018 0.311855 1.00000e-08
##   7:     27.607108: 0.279425 0.644938 0.466432 1.00000e-08
##   8:     27.483107: 0.272793 0.662584 0.224100 1.00000e-08
##   9:     27.470938: 0.273515 0.713212 0.280561 1.00000e-08
##  10:     27.461805: 0.275519 0.668416 0.285442 1.00000e-08
##  11:     27.461128: 0.275092 0.688324 0.253706 1.00000e-08
##  12:     27.460360: 0.275293 0.681762 0.267514 1.00000e-08
##  13:     27.460357: 0.275357 0.681450 0.267138 1.00000e-08
##  14:     27.460357: 0.275367 0.681502 0.267110 1.00000e-08
##  15:     27.460357: 0.275370 0.681501 0.267112 1.00000e-08
## 
## Final Estimate of the Negative LLH:
##  LLH:  -57.38982    norm LLH:  -2.869491 
##           mu        omega       alpha1       alpha2 
## 0.0039574800 0.0001407569 0.2671119120 0.0000000100 
## 
## R-optimhess Difference Approximated Hessian Matrix:
##                  mu         omega        alpha1        alpha2
## mu     -109295.5440    -387099.96    116.902443   -138.231992
## omega  -387099.9622 -317762870.77 -37126.620643 -34668.489288
## alpha1     116.9024     -37126.62    -12.791440     -3.714645
## alpha2    -138.2320     -34668.49     -3.714645     -5.585590
## attr(,"time")
## Time difference of 0.006150961 secs
## 
## --- END OF TRACE ---
## 
## 
## Time to Estimate Parameters:
##  Time difference of 0.02767301 secs
## Warning: Using formula(x) is deprecated when x is a character vector of length > 1.
##   Consider formula(paste(x, collapse = " ")) instead.

cuyos coeficientes

coef(modelo2)
##           mu        omega       alpha1       alpha2 
## 0.0039574800 0.0001407569 0.2671119120 0.0000000100

cuyo modelo tendrá la siguiente forma

\[\begin{align*} h_t^2&=\alpha_0+\alpha_1u^2_{t-1}+\alpha_2u^2_{t-2}\\ &=0.0001407569+ 0.2671119120 u^2_{t-1}0.0000000100u^2_{t-2} \end{align*}\]

veamoslo en un su respectivo gráfico en 5 tiempos después

predict(modelo2,n.ahead=5,plot=TRUE)

ahora hagamos un ARCH(3)

modelo3<-garchFit(~garch(3,0),data=rendimiento)
## 
## Series Initialization:
##  ARMA Model:                arma
##  Formula Mean:              ~ arma(0, 0)
##  GARCH Model:               garch
##  Formula Variance:          ~ garch(3, 0)
##  ARMA Order:                0 0
##  Max ARMA Order:            0
##  GARCH Order:               3 0
##  Max GARCH Order:           3
##  Maximum Order:             3
##  Conditional Dist:          norm
##  h.start:                   4
##  llh.start:                 1
##  Length of Series:          20
##  Recursion Init:            mci
##  Series Scale:              0.01437149
## 
## Parameter Initialization:
##  Initial Parameters:          $params
##  Limits of Transformations:   $U, $V
##  Which Parameters are Fixed?  $includes
##  Parameter Matrix:
##                      U          V     params includes
##     mu     -2.18314636   2.183146 0.21831464     TRUE
##     omega   0.00000100 100.000000 0.10000000     TRUE
##     alpha1  0.00000001   1.000000 0.03333333     TRUE
##     alpha2  0.00000001   1.000000 0.03333333     TRUE
##     alpha3  0.00000001   1.000000 0.03333333     TRUE
##     gamma1 -0.99999999   1.000000 0.10000000    FALSE
##     gamma2 -0.99999999   1.000000 0.10000000    FALSE
##     gamma3 -0.99999999   1.000000 0.10000000    FALSE
##     delta   0.00000000   2.000000 2.00000000    FALSE
##     skew    0.10000000  10.000000 1.00000000    FALSE
##     shape   1.00000000  10.000000 4.00000000    FALSE
##  Index List of Parameters to be Optimized:
##     mu  omega alpha1 alpha2 alpha3 
##      1      2      3      4      5 
##  Persistence:                  0.1 
## 
## 
## --- START OF TRACE ---
## Selected Algorithm: nlminb 
## 
## R coded nlminb Solver: 
## 
##   0:     52.429809: 0.218315 0.100000 0.0333333 0.0333333 0.0333333
##   1:     30.395343: 0.217365 0.731015 0.537153 0.389978 0.503196
##   2:     29.760211: 0.195448  1.01496 0.558788 0.0119617 0.376980
##   3:     29.244896: 0.174659  1.47023 0.302273 1.00000e-08 1.00000e-08
##   4:     29.155552: 0.203039  1.53633 0.180537 1.00000e-08 1.00000e-08
##   5:     28.641369: 0.296252  1.43001 1.00000e-08 1.00000e-08 1.00000e-08
##   6:     28.021278: 0.322309  1.05482 1.00000e-08 1.00000e-08 1.00000e-08
##   7:     27.518165: 0.327025 0.719609 0.205358 1.00000e-08 1.00000e-08
##   8:     27.457419: 0.282959 0.679106 0.259131 1.00000e-08 1.00000e-08
##   9:     27.455516: 0.271914 0.679761 0.266316 1.00000e-08 1.00000e-08
##  10:     27.455476: 0.271858 0.676161 0.264568 1.00000e-08 1.00000e-08
##  11:     27.455442: 0.271261 0.676023 0.267078 1.00000e-08 1.00000e-08
##  12:     27.455441: 0.271526 0.676377 0.266679 1.00000e-08 1.00000e-08
##  13:     27.455441: 0.271516 0.676302 0.266706 1.00000e-08 1.00000e-08
##  14:     27.455441: 0.271519 0.676300 0.266712 1.00000e-08 1.00000e-08
## 
## Final Estimate of the Negative LLH:
##  LLH:  -57.39474    norm LLH:  -2.869737 
##           mu        omega       alpha1       alpha2       alpha3 
## 0.0039021267 0.0001396828 0.2667119485 0.0000000100 0.0000000100 
## 
## R-optimhess Difference Approximated Hessian Matrix:
##                   mu         omega        alpha1        alpha2        alpha3
## mu     -114451.29313    -455593.49     87.967831   -120.777053 -1.236204e+03
## omega  -455593.48672 -318591907.78 -38154.127899 -33221.617445 -2.782825e+04
## alpha1      87.96783     -38154.13    -13.229695     -3.728124 -1.054499e+01
## alpha2    -120.77705     -33221.62     -3.728124     -5.154648  3.580310e+00
## alpha3   -1236.20426     -27828.25    -10.544993      3.580310 -2.878472e-01
## attr(,"time")
## Time difference of 0.01704001 secs
## 
## --- END OF TRACE ---
## Warning in sqrt(diag(fit$cvar)): NaNs produced
## 
## Time to Estimate Parameters:
##  Time difference of 0.04763079 secs
## Warning: Using formula(x) is deprecated when x is a character vector of length > 1.
##   Consider formula(paste(x, collapse = " ")) instead.

con sus respectivos coeficientes

coef(modelo3)
##           mu        omega       alpha1       alpha2       alpha3 
## 0.0039021267 0.0001396828 0.2667119485 0.0000000100 0.0000000100

luego su modelo tendrá la forma

\[\begin{align*} h_t^2&=\alpha_0+\alpha_1u^2_{t-1}+\alpha_2u^2_{t-2}+\alpha_3u^2_{t-3}\\ &=0.0001396828+ 0.2667119485 u^2_{t-1}0.0000000100 u^2_{t-2}+0.0000000100 u_{t-3}^2 \end{align*}\]

cuyo gráfico tendrá la siguiente forma

predict(modelo3,n.ahead=5,plot=TRUE)

Modelo GARCH

En muchas series, como por ejemplo en las financieras, el número de retardos a utilizar es muy elevado, lo cual dificulta su estimación, ya que ello llevaría a un considerable número de iteraciones para alcanzar la solución del sistema planteado, o incluso llegar a no encontrarla debido a las condiciones que presentan sus parámetros. Engle ya propuso ciertas restricciones, como la de añadir pesos a los retardos con tal de reducir el número de parámetros a estimar, aunque estas no recogen cualquier caso.

Es por eso que Tim Peter Bollerslev, un economista danés que actualmente ejerce de profesor en la Universidad Duke en Estados Unidos, propone la generalizaci´on del modelo con heterocedasticidad condicional autoregresiva (GARCH), en el cual la varianza condicional no solo depende de los cuadrados de las perturbaciones, como en el modelo ARCH de Engle, sino que además también depende de las varianzas condicionales de períodos anteriores, es decir, de \(\sigma^2_t\) pasados.

GARCH(p,q): Definición y propiedades

Definición: Un proceso \(R_t, t \in \mathbb{Z}\) es un GARCH\((p,q)\) (Modelo con heterocedasticidad condicional autorregresiva generalizado) si cumple:

  • \(E(R_t|F_{t-1})=0\)

  • La varianza condicional \[\begin{align*} \text{Var}(R_t|F_{t-1})&=\sigma^2_t\\ &=w+\alpha_1R^2_{t-1}+...+\alpha_qR^2_{t-q}+\beta_1\sigma^2_{t-1}+...+\beta_p\sigma^2_{t-p}\\ &=w+\sum_{i=1}^{q}\alpha_iR^2_{t-i}+\sum_{j=1}^{p}\beta_jR^2_{t-j} \end{align*}\]

Donde las condiciones suficientes, aunque no necesarias, para que se cumpla que \(\sigma^2_t>0\), es decir, \(P(\sigma^2_t>0)=1\) son que \(w>0,\alpha_1\geq 0,...,\alpha_1\geq 0\), \(\beta_1\geq0,...,\beta_p\geq0\) y \(\displaystyle Z_t=\frac{R_t}{\sigma_t}\) con \(E(Z_t)=0\) y Var\((Z_t)=1\).

veamos un ejemplo, para el modelo GARCH(1,2):

\[\begin{align*} \sigma^2_t&=w+\alpha_1R^2_{t-1}+\alpha_2R^2_{t-2}+\beta_1\sigma^2_{t-1}\\ &=w+\alpha_1R^2_{t-1}+\alpha_2R^2_{t-2}+\beta_1\left(w+\alpha_1R^2_{t-1}+\alpha_2R^2_{t-2}+\beta_1\sigma^2_{t-2}\right)\\ &=w+\alpha_1R^2_{t-1}+\alpha_2R^2_{t-2}+\beta_1w+\alpha_1\beta_1R^2_{t-2}+\alpha_2\beta_1R^2_{t-3}+\beta^2_1\sigma^2_{t-2}\\ &=...\\ &=w(1+\beta_1+\beta_1^2+...)+\alpha_1\sum_{i=0}^{\infty}\beta_1^iR^2_{t-i-1}+\alpha_2\sum_{j=0}^{\infty}\beta_1^jR^2_{t-j-2}\\ &=\frac{w}{1-\beta_1}+\alpha_1R^2_{t-1}+(\alpha_1\beta_1+\alpha_2)\sum_{j=0}^{\infty}\beta_1^jR^2_{t-j-2} \end{align*}\]

Por lo tanto, para el modelo GARCH(1,2), \(0\leq \beta_1<1,w>0\),\(\alpha_1\geq 0\) y \(\alpha_1\beta_1+\alpha_2\geq0\) son condiciones necesarias y suficientes siempre y cuando \(\displaystyle \sum_{j=0}^{\infty}\beta_1^jR^2_{t-j-2}\) sea convergente.

Como \(\displaystyle Z_t=\frac{R_t}{\sigma_t}\) el modelo se representaría de la forma:

\[R_t=Z_t\sigma_t=Z_t\left(w+\alpha_1R^2_{t-1}+...+\alpha_qR^2_{t-q}+\beta_1\sigma^2_{t-1}+...+\beta_p\sigma^2_{t-p}\right)^{1/2}\]

A este modelo se le llama semi-strong GARCH.

Nota: Si además se cumple que \(\displaystyle Z_t=\frac{R_t}{\sigma_t}\) es i.i.d, se le llama strong GARCH.

Ahora, la varianza condicional \(\sigma^2_t\) es una función lineal del cuadrado de los \(q\) retardos y de las \(p\) varianzas condicionales anteriores.

Propiedades del modelo GARCH(p,q)

Teorema: Sea \(R_t\) un proceso semi-strong GARCH\((p,q)\) con Var\((R_t) =\sigma^2<\infty\), entonces:

\[\sigma^2=\frac{w}{1-\sum_{i=1}^{q}\alpha_i-\sum_{j=1}^{p}\beta_j}\]

con \(\displaystyle \sum_{i=1}^{q}\alpha_i+\sum_{j=1}^{p}\beta_j<1\).

Teorema: Sea \(R_t\) un proceso semi-strong GARCH(1,1) con Var\((R_t)=\sigma^2<\infty\) y \(Z_t\sim N(0,1).\) Entonces se cumple que el momento de cuarto orden de \(R_t\) es constante \((E(R^4_t)<\infty)\) siempre y cuando \(3\alpha_1+2\alpha_1\beta_1+\beta_1^2<1.\) En ese caso, su curtosis se define como:

\[\text{Kurt}(R_t)=\frac{E(R^4_t)}{(E(R^2_t))^2}=3+\frac{6\alpha_1^2}{1-\beta_1^2-2\alpha_1\beta_1-3\alpha_1^2}\]

Ejemplo:

Tomemos como ejemplo las mismas acciones que realizamos en el modelo ARCH, entonces

library(tidyquant)
library(TTR)
library(ggplot2)
library(knitr)
library(nortest)
library(MASS)
library(stats)
library(moments)
library(fGarch)
library(fBasics)


AAPL <- tq_get("AAPL", get="stock.prices", from= "2021-07-01", to= "2021-07-31") 
AAPL <- data.frame(AAPL) 
precio<-AAPL$close
n<-length(precio)
rendimiento<-c()
for (i in 1:n) {
  rendimiento[i]<-(precio[i+1]-precio[i])/(precio[i])
}
rendimiento<-na.omit(rendimiento)



modeloo<-garchFit(~garch(1,2),data=rendimiento)
## 
## Series Initialization:
##  ARMA Model:                arma
##  Formula Mean:              ~ arma(0, 0)
##  GARCH Model:               garch
##  Formula Variance:          ~ garch(1, 2)
##  ARMA Order:                0 0
##  Max ARMA Order:            0
##  GARCH Order:               1 2
##  Max GARCH Order:           2
##  Maximum Order:             2
##  Conditional Dist:          norm
##  h.start:                   3
##  llh.start:                 1
##  Length of Series:          20
##  Recursion Init:            mci
##  Series Scale:              0.01437149
## 
## Parameter Initialization:
##  Initial Parameters:          $params
##  Limits of Transformations:   $U, $V
##  Which Parameters are Fixed?  $includes
##  Parameter Matrix:
##                      U          V    params includes
##     mu     -2.18314636   2.183146 0.2183146     TRUE
##     omega   0.00000100 100.000000 0.1000000     TRUE
##     alpha1  0.00000001   1.000000 0.1000000     TRUE
##     gamma1 -0.99999999   1.000000 0.1000000    FALSE
##     beta1   0.00000001   1.000000 0.4000000     TRUE
##     beta2   0.00000001   1.000000 0.4000000     TRUE
##     delta   0.00000000   2.000000 2.0000000    FALSE
##     skew    0.10000000  10.000000 1.0000000    FALSE
##     shape   1.00000000  10.000000 4.0000000    FALSE
##  Index List of Parameters to be Optimized:
##     mu  omega alpha1  beta1  beta2 
##      1      2      3      5      6 
##  Persistence:                  0.9 
## 
## 
## --- START OF TRACE ---
## Selected Algorithm: nlminb 
## 
## R coded nlminb Solver: 
## 
##   0:     27.974431: 0.218315 0.100000 0.100000 0.400000 0.400000
##   1:     27.963261: 0.218306 0.0984128 0.0925684 0.394741 0.395780
##   2:     27.949471: 0.218248 0.107334 0.0894440 0.396070 0.399248
##   3:     27.928169: 0.218151 0.117980 0.0731233 0.390408 0.398297
##   4:     27.890153: 0.217961 0.149372 0.0485058 0.388109 0.405643
##   5:     27.869473: 0.218131 0.172259 0.0152574 0.386928 0.410137
##   6:     27.865757: 0.218255 0.189488 1.00000e-08 0.384710 0.414745
##   7:     27.865748: 0.217832 0.188354 1.00000e-08 0.384284 0.416287
##   8:     27.865702: 0.216424 0.180617 1.00000e-08 0.381283 0.427025
##   9:     27.865408: 0.211807 0.148607 1.00000e-08 0.368823 0.471576
##  10:     27.853406: 0.194260 0.0185893 1.00000e-08 0.319954 0.650566
##  11:     27.848538: 0.192854 1.00000e-06 1.00000e-08 0.311698 0.678974
##  12:     27.847976: 0.194193 1.00000e-06 1.00000e-08 0.309099 0.681123
##  13:     27.847154: 0.198642 1.00000e-06 1.00000e-08 0.302265 0.686626
##  14:     27.847079: 0.199179 1.00000e-06 1.00000e-08 0.302931 0.685954
##  15:     27.846367: 0.204428 1.00000e-06 1.00000e-08 0.315827 0.673622
##  16:     27.846114: 0.207954 1.00000e-06 1.00000e-08 0.330033 0.660255
##  17:     27.846096: 0.207903 1.00000e-06 1.00000e-08 0.333264 0.657223
##  18:     27.846089: 0.207671 1.00000e-06 1.00000e-08 0.336024 0.654557
##  19:     27.846076: 0.207397 1.00000e-06 1.00000e-08 0.341756 0.648926
##  20:     27.846033: 0.206793 1.00000e-06 1.00000e-08 0.363075 0.627856
##  21:     27.845912: 0.205537 1.00000e-06 1.00000e-08 0.429564 0.561954
##  22:     27.845619: 0.203048 1.00000e-06 1.00000e-08 0.596846 0.395823
##  23:     27.844859: 0.198144 1.00000e-06 1.00000e-08 0.949578 0.0449160
##  24:     27.844042: 0.197826 1.00000e-06 1.00000e-08 0.994378 1.00000e-08
##  25:     27.842566: 0.202157 1.00000e-06 1.00000e-08 0.993745 1.00000e-08
##  26:     27.842136: 0.205236 1.00000e-06 1.00000e-08 0.993071 1.00000e-08
##  27:     27.842136: 0.205174 1.00000e-06 1.00000e-08 0.993073 1.00000e-08
##  28:     27.842136: 0.205171 1.00000e-06 1.00000e-08 0.993072 1.00000e-08
## 
## Final Estimate of the Negative LLH:
##  LLH:  -57.00804    norm LLH:  -2.850402 
##           mu        omega       alpha1        beta1        beta2 
## 2.948618e-03 2.065397e-10 1.000000e-08 9.930725e-01 1.000000e-08 
## 
## R-optimhess Difference Approximated Hessian Matrix:
##                  mu        omega        alpha1        beta1        beta2
## mu      -102523.274     15878295     2606.8522     2951.804     2971.166
## omega  15878295.349 -26388600894 -4828633.9634 -5420412.575 -5442424.328
## alpha1     2606.852     -4828634     -910.9953    -1027.417    -1020.911
## beta1      2951.804     -5420413    -1027.4167    -1102.926    -1108.067
## beta2      2971.166     -5442424    -1020.9112    -1108.067    -1113.227
## attr(,"time")
## Time difference of 0.011554 secs
## 
## --- END OF TRACE ---
## 
## 
## Time to Estimate Parameters:
##  Time difference of 0.05629396 secs

ahora observemos los coeficientes del modelo

coef(modeloo)
##           mu        omega       alpha1        beta1        beta2 
## 2.948618e-03 2.065397e-10 1.000000e-08 9.930725e-01 1.000000e-08

luego el modelo tiene la siguiente forma

\[\begin{align*} \sigma_t^2&=\alpha_0+\alpha_1R^2_{t-1}+\alpha_2R^2_{t-2}+\beta_1\sigma^2_{t-1}\\ &=2.065397e-10+9.930725e-01 R^2_{t-1}+1.000000e-08 R^2_{t-2}+1.000000e-08\sigma^2_{t-1} \end{align*}\]

podemos predecir que ocurre en 5 tiempos depués, por ejemplo

cuyo gráfico tendrá la siguiente forma

predict(modeloo,n.ahead=5,plot=TRUE)