Analisis Series Temporales
Edison Duvan Garzon Gonzalez - Jhonatan Alejandro Gordillo Silva
10 de Julio de 2019
Para comenzar cargamos la serie de tiempo
Colcap <- read_excel("Colcap.xlsx", sheet = "Colcap",col_types=c("date","numeric"))
fechas=as.Date(Colcap$Fecha)
colcap=xts(Colcap$ValorCOLCAP,order.by = fechas)
plot(colcap)
A continuacion presentamos los retornos
retornos=diff(log(colcap))[2:length(colcap)]
plot(retornos)
Dado que esta serie fue diferenciada y transformada por logaritmo (definicion de retornos) podemos afirmar que esta serie es estacionaria, ahora ajustaremos un modelo ARMA como primera medida.
Primero veremos su funcin de autocorrelacion y autocorrelacion parcial.
acf(retornos, lag=sqrt(length(retornos)))
pacf(retornos,lag=sqrt(length(retornos)))
Dado que la funcion de autocorrelacion y autocorrelacion parcial decidimos que a lo mas es un ARMA(3,4) a continuacion ajustaremos un modelo por medio de la funcion auto.arima con el mejor .
modelo=auto.arima(retornos,max.p =3 ,max.q = 4, start.p = 0, start.q = 0, ic="aic")
coeftest(modelo)##
## z test of coefficients:
##
## Estimate Std. Error z value Pr(>|z|)
## ar1 1.259120 0.220683 5.7056 1.160e-08 ***
## ar2 -0.726117 0.233336 -3.1119 0.001859 **
## ar3 0.010325 0.039906 0.2587 0.795843
## ma1 -1.163300 0.219597 -5.2974 1.174e-07 ***
## ma2 0.640592 0.200908 3.1885 0.001430 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1El modelo que nos presenta es un ARMA(3,2) donde \(\phi_{3}\) no es significativo.
modelo1=arima(retornos, order=c(2,0,2), include.mean = FALSE)
coeftest(modelo1)##
## z test of coefficients:
##
## Estimate Std. Error z value Pr(>|z|)
## ar1 -0.59871 0.27161 -2.2043 0.02750 *
## ar2 0.21884 0.12662 1.7283 0.08394 .
## ma1 0.70018 0.27204 2.5738 0.01006 *
## ma2 -0.13713 0.12309 -1.1141 0.26525
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Ahora quitaremos la componente \(\theta_2\).
modelo2=arima(retornos, order=c(2,0,1), include.mean = FALSE)
coeftest(modelo2)##
## z test of coefficients:
##
## Estimate Std. Error z value Pr(>|z|)
## ar1 0.842050 0.184273 4.5696 4.887e-06 ***
## ar2 -0.087263 0.026094 -3.3442 0.0008253 ***
## ma1 -0.739621 0.183650 -4.0274 5.641e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Dado que el modelo ya se ajusto y todos sus parametros son significativos procedemos a un analisis de residuales.
residuales=modelo2$residuals
plot(residuales)
acf(residuales,ci.type="ma",lag.max = sqrt(length(residuales)))
pacf(residuales,lag.max = sqrt(length(residuales)))
Dado que no existe estructura de correlacion entramos a probar si existe evidencia de volatilidad:
arch.test(modelo2)## ARCH heteroscedasticity test for residuals
## alternative: heteroscedastic
##
## Portmanteau-Q test:
## order PQ p.value
## [1,] 4 1195 0
## [2,] 8 1486 0
## [3,] 12 1720 0
## [4,] 16 1874 0
## [5,] 20 1970 0
## [6,] 24 1999 0
## Lagrange-Multiplier test:
## order LM p.value
## [1,] 4 1100 0
## [2,] 8 469 0
## [3,] 12 289 0
## [4,] 16 207 0
## [5,] 20 163 0
## [6,] 24 133 0
Dado que todos los P-valores en todos los ordenes dieron exactamente 0 podemos rechazar la hipotesis nula de homocedasticidad, por tanto el modelo es heterocedastico y por tanto hay que modelar la volatilidad.
res <- residuales^2
plot(res,type='l')
acf(res,ci.type="ma",lag.max = sqrt(length(residuales)))
pacf(res,lag.max = sqrt(length(residuales)))
Se puede ver una alta correlacion hasta 12 en ambos asi que ajustaremos un auto.arima con p=12 y q=12.
mod=auto.arima(res,max.p =12 ,max.q =12, start.p = 0, start.q = 0, ic="aic")
coeftest(mod)##
## z test of coefficients:
##
## Estimate Std. Error z value Pr(>|z|)
## ar1 -0.767625 0.178928 -4.2901 1.786e-05 ***
## ar2 -0.143960 0.069476 -2.0721 0.03826 *
## ar3 -0.897861 0.056891 -15.7820 < 2.2e-16 ***
## ar4 -0.336255 0.105723 -3.1805 0.00147 **
## ar5 0.294241 0.042248 6.9646 3.293e-12 ***
## ma1 0.300789 0.180191 1.6693 0.09506 .
## ma2 -0.485088 0.090707 -5.3479 8.899e-08 ***
## ma3 0.519856 0.052074 9.9830 < 2.2e-16 ***
## ma4 -0.318924 0.076851 -4.1499 3.327e-05 ***
## ma5 -0.719362 0.120316 -5.9790 2.246e-09 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1ahora ajustaremos el modelo ARMA+GARCH.
modelosp5_1=garchFit(garch~arma(2,1)+garch(5,5),data=residuales,trace=F,include.mean = FALSE)## Warning in sqrt(diag(fit$cvar)): Se han producido NaNssummary(modelosp5_1)##
## Title:
## GARCH Modelling
##
## Call:
## garchFit(formula = garch ~ arma(2, 1) + garch(5, 5), data = residuales,
## include.mean = FALSE, trace = F)
##
## Mean and Variance Equation:
## data ~ arma(2, 1) + garch(5, 5)
## <environment: 0x000000001df73b80>
## [data = residuales]
##
## Conditional Distribution:
## norm
##
## Coefficient(s):
## ar1 ar2 ma1 omega alpha1
## 5.9989e-01 1.4602e-02 -6.0624e-01 5.3867e-06 1.5864e-01
## alpha2 alpha3 alpha4 alpha5 beta1
## 4.9759e-02 1.0000e-08 1.0000e-08 1.0000e-08 6.0409e-01
## beta2 beta3 beta4 beta5
## 1.0000e-08 1.0000e-08 1.0000e-08 1.3618e-01
##
## Std. Errors:
## based on Hessian
##
## Error Analysis:
## Estimate Std. Error t value Pr(>|t|)
## ar1 5.999e-01 1.175e-01 5.106 3.29e-07 ***
## ar2 1.460e-02 2.152e-02 0.679 0.497
## ma1 -6.062e-01 1.170e-01 -5.184 2.17e-07 ***
## omega 5.387e-06 NA NA NA
## alpha1 1.586e-01 2.195e-02 7.228 4.92e-13 ***
## alpha2 4.976e-02 1.122e-01 0.444 0.657
## alpha3 1.000e-08 NA NA NA
## alpha4 1.000e-08 1.196e-01 0.000 1.000
## alpha5 1.000e-08 7.069e-02 0.000 1.000
## beta1 6.041e-01 8.555e-01 0.706 0.480
## beta2 1.000e-08 7.781e-01 0.000 1.000
## beta3 1.000e-08 NA NA NA
## beta4 1.000e-08 6.852e-01 0.000 1.000
## beta5 1.362e-01 3.450e-01 0.395 0.693
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Log Likelihood:
## 8954.959 normalized: 3.300759
##
## Description:
## Fri Apr 24 12:26:01 2020 by user: edisonduvan
##
##
## Standardised Residuals Tests:
## Statistic p-Value
## Jarque-Bera Test R Chi^2 751.2317 0
## Shapiro-Wilk Test R W 0.9769415 0
## Ljung-Box Test R Q(10) 13.3206 0.2062948
## Ljung-Box Test R Q(15) 14.96094 0.4542345
## Ljung-Box Test R Q(20) 22.43903 0.3171703
## Ljung-Box Test R^2 Q(10) 70.80817 3.094758e-11
## Ljung-Box Test R^2 Q(15) 72.02754 1.939908e-09
## Ljung-Box Test R^2 Q(20) 73.01429 5.821152e-08
## LM Arch Test R TR^2 6.99188 0.8581497
##
## Information Criterion Statistics:
## AIC BIC SIC HQIC
## -6.591197 -6.560721 -6.591250 -6.580179ahora sus predicciones
prediccion=predict(n.ahead=5,modelosp5_1)## Warning in a_vec[i] <- ar[1:min(u2, i - 1)] * a_vec[(i - 1):(i - u2)] + :
## número de items para para sustituir no es un múltiplo de la longitud del
## reemplazo
## Warning in a_vec[i] <- ar[1:min(u2, i - 1)] * a_vec[(i - 1):(i - u2)] + :
## número de items para para sustituir no es un múltiplo de la longitud del
## reemplazopredict(n.ahead=5,modelosp5_1,plot=TRUE)## Warning in a_vec[i] <- ar[1:min(u2, i - 1)] * a_vec[(i - 1):(i - u2)] + :
## número de items para para sustituir no es un múltiplo de la longitud del
## reemplazo
## Warning in a_vec[i] <- ar[1:min(u2, i - 1)] * a_vec[(i - 1):(i - u2)] + :
## número de items para para sustituir no es un múltiplo de la longitud del
## reemplazo
## meanForecast meanError standardDeviation lowerInterval upperInterval
## 1 -4.486047e-05 0.006690745 0.006690745 -0.01315848 0.01306876
## 2 5.193726e-05 0.007043634 0.007043506 -0.01375333 0.01385721
## 3 3.050177e-05 0.007332725 0.007332233 -0.01434137 0.01440238
## 4 1.905623e-05 0.007464064 0.007463406 -0.01461024 0.01464835
## 5 1.187711e-05 0.007531123 0.007530375 -0.01474885 0.01477261