Predictor de precios de mercado - General Motors - Parte II
Como hemos comentado en los posts anteriores, en este set de diferentes notebooks trataremos de hacer un análisis de la serie temporal de los precios de cierre de mercado de la compañía General Motors en la bolsa de NY. Como objetivo se fijará crear un modelo de referencia y tratar de mejorarlo usando redes secuenciales. El objeto de estudio de hoy es:
El segundo punto, Modelos de referencia, estudiaremos principalmente un modelo autorregresivo integrado de promedio móvil o ARIMA, dado que será el que mejor se ajuste a nuestra serie en función del pre-análisis del punto uno. También, se probarán otros métodos como Naïve method y random-forest para ver si nuestras premisas eran correctas.
Forecast: Modelos de Referencia¶
En este punto estudiaremos principalmente dos modelos, uno autorregresivo integrado de promedio móvil o ARIMA, y otro más simple, naïve methods que aparentemente no sea el que mejor se ajuste a las caracteristicas de nuestra serie en función del pre-análisis del punto uno. También, se probarán otros métodos como random-forest para ver si nuestras premisas eran correctas. EL modelo Naive es teóricamente el que mejor se puede ajustar a este tipo de series.
En el post de hoy trataremos:
1.- Modelo ARIMA
2.- Modelos Naïve methods
3.- Modelo random Forest
Modelo Arima¶
Vamos en este punto a realizar un estudio de los datos para pasar despues a entrenar un modelo univariante tradicional conocido como el modelo ARIMA. De manera teorica para el entrenamiento de un Modelo de este tipo tenemos que partir de un conocimiento de la tipologia de nuestra Secuencia, en cuanto a su tendencia (Trend) y su estacionalidad (Seasonal).
Teoricamente este tipo de modelos sólo se puede aplicar a series estacionarias por ellos tendremos que transformar la serie en ello.
En estadística y econometría, en particular en series temporales, ARIMA es un modelo autorregresivo integrado de promedio móvil o ARIMA (acrónimo del inglés autoregressive integrated moving average) es un modelo estadístico que utiliza variaciones y regresiones de datos estadísticos con el fin de encontrar patrones para la realizacion de predicciones en 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, es por tanto un modelo Univariante.
Fue desarrollado a finales de los sesenta del siglo XX. Box y Jenkins (1976) lo sistematizaron.
El modelo ARIMA necesita identificar los coeficientes y número de regresiones que se utilizarán. La metodología de Box y Jenkins para entrenar un modelo ARIMA o SARIMA( dependiendo de si nuestra serie presenta estacionaridad) se resume en 4 pasos que nos han parecido interesante enumerar y que en la medida de los posible hemos ido siguiendo.
- Identificar el posible modelo ARIMA a seguir. Para ello primeramente deberemos averiguar las transformaciones a realizar para que nuestra serie sea estacionaria. Una vez hecho esto, deberemos encontrar los posibles valores de p y q correspondientes a su estructura autorregresiva (AR) y de media móvil (MA)
- Seleccionado provisionalmente un modelo para la serie estacionaria, se pasa a la segunda etapa de estimación, donde los parámetros AR y MA del modelo se estiman por máxima verosimilitud y se obtienen sus errores estándar y los residuos del modelo
- Se comprueba que los residuos contienen un formato de ruido blanco
- Realización de predicciones con el modelo
1.- Identificar el posible modelo ARIMA a seguir.¶
Ya tenemos identificado el primer parametro referente a la diferenciacion (según nuestro punto anterior de análisis de datos), en nuestro caso de primer orden. Para definir el modeo arima necesitamos establecer el parametro del modelo AR(p)- Modelo Autorregresivo. Como vimos en las graficas anteriores no presentaba autocorrelaciones,todos los valores caian dentro de la zona de confianza, hemos querido establecer un modelo arima inicial .Igualmente formas establecemos unos valores inciales y deberemos de ir probando, hasta encontrar los parametros optimos. Y vemos como se comporta.
from statsmodels.tsa.arima_model import ARIMA
model_ar1 = ARIMA(data['4. close'].values, order=(1,1,0))
model_fit_ar1 = model_ar1.fit(disp=0)
print(model_fit_ar1.summary())
2 y 3.- Etapa de estimación¶
En esta etapa comprobamos que no existan autocorrelaciones con el Tests Ljung-Box y Box-Pierce y que nt se trata de ruido blanco.
import numpy as np
from statsmodels.stats.diagnostic import acorr_ljungbox
LBBP_df = pd.DataFrame(np.column_stack(acorr_ljungbox(model_fit_ar1.resid,boxpierce=True, lags=14)),
columns=['LB_statistic','LB_pvalue', 'BP_statistic','BP_pvalue'])
display(LBBP_df)
Podemos observar que los demás valores presentan no decaen tan rápidamente después del primero como cabría esperar, por ello pobraremos a modificar estos valores.
residuals = pd.DataFrame(model_fit_ar1.resid)
residuals.plot()
plt.show()
residuals.plot(kind='kde')
plt.show()
print(residuals.describe())
Como podemos ver, en la componente estacionaria, esta cercana al ruido blanco, con una media cercana a 0 y una distribucion mas o menos normal.
3B.- Búsqueda de metaparámetros.¶
Una vez tenemos definido que cumple con los puntos dos y tres vamos a probar a ajustar los parámetros de AR(p) y AM(q) para ver si obtenemos un modelo que ajuste mejor.
Para ellos vamos a seguir dos procesos, uno manual y otro automático. En el proceso manual iremos probando a ajustar cada vez un modelo, ver que cumple los requisitos y probar una predicción. En el proceso automático haremos uso de autoArima.
Proceso manual¶
Dado que el proceso es repetitivo, simplemente ejecutar todas las celdas anteriores y hacer el forecast, no pondremos en código todo el proceso de búsqueda de metaparámetros de forma manual. La lista de los modelos probados, cumpliendo los requisitos son:
Con el modelo (5,1,0) error (mse) de 0.4614938000026067
Con el modelo (5,1,1) error (mse) de 0.3895853732933729
Con el modelo (7,2,1) error (mse) de 0.4661530933833349
Con el modelo (1,1,0) error (mse) de 0.4601377978900837
Con el modelo (8,1,1) error (mse) de 0.4634372479003837
Con el modelo (8,1,0) error (mse) de 0.4940168038836873
Por ello nos quedaremos con el modelo ARIMA (5,1,1)
from statsmodels.tsa.arima_model import ARIMA
model_ar1 = ARIMA(data['4. close'].values, order=(5,1,1))
model_fit_ar1 = model_ar1.fit(disp=0)
print(model_fit_ar1.summary())
Proceso automático¶
Dado que la búsqueda de estos parámetros puede ser algo tedioso, se creo una librería para buscar por fuerza bruta y algunos indicadores el mejor modelo. Aunque en nuestro caso cabe decir que no ha sido capaz de encontrarlo, nos ha sugerido modelos con peores resultados que el encontrado de forma manual.
# Ejercicio Construir una función que divida la serie original
#en training y test sets, con un percentaje dado como input
def training_test_split(series, perc):
#print("samples total=", len(series))
n_train = int(len(series) * perc)
#print("samples train=", n_train)
n_test = len(series) - n_train
#print("samples test", int(n_test))
train_s = list(series[:n_train])
test_s = list(series[n_train:])
return train_s, test_s
Utilizaremos la librería de pmdarima.
!pip install -q pmdarima
Probamos sin componente estacional
import pmdarima as pm
train, test = training_test_split(data['4. close'].values, 0.90)
model = pm.auto_arima(train, test='adf',max_p=30, max_q=10 , seasonal=False,
trace=True,error_action='ignore', suppress_warnings=True, stepwise=True)
print(model.summary())
Si probamos con componente estacional
model = pm.auto_arima(train, test='adf',max_p=30, max_q=10 , seasonal=True,
trace=True,error_action='ignore', suppress_warnings=True, stepwise=True)
print(model.summary())
Como los modelos sugerido ARIMA(0,1,0) (serie random walk) y ARIMA(2,1,2) nos da peores resultados (en predicción error MSE) que el encontrado manualmente, ARIMA(5,1,1) nos quedaremos con este.
Aunque estos modelos de Autoarima mejoran el AIC no significa que siempre sea un mejor modelo. [https://otexts.com/fpp2/non-seasonal-arima.html]
Por ello nos quedaremos con el modelo encontrado de forma manual ARIMA(5,1,1)
4.- Forecast¶
model = ARIMA(data['4. close'].values, order=(1,1,0))
#order=(5,1,1))
model_fit = model.fit(disp=0)
model_fit.plot_predict(dynamic=False)
plt.show()
Cuando establece dynamic = False, los valores rezagados en la muestra se utilizan para la predicción.
Es decir, el modelo se entrena hasta el valor anterior para hacer la siguiente predicción. Esto puede hacer que el pronóstico ajustado y los datos reales se vean artificialmente buenos.
from sklearn.metrics import mean_squared_error
train, test = training_test_split(data['4. close'].values, 0.90)
forecast_l =[]
history = [x for x in train]
predictions = list()
for t in range(len(test)):
model = ARIMA(history, order=(5,1,1))
model_fit = model.fit(disp=0)
output = model_fit.forecast()
yhat = output[0]
predictions.append(yhat)
obs = test[t]
history.append(obs)
if t%10==0:
print('predicted=%f, expected=%f' % (yhat, obs))
from matplotlib import pyplot
from math import sqrt
error = mean_squared_error(test, predictions)
RMSE = sqrt(error)
print('Test MSE: %.3f' % error)
print('Test RMSE: %.3f' % RMSE)
# plot
pyplot.plot(test)
pyplot.plot(predictions, color='red')
pyplot.show()
Naïve methods¶
En algunas series temporales de precios, la mejor predicción es la última que tenemos; esto se produce por ejemplo cuando la serie sigue un modelo de Random Walk, como probablemente sea nuestro caso.[https://otexts.com/fpp2/, Capítulo 3.1]. Este tipo de algoritmo, aún a pesar de ser muy simples son extremadamente eficaces en algunas ocasiones, también de ahí su nombre de Naïve methods.
Para los modelos simples se puede hacer una persistencia diaria, semanal o mensual.
from sklearn.metrics import *
from math import sqrt
def smape_kun(y_true, y_pred):
mse = mean_squared_error(y_true, y_pred)
rmse = sqrt(mean_squared_error(y_true, y_pred))
mape = mean_absolute_error(y_true, y_pred)
smape = mean_squared_log_error(y_true, y_pred)
print('MSE: %.2f'% (mse))
print('RMSE: %.2f'% (rmse))
print('MAPE: %.2f %% \nSMAPE: %.2f'% (mape,smape), "%")
def evaluate_model(model_func, train, test):
history = [x for x in train]
predictions = list()
for i in range(len(test)):
yhat_sequence = model_func(history)
predictions.append(yhat_sequence)
history.append(test[i])
smape_kun(test, predictions)
return predictions
def daily_persistence(history):
return history[-1]
def weekly_persistence(history):
return history[-7]
def month_persistence(history):
return history[-30]
models = dict()
models['daily'] = daily_persistence
models['weekly'] = weekly_persistence
models['monthly'] = month_persistence
predictions = dict()
train_s, test_s = training_test_split(data['4. close'], 0.8)
for name, func in models.items():
print('---->',name)
prediction = evaluate_model(func, train_s, test_s)
predictions[name] = prediction
Podemos observar que en la persistencia diaría, el error de MSE y RMSE no son tan bajos como en ARIMA pero no estan tan alejados.
Esto significa que la serie probablemente se acerque a Random Walk pero bajo algunos ARIMAS concretos es capaz de mejorarla. Si siguierámos ejecutando el modelo ARIMA durante una ventana de tiempo, su error probablemente se acabaría acercando a este o empeorándolo.
Random Forest¶
Este algoritmo desarrollado por Leo Breiman y Adele Cutler consiste en una combinación de árboles predictores con una distribución común, en el cual cada uno de estos depende de un vector aleatorio. El error general de este modelo converge hasta conseguir el número de árboles máximo o un conjunto demasiado grande.
El uso de una selección aleatoria de características para dividir cada nodo produce tasas de error que son usualmente muy favorables frente a otros modelos, por ello se ha convertido en uno de los modelos de referencia y más utilizados actualmente. Estas ideas son aplicables tanto a clasificación como a regresión.
Este algoritmo ha demostrado ser robusto a la hora de hacer predicciones en precios de mercado, sobre todo cuando estas tienen diferentes ruidos, por ello la utilizaremos como modelo de referencia para nuestra comparativa.
Random Forest sin normalizar¶
Para el modelo de random forest, dado que no es modelo que trabaja sólo con series temporales conservaremos los datos de entrada de volumen, precio min. La variable objetivo será predecir el precio del día siguiente.
from sklearn.model_selection import train_test_split
from sklearn.metrics import *
from math import sqrt
#Crear la variable objetivo del precio del día siguiente
x_data = data
x_data = x_data.drop('date', axis=1)
y_data = pd.DataFrame(x_data['4. close']).shift(-1)
y_data = y_data.ffill(axis = 0)
x_train, x_valid, y_train, y_valid = train_test_split(x_data, y_data, test_size=0.33, shuffle=False)
x_train = x_train.to_numpy()
x_valid = x_valid.to_numpy()
y_train = y_train.to_numpy()
y_valid = y_valid.to_numpy()
#Usamos en la entrada X todas las variables de entrada, como volumen precio máximo, mínimo.
print('Xtrain_dim:', x_train.shape)
print('Ytrain_dim:', y_train.shape)
print(x_train)
from sklearn.ensemble import RandomForestRegressor
from sklearn.datasets import make_regression
regressor = RandomForestRegressor(
random_state=0,
n_estimators=300,
max_features='auto',
max_depth=20,
min_samples_split=2,
min_samples_leaf=4,
bootstrap=True
)
regressor.fit(x_train, y_train)
print('MSE en entrenamiento es: ', mean_squared_error(y_train, regressor.predict(x_train)))
print('MSE en entrenamiento es: ', mean_squared_error(y_valid, regressor.predict(x_valid)))
y_predict = regressor.predict(x_valid)
print('MSE: ', mean_absolute_error(y_valid, y_predict))
print('RMSE: ', sqrt(mean_squared_error(y_valid, y_predict)))
print('Variance: ', explained_variance_score(y_valid, y_predict))
import numpy as np
plt.figure(figsize=(12, 6))
plt.plot(y_valid, label='real data')
plt.plot(y_predict, label='pred data')
Podemos ver a simple vista que los resultados estan muy lejos de los modelos de referencia anteriores, parace que bajo algunos valores no es capaz de predecir el precio. Probemos a normalizar estos datos en busca de obtener alguna mejoría.
Random Forest normalizado¶
Hacemos uso del Standar Scaler para normalizar estos datos.
from sklearn.preprocessing import *
scaler = StandardScaler()
x_train_est = scaler.fit_transform(x_train)
y_train_est = scaler.fit_transform(y_train)
regressor_standarized = RandomForestRegressor(
random_state=0,
n_estimators=300,
max_features='auto',
max_depth=20,
min_samples_split=2,
min_samples_leaf=4,
bootstrap=True
)
regressor_standarized.fit(x_train_est, y_train_est)
x_test_est = scaler.fit_transform(x_valid)
y_test_est = scaler.fit_transform(y_valid)
print('MSE en entrenamiento es: ', mean_squared_error(y_train_est, regressor.predict(x_train_est)))
print('MSE en validación es: ', mean_squared_error(y_test_est, regressor.predict(x_test_est)))
x_predict = regressor_standarized.predict(x_test_est)
y_predict_trans = scaler.inverse_transform(x_predict)
y_valid_trans = scaler.inverse_transform(y_test_est)
print('MSE: ', mean_squared_error(y_valid_trans, y_predict_trans))
print('RMSE: ', sqrt(mean_squared_error(y_valid_trans, y_predict_trans)))
print('Variance: ', explained_variance_score(y_valid_trans, y_predict_trans))
plt.figure(figsize=(12, 6))
x_predict_trans = scaler.inverse_transform(x_predict)
y_predict_trans = scaler.inverse_transform(y_test_est)
plt.plot(y_predict_trans, label='real data')
plt.plot(x_predict_trans, label='pred data')
Los resultado mejoran y nos arrojan unos errores de MSE 0.45 Y RMSE 0.67, que no estan tan alejados de aquellos errores el ARIMA de 0.36 y 0.63. El modelo de Random Forest no mejora el resultado, pero queda claro que es un modelo a tener en cuenta en el estudio de series temporales. Como antes comentabámos este tipo de modelos son muy robustos frente series temporales ya que muchas veces presentan diferentes ruidos.
Resumen¶
Entre los diferentes modelos probados sin duda el ARIMA ha sido el que mejores resultados ha tenido, pero los demás tampoco han resultado ser mucho peores a la hora de esta predicción, probablemente dado porque la serie temporal se parece mucho a un random - walk, y por ello NaÏve methods y random forest tienen tan buen desempeño.
- ARIMA (5,1,1) => 0.36 (MSE), 0.63 (RMSE)
- NAÏVE DAILY => 0.46 (MSE), 0.68 (RMSE)
- RANDOM FOREST ESTANDARIZADO => 0.45 (MSE), 0.67 (RMSE)
Nos quedaremos entonces con el modelos de ARIMA(5,1,1) para intentar mejorarlo con redes secuenciales
Comentar