ARIMA

Resumen para ajuste de un modelo ARIMA

• Preprocesar y explorar la serie temporal.
• Verificar estacionariedad y aplicar transformaciones si es necesario.
• Identificar los parámetros del modelo (\(p\), \(d\), \(q\)) usando ACF y PACF.
• Ajustar el modelo ARIMA.
• Evaluar el modelo con métricas como RMSE y MAE.
• Diagnóstico revisando residuos.
• Realizar pronósticos.

Este procedimiento es esencial para construir un modelo ARIMA robusto y hacer predicciones confiables a partir de series temporales.

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from statsmodels.graphics.tsaplots import quarter_plot
from statsmodels.tsa.stattools import acf, pacf
from statsmodels.graphics.tsaplots import plot_acf
from statsmodels.graphics.tsaplots import plot_pacf
import statsmodels.tsa.stattools as sts
from statsmodels.tsa.seasonal import seasonal_decompose
from statsmodels.tsa.stattools import adfuller
from statsmodels.tsa.arima_model import ARMA
from statsmodels.tsa.arima_model import ARIMA
from statsmodels.tsa.statespace.sarimax import SARIMAX
from statsmodels.tsa.statespace.sarimax import SARIMAXResults
from scipy.stats.distributions import chi2
#!pip install pmdarima
import pmdarima as pm
import matplotlib.pyplot as plt
import seaborn as sns
import glob
sns.set()
import warnings
warnings.filterwarnings('ignore')

data =  pd.read_csv("./JohnsonJohnson.csv")
data.head()

##          date  sales
## 0  01/01/1960   0.71
## 1  01/04/1960   0.63
## 2  01/07/1960   0.85
## 3  01/10/1960   0.44
## 4  01/01/1961   0.61

data.date.describe()

## count             84
## unique            84
## top       01/01/1960
## freq               1
## Name: date, dtype: object

Convertimos a formato datetime

data['date'] = pd.to_datetime(data['date'], format="%d/%m/%Y")
data['sales'] = pd.to_numeric(data['sales'])
# fijamos índice
data = data.set_index('date')
data.head()

##             sales
## date             
## 1960-01-01   0.71
## 1960-04-01   0.63
## 1960-07-01   0.85
## 1960-10-01   0.44
## 1961-01-01   0.61

Preprocesamiento

Ajustamos la frecuencia de la serie

Fuente: pandas-docs_time series

yt = data.asfreq(freq='QS') # trimestral
yt

##             sales
## date             
## 1960-01-01   0.71
## 1960-04-01   0.63
## 1960-07-01   0.85
## 1960-10-01   0.44
## 1961-01-01   0.61
## ...           ...
## 1979-10-01   9.99
## 1980-01-01  16.20
## 1980-04-01  14.67
## 1980-07-01  16.02
## 1980-10-01  11.61
## 
## [84 rows x 1 columns]

División de la serie

Dividimos los datos en conjunto de entrenamiento y conjunto de prueba. Usualmente se divide entre 80-90% para entrenamiento (train), y 20-10% respectivamente para prueba (test)

def train_test_split(data, train_size=0.85):
  size = int(len(data)*train_size)
  # split data
  X_train = data.iloc[:size]
  X_test = data.iloc[size:]
  return X_train, X_test  

# split
X_train, X_test = train_test_split(yt, train_size=0.85)

print(f"Train : {len(X_train)} \n  Test : {len(X_test)} \n Total : {len(yt)}")

## Train : 71 
##   Test : 13 
##  Total : 84

Graficamos la serie

plt.figure(figsize=(12,5))
plt.plot(X_train, label='Train')
plt.plot(X_test, label = 'Test')
plt.title('Serie Jhonson & Jhonson')
plt.xlabel('date')
plt.ylabel('Sales')
plt.legend()
plt.show()

Podemos observar que la serie no es estacionaria, es decir, su media no es cero, y su varianza no es constante, de hecho va en aumentando a medida que aumenta el tiempo. Si graficamos la función de autocorrelación muestral los rezagos tendrían que decaer de forma exponencial, característica que rebela la no estacionariedad de la serie.

El siguiente gráfico se gráfican los trimestres por separado, con líneas rojas horizontales, las cuales representan los promedios de cada trimestre.

fig = quarter_plot(X_train.sales);
fig.set_size_inches((10, 5))
plt.show()

