<jupyter_start><jupyter_text>US Road Accidents Forecasting Models for Top 10 States Import libraries and Data<jupyter_code>from datetime import datetime

import pandas as pd
import numpy as np
import seaborn as sns
from colorama import Fore
import seaborn as sn

%matplotlib inline
import matplotlib.pyplot as plt
from matplotlib import rcParams
df = pd.read_csv("US_Accidents_Dec20.csv")
df['Date'] = pd.to_datetime(df['Start_Time']).dt.date #short date
df['Date'] = pd.to_datetime(df['Date']) #convert back to date
#filter for California state to build daily forecast models for each city
df = df.loc[df.State.isin(['CA','TX','FL','SC','NC','NY','PA','MI','IL', 'GA'])]
#select columns of interest
df = df[['Date', 'State', 'Temperature(F)', 'Pressure(in)','Visibility(mi)', 'Weather_Condition']]
#roll up the accidents to day level
df = df.groupby(['Date', 'State']).agg({'Temperature(F)' : ['mean', 'count'], \
                                                  'Pressure(in)' : 'mean', 'Visibility(mi)' : 'mean', \
                                                  'Weather_Condition': 'first'}).reset_index()
#convert multi level index to single level index
df.columns = ['_'.join(col) for col in df.columns.values]
df.rename(columns={'Temperature(F)_count':'DailyAccidents', 'Temperature(F)_mean':'Temperature', \
                  'Pressure(in)_mean':'Pressure', 'Visibility(mi)_mean':'Visibility', \
                  'Weather_Condition_first':'Weather_Condition', \
                  'Date_':'Date', \
                  'State_':'State', \
                  }, inplace=True)
df.head(5)<jupyter_output><empty_output><jupyter_text>Plot Road Accidents time series data by States (Top 10)Data exibits time series characteristics such daily/weekly/yearly seasonality and potentially a good candidate for time-series modeling:<jupyter_code>plt.rcParams["figure.figsize"] = (12,7)
sn.lineplot(data = df[['Date','State', 'DailyAccidents']].loc[df.Date<='2017-01-01'], \
            x='Date', y='DailyAccidents', hue='State').set_title('US Road Accidents over time for top 10 States')<jupyter_output><empty_output><jupyter_text>Split train test for modeling<jupyter_code>train_start_date, train_end_date = '2016-02-09', '2019-01-31'
test_start_date, test_end_date = '2019-02-01', '2019-06-30'

train_df = df.loc[(df.Date>= train_start_date) & (df.Date<= train_end_date)]
train_df = train_df.fillna(method='ffill')
test_df = df.loc[(df.Date>= test_start_date) & (df.Date<= test_end_date)]
test_df = test_df.fillna(method='ffill')<jupyter_output><empty_output><jupyter_text>Generalized Additive Model (Prophet Model)<jupyter_code>from fbprophet import Prophet
from fbprophet.diagnostics import cross_validation, performance_metrics
from fbprophet.plot import plot_cross_validation_metric, add_changepoints_to_plot
import warnings
warnings.filterwarnings ("ignore") #ignore warnings

UniqueStates = df.State.unique()
features = ["Temperature", "Pressure", "Visibility"]#features
test_size=(pd.to_datetime(test_end_date) - pd.to_datetime(train_end_date)).days

dfcolumns = ["Date", "State", "DailyAccidents", "Forecast"]
forecast_all_df = pd.DataFrame(columns= dfcolumns)

UniqueStates = ['CA']
#loop through each state to build time series models
for state in UniqueStates:
    print(state)
    
    localdf = train_df.loc[train_df.State == state].copy()#train df local copy
    dfm  = localdf[["Date", "DailyAccidents"]].copy()
    dfm = dfm.sort_values("Date")
    dfm.columns = ['ds','y'] #required for calling in prophet
    dfm['ds'] = pd.to_datetime(dfm['ds'], format="%Y-%m-%d").dt.date
    #add features
    for f in features:
        dfm[f] = localdf[[f]]        
        
    #instantiate model object
    m = Prophet(daily_seasonality=False, yearly_seasonality=True, weekly_seasonality=True, mcmc_samples=500, n_changepoints=0)
    #for adding regressors to model
    for f in features:
        m.add_regressor(f)
    #fit
    m.fit(dfm) 
    #prepare for predict function
    future = m.make_future_dataframe(periods=test_size, freq='D')
    #add features for predict dataframe
    for f in features:
        future[f] = np.array(df.loc[(df.State == state) & (df.Date >= train_start_date) & (df.Date <= test_end_date)].fillna(method='ffill')[[f]])
        
    forecast = m.predict(future)
    forecast_df = forecast #do not change original df from prophet
    forecast_df['DailyAccidents'] = np.array(df.loc[(df.State == state) & (df.Date >= train_start_date) & (df.Date <= test_end_date)][['DailyAccidents']])
    forecast_df['State'] = state
    forecast_df.rename(columns= {'ds': 'Date', 'yhat': 'Forecast'}, inplace=True)
    forecast_df = forecast_df[dfcolumns]
    forecast_df = forecast_df.loc[(forecast_df.Date>= test_start_date) & (forecast_df.Date<= test_end_date)]
    forecast_all_df = pd.concat([forecast_all_df , forecast_df])
    del forecast, forecast_df<jupyter_output>Importing plotly failed. Interactive plots will not work.<jupyter_text>Plot the Forecast, calcuate RMSE and show forecast components:1. Models are able to pick up weekly seasonality pretty well2. The effect of external signals such as Pressure, Temperature are significant3. Model out of sample RMSE is shown below for each state<jupyter_code>from sklearn.metrics import mean_squared_error
