Venta licores.Predicción modelo ARIMA
Tabla de contenidos
Introducción
- Orden cronológico y equidistancia temporal
Estacionalidad
- Enfoque visual
- Media y desviación estándar móvil
- Enfoque estadístico
- Dicker-Fuller Ampliado (Test ADF)
- Transformación
- Transformación logarítmica
- Extracción de la tendencia
- Tendencia evolutiva
- Diferenciación (integrada de primer orden)
- Descomposición de la serie temporal
- Enfoque visual
Análisis exploratorio temporal
- Función de autocorrelación simple (ACF)
- Función de autocorrelación parcial (PACF)
Modelo ARIMA
- Auto Arima: Identificar parámetros óptimos para modelo ARIMA
- Construcción modelo ARIMA
- Validación
- Transformación inversa & Predicción
Introducción¶
Las instancias que recogen el registro de compras de licores del proveedor están compuestas por la variable temporal de forma secuencial. El registro se realiza a diario, aunque los datos no cumplen una equidistancia temporal. Por ello para el estudio fijamos cambios mensuales.
Antes de comenzar, dado que el tamaño de los ficheros es relativamente grande para trabajar en un entorno local se adjunta Notebook con la implementación multiproceo de lectura y preparación de los datos.
A modo sintétizado, creamos varios procesos que ejecutan el código de nuestro programa en Python de manera paralela utlizando el módulo multiprocessing. El módulo multiprocessing implementa la clase Process que representa un flujo de ejecución que corre en un proceso individual. Para crear diferentes procesos controlados por muestro programa principal, crearemos diferentes instancias de la clase Process, especificando la función a ejecutar y el parámetro.
La publicación completa también se encuentra en el repositorio
Orden cronológico y equidistancia temporal¶
Como hemos hecho mención antes, los datos originales son diarios. Se remuestrean mensualmente por el promedio.
Validamos los intervalos de tiempo , y que se cumpla la equidistancia temporal por mes.
# Librerias
import pandas as pd
import numpy as np
import seaborn as sns
import matplotlib.pyplot as plt
plt.style.use('Solarize_Light2')
%matplotlib inline
# Lectura datos e indexacion por fecha
df = pd.read_csv('monthly_data_arima.csv' ,parse_dates=['Date'], index_col=['Date'])
# Validacion intervalos tiempo mensuales
print("Se cumple la equidistancia temporal: {} \n\n".format(
np.unique(df.index + pd.offsets.MonthEnd(0) == df.index)[0]))
df.head()
Se cumple la equidistancia temporal: True
| state_bottle_cost | state_bottle_retail | sale_bottles | sale_dollars | sale_liters | |
|---|---|---|---|---|---|
| Date | |||||
| 2012-01-31 | 10.881787 | 16.341626 | 8.593687 | 134.666037 | 7.898504 |
| 2012-02-29 | 11.207359 | 16.808497 | 8.701056 | 143.331408 | 8.164135 |
| 2012-03-31 | 11.024545 | 16.542414 | 8.695276 | 138.437188 | 7.919058 |
| 2012-04-30 | 11.127428 | 16.717924 | 9.100950 | 149.534302 | 8.590926 |
| 2012-05-31 | 11.222185 | 16.829469 | 8.888101 | 144.791238 | 8.123422 |
Estacionalidad¶
Comprobamos la periodicidad y si se observa en la serie un patrón que se repite.
Series temporales como ARIMA, asumen que los datos subyacentes presentan estacionalidad. \ Se deberían cumplir tres criterios básicos para poder considerar una serie de tiempo como estacionaria:
- La media de la serie no debe ser una función de tiempo
- La varianza de la serie no debe ser una función del tiempo
- La covarianza de la serie no sebe ser una función del tiempo
Enfoque visual¶
El método de medias móviles puede sernos útil en esta ocasión si queremos observar la tendencia de la serie temporal sin tener que ajustarnos a un función previa. Nos ofrece una visión suavizada de la serie, ya que promediando varios valores se elimina parte de los movimientos irregulares de la serie.
Media y desviación estándar móvil¶
La media móvil simple es el promedio no ponderado de los M puntos previos de la serie. La selección de M (sliding windows o ventana) depende del nivel de suavizado que se desee. Para una secuencia de valores, calculamos la media móvil simple en un período de tiempo t de la siguiente forma:
$$MMS_{t}= \frac{ x_{1} + x_{t-1} +x_{t-2} + x_{M-{t-1)}} }{M}$$La desviación estándar móvil indica la volatilidad de los M puntos previos de la serie.
Se obtiene a partir de un periodo t de la media móvil simple calculando la desviación de la variable respecto a las medias móviles
$$Z_{k} = \sqrt { \frac { \sum_{i=k-m+1} ^{k} (y_{i} - \overline{y}) } {m-1}} \forall k =1,...,n$$donde $y_{i}$ valor de la variable de la i-th instancia, $\overline{y}$ promedio de la variable, m la ventana y n el número de instancias.
# Descartamos variables
df = df.drop(['state_bottle_cost', 'state_bottle_retail', 'sale_bottles','sale_liters'], axis=1)
df['rolling_mean'] = df.sale_dollars.rolling(12).mean()
df['rolling_std'] = df.sale_dollars.rolling(12).std()
plt.figure(figsize=(18,13))
sns.lineplot(x=df.index, y=df['sale_dollars'], label='Original')
sns.lineplot(x=df.index, y=df['rolling_mean'], color='black', label='rolling mean')
sns.lineplot(x=df.index, y=df['rolling_std'], color='orange', label='rolling std')
plt.ylabel('Sale dollars', fontsize=14)
plt.xlabel('Año', fontsize=14)
plt.title('Media y desviación estándar móvil', fontsize=20)
plt.tight_layout()
plt.show()
Enfoque estadístico¶
Test ADCF¶
La prueba de ADF es una prueba de raíz unitaria. Cuando una serie tiene raíz unitaria no es estacionaria.
Dada la ecuación de la serie temporal, si el valor de la raíz de la ecuación es igual a $\alpha$ = 1, el proceso no es estacionario : $$Y_{t} = \alpha Y_{t-1} + \beta X_{e} + \epsilon$$
donde $Y_{t}$ es el valor de la serie temporal en t y $X_{e}$ es una variable exógena.
El test de Dickey-Fuller no es solamente una prueba sino que requiere una estrategía de análisis para su aplicación. Consiste en determinar cuál de los casos determina una mejor aproximación a la serie original $Y_{t}$. La decisión se toma con base en el menor AIC.
Para estimar el efecto estacional, en la implementación del test ADCF se desestacionaliza la serie original mediante una diferenciación estacional.El argumento autolag de la función adfuller selecciona aquella diferenciación que minimice la varianza mediante la medida AIC.
La hipótesis nula en este test es que la serie no es estacionaria, luego si el valor resultante (p-value) es menor de 0.05 (ya que el p-value representa la probabilidad de la hipótesis tomada) indica que la serie es estacionaria con un nivel de confianza del 95%. En caso contrario, no habría evidencia para rechazar la hipótesis de no estacionalidad.
- Hipótesis nula (H0) : Serie temporal tiene raíz unitaria. (Serie temporal no es estacionaria)
- Hipótesis alternativa (H1): Serie temporal no tiene raíz unitaria. (Serie temporal estacionaria)
Hay dos formas para rechazar la hipótesis nula:
Por un lado, H0 puede rechazarse como hemos hecho mención anterior por el valor resultante
- p.value > nivel confianza ( por defecto 0.05) No se rechaza H0, serie temporal tiene raíz unitaria y no es estacionaria
p.value <= nivel confianza ( por defecto 0.05) Se rechaza H0, serie temporal no tiene raíz unitaria y es estacionaria
Por otro lado, H0 puede rechazarse si el test estadístico es menor al valor crítico
- estadístico ADF > valor crítico No se rechaza H0, serie temporal tiene raíz unitaria y no es estacionaria
- estadístico ADF < valor crítico Se rechaza H0, serie temporal no tiene raíz unitaria y es estacionaria
from statsmodels.tsa.stattools import adfuller
def visualize_adfuller_results(df, series, title):
"""
Test de Dicker-Fuller Aumentada y visualización según p-value y valores críticos
Args:
df: dataframe con las variable target y temporal
series: variable de la serie temporal
title: titulo del gráfico
Returns:
dfoutput: pandas con los resultados estadísticos
plot: gráfico a partir de los resultados del test
"""
# test Dickey-FUller Aumentada
dftest = adfuller(df[series], autolag='AIC')
# Dataframe con los resultados del test.Se establece index y tras la iteracion sobre los resultados
# del test se asginan los valores obtenidos
dfoutput = pd.Series(dftest[0:4], index=['Test Statistics', 'p-value', 'No. of lags used',
'Number of observations used'])
for key, value in dftest[4].items():
dfoutput['Critical value (%s)' %key] = value
significance_level = 0.05
# Generamos grafico con los resultados anteriores
plt.figure(figsize=(18,13))
significance_level = 0.05
if (dfoutput[1] < significance_level) & (dfoutput[0] < dfoutput[4]):
linecolor='forestgreen'
elif (dfoutput[1] < significance_level) & ((dfoutput[0] < dfoutput[5])):
linecolor='orange'
elif (dfoutput[1] < significance_level) & ((dfoutput[0] < dfoutput[6])):
linecolor='red'
else:
linecolor = 'purple'
sns.lineplot(x=df.index, y=df[series], color=linecolor)
sns.lineplot(x=df.index, y=df['rolling_mean'], color='black', label='rolling mean')
sns.lineplot(x=df.index, y=df['rolling_std'], color='orange', label='rolling std')
plt.ylabel(title, fontsize=14)
plt.xlabel('Año', fontsize=14)
plt.title(f'Test ADF {dfoutput[0]:0.3f},\
p-value: {dfoutput[1]:0.3f}\
Critical value (1%) {dfoutput[4]:0.3f}\
Critical value (5%) {dfoutput[5]:0.3f}\
Critical value (10%) {dfoutput[6]:0.3f}', fontsize=16)
plt.tight_layout()
plt.show()
print("Resultado test Dicker-Fuller Aumentada")
return dfoutput
# Serie original
visualize_adfuller_results(df,'sale_dollars', 'Sale dollars')
Resultado test Dicker-Fuller Aumentada
Test Statistics 1.665810 p-value 0.998042 No. of lags used 12.000000 Number of observations used 84.000000 Critical value (1%) -3.510712 Critical value (5%) -2.896616 Critical value (10%) -2.585482 dtype: float64
Según los resultados del test ADCF, el p-value es muy amplio (máximo 1). Al mismo tiempo los valores críticos no se encuentran cercanos al valor del test estadístico. Por lo tanto, podemos afirmar que nuestra serie temporal en este punto no es estacionaria.
Si los datos no presentan estacionalidad, pero se quiere usar el modelo ARIMA (que requiere esta característica), los datos han de ser transformados.
Transformación¶
La serie no es constante en varianza (serie heterocedástica), y aumenta ligeramente con el nivel de la serie, por lo tanto no sigue una distribución normal. Es recomendable realizar una transformación para acercarse a una distribución normal.
Existen diversas maneras para conseguir la estacionalidad a través de la transformación de los datos. Por ejemplo, tomar el $log_{10}$, $log_{e}$, raíz cuadrada o cúbica, desplazamiento temporal..
Transformación logarítmica¶
En nuestro caso comenzamos con una transformación logarítmica. Nuestro objetivo es eliminar el componente de tendencia. En la transformación logarítmica los valores pequeños se expandirán y los grandes se contrae, y dado que los valores son positivos la podemos aplicar.
# Estimacion tendencia
# Transformacion logaritmica para convertir la serie en estacionaria
logScale_df = pd.DataFrame(np.log(df.sale_dollars))
logScale_df['rolling_mean'] = logScale_df.sale_dollars.rolling(12).mean()
logScale_df['rolling_std'] = logScale_df.sale_dollars.rolling(12).std()
# Grafico comparacion serie tras la transormacion
f, ax = plt.subplots(nrows=1, ncols=2, figsize=(15,9))
# Serie temporal original
sns.lineplot(x=df.index, y=df['sale_dollars'], label='Original', ax=ax[0])
sns.lineplot(x=df.index, y=df['rolling_mean'], color='black', label='rolling mean', ax=ax[0])
ax[0].set_ylabel('Sale dollars', fontsize=14)
ax[0].set_xlabel('Año', fontsize=14)
ax[0].set_title('ST original', fontsize=20)
# Serie temporal tras la transformación logarítmica
sns.lineplot(x=logScale_df.index, y=logScale_df['sale_dollars'], label='Original', ax=ax[1])
sns.lineplot(x=logScale_df.index, y=logScale_df['rolling_mean'], color='black', label='rolling mean', ax=ax[1])
ax[1].set_ylabel('Sale dollars', fontsize=14)
ax[1].set_xlabel('Año', fontsize=14)
ax[1].set_title('ST transformación logarítmica', fontsize=20)
plt.tight_layout()
plt.show()
A partir del gráfico anterior, vemos que aunque la media móvil no es estacionaria, es ligeramente mejor que el caso original donde no se aplica ninguna transformación.
Extracción de la tendencia¶
Tendencia evolutiva
Sabemos que la serie temporal tanto con la escala logarítmica como con la media móvil tiene componente de tendencia.De esta forma, se deduce que substrayendo uno del otro debería eliminar la componente de tendencia.
$$log(X) = E_{Xt} + T_{Xt}$$$$Media \space móvil \space de\space log(A) = E_{At} + T_{At}$$$$Serie \space final \space R = X - A = (E_{Xt} + T_{Xt}) - (E_{At} + T_{At}) = (E_{Xt} -E_{At}) + ( T_{Xt}- T_{At})$$donde $T_{t}$ es la tendencia, $E_{t}$ es la componente estacional, que constituyen la señal o parte determinística.
Dado que X y A son series y sus medias móviles, su tendencia será más o menso similar, por lo que $( T_{Xt}- T_{At}) \simeq 0$. De este modo la componente de tendencia quedará prácticamente eliminada, y tendremos
$$R = (E_{Xt} -E_{At})$$nuestra curva final sin tendencia
# ELiminar componente de tendencia
noTrend_df= pd.DataFrame({'sale_dollars': logScale_df.sale_dollars.values - \
logScale_df.rolling_mean.values},
index=df.index)
# Eliminar valores NAN
noTrend_df.dropna(inplace=True)
noTrend_df['rolling_mean'] = noTrend_df.sale_dollars.rolling(12).mean()
noTrend_df['rolling_std'] = noTrend_df.sale_dollars.rolling(12).std()
# Test adfuller aumentada
visualize_adfuller_results(noTrend_df,'sale_dollars', 'Sale dollars')
Resultado test Dicker-Fuller Aumentada
Test Statistics -1.325895 p-value 0.617244 No. of lags used 12.000000 Number of observations used 73.000000 Critical value (1%) -3.523284 Critical value (5%) -2.902031 Critical value (10%) -2.588371 dtype: float64
Diferenciación (integrada de primer orden)¶
Otro método para eliminar la tendencia consiste en suponer que la tendencia evoluciona lentamente en el tiempo, de manera que en el instante t la tendencia debe estar próxima a la tendencia en el instante t-1. De esta manera, si restamos a cada valor de la serie el valor anterior, la serie resultante estará aproximadamnete libre de tendencia.
$$Y_{t} = X_{t} - X_{t-1}$$donde $X_{t}$ es la serie original
# ELiminar componente de tendencia
first_diff= pd.DataFrame({'sale_dollars': logScale_df.sale_dollars.values - \
logScale_df.sale_dollars.shift()},
index=logScale_df.index)
# Eliminar valores NAN
first_diff.dropna(inplace=True)
first_diff['rolling_mean'] = first_diff.sale_dollars.rolling(12).mean()
first_diff['rolling_std'] = first_diff.sale_dollars.rolling(12).std()
# Test adfuller aumentada
visualize_adfuller_results(first_diff, 'sale_dollars', 'Sale dollars')
Resultado test Dicker-Fuller Aumentada
Test Statistics -7.306384e+00 p-value 1.298228e-10 No. of lags used 1.200000e+01 Number of observations used 8.300000e+01 Critical value (1%) -3.511712e+00 Critical value (5%) -2.897048e+00 Critical value (10%) -2.585713e+00 dtype: float64
A partir de los dos últimos gráficos, podemos obervar que, visualmente a partir de los valores de la media y desviación estándar móvil la serie presenta estacionalidad y ambos no son funciones de tiempo.
La diferenciación presenta mejor estado que la en la tendencia evolutiva en cuanto a la desviación estándar ya que el incremento es menor en el primer caso.
Sobre los resultados del test Dicker-Fuller Aumentada:
- El p-valor es cercano al 0 en la implementación de la diferenciación de la transformación logarítmica
- Los valores críticos no se acercan al test estadístico en la diferencición, sin embargo en la tendencia evolutiva se consiguen mejores resultados
Descomposición de la serie temporal¶
Mediante la implementación de la función seasonal_decompose en la librería statsmodels obtenemos las 3 componentes de la serie. En este punto la serie es estacional.La tendencia revela que las ventas presentan crecimientos y decrecimientos abruptos y los registros de los meses finales van decreciendo. Las ventas tienen una periodicidad de 12 meses, lo que puede explicar la variabilidad. Nos centraremos en el carácter de la componente residual.
Como presuposición personal, no podemos hacer una lectura del todo estadística o numércia sin tener en cuenta que se tratan de datos de los que se deconoce la frecuencia de registro de ventas, o prácticas relacionados a su recopilación.
import statsmodels.api as sm
from pylab import rcParams
rcParams['figure.figsize'] = 18,9
decomposition = sm.tsa.seasonal_decompose(first_diff['sale_dollars'], model='aditive')
fig = decomposition.plot()
Análisis exploratorio¶
Para poder efectuar inferencias sobre los parámetros de un proceso estocástico, es preciso imponer restricciones a este último. Las restricciones que se imponen habitualmente son que sea estacionario y ergódico. Ya se ha trabajado sobre la estacionalidad en los apartados anteriores, por lo que explicaremos un proceso ergódico a través de la autocorrelación.
Las funciones de autocorrelación simple y parcial constituyen uno de los instrumentos clave para ajustar el model que genera nuestra serie temporal.
Función de autocorrelación simple (ACF)¶
Mide la correlación entre dos variables separadas por k periodos. Proporciona la estructura de dependencia lineal de la misma. Así la ACF va a ser una sucesión e valores ($p_{1},p_{2}...p_{k}$) que representan cómo influye una observación sobre las siguientes.
Función de autocorrelación parcial (PACF)¶
Mide la correlación entre dos variables separadas por k periodos pero cuando no se considera la dependencia creada por los retardos intermedios existentes entre ambas
from statsmodels.graphics.tsaplots import plot_acf, plot_pacf
f, ax = plt.subplots(nrows=2, ncols=1, figsize=(16, 8))
plot_acf(first_diff['sale_dollars'], ax=ax[0])
plot_pacf(first_diff['sale_dollars'], ax=ax[1])
plt.show()
Sobre los resultado de ambas funciones:
Para asegurarnos que no hemos aplicado un diferenciación extra de
la necesitada por la serie temporal para ser estacionaria,
nos fijamos en el primer rezago de la gráfica de autocorrelación
y es diferente al valor -0.5, por lo que no consideramos que esté
sobrediferenciada.
En ambas funciones observamos que a medida que se avanza sobre
los periodos la autocorrelación no se hace más pequeña, sino de lo
contrario el comportamiento es irregular, por lo que no se trata de un
proceso ergódico. Esto muestra cierta dependencia entre los datos en
función a los valores anteriores
Se alternan relaciones positivas y negativas entre las observaciones.
Modelar la serie con ARIMA nos va a permitir tener en consideración esta
dependencia que existe entre los datos.Modelo ARIMA¶
Los modelos ARIMA (AutoRegresive Integrated Moving Average) están diseñados para modelizar datos de carácter estacionario. Se derivan de las siguientes componentes:
- Procesos Autoregresivos AR(p) :
el valor actual de la serie puede explicarse o
predecirse en función de p valores pasados, más
un término de error
- Proceso de Medias Móviles MA(q):
el valor de la variable en un período puede explicarse
en función de un término independiente y una sucesión
de errores correspondientes a períodos precedentes
- Proceso Integrado I(d) :
la serie se diferencia d veces para hacerla estacionaria
cuando presenta cambios de nivel en el tiempo o
la varianza no es constanteAuto Arima: Identificar parámetros óptimos para modelo ARIMA¶
Previo a la construcción del modelo, usarmos la función auto_arima que evalúa entre todos los posibles modelos, para identificar los parámetros más óptimos.Conduce diferentes test, en nuestro caso seguimos usando ADF (Augmented Dickey-Fuller por sus siglas en inglés).
En los parámetros start_p y start_q indicamos el orden del proceso autoregresivo y de medias móviles. El período de nuestra serie diferenciada la indicamos en el parámetro m, en nuestro caso es mensual (m=12).
Aunque hemos diferenciado la serie logarítmica manualmente, pasamos la serie logarítmica pues en la transformación inversa nos será más manejable e indicamos d=1 a fin de que sea diferenciada en primer orden.
El ajuste es estacional (seasonal) y por último, se usa un algoritmo gradual (stepwise) que permite la parametrización más óptima al mismo tiempo que reduce el sobre-ajuste del modelo.
from pmdarima.arima import auto_arima
model_autoArima = auto_arima(first_diff.sale_dollars, start_p=0, start_q=0,
test='adf',
m=12,
seasonal=True,
trace=True,
suppress_warnings=True,
stepwise=True)
print(model_autoArima.summary())
Performing stepwise search to minimize aic
ARIMA(0,0,0)(1,0,1)[12] intercept : AIC=-71.211, Time=0.43 sec
ARIMA(0,0,0)(0,0,0)[12] intercept : AIC=-55.852, Time=0.04 sec
ARIMA(1,0,0)(1,0,0)[12] intercept : AIC=-107.416, Time=0.24 sec
ARIMA(0,0,1)(0,0,1)[12] intercept : AIC=inf, Time=0.60 sec
ARIMA(0,0,0)(0,0,0)[12] : AIC=-57.840, Time=0.02 sec
ARIMA(1,0,0)(0,0,0)[12] intercept : AIC=-93.409, Time=0.05 sec
ARIMA(1,0,0)(2,0,0)[12] intercept : AIC=inf, Time=0.73 sec
ARIMA(1,0,0)(1,0,1)[12] intercept : AIC=-126.633, Time=0.46 sec
ARIMA(1,0,0)(0,0,1)[12] intercept : AIC=-98.923, Time=0.16 sec
ARIMA(1,0,0)(2,0,1)[12] intercept : AIC=-137.749, Time=1.64 sec
ARIMA(1,0,0)(2,0,2)[12] intercept : AIC=-136.279, Time=5.58 sec
ARIMA(1,0,0)(1,0,2)[12] intercept : AIC=inf, Time=1.55 sec
ARIMA(0,0,0)(2,0,1)[12] intercept : AIC=-76.804, Time=0.66 sec
ARIMA(2,0,0)(2,0,1)[12] intercept : AIC=-155.950, Time=1.77 sec
ARIMA(2,0,0)(1,0,1)[12] intercept : AIC=-149.464, Time=0.83 sec
ARIMA(2,0,0)(2,0,0)[12] intercept : AIC=inf, Time=1.07 sec
ARIMA(2,0,0)(2,0,2)[12] intercept : AIC=-153.950, Time=2.15 sec
ARIMA(2,0,0)(1,0,0)[12] intercept : AIC=-131.616, Time=0.42 sec
ARIMA(2,0,0)(1,0,2)[12] intercept : AIC=-152.018, Time=1.78 sec
ARIMA(3,0,0)(2,0,1)[12] intercept : AIC=-158.566, Time=4.87 sec
ARIMA(3,0,0)(1,0,1)[12] intercept : AIC=-153.269, Time=2.94 sec
ARIMA(3,0,0)(2,0,0)[12] intercept : AIC=inf, Time=2.40 sec
ARIMA(3,0,0)(2,0,2)[12] intercept : AIC=-156.236, Time=2.39 sec
ARIMA(3,0,0)(1,0,0)[12] intercept : AIC=-134.783, Time=0.67 sec
ARIMA(3,0,0)(1,0,2)[12] intercept : AIC=-155.171, Time=1.87 sec
ARIMA(4,0,0)(2,0,1)[12] intercept : AIC=-160.888, Time=2.30 sec
ARIMA(4,0,0)(1,0,1)[12] intercept : AIC=-157.728, Time=1.07 sec
ARIMA(4,0,0)(2,0,0)[12] intercept : AIC=inf, Time=2.68 sec
ARIMA(4,0,0)(2,0,2)[12] intercept : AIC=-158.566, Time=3.73 sec
ARIMA(4,0,0)(1,0,0)[12] intercept : AIC=-141.833, Time=1.05 sec
ARIMA(4,0,0)(1,0,2)[12] intercept : AIC=-157.992, Time=3.01 sec
ARIMA(5,0,0)(2,0,1)[12] intercept : AIC=-163.639, Time=5.45 sec
ARIMA(5,0,0)(1,0,1)[12] intercept : AIC=-161.874, Time=1.39 sec
ARIMA(5,0,0)(2,0,0)[12] intercept : AIC=inf, Time=2.51 sec
ARIMA(5,0,0)(2,0,2)[12] intercept : AIC=-161.198, Time=4.26 sec
ARIMA(5,0,0)(1,0,0)[12] intercept : AIC=-150.408, Time=1.61 sec
ARIMA(5,0,0)(1,0,2)[12] intercept : AIC=-161.856, Time=3.39 sec
ARIMA(5,0,1)(2,0,1)[12] intercept : AIC=-163.695, Time=3.51 sec
ARIMA(5,0,1)(1,0,1)[12] intercept : AIC=-160.292, Time=1.55 sec
ARIMA(5,0,1)(2,0,0)[12] intercept : AIC=-165.719, Time=2.94 sec
ARIMA(5,0,1)(1,0,0)[12] intercept : AIC=-149.507, Time=1.32 sec
ARIMA(4,0,1)(2,0,0)[12] intercept : AIC=-168.075, Time=2.88 sec
ARIMA(4,0,1)(1,0,0)[12] intercept : AIC=-151.807, Time=2.89 sec
ARIMA(4,0,1)(2,0,1)[12] intercept : AIC=-166.372, Time=5.13 sec
ARIMA(4,0,1)(1,0,1)[12] intercept : AIC=-161.327, Time=1.26 sec
ARIMA(3,0,1)(2,0,0)[12] intercept : AIC=-169.641, Time=2.22 sec
ARIMA(3,0,1)(1,0,0)[12] intercept : AIC=-153.939, Time=1.15 sec
ARIMA(3,0,1)(2,0,1)[12] intercept : AIC=-167.279, Time=2.32 sec
ARIMA(3,0,1)(1,0,1)[12] intercept : AIC=-165.975, Time=1.15 sec
ARIMA(2,0,1)(2,0,0)[12] intercept : AIC=-171.770, Time=2.38 sec
ARIMA(2,0,1)(1,0,0)[12] intercept : AIC=-155.701, Time=1.89 sec
ARIMA(2,0,1)(2,0,1)[12] intercept : AIC=-169.919, Time=4.15 sec
ARIMA(2,0,1)(1,0,1)[12] intercept : AIC=-167.562, Time=1.18 sec
ARIMA(1,0,1)(2,0,0)[12] intercept : AIC=-173.436, Time=1.80 sec
ARIMA(1,0,1)(1,0,0)[12] intercept : AIC=inf, Time=0.72 sec
ARIMA(1,0,1)(2,0,1)[12] intercept : AIC=-170.462, Time=1.50 sec
ARIMA(1,0,1)(1,0,1)[12] intercept : AIC=-168.949, Time=0.91 sec
ARIMA(0,0,1)(2,0,0)[12] intercept : AIC=-161.106, Time=1.25 sec
ARIMA(1,0,2)(2,0,0)[12] intercept : AIC=-171.779, Time=1.65 sec
ARIMA(0,0,0)(2,0,0)[12] intercept : AIC=-77.935, Time=0.47 sec
ARIMA(0,0,2)(2,0,0)[12] intercept : AIC=-171.626, Time=1.29 sec
ARIMA(2,0,2)(2,0,0)[12] intercept : AIC=-169.909, Time=2.25 sec
ARIMA(1,0,1)(2,0,0)[12] : AIC=-174.693, Time=1.74 sec
ARIMA(1,0,1)(1,0,0)[12] : AIC=-154.666, Time=0.55 sec
ARIMA(1,0,1)(2,0,1)[12] : AIC=-172.868, Time=3.78 sec
ARIMA(1,0,1)(1,0,1)[12] : AIC=-170.665, Time=0.69 sec
ARIMA(0,0,1)(2,0,0)[12] : AIC=-161.027, Time=0.82 sec
ARIMA(1,0,0)(2,0,0)[12] : AIC=inf, Time=0.51 sec
ARIMA(2,0,1)(2,0,0)[12] : AIC=-173.145, Time=1.79 sec
ARIMA(1,0,2)(2,0,0)[12] : AIC=-173.038, Time=1.00 sec
ARIMA(0,0,0)(2,0,0)[12] : AIC=-79.920, Time=0.18 sec
ARIMA(0,0,2)(2,0,0)[12] : AIC=-173.388, Time=0.69 sec
ARIMA(2,0,0)(2,0,0)[12] : AIC=inf, Time=0.60 sec
ARIMA(2,0,2)(2,0,0)[12] : AIC=-171.238, Time=1.71 sec
Best model: ARIMA(1,0,1)(2,0,0)[12]
Total fit time: 131.708 seconds
SARIMAX Results
===========================================================================================
Dep. Variable: y No. Observations: 96
Model: SARIMAX(1, 0, 1)x(2, 0, [], 12) Log Likelihood 92.347
Date: Thu, 18 Feb 2021 AIC -174.693
Time: 21:54:14 BIC -161.872
Sample: 0 HQIC -169.511
- 96
Covariance Type: opg
==============================================================================
coef std err z P>|z| [0.025 0.975]
------------------------------------------------------------------------------
ar.L1 -0.4341 0.100 -4.337 0.000 -0.630 -0.238
ma.L1 -0.8534 0.071 -12.012 0.000 -0.993 -0.714
ar.S.L12 0.2668 0.075 3.539 0.000 0.119 0.415
ar.S.L24 0.5919 0.093 6.362 0.000 0.410 0.774
sigma2 0.0070 0.001 7.399 0.000 0.005 0.009
===================================================================================
Ljung-Box (L1) (Q): 0.11 Jarque-Bera (JB): 7.34
Prob(Q): 0.74 Prob(JB): 0.03
Heteroskedasticity (H): 0.97 Skew: 0.57
Prob(H) (two-sided): 0.94 Kurtosis: 3.72
===================================================================================
Warnings:
[1] Covariance matrix calculated using the outer product of gradients (complex-step).
Construcción modelo ARIMA¶
Generamos el modelo ARIMA con órdenes p,d,q con valores 1,0 y 1 respectivamente.También indicamos las órdenes de estacionalidad.
import warnings
warnings.filterwarnings('ignore')
from statsmodels.tsa.statespace.sarimax import SARIMAX
# Modelo Arima con los parametros obtenidos como optimos
# Parametros del modelo y parametros del orden la componente estacional
model = SARIMAX(first_diff.sale_dollars, order=(1,0,1), seasonal_order=(2,0,0,12),
enforce_stationary=False,
enforce_invertibility=False)
fitted = model.fit(disp=-1)
Validación¶
A fin de validar el modelo, generamos el pronóstico desde mediados del último año (2019-Junio hasta 2020-Enero) y lo comparamos con los valores reales de la serie en el mismo perído.
pred = fitted.get_prediction(start=pd.to_datetime('2019-06-30'), dynamic=False)
pred_ci = pred.conf_int()
ax = first_diff.sale_dollars.loc['2018':].plot(label='Sale Dollars')
pred.predicted_mean.plot(ax=ax, label='Pronóstico un paso adelante ', alpha=.7,
figsize=(14,7), color='red')
ax.fill_between(pred_ci.index,
pred_ci.iloc[:,0],
pred_ci.iloc[:,1], color='k', alpha=.10)
plt.ylabel('Sale dollars', fontsize=14)
plt.xlabel('Año', fontsize=14)
plt.title('Validación pronóstico del modelo', fontsize=20)
plt.legend()
plt.tight_layout()
plt.show()
y_forecasted = pred.predicted_mean
y_truth = first_diff.sale_dollars.loc['2019-01':]
mse = ((y_forecasted - y_truth)**2).mean()
print('Error cuadrático medio del pronóstico es {}'.format(round(mse,2)))
print('Raíz del error cuadrático medio del pronóstico {}'.format(round(np.sqrt(mse), 2)))
Error cuadrático medio del pronóstico es 0.01 Raíz del error cuadrático medio del pronóstico 0.09
La RDCM es considerablemente baja, lo que nos sugiere que el ajuste del modelo es bueno y puede sernos útil en la estimación
Transformación inversa & Predicción¶
Visualizamos los resultados del modelo tras realizar la transformación inversa de la serie logarítmica diferenciada. Las predicciones del modelo presentan underfitting (ajuste por debajo).Los motivos pueden ser desde la propia transformación incial de la serie, en la que se puede optar por métodos más adecuados por el tipo de serie hasta por los datos de entrenamiento, que son muy pocos y nuestro modelo no será capaz de generalizar el conocimiento.
# Valores ajustados
fitted_values = pd.Series(fitted.fittedvalues, copy=True)
# Conversion a suma cumulativa
diff_cumulutive_sum = fitted_values.cumsum()
log_predictions = pd.Series(first_diff.sale_dollars.iloc[0], index=first_diff.index)
log_predictions = b.add(a, fill_value=0)
# Transfomracion inversa de logaritmo es exponenciar la serie
predictions = np.exp(b)
ax = df.sale_dollars.loc['2014':].plot(label='Sale Dollars')
predictions.loc['2014':].plot(ax=ax, label='Forecast ', alpha=.7,
figsize=(14,7), color='red')
plt.ylabel('Sale dollars', fontsize=14)
plt.xlabel('Año', fontsize=14)
plt.title('Predicción mensual', fontsize=20)
plt.legend()
plt.tight_layout()
plt.show()
