Arima Estacionalidad
Ricardo R. Palma
21 de marzo de 2022
ARIMA S-ARIMA
Introducción
ARIMA es una técnica básica para el pronóstico de series de tiempo complejas , si son series estacionales se usa S-ARIMA (Seasonal). Este artículo provee una descripción paso a paso para hacer el ajuste y parametrización la biblioteca arima usando R.
Arima es modelo muy popular para simulado y pronóstico de escenarios para empresas de base tecnológica, es flexible ideal para trabajar modelos con variaciones estacionales a partir de información histórica de datos que no son relativos al producto innovador, pero que sí son complementarios o sustitutos. Especialmente para hacer predicciones en los escenarios futuros. Este modelo básico de pronóstico y sus técnicas asociadas pueden ser usados como base para modelos más complejos, como los casos en los que la empresas innovadoras lanzan acciones o fondos de crowdfunding. En este caso de estudio nosotros haremos sobre un ejemplo que examina una serie de tiempo de la demanda de bicicletas en una zona urbana (bicicletas públicas de alquiler) que están relacionadas con el lanzamiento al mercado de una app para proveer información a los dueños de negocios sobre recorridos frecuentes de potenciales clientes que pasan por sus locales. Trataremos de ajustar los parámetros de arima para que se adecuen a este modelo tomando los datos del pasado y proyectando los datos que tenemos en la actualidad (entrenamiento de modelo) para que se ajusten más tarde a la demanda real en el futuro. También generaremos una lista de revisión de los pasos a seguir para ajustar otros tipos de modelos.
ARIMA por sus siglas en inglés , acrónimo del inglés autoregressive integrated moving average , o Modelo autorregresivo integrado de media móvil es muy utilizado en estadística y econometría, en particular en series temporales, un modelo autorregresivo integrado de promedio móvil o simplemente un sistema que se comporta como ARIMA es un modelo estadístico que utiliza variaciones y regresiones de datos estadísticos con el fin de encontrar patrones para una predicción hacia el futuro. Se trata de un modelo dinámico de series temporales, es decir, las estimaciones futuras vienen explicadas por los datos del pasado y no por variables independientes. Fue desarrollado a finales de los sesenta del siglo XX. Box y Jenkins (1976)
Business as Usual
El modelo ARIMA necesita identificar los coeficientes y número de regresiones que se utilizarán. Este modelo es muy sensible a la precisión con que se determinen sus coeficientes.
Se suele expresar como ARIMA(p,d,q) donde los parámetros p, d y q son números enteros no negativos que indican el orden de las distintas componentes del modelo — respectivamente, las componentes autorregresiva, integrada y de media móvil. Cuando alguno de los tres parámetros es cero, es común omitir las letras correspondientes del acrónimo — AR para la componente autorregresiva, I para la integrada y MA para la media móvil. Por ejemplo, ARIMA(0,1,0) se puede expresar como I(1) y ARIMA(0,0,1) como MA(1).
El modelo ARIMA puede generalizarse aún más para considerar el efecto de la estacionalidad. En ese caso, se habla de un modelo SARIMA (seasonal autoregressive integrated moving average).
El modelo ARIMA (p,d,q) se puede representar como:
\[ Y_t= -(\Delta^{d} Y_{t}- Y{t}) + \phi_{0} +\sum _{i=1}^{p}\phi _{i} \Delta ^{d}Y_{t-i} -\sum _{i=1}^{q}\theta _{i}\varepsilon _{t-i}+\varepsilon _{t} \displaystyle \]
en donde \(d\) corresponde a las d diferencias que son necesarias para convertir la serie original en estacionaria, \(\phi_{1},\dot,\dot,\dot,\phi{p}\) son los parámetros pertenecientes a la parte “autorregresiva” del modelo, \({\displaystyle \theta _{1},\ldots ,\theta _{q}}{\displaystyle \theta _{1},\ldots ,\theta _{q}}\) los parámetros pertenecientes a la parte “medias móviles” del modelo.
\(\phi_{0}\) es una constante para cada modelo y \(\varepsilon _{t}\) es el término de error (llamado también innovación o perturbación estocástica esta última asociada más para modelos econométricos uniecuacionales o multiecuacionales).
Se debe tomar en cuenta que:
\[{\displaystyle \Delta Y_{t}=Y_{t}-Y_{t-1}}\]
Caso de estudio
Comenzaremos modelando los datos de una serie temporal. Usando R, convertiremos los datos disponibles en formato de datos de serie temporal. Para hacer esto, ejecutaremos el siguiente comando:
RawData=scan('https://themys.sid.uncu.edu.ar/~rpalma/R-cran/Routers.txt')
plot(RawData, main="Instalacion de Routers ITC (fibra óptica)",ylab="routers",xlab="mes")
Venta de routers de uso hogareños en América Latina
Para el análisis arima es necesario que especifiquemos la duración de la estacionalidad. ¿Se trata de una estacionalidad semanal?, ¿es por horas?, etc. Para ello especificaremos la frecuencia. En nuestro caso se trata de una estacionalidad anula, por ello la frecuencia será 12 meses. En este caso suponemos (tenemos la hypótesis) que se trata de un fenómeno de repetición anual.
tsData = ts(RawData, start = c(2005,1), frequency = 12)
print(tsData)## Jan Feb Mar Apr May Jun
## 2005 927180000 913800000 1116430000 1188880000 1194320000 1277960000
## 2006 1136610000 1082240000 1422560000 1298350000 1507350000 1495540000
## 2007 1359510000 1266150000 1466470000 1658220000 1633650000 1692940000
## 2008 154844000 143552000 171573000 188322000 192756000 195296000
## 2009 1797590000 1738210000 2113870000 2105510000 2183710000 2320570000
## 2010 1939160000 1883750000 2361870000 2490370000 2359570000 2589800000
## 2011 2250100000 2257420000 2651590000 2719860000 2909530000 2851080000
## 2012 2685780000 2560630000 3120410000 3267410000 3151570000 3530160000
## 2013 2891860000 2968810000 3025890000 3340910000 3257900000 3377820000
## 2014 3177600000 2981880000 3634290000 3502030000 3721490000 3718770000
## 2015 3522000000 3349380000 3728910000 3973880000 3856570000 4169610000
## 2016 3633670000 3429790000 3849360000 4217180000 4028770000 4276150000
## Jul Aug Sep Oct Nov Dec
## 2005 1589430000 1780130000 1433850000 1271790000 1144030000 1249000000
## 2006 1857920000 2017580000 1665650000 1480480000 1315810000 1413150000
## 2007 2155380000 2334270000 1844020000 1784320000 1551790000 1633550000
## 2008 252288000 2683790000 2188100000 2035450000 1721480000 1983810000
## 2009 2941730000 3127000000 2518910000 235560000 202876000 224383000
## 2010 3210850000 3345620000 2769320000 2582690000 2335320000 2517550000
## 2011 3626870000 3863470000 3142050000 2921240000 2617400000 2918100000
## 2012 4036620000 4510980000 3568110000 3525660000 3055800000 4106140000
## 2013 4232970000 4541720000 3537270000 3534130000 3152720000 3419020000
## 2014 4724580000 4855170000 4062230000 3772620000 3297940000 3843500000
## 2015 4924800000 5122090000 4115140000 3923800000 3696710000 4002430000
## 2016 5382540000 5280070000plot(tsData,main="Comienzo de contrato con Huawei",ylab="Equipos OptiX PTN 905E", xlab="Año")
Podemos deducir del gráfico en sí mismo que los puntos de datos siguen un patrón de subidas y caídas en una tendencia ascendente. Ahora necesitamos hacer un análisis para averiguar la no estacionariedad exacta y la estacionalidad en los datos.
Pre-Análisis Exploratorio de Datos (PAED)}
Antes de realizar cualquiera de estos pasos en los datos, debemos comprender los tres componentes de los datos de una serie temporal:
En R podemos averiguar estos componentes de la serie temporal con los siguientes comandos:
componentes.ts = decompose(tsData)
plot(componentes.ts)
A continuación, debemos eliminar la parte no estacionaria de ARIMA. Para facilitar el análisis aquí, también eliminaremos la parte estacional de los datos. La parte estacional puede eliminarse del análisis y agregarse más tarde, o puede tratarse en el propio modelo ARIMA.
Para conseguir estacionariedad:
diferencia los datos – calcula las diferencias entre observaciones consecutivas, aplica el logaritmo o raíz cuadrada a los datos de la serie para estabilizar la varianza no constante.
Si los datos contienen una tendencia, ajusta algún tipo de curva a los datos y luego modela los residuales de ese ajuste.
Prueba de raíz unitaria: esta prueba se usa para descubrir la primera diferencia o regresión que se debe usar en los datos de tendencias para hacerla estacionaria. En la prueba, llamada de Kwiatkowski-Phillips-Schmidt-Shin (KPSS), los valores pequeños de p sugieren que se requiera una diferenciación.
El código R para la prueba de raíz unitaria:
#install.packages("fUnitRoots")
library("fUnitRoots")## Loading required package: timeDate## Loading required package: timeSeries## Loading required package: fBasicsurkpssTest(tsData, type = c("tau"), lags = c("short"),use.lag = NULL, doplot = TRUE)
##
## Title:
## KPSS Unit Root Test
##
## Test Results:
## NA
##
## Description:
## Thu Mar 24 13:30:49 2022 by user:Las barras largas (piernas) en los gráficos generales y parciales de ACF nos dicen en que lugares tenemos que tener cuidado si utilizamos el podelo para pronosticar.
Después eliminamos la no estacionariedad: (tendencia lineal creciente). Esto nos permite encontrar que causa o acciones (como campañas de marketing) han tenido efectos expansivos, pero no sostenibles en el tiempo. En nuestro caso de estudio aparece claramente la crisis financiera mundial del 2008/09 y el cierre de las importaciones.
tsstationary = diff(tsData, differences=1)
plot(tsstationary)
Hay diversos gráficos y funciones que nos pueden ayudar a detectar la estacionalidad:
Codigo R para calcular la autocorrelación:
acf(tsData,lag.max=140)
Este gráfico nos muestra las acciones que AFC tuvo que realizar para poder hacer el gráfico sin tendencia de crecimiento.
La función de autocorrelación proporciona la autocorrelación en todos los retrasos posibles. La autocorrelación en el retraso 0 se incluye por defecto, lo que siempre toma el valor 1, ya que representa la correlación entre los datos y ellos mismos. Como podemos observar en el gráfico anterior, la autocorrelación continúa disminuyendo a medida que aumenta el retraso, lo que confirma que no existe una asociación lineal entre las observaciones separadas por retrasos mayores.
timeseriesseasonallyadjusted <- tsData
tsstationary <- diff(timeseriesseasonallyadjusted, differences=1)
tsstationary## Jan Feb Mar Apr May Jun
## 2005 -13380000 202630000 72450000 5440000 83640000
## 2006 -112390000 -54370000 340320000 -124210000 209000000 -11810000
## 2007 -53640000 -93360000 200320000 191750000 -24570000 59290000
## 2008 -1478706000 -11292000 28021000 16749000 4434000 2540000
## 2009 -186220000 -59380000 375660000 -8360000 78200000 136860000
## 2010 1714777000 -55410000 478120000 128500000 -130800000 230230000
## 2011 -267450000 7320000 394170000 68270000 189670000 -58450000
## 2012 -232320000 -125150000 559780000 147000000 -115840000 378590000
## 2013 -1214280000 76950000 57080000 315020000 -83010000 119920000
## 2014 -241420000 -195720000 652410000 -132260000 219460000 -2720000
## 2015 -321500000 -172620000 379530000 244970000 -117310000 313040000
## 2016 -368760000 -203880000 419570000 367820000 -188410000 247380000
## Jul Aug Sep Oct Nov Dec
## 2005 311470000 190700000 -346280000 -162060000 -127760000 104970000
## 2006 362380000 159660000 -351930000 -185170000 -164670000 97340000
## 2007 462440000 178890000 -490250000 -59700000 -232530000 81760000
## 2008 56992000 2431502000 -495690000 -152650000 -313970000 262330000
## 2009 621160000 185270000 -608090000 -2283350000 -32684000 21507000
## 2010 621050000 134770000 -576300000 -186630000 -247370000 182230000
## 2011 775790000 236600000 -721420000 -220810000 -303840000 300700000
## 2012 506460000 474360000 -942870000 -42450000 -469860000 1050340000
## 2013 855150000 308750000 -1004450000 -3140000 -381410000 266300000
## 2014 1005810000 130590000 -792940000 -289610000 -474680000 545560000
## 2015 755190000 197290000 -1006950000 -191340000 -227090000 305720000
## 2016 1106390000 -102470000Para eliminar la estacionalidad de los datos, restamos el componente estacional de la serie original(media cero) y luego lo diferenciamos para que sea estacionario (varianza constante e independiente).
Después de eliminar la estacionalidad y hacer que los datos sean estacionarios, se verá así:
plot(tsstationary,xlab="Mes/año",ylab="Desvio Facturación")
El suavizado se hace generalmente para ayudarnos a ver mejor los patrones, las tendencias en las series temporales. Generalmente alisa las irregularidades para ver una señal más clara. Para datos estacionales, podríamos suavizar la estacionalidad para que podamos identificar la tendencia. El suavizado no nos proporciona un modelo, pero puede ser un buen primer paso para describir varios componentes de la serie.
Para suavizar las series de tiempo:
Ajuste del modelo
Es posible ajustar el modelo para que mejore los pronósticos. El modelo está fuertemente influenciado por la calidad de los datos de entrada, pero aún siendo malos la estacionalidad permite mejorar al modelo para eliminar tanto como se a posible la influencia de los residuals. Una vez que tenemos los datos listos y han satisfecho todas las suposiciones del modelo, para determinar el orden del modelo que se ajustará a los datos, necesitamos tres variables: p, d y q que son enteros positivos que se refieren al orden de las partes medias autorregresivas, integradas y móviles del modelo, respectivamente.
Para identificar que valores de p y q serán apropiados, necesitamos ejecutar las funciones acf() y pacf()
en el retardo k-esimo es la función de autocorrelación que describe la correlación entre todos los puntos de datos que están exactamente k pasos antes, después de tener en cuenta su correlación con los datos entre esos k pasos. Ayuda a identificar el número de coeficientes de autorregresión (AR) (valor p) en un modelo .
El código R para ejecutar los comandos :
acf(tsstationary, lag.max=70)
pacf(tsstationary, lag.max=140)
La forma de acf () para definir valores de p y q es mirar los gráficos y repasando la tabla podemos determinar qué tipo de modelo seleccionar y cuáles serán los valores de p, d y q.
Presentación de tablas
Como tenemos gran cantidad de datos en tablas exploraremos formas de mejorar la apariencia, tanto en html como en pdf. Es por ello que utilizaremos la biblioteca kableExtra.
En rigos esto no pertenece al modelo ARIMA, por lo que puedes no utilizar estas formas de embellecer las tablas o ajustar para que se encuadren dentro de una sola página en los PDF.
library(kableExtra)Contendio de los titulos y mensajes de la tabla.
# library(kableExtra)
t <- c(
"Forma" ,
"Modelo",
"Serie exponencial decayendo a 0 "
,
"Modelo Auto Regresivo (AR). Función pacf () que se
utilizará para identificar el orden del modelo"
,
"Picos alternativos positivos y negativos, decayendo a 0"
,
"Modelo Auto Regresivo (AR). Función pacf () que se
utilizará para identificar el orden del modelo"
,
"Uno o más picos en serie, resto todos son 0"
,
"Modelo de media móvil (MA), identifica el orden donde el gráfico se convierte en 0"
,
"Después de algunos retrasos en general la serie va decayendo."
,
"Modelo mezclado AR & MA"
,
"La serie total es 0 o casi 0"
,
"Datos aleatorios "
,
"Valores medios a intervalos fijos"
,
"Necesitamos incluir el término AR de estacionalidad"
,
"Picos visibles que no decaen a 0"
,
"Series no son estacionarias"
)
Tt <- matrix(t,8,2,byrow = TRUE)Impresión de talba con escala normal
kbl(Tt)| Forma | Modelo |
| Serie exponencial decayendo a 0 | Modelo Auto Regresivo (AR). Función pacf () que se utilizará para identificar el orden del modelo |
| Picos alternativos positivos y negativos, decayendo a 0 | Modelo Auto Regresivo (AR). Función pacf () que se utilizará para identificar el orden del modelo |
| Uno o más picos en serie, resto todos son 0 | Modelo de media móvil (MA), identifica el orden donde el gráfico se convierte en 0 |
| Después de algunos retrasos en general la serie va decayendo. | Modelo mezclado AR & MA |
| La serie total es 0 o casi 0 | Datos aleatorios |
| Valores medios a intervalos fijos | Necesitamos incluir el término AR de estacionalidad |
| Picos visibles que no decaen a 0 | Series no son estacionarias |
kbl(Tt, booktabs = T, linesep = "") %>%
kable_styling(latex_options = "striped", stripe_index = c(1,2, 5:6))| Forma | Modelo |
| Serie exponencial decayendo a 0 | Modelo Auto Regresivo (AR). Función pacf () que se utilizará para identificar el orden del modelo |
| Picos alternativos positivos y negativos, decayendo a 0 | Modelo Auto Regresivo (AR). Función pacf () que se utilizará para identificar el orden del modelo |
| Uno o más picos en serie, resto todos son 0 | Modelo de media móvil (MA), identifica el orden donde el gráfico se convierte en 0 |
| Después de algunos retrasos en general la serie va decayendo. | Modelo mezclado AR & MA |
| La serie total es 0 o casi 0 | Datos aleatorios |
| Valores medios a intervalos fijos | Necesitamos incluir el término AR de estacionalidad |
| Picos visibles que no decaen a 0 | Series no son estacionarias |
Tabla con caption
Como regla general las tablas deben llevar un título (caption o copete) que indique que es lo que nos están mostrando.
kbl(Tt, caption = "¿Qué paramétrica usar?", booktabs = T) %>%
kable_styling(latex_options = c("striped", "hold_position"))| Forma | Modelo |
| Serie exponencial decayendo a 0 | Modelo Auto Regresivo (AR). Función pacf () que se utilizará para identificar el orden del modelo |
| Picos alternativos positivos y negativos, decayendo a 0 | Modelo Auto Regresivo (AR). Función pacf () que se utilizará para identificar el orden del modelo |
| Uno o más picos en serie, resto todos son 0 | Modelo de media móvil (MA), identifica el orden donde el gráfico se convierte en 0 |
| Después de algunos retrasos en general la serie va decayendo. | Modelo mezclado AR & MA |
| La serie total es 0 o casi 0 | Datos aleatorios |
| Valores medios a intervalos fijos | Necesitamos incluir el término AR de estacionalidad |
| Picos visibles que no decaen a 0 | Series no son estacionarias |
Tabla escalada
Este efecto dde escalar las tablas es inapropiado para el formato html, pero se observan muy bien en PDF. La tipografía no pierde nitidez.
kbl(rbind(Tt,Tt,Tt), booktabs = T) %>%
kable_styling(latex_options = c("striped", "scale_down")) %>%
kable_styling(font_size = 7)| Forma | Modelo |
| Serie exponencial decayendo a 0 | Modelo Auto Regresivo (AR). Función pacf () que se utilizará para identificar el orden del modelo |
| Picos alternativos positivos y negativos, decayendo a 0 | Modelo Auto Regresivo (AR). Función pacf () que se utilizará para identificar el orden del modelo |
| Uno o más picos en serie, resto todos son 0 | Modelo de media móvil (MA), identifica el orden donde el gráfico se convierte en 0 |
| Después de algunos retrasos en general la serie va decayendo. | Modelo mezclado AR & MA |
| La serie total es 0 o casi 0 | Datos aleatorios |
| Valores medios a intervalos fijos | Necesitamos incluir el término AR de estacionalidad |
| Picos visibles que no decaen a 0 | Series no son estacionarias |
| Forma | Modelo |
| Serie exponencial decayendo a 0 | Modelo Auto Regresivo (AR). Función pacf () que se utilizará para identificar el orden del modelo |
| Picos alternativos positivos y negativos, decayendo a 0 | Modelo Auto Regresivo (AR). Función pacf () que se utilizará para identificar el orden del modelo |
| Uno o más picos en serie, resto todos son 0 | Modelo de media móvil (MA), identifica el orden donde el gráfico se convierte en 0 |
| Después de algunos retrasos en general la serie va decayendo. | Modelo mezclado AR & MA |
| La serie total es 0 o casi 0 | Datos aleatorios |
| Valores medios a intervalos fijos | Necesitamos incluir el término AR de estacionalidad |
| Picos visibles que no decaen a 0 | Series no son estacionarias |
| Forma | Modelo |
| Serie exponencial decayendo a 0 | Modelo Auto Regresivo (AR). Función pacf () que se utilizará para identificar el orden del modelo |
| Picos alternativos positivos y negativos, decayendo a 0 | Modelo Auto Regresivo (AR). Función pacf () que se utilizará para identificar el orden del modelo |
| Uno o más picos en serie, resto todos son 0 | Modelo de media móvil (MA), identifica el orden donde el gráfico se convierte en 0 |
| Después de algunos retrasos en general la serie va decayendo. | Modelo mezclado AR & MA |
| La serie total es 0 o casi 0 | Datos aleatorios |
| Valores medios a intervalos fijos | Necesitamos incluir el término AR de estacionalidad |
| Picos visibles que no decaen a 0 | Series no son estacionarias |
Tabla a lo ancho
Del mismo modo no se aprecia en html, pero este comando puede compactar una tabla muy grande para que se imprima en PDF en una sola página en forma apaisada.
kbl(Tt, booktabs = T) %>%
kable_styling(full_width = T) %>%
column_spec(1, width = "8cm")| Forma | Modelo |
| Serie exponencial decayendo a 0 | Modelo Auto Regresivo (AR). Función pacf () que se utilizará para identificar el orden del modelo |
| Picos alternativos positivos y negativos, decayendo a 0 | Modelo Auto Regresivo (AR). Función pacf () que se utilizará para identificar el orden del modelo |
| Uno o más picos en serie, resto todos son 0 | Modelo de media móvil (MA), identifica el orden donde el gráfico se convierte en 0 |
| Después de algunos retrasos en general la serie va decayendo. | Modelo mezclado AR & MA |
| La serie total es 0 o casi 0 | Datos aleatorios |
| Valores medios a intervalos fijos | Necesitamos incluir el término AR de estacionalidad |
| Picos visibles que no decaen a 0 | Series no son estacionarias |
Tablas centradas
El siguiente comando tratará de incrustar la tabla PDF a una distancia equidistante de encabezado y el pie de página y distribuirá el texto arriba y abajo de la misma.
kbl(Tt, booktabs = T) %>%
kable_styling(position = "center")| Forma | Modelo |
| Serie exponencial decayendo a 0 | Modelo Auto Regresivo (AR). Función pacf () que se utilizará para identificar el orden del modelo |
| Picos alternativos positivos y negativos, decayendo a 0 | Modelo Auto Regresivo (AR). Función pacf () que se utilizará para identificar el orden del modelo |
| Uno o más picos en serie, resto todos son 0 | Modelo de media móvil (MA), identifica el orden donde el gráfico se convierte en 0 |
| Después de algunos retrasos en general la serie va decayendo. | Modelo mezclado AR & MA |
| La serie total es 0 o casi 0 | Datos aleatorios |
| Valores medios a intervalos fijos | Necesitamos incluir el término AR de estacionalidad |
| Picos visibles que no decaen a 0 | Series no son estacionarias |
Aplicaremos estos modelos de tablas al caso de ARIMA que venimos desarrollando.
Ttt <- data.frame(Tt)
kbl(Ttt, booktabs = T) %>%
kable_styling(full_width = F) %>%
column_spec(1, bold = T, color = "red") %>%
column_spec(2, width = "30em")| X1 | X2 |
|---|---|
| Forma | Modelo |
| Serie exponencial decayendo a 0 | Modelo Auto Regresivo (AR). Función pacf () que se utilizará para identificar el orden del modelo |
| Picos alternativos positivos y negativos, decayendo a 0 | Modelo Auto Regresivo (AR). Función pacf () que se utilizará para identificar el orden del modelo |
| Uno o más picos en serie, resto todos son 0 | Modelo de media móvil (MA), identifica el orden donde el gráfico se convierte en 0 |
| Después de algunos retrasos en general la serie va decayendo. | Modelo mezclado AR & MA |
| La serie total es 0 o casi 0 | Datos aleatorios |
| Valores medios a intervalos fijos | Necesitamos incluir el término AR de estacionalidad |
| Picos visibles que no decaen a 0 | Series no son estacionarias |
El modelo que hemos aplanado con AFC, ahora nos permitiria trabajar con proyecciones como si se tratase de una regresión lineal. Luego Arima aplicará los parámetros \(p,q,d\) para regenerar los datos pero ahora fuera del intervalo en el con el que construimos el modelo. Esto se denomina extrapolar el modelo linealmente y reconstruir su no linealidad.
fitARIMA <- arima(tsData, order=c(1,1,1),seasonal =list(order = c(1,0,0), period = 12),method="ML")
#install.packages("lmtest")
library(lmtest)## Loading required package: zoo##
## Attaching package: 'zoo'## The following object is masked from 'package:timeSeries':
##
## time<-## The following objects are masked from 'package:base':
##
## as.Date, as.Date.numericcoeftest(fitARIMA)##
## z test of coefficients:
##
## Estimate Std. Error z value Pr(>|z|)
## ar1 -0.795254 0.108606 -7.3224 2.436e-13 ***
## ma1 0.648118 0.129662 4.9985 5.777e-07 ***
## sar1 0.484253 0.074262 6.5209 6.991e-11 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1fitARIMA <- arima(tsData, order=c(1,1,1), seasonal =list(order = c(1,0,0), period = 12),method="ML")El proceso de ajuste es un proceso recursivo y necesitamos ejecutar esta función arima () con diferentes valores (p, d, q) para encontrar el modelo más optimizado y eficiente.
La salida de fitarima () incluye los coeficientes ajustados y el error estándar para cada coeficiente. Observando los coeficientes podemos excluir los insignificantes. Podemos usar una función confint () para este propósito. Podemos usar una función confint () para este propósito.
Intervalos de confianza de los modelos ajustados
confint(fitARIMA)## 2.5 % 97.5 %
## ar1 -1.0081178 -0.5823911
## ma1 0.3939846 0.9022509
## sar1 0.3387016 0.6298038Elección del mejor modelo
R utiliza la estimación de máxima verosimilitud (MLE) para estimar el modelo ARIMA. Intenta maximizar la probabilidad logarítmica para valores dados de p, d y q al encontrar estimaciones de parámetros para maximizar la probabilidad de obtener los datos que hemos observado. Averigüa el Criterio de Información de Akaike (AIC) para un conjunto de modelos e investigarlos modelos con los valores AIC más bajos. Prueba el criterio de información bayesiano (BIC) de Schwarz e investiga los modelos con los valores BIC más bajos. Al estimar los parámetros del modelo utilizando la estimación de máxima verosimilitud, es posible aumentar la probabilidad agregando parámetros adicionales, lo que puede resultar en un ajuste excesivo. El BIC resuelve este problema introduciendo un término de penalización para el número de parámetros en el modelo. Junto con AIC y BIC, también debemos observar de cerca los valores de los coeficientes y debemos decidir si incluir ese componente o no de acuerdo con su nivel de significación.
Medidas de diagnostico
Deben observarse los coeficientes de autocorrelación muestral de los residuos y comprobar que ninguno
de ellos supera el valor de las bandas de significatividad al \(5\% (±1,96(1/T ))^{1/2}\). El valor T es una aproximación de la varianza asintótica pero resulta sólo adecuada para valores grandes de «j». Se
aconseja,por tanto, utilizar distinta amplitud de bandas como por ejemplo \(±(1/T )^{1/2}\) para los términos más cercanos a cero.Intenta averiguar el patrón en los residuos del modelo elegido trazando el ACF de los residuos y haciendo una prueba de portmanteau. Necesitamos probar modelos modificados si la gráfica no parece un ruido blanco. Una vez que los residuos parezcan ruido blanco, calcularemos los pronósticos.
Prueba de Box-Ljung
Es una prueba para comprobar si una serie de observaciones en un período de tiempo especifico son aleatorias o independientes en todos los retardos hasta el especificado. En lugar de probar la aleatoriedad en cada retardo distinto, prueba la aleatoriedad «general» basada en un número de retardos y, por lo tanto, es una prueba de comparación. Se aplica a los residuos de un modelo ARIMA ajustado, no a la serie original, y en tales aplicaciones, la hipótesis que se está probando realmente es que los residuos del modelo ARIMA no tienen autocorrelación. Si los retardos no son independientes, un retardo puede estar relacionado con otros retardos k unidades de tiempo después por lo que la autocorrelacion puede reducir la exactitud de un modelo predictivo basado en el tiempo y conducir a una interpretación errónea de los datos. El código R para obtener los resultados de la prueba Box (pruebas de muestreo Montecarlo de cuantil vs cuantiles muestreados):
#install.packages("FitAR")
library(FitAR)## Loading required package: lattice## Loading required package: leaps## Loading required package: ltsa## Loading required package: bestglmacf(fitARIMA$residuals,lag.max=140)
boxresult=LjungBoxTest (fitARIMA$residuals,k=2,StartLag=1)
plot(boxresult[,3],main= "Ljung-Box Q Test", ylab= "P-values", xlab= "Lag")
qqnorm(fitARIMA$residuals)
qqline(fitARIMA$residuals)
Como todos los gráficos apoyan el supuesto de que no hay un patrón en los residuales, podemos seguir adelante y calcular el pronóstico.
Función auto.arima ():
El paquete de pronóstico proporciona dos funciones: ets () y auto.arima () para la selección automática de modelos exponenciales y modelos ARIMA. La función auto.arima () en R utiliza una combinación de pruebas de raíz unitaria, minimización del AIC y MLE para obtener un modelo ARIMA. La prueba KPSS se usa para determinar el número de diferencias (d) En el algoritmo Hyndman- Khandakar para el modelado ARIMA automático. Los p, d y q luego se eligen minimizando el AICc. El algoritmo utiliza una búsqueda por pasos para recorrer el espacio del modelo para seleccionar el mejor modelo con el AICc más pequeño. Si d = 0, entonces se incluye la constante c; si d≥1 entonces la constante c se pone a cero. Las variaciones en el modelo actual se consideran variando p y / o q del modelo actual en ± 1 e incluyendo / excluyendo c del modelo actual. El mejor modelo considerado hasta ahora (ya sea el modelo actual o una de estas variaciones) se convierte en el nuevo modelo actual. Ahora, este proceso se repite hasta que no se pueda encontrar un AIC inferior.
#install.packages("forecast")
library(forecast)## Warning: replacing previous import 'lifecycle::last_warnings' by
## 'rlang::last_warnings' when loading 'tibble'## Warning: replacing previous import 'lifecycle::last_warnings' by
## 'rlang::last_warnings' when loading 'pillar'## Registered S3 method overwritten by 'quantmod':
## method from
## as.zoo.data.frame zoo##
## Attaching package: 'forecast'## The following object is masked from 'package:FitAR':
##
## BoxCoxauto.arima(tsData, trace=TRUE)##
## ARIMA(2,1,2)(1,0,1)[12] with drift : Inf
## ARIMA(0,1,0) with drift : 5971.076
## ARIMA(1,1,0)(1,0,0)[12] with drift : Inf
## ARIMA(0,1,1)(0,0,1)[12] with drift : Inf
## ARIMA(0,1,0) : 5969.545
## ARIMA(0,1,0)(1,0,0)[12] with drift : Inf
## ARIMA(0,1,0)(0,0,1)[12] with drift : Inf
## ARIMA(0,1,0)(1,0,1)[12] with drift : Inf
## ARIMA(1,1,0) with drift : Inf
## ARIMA(0,1,1) with drift : Inf
## ARIMA(1,1,1) with drift : Inf
##
## Best model: ARIMA(0,1,0)## Series: tsData
## ARIMA(0,1,0)
##
## sigma^2 estimated as 2.586e+17: log likelihood=-2983.76
## AIC=5969.52 AICc=5969.55 BIC=5972.45Predecir utilizando un modelo ARIMA
Los parámetros del modelo ARIMA se pueden usar como un modelo predictivo para hacer pronósticos de valores futuros de las series temporales, una vez que se selecciona el modelo más adecuado para los datos de la serie temporal. El valor d afecta a los intervalos de predicción; los intervalos de predicción aumentan de tamaño con valores más altos de «d». Todos los intervalos de predicción serán esencialmente los mismos cuando d = 0 porque la desviación estándar de la previsión a largo plazo irá a la desviación estándar de los datos históricos. Existe una función llamada predict () que se usa para predicciones a partir de los resultados de varias funciones de ajuste de modelos. Se necesita un argumento n.ahead () que especifica la cantidad de tiempo que se adelanta para predecir.
predict(fitARIMA,n.ahead = 5)## $pred
## Jan Feb Mar Apr May Jun Jul Aug Sep Oct Nov
## 2016 4877293715 4717165869 4660853489
## 2017 4621589580
## Dec
## 2016 4766228602
## 2017
##
## $se
## Jan Feb Mar Apr May Jun Jul Aug Sep Oct Nov
## 2016 442244589 581240617 722367284
## 2017 921409972
## Dec
## 2016 819865541
## 2017La función forecast() en el paquete R de forecast también se puede utilizar para pronosticar valores futuros de la serie temporal. Aquí también podemos especificar el nivel de confianza para los intervalos de predicción utilizando el argumento de nivel.
futurVal <- forecast(fitARIMA,h=10, level=c(99.5))
plot(futurVal)
Los pronósticos se muestran como una línea azul, con los intervalos de predicción del 80% como un área sombreada oscura, y los intervalos de predicción del 95% como un área sombreada clara. Este es el proceso general analizar para analizar datos de series temporales y pronosticar valores de las series existentes utilizando ARIMA.
Referencias:
Time Series Analysis Using ARIMA Model In R https://datascienceplus.com/time-series-analysis-using-arima-model-in-r/