En el gráfico anterior podemos apreciar que en el cuatrimestre dos \(q2\) y tres \(q3\), los beneficios trimestrales son mayores en comparación con los cuatrimestre uno \(q1\) y dos \(q2\), lo cual muestran un patrón estacional.

Los siguientes gráficos puede ser útiles para darnos una idea de la existencia de tendencia, y estacionalidad principalmente. En la parte superior observamos la serie original, el segundo gráfico nos muestra la tendencia de la serie original, el siguiente nos muestra la estacionalidad de la serie original, y finalmente el último nos muestra los residuales, que es el resultados de quitar la tendencia y la parte estacional a la serie original.

dec = seasonal_decompose(X_train, model="additive")

series = [dec.observed,dec.trend, dec.seasonal, dec.resid]
titulos= ['Serie observada', 'Componente de tendencia', 'Componente estacional', 'Residuo']
fig, axes = plt.subplots(2,2, sharex=False, figsize=(16, 8))
k=0
for i in range(2):
  for j in range(2):
    serie = series[k]
    axes[i, j].plot(serie)
    axes[i, j].set_title(titulos[k])
    axes[i, j].set_xlabel('date')
    axes[i, j].set_ylabel('Valor')
    k +=1
plt.tight_layout()
plt.show()

Verificación de estacionariedad

Procedimiento de la prueba de Dickey-Fuller Aumentada (ADF)

La prueba ADF ajusta un modelo de regresión a la serie temporal de la siguiente forma: \[ \Delta y_t = \alpha + \beta t + \gamma y_{t-1} + \sum_{i=1}^{p} \phi_i \Delta y_{t-i} + \epsilon_t \] donde:

• \(\Delta y_t = y_t - y_{t-1}\) es la diferencia primera.
• \(\alpha\) es la constante.
• \(\beta t\) es un término de tendencia si se incluye (opcional).
• \(\gamma\) es el coeficiente de la raíz unitaria.
• \(\phi_i\) son los coeficientes de los rezagos.
• \(\epsilon_t\) es el término de error.

Posteriormente, se realiza una prueba de significancia sobre el parámetro \(\gamma\) (coeficiente de \(y_{t-1}\)) para determinar si existe o no una raíz unitaria.

Si el valor de \(\gamma\) es significativamente distinto de cero, se rechaza la hipótesis nula y se concluye que la serie es estacionaria. Si \(\gamma\) es cercano a cero o no es significativo, la serie se considera no estacionaria.

La hipótesis nula de la prueba de Dickey-Fuller aumentada es:

\(H_0 : \gamma = 0\) (significa que la serie no es estacionaria)

\(H_1 : \gamma < 0\) (significa que la serie es estacionaria)

# ADF
def adftest(data):
  result = adfuller(data, autolag='AIC')
  print(f'ADF Statistic: {result[0]}')
  print(f'p-value: {result[1]}')
  for key, value in result[4].items():
    print('Critial Values:')
    print(f'   {key}, {value}')
  # result
  for key, value in result[4].items():
    print('\t%s: %.3f' % (key, value))
  if result[0] < result[4]["5%"]:
    print ("Se rechaza Ho - La serie es estacionaria")
  else:
    print ("No se rechaza Ho - La serie no es estacionaria")

adftest(X_train.sales)

## ADF Statistic: 3.481780217160264
## p-value: 1.0
## Critial Values:
##    1%, -3.5335601309235605
## Critial Values:
##    5%, -2.9064436883991434
## Critial Values:
##    10%, -2.590723948576676
##  1%: -3.534
##  5%: -2.906
##  10%: -2.591
## No se rechaza Ho - La serie no es estacionaria

La prueba Dickey-Fuller, indica que la serie no es estacionaria, entonces, probamos con una diferencia (para quitar la tendencia) y después volvemos a realizar el test, en caso de no funcionar probamos con \(2\) diferencias, y nuevamente aplicamos el test.

Identificación del orden del modelo

Además observe que es necesario aplicar una transformación a la serie debido a que la varianza no es homogénea, vamos a gráficar algunos casos con diferencias y transformación logarítmica.

# Serie Original
fig, axes = plt.subplots(4,3, sharex=False, figsize=(16, 12))
axes[0, 0].plot(X_train.sales);
axes[0, 0].set_title('Serie original')
plot_acf(X_train.sales, ax=axes[0, 1]);
plot_pacf(X_train.sales, ax=axes[0, 2]);