rmse_list =[]
for state in UniqueStates:
    testdf = forecast_all_df.loc[forecast_all_df.State==state]
    rmse = np.round(mean_squared_error(testdf.DailyAccidents, testdf.Forecast, squared=False),2)
    mean_accidents_actual = np.median(testdf[["DailyAccidents"]])
    rmse_list.append([state, rmse, 'Prophet', mean_accidents_actual])

prophet_rmse = pd.DataFrame (rmse_list, columns=['State', 'RSME (test)', 'Model', 'ActualDailyAccidentsMedian'])
prophet_rmse.to_csv('Prophet_RMSE.csv')
prophet_rmse.sort_values('ActualDailyAccidentsMedian', ascending=False)
#plot california Forecast
forecast_all_df['DailyAccidents'] = forecast_all_df['DailyAccidents'].astype(int)
forecast_all_df['Forecast'] = forecast_all_df['Forecast'].astype(int)
plt.rcParams["figure.figsize"] = (12,7)
sn.lineplot(data = forecast_all_df.loc[forecast_all_df.State=="CA"], \
            x='Date', y='DailyAccidents')
sn.lineplot(data = forecast_all_df.loc[forecast_all_df.State=="CA"], \
            x='Date', y='Forecast').set_title('California Road Accident Forecast (Out of sample)');<jupyter_output><empty_output><jupyter_text>Seasonal ARIMA model<jupyter_code>from pmdarima import auto_arima
import warnings
warnings.filterwarnings ("ignore") #ignore warnings

UniqueStates = df.State.unique()
features = ["Temperature", "Pressure", "Visibility"]#features
test_size=(pd.to_datetime(test_end_date) - pd.to_datetime(train_end_date)).days

dfcolumns = ["Date", "State", "DailyAccidents", "Forecast"]
forecast_all_df = pd.DataFrame(columns= dfcolumns)

#UniqueStates = ['CA']
#loop through each state to build time series models
for state in UniqueStates:
    print(state)
    
    localdf = train_df.loc[train_df.State == state].copy()#train df local copy
    
    #prepare y's
    y_to_train = train_df[['Date', 'DailyAccidents']]
    y_to_train.set_index('Date', inplace=True)
    y_to_test = test_df[['Date', 'DailyAccidents']]
    y_to_test.set_index('Date', inplace=True)
    
     # Fit model
    arima_model = auto_arima(y=y_to_train, seasonal=True, m=7, stepwise =True, d=1, \
                         start_p=1, max_p=2 , \
                         start_q =0, max_q=2 , \
                         start_d =1, max_d=1 , \
                         start_P =0, max_P=2 , \
                         start_Q =0, max_Q=2 , \
                         start_D =1, max_D=1 , \
                         D=1 )
    #forecast
    arima_forecast = arima_model.predict(n_periods=(pd.to_datetime(test_end_date)-pd.to_datetime(test_start_date)).days+1, \
                                         return_conf_int=True, alpha = .1)
    

    forecast_df = test_df.loc[test_df.State==state][['Date', 'State', 'DailyAccidents']].copy()
    forecast_df['Forecast'] = arima_forecast[0]    
    forecast_all_df = pd.concat([forecast_all_df , forecast_df])
    del forecast_df<jupyter_output>MI
PA
CA
NY
IL
FL
SC
TX
GA
NC<jupyter_text>Plot the Forecast & calcuate RMSE:1. Models are able to pick up weekly seasonality pretty well2. Model out of sample RMSE is shown below for each state<jupyter_code>#plot california Forecast
forecast_all_df['DailyAccidents'] = forecast_all_df['DailyAccidents'].astype(int)
forecast_all_df['Forecast'] = forecast_all_df['Forecast'].astype(int)
plt.rcParams["figure.figsize"] = (20,5)
sn.lineplot(data = forecast_all_df.loc[forecast_all_df.State=="CA"], \
            x='Date', y='DailyAccidents')
