Sesión 5: Introducción a las Series de Tiempo

Curso: Análisis de datos

Author
Affiliation

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

Magíster en Data Science - Universidad del Desarrollo

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 (o longitudinal)
  • panel

Corte transversal

Cuando se empieza a estudiar modelamiento de datos, la aplicación natural de estos es en lo que llamamos datos de corte transversal. En estos, la unidad de observación son individuos separados unos de otros.

Usualmente, cada fila reprsenta un individuo (personas, familias, empresas) en un momento del tiempo específico (mismo año, mismo mes, etc). 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) o Encuesta de Caracterización Socioeconómica Nacional.

Esta es una encuesta de caracter nacional realizada a hogares, en la cual se busca levantar información para conocer la situación de los hogares y de la población, especialmente de aquella en situación de pobreza y de aquellos grupos definidos como prioritarios y permite evaluar impacto de política social. Es representativa a nivel regional.

Es realizada por el Ministerio de Desarrollo Social desde el año 1990 con una periodicidad; bianual o trianual. Hasta ahora, las encuestas aplicadas corresponden a los años 1990, 1992, 1994, 1996, 1998, 2000, 2003, 2006, 2009, 2011, 2013, 2105 y 2017.

# 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>
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
import matplotlib.pyplot as plt

# Supongamos que ya tienes un DataFrame df_casen2020 y quieres hacer un scatterplot

# Define el tamaño del gráfico
plt.figure(figsize=(10, 6))  # Ajusta el tamaño a tus preferencias (ancho, alto)

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

# Personaliza el gráfico si es necesario (títulos, etiquetas, etc.)

# Muestra el gráfico
plt.show()

Series de tiempo

En cambio, un tipo diferente de datos son las Series Temporales o tambien longitudinales.

En este tipo de información, se tiene diferentes momentos del tiempo (día, semana, mes, año, etc.) para una misma unidad de análisis (individuo, país, empresa, etc)

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

#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()
/Users/melita/anaconda3/lib/python3.10/site-packages/yahoo_fin/stock_info.py:29: UserWarning:

Parsing dates in DD/MM/YYYY format when dayfirst=False (the default) was specified. This may lead to inconsistently parsed dates! Specify a format to ensure consistent parsing.
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(de series de tiempo) 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 Xt sobre Yt?
  • ¿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

Ejemplo idotamente peligroso

