Pronóstico de una serie de tiempo usando redes neuronales — 41:42 min#

41:42 min | Última modificación: Abril 14, 2021 | YouTube

Uno de los casos más comunes de aplicación de las redes neuronales es el pronóstico de series de tiempo. En esta lección, se realiza el pronóstico de la serie correspondiente al número total de personas empleadas mensualmente como mano de obra en el condado de Sutter (California) entre 1946:1 y 1966:12 usando un perceptrón multicapa y se evalua el efecto de la transformación de la serie en la precisión de los resultados.

[1]:

import warnings

warnings.filterwarnings("ignore")

Definición del problema#

La siguiente serie registra el número total de personas empleadas mensualmente como mano de obra en el condado de Sutter (California) entre 1946:1 y 1966:12. Tiene un total de 252 observaciones. Se desea contruir un modelo de pronóstico usando redes neuronales artificiales. Se deben usar los primeros 228 (de 1946:1 a 1964:12) datos para especificar el modelo y los 24 datos restantes (1965:1 a 1966:12) para evaluar el pronóstico un mes adelante por fuera de la muestra de calibración de los parámetros.

[2]:

[3]:

import matplotlib.pyplot as plt

plt.figure(figsize=(14, 5))
plt.plot(data, ".-k")
plt.grid()

# linea vertical para dividir el entrenamiento
# del pronóstico
plt.plot([len(data) - 24, len(data) - 24], [min(data), max(data)], "--", linewidth=2);

../_images/43_modelos_de_redes_neuronales_08_pronostico_series_de_tiempo_6_0.png

Análisis exploratorio#

La inspección visual revela que la serie tiene un componente cíclico con un año de periodo, lo cual se explica porque la contratación de mano de obra es estacionaria (por periodos en el año); la serie presenta también una tendencia creciente en el largo plazo, cuya pendiente cambia aproximadamente sobre la mitad derecha de la gráfica.

La función de autocorrelación (ACF) simple, definida como:

\text{acf}(k) = \frac{1}{T-k} \sum (z_{t+k} - \bar{z})(z_t - \bar{z})

donde:

\bar{z} = \frac{1}{T} \sum_t z_t

presenta el siguiente comportamiento:

[4]:

from statsmodels.tsa.stattools import acf, pacf


def acf_plot(z):
    acf_data = acf(z, fft=False)
    plt.stem(range(len(acf_data)), acf_data, use_line_collection=True)
    plt.ylim(-1, 1),
    plt.grid()


def pacf_plot(z):
    pacf_data = pacf(z)
    plt.stem(range(len(pacf_data)), pacf_data, use_line_collection=True)
    plt.ylim(-1, 1)
    plt.grid()


acf_plot(data)

../_images/43_modelos_de_redes_neuronales_08_pronostico_series_de_tiempo_14_0.png

que se interpreta de la siguiente manera:

El comportamiento cíclico corresponde a la componente periodica dentro del año, y tiene un periodo de 12 meses.
La tendencia decreciente que muere en el largo plazo representa la tendencia de largo plazo.

Existen dos tendencias para modelar la serie usando una red neuronal artificial:

Usar directamente los datos de la serie y dejar que el modelo capture tanto la tendencia como el patrón cíclico.
Remover la tendencia y el ciclo, y la serie que queda es pronósticada usando el modelo de redes neuronales artificiales.

En el pronóstico de series de tiempo, el valor actual de la serie z_t es función de sus valores pasados z_{t-1}, …, z_{t-P}; estos últimos son la entrada a la red neuronal.

Modelo 1: Pronostico de la serie sin transformar#

En este caso se tomará a P = 13. En la práctica deberían tantearse diferentes valores de P.

[5]:

#
# Como primer paso se escala la serie al intervalo [0, 1]
# ya que esto facilita el entrenamiento del modelo
#
import numpy as np
from sklearn.preprocessing import MinMaxScaler

# crea el transformador
scaler = MinMaxScaler()

# escala la serie
data_scaled = scaler.fit_transform(np.array(data).reshape(-1, 1))

# z es un array de listas como efecto
# del escalamiento
data_scaled = [u[0] for u in data_scaled]