sn.lineplot(data = forecast_all_df.loc[forecast_all_df.State=="CA"], \
            x='Date', y='Forecast').set_title('ARIMA: California Road Accident Forecast (Out of sample)');
from sklearn.metrics import mean_squared_error
rmse_list =[]
for state in UniqueStates:
    testdf = forecast_all_df.loc[forecast_all_df.State==state].copy()
    rmse = np.round(mean_squared_error(testdf.DailyAccidents, testdf.Forecast, squared=False),2)
    mean_accidents_actual = np.median(testdf[["DailyAccidents"]])
    rmse_list.append([state, rmse, 'SARIMA', mean_accidents_actual])

prophet_rmse = pd.DataFrame (rmse_list, columns=['State', 'RSME (test)', 'Model', 'ActualDailyAccidentsMedian'])
prophet_rmse.to_csv('SARIMA_RMSE.csv')
prophet_rmse.sort_values('ActualDailyAccidentsMedian', ascending=False)<jupyter_output><empty_output><jupyter_text>RNN-LSTM Models First create utility functions for running LSTM model:<jupyter_code>from sklearn.preprocessing import OneHotEncoder
from sklearn.preprocessing import MinMaxScaler
 
train_start_date, train_end_date = '2016-02-09', '2018-12-31'
test_start_date, test_end_date = '2019-01-01', '2019-01-31'
val_start_date, val_end_date = '2019-02-01', '2019-06-30'#keep same as kept for ARIMA and prophet

#create static list of feature names
CAT_COLS = ['State']
NUM_COLS = ['Temperature', 'Visibility', 'DailyAccidents_Weekly', 'DailyAccidents_Yearly'] 
 
def create_supervise_dataset(df, ohe, minmaxscaler):
    train_df = df.loc[(df.Date>= train_start_date) & (df.Date <= train_end_date)]
    X_train = np.hstack((ohe.transform(train_df[CAT_COLS]), train_df[NUM_COLS]))
    X_train = minmaxscaler.fit(X_train).transform(X_train)
    Y_train = train_df[['DailyAccidents']]
  
    test_df = df.loc[(df.Date>= test_start_date) & (df.Date <= test_end_date)]
    X_test = np.hstack((ohe.transform(test_df[CAT_COLS]), test_df[NUM_COLS]))
    X_test = minmaxscaler.transform(X_test)
    Y_test = test_df[['DailyAccidents']]
  
    val_df = df.loc[(df.Date>= val_start_date) & (df.Date <= val_end_date)]
    X_val = np.hstack((ohe.transform(val_df[CAT_COLS]), val_df[NUM_COLS]))
    X_val = minmaxscaler.transform(X_val)
    Y_val = val_df[['DailyAccidents']]
  
    return X_train, X_test, X_val, Y_train, Y_test, Y_val<jupyter_output><empty_output><jupyter_text>Add seasonality features explicitly in the model:<jupyter_code>TimeSeriesValues = df.State.unique()
 
lead_df_cols = ["Date_Yearly", "Date_Weekly", "DailyAccidents_Yearly", "DailyAccidents_Weekly", "State"] 
lead_df = pd.DataFrame(columns=lead_df_cols)
 
for timeseries in TimeSeriesValues:    
    localdf = df.loc[df.State==timeseries][["Date", "DailyAccidents", "State"]]
    localdf['DailyAccidents_Yearly'] = localdf['DailyAccidents']
    localdf['DailyAccidents_Weekly'] = localdf['DailyAccidents']    
    localdf["Date_Yearly"] = localdf.Date + pd.DateOffset(days=365)
    localdf["Date_Weekly"] = localdf.Date + pd.DateOffset(days=7)
    localdf.drop(['Date', 'DailyAccidents' ], inplace=True, axis=1)
    lead_df = pd.concat([lead_df,localdf])

lead_df.reset_index().drop('index',  inplace=True, axis=1)

df = (pd.merge(df, lead_df[['Date_Weekly', 'State', 'DailyAccidents_Weekly' ]], how = 'left', left_on=['Date', 'State'], right_on=['Date_Weekly', 'State']))
df = (pd.merge(df, lead_df[['Date_Yearly', 'State', 'DailyAccidents_Yearly' ]], how = 'left', left_on=['Date', 'State'], right_on=['Date_Yearly', 'State']))
df.drop(['Date_Weekly', 'Date_Yearly'], inplace=True, axis=1)
df.tail(20)<jupyter_output><empty_output><jupyter_text>Train a single LSTM model for all States:<jupyter_code>from tensorflow.keras.models import Sequential
from tensorflow.keras.layers import Dense
from tensorflow.keras.layers import LSTM
from keras.callbacks import EarlyStopping
from matplotlib import pyplot

