1edam6np06zloswt3kbc Ra.png

Cómo hacer un análisis de datos exploratorio de una serie temporal

Imagen por autor.

Uno de los tipos de datos más populares son las series temporales. Vídeos, imágenes, píxeles, señales, literalmente cualquier cosa que tenga un componente de tiempo podría convertirse en eso. Formalmente, una serie de tiempo es una secuencia de histórico mediciones de una variable observable en intervalos de tiempo iguales.

En este artículo quiero sugerir un pequeño canal que cualquiera puede utilizar al analizar una serie temporal. Podría ayudarle a obtener información significativa sobre los datos en sí, prepararlos para el modelado y sacar algunas conclusiones preliminares.

Como mi palabra favorita es geoespacial 🌏, hoy analizaremos una serie temporal meteorológica. En particular, exploraremos la temperatura del aire a 2 m, la precipitación total, la radiación solar neta en la superficie y la presión en la superficie en un punto en el sudeste de Siberia durante 2023, derivadas de horas Tierra ERA5 [1] Reanálisis climático.

Como siempre, para seguir el tutorial, puedes descargar y ejecutar el cuaderno. aquí.

Para realizar el análisis necesitamos importar varias bibliotecas:

import pandas as pd
import seaborn as sns
import numpy as np

import matplotlib.pyplot as plt
import xarray as xr

import statsmodels.api as sm
from statsmodels.graphics.tsaplots import plot_acf, plot_pacf
from scipy import stats

También decidí probar un nuevo estilo matplotlib de dos bibliotecas, dogmático y ambivalente:

from ambivalent import STYLES
import opinionated
plt.style.use(STYLES['ambivalent'])
plt.style.use("dark_background")

Primero que nada, carguemos y visualicemos los datos que tenemos. Para manejar matrices multidimensionales geoespaciales usaremos matriz x biblioteca.

data = xr.open_dataset('Medium_data.nc')
data
Imagen por autor.

Ahora necesitamos dividir los datos específicamente para la ubicación elegida y convertirlos a un marco de datos de pandas y crear un diagrama de líneas:

df = data.sel(latitude=52.53, longitude=101.63, method='pad').to_pandas().drop(['latitude', 'longitude'], axis=1)
fig, ax = plt.subplots(ncols = 2, nrows = 2, figsize=(16,9))
df['t2m'].plot(ax=ax[0,0])
ax[0,0].set_title('Air Temperature')
df['ssr'].plot(ax=ax[0,1])
ax[0,1].set_title('Surface Net Solar Radiation')
df['sp'].plot(ax=ax[1,0])
ax[1,0].set_title('Surface Pressure')
df['tp'].plot(ax=ax[1,1])
ax[1,1].set_title('Total Precipitation')
plt.tight_layout()
plt.show()
Imagen por autor.

Ya está claro en los gráficos lineales que las cuatro series temporales tienen características diferentes, así que investigémoslas utilizando herramientas matemáticas.

Cualquier serie temporal tiene tres atributos importantes a considerar:

  • Tendenciaque es un cambio fluido a largo plazo en una serie temporal;
  • Estacionalidadque se refiere a una serie de tiempo con un cambio periódico regular en la media de la serie;
  • Ruido (residuales), cual es componente aleatorio de una señal con media igual a cero.

Para obtener cada uno de estos componentes por separado, es posible producir descomposición clásica (aditivo o multiplicativo). Esta operación se produce aplicando un filtro convolucional, por lo que cada componente de la serie temporal se define como

Descomposición aditiva. Imagen por autor.

o

Descomposición multiplicativa. Imagen por autor.

dónde y – es un valor de una serie de tiempo, S — componentes estacionales, t — componente de tendencia y norte – ruido.

Para producir la descomposición, además de seleccionar el tipo de descomposición, es necesario establecer un período estacional (por ejemplo, p=1 para datos anuales, p=4 para datos trimestrales, p=12 para datos mensuales, etc.).

Es importante mencionar que la descomposición clásica antes mencionada es un método realmente ingenuo y simple. Tiene limitaciones importantes, como su linealidad, la incapacidad de capturar la estacionalidad dinámica y la dificultad para manejar la no estacionariedad en una serie temporal. Sin embargo, a los efectos de este artículo, este método es más que suficiente.

