# !pip install arch
# !pip install tabulate
# !pip install xgboost

# Librerias para el desarrollo del proyecto
import math
import datetime
import warnings
import numpy as np
import pandas as pd
import plotly.express as px
import plotly.figure_factory as ff
import matplotlib.pyplot as plt
import seaborn as sns
from arch import arch_model
from statsmodels.graphics.tsaplots import plot_acf, plot_pacf
from ipywidgets import HBox, VBox
from tabulate import tabulate
from sklearn.metrics import mean_absolute_percentage_error
from sklearn.metrics import mean_squared_error

plt.style.use('fivethirtyeight')
# Ignorando algunas advertencias
warnings.filterwarnings('ignore')

import pandas as pd
# Leer los datos como un DataFrame de pandas
repsol_stock_raw_data = pd.read_csv("RepsolStockData20Years.csv")
# Calcular el cambio porcentual por día
returns_repsol = 100 * repsol_stock_raw_data.close.pct_change().dropna()
# Eliminar 0 resultados, hay un error en el conjunto de datos y durante 74 días, el mercado de valores estuvo cerrado, por lo que el retorno es 0
returns_repsol = returns_repsol.drop(returns_repsol[returns_repsol == 0].index)
#Mostrar las 3 primeras filas de la Serie
returns_repsol.head(3)

#Columna de rentabilidad media
returns_repsol.mean()

#Usando los precios de cierre brutos para representar gráficamente la evolución de las acciones
close_prices = pd.DataFrame(repsol_stock_raw_data["close"])
#Al dividir cada precio de cierre por el primer precio de nuestro conjunto de datos, calculamos la rentabilidad acumulada de cada día
cum_rets = close_prices / close_prices.iloc[0,:]
#Usando el módulo plotly.express podemos representar gráficamente nuestra nueva curva cum_rets
fig = px.line(cum_rets.iloc[:,:], width=1000, height=500)
#Añadiendo título
fig.update_layout(title_text='Rentabilidad acumulada de las acciones de Repsol (sin incluir dividendos) desde el 28/10/2002 hasta el 28/10/2022')

#La volatilidad diaria es la desviación estándar de los retornos
daily_volatility = returns_repsol.std()
#La volatilidad mensual es el resultado de multiplicar el vol diario * la raíz cuadrada de 21, esto se debe a que hay 21 días de negociación en un mes
monthly_volatility = math.sqrt(21) * daily_volatility
#La volatilidad anual es el resultado de multiplicar el vol diario * la raíz cuadrada de 252, esto se debe a que hay 252 días de negociación en un año
annual_volatility = math.sqrt(252) * daily_volatility
#Usando el paquete tabulate podemos imprimir una bonita tabla
print(tabulate([['Repsol',daily_volatility,monthly_volatility,annual_volatility]],headers = ['Volatilidad Diaria %', 'Volatilidad Mensual %', 'Volatilidad Anual %'],tablefmt = 'fancy_grid',stralign='centro',numalign='centro',floatfmt=".2f"))

#Podemos representar gráficamente los retornos diarios de Repsol mediante un gráfico lineal usando .plot de pandas
returns_repsol.plot(figsize =(16,5), title = 'Retornos diarios de Repsol');

#Creamos un gráfico de distribución utilizando plotly.figure_factory, rediseñamos los datos para tenerlos en un vector
return_dist_plot = ff.create_distplot([returns_repsol.values.reshape(-1)], group_labels = [' '])
#Especificamos el diseño del gráfico
return_dist_plot.update_layout(showlegend=False, title_text='Distribución de los rendimientos diarios de Repsol', width=1000, height=500)

#Imprimimos el gráfico
return_dist_plot.show()

#Defina un modelo GARCH (4,4) que utilice una distribución GED
model = arch_model(returns_repsol,dist="ged", vol = 'GARCH', p=4, q=4)
#Ajuste del modelo
model_fit = model.fit(disp='off')
#Resumen del modelo
model_fit.summary()

#Definir la serie completa como el modelo definido previamente
full_serie_garch = arch_model(returns_repsol,dist="ged", vol = 'GARCH', p=4, q=4)
#Ajustar el modelo para la serie completa
model_fit_full_serie = full_serie_garch.fit(disp='off')
#Gráficaremos contra la volatilidad móvil
rolling_vol = abs(returns_repsol.rolling(window=22, min_periods=22).std().dropna())
#Concatenar los valores verdaderos y los valores entrenados en nuestro modelo
garch_and_rolling_std = pd.concat([pd.DataFrame(model_fit_full_serie.conditional_volatility),rolling_vol.dropna()], axis=1).dropna()
#Gráficar
garch_and_rolling_std_plot = px.line(garch_and_rolling_std, title = 'GARCH vs. volatilidad continua de los retornos diarios TRAIN', width=1000, height=500)
#Impresión del gráfico
garch_and_rolling_std_plot.show()

#Utilizando un rango numérico de 251 para completar una lista de valores predichos, para cada día ajustamos un nuevo modelo con los mismos parámetros, pero agregando el último día.
test_size = 251
rolling_predictions = []

for i in range(test_size):
  train = returns_repsol[:-(test_size-i)]
  model = arch_model(train,dist="ged", vol = 'GARCH', p=4, q=4)
  model_fit = model.fit(disp='off')
  pred = model_fit.forecast(horizon=1, reindex = False)
  rolling_predictions.append(np.sqrt(pred.variance.values[-1,:][0]))

#Transformándolo en una serie
rolling_predictions = pd.Series(rolling_predictions, index= returns_repsol.dropna().index[-test_size:])
#Establecer parámetros del gráfico
plt.figure(figsize=(10,4))
#Datos verdaderos
true, = plt.plot((rolling_vol)[-test_size:])
#Datos predichos
preds, = plt.plot(rolling_predictions)
#Gráfico de los datos
plt.title('Predicción de volatilidad para los próximos 251 días de negociación - Pronóstico continuo TEST con GARCH(4,4)', fontsize=10)
#Añadir leyenda
plt.legend(['Volatilidad verdadera', 'Volatilidad predicha'], fontsize=16)

RMSE_Serie = mean_squared_error(garch_and_rolling_std["close"],garch_and_rolling_std["cond_vol"],squared=False)
MAPE_Serie = mean_absolute_percentage_error(garch_and_rolling_std["close"], garch_and_rolling_std["cond_vol"])
print(f"El RMSE de nuestro modelo GARCH en los datos de la serie completa es {round(RMSE_Serie,4)}")

print(f"El MAPE de nuestro modelo GARCH en los datos de la serie completa es {round(MAPE_Serie*100,2)}%")

true_vol = rolling_vol[-test_size:]
pred_vol = rolling_predictions
RMSE = mean_squared_error(true_vol, pred_vol,squared=False)
MAPE = mean_absolute_percentage_error(true_vol, pred_vol)
print(f"El RMSE de nuestro modelo GARCH en los datos predichos es {round(RMSE,4)}")

print(f"El MAPE de nuestro modelo GARCH en los datos predichos es {round(MAPE*100,2)}%")

print(tabulate([['RMSE',round(RMSE_Serie,6),round(RMSE,6)]],
headers = ['GARCH ENTRENAMIENTO',  'GARCH PREDICCIONES'],tablefmt = 'fancy_grid',stralign='center',numalign='center',floatfmt=".2f"))

print(tabulate([['MAPE',round(MAPE_Serie*100,2),round(MAPE*100,2)]],
headers = ['GARCH ENTRENAMIENTO %',  'GARCH PREDICCIONES %'],tablefmt = 'fancy_grid',stralign='center',numalign='center',floatfmt=".2f"))