<jupyter_start><jupyter_text>Feature Selection This Notebook Attempts to Apply the Approach Described in "Short-term load forecasting using a two-stage sarimax model" (Tarsitano & Amerise) with the goal of efficiently producing parsimonious model parameters through backwards stepwise regression.We do not attempt to rigorously hypothesize that short term electricity loads, and bicycle counts behave in an analagous manner that makes the model approach generalizable, but rather explore the intuition that the model selection approach described may be similarly useful and specifically extensible to the bicycle count data, as both represent similar computational issues due to the relatively high frequency seasonal components in both model groups.<jupyter_code>import pmdarima as pm
from pmdarima import arima
from pmdarima import model_selection
from pmdarima import pipeline
from pmdarima import preprocessing as ppc
from pmdarima.arima import ADFTest 


from sklearn.compose import ColumnTransformer
from sklearn.preprocessing import PowerTransformer
from sklearn.metrics import mean_squared_log_error, mean_squared_error

from statsmodels.tsa.deterministic import CalendarSeasonality


import cabi.etl.load as l
import cabi.etl.transform as t

import pandas as pd
import numpy as np
from matplotlib import pyplot as plt

print("pmdarima version: %s" % pm.__version__)


def RMSE(y_true, y_pred):
    return mean_squared_error(y_true, y_pred, squared=False)<jupyter_output><empty_output><jupyter_text>Load the Data Select Top Five Most Active ANCs in Either Direction Follow Up Ideas **FLAGGED FOLLOW UP**- Model Checkins/Checkouts by selecting start/ends from trips long- Model Poisson instead of SARIMA?<jupyter_code>counts = l.load_counts_full()
pd.set_option('display.float_format', lambda x: '%.5f' % x)

# 1A/1C has most outflow, 2E/2B/6D have most inflow on average
counts.mean().sort_values()
bot_five = counts.sum().sort_values().head(5).index
top_five = counts.sum().sort_values().tail(5).index
print(top_five, bot_five)
model_groups = list(bot_five) + list(top_five)
model_groups
bot_five = ['1A', '1C', '3C', '5E', '4C']
top_five = ['6C', '2C', '6D', '2E', '2B']
hourly_groups = counts[model_groups].resample('1H').sum()
hourly_groups = hourly_groups[hourly_groups.index > '2020-06-15']
hourly_groups<jupyter_output><empty_output><jupyter_text>Create Weekday/Hourly Dummies, Weekly Fourier Features to Backwards Eliminate<jupyter_code>def get_seasonal_dummies(df):
    """Accepts a time-indexed df of hourly data, returns hourly and weekday dummies as a df
    to passed as exogenous variables in a SARIMAX model"""
    columns = df.columns
    new_df = df.copy()
    new_df['time'] = new_df.index
    
    # create weekday dummy generator
    wday_dumgen = ppc.DateFeaturizer(column_name='time', with_day_of_month=False)

    # since all have the same index, we can use any column in the df to generate the day_dums
    _, wday_dums = wday_dumgen.fit_transform(new_df[columns[0]], new_df)
    
    # drop the columns that aren't dummies
    wday_dums = wday_dums[wday_dums.columns[-7:]]
    
    # set the index for easy merging
    wday_dums.set_index(new_df.index, inplace=True)
    
    # create hourly dummy generator
    hourly_dumgen = CalendarSeasonality('H', 'D')
    
    # generate dummies
    hourly_dums = hourly_dumgen.in_sample(new_df.index)
    
    # merge results
    full_dums = wday_dums.merge(hourly_dums, on='time')
    
    return full_dums
    
# for use with pmdarima, the timestamps must be in a column instead of the index
hourly_groups['time'] = hourly_groups.index
hourly_groups

wday_dums = ppc.DateFeaturizer(column_name='time', with_day_of_month=False)

# since all have the same index, we can use any column in the df to generate the day_dums
_, day_dums = wday_dums.fit_transform(hourly_groups['1A'], hourly_groups)
day_dums
# drop the columns that aren't dummies
day_dums = day_dums[day_dums.columns[-7:]]
day_dums.set_index(hourly_groups.index, inplace=True)
day_dums.columns
day_dums
hourly_dumgen = CalendarSeasonality('H', 'D')
hourly_dummies = hourly_dumgen.in_sample(hourly_groups.index)<jupyter_output><empty_output><jupyter_text>Note on SparsitySee below representation of each hour of the day for each timestamp in the index.Note the large number of zero values this results in (for each row only one of 24 columns will have a non-zero value).We will first attempt to fit the data in this manner, but if it proves inefficient, it may be worth converting these columns to binary data, as in (23 = 10111, instead of 00000...1)<jupyter_code>hourly_dummies
full_dums = day_dums.merge(hourly_dummies, on='time')
full_dums
hourly_groups = hourly_groups.drop('time', axis=1)
hourly_groups
get_seasonal_dummies(hourly_groups)<jupyter_output><empty_output>