Введение
Временные ряды — это тип последовательных данных, содержащих точки данных, отображаемые в определенный последовательный период времени. Это может варьироваться в единицах времени, от секунд до часов, от дней до недель и от месяцев до лет. Временные ряды используются в разных секторах, чтобы увидеть, как определенная переменная изменяется со временем. Например, наиболее популярно использование временных рядов на финансовых рынках для отслеживания изменений в цене и объеме определенного товара и получения информации о закономерностях или тенденциях этого товара.
Существует два типа анализа временных рядов: одномерный и многомерный. В то время как одномерный анализ временных рядов связан с прогнозированием только одной переменной, многомерный анализ шире и более применим к реальному миру, поскольку он учитывает другие переменные, которые могут повлиять на значение зависимой переменной.
В нашем случае нас интересует мониторинг уровней содержания загрязняющих веществ в воздухе на различных лабораторных станциях в Мадриде. Набор данных был получен из набора данных Качество воздуха на Kaggle в рамках моих стажировок. С этим набором данных было немного сложно работать, так как я никогда раньше не имел дела с анализом временных рядов.
На самом деле я избегал временных рядов из-за их сложной репутации в сообществе специалистов по данным. Однако необходимость снова оказалась матерью изобретения, и мне пришлось завершить этот проект в рамках моих обязанностей на вышеупомянутой стажировке. Это расширило мои знания в области функционального кодирования и заставило меня делать больше, поэтому я благодарен за опыт и крутую кривую обучения.
Цель
Я долго думал над постановкой задачи, которая, как я в конце концов решил, будет заключаться в прогнозировании уровней окиси углерода (CO), поскольку это один из наиболее известных опасных загрязнителей воздуха в мире, с использованием уровней других загрязнителей.
Импорт пакетов
import pandas as pd import glob from sklearn.experimental import enable_iterative_imputer from sklearn.impute import IterativeImputer from sklearn import linear_model import matplotlib.pyplot as plt import seaborn as sns from sktime.utils.plotting import plot_series from sklearn.feature_selection import RFE from sklearn.ensemble import RandomForestRegressor import numpy as np from statsmodels.tsa.stattools import adfuller from numpy import log from sklearn.metrics import mean_absolute_error from sklearn.metrics import mean_squared_error
Предварительная обработка данных
Этот набор данных был разбит на несколько файлов, по одному на каждый год с 2001 по 2018 год. Каждый файл содержал столбец даты, а также столбец для каждого из 17 зарегистрированных загрязняющих веществ. Столбец также был посвящен коду станции.
Первое, что нужно было сделать, это объединить файлы за все годы, чтобы я мог лучше понять свой следующий шаг. Я использовал для этого пакет glob на Python, который используется для возврата всех путей к файлам, соответствующих определенному шаблону в указанном каталоге. Затем я передал «глобальные» файлы в один фрейм данных.
#Using the glob package to read all files with datatype csv contained in the folder
files = glob.glob("C:/Users/uchei/OneDrive/Desktop/csvs_per_year/csvs_per_year**/*.csv")
df = pd.DataFrame()
for f in files:
csv = pd.read_csv(f)
df = df.append(csv)
Давайте посмотрим, как выглядел фрейм данных с помощью базовой команды.
Следует отметить, что одна и та же дата и время появляются в этом выводе несколько раз. Это связано с разными кодами станций, как показано в столбце «Станция». В этом коротком простом фрагменте также содержится довольно много нулевых значений, но мы немного разберемся с ними. Продолжая основное исследование, я заметил, что тип данных столбца «дата» является объектом и его необходимо преобразовать в дату и время, как показано ниже.
#date shouldn't be an object dtype, change to datetime df['date'] = pd.to_datetime(df.date, infer_datetime_format = True)
Затем я должен был убедиться в единообразии измерений уровней загрязнения. Я сделал это на веб-сайте Kaggle, указанном выше, и увидел, что есть несколько несоответствий, которые необходимо устранить.
#Divide columns by 1000 to get desired unit (μgm-3) df['CO'] = df['CO']/1000 df['TCH'] = df['TCH']/1000 df['CH4'] = df['CH4']/1000
Выполнив эту задачу, мне нужно было выяснить, как решить проблему с нулевыми значениями. После небольшого размышления и холостых манипуляций с данными я обнаружил, что после группировки набора данных по станциям на некоторых станциях не было наблюдений для определенных газов. Следовательно, единственный способ обработки этих значений для получения надлежащей информации — после разделения по станциям.
df_stations = df.groupby('station')
casa_de_campo = df_stations.get_group(28079024)
casa_de_campo.index = casa_de_campo.date
casa_de_campo = casa_de_campo.drop(['date'], axis = 1)
casa_de_campo.head()
Я выбрал код станции наугад из 39 уникальных кодов станций, содержащихся во всем наборе данных. Кроме того, я передал столбец даты, чтобы он стал индексом, чтобы упростить визуализацию, которую можно будет обработать позже. Эта станция под названием Casa de Campo содержит около 150 000 наблюдений и 17 столбцов, не считая столбца «Станция», который был одинаковым для всей группы. Индекс даты теперь также имеет смысл последовательно, так как значения записываются каждый час в каждом столбце.
Обработка отсутствующих значений
Я обрабатывал нулевые значения, используя исключение по величине процента отсутствующих значений и асимметрии. В некоторых столбцах просто отсутствовало слишком много значений, чтобы их можно было считать полезными, поскольку я установил контрольный показатель сохранения столбцов на уровне 50 % записанных данных.
casa_de_campo.isnull().sum()/151416 * 100
casa_de_campo = casa_de_campo.drop(['MXY', 'OXY', 'PXY', 'NO', 'CH4', 'station'], axis = 1)
Столбец «Станция» также был удален, поскольку он был бы излишним при моделировании, так как все значения одинаковы. В случае асимметрии столбцам со значением асимметрии 0,5 или ниже отсутствующие значения вменялись в их среднее значение.
casa_de_campo.skew()
Тем, кто выше этого порога асимметрии, приписывались медианы в зависимости от их общего процента отсутствующих данных.
#Skewness under 0.5 casa_de_campo['O_3'].fillna(casa_de_campo['O_3'].mean(), inplace = True) #Skewness above 0.5 casa_de_campo['CO'].fillna(casa_de_campo['CO'].median(), inplace = True) casa_de_campo['NMHC'].fillna(casa_de_campo['NMHC'].median(), inplace = True) casa_de_campo['NO_2'].fillna(casa_de_campo['NO_2'].median(), inplace = True) casa_de_campo['PM10'].fillna(casa_de_campo['PM10'].median(), inplace = True) casa_de_campo['SO_2'].fillna(casa_de_campo['SO_2'].median(), inplace = True) casa_de_campo['TCH'].fillna(casa_de_campo['TCH'].median(), inplace = True)
Третья группа, содержащая более высокий процент пропущенных значений и асимметрию более 4, была импутирована с использованием линейной оценки байесовского хребта в MICE Iterative Imputer при обратном заполнении.
# Define MICE Imputer and fill missing values for all other columns mice_imputer = IterativeImputer(estimator=linear_model.BayesianRidge(), n_nearest_features=None, imputation_order='ascending') campo_imputed = pd.DataFrame(mice_imputer.fit_transform(casa_de_campo), columns=casa_de_campo.columns, index = casa_de_campo.index)
Следующим шагом было определение стационарности нового набора данных, который теперь свободен от нулевых значений.
#Checking stationarity of data
result = adfuller(campo_imputed['CO'].values)
print('p-value: %f' % result[1])
p-value: 0.000000
Значение p появляется очень мало. Поскольку значение ниже 0,05, мы можем предположить, что столбец CO является стационарным, что означает, что данные не нуждаются в дальнейшем преобразовании.
После этого короткого теста, с целью максимизации точности, данные еженедельно передискретизировались. Это было сделано с использованием среднего уровня загрязняющих веществ за неделю из предоставленных данных.
#Get the mean of values in the dataset on a weekly basis
weekly_campo = campo_imputed.resample("W").mean()
weekly_campo.head()
Мы должны помнить, что визуализации являются чрезвычайно важной частью науки о данных, чтобы дать понимание с одного взгляда. Следуя этому принципу, я использовал пакет plot_series, содержащийся в sktime.
plot_series(weekly_campo['CO'], x_label = 'Date', y_label = 'Levls of Carbon Monoxide (μgm-3) ')
Из этого графика видно, что уровни выбросов угарного газа на станции Каса-де-Кампо значительно снизились по сравнению с высокими уровнями, наблюдавшимися в начале 2000-х годов.
Выбор функции
Как всегда, при выборе функций для обучения модели вам нужны независимые переменные, которые действительно влияют на зависимую переменную. Распространенным способом определения независимых переменных для выбора является тепловая карта корреляции. Давайте посмотрим один для этого набора данных.
fig, ax = plt.subplots(figsize=(10,5)) sns.heatmap(weekly_campo.corr(), annot = True, ax = ax) plt.show()
Наша целевая переменная — CO, поэтому мы можем сосредоточиться только на этом столбце и игнорировать остальные. Из этой тепловой карты мы можем сказать, что загрязнитель NMHC почти не имеет корреляции с CO. Кроме того, PM10 (твердые частицы) находится на нижней стороне коэффициента Пирсона. O_3 (озон), несмотря на отрицательную корреляцию, выбран потому, что его влияние значительно на уровни CO (угарный газ).
corr_weekly = weekly_campo.drop(['NMHC', 'PM10'], axis = 1) corr_weekly.head()
Впоследствии, чтобы получить еще более высокую точность модели, я создал новый кадр данных, чтобы ввести переменные задержки. Поскольку предыдущий фрейм данных был передискретизирован на недельном таймфрейме, я выбрал 54 лаговых переменных для целевых данных, CO, для работы с ними, так как это идеальное количество лагов для получения недельных данных. Запаздывающие переменные основаны на прошлых значениях зависимой переменной.
#Creating lag variables
dataframe = pd.DataFrame()
#On a weekly basis, the ideal number of lags should be 54.
for i in range(54, 0, -1):
dataframe['t-' + str(i)] = corr_weekly.CO.shift(i)
#Combine the lagged dataframe with the original
lagged_data = pd.concat([corr_weekly, dataframe], axis = 1)
lagged_data.dropna(inplace = True)
Моделирование
Я вручную разделил данные (да, это можно сделать и с помощью временного разделения) на наборы для обучения и тестирования, используя последние 26 недель набора данных в качестве набора для тестирования.
#Slice the data into train and test #The testing slice is for the final 26 weeks of the dataset train = lagged_data.loc['2002-01-20' : '2017-11-05'] test = lagged_data.loc['2017-11-12' : '2018-05-06'] x_test = test.loc[:, test.columns != 'CO'] y_test = test['CO'] x_train = train.loc[:, train.columns != 'CO'] y_train = train['CO']
Я выбрал случайный лесной регрессор в качестве метода ансамбля, выбрав только 13 наиболее важных функций с помощью рекурсивного исключения признаков (RFE).
print("Training model...")
reg_forest = RandomForestRegressor(n_estimators = 100, random_state = 42)
rfe_weekly = RFE(estimator = reg_forest, n_features_to_select = 13)
fit_weekly = rfe_weekly.fit(x_train, y_train)
predictions1 = fit_weekly.predict(x_test)
Оценка модели
Чтобы оценить производительность этой модели, я использовал MAPE, чтобы получить значение точности в процентах. Также использовалась предвзятость.
correct = np.array(y_test)
sumvalue = np.sum(correct)
mape = np.sum(np.abs((correct - predictions1))) / sumvalue * 100
accuracy = 100 - mape
print('Accuracy:', round(accuracy, 2), '%.')
Accuracy: 90.26 %.
Согласно расчетам с использованием MAPE, наша модель имеет ошибку 9,74% при прогнозировании фактических значений. В прогнозировании временных рядов, согласно Монтано и Палмеру и др., это можно интерпретировать как высокоточный прогноз.
#calculating bias
errors = [correct[i] - predictions1[i] for i in range(len(correct))]
bias = sum(errors) * 1.0/len(correct)
print('Bias: %f' % bias)
Bias: -0.000008
Наше смещение очень низкое, однако отрицательный знак перед значением смещения говорит нам, что наша модель превзошла прогноз. Это означает, что наша модель предсказывает значения, которые выше фактических значений.
Чтобы построить график sktime, содержащий прогнозы, необходимо было выполнить некоторую обработку, поскольку функция plot_series не принимает массивы и требует совпадения типов индексов двух кадров данных.
#As the predictions are in array, we need to convert them to a DataFrame in order for the plot_series function to read it predicted_df = pd.DataFrame(predictions1) predicted_df['date'] = y_test.index predicted_df['date'] = pd.to_datetime(predicted_df.date, infer_datetime_format = True) predicted_df.index = predicted_df['date'] predicted_df = predicted_df.drop(['date'], axis = 1) predicted_df.index = pd.to_datetime(predicted_df.index)
Теперь мы можем успешно использовать функцию plot_series.
plot_series(y_train, y_test, predicted_df, x_label = 'Date', y_label = 'CO Levels (μgm-3)', labels = ['Weeks', 'Actual CO Levels (μgm-3)', 'Predicted Levels (μgm-3)']
Хотя у нас есть показатель точности, этот график выглядит немного гранулированным, и его трудно разобрать из-за большого количества содержащихся точек данных. Может быть трудно рассказать историю наших открытий, если визуализация не передает ясности.
Чтобы решить эту проблему, мы можем использовать менее детальный метод повторной выборки. Если данные пересчитываются ежемесячно, график становится более четким, но за счет точности, которая падает до 86,66%, как показано ниже.
monthly_campo = campo_imputed.resample('M').mean()
relevant_monthly = monthly_campo.drop(['NMHC', 'PM10'], axis = 1)
#Creating lag variables
dataframe = pd.DataFrame()
for i in range(12, 0, -1):
dataframe['t-' + str(i)] = relevant_monthly.CO.shift(i)
lagged_monthly = pd.concat([relevant_monthly, dataframe], axis = 1)
lagged_monthly.dropna(inplace = True)
#Adding seasonal variables
train1 = lagged_monthly.loc['2002-01-31' : '2017-11-30']
test1 = lagged_monthly.loc['2017-12-31' : '2018-05-31']
x_test1 = test1.loc[:, test1.columns != 'CO']
y_test1 = test1['CO']
x_train1 = train1.loc[:, train1.columns != 'CO']
y_train1 = train1['CO']
print("Training model...")
monthly_reg = RandomForestRegressor(n_estimators = 100, random_state = 42)
monthly_rfe = RFE(estimator = monthly_reg, n_features_to_select = 9)
monthly_fit = monthly_rfe.fit(x_train1, y_train1)
monthly_predictions = monthly_fit.predict(x_train1)
Создание распознаваемого формата для plot_series.
df1 = pd.DataFrame(monthly_predictions) df1['date'] = y_test1.index df1['date'] = pd.to_datetime(df1.date, infer_datetime_format = True) df1.index = df1['date'] df1 = df1.drop(['date'], axis = 1) df1.index = pd.to_datetime(df1.index) plot_series(y_train1, y_test1, df1, labels = ['Months', 'Actual CO Levels (μgm-3)', 'Predicted Levels'])
На этом графике легче увидеть, как модель укладывается в график, чем на недельных данных, хотя он и менее точен. Прогнозируемые значения падают почти одновременно с фактическими уровнями. Действительно аккуратно!
Заключительные мысли
Много раз я рассматривал возможность проведения базового одномерного анализа временных рядов, особенно когда сталкивался с ошибками (их было так много на этом пути), но я бросал себе вызов, чтобы раздвинуть свои пределы и посмотреть, как далеко я могу зайти. Вы можете найти мой репозиторий GitHub здесь.
Спасибо за прочтение!
