Sesión 5: Introducción a las series de tiempo

Curso: Análisis de datos

Phd (c) Melanie Oyarzún - melanie.oyarzun@udd.cl

Magíster en Data Science - Universidad del Desarrollo

# Paquetes y settings

from dateutil.parser import parse 
import matplotlib as mpl
import matplotlib.pyplot as plt
import seaborn as sns
import numpy as np
import pandas as pd

# setting de graficos
plt.figure(figsize=(5,3), dpi= 200, facecolor='w', edgecolor='k')

<Figure size 1000x600 with 0 Axes>

<Figure size 1000x600 with 0 Axes>

Overview

Resultado de aprendizaje esperado:

Identificar datos de series temporales, sus particularidades y riesgos, en el contexto de posibles aplicaciones profesionales.

Bibliografía recomendada:

Stock & Watson, C.14 link ; Wooldridge, c.12 link, Gujarati, c.12 link

Tipos de datos

Recordemos

Los datos que analizamos con modelos se pueden categorizar en tres grandes tipos:

Corte transversal
Series de tiempo
Panel

Corte transversal

La unidad de observación son individuos separados unos de otros.
Cada fila reprsenta un individuo (personas, familias, empresas) en un momento del tiempo específico (mismo año, mismo mes, etc). (usualmente)
Es análogo a una fotografía de un grupo.

