Introducción Modelo Sarimax
Wilson Sandoval
02/04/2022
Series de Tiempo Multivariadas
La modelización univariante de una serie temporal (útiles para conocer las características de una serie temporal de tendencia, ciclo, estacionalidad, o para predecir), pero no tienen en cuenta la interrelación entre distintas variables (variable objeto de estudio con otras variables relevantes), y por lo tanto su utilidad puede verse limitada
- Por ejemplo las ventas de una empresa están en relación con la renta de los consumidores, o la demanda de turistas está en relación con el crecimiento del PIB o de los precios relativos, etc.
Para resolver estas limitaciones se presentan los modelos multivariantes, que son modelos que sí tienen en cuenta la interrelación entre variables de distintas series temporales, luego estos modelos contienen varias series temporales.
Concretamente estudian la relación entre una variable de interés y una serie de variables explicativas (que influyen en la variable de interés). En el marco multivariante se considera el pasado tanto de la variable que se quiere explicar, como el de las variables que están relacionadas con dicha variable.
Los modelos multivariantes de series temporales son una generalización de los modelos univariantes, la diferencia está en que en vez de una sola variable hay n variables (en vez de una serie, hay varias series).
Series de Tiempo Multivariadas
El análisis de series temporales multivariadas considera simultáneamente series temporales múltiples.
En general, mucho más complicado que el análisis de series de tiempo univariante, especialmente cuando el número de series consideradas es grande.
Toma de decisiones. Los objetivos del análisis de series temporales multivariadas incluyen, por lo tanto:
- Mejorar la precisión de la predicción.
- Estudiar las relaciones dinámicas entre variables.
Sea \(\mathbf{z_t} = (z_{1t},\cdots, z_{kt})\) una serie de tiempo \(k\)-dimensional observada igualmente puntos de tiempo espaciados.
Por ejemplo, dejemos que
\(z_{1t}\) sea el real bruto interno trimestral de EE. UU. (GDP)
\(z_{2t}\) la tasa de desempleo civil trimestral de EE. UU.
Mediante el estudio de \(z_{1t}\) y \(z_{2t}\) conjuntamente, Se puede evaluar la dependencia temporal y contemporánea entre el PIB y la tasa de desempleo.
Cargar las librerias
library(tseries)
library(vars)
library(highcharter)
library(dplyr)
library(ggfortify)
library(plotly)data("USeconomic")
class(USeconomic)## [1] "mts" "ts"USeconomic%>%DT::datatable()plot(USeconomic)
hchart(USeconomic)%>%hc_add_theme(hc_theme_538())s<-autoplot(USeconomic, main = "Series", xlab = "Tiempo", ylab="Useconomic", ts.colour="blue")
ggplotly(s)## Warning: `group_by_()` is deprecated as of dplyr 0.7.0.
## Please use `group_by()` instead.
## See vignette('programming') for more help
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_warnings()` to see where this warning was generated.data(Canada)
plot(Canada)
hchart(Canada) %>% hc_add_theme(hc_theme_economist()) %>% hc_title(text="Canada")autoplot(Canada)
b<-autoplot(Canada)+
xlab("Tiempo")+
ggtitle("Series de tiempo Canada")
ggplotly(b)Estadísticamente hablando, una serie temporal \(k\)-dimensional \(\textbf{z}_t = (z_{1t},..., z_{kt})'\) es un vector aleatorio que consta de \(k\) variables aleatorias.
Modelo ARMAX
Los modelos SARIMAX (auto regresivos integrados de medias móviles estacional con variables explicativas) son una estructura matemática que contiene todas las bondades de un modelo SARIMA adicionalmente, pueden capturar información sobre variables exógenas que permiten a entender y pronosticar la variable de interés. Para una única variable de interés \(y\), y una única variable independiente \(x\).
El modelo que describe su relación es\[\begin{align}y_t = \beta_1 & x_t+\cdots+ \beta_{k} x_{t_{k-1}}+ \phi_1 y_{t-1}+\cdots+\phi_p y_{t-p}+\\ &\theta_1 z_{t-1}+\cdots \theta_{q} z_{t-q} + z_t \end{align}\]
- \(x\) es una variable exógena
- \(z\) es un ruido blanco
Si bien esto parece sencillo, una desventaja es que el coeficiente covariable es difícil de interpretar.
El valor de \(\beta\) no es el efecto en \(y_t\) cuando la \(x_t\) se incrementa en uno (ya que está en regresión).
La presencia de valores rezagados de la variable de respuesta en el lado derecho de la ecuación significa que \(\beta\) solo puede interpretarse condicionalmente en el valor de los valores anteriores de la variable de respuesta, lo cual es poco intuitivo.
SOI Y Recruitment
Cargamos las dos series de tiempo
rec = scan("SARIMAX/recruit.dat")
soi=scan("SARIMAX/soi.dat")También podemos estar interesados en analizar varias series de tiempo a la vez. La Figura muestra los valores mensuales de una serie ambiental llamada Southern Índice de oscilación (SOI) y recluitment asociado (número de peces nuevos)
rec1<-ts(rec, start = c(1950,1), end = c(1986,9), frequency = 12)
soi1<-ts(soi, start = c(1950,1), end = c(1986,9), frequency = 12)
soi.rec<-cbind(rec1, soi1)
autoplot(soi.rec)
hchart(soi.rec)Ambas series tienen una duración de 453 meses.que van a lo largo de los años 1950-1987. El SOI mide los cambios en la presión del aire, relacionado con las temperaturas de la superficie del mar en el Océano Pacífico central. El central El Pacífico se calienta cada tres a siete años debido al efecto de El Niño, que tiene sido culpado, en particular, por las inundaciones de 1997 en las porciones del medio oeste de los Estados Unidos.
Ambas series tienden a exhibir repetitivos comportamiento, con ciclos que se repiten regularmente y que son fácilmente visibles. Este comportamiento periódico es de interés porque los procesos subyacentes de interés pueden ser regular y la tasa o frecuencia de oscilación que caracteriza el comportamiento de la serie subyacente ayudaría a identificarlos.
También se puede comentar que los ciclos del SOI se repiten a un ritmo más rápido que los de la serie Serie de recruitmen. La serie Recruitment también muestra varios tipos de oscilaciones, una frecuencia más rápida que parece repetirse cada 12 meses y una frecuencia más lenta que parece repetirse aproximadamente cada 50 meses.
Los dos las series también tienden a estar algo relacionadas; es fácil imaginar que de alguna manera La población de peces depende del SOI. Quizás incluso una relación rezagada existe, con los cambios de señalización SOI en la población de peces. Esta posibilidad
autoplot(stl(soi1, s.window = "periodic"))
autoplot(stl(rec1, s.window = "periodic"))
Serie Univariada
z<-ts(rec, start = c(1950,1), end = c(1986,9), frequency = 12)
plot(z)
hchart(z)Descomposición de la serie Recruitment
autoplot(stl(z, s.window = "periodic"))
par(mfrow=c(3,3))
plot(z)
acf(z,lag.max = 40)
pacf(z, lag.max = 40)
plot(diff(z))
acf(diff(z), lag.max = 40)
pacf(diff(z), lag.max = 40)
plot(diff(diff(z), 12))
acf(diff(diff(z), 12), lag.max = 40)
pacf(diff(diff(z), 12),lag.max = 40)
par(mfrow=c(1,3))
plot(diff(diff(z), 12))
acf(diff(diff(z), 12), lag.max = 40)
pacf(diff(diff(z), 12),lag.max = 40)
- Dickey- fULLER y Phillis Perron
adf.test(diff(z)) ## Warning in adf.test(diff(z)): p-value smaller than printed p-value##
## Augmented Dickey-Fuller Test
##
## data: diff(z)
## Dickey-Fuller = -9.6547, Lag order = 7, p-value = 0.01
## alternative hypothesis: stationarypp.test(diff(z)) ## Warning in pp.test(diff(z)): p-value smaller than printed p-value##
## Phillips-Perron Unit Root Test
##
## data: diff(z)
## Dickey-Fuller Z(alpha) = -239.6, Truncation lag parameter = 5, p-value
## = 0.01
## alternative hypothesis: stationaryadf.test(diff(diff(z), 12)) ## Warning in adf.test(diff(diff(z), 12)): p-value smaller than printed p-value##
## Augmented Dickey-Fuller Test
##
## data: diff(diff(z), 12)
## Dickey-Fuller = -6.5087, Lag order = 7, p-value = 0.01
## alternative hypothesis: stationarypp.test(diff(diff(z), 12)) ## Warning in pp.test(diff(diff(z), 12)): p-value smaller than printed p-value##
## Phillips-Perron Unit Root Test
##
## data: diff(diff(z), 12)
## Dickey-Fuller Z(alpha) = -275.66, Truncation lag parameter = 5, p-value
## = 0.01
## alternative hypothesis: stationarysd(diff(z))## [1] 10.87691sd(diff(diff(z), 12))## [1] 14.0919library(forecast)## Registered S3 methods overwritten by 'forecast':
## method from
## autoplot.Arima ggfortify
## autoplot.acf ggfortify
## autoplot.ar ggfortify
## autoplot.bats ggfortify
## autoplot.decomposed.ts ggfortify
## autoplot.ets ggfortify
## autoplot.forecast ggfortify
## autoplot.stl ggfortify
## autoplot.ts ggfortify
## fitted.ar ggfortify
## fortify.ts ggfortify
## residuals.ar ggfortifyauto.arima(z)## Series: z
## ARIMA(1,1,0)(0,0,2)[12]
##
## Coefficients:
## ar1 sma1 sma2
## 0.3698 0.0733 0.1644
## s.e. 0.0447 0.0495 0.0451
##
## sigma^2 estimated as 95.78: log likelihood=-1626.92
## AIC=3261.83 AICc=3261.92 BIC=3278.18- SARIMA(1,1,0)(0,1,1)[12] BIC=3278.18
- SARIMA(0,1,1)(0,1,1)[12] BIC= 3156.718
- SARIMA(1,1,0)(0,1,1)[12] BIC= 3155.229
- SARIMA(4,1,0)(0,1,1)[12] BIC= 3157.041
- SARIMA(0,1,11)(0,1,1)[12] BIC= 3172.221 LAGS 1,4,5,8,9,10,11
modelo1<-stats::arima(z,
order=c(1,1,0),
seasonal=list(order=c(0,1,1),
period=12),fixed = c(NA,NA))
modelo1##
## Call:
## stats::arima(x = z, order = c(1, 1, 0), seasonal = list(order = c(0, 1, 1),
## period = 12), fixed = c(NA, NA))
##
## Coefficients:
## ar1 sma1
## 0.2869 -1.000
## s.e. 0.0463 0.059
##
## sigma^2 estimated as 80.68: log likelihood = -1568.53, aic = 3143.05tt <- modelo1$coef[which(modelo1$coef!=0)]/sqrt(diag(modelo1$var.coef))
1 - pt(abs(tt),(modelo1$nobs - length(modelo1$coef[which(modelo1$coef!=0)])))## ar1 sma1
## 6.951388e-10 0.000000e+00BIC del modelo
BIC(modelo1)## [1] 3155.229Residuales
et<-residuals(modelo1)
x1.fit <- fitted(modelo1)
par(mfrow=c(3,2))
plot(z,type="l",lty=2)
lines(x1.fit,type="l",col="red")
plot(scale(et),type="l",main="Residuales")
abline(h=2*sqrt(var(scale(et))),col="red",lty=2)
abline(h=-2*sqrt(var(scale(et))),col="red",lty=2)
acf(et)
pacf(et)
#qqPlot(scale(et))
acf(abs(et)) #Mide Estructura Heterocedástica
Test de Autocorrelacion de Ljung-Box
\(H_o: r1=r2=r3=...=rlag=0\)
$H_a: $
Box.test(et,lag=12,type="Ljung-Box")##
## Box-Ljung test
##
## data: et
## X-squared = 15.768, df = 12, p-value = 0.2021Test de Normalidad basado en Sesgo y Curtosis
\(H_o:\) Datos provienen de una Dist. Normal \(H_a\): Los datos no provienen de una Dist. Normal
jarque.bera.test(et)##
## Jarque Bera Test
##
## data: et
## X-squared = 107.44, df = 2, p-value < 2.2e-16Test de Aleatoriedad
\(H_o:\) Residuales exhiben un comport. de Aleatoriedad \(H_a:\) Residuales no exhiben estructura (Tendencia, o cualquier otro comportamiento predecible)
runs.test(as.factor(sign(et)), alternative="two.sided")##
## Runs Test
##
## data: as.factor(sign(et))
## Standard Normal = -1.5563, p-value = 0.1196
## alternative hypothesis: two.sidedautoplot(modelo1)
Pronóstico Fuera Muestra
plot(forecast(modelo1,h=12, fan=F))
lines(fitted(modelo1), col="red")
lines(ts(rec, start = c(1950,1),frequency = 12),col="green")
Funciones de correlación cruzada y relaciones entre 2 series temporales
El problema básico que estamos considerando es la descripción y el modelado de la relación entre dos series de tiempo.\
En la relación entre dos series de tiempo \(y_t\) y \(x_t\) La serie \(y_t\) puede estar relacionado con retrasos pasados de la serie \(x\) .
La función de correlación cruzada de muestra (CCF) es útil para identificar retrasos de la variable \(x\) que podrían ser predictores útiles de \(y_t\)
En R, el CCF de muestra se define como el conjunto de correlaciones de muestra entre \(x_{t+h}\) y \(y_t\) para \(h = 0, ± 1, ± 2, ± 3\), y así sucesivamente. Un valor negativo para \(h\) es una correlación entre la \(x-\)variable en un momento anterior y la variable y en el tiempo \(t\). Por ejemplo, considere \(h= −2\). El valor CCF daría la correlación entre \(x_{t-2}\) y \(y_t\) .
Cuando uno o más \(x_{t+h}\) , con \(h\) negativo , son predictores de \(y_t\), a veces se dice que \(x\) lleva a \(y\)} .
Cuando uno o más \(x_{t+h}\) , con \(h\) positivo , son predictores de \(y_t\), a veces se dice que \(x\) se retrasa \(y\)
En algunos problemas, el objetivo puede ser identificar qué variable está liderando y cuál está retrasada.
- Sin embargo, en muchos problemas que consideramos, examinaremos las variables \(x\) como una variable principal de la variable \(y\) porque queremos utilizar los valores de la variable \(x\) para predecir los valores futuros de \(y\).
Por lo tanto, generalmente veremos qué sucede con los valores negativos de en la parcela CCF.
El comando CCF es
\(\texttt{ccf(x-variable name, y-variable name)}\).
Si desea especificar cuántos retrasos mostrar, agregue ese número como argumento del comando. Por ejemplo,
\(\texttt {ccf(x, y, 50) dará el CCF para valores de = 0, ± 1, ..., ± 50.}\)
Ejemplo
Índice de oscilación meridional y poblaciones de peces en el hemisferio sur.\
Se describe la relación entre una medida del clima llamada y “recruitment”, una medida de la población de peces en el hemisferio sur. Los datos son estimaciones mensuales para \(n = 453\) meses.
hchart(soi.rec)funcion de correlación cruzada
ccfvalues<-ccf(soi,rec,20)
ccfvalues##
## Autocorrelations of series 'X', by lag
##
## -20 -19 -18 -17 -16 -15 -14 -13 -12 -11 -10
## -0.108 -0.198 -0.253 -0.222 -0.149 -0.092 -0.076 -0.103 -0.175 -0.267 -0.369
## -9 -8 -7 -6 -5 -4 -3 -2 -1 0 1
## -0.476 -0.560 -0.598 -0.599 -0.527 -0.297 -0.146 -0.042 0.011 0.025 -0.013
## 2 3 4 5 6 7 8 9 10 11 12
## -0.086 -0.154 -0.228 -0.259 -0.232 -0.144 -0.017 0.094 0.154 0.174 0.162
## 13 14 15 16 17 18 19 20
## 0.118 0.043 -0.057 -0.129 -0.156 -0.131 -0.049 0.060