(Fuente: https://www.tylervigen.com/spurious-correlations)

Trabajar con series de tiempo no es para nada sencillo: necesitaremos técnicas nuevas.

¿Por qué?

Necesitaremos modelos sofisticados de errores y técnicas para estimarlos, por ciertas dieferencias en la naturaleza del proceso generador de datos

Panel

El último tipo de datos, es lo que se llama un panel. Este es la combinacion de los dos tipos anteriores, donde para un conjunto de individuos tenemos varias observaciones en el tiempo.

Suele ser el tipo de datos más completo, pero tabien 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 etsos 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.

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

Corte transversal

Estimamos modelos de la forma: yi=β0+β1x1i+…+β1xki+ui y_{i}=\beta_{0}+\beta_{1}x_{1i}+\dots+\beta_{1}x_{ki}+u_{i}

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[ui|Xi]=0 E[u_{i}|X_{i}]=0

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

E[ui|Xj]=0 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.

La imposibilidad de la muestra aleatoria en series de tiempo

En serie de tiempo, nuestro universo es un proceso estocástico:

Corte transversal

Estimamos modelos de la forma:

yi=β0+β1x1i+…+β1xki+ui y_{i}=\beta_{0}+\beta_{1}x_{1i}+\dots+\beta_{1}x_{ki}+u_{i}

  • Estamos en presencia de un proceso estocástico
  • En cada momento se observa un posible resultado (o realización) del proceso estocásticos.

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[ut|Xs]=0∀sE[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.

¿Qué hacer entonces?

Reconocer los datos como procesos estocásticos e incluir sus particularidades en la modelación. De eso se tratarán los dos siguientes sesiones - Sesion 5: Peculiaridades de los procesos estocasticos y su exploración en la data - Sesion 6: Modelamiento

Conceptos propios 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:

Notación y transformaciones:

Notación

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

Rezagos

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

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
print(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: Δyt=yt−yt−1\Delta y_t = y_t- y_{t-1}

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: Δln(yt)=ln(yt)−ln(yt−1) \Delta ln(y_{t})=ln(y_{t})-ln(y_{t-1})

  • Cambio porcentual de yty_t entre t−1t -1 y t≈100×Δln(yt)t\approx 100\times \Delta ln(y_{t})

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

yt=β0+β1zt+uty_{t}=\beta_{0}+\beta_{1}z_{t}+u_{t}

Ejemplos: * Modelo Simple: inflaciónt=β0+β1desempleot+ut inflación_{t}=\beta_{0}+\beta_{1}desempleo_{t}+u_{t}

  • Modelo múltiple: homicidiot=β0+β1condenat+β2desmepleot+β3hombrest+ut homicidio_{t}=\beta_{0}+\beta_{1}condena_{t}+\beta_{2}desmepleo_{t}+\beta_{3}hombres_{t}+u_{t}

Modelo dinámico

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

inflaciónt=β0+β1desempleot+β2inflaciónt−1+ut inflación_{t}=\beta_{0}+\beta_{1}desempleo_{t}+\beta_{2} inflación_{t-1} +u_{t}

Pandas y series de tiempo

  • Pandas fue desarrollado originalmente para trabajar con datos financieros.

  • Como las series temporales son un tipo de datos comunes encontrados en aplicaiones de finanzas, naturalmente pandas tiene muy buen soporte para la mayoria de operaciones comunes.

  • pd.to_datetime() : conviere una serie o valor en un timestamp

  • Este formato permite un mejor manejo de las series.

#pandas remote data access support for calls to the World Bank Indicators API
from pandas_datareader import data, wb #instalar conda install pandas-datareader  o  pip installpandas-datareader
#Revisemos que indicadores hay disponibles. En este caso revisare de PIB (GDP en ingés), pero se pueden explorar muchas más opciones.

wb.search('gdp.*capita.*const')
id name unit source sourceNote sourceOrganization topics
691 6.0.GDPpc_constant GDP per capita, PPP (constant 2011 internation... LAC Equity Lab GDP per capita based on purchasing power parit... b'World Development Indicators (World Bank)' Economy & Growth
10978 NY.GDP.PCAP.KD GDP per capita (constant 2015 US$) World Development Indicators GDP per capita is gross domestic product divid... b'World Bank national accounts data, and OECD ... Economy & Growth
10980 NY.GDP.PCAP.KN GDP per capita (constant LCU) World Development Indicators GDP per capita is gross domestic product divid... b'World Bank national accounts data, and OECD ... Economy & Growth
10982 NY.GDP.PCAP.PP.KD GDP per capita, PPP (constant 2017 internation... World Development Indicators GDP per capita based on purchasing power parit... b'International Comparison Program, World Bank... Economy & Growth
10983 NY.GDP.PCAP.PP.KD.87 GDP per capita, PPP (constant 1987 internation... WDI Database Archives b'' 
# Obtengamos la lista de paises disponibles
countries=wb.get_countries()

#Preview primeras filas lista de paises
countries[:5]
iso3c iso2c name region adminregion incomeLevel lendingType capitalCity longitude latitude
0 ABW AW Aruba Latin America & Caribbean High income Not classified Oranjestad -70.0167 12.5167
1 AFE ZH Africa Eastern and Southern Aggregates Aggregates Aggregates NaN NaN
2 AFG AF Afghanistan South Asia South Asia Low income IDA Kabul 69.1761 34.5228
3 AFR A9 Africa Aggregates Aggregates Aggregates NaN NaN
4 AFW ZI Africa Western and Central Aggregates Aggregates Aggregates NaN NaN 
#sabemos que queremos Chile, asi que busquemos su info

countries[ countries['name'] == 'Chile' ]
iso3c iso2c name region adminregion incomeLevel lendingType capitalCity longitude latitude
49 CHL CL Chile Latin America & Caribbean High income IBRD Santiago -70.6475 -33.475 
# Descarguemos la data desde la API del banco muncial a un dataframe

df_GPDpc_Chile = wb.download(
                    #Use the indicator attribute to identify which indicator or indicators to download
                    indicator='NY.GDP.PCAP.KD',
                    #Use the country attribute to identify the countries you want data for
                    country=['CL'],
                    #Identify the first year for which you want the data, as an integer or a string
                    start='1980',
                    #Identify the last year for which you want the data, as an integer or a string
                    end=2020,
                )

#Veamos el data frame
df_GPDpc_Chile.head(10)
NY.GDP.PCAP.KD
country year
Chile 2020 12741.157507
2019 13761.374474
2018 13906.770558
2017 13615.523858
2016 13644.623261
2015 13569.948127
2014 13421.538342
2013 13318.595215
2012 13017.067840
2011 12382.379464 
wb.download( indicator=['NY.GDP.PCAP.PP.KD','NY.GDP.PCAP.KD'], country=['CL'], start=2008, end=2010 )
NY.GDP.PCAP.PP.KD NY.GDP.PCAP.KD
country year
Chile 2010 21225.103857 11773.004400
2009 20255.098784 11234.968211
2008 20695.562561 11479.281832 
sns.scatterplot(data=df_GPDpc_Chile, x="year", y="NY.GDP.PCAP.KD")
<Axes: xlabel='year', ylabel='NY.GDP.PCAP.KD'>

ax = df_GPDpc_Chile.plot(legend=True)

Ejemplos

1. PIB de EEUU

Empecemos con un ejemplo muy clásico de series de tiempo, con datos del PIB de Estados Unidos

GDP = pd.read_csv('data_sesion5/GDP.csv', parse_dates=['DATE'])
GDP.info()
<class 'pandas.core.frame.DataFrame'>
RangeIndex: 290 entries, 0 to 289
Data columns (total 2 columns):
 #   Column  Non-Null Count  Dtype         
---  ------  --------------  -----         
 0   DATE    290 non-null    datetime64[ns]
 1   GDP     290 non-null    float64       
dtypes: datetime64[ns](1), float64(1)
memory usage: 4.7 KB
GDP.set_index('DATE', inplace=True)
GDP.info()
<class 'pandas.core.frame.DataFrame'>
DatetimeIndex: 290 entries, 1947-01-01 to 2019-04-01
Data columns (total 1 columns):
 #   Column  Non-Null Count  Dtype  
---  ------  --------------  -----  
 0   GDP     290 non-null    float64
dtypes: float64(1)
memory usage: 4.5 KB
GDP.head()
GDP
DATE
1947-01-01 2033.061
1947-04-01 2027.639
1947-07-01 2023.452
1947-10-01 2055.103
1948-01-01 2086.017 
ax = GDP['2009':].plot(legend=False)
ax.set_ylabel(r'GDP ($\$B$)')
Text(0, 0.5, 'GDP ($\\$B$)')

Observamos claramente una tendencia al alza.

2. Mortalidad por influenza

Vemos un segundo ejemplo

ILI = pd.read_csv('data_sesion5/CDC.csv')
ILI.head()
Year Week Percent of Deaths Due to Pneumonia and Influenza Expected Threshold All Deaths Pneumonia Deaths Influenza Deaths
0 2012 1 8.479120 8.15718 8.49104 51102 4323 10
1 2012 2 8.343472 8.22181 8.55556 50962 4245 7
2 2012 3 8.370908 8.27534 8.60898 51010 4261 9
3 2012 4 8.448458 8.31696 8.65049 50163 4227 11
4 2012 5 8.140332 8.34602 8.67945 49568 4026 9 
ILI['date'] = ILI['Year']+ILI['Week']/52.
ILI.head()
Year Week Percent of Deaths Due to Pneumonia and Influenza Expected Threshold All Deaths Pneumonia Deaths Influenza Deaths date
0 2012 1 8.479120 8.15718 8.49104 51102 4323 10 2012.019231
1 2012 2 8.343472 8.22181 8.55556 50962 4245 7 2012.038462
2 2012 3 8.370908 8.27534 8.60898 51010 4261 9 2012.057692
3 2012 4 8.448458 8.31696 8.65049 50163 4227 11 2012.076923
4 2012 5 8.140332 8.34602 8.67945 49568 4026 9 2012.096154 
ILI.plot(x='date', y=['Percent of Deaths Due to Pneumonia and Influenza', 'Expected', 'Threshold'])
ax = plt.gca()
ax.legend(['Mortality', 'Expected', 'Threshold'])
ax.set_xlabel('Date')
ax.set_ylabel('% Mortality')
Text(0, 0.5, '% Mortality')

El comportamiento es muy diferente, podemos ver una estacionalidad.

Es decir hay un cambio periodico asociado al tiempo de la variable.

Transformaciones

En base a estos conceptos, entonces, es muy comun realizar transformaciones a las series de tiempo ya sea para obtener rezagos, diferencias, tasas de crecimiento, etc.

Ilustraremos varios de estos, con el dataset del DOw-Jonses indsutrial Average.

DJIA = pd.read_csv('data_sesion5/DJIA.csv', parse_dates=['DATE'], na_values='.').dropna()
DJIA.plot(x='DATE', legend=False)
ax = plt.gca()
ax.set_ylabel('DJIA')
ax.set_xlabel('Date')
Text(0.5, 0, 'Date')

Diferencias

Una manera de remover tendencias de un dataset es diferenciandolo. Como nuestros datasets son discretos, usaremos diferencias finitas.

def differentiate(values, d=1):
    # First value is required so that we can recover the original values with np.cumsum
    x = np.concatenate([[values[0]], values[1:]-values[:-1]])

    if d == 1:
        return x
    else:    
        return difference(x, d - 1)
values = DJIA['DJIA'].values
differences = differentiate(values)

Como podemos ver, el grafico se ve mucho más estacionario.

plt.plot(DJIA['DATE'].iloc[1:], differences[1:])
plt.xlabel('Date')
plt.ylabel('Differences')
Text(0, 0.5, 'Differences')

Para recuperar la data original, basta con integrar los puntos diferenciados.

def integrate(values, d=1):
    x = np.cumsum(values)
    
    if d == 1:
        return x
    else:
        return integrate(x, d-1)
rebuilt = integrate(differences)

Y con un chequeo rapido, podemos verificar que son los mismos valores.

np.mean(rebuilt-values)
0.0

Windowing

We also often want to calculate running values of some quantity. This requires the use of windowing functions that return the proper element at each step. def rolling(x, order): npoints = x.shape[0] running = []

for i in range(npoints-order+1):
    running.append(x[i:i+order])
    
return np.array(running)

And a simple example values = np.arange(11) values rolling(values, 6) Since we return a numpy array with all the individual windows, this also provides us with a simple way to take running averages by chaining methods rolling(values, 2) rolling(values, 2).mean(axis=1) Or the running maximum, etc. rolling(values, 2).max(axis=1)

Exponential Smoothing

Another form of smoothing a noise time series is called exponential smoothing. This is equivalent to an exponentially weighted running average where past values get exponentially reduced. def ES(values, alpha= 0.05): N = len(values) S = [values[0]*alpha]

for i in range(1, N):
    S.append(alpha*values[i]+(1-alpha)*S[-1])
    
return np.array(S)

As we can see with a few quick examples, the smaller the value of alpha the smoother (less noisy) the result smooth = [] smooth.append(ES(differences[1:], 0.01)) smooth.append(ES(differences[1:], 0.1)) smooth.append(ES(differences[1:], 0.5)) plt.plot(DJIA[‘DATE’].iloc[1:100], differences[1:100], label=‘Differences’) plt.plot(DJIA[‘DATE’].iloc[1:100], smooth[2][:99], label=r’α=0.5\alpha=0.5‘) plt.plot(DJIA[’DATE’].iloc[1:100], smooth[1][:99], label=r’α=0.1\alpha=0.1‘) plt.plot(DJIA[’DATE’].iloc[1:100], smooth[0][:99], label=r’α=0.01\alpha=0.01‘) plt.xlabel(’Date’) plt.ylabel(‘Differences’) plt.legend()

Missing Data

Desafortunadamente, los datos no siempre están limpios o completos, lo que nos obliga a lidiar datos faltantes. Aquí ilustramos varios enfoques para introducir valores perdidos. Comenzamos generando un conjunto de datos con valores perdidos.

x = np.linspace(-np.pi, np.pi, 100)
y = np.cos(x)
y_missing = y.copy()
y_missing[40:55] = np.nan

#This is simply a cosine function with a few missing values at the peak.

plt.plot(x, y, '*')
plt.plot(x, y_missing)

Quizás la estrategia más común es simplemente mantener el último valor ‘bueno’ conocido y usarlo para completar los puntos de datos faltantes. Este enfoque no puede lidiar con los valores faltantes al comienzo del conjunto de datos.

def ffill(y):
    y0 = y.copy()
    N = len(y0)
    
    current = None
    for i in range(1, N):
        if np.isnan(y0[i]):
            y0[i] = current
        else:
            current = y0[i]
    
    return y0

Naturalmente, el enfoque opuesto también es común cuando usamos el siguiente valor bueno. De esta manera podemos manejar fácilmente los valores perdidos iniciales, pero no podemos hacer nada con los valores perdidos al final de la serie de tiempo.

def bfill(y):
    y0 = y.copy()
    N = len(y0)
    
    current = None
    for i in range(N-1, 0, -1):
        if np.isnan(y0[i]):
            y0[i] = current
        else:
            current = y0[i]
    
    return y0

Back-fill y Forward-fill son enfoques simples pero poderosos para lidiar con los datos faltantes. Sin embargo, a menudo queremos tener más cuidado con el valor que atribuimos. Un enfoque común es interpolar entre el valor anterior y el siguiente y conectarlos con una línea recta.

def interpolate(y):
    y0 = y.copy()
    N = len(y0)
    
    pos = 0
    while pos < N:
        if np.isnan(y0[pos]):
            count = 0
            
            while np.isnan(y0[pos+count]):
                count += 1
            
            current = y0[pos-1]
            future = y0[pos+count]
            slope = (future-current)/count
            
            y0[pos:pos+count] = current + np.arange(1, count+1)*slope
            
            pos += count
        else:
            pos += 1
            
    return y0

La imputación de datos (el cálculo de los valores perdidos esperados) es un gran subcampo de estadísticas con una amplia gama de técnicas y enfoques.

y_bfill = bfill(y_missing)
y_ffill = ffill(y_missing)
y_inter = interpolate(y_missing)
# And a quick plot to visualize the differences
plt.plot(x, y_bfill, label='back fill')
plt.plot(x, y_ffill, label='forward fill')
plt.plot(x, y_inter, label='interpolate')
plt.plot(x, y_missing, label='Data')
plt.legend()

Resampling

En muchos casos, también necesitamos cambiar la frecuencia a la que estamos operando. Por ejemplo, nuestro conjunto de datos DJIA tiene valores al final del día, pero podríamos estar interesados en puntos de datos semanales o mensuales. El remuestreo es una serie de técnicas diseñadas para lidiar con esta situación y es similar en espíritu a las técnicas de ventanas que vimos anteriormente. La principal diferencia es que en lugar de simplemente mover la ventana en un paso fijo, cada ventana corresponde a nuestro período de interés.

mapping = DJIA['DATE'].dt.year
values = DJIA['DJIA'].values

En el caso más simple, simplemente calculamos cuál es la ventana correcta para cada punto de datos y la agregamos en consecuencia.

def groupBy(values, mapping, func = None):
    agg = {}
    pos = {}
    
    for i in range(values.shape[0]):
        key = mapping.iloc[i]
        
        if key not in agg:
            agg[key] = []
        
        pos[key] = i
        
        if not np.isnan(values[i]):
            agg[key].append(values[i])
        
    order = sorted(agg.keys())
    
    if func is not None:
        for key in agg:
            agg[key] = func(np.array(agg[key]).astype('float'))
            
    return agg, pos

Naturalmente, esta función groupBy es útil no solo para remuestrear sino también para una amplia gama de análisis estadísticos. Además de un mapeo, también debemos especificar qué función de agregación queremos usar. ¿Nos interesa el valor medio? ¿el maximo? ¿Desviación Estándar?

agg, pos = groupBy(values, mapping, np.mean)

Aca simplemente calculamos la media por año

agg
{2009: 10232.334385964914,
 2010: 10668.589087301589,
 2011: 11957.570238095239,
 2012: 12965.28744,
 2013: 15009.523492063492,
 2014: 16777.690912698414,
 2015: 17587.029166666664,
 2016: 17927.107341269842,
 2017: 21750.20374501992,
 2018: 25046.85734939759,
 2019: 26000.832040816327}

Como nuestra función groupBy también devuelve las posiciones de índice de la última vez que se vio cada bin, podemos comparar fácilmente los datos originales con los muestreados nuevamente.

aggregated = []

for key in pos:
    aggregated.append([pos[key], agg[key]])

aggregated = np.array(aggregated)
aggregated
plt.plot(DJIA['DATE'], DJIA['DJIA'])
ax = plt.gca()
ax.plot(DJIA.set_index('DATE').index[aggregated.T[0].astype('int')], aggregated.T[1])
ax.plot(DJIA.set_index('DATE').index[aggregated.T[0].astype('int')], aggregated.T[1], 'ro', markersize=10,)
ax.set_ylabel('DJIA')
ax.set_xlabel('Date')
Text(0.5, 0, 'Date')

Jackknife estimators

Finalmente, en muchos casos queremos estimar cantidades estadísticas de series de tiempo. Un ejemplo obvio podría ser estimar la media móvil de una serie sin tendencia para verificar si de hecho está lo suficientemente cerca de cero para ser considerada estacionaria.

El estimador JackKnife nos permite obtener no solo el valor esperado en cuestión, sino también una medida de su varianza. Esto se logra mediante el uso de un enfoque de omisión para calcular N estimaciones de nuestra métrica. A partir de esta población de estimación, podemos obtener los promedios como la mejor estimación posible y la desviación estándar como una medida de las barras de error involucradas.

def jackknife(x, func, variance = False):
    N = len(x)
    pos = np.arange(N)
    values = [func(x[pos != i]) for i in pos]
    jack = np.sum(values)/N
    
    if variance:
        values = [np.power(func(x[pos != i]) - jack, 2.0) for i in pos]
        var = (N-1)/N * np.sum(values)
        return jack, var
    else:
        return jack
x = np.random.normal(0, 2, 100)
print(x.std())
jackknife(x, np.std, True)
1.8992867419484192
(1.8991395820804744, 0.01890134168632936)

Con jackknife obtenemos no solo una estimación del valor sino también una medida del error

Bootstrapping

Otra técnica común para estimar propiedades estadísticas se conoce como bootstrapping y está estrechamente relacionada con Jackknife. En este enfoque, simplemente tomamos muestras (con reemplazo) de la población original para obtener una medida de cuánta variabilidad se puede esperar.

def bootstrapping(x, n_samples, func=np.mean):
    y = x.copy()
    N = len(y)
    population = []
    
    for i in range(n_samples):
        population.append(func(np.random.choice(y, N, replace=True)))
        
    return np.array(population)

Podemos generar fácilmente un histograma de las muestras bostrappeadas

def histogram(values, n_bins=100):
    xmax = values.max()
    xmin = values.min()
    delta  = (xmax-xmin)/n_bins
    
    counts = np.zeros(n_bins+1, dtype='int')
    
    for value in values:
        val_bin = np.around((value-xmin)/delta).astype('int')
        counts[val_bin] += 1.0
    
    bins = xmin+delta*np.arange(n_bins+1)
    
    return bins, counts/values.shape[0]
x = np.random.normal(0, 2, size=100)
boot = bootstrapping(x, 1000)
x.mean()
-0.08580476619590074

Actividad:

  1. Siga este ejemplo practico de importar datos desde la API del Banco Mundial y preparar la base para su análisis de series de tiempo.
  2. Importe usted la serie de GDP total Y Percapita para otro país serie desde la API del Banco mundial, muestre sus principales características y realice un grafico.
  3. Obtenga las primeras diferencias y compare, ¿pareciera haber tendencias?

Desarrollo

#pandas remote data access support for calls to the World Bank Indicators API
from pandas_datareader import data, wb #instalar conda install pandas-datareader  o  pip installpandas-datareader

#Revisemos que indicadores hay disponibles. En este caso revisare de PIB (GDP en ingés), pero se pueden explorar muchas más opciones.

wb.search('gdp.*capita.*const')
id name unit source sourceNote sourceOrganization topics
691 6.0.GDPpc_constant GDP per capita, PPP (constant 2011 internation... LAC Equity Lab GDP per capita based on purchasing power parit... b'World Development Indicators (World Bank)' Economy & Growth
10978 NY.GDP.PCAP.KD GDP per capita (constant 2015 US$) World Development Indicators GDP per capita is gross domestic product divid... b'World Bank national accounts data, and OECD ... Economy & Growth
10980 NY.GDP.PCAP.KN GDP per capita (constant LCU) World Development Indicators GDP per capita is gross domestic product divid... b'World Bank national accounts data, and OECD ... Economy & Growth
10982 NY.GDP.PCAP.PP.KD GDP per capita, PPP (constant 2017 internation... World Development Indicators GDP per capita based on purchasing power parit... b'International Comparison Program, World Bank... Economy & Growth
10983 NY.GDP.PCAP.PP.KD.87 GDP per capita, PPP (constant 1987 internation... WDI Database Archives b'' 
# Obtengamos la lista de paises disponibles
countries=wb.get_countries()

#Preview primeras filas lista de paises
countries[:5]
iso3c iso2c name region adminregion incomeLevel lendingType capitalCity longitude latitude
0 ABW AW Aruba Latin America & Caribbean High income Not classified Oranjestad -70.0167 12.5167
1 AFE ZH Africa Eastern and Southern Aggregates Aggregates Aggregates NaN NaN
2 AFG AF Afghanistan South Asia South Asia Low income IDA Kabul 69.1761 34.5228
3 AFR A9 Africa Aggregates Aggregates Aggregates NaN NaN
4 AFW ZI Africa Western and Central Aggregates Aggregates Aggregates NaN NaN 
#sabemos que queremos Chile, asi que busquemos su info

countries[ countries['name'] == 'Chile' ]
iso3c iso2c name region adminregion incomeLevel lendingType capitalCity longitude latitude
49 CHL CL Chile Latin America & Caribbean High income IBRD Santiago -70.6475 -33.475 
# Descarguemos la data desde la API del banco muncial a un dataframe

df_GPDpc_Chile = wb.download(
                    #Use the indicator attribute to identify which indicator or indicators to download
                    indicator='NY.GDP.PCAP.KD',
                    #Use the country attribute to identify the countries you want data for
                    country=['CL'],
                    #Identify the first year for which you want the data, as an integer or a string
                    start='1980',
                    #Identify the last year for which you want the data, as an integer or a string
                    end=2020
                )

#Veamos el data frame
df_GPDpc_Chile.head(10)
NY.GDP.PCAP.KD
country year
Chile 2020 12741.157507
2019 13761.374474
2018 13906.770558
2017 13615.523858
2016 13644.623261
2015 13569.948127
2014 13421.538342
2013 13318.595215
2012 13017.067840
2011 12382.379464 
wb.download( indicator=['NY.GDP.PCAP.PP.KD','NY.GDP.PCAP.KD'], country=['CL'], start=2008, end=2010 )
NY.GDP.PCAP.PP.KD NY.GDP.PCAP.KD
country year
Chile 2010 21225.103857 11773.004400
2009 20255.098784 11234.968211
2008 20695.562561 11479.281832 
sns.scatterplot(data=df_GPDpc_Chile, x="year", y="NY.GDP.PCAP.KD")
<Axes: xlabel='year', ylabel='NY.GDP.PCAP.KD'>

Transformemos la columna de tiempo a un indice en el dataframe

#df_GPDpc_Chile.index=pd.to_datetime(df_GPDpc_Chile['year'], format="%Y" )

Supuestos no cumplidos y opciones

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:

  1. Modelo lineal en parámetros yt=β0+β1x1t+…+β1xkt+uty_{t}=\beta_{0}+\beta_{1}x_{1t}+\dots+\beta_{1}x_{kt}+u_{t}

  2. Media condicionada nula o exogeneidad (que ahora es estricta”) E[ut|Xjs]=0∀s E[u_{t}|X_{js}]=0 \quad \forall s

  3. No hay colinealidad perfecta.

  4. Homoscedasticidad Var[ut|Xjs]=σ∀s Var[u_{t}|X_{js}]=\sigma \quad \forall s

5 No hay correlación serial. Cov[ut,us|Xjs]=σ∀s≠t Cov[u_{t}, u_{s}|X_{js}]=\sigma \quad \forall s\ne t

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 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 este bloque nos enfocaremos en la primera aplicación.

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.

Patrones de dependencia intertemporal

Primera aproximación a series de tiempo desde la intuiciónn - Con las series Tendremos dos casos: - La serie es estacionaria: bien comportada y el pasado describe bien el futuro. - Tenemos particularidades temporales actuando en el proceso estocástico:

* Efecto rezagados (variable independiente rezagada)
* Retroalimentación entre variables
* Inercia o auto-dependencia la variable depende de si misma en el pasado

Estacionareidad

  • Entenderemos estacionareidad como el futuro se parece al pasado, al menos en un sentido probabilistico.
  • En una serie estacionaria:
    • Nunca se aleja mucho de la media (media constante)
    • El efecto de los errores caen a lo largo del tiempo (varianza constante)
    • Lo que ocurre recientemente es relativamente más importante que lo que pasó en momentos más alejados (covarianza constante)
  • Si las series son estacionarias, no necesitamos muestreo aleatorio y podemos seguir usando OLS como siempre.
  • (En términos técnicos: Podemos reemplazar el supuesto de exogeneidad extricta por exogeneidad debil SI la serie es estacionaria.)

Estacionareidad

::: {.callout-note collapse=“true” incremental=“false”} ## Condiciones matemáticas:

  1. La media es constante a lo largo del tiempo, es decir, μt=μ\mu_t = \mu para todo tt.
  2. La varianza es constante a lo largo del tiempo, es decir, σt2=σ2\sigma^2_t = \sigma^2 para todo tt.
  3. La covarianza no depende del tiempo, lo que implica que Cov(Xt,Xt+k)=Cov(Xs,Xs+k)Cov(X_t, X_{t+k}) = Cov(X_s, X_{s+k}) para cualquier par de instantes tt y ss, y para cualquier desfase (lag)

Ejemplo

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 estacionareidad 1: Tendencias

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

# 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)

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

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')

Podemos identificar la tendencia en los datos, al calcular la media movil.

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. Ahora podemos agregar esta línea de tendencia al gráfico anterior.

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()

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

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

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 disvision 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.

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.

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()
/Users/melita/anaconda3/lib/python3.10/site-packages/yahoo_fin/stock_info.py:29: UserWarning:

Parsing dates in DD/MM/YYYY format when dayfirst=False (the default) was specified. This may lead to inconsistently parsed dates! Specify a format to ensure consistent parsing.
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

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()

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()

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