Un ejemplo de esto es la encuesta Casen (http://observatorio.ministeriodesarrollosocial.gob.cl/encuesta-casen).

import pandas as pd

df_casen2020= pd.read_stata("data_sesion5/casen_2020_ingresos.dta")
df_casen2020.head(5)

	folio	o	id_persona	region	comuna	zona	expr	edad	sexo	tot_per	...	esc2	educ	o1	yaut	yauth	yautcor	yautcorh	ytrabajocor	ytrabajocorh	yae
0	1.101100e+11	1	5	Región de Tarapacá	Iquique	Urbano	67	34	Mujer	2	...	12.0	Media humanista completa	No	220000.0	300000	220000.0	300000	150000.0	150000.0	240586.0
1	1.101100e+11	2	6	Región de Tarapacá	Iquique	Urbano	67	4	Mujer	2	...	NaN	Sin educación formal	NaN	80000.0	300000	80000.0	300000	NaN	150000.0	240586.0
2	1.101100e+11	2	31	Región de Tarapacá	Iquique	Urbano	67	5	Mujer	3	...	NaN	Básica incompleta	NaN	25000.0	941583	25000.0	941583	NaN	891583.0	439170.0
3	1.101100e+11	1	32	Región de Tarapacá	Iquique	Urbano	67	45	Hombre	3	...	15.0	Técnico nivel superior incompleta	Sí	889500.0	941583	889500.0	941583	889500.0	891583.0	439170.0
4	1.101100e+11	3	30	Región de Tarapacá	Iquique	Urbano	67	19	Mujer	3	...	NaN	No sabe	No	27083.0	941583	27083.0	941583	2083.0	891583.0	439170.0

5 rows × 22 columns

df_casen2020.info()

<class 'pandas.core.frame.DataFrame'>
Int64Index: 185437 entries, 0 to 185436
Data columns (total 22 columns):
 #   Column        Non-Null Count   Dtype   
---  ------        --------------   -----   
 0   folio         185437 non-null  float64 
 1   o             185437 non-null  int8    
 2   id_persona    185437 non-null  int16   
 3   region        185437 non-null  category
 4   comuna        185437 non-null  category
 5   zona          185437 non-null  category
 6   expr          185437 non-null  int16   
 7   edad          185437 non-null  int16   
 8   sexo          185437 non-null  category
 9   tot_per       185437 non-null  int8    
 10  ecivil        153892 non-null  category
 11  esc           148886 non-null  float64 
 12  esc2          148886 non-null  float64 
 13  educ          185437 non-null  category
 14  o1            151315 non-null  category
 15  yaut          95399 non-null   float64 
 16  yauth         185437 non-null  int32   
 17  yautcor       102165 non-null  float64 
 18  yautcorh      185437 non-null  int32   
 19  ytrabajocor   73509 non-null   float32 
 20  ytrabajocorh  185437 non-null  float32 
 21  yae           185339 non-null  float64 
dtypes: category(7), float32(2), float64(6), int16(3), int32(2), int8(2)
memory usage: 15.6 MB

import seaborn as sns
# Crea el scatterplot
sns.scatterplot(data=df_casen2020, x="esc", y="yaut", alpha=0.5)

<Axes: xlabel='esc', ylabel='yaut'>

Series de tiempo

Se tienen diferentes momentos del tiempo (día, semana, mes, año, etc.) para una misma unidad de análisis
- individuo, país, empresa, etc

Serie Temporal

Ejemplos típicos de series de tiempo son:
- Datos macroeconómicos (PIB. Inglación, empleo, etc.)
- Financieros (precios de acciones)
- Empresariales (Ventas, costos, etc…)

Ejemplo datos de acciones

Cargar datos
Plot

#Usemos la api de Yahoo finance: instalar: pip install yahoo_fin y pip install requests_html

from yahoo_fin.stock_info import get_data

amazon_weekly= get_data("amzn", start_date="12/04/2009", end_date="25/09/2024", index_as_date = False, interval="1wk")

amazon_weekly.head()

	date	open	high	low	close	adjclose	volume	ticker
0	2009-11-30	7.1810	7.2955	6.7555	6.8790	6.8790	627018000	AMZN
1	2009-12-07	6.9000	6.9500	6.4910	6.7075	6.7075	957260000	AMZN
2	2009-12-14	6.6250	6.6305	6.2825	6.4240	6.4240	915898000	AMZN
3	2009-12-21	6.5240	6.9850	6.5095	6.9235	6.9235	648120000	AMZN
4	2009-12-28	6.9875	7.1290	6.7260	6.7260	6.7260	572014000	AMZN

sns.scatterplot(data=amazon_weekly, y="open", x="date")

<Axes: xlabel='date', ylabel='open'>

Uso series de tiempo

La información temporal permite responder cómo una variable (o más) responde a cambios a través del tiempo.

¿Cuál es el efecto causal dinámico de \(X_t\) sobre \(Y_t\)?
¿Cuál es la mejor predicción del valor de Y en el futuro?

PERO su uso sin cuidado puede llevarnos a conclusiones MUY equivocadas

Panel

Este es la combinacion de los dos tipos anteriores
Para un conjunto de individuos tenemos varias observaciones en el tiempo.

Panel

Suele ser el tipo de datos más completo,
- pero también más dificil de conseguir.
Además, enfrenta los problemas típicos de ambos tipos anteriores, dependiendo si es un panel corto o largo.
Una alternativa a estos son los llamados pooled cross section.´
El análisis de este tipo de datos escapa a los objetivos del taller, pero una introducción pueden revisarla en Stock y Watson Cap. 10.

VERSUS

Cortes transversal vs series de tiempo

Existe una diferencia clave, que se desprende de trabajar con corte tranvseral vs series de tiempo y es el concepto de la muestra aleatoria.

Muestras en corte transversal

En corte transversal solemos trabajar con MUESTRAS

Estimamos modelos de la forma:

\[ y_{i}=\beta_{0}+\beta_{1}x_{1i}+\dots+\beta_{1}x_{ki}+u_{i} \]

Muestras en corte transversal

Estimamos modelos de la forma:

\[ y_{i}=\beta_{0}+\beta_{1}x_{1i}+\dots+\beta_{k}x_{ki}+u_{i} <-> Y= \beta_{0}+\beta_{1}X_{1}+\dots+\beta_{k}X_{k}+U \]

Estos modelos representan una correlación marginal en las observaciones entre y y x’s (escalada por la varianza de x)
y para que podamos interpretarlas causalmente tenemos varias condiciones o supuestos que de se deben cumplir.
Uno de estos, es el supuesto de exogeneidad: \[ E[u_{i}|X_{i}]=0 \]

Muestras en corte transversal

Sin embargo, esta forma de ver este supuesto es una simplificación
Como estamos en una muestra aleatoria, no tenemos que verificar que los efectos cruzados tambien sean exogenos:

\[ E[u_{i}|X_{j}]=0 \]

Esto se daba por cumplido, como consecuencia de que era una muestra aletaoria.
Lo cual, generalente, ocurre en corte transversal por lo cual cada observación es i.i.d.

Series de tiempo y procesos estocásticos

En serie de tiempo, nuestro universo es un proceso estocástico
En cada momento se observa un posible resultado (o realización) del proceso estocásticos.
Estimamos modelos de la forma: \[ y_{i}=\beta_{0}+\beta_{1}x_{1i}+\dots+\beta_{1}x_{ki}+u_{i} \]

La imposibilidad de la muestra aleatoria en series de tiempo

Esto tiene vaias implicancias:

NO SON INDEPENDIENTES ya que por construcción, viene del mismo proceso.
Por lo tanto, el supuesto de exogeneidad \(E[u_{i}|X_{i}]=0\) NO ES SUFICIENTE
Requerimos su versión más exigente:
(Exogeneidad estricta): \[E[u_{t}|X_{s}]=0 \qquad \forall s \]
Este supuesto, generalmente NO SE CUMPLE.
Si se cumpliese, seguiruíamos operando como siempre con modelos de regresión multiple estándar.

Trabajando con series de tiempo:

¿Qué hacer entonces?

Reconocer los datos como procesos estocásticos e incluir sus particularidades en la modelación. De eso se tratarán nuestras dos sesiones
- Peculiaridades de los procesos estocasticos y su exploración en la data (Sesion 5)
  - Dependencia temporal, inercia, estacionaerieada y estacionalidad.
- Modelos específicos para series de tiempo (Sesión 6)
  - ARIMA, SARIMAX, GARCH, entre otros

El lenguaje de las series de tiempo

Notación

Variables de series de tiempo se denominan con sub-índice “t” para indicar el perído en el tiempo: \(y_t\)
El total de periodos se suele referir como T.

Transformaciones temporales

Revisemos alginas peculiaridades de las series de tiempo

Dado que una serie de tiempo tiene un orden específico, que en si mismo es importante.
En base a este orden se suelen crear nuevos indicadores o variables.
- Revisemos los más comunes:
  - Rezagos
  - Diferencias
  - Tazas de crecimiento

Rezagos

Un rezago es el valor de la variable en períodos anteriores:
- Primer rezago: \(y_{t-1}\) es el valor 1 período anterior
- Segundo rezago: \(y_{t-2}\) es el valor 2 períodos atrás
- j-ésimo rezago: \(y_{t-j}\) es el valor j períodos atrás

Ejemplo rezago

Calculemos el primer y segundo rezago en los datos de amazon:

# Calcular el primer rezago (Lag 1) para la columna 'open'
amazon_weekly['open_lag1'] = amazon_weekly['open'].shift(1)

# Calcular el segundo rezago (Lag 2) para la columna 'open'
amazon_weekly['open_lag2'] = amazon_weekly['open'].shift(2)

# Ver los primeros registros del DataFrame con las columnas de rezago
amazon_weekly[['date', 'open', 'open_lag1', 'open_lag2']].head()

	date	open	open_lag1	open_lag2
0	2009-11-30	7.1810	NaN	NaN
1	2009-12-07	6.9000	7.181	NaN
2	2009-12-14	6.6250	6.900	7.181
3	2009-12-21	6.5240	6.625	6.900
4	2009-12-28	6.9875	6.524	6.625

Diferencias

Una diferencia corresponde al cambio en una variable entre dos periodos específicos y se usa la notación \(\Delta\)

Primera diferencia: \[\Delta y_t = y_t- y_{t-1}\]
Ejemplo: Calculemos la primera y segunda diferencia en los datos de Amazon:

# Calcular la primera diferencia para la columna 'open'
amazon_weekly['open_diff1'] = amazon_weekly['open'] - amazon_weekly['open'].shift(1)

# Calcular la segunda diferencia para la columna 'open'
amazon_weekly['open_diff2'] = amazon_weekly['open_diff1'] - amazon_weekly['open_diff1'].shift(1)

# Ver los primeros registros del DataFrame con las columnas de diferencia
print(amazon_weekly[['date', 'open', 'open_diff1', 'open_diff2']].head())

        date    open  open_diff1  open_diff2
0 2009-11-30  7.1810         NaN         NaN
1 2009-12-07  6.9000     -0.2810         NaN
2 2009-12-14  6.6250     -0.2750      0.0060
3 2009-12-21  6.5240     -0.1010      0.1740
4 2009-12-28  6.9875      0.4635      0.5645

Tasas de crecimiento

Si calculamos la primera diferencia el logaritmo natural, podemos obtener incorporar la tasa de crecimiento en una regresión:

Primera diferencia en logs: \[ \Delta ln(y_{t})=ln(y_{t})-ln(y_{t-1}) \]
Cambio porcentual de \(y_t\) entre \(t -1\) y \(t\approx 100\times \Delta ln(y_{t})\)

Ejemplo Tasa de crecimiento

Calculemos la primera diferencia en logaritmo natural y la taza de crecimiento para la serie de Amazon, en su precio de apertura:

import numpy as np

# Calcular el logaritmo natural de la columna 'open' con NumPy
amazon_weekly['open_ln'] = np.log(amazon_weekly['open'])

# Calcular la primera diferencia en los logaritmos naturales ('open_ln_diff1')
amazon_weekly['open_ln_diff1'] = amazon_weekly['open_ln'] - amazon_weekly['open_ln'].shift(1)

# Calcular el cambio porcentual entre t-1 y t para 'open' ('open_percent_change')
amazon_weekly['open_percent_change'] = (amazon_weekly['open'] - amazon_weekly['open'].shift(1)) / amazon_weekly['open'].shift(1) * 100

# Ver los primeros registros del DataFrame con las columnas calculadas
print(amazon_weekly[['date', 'open', 'open_ln', 'open_ln_diff1', 'open_percent_change']].head())

        date    open   open_ln  open_ln_diff1  open_percent_change
0 2009-11-30  7.1810  1.971439            NaN                  NaN
1 2009-12-07  6.9000  1.931521      -0.039917            -3.913106
2 2009-12-14  6.6250  1.890850      -0.040671            -3.985509
3 2009-12-21  6.5240  1.875488      -0.015363            -1.524526
4 2009-12-28  6.9875  1.944123       0.068635             7.104537

Tipos de modelos:

Modelo estático

Se modela la relación contemporánea entre dos variables, i.e., relación en el mismo momento en el tiempo

\[y_{t}=\beta_{0}+\beta_{1}z_{t}+u_{t}\]

Ejemplos:

Modelo Simple: \[ inflación_{t}=\beta_{0}+\beta_{1}desempleo_{t}+u_{t}\]
Modelo múltiple: \[ homicidio_{t}=\beta_{0}+\beta_{1}condena_{t}+\beta_{2}desmepleo_{t}+\beta_{3}hombres_{t}+u_{t} \]

Tipos de modelos:

Modelo dinámico

Se incluyen efectos temporales, que se piensan que tienen que ver con tendencias, inercia o estacionalidades

Un modelo con inercia en la inflación (por conceptos): \[ inflación_{t}=\beta_{0}+\beta_{1}desempleo_{t}+\beta_{2} inflación_{t-1} +u_{t}\]
Hay modelos conocidos:
- Ejemplo: ARIMA de primer órden \[ \Delta inflacion_t = \phi \cdot \Delta inflacion_{t-1} + \varepsilon_t\]

Supuestos no cumplidos y opciones de modelamiento

Para poder interpretar causalmente un modelo de regresión lineal multiple en series de tiempo necesitamos que se cumplan los siguientes supuestos de Gauss Markov:

Modelo lineal en parámetros \[y_{t}=\beta_{0}+\beta_{1}x_{1t}+\dots+\beta_{1}x_{kt}+u_{t}\]
Media condicionada nula o exogeneidad (que ahora es estricta) \[ E[u_{t}|X_{js}]=0 \quad \forall s\]
No hay colinealidad perfecta.
Homoscedasticidad \[ Var[u_{t}|X_{js}]=\sigma \quad \forall s \]
No hay correlación serial. \[ Cov[u_{t}, u_{s}|X_{js}]=\sigma \quad \forall s\ne t\]

Supuestos no cumplidos y opciones

Cada uno de estos supuestos tiene diferentes implicancias:
1-3 es que el estimador de MCO \(\hat{\beta}\) es insesgado.
4-5 que tiene minima varianza de los estimadores lineales (eficiente)
1-5 -> MELI
Vemos que hay dos elementos diferentes a corte transversal:
- cambia la exogeneidad a uno más estricto (ya que no hay muestreo aleatorio)
- y aparece un nuevo supuesto, correlación serial.

Objetivos diferentes de los modelos

Los modelos de regresión se pueden usar con dos objetivos principales: predicción y explicación.
En series de tiempo, como no hay muestreo aleatorio, nos enfocaremos en la predicciónn.
- Un modelo de regresión puede ser útil para la predicción, aún cuando ninguno de sus coeficientes tenga interpretación causal.
- Desde el punto de vista de la predicción, lo que es importante es que el modelo entregue una predicción lo más precisa posible.
La idea base es aprovechar las peculiaridades de las series temporales: la dependencia temporal y tendencias, para identificar patrones en la data que enriquezcan la predicción aun cuando no tengan valor explicativo. Estos elementos los incluiremos DENTRO de la regresión, para que represente de mejor manera el fenómeno a predecir.

Modelos explicativos en series temporales:

Los modelos que se centran en la explicación: Se basen en modelos teoricos estructurales que se corroboran en la data
Se usan otras aproximaciones empíricas como Causalidad de Granger
- un concepto estadístico que evalúa si una serie temporal puede predecir otra serie temporal.
- Ayuda a determinar si un conjunto de datos precedentes es útil para predecir un conjunto de datos posterior.
Un Concepto relacionado es la cointegración:
- La cointegración identifica relaciones estadísticas a largo plazo entre series temporales, lo que puede ser relevante al evaluar la causalidad de Granger.
- A pesar de tendencias o patrones a largo plazo, existe una combinación lineal que es estacionaria.

Patrones de dependencia intertemporal

Modelando series de tiempo desde la intuiciónn

¿Cuándo falla supuesto de exogeneidad estricta?
Cuando tenemos particularidades temporales actuando en el proceso estocástico:
- Efecto rezagados (variable independiente rezagada)
- Retroalimentación entre variables
- Auto-dependencia, la variable depende de si misma en el pasado (variable dependiente rezagada, AR(1) )
Todas estas pueden evitar que nuestras series sean estacionarias

Estacionareidad

Entenderemos estacionareidad si el futuro se parece al pasado, al menos en un sentido probabilistico.
Una serie es estacionaria si:

Su distribución de probabilidad no varia en el tiempo,
Su media, varianza y covarianza es constante

Podemos reemplazar el supuesto de exogeneidad extricta por exogeneidad debil SI la serie es estacionaria.

Ejemplo estacionaereidad

Hay muchos ejemplos de tipos de no estacionareidad.
Ilustremos los más tipicos con unos ejemplos ficticios.

x = np.linspace(0, np.pi*10, 360)
y = np.sin(x)

Podemos generar los diferentes tipos con algunas manipulaciones algebráicas.

import matplotlib as plt
import matplotlib.pyplot as plt

fig, axs = plt.subplots(2, 2, figsize=(18, 18))
axs[0][0].plot(x, y)
axs[0][0].set_title('Serie Estacionaria')
axs[0][0].set_xlabel('tiempo')
axs[0][0].set_ylabel('Amplitud')

axs[0][1].plot(x, y+x/10)
axs[0][1].set_title('Media cambiante')
axs[0][1].set_xlabel('time')
axs[0][1].set_ylabel('Amplitud')


axs[1][0].plot(x, y*x/10)
axs[1][0].set_title('Varianza cambiante')
axs[1][0].set_xlabel('time')
axs[1][0].set_ylabel('Amplitud')

axs[1][1].plot(np.sin(x+x*x/30))
axs[1][1].set_title('Co-variance cambiante')
axs[1][1].set_xlabel('time')
axs[1][1].set_ylabel('Amplitud')

plt.tight_layout()

Ausencia de estacionaeridad:

Revisaremos los casos más típicos en los cuales no hay estacionaeridad:

Tendencias
Estacionalidad
Quiebre estructural

Ausencia de estacionareidad 1: Tendencias

Un claro ejemplo de series no estacionarias es cuando hay tendencias.
Podemos identificar la tendencia a través de una media movil
Revisemoslo con los datos de pasajeros de aerolineas.

import pandas as pd

# Cargar datos

airline = pd.read_csv('data_sesion5/international-airline-passengers.csv', sep=';')

# Nota: si no les reconoce bien la dependencia de la carpeta, pueden usar también
#airline= pd.read_csv("https://github.com/melanieoyarzun/taller_seriestiempo_IDS/blob/8c0b9774be8d4103da3801d3069d82b4fe006461/Data/international-airline-passengers.csv?raw=true", sep=';')

airline['Month'] = pd.to_datetime(airline['Month']+'-01')
airline.set_index('Month', inplace=True)

airline.head()

	Passengers
Month
1949-01-01	112
1949-02-01	118
1949-03-01	132
1949-04-01	129
1949-05-01	121

Ausencia de estacionareidad 1: Tendencias

Ejemplo Aerolinea

Con un grafico rápido, identificamos que está todo bien y que efectivamente se observa que la serie no es estacionaria.

. . .

import matplotlib as mpl
import matplotlib.pyplot as plt

fig, ax = plt.subplots(1, 1)
ax.plot(airline.index, airline['Passengers'])
ax.set_xlabel('Año')
ax.set_ylabel('Pasajeros')

Text(0, 0.5, 'Pasajeros')

Ausencia de estacionareidad 1: Tendencias

Ejemplo Aerolinea

Podemos identificar la tendencia en los datos, al calcular la media movil.
Definimos una función para calcular la media móvil:

. . .

def running_average(x, order):
    current = x[:order].sum()
    running = []
    
    for i in range(order, x.shape[0]):
        current += x[i]
        current -= x[i-order]
        running.append(current/order)
    
    return np.array(running)

Esta función es autoexplicativa.
Simplemente corre en el dataset paso a paso y calcula la media en una ventana específica.

Ausencia de estacionareidad 1: Tendencias

Ejemplo Aerolinea

Ahora podemos agregar esta línea de tendencia al gráfico anterior.

import numpy as np

trend = running_average(airline['Passengers'], 12)

fig, ax = plt.subplots(1, 1)
ax.plot(airline.index, airline['Passengers'])
ax.set_xlabel('Año')
ax.set_ylabel('Pasajeros')
ax.plot(airline.index[12:], trend, label='Tendencia')
ax.legend()

Ausencia de estacionareidad 1: Tendencias

Ejemplo Aerolinea

A la serie, entonces, le podemos sacar esta tendencia, al dividr.

detrended = airline.iloc[12:].values.flatten()/trend

. . .

Y graficamente:

fig, ax = plt.subplots(1, 1)
ax.plot(airline.index[12:], detrended)
ax.set_xlabel('Date')
ax.set_ylabel('Detrended value')

Text(0, 0.5, 'Detrended value')

Ausencia de estacionareidad 2: Estacionalidad

(no confundir con estacionareidad)

Es cuando hay un patron claro de ciclos que se repite en el tiempo.
Generalmente los identificamos a priori con una inspección del grafico.
Podemos usar una función para identificarlos.

. . .

def plot_seasons(detrended, order, plot_mean = True):
    colors = plt.rcParams['axes.prop_cycle'].by_key()['color']

    N = len(detrended)

    data = np.array([detrended[i::order] for i in range(order)])
    
    means = np.mean(data, axis=1)
    medians = np.median(data, axis=1)
    
    counts = [0]
    counts.extend([len(data[i]) for i in range(order)])
    counts = np.cumsum(counts)

    ticks = (counts[:-1]+counts[1]/2)
    
    for i in range(order):
        values = data[i, :]
        npoints = len(values)

        plt.plot(range(counts[i], counts[i+1]), values, c=colors[0])
        plt.plot(range(counts[i], counts[i+1]), np.ones(npoints)*means[i], c=colors[1])
        plt.plot(range(counts[i], counts[i+1]), np.ones(npoints)*medians[i], c=colors[2])

    plt.legend(['data', 'mean', 'median'])
    plt.xlabel('season')
    plt.ylabel('values')
    plt.xticks(ticks, np.arange(order));
    
    if plot_mean:
        plt.plot(ticks, means, c=colors[3])
    
    return means

Ausencia de estacionareidad 2: Estacionalidad

Acá, simplemente vamos a graficar la curva para diferentes periodos en la temporada.
Esto se hace yendo a traves del dataset con un stride igual al periodo estacional.
La figura tambien nos provee una forma arquetípica de comportamiento estacional.
Con esto, podememos terminar la descomposición de la data en tres componentes.

# Descomposición multiplicativa
def decomposition(data, order, plot=True):
    values = data.values.flatten()
    trend = running_average(values, order)
    detrended = values[order:]/trend
    
    season = [detrended[i::order].mean() for i in range(order)]
    seasonality = np.array(season*(detrended.shape[0]//order+1))[:detrended.shape[0]]
    residuals = values[order:]/(trend*seasonality)

    if plot:
        fig, axs = plt.subplots(4, 1, figsize=(22, 16), sharex=True)
        index = data.index

        axs[0].plot(index, values)
        axs[0].set_title('Data original')
        
        axs[1].plot(index[order:], trend)
        axs[1].set_title('Tendencia')

        axs[2].plot(index[order:], detrended)
        axs[2].set_title('Estacionalidad')

        axs[3].plot(index[order:], residuals)
        axs[3].set_title('Residuos')
        
    return values, trend, seasonality, residuals

Ausencia de estacionareidad 2: Estacionalidad

Con esto , el patron estacional se remueve al repetir el patron medio estacional identificado anteriormente y dividiendo, desde la data sin tendencia.
El resultado de esta division son simpemente los residuos.
Lo que nos muestra que esta descomposición es demasiado sencilla aun.

values, trend, seasonality, residuals = decomposition(airline, 12)

Acá usamos una descomposición multiplicativa, pero este mismo proneso se podría haber hecho siguiendo una descomposición aditiva, con pequeños cambios al codigo.

Ausencia de estacionareidad 3: Quiebre estructural

Otro tipo de no estacionaeridad se presenta cuando la función de regresión poblacional cambia en el transcurso de la las observaciones.
Esto puede ocurrir por varios motivos, por ejemplo cambios en una politica económica, cambios en la estructura de la economía, un nuevo invento o disrupción tecnológica, etc.
Si ocurren tales “cambios estructurales” o “rupturas”, entonces un modelo de regresión que no tenga en cuenta esos cambios puede proporcionar una base engañosa para la inferencia ya la predicción.

Ausencia de estacionareidad 3: Quiebre estructural

Identificación

Las estrategias para identificar un cambio estructurale son varias revisaremos dos:
Contrastes de hipótesis comparando cambios en los coeficientes de regresión mediante estadístico F o Test de Chow.
La segunda es biuscar potenciales cambios estructurales desde la predicción: se simula que la mmuestra termina antes de lo que realmente lo hace y se comparan las predicciones.
- Los cambios estructurales se detectan cuando la capacidad de predicción es sustancialmente peor de lo esperado.

Ausencia de estacionareidad 3: Quiebre estructural

Ejemplo datos de amazon

Podemos utilizar datos de acciones y verificar si hay un quiebre estructural antes y después de la crisis financiera de 2008.
Para hacerlo, necesitaremos datos históricos de acciones y realizaremos análisis de series temporales.

. . .

# Importar las bibliotecas necesarias
import numpy as np
import pandas as pd
import statsmodels.api as sm
import matplotlib.pyplot as plt

# Obtener datos históricos de acciones de Amazon
from yahoo_fin.stock_info import get_data
amazon_subprime= get_data("amzn", start_date="12/04/2000", end_date="25/09/2015", index_as_date = False, interval="1wk")
amazon_subprime.head()

	date	open	high	low	close	adjclose	volume	ticker
0	2000-12-04	1.259375	1.381250	1.006250	1.171875	1.171875	1013370000	AMZN
1	2000-12-11	1.143750	1.375000	1.087500	1.143750	1.143750	804546000	AMZN
2	2000-12-18	1.037500	1.059375	0.743750	0.778125	0.778125	1390856000	AMZN
3	2000-12-25	0.815625	0.925000	0.750000	0.778125	0.778125	678338000	AMZN
4	2001-01-01	0.790625	0.893750	0.678125	0.728125	0.728125	866064000	AMZN

Ausencia de estacionareidad 3: Quiebre estructural

Ejemplo datos de amazon

Graficamos los datos

import seaborn as sns
# Crea el scatterplot
sns.scatterplot(data=amazon_subprime, y="open", x="date")
plt.xlabel("Fecha")
plt.ylabel("Precio de Apertura")
plt.title("Precios de Apertura de Acciones de Amazon (2000-2015)")
plt.show()

Ausencia de estacionareidad 3: Quiebre estructural

Ejemplo datos de amazon

Podriamos ajustar un modelo a todos los datos o separar anres y después de la crisis subprime.

. . .

# Dividir los datos en tres conjuntos: antes de 2008, después de 2008 y todos los datos
data_before_2008 = amazon_subprime[amazon_subprime['date'] < '2008-01-01']
data_after_2008 = amazon_subprime[amazon_subprime['date'] >= '2008-01-01']

# Ajustar modelos de regresión lineal para cada conjunto de datos
model_all_data = sm.OLS(amazon_subprime['open'], sm.add_constant(np.arange(len(amazon_subprime)))).fit()
model_before_2008 = sm.OLS(data_before_2008['open'], sm.add_constant(np.arange(len(data_before_2008)))).fit()
model_after_2008 = sm.OLS(data_after_2008['open'], sm.add_constant(np.arange(len(data_after_2008)))).fit()

# Graficar los datos y las líneas de regresión
plt.figure(figsize=(12, 6))

# Datos de todos los datos (en negro)
plt.scatter(amazon_subprime['date'], amazon_subprime['open'], color='lightgray', label='Todos los datos')

# Línea de regresión para todos los datos (en negro)
plt.plot(amazon_subprime['date'], model_all_data.predict(sm.add_constant(np.arange(len(amazon_subprime)))), color='black', linewidth=2)

# Datos hasta 2008 (en azul)
plt.scatter(data_before_2008['date'], data_before_2008['open'], color='lightblue', label='Hasta 2008')

# Línea de regresión hasta 2008 (en azul)
plt.plot(data_before_2008['date'], model_before_2008.predict(sm.add_constant(np.arange(len(data_before_2008)))), color='blue', linewidth=2)

# Datos después de 2008 (en rojo)
plt.scatter(data_after_2008['date'], data_after_2008['open'], color='lightcoral', label='Después de 2008')

# Línea de regresión después de 2008 (en rojo)
plt.plot(data_after_2008['date'], model_after_2008.predict(sm.add_constant(np.arange(len(data_after_2008)))), color='red', linewidth=2)

plt.xlabel("Fecha")
plt.ylabel("Precio de Apertura")
plt.title("Líneas de Regresión para Precios de Apertura de Acciones de Amazon")
plt.legend()
plt.show()

Ausencia de estacionareidad 3: Quiebre estructural

Ejemplo datos de amazon

-Claramente, los modelos por separados oarecen que explican mejor. - Al parecer hay quiebre estructural. - Revisemos usando el test de Chow:

# Dividir los datos en dos conjuntos: antes de 2008 y después de 2008
data_before_2008 = amazon_subprime[amazon_subprime['date'] < '2008-01-01']
data_after_2008 = amazon_subprime[amazon_subprime['date'] >= '2008-01-01']

# Ajustar modelos de regresión lineal para cada conjunto de datos
model_before_2008 = sm.OLS(data_before_2008['open'], sm.add_constant(np.arange(len(data_before_2008)))).fit()
model_after_2008 = sm.OLS(data_after_2008['open'], sm.add_constant(np.arange(len(data_after_2008)))).fit()

# Calcular SSR para cada modelo
ssr_before_2008 = np.sum(model_before_2008.resid ** 2)
ssr_after_2008 = np.sum(model_after_2008.resid ** 2)

# Combinar todos los datos en un solo modelo
all_data = amazon_subprime['open']
model_all_data = sm.OLS(all_data, sm.add_constant(np.arange(len(amazon_subprime)))).fit()

# Calcular SSR para el modelo completo
ssr_all_data = np.sum(model_all_data.resid ** 2)

# Calcular el estadístico F para el Test de Chow
k = 2  # Número de coeficientes (incluyendo el intercepto)
N = len(amazon_subprime)
f_statistic = ((ssr_all_data - (ssr_before_2008 + ssr_after_2008)) / k) / ((ssr_before_2008 + ssr_after_2008) / (N - 2 * k))

# Calcular el valor p asociado al estadístico F
from scipy.stats import f
p_value = 1 - f.cdf(f_statistic, k, N - 2 * k)

print("Valor p del Test de Chow:", p_value)

Valor p del Test de Chow: 1.1102230246251565e-16

Podemos usar paquetes que tengan el test directamente programado. Por ejemplo chowtest

Ausencia de estacionareidad 3: Quiebre estructural

Ejemplo datos de amazon

Podemos ver los tres modelos en una tabla integrada:

from stargazer.stargazer import Stargazer

# Crear una lista de modelos que deseas incluir en la tabla
modelos = [model_before_2008, model_after_2008, model_all_data]

# Crear una lista de etiquetas para los modelos
etiquetas = ["Modelo Antes de 2008", "Modelo Después de 2008", "Modelo Completo"]

# Configurar Stargazer
stargazer = Stargazer(modelos)
stargazer.custom_columns(etiquetas, [1, 1, 1])  # Asignar etiquetas a las columnas

# Generar la tabla HTML o LaTeX 
#tabla = stargazer.render_html()  # Para HTML
# tabla = stargazer.render_latex()  # Para LaTeX

# Imprimir la tabla
stargazer


	Dependent variable: open

	Modelo Antes de 2008	Modelo Después de 2008	Modelo Completo
	(1)	(2)	(3)

const	0.491^***	1.089^***	-3.132^***
	(0.063)	(0.169)	(0.194)
x1	0.007^***	0.048^***	0.025^***
	(0.000)	(0.001)	(0.000)

Observations	370	403	773
R²	0.605	0.915	0.809
Adjusted R²	0.604	0.915	0.808
Residual Std. Error	0.611 (df=368)	1.701 (df=401)	2.695 (df=771)
F Statistic	563.964^*** (df=1; 368)	4302.783^*** (df=1; 401)	3256.333^*** (df=1; 771)

Note:	^p<0.1; ^p<0.05; ^**p<0.01