Введение
Для этого проекта я решил проанализировать и построить модель временных рядов, предсказывающую количество преступлений в ближайшие 5 лет, разделенных по районам. Я получил набор данных с веб-сайта полиции Нью-Йорка со всеми типами информации о преступлениях с 2006 по 2019 год. После очистки описаний преступлений, мест, типов данных и нулевых значений я получил около 5 миллионов инцидентов с начала 2006 года по конец 2019 года. конец 2019 года. Затем данные были сгруппированы по месяцам и районам. Каждый набор районов был разделен на обучающий набор с 2006 по 2016 год и тестовый набор с 2017 по 2019 год.
Чтобы создать наилучшую модель, я решил использовать среднеквадратичную ошибку (RMSE), чтобы показать, насколько неточными будут мои прогнозы. Эта метрика была зафиксирована для каждой модели, а окончательные результаты были получены в конце. Цель этого поста — показать краткий обзор экспоненциального сглаживания Холта-Уинтерса, SARIMAX и Facebook Prophet. Вы можете проверить мой репозиторий здесь для получения более подробной информации о том, как я очищал и исследовал данные, разлагал временные ряды и проверял компонент стационарности.
САРИМАКС
Прошлые временные точки данных временных рядов могут влиять на текущие и будущие временные точки. Модели ARIMA учитывают эту концепцию при прогнозировании текущих и будущих значений. ARIMA использует ряд запаздывающих наблюдений временных рядов для прогнозирования наблюдений. Вес применяется к каждому прошедшему термину, и веса могут варьироваться в зависимости от того, насколько они недавние. SARIMA(X) — это расширение ARIMA, которое явно поддерживает одномерные данные временных рядов с сезонным компонентом. Он фиксирует скользящие средние по временным рядам, принимая во внимание тенденцию данных, сезонность и шум.
Ниже приведена функция для поиска параметров для заказа SARIMAX и Season_order.
def best_parameters(ts):
'''Enter time series. This function will let auto_arima search for hyperparameters, then SARIMAX will fit these
hyperparameters to get the best model with lowest AIC'''
best_orders = pm.auto_arima(ts, start_p = 0, start_q = 0, max_p = 2, max_q = 2, m = 12, seasonal = True,
stationary = False, stepwise = True, trend = 'ct', suppress_warnings = True,
trace = False, error_action = 'ignore')
best_model = SARIMAX(ts, order = best_orders.order, seasonal_order = best_orders.seasonal_order).fit()
best_parameters = []
best_parameters.append([best_orders.order, best_orders.seasonal_order, best_model.aic])
print('ARIMA {} x {}, AIC Calculated: {}'.format(best_orders.order, best_orders.seasonal_order, best_model.aic))
print(best_model.summary())
return best_modelСледующая функция получает прогноз на один шаг вперед и сравнивает его с фактическими точками данных. Кроме того, он также делает будущие прогнозы для количества шагов, переданных в качестве параметров, отображает результаты с доверительным интервалом 95% и вычисляет RMSE и процентное изменение.
def get_predictions(ts, model, steps = 84, plot=True, show=True): '''Parameters: time series dataframe, model, steps, plot, and show.''' # Get preditions from model for the last 3 years of data period pred = model.get_prediction(start='2017-01-31', dynamic=False) conf = pred.conf_int()# Plot observed and predicted values with confidence interval if plot: ax = ts['2006':].plot(label='Observed', figsize=(10, 8)) pred.predicted_mean.plot(ax=ax, label='One-step-ahead Forecast', alpha=.5) ax.fill_between(conf.index, conf.iloc[:, 0], conf.iloc[:, 1], color='g', alpha=.3, label='Confidence Interval') ax.set_ylabel('Value') ax.set_xlabel('Year') plt.title('Observations vs Predictions') ax.legend() plt.show() # Compare real and predicted values to validade model and compute the rmse predicted = pred.predicted_mean real = ts['2017-01-31':].Crime_Number mse = np.square(np.subtract(real,predicted)).mean() rmse = math.sqrt(mse) # Get forecast and confidence interval for steps ahead in future future = model.get_forecast(steps=steps, dynamic=True) future_conf = future.conf_int()# Plot future forecast with confidence interval if plot: ax = ts['2006':].plot(label='Observed', figsize=(10, 8)) future.predicted_mean.plot(ax=ax, label='Forecast') ax.fill_between(future_conf.index, future_conf.iloc[:, 0], future_conf.iloc[:, 1], color='g', alpha=.3) ax.fill_betweenx(ax.get_ylim(), pd.to_datetime('2026-12-31'), predicted.index[-1], alpha=.1, zorder=-1) ax.set_xlabel('Year') ax.set_ylabel('Number of Crimes') plt.title('Future Forecast') ax.legend() plt.show() # show prediction for end of step-period forecast = future.predicted_mean[-1] upper = future_conf.iloc[-1,1] lower = future_conf.iloc[-1,0] predictions = {} predictions['Upper Bound'] = upper predictions['Expected Forecast'] = forecast predictions['Lower Bound'] = lower predictions = pd.DataFrame.from_dict(predictions, orient='index', columns=['Prediction']) # calculate return percentages crime_2019 = ts.loc['2019-12-31'] forecast_2026 = forecast forecast_lower = lower forecast_upper = upper return_percentage_predictions = {} predicted_percent_change = ((forecast_2026- crime_2019) / crime_2019)*100 upper_percent_change = ((forecast_upper - crime_2019) / crime_2019)*100 lower_percent_change = ((forecast_lower - crime_2019) / crime_2019)*100 return_percentage_predictions['Predicted % Change'] = predicted_percent_change return_percentage_predictions['Upper % Change'] = upper_percent_change return_percentage_predictions['Lower % Change'] = lower_percent_change return_percentage_predictions = pd.DataFrame.from_dict(return_percentage_predictions,orient='index') if show: print(predictions) print('\n' + f'The RMSE of our forecast is {round(rmse, 2)}' + '\n') print(return_percentage_predictions)
Вот сюжет предсказания преступности на Манхэттене.
Фейсбук пророк
Facebook Prophet — это процедура прогнозирования данных временных рядов, основанная на аддитивной модели, в которой нелинейные тренды соответствуют годовой, еженедельной и ежедневной сезонности, а также праздничным эффектам. Он лучше всего работает с временными рядами, которые имеют сильные сезонные эффекты и несколько сезонов исторических данных. Пророк устойчив к отсутствующим данным и сдвигам в тренде и обычно хорошо справляется с выбросами.
Вот функция, которую я создал, чтобы получить прогноз с доверительным интервалом 95%, подогнать модель, построить прогнозы и вычислить RMSE и процентное изменение. Пророк также отображает тенденции, годовые, еженедельные и ежедневные компоненты сезонности.
def prophet_model(ts): '''Input time series data to get prediction with 95% confidence interval, fit the model, plot predictions, calculate percentage change and RMSE'''# train-test split train_data = ts.iloc[:len(ts)-36] test_data = ts.iloc[len(ts)-36:] # Set the uncertainty interval to 95% and yearly seasonality model = Prophet(interval_width=0.95, yearly_seasonality = True, weekly_seasonality=True, daily_seasonality=True) # Fit the timeseries to model model.fit(train_data) # Use make_future_dataframe() with a monthly frequency and periods = 84 future = model.make_future_dataframe(periods=120, freq='M') # Predict the values for future dates and take the head of forecast forecast = model.predict(future) # Use Prophet's plot method to plot the predictions model.plot(forecast, uncertainty=True) plt.show() # Plot model components model.plot_components(forecast) plt.show() # Calculate RMSE y_true = ts['y'].values y_pred = forecast['yhat'][:168].values RMSE = np.sqrt(mean_squared_error(y_true,y_pred)) # show prediction for end of step-period forecast_2026 = forecast['yhat'][251] forecast_lower = forecast['yhat_lower'][251] forecast_upper = forecast['yhat_upper'][251] print(f'Predicted Number of Crimes: {round(forecast_2026, 2)}' + '\n') print(f'Upper Number of Crimes: {round(forecast_upper, 2)}' + '\n') print(f'Lower Number of Crimes: {round(forecast_lower, 2)}' + '\n') # calculate return percentages crime_2019 = ts['y'][167] forecast_2026 = forecast['yhat'][251] forecast_lower = forecast['yhat_lower'][251] forecast_upper = forecast['yhat_upper'][251] return_percentage_predictions = {} predicted_percent_change = ((forecast_2026- crime_2019) / crime_2019)*100 upper_percent_change = ((forecast_upper - crime_2019) / crime_2019)*100 lower_percent_change = ((forecast_lower - crime_2019) / crime_2019)*100 return_percentage_predictions['Predicted % Change'] = predicted_percent_change return_percentage_predictions['Upper % Change'] = upper_percent_change return_percentage_predictions['Lower % Change'] = lower_percent_change print(f'The RMSE of our forecast is {round(RMSE, 2)}' + '\n') print(f'Predicted % Change: {round(predicted_percent_change, 2)}' + '\n') print(f'Upper % Change: {round(upper_percent_change, 2)}' + '\n') print(f'Lower % Change: {round(lower_percent_change, 2)}' + '\n')
График предсказания Манхэттена должен выглядеть примерно так, как показано ниже.
Экспоненциальное сглаживание Холта-Уинтерса
Экспоненциальное сглаживание назначает экспоненциально уменьшающиеся веса и значения по сравнению с историческими данными, чтобы уменьшить значение веса для более старых данных. Он используется для прогнозирования данных временных рядов, которые демонстрируют как тенденцию, так и сезонные колебания.
Подобно моделям выше, я также создал функцию для обучения и тестирования временных рядов. Он также отображает прошлые данные и прогнозы на будущее, а также вычисляет RMSE и процентное изменение.
def hwes (ts, trend, seasonal): '''Enter time series data frame, trend (mul or add), and seasonal (mul or add) to train and test Holt-Winters Exponential Smoothing. This function also plots future predictions and calculate percentage change & RMSE''' # split between the training and the test data sets. The last 36 periods form the test data ts_train = ts.iloc[:-36] ts_test = ts.iloc[-36:] #build and train the model on the training data hwmodel = ExponentialSmoothing(ts_train,trend=trend, seasonal=seasonal, seasonal_periods=12).fit()#create an out of sample forcast for the next 120 steps beyond the final data point in the training data set test_pred = hwmodel.forecast(steps=120) test_pred.to_frame()#plot the training data, the test data and the forecast on the same plot ts_train['Crime_Number'].plot(legend=True, label='Train', figsize=(10,6)) ts_test['Crime_Number'].plot(legend=True, label='Test') test_pred.plot(legend=True, label='Predicted') # calculate RMSE RMSE = np.sqrt(mean_squared_error(ts_test,test_pred[:36])) print('\n' + f'The RMSE of our forecast is {round(RMSE, 2)}' + '\n')# Number of crime by end of 2026 crime_2026 = test_pred.iloc[-1] print(f'Predicted Number of Crimes by End of 2026 is {round(crime_2026, 2)}' + '\n')# Percent change crime_2019 = ts_test['Crime_Number'][-1] return_percentage_predictions = {} predicted_percent_change = ((crime_2026- crime_2019) / crime_2019)*100 return_percentage_predictions['Predicted % Change'] = predicted_percent_change print(f'Predicted % Change: {round(predicted_percent_change, 2)}' + '\n')
А вот и график предсказания количества преступлений на Манхэттене.
Заключение
Для этого конкретного проекта прогнозирования я обнаружил, что SARIMAX вернул самое низкое среднеквадратичное отклонение. Однако я думаю, что у каждого метода есть свои преимущества. SARIMAX намного лучше справляется с сезонностью и компонентами, но поиск гиперпараметров и лучшей модели занимает гораздо больше времени. Экспоненциальное сглаживание Холта-Винтерса — самая быстрая модель для обучения. А Facebook Prophet четко показывает компоненты сезонности (годовые, еженедельные и ежедневные) и «выбросы».
Попробуйте эти модели и дайте мне знать, какая вам нравится :)