plt.figure(figsize=(14, 5))
plt.plot(data_scaled, ".-k")
plt.grid()
plt.plot(
    [len(data_scaled) - 24, len(data_scaled) - 24],
    [min(data_scaled), max(data_scaled)],
    "--",
    linewidth=2,
);

../_images/43_modelos_de_redes_neuronales_08_pronostico_series_de_tiempo_20_0.png

Ya que la implementación dispobible en sklearn es para modelos de regresión, se debe armar una matrix donde las variables independientes son z_{t-1}, …, z_{t-P} y la variable dependiente es z_t.

[6]:

P = 13

X = []
for t in range(P - 1, 252 - 1):
    X.append([data_scaled[t - n] for n in range(P)])

observed_scaled = data_scaled[P:]

[7]:

#
# Numero total de patrones al convertir
# los datos a un modelo de regresión
#
len(X)

[7]:

[8]:

from sklearn.neural_network import MLPRegressor

np.random.seed(123456)

H = 1  # Se escoge arbitrariamente

mlp = MLPRegressor(
    hidden_layer_sizes=(H,),
    activation="logistic",
    learning_rate="adaptive",
    momentum=0.0,
    learning_rate_init=0.1,
    max_iter=10000,
)

# Entrenamiento
mlp.fit(X[0:215], observed_scaled[0:215])  # 239 - 24 = 215

# Pronostico
y_scaled_m1 = mlp.predict(X)

plt.figure(figsize=(14, 5))
plt.plot(data_scaled, ".-k")
plt.grid()

#
# No hay pronóstico para los primeros 13 valores
# de la serie
#
plt.plot([None] * P + y_scaled_m1.tolist(), "-r")

#
# linea vertical para dividir el entrenamiento
# del pronóstico. Se ubica en el ultimo dato
# usando para entrenamiento
#
plt.plot(
    [len(data_scaled) - 24, len(data_scaled) - 24],
    [min(data_scaled), max(data_scaled)],
    "--",
    linewidth=2,
);

../_images/43_modelos_de_redes_neuronales_08_pronostico_series_de_tiempo_24_0.png

[9]:

#
# Se desescala para llevar los valores a la escala de los datos originales
#
y_m1 = scaler.inverse_transform([[u] for u in y_scaled_m1])
y_m1 = [u[0] for u in y_m1]

plt.figure(figsize=(14, 5))
plt.plot(data, ".-k")
plt.grid()
plt.plot([None] * P + y_m1, "-r")
plt.plot([len(data) - 24, len(data) - 24], [min(data), max(data)], "--", linewidth=2);

../_images/43_modelos_de_redes_neuronales_08_pronostico_series_de_tiempo_25_0.png

Ejercicio.— Determine el número óptimo de neuronas en la capa oculta y la cantidad de rezagos P. Recuerde que los resultados son sensibles a los parámetros de la optimización y el punto de arranque.

Modelo 2: Pronóstico de la serie transformada#

Este es el procedimiento común usado en Estadística Aplicada y Econometría. Acá se juzga el modelo por la precisión, no por el cumplimiento de criterios estadísticos.

[10]:

# Se remueve la tendencia
data_d1 = [data[t] - data[t - 1] for t in range(1, len(data))]

# En la gráfica queda la componente cíclica
plt.figure(figsize=(14, 5))
plt.plot(data_d1, ".-k")
plt.grid()
plt.plot(
    [len(data_d1) - 24, len(data_d1) - 24],
    [min(data_d1), max(data_d1)],
    "--",
    linewidth=2,
);

../_images/43_modelos_de_redes_neuronales_08_pronostico_series_de_tiempo_29_0.png

[11]:

#
# Se remueve la componente cíclica
# restando z[t] - z[t-12]
#
data_d1d12 = [data_d1[t] - data_d1[t - 12] for t in range(12, len(data_d1))]

plt.figure(figsize=(14, 5))
plt.plot(data_d1d12, ".-k")
plt.grid()
plt.plot(
    [len(data_d1d12) - 24, len(data_d1d12) - 24],
    [min(data_d1d12), max(data_d1d12)],
    "--",
    linewidth=2,
);

../_images/43_modelos_de_redes_neuronales_08_pronostico_series_de_tiempo_30_0.png

[12]:

acf_plot(data_d1d12)

