Implementación de Python
Antes de comenzar, crea una cuenta con Datos económicos de la Reserva Federal (FRED) y obtener una clave API usando Este enlaceTenga en cuenta que este producto utiliza la API FRED® pero no está avalado ni certificado por el Banco de la Reserva Federal de St. Louis.
Comenzamos instalando y cargando los módulos necesarios.
import requests
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
from datetime import datetime, timedelta
from pandas.tseries.offsets import MonthEnd
from sklearn.decomposition import PCA
from sklearn.preprocessing import StandardScaler
from sklearn.metrics import mean_absolute_percentage_error
from statsmodels.tsa.stattools import adfuller
from statsmodels.tsa.stattools import grangercausalitytestsfrom pmdarima import auto_arima
A continuación, crearemos una función personalizada para leer los datos utilizando la API de FRED.
FRED_API_KEY = '__YOUR_API_KEY__'# Function to read data from FRED API
def get_fred_data(data_id, data_name):
response = requests.get(f'https://api.stlouisfed.org/fred/series/observations?series_id={data_id}&api_key={FRED_API_KEY}&file_type=json')
df = pd.DataFrame(response.json()['observations'])[['date', 'value']].rename(columns={'value': data_name})
df[data_name] = pd.to_numeric(df[data_name], errors='coerce')
df['date'] = pd.to_datetime(df['date']) + MonthEnd(1)
df.set_index('date', inplace=True)
df.index.freq='M'
return df
Ahora, leamos nuestros datos y almacenémoslos en un marco de datos de pandas.
dependent_timeseries_id = 'MRTSSM4453USN'
dependent_timeseries_name = 'liquor_sales'potential_leading_indicators = {
'USACPIALLMINMEI': 'consumer_price_index',
'PCU44534453': 'liquor_ppi',
'DSPIC96': 'real_income',
'LFWA25TTUSM647S': 'working_age_population',
'UNRATENSA': 'unemployment_rate',
'LNU01300000': 'labor_force_participation',
'A063RC1': 'social_benefits',
'CCLACBM027NBOG': 'consumer_loans',
}
# Read dependent time series
timeseries = get_fred_data(dependent_timeseries_id, dependent_timeseries_name)# Join timeseries with potential leading indicators
for data_id, data_name in potential_leading_indicators.items():
df = get_fred_data(data_id, data_name)
timeseries = timeseries.join(df)
# We will start our analysis from Jan-2010
timeseries = timeseries['2010':]
# add month we want to predict liquor_sales
timeseries = timeseries.reindex(timeseries.index.union([timeseries.index[-1] + timeseries.index.freq]))
timeseries
Un análisis visual rápido de nuestros datos muestra que nuestra serie temporal dependiente (ventas de licores) sigue más o menos el mismo ciclo cada 12 meses. Usaremos este ciclo de 12 meses como parámetro más adelante en nuestro pronóstico de series temporales.
timeseries[dependent_timeseries_name].plot(figsize=(20,8));
Antes de comprobar la causalidad, debemos confirmar la estacionariedad de los datos de nuestra serie temporal. Para lograrlo, utilizaremos la prueba de Dickey-Fuller aumentada. Si nuestro conjunto de datos no supera esta prueba de estacionariedad, debemos emplear un enfoque de diferenciación recursiva hasta que satisfaga los criterios de la prueba.
# create a copy of the timeseries to use for tests. Be sure to exclude the additional row we added in the previous task
timeseries_for_gc_tests = timeseries[:-1]
all_cols = timeseries_for_gc_tests.columnsstationary_cols = []
diff_times = 0
while True:
# Test for stationarity
for col in all_cols:
adf, pvalue, lagsused, observations, critical_values, icbest = adfuller(timeseries_for_gc_tests[col])
if pvalue <= 0.05:
stationary_cols.append(col)
# Difference the time series if at least one column fails the stationary test
if set(stationary_cols) != set(all_cols):
timeseries_for_gc_tests = timeseries_for_gc_tests.diff().dropna()
diff_times += 1
stationary_cols = []
else:
print(f'No of Differencing: {diff_times}')
break
Ahora que hemos cargado nuestros datos de series de tiempo en un marco de datos de pandas, y pasa la prueba de estacionariedad, probamos la causalidad utilizando la prueba de causalidad de Granger.
maxlag = 6 # represents the maximum number of past time periods to look for potential causality. We cap ours at 6 months
leading_indicators = []for x in all_cols[1:]:
gc_res = grangercausalitytests(timeseries_for_gc_tests[[dependent_timeseries_name, x]], maxlag=maxlag, verbose=0)
leading_indicators_tmp = []
for lag in range(1, maxlag+1):
ftest_stat = gc_res[lag][0]['ssr_ftest'][0]
ftest_pvalue = gc_res[lag][0]['ssr_ftest'][1]
if ftest_pvalue <= 0.05:
leading_indicators_tmp.append({'x': x, 'lag': lag, 'ftest_pvalue': ftest_pvalue, 'ftest_stat': ftest_stat, 'xlabel': f'{x}__{lag}_mths_ago'})
if leading_indicators_tmp:
leading_indicators.append(max(leading_indicators_tmp, key=lambda x:x['ftest_stat']))
# Display leading indicators as a dataframe
pd.DataFrame(leading_indicators).reset_index(drop=True).reset_index(drop=True)
De nuestras pruebas, podemos ver que las ventas de licor del mes actual se ven afectadas por el Índice de Precios al Consumidorᵈ² y los Préstamos al Consumidorᵈ¹⁰ de hace 2 meses; y la Participación de la Fuerza Laboralᵈ⁷ de hace 6 meses.
Ahora que hemos establecido nuestros indicadores principales, cambiaremos sus registros para que sus cifras rezagadas estén en la misma fila que los datos actuales de ventas de licor que ellos “causan”.
# shift the leading indicators by their corresponding lag periods
for i in leading_indicators:
timeseries[i['xlabel']] = timeseries[i['x']].shift(periods=i['lag'], freq='M')# select only the dependent_timeseries_name and leading indicators for further analysis
timeseries = timeseries[[dependent_timeseries_name, *[i['xlabel'] for i in leading_indicators]]].dropna(subset=[i['xlabel'] for i in leading_indicators], axis=0)
timeseries
A continuación, escalamos nuestros datos para que todas las características estén dentro del mismo rango y luego aplicamos PCA para eliminar la multicolinealidad entre nuestros indicadores principales.
# Scale dependent timeseries
y_scaler = StandardScaler()
dependent_timeseries_scaled = y_scaler.fit_transform(timeseries[[dependent_timeseries_name]])# Scale leading indicators
X_scaler = StandardScaler()
leading_indicators_scaled = X_scaler.fit_transform(timeseries[[i['xlabel'] for i in leading_indicators]])
# Reduce dimensionality of the leading indicators
pca = PCA(n_components=0.90)
leading_indicators_scaled_components = pca.fit_transform(leading_indicators_scaled)leading_indicators_scaled_components.shape
Finalmente, podemos construir nuestro modelo SARIMAX con la ayuda de auto_arima. Para el propósito de esta implementación, dejaremos todos los parámetros con sus valores predeterminados, con la excepción del indicador de estacionalidad y el número de períodos en cada ciclo (m).
Entrenaremos nuestro modelo utilizando los datos de series temporales hasta ‘2024–05–31’, probaremos con los datos ‘2024–06–30’ y luego predeciremos las ventas de licor ‘2024–07–31’.
# Build SARIMAX model
periods_in_cycle = 12 # number of periods per cycle. In our case, its 12 months
model = auto_arima(y=dependent_timeseries_scaled[:-2], X=leading_indicators_scaled_components[:-2], seasonal=True, m=periods_in_cycle)
model.summary()
# Forecast the next two periods
preds_scaled = model.predict(n_periods=2, X=leading_indicators_scaled_components[-2:])
pred_2024_06_30, pred_2024_07_31 = np.round(y_scaler.inverse_transform([preds_scaled]))[0]print("TEST\n----")
print(f"Actual Liquor Sales for 2024-06-30: {timeseries[dependent_timeseries_name]['2024-06-30']}")
print(f"Predicted Liquor Sales for 2024-06-30: {pred_2024_06_30}")
print(f"MAPE: {mean_absolute_percentage_error([timeseries[dependent_timeseries_name]['2024-06-30']], [pred_2024_06_30]):.1%}")
print("\nFORECAST\n--------")
print(f"Forecasted Liquor Sales for 2024-07-31: {pred_2024_07_31}")
Siguiendo el proceso paso a paso, pronosticamos la cifra de ventas de licores para julio de 2024 con un MAPE estimado de solo 0,4%.
Para mejorar aún más la precisión de nuestra predicción, podemos explorar la posibilidad de agregar más indicadores líderes potenciales y ajustar los modelos utilizados.
Conclusión
Los indicadores adelantados, como hemos visto, sirven como señales tempranas de tendencias futuras, lo que proporciona una ventaja crucial para anticipar los cambios antes de que se materialicen por completo. Al aprovechar técnicas como las pruebas de causalidad de Granger para identificar series de indicadores adelantados e incorporarlas a los modelos de pronóstico, podemos mejorar significativamente la precisión y la solidez de nuestras predicciones.