# Logaritmo natural de la serie
axes[1, 0].plot(np.log(X_train.sales));
axes[1, 0].set_title('Logaritmo natural de la serie')
plot_acf(np.log(X_train.sales), ax=axes[1, 1]);
plot_pacf(np.log(X_train.sales), ax=axes[1, 2]);

# 1ra diferenciación
axes[2, 0].plot(np.log(X_train.sales).diff().diff(4));
axes[2, 0].set_title('1ra diferenciación')
plot_acf(np.log(X_train.sales).diff().diff(4).dropna(), ax=axes[2, 1]);
plot_pacf(np.log(X_train.sales).diff().diff(4).dropna(), ax=axes[2, 2]);

# 2da diferenciación
axes[3, 0].plot(np.log(X_train).diff().diff().diff(4));
axes[3, 0].set_title('2da diferenciación')
plot_acf(np.log(X_train.sales).diff().diff().diff(4).dropna(), ax=axes[3, 1]);
plot_pacf(np.log(X_train.sales).diff().diff().diff(4).dropna(), ax=axes[3, 2]);
plt.tight_layout()
plt.show()

Los casos anteriores muestran que con una diferenciación en la parte regular y una diferenciación en la parte estacional, y transformación logarítmica la serie es casi estacionarias, es decir, que no existe eun decaimiento exponencial tanto en la ACF, como en la PACF. Aplicamos la transformación logarítmica a la serie original, hacemos la diferenciación, y eliminamos el valor NA que se genere al diferenciar la serie.

X_train_log_diff =  np.log(X_train).diff().diff(4).dropna()

# gráficamos serie con una diferencia
plt.figure(figsize=(10,4))
plt.plot(X_train_log_diff, label='Train')
plt.title('Serie Jhonson')
plt.xlabel('date')
plt.ylabel('Sales')
plt.legend()
plt.show()

Aplicamos nuevamente el test de Dickey-Fuller

adftest(X_train_log_diff.sales)

## ADF Statistic: -6.082457671940423
## p-value: 1.0839818085936324e-07
## Critial Values:
##    1%, -3.5443688564814813
## Critial Values:
##    5%, -2.9110731481481484
## Critial Values:
##    10%, -2.5931902777777776
##  1%: -3.544
##  5%: -2.911
##  10%: -2.593
## Se rechaza Ho - La serie es estacionaria

Obtenemos una serie estacionaria, por lo tanto podemos contienuar con ka estimación de los parámetros.

Estimación

SARIMAX(endog, exog=None, order=(0, 0, 0), seasonal_order=(0, 0, 0, 0), trend=None, enforce_stationarity=True, enforce_invertibility=True, …)

Parámetros:

• endog : El proceso de la serie temporal observada.
• exog : Conjunto de regresores exógenos.
• order : El orden \((p,d,q)\) del modelo para los componentes autorregresivos \((p)\), diferencias \((d)\) y el orden de media móvil \((q)\). \(d\) es siempre un entero, mientras que \(p\) y \(q\) pueden ser enteros o listas de enteros.
• seasonal_order : El orden \((P,D,Q,s)\) del componente estacional del modelo para los parámetros AR(P), las diferencias I(D), los parámetros de media móvil MA(Q) y la periodicidad \((s)\). Por defecto es \((0, 0, 0, 0)\). \(D\) y \(s\) son siempre enteros, mientras que \(P\) y \(Q\) pueden ser enteros o listas de enteros positivos.
• trend : Parámetro que controla la tendencia determinista. Puede especificarse como una cadena en la que c indica un término constante, t indica una tendencia lineal en el tiempo, y ct incluye ambos.
• enforce_stationarity : Requerir o no que los parámetros autorregresivos correspondan a un proceso estacionario.
• enforce_stationarity : Requerir o no que los parámetros de media móvil correspondan a un proceso invertible.

model = SARIMAX(np.log(X_train), order=(0,1,1), seasonal_order=(1,1, 0, 4)).fit()
model.summary()

Evaluación del modelo

Una vez estimado el modelo procedemos a validar que los residuos cumplan con lo siguiente: - Se distribuyan normalmente. - Media cero y varianza constante. - No esten correlacionados

model.plot_diagnostics(figsize=(16,8));
plt.show()

Predicción

# Predicciones conjunto de entrenamiento
pred = model.get_prediction(step=1)

conficence_int = pred.conf_int()
predicted = pred.prediction_results
print(conficence_int)

