MAT281 - Laboratorio N°10¶

I.- Problema 01¶

Lista de actos delictivos registrados por el Service de police de la Ville de Montréal (SPVM).

El conjunto de datos en estudio interventionscitoyendo.csv corresponde a todos los delitos entre 2015 y agosto de 2020en Montreal. Cada delito está asociado en grandes categorías, y hay información sobre la ubicación, el momento del día, etc.

Nota: Para más información seguir el siguiente el link.

In [ ]:

# librerias 
import pandas as pd
import numpy as np 
import matplotlib.pyplot as plt 
import seaborn as sns 
import statsmodels.api as sm
from prophet import Prophet
from statsmodels.tsa.statespace.sarimax import SARIMAX

# graficos incrustados
sns.set_style('whitegrid')
%matplotlib inline

# parametros esteticos de seaborn
sns.set_palette("deep", desat=.6)
sns.set_context(rc={"figure.figsize": (12, 4)})

# librerias import pandas as pd import numpy as np import matplotlib.pyplot as plt import seaborn as sns import statsmodels.api as sm from prophet import Prophet from statsmodels.tsa.statespace.sarimax import SARIMAX # graficos incrustados sns.set_style('whitegrid') %matplotlib inline # parametros esteticos de seaborn sns.set_palette("deep", desat=.6) sns.set_context(rc={"figure.figsize": (12, 4)})

In [ ]:

from sklearn.metrics import mean_absolute_error, mean_squared_error

def mean_absolute_percentage_error(y_true, y_pred):
    return np.mean(np.abs((y_true - y_pred) / y_true)) * 100

def regression_metrics(df):
    """
    Aplicar las distintas métricas definidas
    :param df: DataFrame con las columnas: ['y', 'yhat']
    :return: DataFrame con las métricas especificadas
    """
    df_result = pd.DataFrame()

    y_true = df['y']
    y_pred = df['yhat']

    df_result['mae'] = [round(mean_absolute_error(y_true, y_pred), 4)]
    df_result['mse'] = [round(mean_squared_error(y_true, y_pred), 4)]
    df_result['rmse'] = [round(np.sqrt(mean_squared_error(y_true, y_pred)), 4)]
    df_result['mape'] = [round(mean_absolute_percentage_error(y_true, y_pred), 4)]
    df_result['smape'] = [round(2 * mean_absolute_percentage_error(y_true, y_pred) / (mean_absolute_percentage_error(y_true, y_pred) + 100), 4)]

    return df_result

from sklearn.metrics import mean_absolute_error, mean_squared_error def mean_absolute_percentage_error(y_true, y_pred): return np.mean(np.abs((y_true - y_pred) / y_true)) * 100 def regression_metrics(df): """ Aplicar las distintas métricas definidas :param df: DataFrame con las columnas: ['y', 'yhat'] :return: DataFrame con las métricas especificadas """ df_result = pd.DataFrame() y_true = df['y'] y_pred = df['yhat'] df_result['mae'] = [round(mean_absolute_error(y_true, y_pred), 4)] df_result['mse'] = [round(mean_squared_error(y_true, y_pred), 4)] df_result['rmse'] = [round(np.sqrt(mean_squared_error(y_true, y_pred)), 4)] df_result['mape'] = [round(mean_absolute_percentage_error(y_true, y_pred), 4)] df_result['smape'] = [round(2 * mean_absolute_percentage_error(y_true, y_pred) / (mean_absolute_percentage_error(y_true, y_pred) + 100), 4)] return df_result

In [ ]:

# read data

validate_categorie = [
  'Introduction', 'Méfait','Vol dans / sur véhicule à moteur', 'Vol de véhicule à moteur',
]

df = pd.read_csv("https://raw.githubusercontent.com/fralfaro/MAT281_2024/main/docs/labs/data/interventionscitoyendo.csv", sep=",", encoding='latin-1')
df.columns = df.columns.str.lower()
df['date'] = pd.to_datetime(df['date'], format='%Y-%m-%d')

df = df.loc[lambda x: x['categorie'].isin(validate_categorie)]
df = df.sort_values(['categorie','date'])
df.head()

# read data validate_categorie = [ 'Introduction', 'Méfait','Vol dans / sur véhicule à moteur', 'Vol de véhicule à moteur', ] df = pd.read_csv("https://raw.githubusercontent.com/fralfaro/MAT281_2024/main/docs/labs/data/interventionscitoyendo.csv", sep=",", encoding='latin-1') df.columns = df.columns.str.lower() df['date'] = pd.to_datetime(df['date'], format='%Y-%m-%d') df = df.loc[lambda x: x['categorie'].isin(validate_categorie)] df = df.sort_values(['categorie','date']) df.head()

Como tenemos muchos datos por categoría a nivel de día, agruparemos a nivel de semanas y separaremos cada serie temporal.

In [ ]:

cols = ['date','pdq']
y_s1 = df.loc[lambda x: x.categorie == validate_categorie[0] ][cols].set_index('date').resample('W').mean()
y_s2 = df.loc[lambda x: x.categorie == validate_categorie[1] ][cols].set_index('date').resample('W').mean()
y_s3 = df.loc[lambda x: x.categorie == validate_categorie[2] ][cols].set_index('date').resample('W').mean()
y_s4 = df.loc[lambda x: x.categorie == validate_categorie[3] ][cols].set_index('date').resample('W').mean()