Para hacer la descomposición clásica usaremos descomposición_estacional función de la biblioteca statsmodels con un período igual a 24, ya que estamos tratando con datos horarios:

vars = {'t2m': 'Air Temperature', 'tp': 'Total Precipitation', 'sp': 'Surface Pressure', 'ssr': 'Surface Net Solar Radiation'}
for var in df.columns:
result = sm.tsa.seasonal_decompose(df[var], model='additive', period = 24)
results_df = pd.DataFrame({'trend': result.trend, 'seasonal': result.seasonal, 'resid': result.resid, 'observed': result.observed})
fig, ax = plt.subplots(ncols = 2, nrows = 2,figsize=(16,9))
ax[0,0].plot(df.index, results_df.trend)
ax[0,0].set_title('Trend')
ax[0,0].set_ylabel('Value')

ax[0,1].plot(df.index, results_df.seasonal)
ax[0,1].set_title('Seasonal')

ax[1,0].plot(df.index, results_df.resid)
ax[1,0].set_title('Residual')
ax[1,0].set_ylabel('Value')
ax[1,0].set_xlabel('time')

ax[1,1].plot(df.index, results_df.observed)
ax[1,1].set_title('Observed')
ax[1,1].set_xlabel('time')

opinionated.set_title_and_suptitle(vars[var], f"Dickey-Fuller test: {round(sm.tsa.stattools.adfuller(df[var])[1],5)}", position_title=[0.45,1],
position_sub_title=[0.95, 1])
plt.tight_layout()
plt.savefig(f'Seasonal_{var}.png')
plt.show()

Imagen por autor.
Imagen por autor.
Imagen por autor.
Imagen por autor.

Puedes ver que para todas las variables el componente estacional parece un desastre. Dado que analizamos datos horarios, estas variaciones estacionales se observan dentro de un día y no son realmente informativas. En este caso, vale la pena intentar volver a muestrear los datos a resolución diaria y realizar la descomposición durante el período de un día.

df_d = df.resample('1d').mean()
Imagen por autor.

A estas alturas, algunos de ustedes probablemente habrán notado la etiqueta de prueba de Dickey-Fuller en la esquina superior derecha del gráfico. Esto es un estacionariedad prueba, que se realizó utilizando el adfuller función de la misma biblioteca. En el caso de una serie temporal, la estacionariedad significa que las propiedades de una serie temporal no cambian con el tiempo. Al decir propiedades me refiero a varianza, estacionalidad, tendencia y autocorrelación.

Al aplicar la prueba aumentada de Dickey-Fuller (ADF) a una serie de tiempo, planteamos la hipótesis nula de que la serie de tiempo no es estacionaria. Luego seleccionamos un nivel de significancia α, que suele ser del 5%. En esencia, α es la probabilidad de rechazar incorrectamente la hipótesis nula cuando en realidad es cierta. Entonces, en nuestro caso, α=5%, hay un 5% de riesgo de concluir que la serie temporal es estacionaria, cuando en realidad no lo es.

El resultado de la prueba nos dará un valor p. Si es inferior a 0,05, podemos rechazar nuestra hipótesis nula.

Como puede ver, las 4 variables son estacionarias según la prueba ADF.

Por lo general, para aplicar algunos modelos de pronóstico de series temporales como ARIMA y otros, la estacionariedad es imprescindible, por lo que tenemos suerte aquí. Generalmente, los datos meteorológicos y climáticos suelen analizarse en diferentes materiales de aprendizaje relacionados con series temporales, ya que en la mayoría de los casos son estacionarios.

Después de concluir que todas nuestras series temporales son estacionarias, echemos un vistazo a cómo se distribuyen. Para ello utilizaremos la conocida biblioteca seaborn y su función. diagrama de paresque permite crear gráficos informativos con hists y kdes.

ax = sns.pairplot(df, diag_kind='kde')
ax.map_upper(sns.histplot, bins=20)
ax.map_lower(sns.kdeplot, levels=5, color='.1')
plt.show()
Imagen por autor.

Consideremos el ejemplo de t2m (1 fila y 1 columna). Al analizar el gráfico de estimación de densidad del kernel (kde), queda claro que la distribución de esta variable es multimodal, es decir, tiene 2 o más “campanas”. Entonces, durante las siguientes etapas del artículo actual intentaremos transformar la variable para que se parezca a una distribución normal.

