<jupyter_start><jupyter_text>Import modules<jupyter_code>import numpy as np
import pandas as pd
import torch
import torch.nn as nn

from sklearn.model_selection import train_test_split, cross_validate, cross_val_score
from sklearn.base import BaseEstimator, TransformerMixin
from sklearn.compose import ColumnTransformer, make_column_transformer
from sklearn.preprocessing import StandardScaler, MinMaxScaler
from sklearn.pipeline import Pipeline, make_pipeline
from sklearn.metrics import make_scorer
from sklearn.dummy import DummyRegressor

from xgboost import XGBRegressor
from sklearn.ensemble import RandomForestRegressor
from sklearn.linear_model import Ridge
import statsmodels.api as sm
import statsmodels.api as sm
import matplotlib.pyplot as plt

import altair as alt
# Save a vega-lite spec and a PNG blob for each plot in the notebook
alt.renderers.enable('mimetype')
# Handle large data sets without embedding them in the notebook
alt.data_transformers.enable('data_server')<jupyter_output><empty_output><jupyter_text>Utility Functions and Classes<jupyter_code># Function to plot the raw data
def plot_data(data, x_axis, y_axis):
    """
    Plots line chart of the time series data (y-axis) against time (x-axis)

    Parameters
    ----------
    data: pd.DataFrame
        input dataframe
    x_axis: str
        the name of the time(datetime) column 
    y_axis: str
        the name of the time series data column

    Returns
    ----------
    lines: altair.vegalite.v4.api.Chart
        altair chart object 
    """
    data = data.reset_index()
    lines = alt.Chart(data).mark_line().encode(
        x=alt.X(x_axis),
        y=alt.Y(y_axis)
    ).properties(
        width=700,
        height=250
    )
    return lines
# Generate lags
def generate_lags(df, num_lags):
    """
    Creates new features which are lags of the target variable
    Parameters
    ----------
    df: pd.DataFrame
        input dataframe
    num_lags: int
        the number of lags

    Returns
    ----------
    df_tmp: pd.DataFrame
        dataframe with lagged features created
    """
    df_tmp = df.copy()
    for n in range(1, num_lags+1):
        df_tmp[f'lag{n}'] = df_tmp['target'].shift(n)
    df_tmp = df_tmp.iloc[num_lags:]
    return df_tmp
# Custom transformer for automated feature extraction
class UniveriateTimeSeriesAddFeatures(BaseEstimator, TransformerMixin):
    def __init__(self, num_feats=4):
        self.num_feats = num_feats
        
    def __extract_features(self, X):
        """
        Creates new features from timestamp
        Parameters
        ----------
        X: pd.DataFrame
            input dataframe
        
        Returns
        ----------
        df_features: pd.DataFrame
            dataframe with new features created
        """
        df_features = X.copy()
        df_features = (
                df_features
                .assign(hour = df_features.index.hour)
                .assign(day = df_features.index.day)
                .assign(month = df_features.index.month)
                .assign(day_of_week = df_features.index.dayofweek)
                # .assign(week_of_year = df.index.isocalendar().week)
        )
        df_features = df_features.reset_index()
        return df_features
    
    def fit(self, X, y=None):
        return self

    def transform(self, X, y=None):
        df_features = X.copy()
        df_features = self.__extract_features(df_features)
        df_features = df_features.iloc[:,1:]
        return df_features
# Function to collate results
def mean_std_cross_val_scores(model, X_train, y_train, **kwargs):
    """
    Returns mean and std of cross validation

    Parameters
    ----------
    model :
        scikit-learn model
    X_train : numpy array or pandas DataFrame
        X in the training data
    y_train :
        y in the training data

    Returns
    ----------
        pandas Series with mean scores from cross_validation
    """

    scores = cross_validate(model, X_train, y_train, **kwargs)

    mean_scores = pd.DataFrame(scores).mean()
    std_scores = pd.DataFrame(scores).std()
    out_col = []

    for i in range(len(mean_scores)):
        out_col.append((f"%0.3f (+/- %0.3f)" % (mean_scores[i], std_scores[i])))

    return pd.Series(data=out_col, index=mean_scores.index)
# Function to calculate MAPE
def mape(true, pred):
    return 100.0 * np.mean(np.abs((pred - true) / true))<jupyter_output><empty_output><jupyter_text>Read the data<jupyter_code># Read the data
df = pd.read_csv('data/PJME_hourly.csv')

# Remove spaces from column names
df.columns = df.columns.str.replace(' ', '')

# Convert Datetime column to datetime data type
df['Datetime'] = pd.to_datetime(df['Datetime'])