##             lower sales  upper sales
## date                                
## 1960-01-01 -2771.807658  2771.807658
## 1960-04-01 -1960.306488  1959.621508
## 1960-07-01 -1960.426033  1959.501962
## 1960-10-01 -1960.126517  1959.801479
## 1961-01-01 -1386.896067  1384.911615
## ...                 ...          ...
## 1976-07-01     2.004123     2.361525
## 1976-10-01     1.761665     2.119067
## 1977-01-01     1.966533     2.323934
## 1977-04-01     2.134814     2.492216
## 1977-07-01     2.097555     2.454957
## 
## [71 rows x 2 columns]

# Predicción del conjunto de prueba
n = len(X_test)
forec = model.get_forecast(n)

conficence_int = np.exp(forec.conf_int())
forecast = np.exp(forec.predicted_mean)
print(conficence_int)

##             lower sales  upper sales
## 1977-10-01     6.547517     9.360395
## 1978-01-01     8.719397    12.732220
## 1978-04-01     9.516761    14.178123
## 1978-07-01     8.762997    13.306882
## 1978-10-01     6.759609    11.673095
## 1979-01-01     9.047982    16.212690
## 1979-04-01     9.715160    18.023344
## 1979-07-01     8.884934    17.034057
## 1979-10-01     6.750355    15.144596
## 1980-01-01     8.916708    21.113987
## 1980-04-01     9.528964    23.739034
## 1980-07-01     8.648348    22.606771
## 1980-10-01     6.596243    20.033938

plt.figure(figsize=(16,8))
plt.plot(X_train, label='Train')
plt.plot(X_test, label = 'Test')
plt.plot(forecast, label = 'Predicted')
plt.fill_between(forecast.index, conficence_int.iloc[:, 0], conficence_int.iloc[:, 1], color='b', alpha=.2)
plt.title('Serie Jhonson')
plt.xlabel('Quatrimestre')
plt.ylabel('Ventas')
plt.legend()
plt.show()

Métricas de error

Para obtener la métricas de error, empleamos la siguiente función:

def metrics_eval(model):
  residuals = model.resid
  aic = model.aic
  bic = model.bic
  mse = np.mean(residuals**2)
  rmse = np.sqrt(mse)
  mae = np.mean(np.abs(residuals))
  result = {'Metric': ['AIC', 'BIC', 'MSE', 'RMSE', 'MAE'], 'Values': [aic, bic, mse, rmse, mae]}
  return pd.DataFrame(result)

metrics = metrics_eval(model)
metrics

##   Metric      Values
## 0    AIC -121.581843
## 1    BIC -115.012879
## 2    MSE    0.020905
## 3   RMSE    0.144585
## 4    MAE    0.098418

Lo ideal sería probar otro modelo, y comparar las métricas de error, para ver cual es mejor, y sobre ello realizar pronósticos futuros.

Pronósticos (Forecasting)

Ya identificamos un modelo aceptable, ahora utilizamos todo el conjunto de datos, y realizamos un pronóstico de valores futuros.

model = SARIMAX(np.log(yt), order=(0,1,1), seasonal_order=(1,1, 0, 4)).fit()
print(model.summary())

##                                      SARIMAX Results                                      
## ==========================================================================================
## Dep. Variable:                              sales   No. Observations:                   84
## Model:             SARIMAX(0, 1, 1)x(1, 1, [], 4)   Log Likelihood                  78.457
## Date:                          lun., 11 nov. 2024   AIC                           -150.914
## Time:                                    05:20:25   BIC                           -143.805
## Sample:                                01-01-1960   HQIC                          -148.066
##                                      - 10-01-1980                                         
## Covariance Type:                              opg                                         
## ==============================================================================
##                  coef    std err          z      P>|z|      [0.025      0.975]
## ------------------------------------------------------------------------------
## ma.L1         -0.6800      0.096     -7.115      0.000      -0.867      -0.493
## ar.S.L4       -0.3216      0.108     -2.967      0.003      -0.534      -0.109
## sigma2         0.0079      0.002      5.150      0.000       0.005       0.011
## ===================================================================================
## Ljung-Box (L1) (Q):                   0.05   Jarque-Bera (JB):                 1.59
## Prob(Q):                              0.83   Prob(JB):                         0.45
## Heteroskedasticity (H):               0.45   Skew:                             0.08
## Prob(H) (two-sided):                  0.05   Kurtosis:                         2.33
## ===================================================================================
## 
## Warnings:
## [1] Covariance matrix calculated using the outer product of gradients (complex-step).

Consideramos un horizonte de predicción de 10 trimestres \(h=10\)