Otros gráficos en la primera columna y fila son idénticos en términos de la información que brindan, pero se visualizan de manera diferente. Básicamente, se trata de diagramas de dispersión, que permiten identificar cómo se correlacionan dos variables. Entonces, cuanto más oscuro sea el color de un punto o cuanto más cerca esté del círculo central, mayor será la densidad de puntos en esta área.

Como hemos descubierto que la serie temporal de la temperatura del aire es estacionaria, pero no está distribuida normalmente, intentemos solucionarlo mediante la transformación de Box-Cox. Para hacer eso usaremos el paquete scipy y su función. boxcox.

df_d['t2m_box'], _ = stats.boxcox(df_d.t2m)
fig, ax = plt.subplots(nrows=1, ncols=2, figsize=(15,7))
sns.histplot(df_d.t2m_box, kde=True, ax=ax[0])
sns.histplot(df_d.t2m, kde=True, ax=ax[1])
Imagen por autor.

La parte izquierda de la gráfica es nuestra distribución de series temporales después de la transformación BoxCox y, como puede ver, todavía está lejos de ser llamada distribución «normal». Pero si lo comparamos con el de la derecha, podemos decir que definitivamente se acercó más a uno normal.

Otra cosa que podemos hacer para asegurarnos de que la transformación realizada fue útil es crear un gráfico de probabilidad. En esencia, trazamos cuantiles de una distribución teórica (en nuestro caso normal) contra muestras de nuestros datos empíricos (es decir, las series de tiempo que consideramos). Cuanto más cerca estén los puntos de la línea blanca, mejor.

fig = plt.figure()

ax1 = fig.add_subplot(211)
prob = stats.probplot(df_d.t2m, dist=stats.norm, plot=ax1)
ax1.get_lines()[1].set_color('w')
ax1.get_lines()[0].set_color('#8dd3c7')
ax1.set_title('Probplot against normal distribution')

ax2 = fig.add_subplot(212)
prob = stats.probplot(df_d.t2m_box, dist=stats.norm, plot=ax2)
ax2.get_lines()[1].set_color('w')
ax2.get_lines()[0].set_color('#8dd3c7')
ax2.set_title('Probplot after Box-Cox transformation')
plt.tight_layout()fig = plt.figure()

ax1 = fig.add_subplot(211)
prob = stats.probplot(df_d.t2m, dist=stats.norm, plot=ax1)
ax1.set_title('Probplot against normal distribution')

ax2 = fig.add_subplot(212)
prob = stats.probplot(df_d.t2m_box, dist=stats.norm, plot=ax2)
ax2.set_title('Probplot after Box-Cox transformation')
plt.tight_layout()

Imagen por autor.

Si va a utilizar la serie temporal transformada para el modelado de ML, no olvide aplicar la transformación BoxCox inversa; de lo contrario, tendrá que lidiar con números inadecuados.

Y el último paso de nuestro análisis es la autocorrelación. La función de autocorrelación (ACF) estima la correlación entre una serie temporal y una versión retrasada de la misma. O en otras palabras, cómo se correlaciona un valor específico de una serie temporal con otros valores anteriores en diferentes intervalos de tiempo.
También podría resultar útil trazar la función de autocorrelación parcial (PACF), que es lo mismo que la autocorrelación, pero sin la correlación en retrasos más cortos. Por tanto, estima la correlación entre valores dentro de una determinada marca de tiempo, pero controlando la influencia de otros valores.

for var in df.columns[:-1]:
fig, (ax1, ax2) = plt.subplots(2,1,figsize=(10,8))
plot_acf(df_d.t2m, ax = ax1)
plot_pacf(df_d.t2m, ax = ax2)
opinionated.set_title_and_suptitle(vars[var], '',position_title=[0.38,1],
position_sub_title=[0.95, 1])
plt.tight_layout()
plt.show()
Imagen por autor.

Como puede ver, existe una autocorrelación parcial realmente fuerte en la serie temporal de presión superficial con un retraso de 1 día. Luego se debilita significativamente y después de un retraso de 3 días casi desaparece. Un análisis de este tipo podría ayudarle a comprender mejor la naturaleza de los datos con los que está tratando y, por tanto, a extraer conclusiones más significativas.