# Set index to datetime column
df = df.set_index('Datetime')

# Sort the index to ascending order
df = df.sort_index()

# Rename the PJME_MW column to target column
df = df.rename(columns={'PJME_MW': 'target'})

# Check head of the resulting dataframe
df.head()<jupyter_output><empty_output><jupyter_text>Split into train and test<jupyter_code>df_train, df_test = train_test_split(df, test_size=0.2, shuffle=False, random_state=123)
print(df_train.shape)
print(df_test.shape)<jupyter_output>(116292, 1)
(29074, 1)<jupyter_text>Visualize the data - Basic EDA<jupyter_code>plot_data(df_train, x_axis='Datetime', y_axis='target')<jupyter_output><empty_output><jupyter_text>ACF and PACF plots<jupyter_code>ts = df_train.target

fig,ax = plt.subplots(2,1,figsize=(8,8))
fig = sm.graphics.tsa.plot_acf(ts, lags=100, ax=ax[0])
fig = sm.graphics.tsa.plot_pacf(ts, lags=100, ax=ax[1])
plt.show()<jupyter_output>/Users/arijeetchatterjee/opt/miniconda3/envs/571/lib/python3.9/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'.
  warnings.warn(<jupyter_text>Separate train and test   Total Lags = 10<jupyter_code># Total lags defined
tot_lags = 10
df_train = generate_lags(df=df_train, num_lags=tot_lags)
df_test = generate_lags(df=df_test, num_lags=tot_lags)
df_test.head()
X_train, y_train = (
    df_train.drop(columns=['target']),
    df_train['target']
)
X_test, y_test = (
    df_test.drop(columns=['target']),
    df_test['target']
)

# Print shapes of the resulting train and test data sets
print(f'Shape of X_train and y_train: {X_train.shape, y_train.shape}')
print(f'Shape of X_test and y_test: {X_test.shape, y_test.shape}')<jupyter_output>Shape of X_train and y_train: ((116282, 10), (116282,))
Shape of X_test and y_test: ((29064, 10), (29064,))<jupyter_text>Feature Engineering Using Pipeline Transformations<jupyter_code>preprocessor = Pipeline(
    [
        ('add_lag_step', UniveriateTimeSeriesAddFeatures()),
        ('stdscaler', StandardScaler())
    ]
)<jupyter_output><empty_output><jupyter_text>Models<jupyter_code># make a scorer function that we can pass into cross-validation
mape_scorer = make_scorer(mape, greater_is_better=False)

scoring_metrics = {
    "neg RMSE": "neg_root_mean_squared_error",
    "r2": "r2",
    "mape": mape_scorer,
}
results = {}<jupyter_output><empty_output><jupyter_text>Dummy Regressor<jupyter_code>dummy = DummyRegressor()

results["Dummy"] = mean_std_cross_val_scores(
    dummy, X_train, y_train, return_train_score=True, scoring=scoring_metrics
)

pd.DataFrame(results)<jupyter_output><empty_output><jupyter_text>Other Models<jupyter_code>pipe_ridge = make_pipeline(preprocessor, Ridge(random_state=123))
pipe_xgb = make_pipeline(preprocessor, XGBRegressor(random_state=123))
pipe_rf = make_pipeline(preprocessor, RandomForestRegressor(random_state=123))

models = {
    'ridge': pipe_ridge,
    'xgboost': pipe_xgb,
    'rf_regressor': pipe_rf
}
for (name, model) in models.items():
    results[name] = mean_std_cross_val_scores(
        model, X_train, y_train, return_train_score=True, scoring=scoring_metrics
    )

pd.DataFrame(results)<jupyter_output><empty_output><jupyter_text>Best Model (Example of XGBRegressor provided below)<jupyter_code>pipe_xgb.fit(X_train, y_train);
preds = pipe_xgb.predict(X_test)
actual_vs_pred_df = y_test.reset_index()
actual_vs_pred_df['predicted'] = preds
actual_vs_pred_df.head()
actual_vs_pred_df_long = pd.melt(actual_vs_pred_df, id_vars='Datetime')

actual_vs_pred_df_long.tail()
# Visualize the actual vs predicted chart for the last 2 days of the test data period
alt.Chart(actual_vs_pred_df_long[actual_vs_pred_df_long.Datetime > '2018-07-24'], 
          title='Actual vs Predicted on Test Data').mark_line().encode(
    x=alt.X('Datetime', type='temporal', title='Datetime'),
    y=alt.Y('value', type='quantitative', title='Actual vs Predicted'),
    color=alt.Color('variable', title='Type')
).properties(
    width=700
).interactive()
print(f'MAPE for test set: {mape(actual_vs_pred_df["target"], actual_vs_pred_df["predicted"]):.2f}%')<jupyter_output>MAPE for test set: 0.82%<jupyter_text>Total Lags = 50<jupyter_code>df_train, df_test = train_test_split(df, test_size=0.2, shuffle=False, random_state=123)

# Total lags defined
tot_lags = 50
df_train = generate_lags(df=df_train, num_lags=tot_lags)
df_test = generate_lags(df=df_test, num_lags=tot_lags)

X_train, y_train = (
    df_train.drop(columns=['target']),
    df_train['target']
)
X_test, y_test = (
    df_test.drop(columns=['target']),
    df_test['target']
)

# Print shapes of the resulting train and test data sets
print(f'Shape of X_train and y_train: {X_train.shape, y_train.shape}')
print(f'Shape of X_test and y_test: {X_test.shape, y_test.shape}')
preprocessor = Pipeline(
    [
        ('add_lag_step', UniveriateTimeSeriesAddFeatures()),
        ('stdscaler', StandardScaler())
    ]
)
results = {}
dummy = DummyRegressor()

results["Dummy"] = mean_std_cross_val_scores(
    dummy, X_train, y_train, return_train_score=True, scoring=scoring_metrics
)

pd.DataFrame(results)
pipe_ridge = make_pipeline(preprocessor, Ridge(random_state=123))
pipe_xgb = make_pipeline(preprocessor, XGBRegressor(random_state=123))
pipe_rf = make_pipeline(preprocessor, RandomForestRegressor(random_state=123))

models = {
    'ridge': pipe_ridge,
    'xgboost': pipe_xgb,
    'rf_regressor': pipe_rf
}
for (name, model) in models.items():
    results[name] = mean_std_cross_val_scores(
        model, X_train, y_train, return_train_score=True, scoring=scoring_metrics
    )

pd.DataFrame(results)
pipe_xgb.fit(X_train, y_train);

preds = pipe_xgb.predict(X_test)

actual_vs_pred_df = y_test.reset_index()
actual_vs_pred_df['predicted'] = preds
# actual_vs_pred_df.head()

actual_vs_pred_df_long = pd.melt(actual_vs_pred_df, id_vars='Datetime')

# actual_vs_pred_df_long.tail()
# Visualize the actual vs predicted chart for the last 2 days of the test data period
alt.Chart(actual_vs_pred_df_long[actual_vs_pred_df_long.Datetime > '2018-07-24'], 
          title='Actual vs Predicted on Test Data').mark_line().encode(
    x=alt.X('Datetime', type='temporal', title='Datetime'),
    y=alt.Y('value', type='quantitative', title='Actual vs Predicted'),
    color=alt.Color('variable', title='Type')
).properties(
    width=700
).interactive()
print(f'MAPE for test set: {mape(actual_vs_pred_df["target"], actual_vs_pred_df["predicted"]):.2f}%')<jupyter_output>MAPE for test set: 0.89%<jupyter_text>Total Lags = 5<jupyter_code>df_train, df_test = train_test_split(df, test_size=0.2, shuffle=False, random_state=123)

# Total lags defined
tot_lags = 5
df_train = generate_lags(df=df_train, num_lags=tot_lags)
df_test = generate_lags(df=df_test, num_lags=tot_lags)

X_train, y_train = (
    df_train.drop(columns=['target']),
    df_train['target']
)
X_test, y_test = (
    df_test.drop(columns=['target']),
    df_test['target']
)

# Print shapes of the resulting train and test data sets
print(f'Shape of X_train and y_train: {X_train.shape, y_train.shape}')
print(f'Shape of X_test and y_test: {X_test.shape, y_test.shape}')
preprocessor = Pipeline(
    [
        ('add_lag_step', UniveriateTimeSeriesAddFeatures()),
        ('stdscaler', StandardScaler())
    ]
)
results = {}
dummy = DummyRegressor()

results["Dummy"] = mean_std_cross_val_scores(
    dummy, X_train, y_train, return_train_score=True, scoring=scoring_metrics
)

pd.DataFrame(results)
pipe_ridge = make_pipeline(preprocessor, Ridge(random_state=123))
pipe_xgb = make_pipeline(preprocessor, XGBRegressor(random_state=123))
pipe_rf = make_pipeline(preprocessor, RandomForestRegressor(random_state=123))

models = {
    'ridge': pipe_ridge,
    'xgboost': pipe_xgb,
    'rf_regressor': pipe_rf
}
for (name, model) in models.items():
    results[name] = mean_std_cross_val_scores(
        model, X_train, y_train, return_train_score=True, scoring=scoring_metrics
    )

pd.DataFrame(results)
pipe_xgb.fit(X_train, y_train);

preds = pipe_xgb.predict(X_test)

actual_vs_pred_df = y_test.reset_index()
actual_vs_pred_df['predicted'] = preds
# actual_vs_pred_df.head()

actual_vs_pred_df_long = pd.melt(actual_vs_pred_df, id_vars='Datetime')

# actual_vs_pred_df_long.tail()
# Visualize the actual vs predicted chart for the last 2 days of the test data period
alt.Chart(actual_vs_pred_df_long[actual_vs_pred_df_long.Datetime > '2018-07-24'], 
          title='Actual vs Predicted on Test Data').mark_line().encode(
    x=alt.X('Datetime', type='temporal', title='Datetime'),
    y=alt.Y('value', type='quantitative', title='Actual vs Predicted'),
    color=alt.Color('variable', title='Type')
).properties(
    width=700
).interactive()
print(f'MAPE for test set: {mape(actual_vs_pred_df["target"], actual_vs_pred_df["predicted"]):.2f}%')<jupyter_output>MAPE for test set: 0.82%<jupyter_text>Total Lags = 30<jupyter_code>df_train, df_test = train_test_split(df, test_size=0.2, shuffle=False, random_state=123)

# Total lags defined
tot_lags = 30
df_train = generate_lags(df=df_train, num_lags=tot_lags)
df_test = generate_lags(df=df_test, num_lags=tot_lags)

X_train, y_train = (
    df_train.drop(columns=['target']),
    df_train['target']
)
X_test, y_test = (
    df_test.drop(columns=['target']),
    df_test['target']
)

# Print shapes of the resulting train and test data sets
print(f'Shape of X_train and y_train: {X_train.shape, y_train.shape}')
print(f'Shape of X_test and y_test: {X_test.shape, y_test.shape}')
preprocessor = Pipeline(
    [
        ('add_lag_step', UniveriateTimeSeriesAddFeatures()),
        ('stdscaler', StandardScaler())
    ]
)
results = {}
dummy = DummyRegressor()

results["Dummy"] = mean_std_cross_val_scores(
    dummy, X_train, y_train, return_train_score=True, scoring=scoring_metrics
)

pd.DataFrame(results)
pipe_ridge = make_pipeline(preprocessor, Ridge(random_state=123))
pipe_xgb = make_pipeline(preprocessor, XGBRegressor(random_state=123))
pipe_rf = make_pipeline(preprocessor, RandomForestRegressor(random_state=123))

models = {
    'ridge': pipe_ridge,
    'xgboost': pipe_xgb,
    'rf_regressor': pipe_rf
}
for (name, model) in models.items():
    results[name] = mean_std_cross_val_scores(
        model, X_train, y_train, return_train_score=True, scoring=scoring_metrics
    )

pd.DataFrame(results)
pipe_xgb.fit(X_train, y_train);

preds = pipe_xgb.predict(X_test)

actual_vs_pred_df = y_test.reset_index()
actual_vs_pred_df['predicted'] = preds
# actual_vs_pred_df.head()

actual_vs_pred_df_long = pd.melt(actual_vs_pred_df, id_vars='Datetime')

# actual_vs_pred_df_long.tail()
# Visualize the actual vs predicted chart for the last 2 days of the test data period
alt.Chart(actual_vs_pred_df_long[actual_vs_pred_df_long.Datetime > '2018-07-24'], 
          title='Actual vs Predicted on Test Data').mark_line().encode(
    x=alt.X('Datetime', type='temporal', title='Datetime'),
    y=alt.Y('value', type='quantitative', title='Actual vs Predicted'),
    color=alt.Color('variable', title='Type')
).properties(
    width=700
).interactive()
print(f'MAPE for test set: {mape(actual_vs_pred_df["target"], actual_vs_pred_df["predicted"]):.2f}%')<jupyter_output>MAPE for test set: 0.89%
{'fraction_non_alphanumeric': 0.12175102599179206, 'fraction_numerical': 0.013679890560875513, 'mean_word_length': 4.699149539333805, 'pattern_counts': {'":': 3, '<': 54, '<?xml version=': 0, '>': 58, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '15', 'dir': 'jupyter-structured-clean-dedup', 'id': '64746', 'license': 'MIT', 'n_tokens_mistral': 5877, 'n_tokens_neox': 5657, 'n_words': 1157, 'path': 'timeseries_using_lags.ipynb', 'repo_name': 'arijc76/stock-lstm', 'size': '16082'}