Series temporales y combinación de clasificadores¶
Esta práctica está dividida en dos partes:
- En el primer ejercicio veremos como descomponer y componer series temporales para realizar predicciones a futuro.
- En el segundo ejercicio estudiaremos diferentes métodos de combinación de clasificadores.
import pickle
import scipy.stats
import numpy as np
import pandas as pd
import matplotlib as mpl
from sklearn import svm
from sklearn import ensemble
from matplotlib import pyplot as plt
from sklearn.decomposition import PCA
from sklearn.metrics import accuracy_score
from sklearn.preprocessing import OneHotEncoder
from sklearn.tree import DecisionTreeClassifier
from sklearn.linear_model import LinearRegression
from sklearn.neighbors import KNeighborsClassifier
from sklearn.model_selection import cross_val_score
from sklearn.model_selection import train_test_split
!pip install statsmodels
from statsmodels.tsa.statespace.sarimax import SARIMAX
from statsmodels.graphics.tsaplots import plot_pacf, plot_acf
%matplotlib inline
1. Series temporales¶
En este primer ejercicio trabajaremos las series temporales. Para ello usaremos el dataset AirPassangers que contiene información del número de vuelos que se realizaron a lo largo de muchos años.
Empezaremos leyendo los datos y observando gráficamente su distribución. Como se puede apreciar es un claro caso de serie temporal, con heterocedasticidad, tendencia, periodo y ruido. A lo largo de este ejercicio trataremos cada uno de estos puntos.
data = pd.read_csv('AirPassengers.csv', parse_dates=['Month'], index_col='Month', header=0)
data.head()
data.plot(figsize=(15, 5))
Antes de empezar a tratar las diferentes componenetes de una serie temporal, eliminaremos del dataset original un par de años de datos. Así cuando hagamos una predicción a futuro podremos comprobar si se ajusta a los datos reales.
TEST_SIZE = 24
train, test = data.iloc[:-TEST_SIZE], data.iloc[-TEST_SIZE:]
x_train, x_test = np.array(range(train.shape[0])), np.array(range(train.shape[0], data.shape[0]))
train.shape, x_train.shape, test.shape, x_test.shape
fig, ax = plt.subplots(1, 1, figsize=(15, 5))
ax.plot(x_train, train)
ax.plot(x_test, test)
1.1 Componentes de la serie temporal¶
1.1.a. Heterocedasticidad¶
Se dice que una serie temporal tiene heterocedasticidad cuando la variancia varía con el tiempo (https://es.wikipedia.org/wiki/Heterocedasticidad). En nuestro caso, observamos que tenemos heterocedasticidad, ya que la amplitud de onda varía con el tiempo. En este primer apartado debéis eliminar la heterocedasticidad de la serie temporal. Es decir, que la diferencia entre el mínimo y el máximo de la estacionalidad (anual) sea más o menos la misma a lo largo del tiempo.
train_log = np.log(train)
test_log = np.log(test)
print('logarítmica')
fig, ax = plt.subplots(1, 1, figsize=(15, 5))
ax.plot(x_train, train_log)
ax.plot(x_test, test_log, label='test')
ax.set_title('logarítmica')
plt.legend([l1, l2],["train", 'test'])
print('exponencial')
train_exp = np.exp(train)
fig, ax = plt.subplots(1, 1, figsize=(15, 5))
ax.plot(x_train, train_exp)
ax.set_title('exponencial')
print('raiz cuadrada')
train_sqrt = np.sqrt(train)
fig, ax = plt.subplots(1, 1, figsize=(15, 5))
ax.plot(x_train, train_sqrt)
ax.set_title('raiz cuadrada')
Con el resultado obtenido, considero que la mejor opción para eliminarla heterocedasticidad en el dataset es la función logaritmica.
1.1.b. Tendencia¶
La tendencia es el comportamiento que tienen los datos a largo plazo (https://miro.medium.com/max/1872/1*rDQL2fAp_X_dgAHNZuwRfw.png). En nuestra serie temporal tenemos una tendencia lineal creciente. En este apartado debéis eliminar la tendencia, quedando una serie temporal con tendencia constante.
model = LinearRegression()
model.fit(x_train.reshape(-1, 1), train_log.values)
train_log_trend = model.predict(x_train.reshape(-1, 1))
fig, ax = plt.subplots(1, 1, figsize=(15, 5))
ax.plot(x_train, train_log)
ax.plot(x_test, test_log)
ax.plot(x_train, train_log_trend)
plt.legend([l1, l2],["train", 'test'])
train_log_trend = (train_log - model.predict(x_train.reshape(-1, 1))).values.squeeze()
test_log_trend = (test_log - model.predict(x_test.reshape(-1, 1))).values.squeeze()
fig, ax = plt.subplots(1, 1, figsize=(15, 5))
ax.plot(x_train, train_log_trend)
ax.plot(x_test, test_log_trend)
plt.legend([l1, l2],["train", 'test'])
1.1.c. Estacionalidad¶
Definimos estacionalidad como la variación cíclica que se produce en los datos (https://es.wikipedia.org/wiki/Estacionalidad). En este apartado se debe encontrar y eliminar la estacionalidad.
1.1.c.a. Encontrar el periodo de la estacionalidad¶
corr_coefficients = []
for i in range(1, 24):
corr_coefficients.append(np.corrcoef(train_log_trend[i:], train_log_trend[:-i])[0,1])
fig, ax = plt.subplots(1, 1, figsize=(15, 5))
ax.plot(range(len(corr_coefficients)), corr_coefficients)
1.1.c.b. Aplicar un modelo SARIMA¶
Para tratar la estacionalidad aplicaremos un modelo SARIMA. Las siglas corresponden a "stationality ARIMA", es decir, un modelo ARIMA con estacionalidad (la que acabamos de encontrar).
El modelo ARIMA nos va a permitir tratar el ruido que queda al eliminar la heterocedasticidad, tendencia y estacionalidad. Para hacerlo tiene en cuenta las siguientes componentes:
- AR: auto-regressive, se denomina p. Tiene en cuenta la correlación con sus lags, es decir, mira si las observaciones pasadas afectan para calcular el siguiente punto.
- I: integrated, se denomina d. Es el orden de diferenciación, en nuestro caso no es necesario, puesto que la serie ya no tiene heterocedasticidad ni tendencia.
- MA: moving-average, se denomina q. Tiene en cuenta la correlación con los lags de los errores, es decir, una vez se ha aplicado el modelo, mira los errores del modelo versus los datos reales.
Para más información sobre los ARIMA podéis consultar este link: https://otexts.com/fpp2/arima.html
Los modelos SARIMA dependen de varios parámetros (p,d,q)(P,D,Q)s, donde los primeros (minúsculas) corresponden a AR, I, MA de la serie normal, y los segundos (mayúsculas) a AR, I, MA con estacionalidad.
Para encontrar estos parámetros tenemos que mirar los gráficos PACF (partial autocorrelation function) y ACF (autocorrelation function).
El gráfico PACF nos determina el parámetro p, es decir, la AR.
# Supongamos que los datos sin heterocedasticidad ni tendencia se llaman train_log_trend
plot_pacf(train_log_trend)
plt.show()
Para leer este gráfico, simplemente nos tenemos que fijar en los valores que salen fuera del intervalo de confianza (zona azul).
Nada mas empezar hay dos valores fuera del intervalo. De todos modos el primero no se debe tener en cuenta, puesto que mira la correlación de un valor consigo mismo, y éste siempre será 1. Si no tenemos en cuenta este primer valor, solo hay un valor fuera del intervalo de confianza, con lo cual p = 1.
Cuando se repite el ciclo, es decir, a partir del valor 11, vuelven a haber dos valores fuera del intervalo, con lo cual P = 2.
Veamos ahora el gráfico ACF, este determinará el valor de q.
# Supongamos que los datos sin heterocedasticidad ni tendencia se llaman train_log_trend
plot_acf(train_log_trend)
plt.show()
Este gráfico se lee igual que el anterior. Nada más empezar hay un valor fuera del intervalo, con lo cual q = 1. Cuando se repite el ciclo, es decir, a partir del valor 11, hay tres valores fuera del intervalo, con lo cual Q = 3.
Como los datos no tienen ni tendencia ni heterocedasticidad, d = D = 0.
Como hemos visto en el apartado anterior, el ciclo es 12, con lo cual s = 12.
import warnings
warnings.filterwarnings("ignore") # specify to ignore warning messages
import statsmodels.api as sm
p=1
q=1
d=0
P=2
D=0
Q=3
m=12
#SARIMA(p,d,q)(P,D,Q)m
mod = sm.tsa.statespace.SARIMAX(train_log_trend,
order=(p,d,q),
seasonal_order=(P,D,Q, m))
results = mod.fit()
print('----Resultado del modelo----')
print(results.summary())
results.plot_diagnostics(figsize=(15, 12))
plt.show()
print('En este punto, únicamente utilizo la parte test para ver si el modelo es fiel a la realidad')
start = len(train)
end = len(train) + len(test) - 1
predictions = results.predict(start, end,
typ = 'levels')
predictions= results.predict(start =start, end=end,exog=x_test)
forecast_1= results.forecast(steps=len(test) , exog=x_test)
train_size=len(train)
act= pd.DataFrame(test_log_trend)
predictions=pd.DataFrame(predictions)
predictions.reset_index(drop=True, inplace=True)
predictions.index=pd.DataFrame(x_test).index
predictions['Actual'] = act
predictions.rename(columns={0:'Pred'}, inplace=True)
predictions['Actual'].plot(figsize=(20,8), legend=True, color='blue')
predictions['Pred'].plot(legend=True, color='red', figsize=(20,8),title='Predición sobre el test')