h = 10
forec = model.get_forecast(h)

conficence_int = np.exp(forec.conf_int()) # intervalos de confianza
forecast = np.exp(forec.predicted_mean) # predicción
print(conficence_int)

##             lower sales  upper sales
## 1981-01-01    15.423928    21.859021
## 1981-04-01    13.931732    20.091135
## 1981-07-01    15.323082    22.468047
## 1981-10-01    10.756174    16.024926
## 1982-01-01    16.064791    27.272916
## 1982-04-01    14.375544    25.202541
## 1982-07-01    15.624687    28.238227
## 1982-10-01    10.981130    20.428339
## 1983-01-01    16.185943    35.046623
## 1983-04-01    14.383716    32.653437

plt.figure(figsize=(12,5))
plt.plot(yt, label='Serie')
plt.plot(forecast, label = 'Pronóstico')
plt.fill_between(forecast.index, conficence_int.iloc[:, 0], conficence_int.iloc[:, 1], color='b', alpha=.2)
plt.title('Serie Jhonson & Jhonson')
plt.xlabel('Trimestre')
plt.ylabel('Ventas')
plt.legend()
plt.show()

Auto ARIMA

Comparemos el resultado anterior con auto arima

import pmdarima as pm
model = pm.auto_arima(np.log(yt), start_p=1, start_q=1,
                      test='adf',       #adftest para encontrar el valor óptimo de 'd'
                      max_p=3, max_q=3, # máx p y q
                      m=4,              # frecuencia de la serie
                      d=1,              # diferenciación de la parte regular (no estacional)
                      seasonal=True,    #
                      start_P=0, start_Q=0,  # valores de inicip de P y Q de la parte estacional
                      D=1,                   # Diferenciación de la parte estacional
                      suppress_warnings=True,
                      stepwise=True)

print(model.summary())

##                                      SARIMAX Results                                      
## ==========================================================================================
## Dep. Variable:                                  y   No. Observations:                   84
## Model:             SARIMAX(0, 1, 1)x(1, 1, [], 4)   Log Likelihood                  78.457
## Date:                          lun., 11 nov. 2024   AIC                           -150.914
## Time:                                    05:20:33   BIC                           -143.805
## Sample:                                01-01-1960   HQIC                          -148.066
##                                      - 10-01-1980                                         
## Covariance Type:                              opg                                         
## ==============================================================================
##                  coef    std err          z      P>|z|      [0.025      0.975]
## ------------------------------------------------------------------------------
## ma.L1         -0.6800      0.096     -7.115      0.000      -0.867      -0.493
## ar.S.L4       -0.3216      0.108     -2.967      0.003      -0.534      -0.109
## sigma2         0.0079      0.002      5.150      0.000       0.005       0.011
## ===================================================================================
## Ljung-Box (L1) (Q):                   0.05   Jarque-Bera (JB):                 1.59
## Prob(Q):                              0.83   Prob(JB):                         0.45
## Heteroskedasticity (H):               0.45   Skew:                             0.08
## Prob(H) (two-sided):                  0.05   Kurtosis:                         2.33
## ===================================================================================
## 
## Warnings:
## [1] Covariance matrix calculated using the outer product of gradients (complex-step).

Observamos que el método \(auto\_arima()\) nos lleva a la misma solución, el modelo explicito sería:

\[\begin{eqnarray*} \Phi_1(B^{s=4}) \nabla^{d=1} \nabla^{D=1} x_t & = & \theta_1(B)u_t \end{eqnarray*}\]

\[\begin{eqnarray*} -0.3216(B^{s=4}) \nabla^{d=1} \nabla^{D=1} x_t & = & -0.68(B)u_t \end{eqnarray*}\]

donde \[\begin{eqnarray*} \phi_p(B) & = & (1- \phi_1 B - \phi_2 B^2 -... - \phi_p B^p)\\ \Phi_P(B^s) & = & (1- \Phi_1 B^{s} - \Phi_2 B^{2s} ... - \Phi_p B^{Ps})\\ \nabla^d & = & (1-B)^d\\ \nabla^D & = & (1-B^s)^D\\ \Theta_Q(B^s) & = & (1- \Theta_1 B^s - \Theta_2 B^{2s} - ... - \Theta_Q B^{Qs})\\ \theta_q(B) & = & (1- \theta_1 B - \theta_2 B^{2} - ... - \theta_q B^{Q})\\ \end{eqnarray*}\]