ARMAX

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 \(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)

## Warning: package 'ggfortify' was built under R version 4.3.1

library(plotly)

data("USeconomic")
class(USeconomic)

## [1] "mts" "ts"

USeconomic%>%DT::datatable()

?USeconomic

## starting httpd help server ... done

plot(USeconomic)

hchart(USeconomic)%>%hc_add_theme(hc_theme_538())

## Warning: Deprecated function. Use the `create_axis` function.

s<-autoplot(USeconomic, main = "Series", xlab = "Tiempo", ylab="Useconomic", ts.colour="blue")
ggplotly(s)

?Canada

data(Canada)
plot(Canada)

hchart(Canada) %>% hc_add_theme(hc_theme_economist()) %>% hc_title(text="Canada")

## Warning: Deprecated function. Use the `create_axis` function.

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("recruit.dat")
soi=scan("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)

## Warning: Deprecated function. Use the `create_axis` function.

• 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 Pacífico se calienta cada tres a siete años debido al efecto de el Niño, que ha 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 recruitment. 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.

-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)

• 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: stationary

pp.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: stationary

adf.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: stationary

pp.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: stationary

sd(diff(z))

## [1] 10.87691

sd(diff(diff(z), 12))

## [1] 14.0919

library(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           ggfortify

auto.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 = 95.78:  log likelihood = -1626.92
## AIC=3261.83   AICc=3261.92   BIC=3278.18

library(TSstudio)

ts_cor(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 = 50)

• 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.05

tt <- 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.951273e-10 0.000000e+00

BIC del modelo

BIC(modelo1)

## [1] 3155.229

Residuales

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:\) Al menos una es dif de cero

Box.test(et,lag=12,type="Ljung-Box")

## 
##  Box-Ljung test
## 
## data:  et
## X-squared = 15.768, df = 12, p-value = 0.2021

tsdiag(modelo1, gof.lag=20)

Test de Normalidad basado en Sesgo y Curtosis

\(H_0:\) 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-16

Test de Aleatoriedad

\(H_0:\) 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.sided

autoplot(modelo1)

Pronóstico Fuera Muestra

plot(forecast(modelo1,h=6, fan=T))
lines(fitted(modelo1), col="red")
lines(ts(rec, start = c(1950,1),frequency = 12),col="green")

forecast(modelo1,h=12, fan=T)