../_images/43_modelos_de_redes_neuronales_08_pronostico_series_de_tiempo_31_0.png

[13]:

pacf_plot(data_d1d12)

../_images/43_modelos_de_redes_neuronales_08_pronostico_series_de_tiempo_32_0.png

A partir de la gráfica del PACF, se deduce que un valor apropiado para P podría ser 12.

[14]:

#
# Se escalan los valores como en el caso anaterior
#
scaler = MinMaxScaler()
data_d1d12_scaled = scaler.fit_transform(np.array(data_d1d12).reshape(-1, 1))
data_d1d12_scaled = [u[0] for u in data_d1d12_scaled]

#
# Se verifica visualmente el escalamiento
#
plt.figure(figsize=(14, 5))
plt.plot(data_d1d12_scaled, ".-k")
plt.grid()
plt.plot(
    [len(data_d1d12_scaled) - 24, len(data_d1d12_scaled) - 24],
    [min(data_d1d12_scaled), max(data_d1d12_scaled)],
    "--",
    linewidth=2,
);

../_images/43_modelos_de_redes_neuronales_08_pronostico_series_de_tiempo_34_0.png

[15]:

#
# La longitud de la serie cambia
#
len(data_d1d12_scaled)

[15]:

[16]:

#
# Se construye la matriz de regresores
#
P = 13

X = []
for t in range(P - 1, len(data_d1d12_scaled) - 1):
    X.append([data_d1d12_scaled[t - n] for n in range(P)])

d = data_d1d12_scaled[P:]

[17]:

#
# Numero total de patrones al convertir
# los datos a un modelo de regresión
#
len(X)

[17]:

[18]:

len(d)

[18]:

[19]:

H = 4  # Se escoge arbitrariamente

np.random.seed(123456)

mlp = MLPRegressor(
    hidden_layer_sizes=(H,),
    activation="logistic",
    learning_rate="adaptive",
    momentum=0.0,
    learning_rate_init=0.002,
    max_iter=100000,
)

# Entrenamiento
mlp.fit(X[0:202], data_d1d12_scaled[0:202])  # 226 - 24 = 202

# Pronostico
y_d1d12_scaled_m2 = mlp.predict(X)


plt.figure(figsize=(14, 5))
plt.plot(data_d1d12_scaled, ".-k")
plt.grid()

# No hay pronóstico para los primeros 13 valores
# de la serie
plt.plot([None] * P + y_d1d12_scaled_m2.tolist(), "-r")

# linea vertical para dividir el entrenamiento
# del pronóstico
plt.plot(
    [len(data_d1d12_scaled) - 24, len(data_d1d12_scaled) - 24],
    [min(data_d1d12_scaled), max(data_d1d12_scaled)],
    "--",
    linewidth=2,
);

../_images/43_modelos_de_redes_neuronales_08_pronostico_series_de_tiempo_39_0.png

[20]:

print(len(y_d1d12_scaled_m2), ",", len(data_d1d12_scaled) - P)

226 , 226

[21]:

#
# La gráfica anterior no es comparativa. Se realizan las
#  transformaciones inversas a las realizadas
#
y_d1d12_scaled_m2 = data_d1d12_scaled[0:P] + y_d1d12_scaled_m2.tolist()

y_d1d12_m2 = scaler.inverse_transform([[u] for u in y_d1d12_scaled_m2])
y_d1d12_m2 = [u[0] for u in y_d1d12_m2.tolist()]

y_d1_m2 = [y_d1d12_m2[t] + data_d1[t] for t in range(len(y_d1d12_m2))]
y_d1_m2 = data_d1[0:12] + y_d1_m2

y_m2 = [y_d1_m2[t] + data[t] for t in range(len(y_d1_m2))]

y_m2 = [data[0]] + y_m2

[22]:

plt.figure(figsize=(14, 5))
plt.plot(data, ".-k")
plt.grid()

plt.plot(y_m2, "-r")

# linea vertical para dividir el entrenamiento
# del pronóstico. Se ubica en el ultimo dato
# usando para entrenamiento
plt.plot([len(data) - 24, len(data) - 24], [min(data), max(data)], "--", linewidth=2);

../_images/43_modelos_de_redes_neuronales_08_pronostico_series_de_tiempo_42_0.png