# Parámetros conocidos que luego vamos a estimar
mu <- 2
w <- 0.1
alpha1 <- 0.5
alpha2 <- 0.3
beta1 <- 0.05
# Simulación de n rendimientos
n <- 1000
set.seed(1234) # fijamos la semilla
et <- rnorm(n=2, mean=0, sd=1)
varianza <- numeric(n)
# La recursividad de presente de t=3 en adelante
for (i in 3:n) {
varianza[i] <- w + alpha1*(et[i-1])^2 + alpha2*(et[i-2])^2 + beta1*varianza[i-1]
et[i] <- rnorm(n=1, mean=0, sd=sqrt(varianza[i]))
}
rt <- mu + etAplicación de simulación Monte Carlo: estimación de un modelo GARCH
Objetivo
Evaluar el efecto del tamaño de muestra en la estimación de los parámetros de un modelo garch.
Introducción
Vamos a simular rendimientos diarios \(r_t\) y la varianza \(\sigma_t^2\) a usar estará condicionada a la información pasada según un modelo GARCH. Vamos a suponer que la media \(\mu\) de los rendimientos, o media de largo plazo, es constante. Cada día hay un choque aleatorio que se suma a la media, y vamos a suponer que esos choques aleatorios \({e_t}\) son i.i.d con distribución \(N(0, \sigma^2_t)\).
La ecuación del rendimiento en el tiempo \(t\) será:
\[ r_t = \mu + e_t \]
La varianza diaria del proceso \(\sigma_t^2\) tiene la siguiente estructura:
\[\sigma_t^2 = \omega + \alpha_1 e^2_{t-1} + \alpha_2 e^2_{t-2} + \beta_1 \sigma_{t-1}^2,\]
donde \(\omega \geq 0\), \(\alpha_1 >0\), \(\alpha_2 >0\), \(\beta_1 > 0\), adicionalmente se debe cumplir que \({\omega + \alpha_1 + \alpha_2 + \beta_1 < 1}\) (No se para que sirve esta condición pero Henry la puso).
Los parámetros del modelo los vamos a suponer conocidos con los siguientes valores:
\[ \boldsymbol{\Theta} = (\mu=2, \omega=0.1, \alpha_1=0.5, \alpha_2=0.3, \beta_1=0.05)^\top \]
El objetivo es simular datos del modelo anterior y luego estimar \(\boldsymbol{\Theta}\).
Simulación de \(n=1000\) retornos
El código para definir \(\boldsymbol{\Theta}\) y simular los datos es el siguiente:
Vamos a mostrar los primeros diez rendimientos simulados.
rt[1:10] [1] 0.7929343 2.2774292 2.8227370 0.3574797 2.5556330 2.5421065 1.6378975
[8] 1.7140902 1.7514576 1.6383507Ahora vamos a dibujar los rendimientos \(r_t\) y las varianzas \(\sigma_t^2\) para las primeras 50 observaciones, no para las mil observaciones.
par(mfrow=c(1, 2))
plot(x=1:50, y=rt[1:50], xlab="Tiempo", type = "b", las=1)
plot(x=1:50, y=varianza[1:50], xlab="Tiempo", type = "b", las=1)
Recordemos que el rendimiento sigue esta estructura \(r_t = \mu + e_t\), es decir, que \(r_t\) oscila alrededor de \(\mu\) pero la varianza de \(e_t\) la que tiene dependencia autoregresiva.
Estimación de los parámetros
Para estimar el vector de parámetros \(\boldsymbol{\Theta}\) vamos a usar el paquete rugarch, para mayor información se puede consultar la página del paquete.
A continuación se muestra el código para definir y estimar el modelo.
library("rugarch")
# Paso 1: definir el modelo
garchspec <- ugarchspec(mean.model = list(armaOrder = c(0,0)),
variance.model = list(garchOrder = c(2,1)),
distribution.model = "norm")
# Paso 2: estimar el modelo
garchfit21 <- ugarchfit(data = rt, spec = garchspec)Para ver los resultados del modelo hacemos lo siguiente.
garchfit21
*---------------------------------*
* GARCH Model Fit *
*---------------------------------*
Conditional Variance Dynamics
-----------------------------------
GARCH Model : sGARCH(2,1)
Mean Model : ARFIMA(0,0,0)
Distribution : norm
Optimal Parameters
------------------------------------
Estimate Std. Error t value Pr(>|t|)
mu 1.989984 0.013798 144.2224 0.000000
omega 0.086276 0.015442 5.5870 0.000000
alpha1 0.453253 0.066676 6.7978 0.000000
alpha2 0.168325 0.088232 1.9077 0.056425
beta1 0.204166 0.098731 2.0679 0.038650
Robust Standard Errors:
Estimate Std. Error t value Pr(>|t|)
mu 1.989984 0.014035 141.7922 0.000000
omega 0.086276 0.015342 5.6236 0.000000
alpha1 0.453253 0.073966 6.1278 0.000000
alpha2 0.168325 0.081455 2.0665 0.038783
beta1 0.204166 0.094473 2.1611 0.030687
LogLikelihood : -810.971
Information Criteria
------------------------------------
Akaike 1.6319
Bayes 1.6565
Shibata 1.6319
Hannan-Quinn 1.6413
Weighted Ljung-Box Test on Standardized Residuals
------------------------------------
statistic p-value
Lag[1] 0.709 0.3998
Lag[2*(p+q)+(p+q)-1][2] 1.169 0.4466
Lag[4*(p+q)+(p+q)-1][5] 1.757 0.6768
d.o.f=0
H0 : No serial correlation
Weighted Ljung-Box Test on Standardized Squared Residuals
------------------------------------
statistic p-value
Lag[1] 0.4667 0.4945
Lag[2*(p+q)+(p+q)-1][8] 1.5244 0.9244
Lag[4*(p+q)+(p+q)-1][14] 4.3491 0.8450
d.o.f=3
Weighted ARCH LM Tests
------------------------------------
Statistic Shape Scale P-Value
ARCH Lag[4] 0.1170 0.500 2.000 0.7323
ARCH Lag[6] 0.2992 1.461 1.711 0.9454
ARCH Lag[8] 0.3745 2.368 1.583 0.9903
Nyblom stability test
------------------------------------
Joint Statistic: 0.6688
Individual Statistics:
mu 0.12885
omega 0.04460
alpha1 0.06781
alpha2 0.15001
beta1 0.18570
Asymptotic Critical Values (10% 5% 1%)
Joint Statistic: 1.28 1.47 1.88
Individual Statistic: 0.35 0.47 0.75
Sign Bias Test
------------------------------------
t-value prob sig
Sign Bias 1.932 0.05369 *
Negative Sign Bias 1.328 0.18435
Positive Sign Bias 1.260 0.20814
Joint Effect 4.450 0.21677
Adjusted Pearson Goodness-of-Fit Test:
------------------------------------
group statistic p-value(g-1)
1 20 13.48 0.8131
2 30 17.30 0.9573
3 40 36.96 0.5633
4 50 34.90 0.9358
Elapsed time : 0.3453538 Con \(n = 1000\) observaciones se obtienen parámetros estimados más o menos cercanos a los verdaderos valores, por lo menos para los tres primeros parámetros. El vector estimado de parámetros es \(\hat{\boldsymbol{\Theta}} = (1.989, 0.086, 0.453, 0.168, 0.204)^\top\) mientras que el verdadero vector de parámetros es \(\boldsymbol{\Theta} = (2, 0.1, 0.5, 0.3, 0.05)^\top\).
¿Cuál será el ECM de esta estimación?
theta <- c(2, 0.1, 0.5, 0.3, 0.05)
theta_hat <- c(1.989, 0.086, 0.453, 0.168, 0.204)
ecm <- mean((theta - theta_hat)^2)
ecm[1] 0.0087332De la salida anterior se observa que todos los parámetros fueron significativos a un nivel del 5%.
Simulación de \(n=10000\) retornos
¿Será que la estimación de los parámetros mejora si aumentamos el número de observaciones?
Vamos a explorar.
n <- 10000
set.seed(1234) # fijamos la semilla
et <- rnorm(n=2, mean=0, sd=1)
varianza <- numeric(n)
# La recursividad de presente de t=3 en adelante
for (i in 3:n) {
varianza[i] <- w + alpha1*(et[i-1])^2 + alpha2*(et[i-2])^2 + beta1*varianza[i-1]
et[i] <- rnorm(n=1, mean=0, sd=sqrt(varianza[i]))
}
rt <- mu + etEstimemos el vector de parámetros.
# Paso 1: definir el modelo
garchspec <- ugarchspec(mean.model = list(armaOrder = c(0,0)),
variance.model = list(garchOrder = c(2,1)),
distribution.model = "norm")
# Paso 2: estimar el modelo
garchfit21 <- ugarchfit(data = rt, spec = garchspec)
garchfit21
*---------------------------------*
* GARCH Model Fit *
*---------------------------------*
Conditional Variance Dynamics
-----------------------------------
GARCH Model : sGARCH(2,1)
Mean Model : ARFIMA(0,0,0)
Distribution : norm
Optimal Parameters
------------------------------------
Estimate Std. Error t value Pr(>|t|)
mu 1.998377 0.004230 472.4764 0.000000
omega 0.098225 0.004392 22.3642 0.000000
alpha1 0.474999 0.020001 23.7488 0.000000
alpha2 0.260740 0.022692 11.4906 0.000000
beta1 0.073792 0.022586 3.2671 0.001086
Robust Standard Errors:
Estimate Std. Error t value Pr(>|t|)
mu 1.998377 0.004047 493.7932 0.000000
omega 0.098225 0.004150 23.6688 0.000000
alpha1 0.474999 0.020349 23.3421 0.000000
alpha2 0.260740 0.021254 12.2679 0.000000
beta1 0.073792 0.021801 3.3848 0.000712
LogLikelihood : -8004.264
Information Criteria
------------------------------------
Akaike 1.6019
Bayes 1.6055
Shibata 1.6019
Hannan-Quinn 1.6031
Weighted Ljung-Box Test on Standardized Residuals
------------------------------------
statistic p-value
Lag[1] 0.04614 0.8299
Lag[2*(p+q)+(p+q)-1][2] 0.63586 0.6331
Lag[4*(p+q)+(p+q)-1][5] 1.83159 0.6587
d.o.f=0
H0 : No serial correlation
Weighted Ljung-Box Test on Standardized Squared Residuals
------------------------------------
statistic p-value
Lag[1] 0.04702 0.8283
Lag[2*(p+q)+(p+q)-1][8] 3.65930 0.5656
Lag[4*(p+q)+(p+q)-1][14] 5.49053 0.7067
d.o.f=3
Weighted ARCH LM Tests
------------------------------------
Statistic Shape Scale P-Value
ARCH Lag[4] 2.025 0.500 2.000 0.1547
ARCH Lag[6] 2.780 1.461 1.711 0.3420
ARCH Lag[8] 3.258 2.368 1.583 0.4963
Nyblom stability test
------------------------------------
Joint Statistic: 0.6601
Individual Statistics:
mu 0.08059
omega 0.21301
alpha1 0.10729
alpha2 0.04307
beta1 0.06519
Asymptotic Critical Values (10% 5% 1%)
Joint Statistic: 1.28 1.47 1.88
Individual Statistic: 0.35 0.47 0.75
Sign Bias Test
------------------------------------
t-value prob sig
Sign Bias 0.6339 0.5261
Negative Sign Bias 0.2274 0.8201
Positive Sign Bias 1.1763 0.2395
Joint Effect 1.4358 0.6972
Adjusted Pearson Goodness-of-Fit Test:
------------------------------------
group statistic p-value(g-1)
1 20 25.42 0.1471
2 30 32.04 0.3180
3 40 36.26 0.5953
4 50 49.58 0.4500
Elapsed time : 1.712335 Con \(n = 10000\) observaciones se obtienen parámetros estimados más cercanas a los verdaderos valores. El vector estimado de parámetros es \(\hat{\boldsymbol{\Theta}} = (1.998, 0.098, 0.474, 0.260, 0.073)^\top\) mientras que el verdadero vector de parámetros es \(\boldsymbol{\Theta} = (2, 0.1, 0.5, 0.3, 0.05)^\top\).
¿Cuál será el ECM de esta estimación?
theta <- c(2, 0.1, 0.5, 0.3, 0.05)
theta_hat <- c(1.998, 0.098, 0.474, 0.260, 0.073)
ecm <- mean((theta - theta_hat)^2)
ecm[1] 0.0005626¿Cómo afecta/beneficia la disminución/aumento de \(n\) la estimación y el ECM?
De los dos ejemplos anteriores vimos que \(ECM_{n=1000}=0.009\) mientras que \(ECM_{n=10000}=0.00056\), claramente el ECM disminuyó al aumentar el tamaño de la muestra \(n\).
¿Cómo afecta \(n\) la estimación de los parámetros?
Tarea
- Realizar un estudio de simulación para ver lo que sucede con el ECM cuando se tienen valores de \(n=50, 100, 200, 400, 800, 1600, 3200, 6400\).
- Para cada valor de \(n\) se deben hacer 50 repeticiones, es decir, tener 50 valores de ECM para luego promediarlos.
- Construir una figura del promedio del ECM versus \(n\).