ARIMA S-ARIMA

Introducción

ARIMA es una técnica básica para el pronóstico de series de tiempo complejas, basado en inteligencia artificial. 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 de 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 adecuén 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

Estudiaremos el caso de un producto con fuerte demanda estacional. En nuestro caso la demanda de cerveza en Australia.

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)

Venta de pallets de cerceza.

tsData = ts(RawData, start = c(2006,1), frequency = 12)
print(tsData)
##             Jan        Feb        Mar        Apr        May        Jun
## 2006  927180000  913800000 1116430000 1188880000 1194320000 1277960000
## 2007 1136610000 1082240000 1422560000 1298350000 1507350000 1495540000
## 2008 1359510000 1266150000 1466470000 1658220000 1633650000 1692940000
## 2009  154844000  143552000  171573000  188322000  192756000  195296000
## 2010 1797590000 1738210000 2113870000 2105510000 2183710000 2320570000
## 2011 1939160000 1883750000 2361870000 2490370000 2359570000 2589800000
## 2012 2250100000 2257420000 2651590000 2719860000 2909530000 2851080000
## 2013 2685780000 2560630000 3120410000 3267410000 3151570000 3530160000
## 2014 2891860000 2968810000 3025890000 3340910000 3257900000 3377820000
## 2015 3177600000 2981880000 3634290000 3502030000 3721490000 3718770000
## 2016 3522000000 3349380000 3728910000 3973880000 3856570000 4169610000
## 2017 3633670000 3429790000 3849360000 4217180000 4028770000 4276150000
##             Jul        Aug        Sep        Oct        Nov        Dec
## 2006 1589430000 1780130000 1433850000 1271790000 1144030000 1249000000
## 2007 1857920000 2017580000 1665650000 1480480000 1315810000 1413150000
## 2008 2155380000 2334270000 1844020000 1784320000 1551790000 1633550000
## 2009  252288000 2683790000 2188100000 2035450000 1721480000 1983810000
## 2010 2941730000 3127000000 2518910000  235560000  202876000  224383000
## 2011 3210850000 3345620000 2769320000 2582690000 2335320000 2517550000
## 2012 3626870000 3863470000 3142050000 2921240000 2617400000 2918100000
## 2013 4036620000 4510980000 3568110000 3525660000 3055800000 4106140000
## 2014 4232970000 4541720000 3537270000 3534130000 3152720000 3419020000
## 2015 4724580000 4855170000 4062230000 3772620000 3297940000 3843500000
## 2016 4924800000 5122090000 4115140000 3923800000 3696710000 4002430000
## 2017 5382540000 5280070000
plot(tsData)

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)

En este gráfico podemos observar 4 sub-graficos:

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 Kwiatkowski-Phillips-Schmidt-Shin (KPSS), los valores pequeños de p sugieren que se requiera una diferenciación.

El codigo R para la prueba de raíz unitaria:

#install.packages("fUnitRoots")
library("fUnitRoots")
urkpssTest(tsData, type = c("tau"), lags = c("short"),use.lag = NULL, doplot = TRUE)

## 
## Title:
##  KPSS Unit Root Test
## 
## Test Results:
##   NA
## 
## Description:
##  Wed Jun  7 16:20:12 2023 by user:

Después eliminamos la no estacionariedad: (tendencia lineal creciente)

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)

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-componentes.ts$seasonal
tsstationary <- diff(timeseriesseasonallyadjusted, differences=1)

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

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

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

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.

library(kableExtra)
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)
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

kbl(Tt, caption = "¿Qué paramétrica usar?", booktabs = T) %>%
kable_styling(latex_options = c("striped", "hold_position"))
¿Qué paramétrica usar?
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

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

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

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 
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 
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 objects are masked from 'package:base':
## 
##     as.Date, as.Date.numeric
coeftest(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 ' ' 1
fitARIMA <- 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.

confint(fitARIMA)
##           2.5 %     97.5 %
## ar1  -1.0081178 -0.5823911
## ma1   0.3939846  0.9022509
## sar1  0.3387016  0.6298038

Elecció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:

#install.packages("FitAR")
#library(FitAR)
#acf(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)
## Registered S3 method overwritten by 'quantmod':
##   method            from
##   as.zoo.data.frame zoo
auto.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 = 2.586e+17:  log likelihood = -2983.76
## AIC=5969.52   AICc=5969.55   BIC=5972.45

Predecir 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 = 25)
## $pred
##             Jan        Feb        Mar        Apr        May        Jun
## 2017                                                                  
## 2018 4621589580 4495873982 4720512743 4881563718 4803898156 4912899005
## 2019 5085922151 5023317346 5133472426 5210369783 5173628451 5225721786
##             Jul        Aug        Sep        Oct        Nov        Dec
## 2017                       4877293715 4717165869 4660853489 4766228602
## 2018 5457254970 5400807420 5211190472 5129331056 5105494810 5153792743
## 2019 5489876839 5462105194 5370630020                                 
## 
## $se
##             Jan        Feb        Mar        Apr        May        Jun
## 2017                                                                  
## 2018  921409972 1002298783 1084880949 1155826294 1226992805 1290977682
## 2019 1957329030 2046140721 2134457223 2216808199 2298088879 2375143085
##             Jul        Aug        Sep        Oct        Nov        Dec
## 2017                        442244589  581240617  722367284  819865541
## 2018 1354418513 1413125990 1544150757 1651889306 1762718038 1859667269
## 2019 2450890506 2523504017 2620829281

La 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/