df.dropna(inplace=True)
df_copy = df.copy()#remove records adding yearly seasonality

#one hot encoding based on train set
ohe = OneHotEncoder(sparse=False, drop='first').fit(df_copy.loc[df_copy.Date <= train_end_date][CAT_COLS])
encoded_X = ohe.transform(df_copy[CAT_COLS]) #later use for decoding

#normalize features
minmaxscaler_X = MinMaxScaler(feature_range=(0, 1))
minmaxscaler_Y = MinMaxScaler(feature_range=(0, 1))
 
#train test prepare
X_train, X_test, X_val, Y_train, Y_test, Y_val = create_supervise_dataset(df_copy, ohe, minmaxscaler_X)
Y_train= minmaxscaler_Y.fit_transform(Y_train)
Y_test = minmaxscaler_Y.transform(Y_test)
Y_val = minmaxscaler_Y.transform(Y_val)
 
print('shape before reshape')
print(X_train.shape, X_test.shape, X_val.shape)
print(Y_train.shape, Y_test.shape, Y_val.shape)
 
# reshape input to be 3D [samples, timesteps, features]
print('shape after reshape')
X_train = X_train.reshape((X_train.shape[0], 1, X_train.shape[1]))
X_test = X_test.reshape((X_test.shape[0], 1, X_test.shape[1]))
X_val = X_val.reshape((X_val.shape[0], 1, X_val.shape[1]))
print(X_train.shape, X_test.shape, X_val.shape)
 
# design network
model = Sequential()
model.add(LSTM(25, input_shape=(X_train.shape[1], X_train.shape[2]), dropout=.1))
model.add(Dense(1))
model.compile(loss='mse', optimizer='adam')
#early stopping
es = EarlyStopping(monitor='val_loss', mode='min', verbose=1, patience=25)
# fit network
history = model.fit(X_train, Y_train, epochs=100, batch_size=64, validation_data=(X_test, Y_test), verbose=2, shuffle=False, callbacks=[es]) #, stateful=True
# plot history
pyplot.plot(history.history['loss'], label='train')
pyplot.plot(history.history['val_loss'], label='test')
pyplot.legend()
pyplot.show()
X_train, X_test, X_val, Y_train, Y_test, Y_val = create_supervise_dataset(df_copy, ohe, minmaxscaler_X)
X_val = X_val.reshape((X_val.shape[0], 1, X_val.shape[1]))
X_test_Forecast = pd.DataFrame(minmaxscaler_Y.inverse_transform(((model.predict(X_val)))), columns=['Forecast']).round(0)

forecast_all_df = df.loc[(df.Date>= val_start_date) & (df.Date <= val_end_date)].reset_index()
forecast_all_df['Forecast'] = X_test_Forecast['Forecast']
forecast_all_df.head(5)
from sklearn.metrics import mean_squared_error
rmse_list =[]
for state in df.State.unique():
    testdf = forecast_all_df.loc[forecast_all_df.State==state].copy()
    rmse = np.round(mean_squared_error(testdf.DailyAccidents, testdf.Forecast, squared=False),2)
    mean_accidents_actual = np.median(testdf[["DailyAccidents"]])
    rmse_list.append([state, rmse, 'LSTM_1', mean_accidents_actual])

lstm_rmse = pd.DataFrame (rmse_list, columns=['State', 'RSME (test)', 'Model', 'ActualDailyAccidentsMedian'])
lstm_rmse.sort_values('ActualDailyAccidentsMedian', ascending=False).to_csv('LSTM_4.csv')
lstm_rmse.sort_values('ActualDailyAccidentsMedian', ascending=False)
test = forecast_all_df
 
lineStart = test.DailyAccidents.min() 
lineEnd = test.DailyAccidents.max()  
 
plt.close("all")
plt.rcParams["figure.figsize"] = [8,8]
 
plt.scatter(test['DailyAccidents'],test['Forecast'])
plt.plot([lineStart, lineEnd], [lineStart, lineEnd], 'k-', color = 'r')
plt.xlim(lineStart, lineEnd)
plt.ylim(lineStart, lineEnd)<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.14085672458786838, 'fraction_numerical': 0.01164789001336643, 'mean_word_length': 5.353416902547513, 'pattern_counts': {'":': 0, '<': 50, '<?xml version=': 0, '>': 48, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 6, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '10', 'dir': 'jupyter-structured-clean-dedup', 'id': '184314', 'license': 'MIT', 'n_tokens_mistral': 5718, 'n_tokens_neox': 5461, 'n_words': 1042, 'path': 'U.S Road Accident Forecasting.ipynb', 'repo_name': 'Arora-Dhruv/U.S-Road-Accidents-Forecasting-using-LSTM', 'size': '15711'}