import pyreadstat
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
from statsmodels.graphics.tsaplots import plot_acf, plot_pacf
from statsmodels.tsa.arima.model import ARIMA
from statsmodels.tsa.stattools import adfuller, acf, pacf
from statsmodels.stats.diagnostic import acorr_ljungbox
import seaborn as sns
import scipy.stats as statsANALISIS DE INTERVENCION DE SERIES TEMPORALES EN PYTHON
import pandas as pd
df=pd.read_spss('ipi.sav')
df.head()| Ipi | YEAR_ | MONTH_ | DATE_ | |
|---|---|---|---|---|
| 0 | 90.2 | 1988.0 | 1.0 | JAN 1988 |
| 1 | 96.5 | 1988.0 | 2.0 | FEB 1988 |
| 2 | 103.4 | 1988.0 | 3.0 | MAR 1988 |
| 3 | 95.0 | 1988.0 | 4.0 | APR 1988 |
| 4 | 101.2 | 1988.0 | 5.0 | MAY 1988 |
# Crear una columna de fechas en formato datetime a partir de YEAR_ y MONTH_
df['Fecha'] = pd.to_datetime(df['YEAR_'].astype(int).astype(str) + '-' + df['MONTH_'].astype(int).astype(str))
# Establecer la nueva columna 'Fecha' como índice del DataFrame
#df.set_index('Fecha', inplace=True)
# Establecer la nueva columna 'Fecha' como índice del DataFrame y agregar la frecuencia
df.set_index(pd.DatetimeIndex(df['Fecha'], freq='MS'), inplace=True)# Asegurarse de que los datos estén ordenados por fecha
df.sort_index(inplace=True)df| Ipi | YEAR_ | MONTH_ | DATE_ | Fecha | |
|---|---|---|---|---|---|
| Fecha | |||||
| 1988-01-01 | 90.2 | 1988.0 | 1.0 | JAN 1988 | 1988-01-01 |
| 1988-02-01 | 96.5 | 1988.0 | 2.0 | FEB 1988 | 1988-02-01 |
| 1988-03-01 | 103.4 | 1988.0 | 3.0 | MAR 1988 | 1988-03-01 |
| 1988-04-01 | 95.0 | 1988.0 | 4.0 | APR 1988 | 1988-04-01 |
| 1988-05-01 | 101.2 | 1988.0 | 5.0 | MAY 1988 | 1988-05-01 |
| ... | ... | ... | ... | ... | ... |
| 2000-08-01 | 86.9 | 2000.0 | 8.0 | AUG 2000 | 2000-08-01 |
| 2000-09-01 | 125.1 | 2000.0 | 9.0 | SEP 2000 | 2000-09-01 |
| 2000-10-01 | 126.8 | 2000.0 | 10.0 | OCT 2000 | 2000-10-01 |
| 2000-11-01 | 133.3 | 2000.0 | 11.0 | NOV 2000 | 2000-11-01 |
| 2000-12-01 | 112.3 | 2000.0 | 12.0 | DEC 2000 | 2000-12-01 |
156 rows × 5 columns
# Extraer la columna 'Ipi' como una serie de tiempo mensual
ipi_series = df['Ipi']# Graficar la serie temporal para verificar su formato
import matplotlib.pyplot as plt
plt.figure(figsize=(10, 6))
plt.plot(ipi_series)
plt.title("Índice de Producción Industrial (IPI) - Enero 1988 a Diciembre 2009")
plt.xlabel("Fecha")
plt.ylabel("IPI")
plt.show()
# Paso 3: Correlogramas (ACF y PACF) de la serie original
plot_acf(ipi_series, lags=40)
plt.title("Correlograma ACF de la serie original")
plt.show()
plot_pacf(ipi_series, lags=40)
plt.title("Correlograma PACF de la serie original")
plt.show()C:\Users\jipm1\AppData\Local\r-miniconda\envs\r-reticulate\lib\site-packages\statsmodels\graphics\tsaplots.py:348: FutureWarning:
The default method 'yw' can produce PACF values outside of the [-1,1] interval. After 0.13, the default will change tounadjusted Yule-Walker ('ywm'). You can use this method now by setting method='ywm'.
# Paso 4: Test de Dickey-Fuller en la serie original
adf_result = adfuller(ipi_series.dropna())
print("Test de Dickey-Fuller para la serie original:")
print(f"ADF Statistic: {adf_result[0]}")
print(f"p-value: {adf_result[1]}")
print(f"Critical Values: {adf_result[4]}")Test de Dickey-Fuller para la serie original:
ADF Statistic: -0.7475805058237533
p-value: 0.8340373085251793
Critical Values: {'1%': -3.4776006742422374, '5%': -2.882265832283648, '10%': -2.5778219289774156}# Paso 5: Diferenciar la serie (regular y estacionalmente)
# Diferenciación regular (d = 1)
ipi_diff = ipi_series.diff().dropna()# Diferenciación estacional (D = 1, s = 12)
ipi_diff_seasonal = ipi_diff.diff(12).dropna()# Graficar la serie diferenciada
plt.figure(figsize=(10, 6))
plt.plot(ipi_diff_seasonal)
plt.grid()
plt.title("Serie diferenciada (regular y estacional)")
plt.xlabel("Fecha")
plt.ylabel("IPI Diferenciado")
plt.show()
# Paso 6: Correlogramas de la serie diferenciada
plot_acf(ipi_diff_seasonal, lags=40)
plt.title("Correlograma ACF de la serie diferenciada")
plt.show()
plot_pacf(ipi_diff_seasonal, lags=40)
plt.title("Correlograma PACF de la serie diferenciada")
plt.show()C:\Users\jipm1\AppData\Local\r-miniconda\envs\r-reticulate\lib\site-packages\statsmodels\graphics\tsaplots.py:348: FutureWarning:
The default method 'yw' can produce PACF values outside of the [-1,1] interval. After 0.13, the default will change tounadjusted Yule-Walker ('ywm'). You can use this method now by setting method='ywm'.
# Paso 7: Test de Dickey-Fuller en la serie diferenciada
adf_result_diff = adfuller(ipi_diff_seasonal.dropna())
print("Test de Dickey-Fuller para la serie diferenciada:")
print(f"ADF Statistic: {adf_result_diff[0]}")
print(f"p-value: {adf_result_diff[1]}")
print(f"Critical Values: {adf_result_diff[4]}")Test de Dickey-Fuller para la serie diferenciada:
ADF Statistic: -3.9251125688389625
p-value: 0.0018576119413281822
Critical Values: {'1%': -3.482087964046026, '5%': -2.8842185101614626, '10%': -2.578864381347275}# Basado en los correlogramas, identificamos un modelo ARIMA
# Vamos a ajustar ARIMA(1, 1, 1)(1, 1, 1, 12) como un posible modelo
# Paso 8: Ajustar el modelo ARIMA
model = ARIMA(df['Ipi'], order=(2, 1, 0), seasonal_order=(0, 1, 1, 12))
model_fit = model.fit()# Mostrar el resumen del modelo ajustado
print(model_fit.summary()) SARIMAX Results
==========================================================================================
Dep. Variable: Ipi No. Observations: 156
Model: ARIMA(2, 1, 0)x(0, 1, [1], 12) Log Likelihood -383.984
Date: Tue, 10 Sep 2024 AIC 775.968
Time: 23:18:50 BIC 787.820
Sample: 01-01-1988 HQIC 780.784
- 12-01-2000
Covariance Type: opg
==============================================================================
coef std err z P>|z| [0.025 0.975]
------------------------------------------------------------------------------
ar.L1 -0.8655 0.082 -10.514 0.000 -1.027 -0.704
ar.L2 -0.5545 0.078 -7.085 0.000 -0.708 -0.401
ma.S.L12 -0.8323 0.112 -7.437 0.000 -1.052 -0.613
sigma2 11.3221 1.333 8.495 0.000 8.710 13.934
===================================================================================
Ljung-Box (L1) (Q): 0.00 Jarque-Bera (JB): 9.13
Prob(Q): 0.97 Prob(JB): 0.01
Heteroskedasticity (H): 1.19 Skew: 0.11
Prob(H) (two-sided): 0.55 Kurtosis: 4.22
===================================================================================
Warnings:
[1] Covariance matrix calculated using the outer product of gradients (complex-step).# Residuos del modelo ajustado
residuals = model_fit.resid# 1. Graficar los residuos
plt.figure(figsize=(10, 6))
plt.plot(residuals)
plt.axhline(0, color='red', linestyle='--', lw=2) # Media cero
plt.title("Gráfico de los residuos del modelo ARIMA")
plt.xlabel("Fecha")
plt.ylabel("Residuos")
plt.show()
#Correlogramas de los residuos (ACF y PACF)
plot_acf(residuals, lags=40)
plt.title("Correlograma ACF de los residuos")
plt.show()
plot_pacf(residuals, lags=40)
plt.title("Correlograma PACF de los residuos")
plt.show()C:\Users\jipm1\AppData\Local\r-miniconda\envs\r-reticulate\lib\site-packages\statsmodels\graphics\tsaplots.py:348: FutureWarning:
The default method 'yw' can produce PACF values outside of the [-1,1] interval. After 0.13, the default will change tounadjusted Yule-Walker ('ywm'). You can use this method now by setting method='ywm'.
# 3. Normalidad de los residuos: Histograma y Gráfico QQ
# Histograma de los residuos con curva KDE (Kernel Density Estimate)
plt.figure(figsize=(10, 6))
sns.histplot(residuals, kde=True, color='blue', stat="density")
plt.title("Histograma de los residuos con KDE")
plt.show()
# Gráfico QQ de los residuos para comprobar la normalidad
plt.figure(figsize=(10, 6))
stats.probplot(residuals, dist="norm", plot=plt)
plt.title("Gráfico QQ de los residuos")
plt.show()
# 4. Test de Ljung-Box para verificar si los residuos son ruido blanco (sin autocorrelación)
ljung_box_test = acorr_ljungbox(residuals, lags=[10], return_df=True)
print("Test de Ljung-Box para los residuos:")
print(ljung_box_test)Test de Ljung-Box para los residuos:
lb_stat lb_pvalue
10 24.180782 0.007135# 5. Comprobar media cero y varianza constante
print(f"Media de los residuos: {np.mean(residuals)}")
print(f"Varianza de los residuos: {np.var(residuals)}")Media de los residuos: 0.27260126919964844
Varianza de los residuos: 100.36891402541742# Paso 9: Pronosticar los próximos 12 meses
forecast = model_fit.get_forecast(steps=12)
forecast_index = pd.date_range(start=ipi_series.index[-1] + pd.DateOffset(months=1), periods=12, freq='M')# Extraer las predicciones y revertir la diferenciación regular y estacional
predicted_diff = forecast.predicted_mean
ci = forecast.conf_int()# Revertir diferenciación regular (cumsum sobre predicciones)
predicted_diff = ipi_series.iloc[-1] + predicted_diff.cumsum()# Revertir la diferenciación estacional: sumar los valores de 12 meses atrás
predictions = predicted_diff.copy()
for i in range(12):
if i < 12:
predictions[i] = predicted_diff[i] + ipi_series.iloc[-12 + i]
else:
predictions[i] = predicted_diff[i] + predictions[i - 12]# Mostrar los pronósticos en la escala original
plt.figure(figsize=(10, 6))
plt.plot(ipi_series, label="Serie original")
plt.plot(forecast_index, predictions, label="Predicciones (revertidas)", color='red')
plt.fill_between(forecast_index, ci.iloc[:, 0], ci.iloc[:, 1], color='pink', alpha=0.3)
plt.title("Pronóstico ARIMA del IPI (escala original)")
plt.xlabel("Fecha")
plt.ylabel("IPI")
plt.legend()
plt.show()
# Mostrar las predicciones revertidas
print("Pronósticos (revertidos a la escala original):")
print(predictions)Pronósticos (revertidos a la escala original):
2001-01-01 351.813894
2001-02-01 482.167792
2001-03-01 620.588591
2001-04-01 720.410870
2001-05-01 867.613029
2001-06-01 996.628027
2001-07-01 1123.167195
2001-08-01 1167.592555
2001-09-01 1332.911908
2001-10-01 1465.260428
2001-11-01 1602.924849
2001-12-01 1700.332116
Freq: MS, Name: predicted_mean, dtype: float64A CONTINUACION EL ANALISIS DE INTERVENCION DE LA SERIE TEMPORAL.
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
from statsmodels.tsa.arima.model import ARIMA
import seaborn as sns
# Paso 1: Leer el archivo .sav
df = pd.read_spss('ipi.sav')
# Crear una columna de fechas en formato datetime a partir de YEAR_ y MONTH_
df['Fecha'] = pd.to_datetime(df['YEAR_'].astype(int).astype(str) + '-' + df['MONTH_'].astype(int).astype(str))
# Establecer la nueva columna 'Fecha' como índice del DataFrame y agregar la frecuencia
df.set_index(pd.DatetimeIndex(df['Fecha'], freq='MS'), inplace=True)
# Asegurarse de que los datos estén ordenados por fecha
df.sort_index(inplace=True)
# Extraer la columna 'Ipi' como una serie de tiempo mensual
ipi_series = df['Ipi']# Paso 2: Identificación del evento de intervención
# el evento ocurrió en abril de 1997
fecha_intervencion = pd.Timestamp('1997-03-01')# Paso 3: Creación de la variable de intervención (dummy)
# La variable será 0 antes de enero 1997 y 1 después de enero 1997 abril
df['intervencion'] = np.where(df.index >= fecha_intervencion, 1, 0)# Graficar la serie con la intervención
plt.figure(figsize=(10, 6))
plt.plot(ipi_series, label='Serie original')
plt.axvline(fecha_intervencion, color='red', linestyle='--', label='Intervención (Abril 1997)')
plt.title("Índice de Producción Industrial (IPI) con intervención")
plt.xlabel("Fecha")
plt.ylabel("IPI")
plt.legend()
plt.show()
# Paso 4: Ajustar el modelo ARIMA con intervención
# El modelo ARIMA incluirá la variable de intervención como exógena
model = ARIMA(ipi_series, order=(2, 1, 0), seasonal_order=(0, 1, 1, 12), exog=df[['intervencion']])
model_fit = model.fit()# Mostrar el resumen del modelo ajustado
print(model_fit.summary()) SARIMAX Results
==========================================================================================
Dep. Variable: Ipi No. Observations: 156
Model: ARIMA(2, 1, 0)x(0, 1, [1], 12) Log Likelihood -382.295
Date: Tue, 10 Sep 2024 AIC 774.589
Time: 23:18:54 BIC 789.403
Sample: 01-01-1988 HQIC 780.609
- 12-01-2000
Covariance Type: opg
================================================================================
coef std err z P>|z| [0.025 0.975]
--------------------------------------------------------------------------------
intervencion 4.1948 1.126 3.724 0.000 1.987 6.403
ar.L1 -0.8900 0.083 -10.748 0.000 -1.052 -0.728
ar.L2 -0.5668 0.076 -7.447 0.000 -0.716 -0.418
ma.S.L12 -0.8363 0.112 -7.480 0.000 -1.055 -0.617
sigma2 11.0331 1.404 7.861 0.000 8.282 13.784
===================================================================================
Ljung-Box (L1) (Q): 0.04 Jarque-Bera (JB): 0.88
Prob(Q): 0.84 Prob(JB): 0.64
Heteroskedasticity (H): 1.12 Skew: -0.14
Prob(H) (two-sided): 0.70 Kurtosis: 3.26
===================================================================================
Warnings:
[1] Covariance matrix calculated using the outer product of gradients (complex-step).# Paso 5: Interpretación de los resultados
# El coeficiente de la variable de intervención indica el impacto del evento en la serie temporal
intervencion_coef = model_fit.params['intervencion']
print(f"Coeficiente de la intervención: {intervencion_coef}")Coeficiente de la intervención: 4.194795887472024# Paso 6: Pronósticos y comparación con el modelo sin intervención
# Hacer pronósticos con la intervención
forecast_with_intervention = model_fit.get_forecast(steps=12, exog=df[['intervencion']].iloc[-12:])
forecast_index = pd.date_range(start=ipi_series.index[-1] + pd.DateOffset(months=1), periods=12, freq='M')# Mostrar las predicciones con intervención
predicted_with_intervention = forecast_with_intervention.predicted_mean
predicted_with_intervention2001-01-01 121.193560
2001-02-01 123.658679
2001-03-01 126.818244
2001-04-01 121.090296
2001-05-01 128.453011
2001-06-01 128.961282
2001-07-01 130.407406
2001-08-01 85.574438
2001-09-01 126.722446
2001-10-01 130.306747
2001-11-01 130.724608
2001-12-01 118.038667
Freq: MS, Name: predicted_mean, dtype: float64# Graficar la serie original con los pronósticos con intervención
plt.figure(figsize=(10, 6))
plt.plot(ipi_series, label='Serie original')
plt.plot(forecast_index, predicted_with_intervention, label='Predicciones con intervención', color='red')
plt.axvline(fecha_intervencion, color='blue', linestyle='--', label='Intervención (Enero 2000)')
plt.title("Pronóstico ARIMA con intervención")
plt.xlabel("Fecha")
plt.ylabel("IPI")
plt.legend()
plt.show()
# ------------------- Reversión de diferenciaciones ------------------- #
# Revertir diferenciación regular (cumsum sobre predicciones)
predicted_diff = ipi_series.iloc[-1] + predicted_diff.cumsum()
# Revertir la diferenciación estacional: sumar los valores de 12 meses atrás
predictions = predicted_diff.copy()
for i in range(12):
if i < 12:
predictions[i] = predicted_diff[i] + ipi_series.iloc[-12 + i]
else:
predictions[i] = predicted_diff[i] + predictions[i - 12]# ------------------- Graficar los pronósticos ------------------- #
# Graficar la serie original con los pronósticos con intervención
plt.figure(figsize=(10, 6))
plt.plot(ipi_series, label='Serie original')
plt.plot(forecast_index, predictions, label='Predicciones (revertidas)', color='red')
plt.fill_between(forecast_index, ci.iloc[:, 0], ci.iloc[:, 1], color='pink', alpha=0.3)
plt.axvline(fecha_intervencion, color='blue', linestyle='--', label='Intervención (Enero 2000)')
plt.title("Pronóstico ARIMA con intervención (escala original)")
plt.xlabel("Fecha")
plt.ylabel("IPI")
plt.legend()
plt.show()
# Mostrar las predicciones revertidas a la escala original
print("Pronósticos (revertidos a la escala original):")
print(predictions)Pronósticos (revertidos a la escala original):
2001-01-01 464.113894
2001-02-01 827.781687
2001-03-01 1323.170278
2001-04-01 1907.281148
2001-05-01 2660.094178
2001-06-01 3523.622204
2001-07-01 4514.089399
2001-08-01 5553.181954
2001-09-01 6799.193862
2001-10-01 8139.354289
2001-11-01 9615.479138
2001-12-01 11182.511254
Freq: MS, Name: predicted_mean, dtype: float64