---
title: Data Screening
description: Visualización y análisis estadístico del feature dataset
authors: Ian Contreras
format:
  html:
    fig-format: svg
---

## Importaciones

### Librerías
```{python}
# | label: importing-libraries
# Manejo de datos y análisis
import numpy as np
import pandas as pd
import scipy.stats as stats

# Modelos estadísticos y econométricos
import statsmodels.api as sm
from statsmodels.graphics.tsaplots import plot_acf, plot_pacf
from statsmodels.tsa.stattools import adfuller
from pmdarima.arima import auto_arima
from arch import arch_model

# Visualización
import matplotlib.pyplot as plt

# Presentación en IPython
from IPython.display import Markdown, HTML
```

### Data
```{python}
# | label: importing-data
df = pd.read_csv(r'..\data\csv\irp.csv', parse_dates=[
                 'date'], index_col='date')
returns = df['price_return']
split_date = '2020-12-31'
R_test = df[df.index >= split_date]['price_return'].rolling(
    window=5).std().dropna()
```

## Análisis Exploratorio. 

### Serie Histórica
```{python}
# | label: fig-price-return-series
# | fig-cap: Serie de retornos diarios IRP-GOBIX.
fig, ax = plt.subplots(figsize=(12, 10))
ax.plot(df.index, df['price_return'])
ax.set_ylabel('Retorno Precio (BPS)')
ax.set_xlim([df.index.min(), df.index.max()])
ax.legend()
plt.show()
```

### Análisis descriptivo de la muestra
```{python}
# | label: return-descriptive-stats
# | fig-cap: Tabla de las estadísticas descriptiva de la serie de retornos
desc_stats = returns.describe()

skewness = returns.skew()
kurtosis = returns.kurtosis()
jb_test = sm.stats.jarque_bera(returns)

descriptive_table = pd.DataFrame({
    'Observations': [int(desc_stats['count'])],
    'Mean': [desc_stats['mean']],
    'Median': [desc_stats['50%']],
    'Std. Dev': [desc_stats['std']],
    'Skewness': [skewness],
    'Kurtosis': [kurtosis],
    'Jarque-Bera': [jb_test[0]],
    'Prob.': [jb_test[1]]
})
Markdown(descriptive_table.to_markdown(index=False))
```

### Ajuste de distribuciones
```{python}
# | label: fig-distribution-fitting
# | fig-cap: Ajuste de distribuciones a los retornos de precio.
distributions = [stats.norm, stats.t]

plt.figure(figsize=(10, 6))
plt.hist(returns, bins=50, density=True, alpha=0.6, label='Resultados Empíricos')

for dist in distributions:
    params = dist.fit(returns)
    x = np.linspace(returns.min(), returns.max(), 100)
    plt.plot(x, dist.pdf(x, *params), label=f'{dist.name} fit')

plt.xlabel('Retorno Precio (BPS)')
plt.ylabel('Verosimilitud')
plt.legend()
plt.show()
```

### Q-Q plot con respecto a distribución t
```{python}
# | label: fig-tqq-plot
# | fig-cap: Q-Q Plot de retornos de precio contra la distribución t.
t_params = stats.t.fit(returns)

fig, ax = plt.subplots(figsize=(12, 10))
stats.probplot(returns, dist="t", sparams=t_params, plot=plt)
plt.grid(True)
plt.show()
```

### Kolmorov Smirnov test para distribución t
```{python}
# | label: kolmorov-smirnov
ks_stat, p_value = stats.kstest(returns, 't', args=t_params)
print(f"Kolmogorov-Smirnov statistic: {ks_stat}")
print(f"P-value: {p_value}")
```

### Test de estacionariedad/ robustez
```{python}
# | label: adf
result = adfuller(returns)
print('ADF Statistic:', result[0])
print('p-value:', result[1])
```

### Test de autocorrelación
```{python}
# | label: fig-acf-pacf
# | fig-cap: Función de autocorrelación (ACF) y función de autocorrelación parcial (PACF) de los retornos.
fig, axes = plt.subplots(2, 1, figsize=(10, 8))
plot_acf(returns, lags=24, ax=axes[0])
axes[0].set_title('Función de Autocorrelación (ACF)')

plot_pacf(returns, lags=12, ax=axes[1])
axes[1].set_title('Función de Autocorrelación Parcial (PACF)')

plt.tight_layout()
plt.show()
```

## Resultados 
### Estimación del modelo ARIMA 
```{python}
# | label: arima-model
# | fig-cap: Modelo ARIMA de maxima verosimilitud para la serie de retornos. 
model_auto = auto_arima(returns)
model_auto.summary()
```

### Estimación del AR-GARCH
```{python}
# | label: garch-model-fitting
# | fig-cap: Ajuste del modelo Zero-Garch ala serie de retornos
ar = arch_model(returns, mean='Zero', vol='GARCH', dist='t')
res = ar.fit(last_obs=split_date)
```

### Resultados del AR-GARCH
```{python}
# | label: garch-model
# | fig-cap: Modelo Zero-Garch de la serie de retornos
HTML(res.summary().as_html())
```

### Residuos del AR-GARCH
```{python}
# | label: fig-garch-residuals
# | fig-cap: Residuos del modelo AR-GARCH.
fig = res.plot()
```

### Residuos contra volatilidad condicional
```{python}
# | label: fig-residuals-vs-volatility
# | fig-cap: Residuos contra la volatilidad condicional del AR-GARCH.
std_resid = res.resid / res.conditional_volatility
unit_var_resid = res.resid / res.resid.std()
df = pd.concat([std_resid, unit_var_resid], axis=1)
df.columns = ["Residuos Estandarizados", "Residuos de Varianza Unitaria"]
subplot = df.plot(kind="kde", xlim=(-4, 4))
```

## Evaluación Retrospectiva (Backtesting)

### Conversion de varianza predicha. 
```{python}
# | label: forecast
forecasts = res.forecast(horizon=1)
forecasted_std = np.sqrt(forecasts.variance)
forecasted_std = forecasted_std[forecasted_std.index >
                                split_date].loc[R_test.index, 'h.1']
```

### MAPE del modelo
```{python}
# | label: mape
def MAPE(actual, predicted):
    actual_std = actual.rolling(window=5).std().dropna()
    return np.mean(np.abs((actual_std - predicted) / actual_std)) * 100

mape_value = MAPE(R_test, forecasted_std)
print(f"MAPE on Test Set (based on moving std): {mape_value:.2f}%")
```

### Predichos vs Actuales Plot-Backtested
```{python}
# | label: fig-predicted-vs-actual-backtest
# | fig-cap: Desviación Estándar Móvil Real vs Pronosticada - Modelo Zero-GARCH
fig, ax = plt.subplots(figsize=(10, 6))
ax.plot(R_test.index, R_test,
        label='Desviación Estándar Móvil Real', color='blue')
ax.plot(forecasted_std.index,
        forecasted_std, label='Desviación Estándar Móvil Pronosticada', color='red', linestyle='--')
ax.set_xlim([forecasted_std.index.min(), forecasted_std.index.max()])
plt.legend()
plt.show()
```