cols = ['date','pdq'] y_s1 = df.loc[lambda x: x.categorie == validate_categorie[0] ][cols].set_index('date').resample('W').mean() y_s2 = df.loc[lambda x: x.categorie == validate_categorie[1] ][cols].set_index('date').resample('W').mean() y_s3 = df.loc[lambda x: x.categorie == validate_categorie[2] ][cols].set_index('date').resample('W').mean() y_s4 = df.loc[lambda x: x.categorie == validate_categorie[3] ][cols].set_index('date').resample('W').mean()

El objetivo de este laboratorio es poder realizar un análisis completo del conjunto de datos en estudio, para eso debe responder las siguientes preguntas:

1. Realizar un gráfico para cada serie temporal $y\_{si}, i =1,2,3,4$.

In [ ]:

# FIXME

# FIXME

2. Escoger alguna serie temporal $y\_{si}, i =1,2,3,4$. Luego:

• Realizar un análisis exploratorio de la serie temporal escogida
• Aplicar el modelo de pronóstico $SARIMA(p,d,q)x(P,D,Q,S)$, probando varias configuraciones de los hiperparámetros. Encuentre la mejor configuración. Concluya.
• Para el mejor modelo encontrado, verificar si el residuo corresponde a un ruido blanco.

Hint: Tome como target_date = '2021-01-01'. Recuerde considerar que su columna de valores se llama pdq.

In [ ]:

# Creando la clase SarimaModels
class SarimaModels:
    def __init__(self, params):
        self.params = params

    @property
    def name_model(self):
        return f"SARIMA_{self.params[0]}X{self.params[1]}".replace(' ', '')

    @staticmethod
    def test_train_model(y, date):
        mask_ds = y.index < date
        y_train = y[mask_ds]
        y_test = y[~mask_ds]
        return y_train, y_test

    def fit_model(self, y, date):
        y_train, y_test = self.test_train_model(y, date)
        model = SARIMAX(y_train,
                        order=self.params[0],
                        seasonal_order=self.params[1],
                        enforce_stationarity=False,
                        enforce_invertibility=False)
        model_fit = model.fit(disp=0)
        return model_fit

    def df_testing(self, y, date):
        y_train, y_test = self.test_train_model(y, date)
        model = SARIMAX(y_train,
                        order=self.params[0],
                        seasonal_order=self.params[1],
                        enforce_stationarity=False,
                        enforce_invertibility=False)
        model_fit = model.fit(disp=0)
        start_index = y_test.index.min()
        end_index = y_test.index.max()
        preds = model_fit.get_prediction(start=start_index, end=end_index, dynamic=False)
        df_temp = pd.DataFrame({
            'y': y_test['pdq'],
            'yhat': preds.predicted_mean
        })
        return df_temp

    def metrics(self, y, date):
        df_temp = self.df_testing(y, date)
        df_metrics = regression_metrics(df_temp)
        df_metrics['model'] = self.name_model
        return df_metrics

# Definir parámetros
import itertools
p = d = q = range(0, 2)
pdq = list(itertools.product(p, d, q))
seasonal_pdq = [(x[0], x[1], x[2], 12) for x in list(itertools.product(p, d, q))]
params = list(itertools.product(pdq, seasonal_pdq))
target_date = '2021-01-01'

# Creando la clase SarimaModels class SarimaModels: def __init__(self, params): self.params = params @property def name_model(self): return f"SARIMA_{self.params[0]}X{self.params[1]}".replace(' ', '') @staticmethod def test_train_model(y, date): mask_ds = y.index < date y_train = y[mask_ds] y_test = y[~mask_ds] return y_train, y_test def fit_model(self, y, date): y_train, y_test = self.test_train_model(y, date) model = SARIMAX(y_train, order=self.params[0], seasonal_order=self.params[1], enforce_stationarity=False, enforce_invertibility=False) model_fit = model.fit(disp=0) return model_fit def df_testing(self, y, date): y_train, y_test = self.test_train_model(y, date) model = SARIMAX(y_train, order=self.params[0], seasonal_order=self.params[1], enforce_stationarity=False, enforce_invertibility=False) model_fit = model.fit(disp=0) start_index = y_test.index.min() end_index = y_test.index.max() preds = model_fit.get_prediction(start=start_index, end=end_index, dynamic=False) df_temp = pd.DataFrame({ 'y': y_test['pdq'], 'yhat': preds.predicted_mean }) return df_temp def metrics(self, y, date): df_temp = self.df_testing(y, date) df_metrics = regression_metrics(df_temp) df_metrics['model'] = self.name_model return df_metrics # Definir parámetros import itertools p = d = q = range(0, 2) pdq = list(itertools.product(p, d, q)) seasonal_pdq = [(x[0], x[1], x[2], 12) for x in list(itertools.product(p, d, q))] params = list(itertools.product(pdq, seasonal_pdq)) target_date = '2021-01-01'

In [ ]:

# FIXME

# FIXME

3. Resuelva el ejercicio anterior utilizando la librería de Prophet.

In [ ]:

# FIXME

# FIXME