##          Point Forecast    Lo 51     Hi 51    Lo 54     Hi 54    Lo 57
## Oct 1986       86.12217 79.83613  92.40822 79.39414  92.85021 78.93569
## Nov 1986       91.71160 81.46679 101.95641 80.74644 102.67675 79.99928
## Dec 1986       96.39389 83.01340 109.77438 82.07257 110.71521 81.09672
## Jan 1987       96.20277 80.21239 112.19316 79.08805 113.31750 77.92186
## Feb 1987       94.77058 76.52110 113.02005 75.23792 114.30323 73.90697
## Mar 1987       92.32968 72.06619 112.59317 70.64139 114.01797 69.16356
## Apr 1987       88.57743 66.48135 110.67351 64.92770 112.22716 63.31622
## May 1987       84.42993 60.64166 108.21821 58.96902 109.89085 57.23412
## Jun 1987       78.23499 52.86705 103.60293 51.08334 105.38664 49.23324
## Jul 1987       70.86156 44.00668  97.71643 42.11842  99.60469 40.15987
## Aug 1987       68.11405 39.85035  96.37775 37.86303  98.36507 35.80174
## Sep 1987       75.93821 46.33266 105.54376 44.25099 107.62543 42.09183
##              Hi 57    Lo 60     Hi 60    Lo 63     Hi 63    Lo 66     Hi 66
## Oct 1986  93.30866 78.45826  93.78609 77.95877  94.28558 77.43342  94.81093
## Nov 1986 103.42392 79.22118 104.20202 78.40713 105.01607 77.55093 105.87227
## Dec 1986 111.69106 80.08046 112.70732 79.01725 113.77053 77.89898 114.88880
## Jan 1987 114.48369 76.70737 115.69818 75.43677 116.96878 74.10039 118.30516
## Feb 1987 115.63418 72.52091 117.02025 71.07080 118.47035 69.54562 119.99554
## Mar 1987 115.49580 67.62453 117.03483 66.01438 118.64497 64.32088 120.33848
## Apr 1987 113.83865 61.63800 115.51687 59.88224 117.27263 58.03558 119.11929
## May 1987 111.62575 55.42738 113.43249 53.53716 115.32271 51.54908 117.31079
## Jun 1987 107.23675 47.30652 109.16347 45.29077 111.17921 43.17067 113.29931
## Jul 1987 101.56324 38.12022 103.60289 35.98632 105.73679 33.74195 107.98116
## Aug 1987 100.42636 33.65509 102.57301 31.40924 104.81886 29.04713 107.18097
## Sep 1987 109.78459 39.84326 112.03316 37.49080 114.38563 35.01654 116.85988
##             Lo 69     Hi 69    Lo 72     Hi 72    Lo 75     Hi 75    Lo 78
## Oct 1986 76.87743  95.36692 76.28464  95.95971 75.64694  96.59741 74.95324
## Nov 1986 76.64479 106.77841 75.67868 107.74452 74.63938 108.78382 73.50881
## Dec 1986 76.71550 116.07228 75.45370 117.33408 74.09628 118.69150 72.61968
## Jan 1987 72.68606 119.71948 71.17814 121.22740 69.55596 122.84959 67.79135
## Feb 1987 67.93148 121.60967 66.21052 123.33063 64.35916 125.18199 62.34525
## Mar 1987 62.52861 122.13075 60.61772 124.04163 58.56205 126.09731 56.32588
## Apr 1987 56.08121 121.07365 53.99752 123.15735 51.75593 125.39894 49.31752
## May 1987 49.44504 119.41483 47.20177 121.65810 44.78851 124.07136 42.16336
## Jun 1987 40.92692 115.54307 38.53468 117.93531 35.96116 120.50882 33.16170
## Jul 1987 31.36668 110.35643 28.83422 112.88889 26.10986 115.61325 23.14630
## Aug 1987 26.54725 109.68085 23.88194 112.34616 21.01466 115.21344 17.89563
## Sep 1987 32.39797 119.47845 29.60612 122.27030 26.60271 125.27371 23.33561
##             Hi 78    Lo 81     Hi 81    Lo 84     Hi 84     Lo 87     Hi 87
## Oct 1986  97.2911 74.18786  98.05648 73.32740  98.91694 72.334559  99.90979
## Nov 1986 109.9144 72.26142 111.16178 70.85907 112.56413 69.240960 114.18224
## Dec 1986 120.1681 70.99049 121.79729 69.15892 123.62886 67.045544 125.74224
## Jan 1987 124.6142 65.84438 126.56117 63.65555 128.75000 61.129960 131.27559
## Feb 1987 127.1959 60.12322 129.41793 57.62515 131.91600 54.742751 134.79840
## Mar 1987 128.3335 53.85862 130.80073 51.08487 133.57449 47.884367 136.77499
## Apr 1987 127.8373 46.62714 130.52773 43.60253 133.55234 40.112581 137.04228
## May 1987 126.6965 39.26694 129.59293 36.01069 132.84918 32.253472 136.60640
## Jun 1987 123.3083 30.07293 126.39705 26.60046 129.86952 22.593738 133.87624
## Jul 1987 118.5768 19.87649 121.84662 16.20048 125.52263 11.958909 129.76420
## Aug 1987 118.3325 14.45428 121.77382 10.58543 125.64267  6.121337 130.10676
## Sep 1987 128.5408 19.73088 132.14555 15.67834 136.19808 11.002314 140.87411
##               Lo 90    Hi 90      Lo 93    Hi 93      Lo 96    Hi 96      Lo 99
## Oct 1986 71.1439162 101.1004 69.6226724 102.6217  67.420461 104.8239  62.666327
## Nov 1986 67.3004857 116.1227 64.8212081 118.6020  61.232110 122.1911  53.483965
## Dec 1986 64.5111391 128.2766 61.2730162 131.5148  56.585384 136.2024  46.465725
## Jan 1987 58.1012129 134.3043 54.2314861 138.1741  48.629520 143.7760  36.535995
## Feb 1987 51.2861096 138.2550 46.8696762 142.6715  40.476277 149.0649  26.674205
## Mar 1987 44.0462487 140.6131 39.1424159 145.5169  32.043438 152.6159  16.718164
## Apr 1987 35.9273516 141.2275 30.5800265 146.5748  22.839032 154.3158   6.127773
## May 1987 27.7477233 141.1121 21.9908811 146.8690  13.657054 155.2028  -4.334013
## Jun 1987 17.7887840 138.6812 11.6496572 144.8203   2.762420 153.7076 -16.423348
## Jul 1987  6.8723136 134.8508  0.3733439 141.3498  -9.034816 150.7579 -29.345151
## Aug 1987  0.7678965 135.4602 -6.0720127 142.3001 -15.973730 152.2018 -37.349556
## Sep 1987  5.3947111 146.4817 -1.7699317 153.6464 -12.141746 164.0182 -34.532418
##             Hi 99
## Oct 1986 109.5780
## Nov 1986 129.9392
## Dec 1986 146.3221
## Jan 1987 155.8696
## Feb 1987 162.8669
## Mar 1987 167.9412
## Apr 1987 171.0271
## May 1987 173.1939
## Jun 1987 172.8933
## Jul 1987 171.0683
## Aug 1987 173.5777
## Sep 1987 186.4088

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)

## Warning: Deprecated function. Use the `create_axis` function.

función de correlación cruzada

ccfvalues<-ccf(soi,rec,20)

Las correlaciones cruzadas más dominantes ocurren en algún lugar entre = −10 y aproximadamente = −4.

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

Las correlaciones en esta región son negativas, lo que indica que es probable que un valor de SOI superior al promedio dé lugar a un valor de “recreuitment”(cantidad de peces) inferior al promedio unos 6 meses después. Y, un SOI por debajo del promedio se asocia con un valor Recruitment probablemente superior al promedio unos 6 meses después.