Прогнозирование временных рядов для фондовых рынков

  • Часть 1
  • Прогнозирование цен и тенденций Dow Jones с помощью evoML
  • Прогнозирование временных рядов с помощью Prophet и API
  • Прогнозирование временных рядов с использованием LSTM
  • Transformers vs. LSTM для прогнозирования временных рядов цен на акции
  • Прогнозирование временных рядов цен с помощью Python
  • Прогнозирование временных рядов с помощью PyCaret
  • Прогнозирование времени и рядов с помощью рекуррентной нейронной сети LSTM
  • Исследовательский анализ данных
  • Скользящая средняя, экспоненциальное сглаживание и SARIMA
  • Прогнозирование временных рядов с помощью AutoGluon
  • Преодоление ограничений древовидных моделей при прогнозировании временных рядов
  • Раскрытие возможностей VAR-моделирования для динамического прогнозирования временных рядов
  • Прогнозирование временных рядов с помощью трансформеров

Часть 1

Всегда полезно знать о будущем, не так ли? Особенно, если речь идет о зарабатывании денег?. Принятие решений на основе данных было бы разумнее в мире с большим количеством данных и достижениями в области методов машинного обучения.

В этой статье мы покажем, как использовать сети долговременной и краткосрочной памяти (LSTM) для прогнозирования будущей цены двух фондовых рынков. Новый подход, который мы собираемся обсудить в этой статье, заключается в том, что одна сеть обучается более чем на одних биржевых данных.

Для этого мы будем использовать данные по ценам на акции Apple и Google с 2010 по 2018 год. Данные включают максимум, минимум, цены открытия, закрытия и скорректированные цены закрытия. Он также включает в себя общий объем рынка за период. Таким образом, у нас есть шесть временных рядов для каждой компании, и мы собираемся обучать нашу сеть на основе полученных данных.

Примечание: Этот пост написан только в образовательных целях; Вы не должны рассматривать это как инвестиционный совет.

Загрузка данных

Для начала импортируем необходимые библиотеки.

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

Теперь давайте посмотрим на данные, которые мы получили! Мы построим график цены закрытия для акций Google и Apple.

df_google = pd.read_csv("./GOOG.csv",sep="\t")
df_apple = pd.read_csv("./AAPL.csv",sep="\t")
plt.plot(df_google['Close'])
plt.title("Google Stock Prices")
plt.xlabel("Time (oldest -> latest)")
plt.ylabel("Stock Closing Price")
plt.show()

plt.plot(df_apple['Close'])
plt.title("Apple Stock Prices")
plt.xlabel("Time (oldest -> latest)")
plt.ylabel("Stock Closing Price")
plt.show()
Цены закрытия акций Google и Apple за период.

Как мы видим, диапазон цен закрытия для акций Google и Apple имеет разные масштабы. Поэтому обучение сети становится немного сложнее.
Подготовим данные обучения и тестирования.

from sklearn.preprocessing import MinMaxScaler
from sklearn.model_selection import train_test_split
from sklearn import preprocessing

goog = df_google[['High', 'Low', 'Open', 'Volume','Close','Adj Close']].to_numpy()
aapl = df_apple[['High', 'Low', 'Open', 'Volume','Close','Adj Close']].to_numpy()
# MinMax
scaler = preprocessing.MinMaxScaler(feature_range=(0.1, 1))
traaain = np.concatenate((goog,aapl),axis=1)
goog = scaler.fit_transform(goog)
aapl = scaler.fit_transform(aapl)
traaain = scaler.fit_transform(traaain)

Теперь нам нужно инициализировать загрузчик данных, чтобы подготовить данные в пакетном режиме. Для этого мы определим функцию ниже.

def data_for_training(data_s1, period1):
a1 = np.arange(len(data_s1) - period1 + 1)
data_list1 = []
for i in range(len(a1)):
sample_list = []
for j in range(period1):
value1 = data_s1[i + j]
sample_list.append(value1)
data_list1.append(sample_list)
data_list1 = np.array(data_list1)
# np.random.shuffle(data_list1)
return data_list1

Эта функция берет данные и период и подготавливает пакеты с периодом размера. Каждая партия будет перекрываться со значениями периода-1 предыдущей партии.
Теперь давайте подготовим наши партии. Мы выберем 30-дневный период. В качестве цели мы выбираем цену закрытия.

period = 30
data_list_goog = data_for_training(goog, period)
data_list_aapl = data_for_training(aapl, period)
data_train = data_for_training(traaain,30)

target_train = np.zeros((len(data_list_aapl), 2))
for i in range(len(data_list_aapl)):
value1 = data_train[i][29][4]
value2 = data_train[i][29][10]
target_train[i][0] = value1
target_train[i][1] = value2

data_train = data_for_training(traaain,30)

Теперь давайте разделим наши данные для обучения и тестирования. Мы будем использовать 10% данных для теста.

my_test_size = 0.1
ratio_train = int((1 - my_test_size) * len(data_list_goog))
X_train = data_train[:ratio_train]
y_train = target_train[:ratio_train]
X_test = data_train[ratio_train:]
y_test = target_train[ratio_train:]

Создание сети LSTM для прогнозирования

Теперь у нас есть готовые обучающие данные. Давайте построим нашу сеть.

from tensorflow.keras import Sequential, Model, Input
from tensorflow.keras.layers import Dense, LSTM, Dropout, Concatenate, concatenate
from tensorflow.keras.callbacks import EarlyStopping, ModelCheckpoint, ReduceLROnPlateau
import tensorflow as tf

data_dim = 12 # we have 6 dimensions for each data. Therefore total dimension is 12 in total as we concatenated the datasets.
timesteps = 30
nb_classes = 2 # we have to predict closing prices for the two companies.
model = Sequential()
model.add(LSTM(64, batch_input_shape=(None, timesteps, data_dim), return_sequences=True, recurrent_dropout=0))
model.add(LSTM(64, return_sequences=False, recurrent_dropout=0))
model.add(Dense(nb_classes, activation='relu'))
model.compile(loss='mse',
optimizer='adam',
metrics=['mae'])

model.summary()
Сетевая архитектура

А теперь пусть начнется обучение!

history = model.fit(X_train,y_train,
batch_size=64, epochs=100)

Построение графиков потерь

results_train = model.evaluate(X_train, y_train)
results_test = model.evaluate(X_test, y_test)

print('train loss: {}, train mae: {}'.format(results_train[0], results_train[1]))
print('test loss: {}, test mae: {}'.format(results_test[0], results_test[1]))

plt.plot(history.history['loss'])
plt.title('Loss plot')
plt.xlabel('Epochs')
plt.ylabel('Loss')
plt.show()

plt.plot(history.history['loss'])
plt.yscale('log')
plt.title('Log Loss plot')
plt.xlabel('Epochs')
plt.ylabel('Loss')
plt.show()
Потери и потери журнала во время тренировки против эпохи

Результаты

Поскольку мы выполнили некоторую предварительную обработку обучающих данных, мы должны выполнить обратное преобразование прогнозируемых данных. Для этого мы определяем функцию ниже

def inv_transform(scaler,data):
ttest = np.zeros((len(data),12))
ttest[:,10]=data[:,1]
ttest[:,4]=data[:,0]
return scaler.inverse_transform(ttest)

Теперь давайте посмотрим, насколько хороша наша модель в прогнозировании цен.

z=model.predict(X_test)
zz = inv_transform(scaler,y_test)
scaled_results = inv_transform(scaler,z)
plt.plot(zz[:,4])
plt.plot(scaled_results[:,4])
plt.title('GOOG')
plt.xlabel('Epochs')
plt.ylabel('Loss')
plt.legend(["Real","Predicted"])
plt.show()

plt.plot(zz[:,10])
plt.plot(scaled_results[:,10])
plt.title('AAPL')
plt.xlabel('Epochs')
plt.ylabel('Loss')
plt.legend(["Real","Predicted"])
plt.show()
Прогнозируемая цена закрытия по сравнению с реальной ценой акций Apple и Google

Заключение

В этой статье мы продемонстрировали, как построить многомерную модель прогнозирования временных рядов на основе LSTM. Мы показали, как можно использовать одну сеть для предсказания двух временных рядов с разными масштабами

Набор данных и полный код доступны на GitHub.

Прогнозирование цен и тенденций Dow Jones с помощью evoML

Прогнозирование является важнейшим компонентом в принятии решений. Для того чтобы принимать правильные решения в настоящем, мы должны уметь, по крайней мере, с некоторой точностью, предвидеть будущие события и последствия наших действий для этих будущих инцидентов.

Прогнозирование временных рядов

Машинное обучение может быть применено к данным временных рядов, чтобы помочь нам прогнозировать будущее и принимать оптимальные решения в настоящем. Временной ряд — это любой набор значений, которые были измерены в последовательное время и размещены в хронологическом порядке. Существует огромное количество примеров таких данных, почерпнутых из инцидентов, которые затрагивают нас ежедневно. Некоторые из этих примеров включают данные о погоде (например, ветер, температура и влажность), солнечную активность, спрос на электроэнергию и финансовые данные (например, цены на акции, обменные курсы и цены на сырьевые товары).

Упорядоченная во времени структура этих наборов данных позволяет нам применять новые инструменты при составлении прогнозов. Например, мы можем изучать корреляции между прошлыми и настоящими значениями ряда, чтобы делать прогнозы на будущее. Тем не менее, мы должны действовать осторожно, так как в этой ситуации очень легко допустить ошибки, такие как обучение модели прогнозирования использованию данных, которые не будут доступны на момент составления прогноза.

В этом блоге мы обсудим эти проблемы со ссылкой на промышленный индекс Доу-Джонса (DJIA) и то, как мы справляемся с этими проблемами с помощью платформы evoML. Читатели могут дополнительно применять фундаментальные концепции и evoML в других проектах временных рядов.

Шаблон данных изменяется с течением времени

Основное предположение контролируемого машинного обучения заключается в том, что шаблоны, полученные из набора обучающих данных, могут быть применены к новому невидимому набору данных. Часто это работает очень хорошо: это предположение лежит в основе науки в целом. Мы делаем наблюдения, находим модель, объясняющую эти наблюдения, и используем эту модель для прогнозирования будущего поведения. Например, используя законы движения и тяготения Ньютона, мы можем довольно точно предсказать движение планет. Хотя теперь у нас есть лучшие модели, эта ньютоновская модель сегодня так же точна, как и в 17 веке.

Но мы не всегда можем воспринимать эту стабильность как нечто само собой разумеющееся. Многие модели машинного обучения относятся к моделям поведения человека, которые могут меняться со временем, и это означает, что модели, обученные на старых данных, не обязательно будут действительны в будущем.

Прогнозирование промышленного индекса Доу-Джонса

Индекс Доу-Джонса является одним из старейших фондовых индексов в мире и имеет историю, восходящую к 1885 году. Он обеспечивает средневзвешенное значение цен на акции 30 крупных публично торгуемых американских компаний. Мы можем получить ежедневные цены закрытия Dow Jones из TradingView.

Рисунок 1. Дневные цены закрытия индекса Доу-Джонса

Теперь мы задаемся простым вопросом: будет ли завтрашняя цена закрытия выше сегодняшней?

Это сводит проблему прогнозирования к двоичной классификации. Чтобы обучить модель, мы должны сгенерировать функции и цель из наших данных. Для прогнозирования тренда цен закрытия мы будем использовать предыдущие пять дней цен закрытия (не включая выходные и праздничные дни, когда торговля остановлена).

Прежде чем мы обучим классификатор предсказывать тренд, мы должны быть осторожны с тем, как измерить точность нашей модели.

Рисунок 2. Предыдущие пять дней цен закрытия

Почему перекрестная проверка K-fold не работает

Обычно одним из самых популярных методов для этой задачи является k-кратная перекрестная проверка. В этой процедуре обучающие данные сначала разбиваются на k отдельных блоков. Затем у нас есть k раундов обучения и проверки, в которых каждый из блоков откладывается, модель обучается на оставшихся данных, а ее производительность оценивается на блоке, который мы исключили. Показатели оценки, такие как точность классификации, могут быть усреднены за k раундов, чтобы дать общую меру производительности модели.

Рисунок 3. К-кратная перекрестная проверка. Изображение предоставлено TurinTech.

Однако, когда мы рассматриваем данные временных рядов, эта процедура должна быть изменена. Каждый раз, когда мы оцениваем модель, мы должны быть уверены, что используемая нами процедура отражает то, как модель будет использоваться на практике. Очевидно, что если мы обучаем модель, чтобы она говорила нам о будущем, мы сможем обучить ее только на данных, которые мы собрали до того, как будет сделан прогноз. Тем не менее, если мы используем перекрестную проверку k-fold, то некоторые из наших данных проверки будут предшествовать соответствующим обучающим данным.

Проверка модели таким образом может дать неточную картину производительности модели, позволяя модели изучать будущую информацию для прогнозирования прошлого поведения временных рядов. Такого рода смещение будет проблемой, особенно если факторы, управляющие динамикой временных рядов (например, такие факторы, как процентные ставки, обменные курсы и т. д.), изменяются с течением времени.

Хронологически упорядоченное обучение и валидация

Альтернативой, которую мы применяем, является использование хронологически упорядоченных блоков обучения и проверки. Как и в k-кратной перекрестной проверке, мы делим данные на k блоков, но теперь мы проходим k-1 раунды обучения и проверки. Мы начинаем с обучения на первом блоке и валидации на втором. В каждом последующем раунде мы расширяем обучающий набор, добавляя предыдущий блок проверки к обучающему набору и проверяя следующий блок, пока не охватим весь набор данных. Это позволяет нам проверять нашу модель во многих различных окнах во времени, всегда обучая ее на прошлых данных.

Рисунок 4. Разделение временных рядов. Изображение предоставлено TurinTech.

Распараллеливание моделей машинного обучения, чтобы найти лучшие из них

Используя платформу evoML, мы можем быстро создать пробную версию для сравнения и оптимизации сотен моделей одновременно.

Мы можем просто перетащить данные Dow Jones через интерфейс без кода на платформе evoML.

Рисунок 5. Импортируйте данные из любого места. Изображение предоставлено TurinTech.

Сначала мы проводим перекрестную проверку, используя данные между 1885 и 1950 годами. Затем мы тестируем данные с 1950 по 1980 год. После сравнения классификаторов LightGBM, случайного леса и логистической регрессии мы обнаружили, что логистическая регрессия работает лучше всего, демонстрируя точность 0,539 во время проверки модели.

Цифра. 6 Генерация и оценка моделей на evoML. Изображение предоставлено TurinTech.

Развертывание модели и прогноза

Мы разрабатываем простую торговую стратегию. Мы начинаем со 100 долларов. Если мы прогнозируем, что индекс Доу-Джонса вырастет на следующем временном шаге, то мы будем либо покупать, либо держать (в зависимости от того, были ли уже вложены наши деньги или нет). Но если мы прогнозируем падение, то мы продаем наши акции (или ничего не делаем, если нас не инвестируют). Мы проверяем, как эта стратегия работала бы в период с 1950 по 1980 год.

Рисунок 7. Эффективность торговой стратегии в период с 1950 по 1980 год. Изображение предоставлено TurinTech.

Если бы мы просто купили акции компаний, входящих в индекс Доу-Джонса, наши инвестиции выросли бы в цене до 419 долларов к 1980 году. Но, используя нашу стратегию, мы могли бы увеличить стоимость наших инвестиций до 5 708 долларов, опередив рынок в 13,6 раза.

Постоянная оптимизация модели для достижения максимальной производительности

Эта базовая стратегия сработала бы на удивление хорошо. Но как это работает в наши дни? К сожалению, не все так хорошо. Если мы повторим то же упражнение, но охватим период с 1980 года до начала 2022 года, мы обнаружим, что наши первоначальные 100 долларов увеличились бы до 864 долларов, если бы мы следовали нашей торговой стратегии, но увеличились бы до 4,333 долларов, если бы мы просто следовали за рынком.

Рисунок 8. Эффективность торговой стратегии в период с 1980 по 2020 год. Изображение предоставлено TurinTech.

Таким образом, наша предыдущая стратегия сейчас работает очень плохо. Почему так? Общеизвестно, что прогнозирование фондового рынка – очень сложная задача. Внимательно изучив его, мы сможем узнать определенные корреляции между прошлым и настоящим поведением, но нет никакой гарантии, что эти корреляции будут сохраняться в будущем.

Законы, регулирующие фондовый рынок, являются продуктом человеческого поведения и могут меняться со временем. Например, фондовый рынок находится под сильным влиянием наших усилий предсказать его таким образом, который явно не применим к другим системам, таким как движение планет или погода. Если мы считаем, что акции вырастут в цене, мы можем их купить, но другие трейдеры также могут быть в процессе покупки этих акций, потому что они пришли к такому же выводу. Эффект всех этих ордеров на покупку заключается в том, чтобы увеличить цену акции и помешать нам купить ее по низкой цене, которую мы предполагали. Каждый трейдер конкурирует с любым другим трейдером, а достижения в области информационных технологий и математики делают старые торговые стратегии устаревшими. Поэтому невозможно просто выбрать стратегию и придерживаться ее. Модели должны постоянно контролироваться и обновляться, чтобы обеспечить будущий успех.

Вот где evoML может помочь. Автоматически оптимизируя десятки потенциальных моделей и сотни их параметров, мы даем пользователям возможность быстро адаптировать и тестировать новые подходы на основе их пользовательских критериев (например, точности, объяснимости и времени прогнозирования). Затем наилучшие результаты могут быть легко развернуты вместе с их полным исходным кодом, чтобы пользователь не только получил лучшую модель для работы, но и понял, что это за модель и как она работает.

Рисунок 9. Объяснение модели с полным кодом на evoML. Изображение предоставлено TurinTech.

Прогнозирование временных рядов с помощью Prophet и API

Один из самых простых, но эффективных методов прогнозирования временных рядов в Python

Прогнозирование временных рядов является мощным инструментом в области предиктивного финансового анализа. Независимо от того, идет ли речь о ценах на акции, географической сегментации или выручке в зависимости от продукта, способность предсказывать будущие тенденции имеет решающее значение для принятия обоснованных решений. В этой статье мы рассмотрим, как мы можем извлечь ценную информацию из API подготовки к финансовому моделированию (FMP), предоставляя дорожную карту как для отдельных лиц, так и для компаний для улучшения их стратегического планирования.

Импорт необходимых пакетов

import requests
import pandas as pd
import matplotlib.pyplot as plt
from statsmodels.tsa.stattools import adfuller
from prophet import Prophet

plt.rcParams[‘figure.figsize’] = (20,10)
plt.style.use(‘fivethirtyeight’)

requests: предоставляет простой API для взаимодействия с HTTP-операциями, такими как GET POST т. д.

pandas — это быстрый, мощный, гибкий и простой в использовании инструмент для анализа и обработки данных с открытым исходным кодом.

matpllotlib для создания статических, анимированных и интерактивных визуализаций на Python.

adfuller для проведения расширенного теста Дики-Фуллера (ADF) для проверки стационарности в данных временных рядов.

Prophet для прогнозирования временных рядов, разработанная Facebook. Класс Prophet используется для создания моделей прогнозирования.

Прежде чем двигаться дальше, для извлечения данных мы будем использовать API FinancialModelingPrep. Таким образом, для бесперебойного выполнения предстоящего блока кода вам необходимо иметь учетную запись разработчика FMP, которую вы можете легко создать по ссылке здесь.

Извлечение информации из API FMP

Следующая функция extract_geographical_revenue помогает в извлечении важной информации о доходах компании в различных географических регионах. Это важно для понимания глобального распределения доходов компании.def extract_geographical_revenue(symbol):
api_key = ‘YOUR API KEY’
url = f’https://financialmodelingprep.com/api/v4/revenue-geographic-segmentation?symbol={symbol}&period=quarter&structure=flat&apikey={api_key}’
raw_df = requests.get(url).json()

# Convert JSON to DataFrame
data = []
for entry in raw_df:
for date, values in entry.items():
row = {‘Date’: date}
for country, value in values.items():
row[country] = value
data.append(row)

df = pd.DataFrame(data)

return df

df_geographical = extract_geographical_revenue(‘TSLA’)
df_geographical.head()

Изображение автора

Следующая функция extract_product_revenue извлечения информации о том, как выручка компании распределяется по различным категориям продуктов. Эта сегментация имеет решающее значение для понимания того, какие продукты вносят наибольший вклад в общий доход компании.def extract_product_revenue(symbol):
api_key = ‘YOUR API KEY’
url = f’https://financialmodelingprep.com/api/v4/revenue-product-segmentation?symbol={symbol}&period=quarter&structure=flat&apikey={api_key}’
raw_df = requests.get(url).json()
# Create an empty DataFrame
df = pd.DataFrame()

# Iterate through the data and append to the DataFrame
for entry in raw_df:
for date, values in entry.items():
df = pd.concat([df, pd.DataFrame.from_dict({date: values}, orient=’index’)])

# Reset the index
df.reset_index(inplace=True)
df.rename(columns={‘index’: ‘Date’}, inplace=True)

return df

df_product = extract_product_revenue(‘TSLA’)
df_product.head()

Изображение автора

Изучение и подготовка набора данных

Перед обучением модели необходимо изменить данные в нужном формате, чтобы модель лучше понимала и выполняла анализ. В этом разделе мы рассмотрим то же самое.# Changing all the NaN values to zero for further modifications.
df_geographical = df_geographical.fillna(0)

# Merging all the United States columns.
df_geographical[‘United States’] = df_geographical[‘UNITED STATES’] + df_geographical[‘U S’] + df_geographical[‘U’]

# Merging all the China colums.
df_geographical[‘China’] = df_geographical[‘CHINA’] + df_geographical[‘C N’]

# Drop the unused columns.
df_geographical = df_geographical.drop([«UNITED STATES» ,»U S», ‘U’,’CHINA’, ‘C N’],axis =1)

# Reversing the rows to make it from earlier to newer.
df_geographical = df_geographical .iloc[::-1]

# Reseting the index
df_geographical = df_geographical.reset_index().drop([‘index’],axis=1)
df_geographical.head()

# Set the ‘DATE’ column as the index
df_geographical = df_geographical.set_index(‘Date’)
df_geographical.index = pd.to_datetime(df_geographical.index)

# Dropping the NA values.
df_geographical.dropna(inplace = True)

df_geographical.head()

Изображение автора

Общие сведения о стационарности

Стационарность является фундаментальным понятием в анализе временных рядов. Стационарный временной ряд — это ряд, статистические свойства, такие как среднее значение и дисперсия, остаются постоянными во времени. Это означает, что ряд имеет последовательное поведение, а его закономерности предсказуемы в разные периоды времени.

Стационарные временные ряды обладают следующими свойствами:

Среднее постоянное: Среднее значение ряда остается неизменным во времени.

Постоянная дисперсия: Дисперсия (или стандартное отклонение) ряда остается постоянной с течением времени.

Постоянная автоковариация: Связь между наблюдениями в разные временные лаги остается постоянной.

Стационарность важна, потому что многие методы анализа временных рядов и модели прогнозирования предполагают стационарность или лучше всего работают со стационарными данными.# Visualize the time series to find stationarity and patterns

plt.plot(df_geographical.index, df_geographical[‘United States’])
plt.xlabel(‘Quarters’)
plt.ylabel(‘United States production’)
plt.title(‘Production from United States each quarter’)
plt.show()

Изображение автора

Производство Tesla в Соединенных Штатах растет с увеличением наклона от квартала к кварталу.

Тест ADF помогает определить, является ли временной ряд стационарным или нет. Он предоставляет статистику ADF, p-значение и критические значения для различных уровней значимости,

Если p-значение меньше выбранного уровня значимости (например, 0,05), то можно отвергнуть нулевую гипотезу и сделать вывод о стационарности временного ряда. В противном случае, если p-значение больше уровня значимости, мы не можем отвергнуть нулевую гипотезу, предполагая, что временной ряд нестационарен.# Perform the Augmented Dickey-Fuller test

result = adfuller(df_geographical[‘United States’])
print(‘ADF Statistic:’, result[0])
print(‘p-value:’, result[1])
print(‘Critical Values:’)
for key, value in result[4].items():
print(key, ‘:’, value)

Изображение автора

Анализ трендов

Тренд помогает выявить основное направление или закономерность во временном ряду, что может быть полезно для прогнозирования и принятия решений.

Для выполнения анализа тренда на временном ряду мы можем использовать различные методы, включая визуальный осмотр, скользящие средние и регрессионный анализ.

Визуальный осмотр: Построение графика данных временных рядов часто может выявить наличие тренда. Четкое восходящее или нисходящее движение с течением времени говорит о наличии трендовой составляющей. Визуальный осмотр позволяет наблюдать за общей картиной и выявлять любые отклонения или изменения в ряду.

Скользящие средние: Скользящие средние широко используются для анализа тренда. Они помогают сглаживать краткосрочные колебания данных, облегчая выявление основного тренда. К распространенным типам скользящих средних относятся простая скользящая средняя (SMA), взвешенная скользящая средняя (WMA) и экспоненциальная скользящая средняя (EMA).

Вот три распространенных типа скользящих средних, используемых для анализа тренда:

Простая скользящая средняя (SMA): Простая скользящая средняя вычисляет среднее значение заданного количества точек данных за определенное окно. Каждая точка данных вносит равный вклад в среднее значение, и более старые наблюдения имеют такой же вес, как и новые. SMA обеспечивает сглаженное представление данных за счет уменьшения случайных флуктуаций.

Взвешенная скользящая средняя (WMA): Взвешенная скользящая средняя присваивает различные веса точкам данных в пределах окна. Весовые коэффициенты могут быть линейными или следовать определенной схеме. WMA уделяет больше внимания недавним наблюдениям, что позволяет ей быстрее реагировать на изменения данных по сравнению с SMA.

Экспоненциальная скользящая средняя (EMA): Экспоненциальная скользящая средняя — это тип взвешенной скользящей средней, которая присваивает экспоненциально убывающие веса точкам данных. Он придает большее значение недавним наблюдениям, постепенно снижая важность более старых наблюдений. EMA более чутко реагирует на последние изменения данных и часто используется в техническом анализе.

Регрессионный анализ: Регрессионный анализ может быть применен для оценки и визуализации трендовой составляющей временного ряда. Он включает в себя подгонку регрессионной модели к данным, где время рассматривается как независимая переменная, а интересующая нас переменная является зависимой переменной. Затем компонент тренда может быть извлечен из регрессионной модели.# Calculate the rolling mean (simple moving average) with a window size of 3 months
rolling_mean = df_geographical[‘United States’].rolling(window=3).mean()

# Visualize the original time series and the trend component
plt.plot(df_geographical.index, df_geographical[‘United States’], label=’Original’)
plt.plot(df_geographical.index, rolling_mean, color=’red’, label=’Trend (Moving Average)’)
plt.xlabel(‘Quarter’)
plt.ylabel(‘Production from United States’)
plt.title(‘Trend Analysis: United States Production’)
plt.legend()
plt.show()

Изображение автора

Обучение модели и прогнозирование

После того, как данные приведены в надлежащий формат, мы готовы к обучению модели и прогнозированию будущего. Вот так:# Format data for prophet model using ds and y
df_geographical = df_geographical.reset_index() \
.rename(columns={‘Date’:’ds’,
‘United States’:’y’})

# Setup and train model and fit
model = Prophet()
model.fit(df_geographical)

# define the period for which we want a prediction
future = [‘2023-09-30′,’2023-12-30′,’2024-03-30′,’2024-06-30′,’2024-09-30′,’2024-12-30′,’2025-03-30′,’2025-06-30′,’2025-09-30’]
future = pd.DataFrame(future)
future.columns = [‘ds’]
future[‘ds’]= pd.to_datetime(future[‘ds’])

forecast = model.predict(future)

# plot forecast
model.plot(forecast)
plt.show()

Мы видим, что модель предсказывает, что производство снизится по сравнению с предыдущим годом, но увеличится в будущем.

Заключение

В заключение следует отметить, что прогнозирование временных рядов меняет правила игры в финансовом мире. Использование API-интерфейсов, таких как Financial Modeling Prep, открывает доступ к сокровищнице данных, позволяя отдельным лицам и компаниям принимать стратегические решения на основе ценной информации. От очистки и подготовки данных до анализа тенденций и прогнозирования — эта статья служит исчерпывающим руководством для всех, кто хочет ориентироваться в сложном ландшафте анализа финансовых временных рядов.

Прогнозирование временных рядов с использованием LSTM

Знакомство

Вот среднемесячные температуры во Флориде с 1743 года, о которых я говорил в предыдущей статье. В этой статье я сначала продемонстрирую операции базового уровня с использованием LSTM с этим набором данных.

Не объясняйте мне, покажите код

Доступ к репозиторию GitHub можно получить здесь.

Что такое LSTM и почему он важен для временных рядов?

Long short-term memory (LSTM) — это архитектура искусственной повторяющейся нейронной сети (RNN), используемая в области глубокого обучения. Хотя он не отличается от RNN с точки зрения рабочей логики, он позволяет работать гораздо более длинным последовательностям. По мнению RNN, это может в значительной степени решить проблему градиента.

Короче говоря, модели LSTM могут хранить информацию в течение определенного периода времени. Благодаря этой функции LSTM использование LSTM чрезвычайно полезно при работе с временными рядами или последовательными данными. Конечно, поскольку эта статья является частью серии статей о временных рядах, подробное объяснение LSTM выходит за рамки этой статьи. Но его использование может дать вам некоторое представление.

Предпочтительные и возможные функции активации

Когда мы знаем, что нужно предсказать, функция активации, которую мы предпочитаем, может появиться спонтанно. Если мы хотим продолжить без нашего ведома на этом этапе, он вернется к нам таким образом, что увеличит время обучения.
Например, сигмоида хорошо работает для классификатора, который мы можем понять, просто взглянув на график для сигмоиды.

График сигмовидной функции из https://en.wikipedia.org/wiki/Sigmoid_function

С другой стороны, если мы хотим пойти легко, ReLU вычислительно проще, чем tanh и sigmoid, поскольку он включает в себя более простые математические операции и требует меньше энергии.

График функций ReLU из https://en.wikipedia.org/wiki/Rectifier_(neural_networks)

Итак, приступим

0. Загрузка и подготовка данных

import numpy as np
import matplotlib.pyplot as plt
import pandas as pd
import datetime as dt

from sklearn.preprocessing import MinMaxScaler
from tensorflow.keras.models import Sequential
from tensorflow.keras.layers import Dense, Dropout, LSTM

# Load Data
florida = pd.read_csv(‘github/florida_file_.csv’)
florida.head()
florida.tail(5)
florida[«Date»] = pd.to_datetime(florida[«Date»])

florida = florida[[«Date», «Avg_Temp»]]
# florida = florida[«Avg_Temp»].resample(‘MS’).mean()
florida = florida.fillna(florida.bfill())
florida.columns = [‘Date’, ‘Avg_Temp’]

train = florida[:-225]
len(train)
test = florida[-225:]
len(test)
train_dates = pd.to_datetime(train[‘Date’])
test_dates = pd.to_datetime(test[‘Date’])

После установки необходимых библиотек и загрузки данных я выполняю разделение. Здесь набор тестовых данных будет представлять собой данные, которые я зарезервировал для использования в конце, чтобы использовать значения, которые я предскажу с помощью модели. В будущем я разделю обучающий набор данных на x_train и y_train, но перед всем этим я масштабирую данные. (Масштабирование данных очень важно при работе с нейронными сетями, эта модель значительно сокращает время обучения.)

1. Подготовьте данные

scaler = MinMaxScaler(feature_range=(0,1))
scaled_data = scaler.fit_transform(train[‘Avg_Temp’].values.reshape(-1,1))

prediction_days = 225

x_train = []
y_train = []

for x in range(prediction_days, len(scaled_data)):
x_train.append(scaled_data[x-prediction_days:x, 0])
y_train.append(scaled_data[x, 0])

x_train, y_train = np.array(x_train), np.array(y_train)
x_train = np.reshape(x_train, (x_train.shape[0], x_train.shape[1], 1))

После масштабирования всех значений в наборе данных от 0 до 1 я указываю количество дней для прогнозирования. После этого значения я разделяю данные как x_train и y_train. При приведенном нами значении 225 модель будет работать таким образом, что она исследует 225 данных и предсказывает следующую, затем снова исследует 225 данных и пытается предсказать следующую.

Заполнив наши списки, мы преобразуем их в массивы NumPy. Затем мы перекраиваем x_train так, чтобы x_train могли работать с нейронной сетью.

3. Постройте модель

model = Sequential()

model.add(LSTM(units =128, activation=’relu’, return_sequences=True, input_shape = (x_train.shape[1],1)))
model.add(Dropout(0.2))
model.add(LSTM(units =128, activation=’relu’, return_sequences=True))
model.add(Dropout(0.2))
model.add(LSTM(units =128, activation=’relu’, return_sequences=False))
model.add(Dropout(0.2))
model.add(Dense(units=1)) # Prediction of the next value

На этом шаге мы настроили нашу модель, поскольку это базовая нейронная сеть, мы настроили последовательную модель напрямую. Здесь, выбрав слои LSTM, я применил Dropout для повышения эффективности и устранения бесполезных узлов, чтобы после каждого слоя LSTM я не переходил с худшими узлами к следующему шагу. Здесь я устранил 20% от общего количества узлов.

Внимание!! Мы должны присвоить input_shape значение нашему первому слою LSTM, потому что наша модель не знает размер данных, которые она будет обрабатывать, нам не нужно вводить это значение на следующих шагах. На последнем шаге я вернул одно значение с плотным слоем, поэтому это оценочное значение будет средней температурой, предсказанной нашей моделью.model.compile(optimizer=’adam’, loss=’mean_squared_error’)
model.summary()

history = model.fit(x_train, y_train, epochs = 25, batch_size=32, validation_split=0.1)

Затем мы компилируем нашу модель, готовим набор тестовых данных, который наша модель увидит вместе с готовой моделью. Наконец, мы перекраиваем x_test и приводим его к формату, соответствующему нашей модели.

4. Прогноз

actual_temp = test[‘Avg_Temp’].values
total_temp = pd.concat((train[‘Avg_Temp’], test[‘Avg_Temp’]),axis=0)

model_inputs = total_temp[len(total_temp)-len(test)-prediction_days:].values
model_inputs = model_inputs.reshape(-1,1)
model_inputs = scaler.transform(model_inputs)

# Make Predictions on Test Data
x_test = []

for x in range(prediction_days, len(model_inputs)):
x_test.append(model_inputs[x-prediction_days:x, 0])

x_test = np.array(x_test)
x_test = np.reshape(x_test, (x_test.shape[0], x_test.shape[1], 1))

pred = model.predict(x_test)
pred = scaler.inverse_transform(pred)

Затем мы выполняем операцию прогнозирования на нашей модели. Затем мы нормализуем эти оценочные значения с помощью операции inverse_transform, чтобы сравнить эти прогнозируемые значения с фактическими значениями, которые у нас есть.

Учитывая, что ценность нашей эпохи является низкой величиной, не будет ошибкой сказать, что успех нашей модели находится на приемлемом уровне.

Ссылки

[1] https://www.veribilimiokulu.com

[2] https://www.analyticsvidhya.com/blog/2020/10/how-to-create-an-arima-model-for-time-series-forecasting-in-python/

[3] https://www.statisticssolutions.com/time-series-analysis/

[4] https://www.tensorflow.org/api_docs/python/tf/keras/layers/LSTM

Transformers vs. LSTM для прогнозирования временных рядов цен на акции

Диаграмма трансформатора

Поэтому я подумал, что должен продолжить обсуждение проблемы биржевых временных рядов, которую я начал в своем первом посте в блоге. В этом первом посте мы использовали архитектуры моделей LSTM и CNN-LSTM для прогнозирования будущих цен на акции, используя только предыдущие цены в качестве входных данных, и у нас это получилось довольно хорошо. Наша средняя абсолютная процентная ошибка была на уровне или ниже трех процентов, и все выглядело хорошо. В этом посте я упоминал, что буду обсуждать способы улучшения этой цифры, и первоначально я намеревался сделать это, создав ансамблевую LSTM-модель с использованием LSTM-сетей с различными гиперпараметрами, а затем усреднив их предсказание, чтобы получить (надеюсь) что-то лучшее. Но должен признаться, что я немного отвлекся, занявшись другим, на мой взгляд, более интересным направлением исследования. Одним словом, этот новый подход – «Трансформеры».

Что такое архитектура модели трансформатора?

Прежде чем мы перейдем к конкретным особенностям проблемы, я подумал, что было бы неплохо немного экспозиции. Однако следует отметить, что трансформеры — это высокоуровневая концепция машинного обучения, в которую я только погружаюсь, и поэтому я далек от эксперта в предметной области. Но при этом преподавание (или попытка) — неотъемлемая часть обучения, так что в духе прогресса вот оно.

На очень высоком уровне трансформатор представляет собой архитектуру сети с прямой связью, которая использует «механизм самовнимания» для сокращения времени обучения и количества параметров, оставаясь при этом высокопрогностическим. Я сформулировал это следующим образом: трансформаторы могут выполнять работу РНС без всех этих повторений, замедляющих работу и добавляющих параметры. Это работает следующим образом: механизм многоголового внимания преобразователя (реализованный tf.keras.layers.multiheadattentionв tensorflow) позволяет модели отслеживать каждую точку данных, относящуюся к конкретной другой точке данных. Если мы пытаемся предсказать завтрашнюю цену акций, то сегодняшняя цена явно имеет значение, но и другие цены тоже. Например, ценовое действие на этой неделе может быть очень похоже на то, что произошло три месяца назад, и поэтому для нас будет полезно, если наша модель сможет «запомнить» это ценовое действие. Эта потребность в «памяти» была одной из причин, по которой мы решили использовать рекуррентные нейронные сети и модели LSTM (длительная краткосрочная память). Позволяя предыдущим точкам данных передаваться в нашу модель в обратном направлении, RNN могут взвешивать предыдущие значения, создавая эффект памяти. Что отличает трансформеров, так это их уникальный механизм самовнимания, означающий, что они могут иметь тот же самый эффект «памяти» без всех повторений: данные подаются один раз, а самовнимание отслеживает соответствующие точки данных на каждом шаге (оно взвешивает то ценовое действие, которое было похоже на то, что наблюдалось в настоящее время, скажем, три месяца назад). Это означает, что трансформаторы обучаются намного быстрее, чем RNN, поскольку мы вводим данные только один раз, и им требуется гораздо меньше параметров, что является беспроигрышным вариантом. Для более убедительного и экспертного объяснения трансформаторов я дам ссылку на оригинальную статью, описывающую их, в конце этого поста. А пока перейдем к нашей проблеме.

Задача биржевых временных рядов

Вопрос, на который пытается ответить это исследование, можно сформулировать так: «Можем ли мы предсказать цены на акции в течение следующих пяти дней, используя только предыдущие цены акций?». Как я уже упоминал выше, мы уже частично ответили на этот вопрос предварительным «да». Тем не менее, модели, которые мы построили, хотя и удивительно искусны, далеки от совершенства или даже частичной полезности: в среднем 3,5% недостаточно, чтобы чувствовать себя уверенно в торговле. Отчасти это связано с тем, что проблема просто сложная, данные зашумлены и почти случайны. Оказывается, фондовый рынок не зря называют случайной прогулкой. Но даже в этом случае нас это не останавливает, и поэтому мы собираемся двигаться вперед и отмечать следующие шаги в поисках лучшей модели:

  1. Создайте надежную базовую модель LSTM для оценки производительности наших трансформаторов.
  2. Построение модели трансформатора
  3. Обучайте и оценивайте наши модели

Часть первая: LSTM

Давайте построим нашу базовую модель. Из предыдущих экспериментов я обнаружил, что наличие трех слоев, слоя LSTM с 200 нуэронами и двух плотных узлов с 50 и 1 нуэроном (нуэронами) соответственно, очень хорошо работает для решения этой задачи. Вот код для реализации этого:from keras.models import Sequential
from keras.layers import Dense
from keras.layers import LSTM

def build_lstm(etl: ETL, epochs=25, batch_size=32) -> tuple[tf.keras.Model, tf.keras.History]]:
«»»
Builds, compiles, and fits our LSTM baseline model.
«»»
n_timesteps, n_features, n_outputs = 5, 1, 5
callbacks = [tf.keras.callbacks.EarlyStopping(patience=10,
restore_best_weights=True)]
model = Sequential()
model.add(LSTM(200, activation=’relu’,
input_shape=(n_timesteps, n_features)))
model.add(Dense(50, activation=’relu’))
model.add(Dense(n_outputs))
print(‘compiling baseline model…’)
model.compile(optimizer=’adam’, loss=’mse’, metrics=[‘mae’, ‘mape’])
print(‘fitting model…’)
history = model.fit(etl.X_train, etl.y_train,
batch_size=batch_size,
epochs=epochs,
validation_data=(etl.X_test, etl.y_test),
verbose=1,
callbacks=callbacks)
return model, history

А вот краткое описание модели, связанное с этой моделью после ее компиляции:

Итак, это базовый уровень, который мы собираемся использовать для оценки нашего трансформатора, мы обсудим его результаты в разделе оценки, а пока давайте поговорим об архитектуре нашего трансформатора.

Часть вторая: Трансформер

Теперь для нашей архитектуры трансформатора мы будем использовать конструкцию, рекомендованную в документации keras (я дам ссылку в конце поста), но эта конфигурация создана для классификации, поэтому мы внесем небольшое изменение; Мы изменим функцию активации конечных выходных слоев с Softmax на Relu, а нашу функцию потерь на Mean squared Error. Кроме того, я установил гиперпараметры на то, что, по моему мнению, лучше всего подходит для этой задачи с помощью экспериментов. Вот что мы получаем:def transformer_encoder(inputs, head_size, num_heads, ff_dim,
dropout=0, attention_axes=None):
«»»
Creates a single transformer block.
«»»
x = layers.LayerNormalization(epsilon=1e-6)(inputs)
x = layers.MultiHeadAttention(
key_dim=head_size, num_heads=num_heads, dropout=dropout,
attention_axes=attention_axes
)(x, x)
x = layers.Dropout(dropout)(x)
res = x + inputs

# Feed Forward Part
x = layers.LayerNormalization(epsilon=1e-6)(res)
x = layers.Conv1D(filters=ff_dim, kernel_size=1, activation=»relu»)(x)
x = layers.Dropout(dropout)(x)
x = layers.Conv1D(filters=inputs.shape[-1], kernel_size=1)(x)
return x + res

def build_transfromer(head_size,
num_heads,
ff_dim,
num_trans_blocks,
mlp_units, dropout=0, mlp_dropout=0) -> tf.keras.Model:
«»»
Creates final model by building many transformer blocks.
«»»
n_timesteps, n_features, n_outputs = 5, 1, 5
inputs = tf.keras.Input(shape=(n_timesteps, n_features))
x = inputs
for _ in range(num_trans_blocks):
x = transformer_encoder(x, head_size, num_heads, ff_dim, dropout)

x = layers.GlobalAveragePooling1D(data_format=»channels_first»)(x)
for dim in mlp_units:
x = layers.Dense(dim, activation=»relu»)(x)
x = layers.Dropout(mlp_dropout)(x)

outputs = layers.Dense(n_outputs, activation=’relu’)(x)
return tf.keras.Model(inputs, outputs)

transformer = build_transfromer(head_size=128, num_heads=4, ff_dim=2,
num_trans_blocks=4, mlp_units=[256],
mlp_dropout=0.10, dropout=0.10,
attention_axes=1)

transformer.compile(
loss=»mse»,
optimizer=tf.keras.optimizers.Adam(learning_rate=1e-3),
metrics=[«mae», ‘mape’],
)

callbacks = [tf.keras.callbacks.EarlyStopping(patience=10,
restore_best_weights=True)]

t_hist = transformer.fit(data.X_train, data.y_train, batch_size=32,
epochs=25, validation_data=(data.X_test, data.y_test),
verbose=1, callbacks=callbacks)

И этот код построит модель, которая будет выглядеть примерно так (повторяется num_trans_blocks раза).

Это один трансформаторный блок, в нашей модели будет 4 таких блока, сложенных друг на друга.

Наконец, в общей сложности наш трансформаторный блок имеет 17 205 общих параметров, что составляет чуть более одной десятой от нашего LSTM. Теперь приступим к тренировкам.

Часть третья: Обучение и оценка

Обучим обе сети за 25 эпох с помощью оптимизатора Адама. Во-первых, давайте обсудим производительность наших базовых моделей на тестовых данных. Прежде всего, я хотел бы отметить, что LSTM удивительно стабилен, то есть, если вы обучите его 10 раз подряд, он даст вам прогнозы, которые будут работать почти одинаково на тестовом наборе во всех десяти случаях. Кроме того, он тренируется быстрее, чем я ожидал, в общей сложности 143 секунды. Он также хорошо справляется со временем вывода, занимая всего двадцать две секунды. Это тем более впечатляюще, что, по крайней мере, в этом сравнении, LSTM имеет вес в 171 000+ параметров, как упоминалось ранее.

Визуализация наших прогнозов LSTM.

С точки зрения его предсказательной способности, мы считаем, что LSTM является искусным (не удивительно, это наша базовая линия не просто так). Он набирает MAPE 2,44% в нашем тестовом наборе. В целом, LSTM является последовательной, легко обучаемой моделью, которая умело предсказывает данные временных рядов акций. Единственное его ограничение заключается в том, что он большой, и поэтому его нельзя легко масштабировать. Однако, когда дело доходит до цен на акции, на самом деле недостаточно данных для построения очень глубоких моделей, если вы это сделаете, вы фактически начнете терять производительность. Таким образом, я не думаю, что большое количество параметров действительно сильно вредит LSTM.

Вот график Трансформера. Он взят из модели с MAPE 2,8%. Мне не удалось получить визуализацию золотых моделей.

Далее поговорим о трансформаторе. Первое, что хотелось бы отметить, это то, что трансформатор более нестабилен, чем LSTM. Что я имею в виду? Я обучал одни и те же модели (читай: одни и те же гиперпараметры) много-много раз. Каков результат? Я обучил лучшую модель, которую я построил для прогнозирования цен на акции, и пришел с MAPE 2,37%. Я построил еще один с MAPE 2,41%. Оба показателя, как вы заметите, лучше, чем наш базовый уровень, который стабильно весил около 2,45%-2,5%. Но, к сожалению, трансформатор не смог последовательно воспроизвести эту производительность, даже с теми же гиперпараметрами, что и эти два золотых тренировочных прогона. Более того, нельзя сказать, что трансформатор работал немного хуже. Были времена, когда его MAPE превышал 3%. Это проблема, если мы пытаемся строить, а затем переобучать модели каждый месяц или, скажем, квартал. У вас может быть мощная модель трансформатора, вы можете переобучить ее и остаться с кучей хлама. Так что в этом отношении трансформатор не был идеальным.

Так где же на самом деле засиял трансформатор? Подсчет параметров. Он весит чуть более одной десятой по количеству параметров, чем LSTM. Это означало, что трансформатор обучался быстрее, чем LSTM, ему требовалось всего 138 секунд против 143 у LSTM. Но, как ни странно, логический вывод на LSTM был быстрее: трансформатор занял в общей сложности 25 секунд.

Наконец, с точки зрения относительной дисперсии прогнозов по сравнению с фактическими значениями, другими словами, насколько много неопределенности заложено в наши прогнозы, LSTM подстраховался: она составила от 2,4% до 2,6%.

Заключение

Прежде всего, следует отметить, что производительность трансформера была действительно впечатляющей. Имея всего ~ 17 000, он смог не отставать от LSTM с более чем 170 000. Это не подлые ноги. Однако она была непредсказуемой и нестабильной в условиях переподготовки. Я не знаю, является ли это общей особенностью трансформаторов или проблемой с моей реализацией и использованием их, но это определенно было очень заметно. Это может не быть проблемой, если вы обучаете большую языковую модель; Английский не меняется так быстро, чтобы каждые пару недель приходилось выпускать новую модель. Но для прогнозирования цен на акции, возможность переобучаться по мере поступления новых данных, не беспокоясь о том, что ваша ранее очень умелая модель теперь просто посредственна, что ж, это становится очень важным. Я думаю, что по этой причине, а также из-за того, что наша экономия параметров не ускорила время обучения, но настолько, я думаю, что за свои деньги я бы все равно выбрал LSTM. Кроме того, нам не нужно масштабировать эти модели до миллиардов параметров, в том режиме, в котором трансформаторы действительно блистают. Параметры 1B по сравнению с 10B — это совсем другой выбор, чем 17 000 против 170 000.

Так что, хотя мне очень понравилось работать с трансформаторами в первый раз, я не думаю, что они очень хорошо подходят для этой конкретной проблемы. Однако, пока я экспериментировал, мне пришла в голову мысль, которую, я думаю, я рассмотрю в следующем посте: объединив архитектуру Transformer и LSTM, можем ли мы получить лучшее из обоих миров? Большая воспроизводимость, меньшее количество параметров и более точные прогнозы? Следите за новостями. Спасибо за внимание!

Дюны:

Оригинальная статья о трансформерах: https://proceedings.neurips.cc/paper_files/paper/2017/file/3f5ee243547dee91fbd053c1c4a845aa-Paper.pdf

Моя записная книжка для совместной работы с исходным кодом: https://github.com/maym5/lstm_vs_transformer/blob/main/lstm_vs__transformer.ipynb

Документация по Keras Transformer:

https://keras.io/examples/timeseries/timeseries_transformer_classification/

Прогнозирование временных рядов цен с помощью Python

Распространенной задачей в финансах является прогнозирование значений. Существует несколько методов создания прогнозов, таких как ARIMA, байесовские структурные временные ряды и моделирование. Метод Монте-Карло — это метод моделирования, основанный на генерации случайных блужданий за определенный период времени.

Мы начинаем с импорта модулей, которые нам понадобятся. Мы собираемся получить доступ к цене на природный газ из Financial Modeling Prep, веб-API. Вам нужен API-ключ от них, чтобы получить данные.# getting historical data for US Natural Gas prices.
# This code calls the API and transforms the result into a DataFrame.

import numpy as np
np.random.seed(3363)
import pandas as pd

from scipy.stats import norm
import datetime

import matplotlib.pyplot as plt
%matplotlib inlineticker = «NGUSD» #US natural gas
base = ‘https://financialmodelingprep.com/api/v3/’
key = YOUR_KEY
target = «{}historical-price-full/{}?apikey={}».format(base, ticker, key)

df = pd.read_json(target)
df = pd.json_normalize(df[‘historical’])
df[‘date’] = pd.to_datetime(df[‘date’])
df.set_index(‘date’, inplace=True)
# saving the file
df.to_csv(‘data/NGUSD data.csv’)
df.head()

Мы извлекаем данные и сохраняем их в формате CSV, чтобы использовать их позже. Давайте построим график данных, чтобы увидеть, как цена менялась с течением времени.#Plot of asset historical closing price
df[‘adjClose’].plot(figsize=(10, 6), title = «Price of {} from {} to {}».format(ticker, df.index.min(), df.index.max()))

Цена на природный газ в США с июля 2018 года по июль 2023 года

Геометрическое броуновское движение

Мы хотим смоделировать будущие значения природных цен на газ на основе исторических значений. Чтобы сделать это более относительным, мы ограничиваем величину, на которую цена может измениться за период, средним значением исторического изменения цены с помощью некоторого случайного шока. Мы предполагаем, что лог возвратов (процентных изменений) нормально распределен. Мы также исходим из того, что рынок эффективен.

Формула изменения цены между периодами представляет собой цену акции в t0, умноженную на ожидаемый дрейф (среднее изменение цены) плюс экзогенный шок.

Формула для броуновского движения

pred_end_date = datetime.datetime(2023, 7, 20)
forecast_dates = [d if d.isoweekday() in range(1, 6) else np.nan for d in pd.date_range(df.index.max(), pred_end_date)]
intervals = len(forecast_dates)
iterations = 10000
#Preparing log returns from data
log_returns = np.log(1 + df[‘adjClose’].pct_change())

#Setting up drift and random component in relation to asset data
u = log_returns.mean()
var = log_returns.var()
drift = u — (0.5 * var)
stdev = log_returns.std()
daily_returns = np.exp(drift + stdev * norm.ppf(np.random.rand(intervals, iterations)))

#Takes last data point as startpoint point for simulation
S0 = df[‘adjClose’].iloc[-1]
price_list = np.zeros_like(daily_returns)
price_list[0] = S0
#Applies Monte Carlo simulation in asset
for t in range(1, intervals):
price_list[t] = price_list[t — 1] * daily_returns[t]

forecast_df = pd.DataFrame(price_list)

forecast_df.plot(figsize=(10,6), legend=False, title = «{} Simulated Future Paths».format(iterations))

Мы запускаем симуляцию 10 000 раз. Мы могли бы запустить его больше, но пока этого достаточно. Построим график конечных значений. Это ожидаемые значения цены на природный газ за 14 периодов (в данном случае 14 дней).# Plotting with a histogram

x = forecast_df.values[-1]
sigma = np.std(x)
mu = np.mean(x)

num_bins = 15

fig, ax = plt.subplots()

# the histogram of the data
n, bins, patches = ax.hist(x, num_bins, density=1, alpha=.75)

# add a ‘best fit’ line
y = ((1 / (np.sqrt(2 * np.pi) * sigma)) *
np.exp(-0.5 * (1 / sigma * (bins — mu))**2))
ax.plot(bins, y, ‘—‘)
ax.axvline(np.mean(x), color=’r’)
ax.axvline(mu+sigma*1.96, color=’g’, ls=’—‘)
ax.axvline(mu-sigma*1.96, color=’g’, ls=’—‘)
ax.axvline(S0)
ax.set_xlabel(‘Predicted Price on {}’.format(pred_end_date))
ax.set_ylabel(‘Probability density’)
ax.set_title(r’Histogram of {ticker}: $\mu={mu:.02f}$, $\sigma={sigma:.02f}$’.format(ticker = ticker, mu=mu, sigma=sigma))

# Tweak spacing to prevent clipping of ylabel
fig.tight_layout()
plt.show()

Эта гистограмма показывает нам среднее математическое ожидание и диапазон ожидаемых значений через 2 недели. Исходя из этого, мы видим, что можно иметь цену в 5 долларов, но это маловероятно.

Расширенные возможности построения графиков

Фрейм данных, который мы создали для моделирования, содержал каждый прогон в виде отдельного столбца. Это называется wide форматом, который полезен для некоторых аналитических задач, но не оптимален для визуализации. Мы используем технику pd.melt для коллажа из 1000 столбцов в три. Этот новый формат, получивший название long, соответствует принципам tidy data.

Сюжет

Plotly — это высокоуровневая библиотека для интерактивной визуализации. Это «декларативная» библиотека, в которой мы говорим то, что мы хотим, а не то, как производить то, что мы хотим. Обратите внимание, что код для этой визуализации намного меньше, чем код для Matplotlib, и тем не менее результат намного лучше.forecast_df[‘date’] = [df.index.max()+pd.Timedelta(days=i) for i in forecast_df.index]
forecast_df.set_index(‘date’, inplace=True)

df[‘Source’] = ‘Actual’
forecast_df[‘Source’] = ‘Forecast’
result = forecast_df.append(df[[‘adjClose’, ‘Source’]], sort=False)

r = result.reset_index()
r = pd.melt(r, id_vars=[‘Source’, ‘date’])

Авторегрессионная интегрированная скользящая средняя (Arima)

ARIMA — это простой способ прогнозирования будущих значений с помощью статистики. В StatsModels встроен хороший инструмент ARIMA, и мы можем использовать его для прогнозирования будущих значений.from sklearn.metrics import mean_squared_error
from math import sqrt
from statsmodels.tsa.arima.model import ARIMA

# fit ARIMA model
model = ARIMA(df[‘adjClose’], order=(5,1,0))
model_fit = model.fit()
# summary of fit model
print(model_fit.summary())
# line plot of residuals
residuals = pd.DataFrame(model_fit.resid)
residuals.plot()
plt.show()
# density plot of residuals
residuals.plot(kind=’kde’)
plt.show()
# summary stats of residuals
print(residuals.describe())

series = df[‘adjClose’]
series.index = df[‘adjClose’].index
# split into train and test sets
X = series.values
size = int(len(X) * 0.66)
train, test = X[0:size], X[size:len(X)]
history = [x for x in train]
predictions = list()
# walk-forward validation
for t in range(len(test)):
model = ARIMA(history, order=(5,1,0))
model_fit = model.fit()
output = model_fit.forecast()
yhat = output[0]
predictions.append(yhat)
obs = test[t]
history.append(obs)
#print(‘predicted=%f, expected=%f’ % (yhat, obs))
# evaluate forecasts
rmse = sqrt(mean_squared_error(test, predictions))
print(‘Test RMSE: %.3f’ % rmse)
# plot forecasts against actual outcomes
plt.plot(test)
plt.plot(predictions, color=’red’)
plt.show()

Красная линия — это прогнозируемые значения, а синяя (едва видимая) — фактические значения.

Prophet — это байесовский метод временных рядов, реализованный Facebook. Это позволяет нам использовать априорную информацию для прогнозирования будущих значений. На графике показаны фактические значения (точки) и прогнозируемый диапазон (светло-голубой).! pip install -q prophet

from prophet import Prophet
# data must be formated for Prophet

[X, Y] = df.index, df[‘adjClose’]
p = pd.DataFrame(Y, X)
p.reset_index(inplace=True)
p.rename(columns={‘date’: ‘ds’, ‘adjClose’: ‘y’}, inplace=True)

m = Prophet()
m.fit(p)

future = m.make_future_dataframe(periods=5)
future.tail()

forecast = m.predict(future)
forecast[[‘ds’, ‘yhat’, ‘yhat_lower’, ‘yhat_upper’]].tail()
fig1 = m.plot(forecast)

Прогнозирование временных рядов с помощью PyCaret

Здесь объясняется, как выполнять прогнозирование с помощью low-code библиотеки AutoML Python PyCaret.

PyCaret

PyCaret — это библиотека Python для машинного обучения с минимумом программирования, которая автоматизирует рабочие процессы и предлагает эффективные сквозные эксперименты. Автоматизация основных этапов оценки и сравнения моделей облегчает нам изучение производительности различных моделей и выбор наилучшей. На этой странице давайте посмотрим, как подойти к задаче прогнозирования временных рядов знакомых данных авиакомпаний с помощью PyCaret.

Установите и импортируйте PyCaret.

pip install pycaret
import pycaret

Подготовьте данные.

Данные, доступные в библиотеке PyCaret, можно импортировать с помощью pycaret.datasets.get_data(). В противном случае импорт с помощью pandas с устройства как обычно.»’load Airline data»’

import pandas as pd
data = pd.read_excel(‘Airlines_data.xlsx’)
data.set_index(‘Month’,inplace=True)
# data = pycaret.datasets.get_data(‘airlines) # to load from pycaret library
data.head()

Данные авиакомпаний

Среда установки.

Функция setup() инициализирует среду для обучения и создает конвейер преобразования. Важно вызывать setup() перед выполнением любой другой функции. Для эксперимента «Временные ряды» в функцию настройки должны быть переданы следующие параметры.

fh: forecasting horizon

fold: number of folds for cross validation

session_id: random seed for experiment»’Import Time Series libraries»’
from pycaret.time_series import *

»’Initialize the training environment»’
s = setup(data, fh = 3, fold = 5, session_id = 101)

Как видно из приведенной выше таблицы, AutoML анализирует данные временных рядов и предоставляет подробную сводку, особенно о сезонности. Взгляните на записи Seasonality, Present и Significant Seasonal Period(s) в приведенной выше таблице, созданной AutoML.

Выберите API и сравните модели.

PyCaret имеет два API. Функционал и ООП. Функциональный API следует процедурному подходу и предлагает простой интерфейс. Поскольку настройки применяются глобально между функциями, этот API удобен для быстрых экспериментов. API OPP следует объектно-ориентированному подходу, обеспечивая настраиваемый и модульный рабочий процесс и, следовательно, больше подходит для независимого выполнения нескольких экспериментов. В этом эксперименте я использую функциональный API (Коды для ООП API закомментированы).

Функция compare_model() возвращает таблицу, содержащую различные модели и метрики производительности. Из таблицы видно, что экспоненциальное сглаживание превосходит все остальные модели.»’Compare models using Functional API»’
best = compare_models()

# »’OOP API»’
# best = s.compare_models()

Построение графика прогноза и диагностика остатков.

Теперь построим прогноз с помощью функции plot_model() В data_kwargs словаря мы также можем указать ключевые слова для сюжета, такие как title, legend, opacity, height, width и т. д. в соответствии с требованиями.»’Plot forecasting using functional API»’
plot_model(best, plot = ‘forecast’, data_kwargs = {‘fh’ : 24})

# OOP API
# s.plot_model(best, plot = ‘forecast’, data_kwargs = {‘fh’ : 24})

Прогноз модели на основе обучающих данных (в выборке) можно отобразить, как показано ниже.»’Insample plot using functional API»’
plot_model(best, plot = ‘insample’)

# OOP API
# s.plot_model(best, plot = ‘insample’)

Чтобы диагностировать остатки, мы можем использовать ту же функцию plot_model() с ключевым словом plot = ‘diagnostics’, как показано ниже.»’Diagnose best model using functional API»’
plot_model(best, plot = ‘diagnostics’)

# OOP API
# s.plot_model(best, plot = ‘diagnostics’)

Прогноз с использованием наилучшей модели

Поскольку у нас есть экспоненциальное сглаживание в качестве лучшей модели, давайте сделаем прогноз и построим график для горизонта 24.»’Forecast using functional API»’
forecast = predict_model(best, fh = 24)

# OOP API
# s.predict_model(best, fh = 24)

»’Plot forecast»’
import matplotlib.pyplot as plt
plt.plot(forecast.index.strftime(‘%Y-%m’), forecast[‘y_pred’].astype(float),label = ‘forecast’)
plt.legend(loc=[1,1])
plt.xticks(rotation=45);

Сохраните и загрузите лучшую модель.

Модель может быть сохранена с помощью функции save_model() которая сохраняет модель и возвращает архитектуру модели.»’Save model using functional API»’

finalModel = finalize_model(best)
save_model(finalModel, ‘bestModel’)

# OOP API
# s.save_model(finalModel, ‘bestModel’)

Сохраненная модель может быть загружена с помощью функции load_model() следующим образом.»’Load saved model using functional API»’
loaded_model = load_model(‘bestModel’)
# OOP API
# loaded_model = s.load_model(‘my_final_best_model’)

Итоги

На этой странице мы заметили, что PyCaret упрощает такие задачи, как предварительная обработка данных, выбор модели, настройка гиперпараметров и оценка модели, преобразуя рабочий процесс для пользователей. Возможности AutoML значительно сокращают затраты на написание кода, необходимые для разработки надежных моделей машинного обучения, что делает его эффективным и удобным инструментом как для начинающих, так и для опытных специалистов по обработке и анализу данных. Кроме того, PyCaret можно использовать для других задач машинного обучения, таких как регрессия, классификация, кластеризация, обнаружение аномалий и т. д.

Ссылки.

  1. PyCaret 3.0 — Документы (gitbook.io)
  2. https://medium.com/python-in-plain-english/pycaret-simplified-mastering-low-code-machine-learning-e0277faed913

Прогнозирование времени и рядов с помощью рекуррентной нейронной сети LSTM

Каждый день люди делают пассивные прогнозы при выполнении таких задач, как пересечение дороги, где они оценивают скорость автомобилей и расстояние до них, или ловля мяча, угадывая его скорость и соответствующим образом располагая руки. Эти навыки приобретаются с опытом и практикой. Однако прогнозирование сложных явлений, таких как погода или экономика, может быть затруднительным из-за множества задействованных переменных. В таких ситуациях используется прогнозирование по времени и рядам, опираясь на исторические данные и математические модели для прогнозирования будущих тенденций и закономерностей. В этой статье мы рассмотрим пример прогнозирования с математической концепцией с использованием набора данных авиакомпаний.

Часть 1: Математические понятия

В контексте алгоритма прогнозирования временных рядов, используемого в этой статье, вместо ручного вычисления наклона и пересечения линии, алгоритм использует нейронную сеть со слоями LSTM для изучения базовых закономерностей и взаимосвязей в данных временных рядов. Нейронная сеть обучается на части данных, а затем используется для прогнозирования оставшейся части. В этом алгоритме предсказание следующего временного шага основано на предыдущих n_inputs временных шагов, что аналогично концепции использования y(t) для предсказания y(T+1) в примере линейной регрессии. Однако вместо использования простого линейного уравнения предсказание в этом алгоритме генерируется с помощью функции активации слоя LSTM. Функция активации позволяет модели фиксировать нелинейные зависимости в данных, что делает ее более эффективной при регистрации сложных закономерностей в данных временных рядов.

Функция активации

Функция активации, используемая в модели LSTM, является функцией активации выпрямленного линейного блока (ReLU). Эта функция активации обычно используется в моделях глубокого обучения из-за ее простоты и эффективности при решении проблемы исчезающего градиента. В модели LSTM функция активации ReLU применяется к выходу каждого блока LSTM, чтобы внести нелинейность в модель и позволить ей изучать сложные закономерности в данных. Функция ReLU имеет простое пороговое поведение, при котором любой отрицательный входной сигнал сопоставляется с нулем, а любой положительный вход пропускается без изменений, что делает ее вычислительно эффективной.

Часть 2 : Реализация

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
df = pd.read_csv(‘airline-passengers.csv’, index_col=’Month’, parse_dates=True)
df.index.freq = ‘MS’
df.shape
df.columns
plt.figure(figsize=(20, 4))
plt.plot(df.Passengers, linewidth=2)
plt.show()

Код импортирует три важные библиотеки: numpy, pandas и matplotlib. Библиотека pandas используется для чтения файла «airline-passengers.csv» и установки столбца «Месяц» в качестве индекса, что позволяет анализировать данные с течением времени. Затем код использует библиотеку matplotlib для создания линейного графика, показывающего количество авиапассажиров с течением времени. Наконец, график отображается с помощью функции ‘plt.show’. Этот код полезен для всех, кто интересуется анализом данных временных рядов, и он демонстрирует, как использовать pandas и matplotlib для визуализации тенденций в данных.

nobs = 12
df_train = df.iloc[:-nobs]
df_test = df.iloc[-nobs:]
df_train.shape
df_test.shape

Этот код создает два новых фрейма данных «df_train» и «df_test» путем разделения существующего фрейма данных временных рядов «df» на обучающий и проверочный наборы. Переменная ‘nobs’ установлена в 12, что означает, что последние 12 наблюдений ‘df’ будут использоваться для тестирования, а остальные данные будут использоваться для обучения. Обучающий набор хранится в ‘df_train’ и состоит из всех строк ‘df’, за исключением последних 12 строк, в то время как проверочный набор хранится в ‘df_test’ и состоит только из последних 12 строк ‘df’. Атрибут shape затем используется для печати количества строк и столбцов в каждом фрейме данных, что подтверждает, что разбиение было выполнено правильно. Этот код полезен для подготовки данных временных рядов для моделирования и тестирования путем разделения их на два набора.

Архитектура модели

Фото предоставлено @learnwithutsav

from keras.preprocessing.sequence import TimeseriesGenerator
from sklearn.preprocessing import MinMaxScaler
scaler = MinMaxScaler()
scaler.fit(df_train)
scaled_train = scaler.transform(df_train)
scaled_test = scaler.transform(df_test)n_inputs = 12
n_features = 1
generator = TimeseriesGenerator(scaled_train, scaled_train, length = n_inputs, batch_size =1)

for i in range(len(generator)):
X, y = generator[i]
print(f’ \n {X.flatten()} and {y}’)

В этом фрагменте кода показано, как использовать класс ‘TimeseriesGenerator’ из Keras и класс ‘MinMaxScaler’ из scikit-learn для создания входных и выходных массивов для модели прогнозирования временных рядов. Сначала код создает экземпляр класса ‘MinMaxScaler’ и подгоняет его к набору обучающих данных (‘df_train’) для масштабирования данных. Затем масштабированные данные хранятся во фреймах данных scaled_train и scaled_test. Количество временных шагов (‘n_inputs’) установлено равным 12, а количество объектов (‘n_features’) установлено равным 1. Создается объект ‘TimeseriesGenerator’ с данными ‘scaled_train’, длиной окна ‘n_inputs’ и размером пакета 1. Наконец, цикл используется для перебора объекта ‘generator’ и вывода входного и выходного массивов для каждого временного шага. Переменные ‘X’ и ‘y’ представляют входной и выходной массивы для каждого временного шага соответственно. Метод ‘flatten()’ используется для преобразования входного массива в 1D-массив для упрощения печати. В целом, этот код полезен для подготовки данных временных рядов для моделей прогнозирования, использующих подход скользящего окна.

X.shape

Этот код возвращает форму массива или матрицы ‘X’. Атрибут ‘shape’ является свойством массивов NumPy и возвращает кортеж, представляющий размерности массива. Код не предоставляет никакого дополнительного контекста, поэтому неясно, какова форма ‘X’. Вывод будет в формате (строки, столбцы).

from keras.models import Sequential
from keras.layers import Dense
from keras.layers import LSTM

model = Sequential()
model.add(LSTM(200, activation=’relu’, input_shape = (n_inputs, n_features)))
model.add(Dense(1))
model.compile(optimizer=’adam’, loss=’mse’)

model.summary()

Этот код демонстрирует, как создать модель нейронной сети LSTM для прогнозирования временных рядов с помощью Keras. Во-первых, импортируются необходимые классы Keras, в том числе ‘Sequential’, ‘Dense’ и ‘LSTM’. Модель создается как «последовательный» объект, и добавляется слой LSTM с 200 нейронами, функцией активации «relu» и входной формой, определяемой «n_inputs» и «n_features». Затем выходные данные слоя LSTM передаются в слой «Dense» с одним выходным нейроном. Модель скомпилирована с оптимизатором adam и функцией потерь среднеквадратичной ошибки (mse). Метод ‘summary()’ используется для отображения сводной информации об архитектуре, включая количество параметров и формы входного и выходного тензоров для каждого слоя. Этот код может быть полезен для создания LSTM-модели для прогнозирования временных рядов, так как он представляет собой простой пример, который можно адаптировать к различным наборам данных и задачам прогнозирования.

Этап обучения

model.fit(generator, epochs = 50)

Этот код обучает модель нейронной сети LSTM с помощью метода fit()’ в Keras для 50 эпох. Объект ‘TimeseriesGenerator’ генерирует пакеты пар ввода-вывода, на основе которых модель будет учиться. Метод fit()’ обновляет параметры модели с помощью обратного распространения ошибки на основе функции потерь и оптимизатора, определенных при компиляции модели. Обучая модель, она учится делать прогнозы на основе новых, невидимых данных на основе закономерностей, полученных в обучающих данных.

plt.plot(model.history.history[‘loss’])

last_train_batch = scaled_train[-12:]

last_train_batch = last_train_batch.reshape(1, 12, 1)

last_train_batch

model.predict(last_train_batch)

Этот код использует обученную модель нейронной сети LSTM для прогнозирования новой точки данных. Последние 12 точек данных из обучающих данных выбираются, масштабируются и преобразуются в формат, соответствующий модели. Метод predict()’ вызывается для модели с измененными данными в качестве входных данных, а выходными данными является прогнозируемое значение для следующего временного шага во временном ряду. Это важный шаг в использовании модели LSTM для прогнозирования временных рядов.

scaled_test[0]

Этот код печатает первый элемент масштабируемого массива тестовых данных. Переменная ‘scaled_test’ представляет собой массив NumPy тестовых данных, которые были преобразованы с помощью объекта ‘MinMaxScaler’. При печати первого элемента этого массива отображается масштабированное значение для первого временного шага в тестовых данных.

Прогнозирование

y_pred = []

first_batch = scaled_train[-n_inputs:]
current_batch = first_batch.reshape(1, n_inputs, n_features)

for i in range(len(scaled_test)):
batch = current_batch
pred = model.predict(batch)[0]
y_pred.append(pred)
current_batch = np.append(current_batch[:,1:, :], [[pred]], axis = 1)

y_pred

scaled_test

Этот код генерирует прогнозы для тестовых данных с помощью обученной модели LSTM. Он использует цикл for для перебора каждого элемента в масштабируемых тестовых данных. В каждой итерации текущий пакет используется для составления прогноза с помощью метода ‘predict()’ модели. Затем прогнозируемое значение добавляется в список y_pred и текущий пакет обновляется. Наконец, список «y_pred» печатается вместе с данными «scaled_test» для сравнения прогнозируемых значений с фактическими значениями. Этот шаг имеет решающее значение для оценки производительности LSTM-модели на тестовых данных.

df_test

y_pred_transformed = scaler.inverse_transform(y_pred)

y_pred_transformed = np.round(y_pred_transformed,0)

y_pred_final = y_pred_transformed.astype(int)

y_pred_final

Этот код преобразует прогнозируемые значения, сгенерированные на предыдущем шаге, обратно в исходную шкалу с помощью метода ‘inverse_transform()’ объекта scaler. Преобразованные значения округляются до ближайшего целого числа с помощью функции ’round()’ и преобразуются в целые числа с помощью метода ‘astype()’. Результирующий массив прогнозируемых значений, ‘y_pred_final’, печатается для отображения окончательных прогнозируемых значений для тестовых данных. Этот шаг важен для оценки точности предсказаний LSTM-модели в исходном масштабе данных.

df_test.values, y_pred_final

df_test[‘Predictions’] = y_pred_final

df_test

В приведенном выше коде показаны прогнозируемые значения, сгенерированные LSTM-моделью, добавляемые к исходному тестовому набору данных. Во-первых, атрибут ‘values’ используется для извлечения значений кадра данных ‘df_test’, которые затем объединяются с прогнозируемыми значениями ‘y_pred_final’. Затем в кадр данных df_test» добавляется новый столбец под названием «Прогнозы» для хранения прогнозируемых значений. Наконец, кадр данных «df_test» печатается вместе с недавно добавленным столбцом «Прогнозы». Этот шаг важен для визуального сравнения фактических значений тестового набора данных с прогнозируемыми значениями и оценки точности модели.

plt.figure(figsize=(15, 6))
plt.plot(df_train.index, df_train.Passengers, linewidth=2, color=’black’, label=’Train Values’)
plt.plot(df_test.index, df_test.Passengers, linewidth=2, color=’green’, label=’True Values’)
plt.plot(df_test.index, df_test.Predictions, linewidth=2, color=’red’, label=’Predicted Values’)
plt.legend()
plt.show()

Этот блок кода генерирует график с помощью библиотеки matplotlib. Сначала он задает размер фигуры, а затем отображает обучающие данные в виде черной линии, истинные тестовые значения — в виде зеленой линии, а прогнозируемые тестовые значения — в виде красной линии. Он также добавляет легенду к графику и отображает ее с помощью метода show()

Среднеквадратичная ошибка

Среднеквадратичная ошибка (MSE) — это мера того, насколько близко линия регрессии находится к набору точек. Он вычисляется путем взятия среднего квадрата разностей между прогнозируемым и фактическим значениями. Квадратный корень из MSE известен как среднеквадратичная ошибка (RMSE), которая является популярной мерой точности прогнозов. В этом блоке кода RMSE вычисляется с помощью функции mean_squared_error из модуля sklearn.metrics и функции sqrt из модуля math. RMSE используется для оценки точности прогнозов модели LSTM по сравнению с истинными значениями в тестовом наборе.from sklearn.metrics import mean_squared_error
from math import sqrt

sqrt(mean_squared_error(df_test.Passengers, df_test.Predictions))

Этот код вычисляет среднеквадратичную ошибку (RMSE) между фактическими пассажирскими значениями в тестовом наборе (df_test.Passengers) и прогнозируемые значения пассажиропотока (df_test.Predictions). RMSE — это широко используемая метрика для оценки производительности регрессионных моделей. Он измеряет среднее расстояние между прогнозируемыми значениями и фактическими значениями с учетом квадрата разностей между ними. RMSE является полезной метрикой, поскольку она наказывает за большие ошибки сильнее, чем за маленькие, что делает ее хорошим индикатором общей точности прогнозов модели.

В заключение мы реализовали модель прогнозирования временных рядов с использованием алгоритма LSTM в Keras. Мы обучили модель на ежемесячном наборе данных о пассажирах авиакомпаний и использовали ее для составления прогнозов на следующие 12 месяцев. Модель показала хорошие результаты со среднеквадратичной ошибкой 30,5. Визуализация истинных, прогнозируемых и обучающих значений показала, что модель способна уловить общую тенденцию и сезонность в данных. Это демонстрирует возможности LSTM в регистрации сложных временных взаимосвязей в данных временных рядов и его потенциал для точных прогнозов.

Исследовательский анализ данных

Часть 1

В этом разделе я расскажу о различных методах, используемых для анализа данных временных рядов с целью их подготовки к дальнейшему прогнозированию с использованием библиотек python. Цель состоит в том, чтобы подготовить данные временных рядов, которые полезны для алгоритмов прогнозирования временных рядов для прогнозирования производства ветровой энергии. Алгоритмы прогнозирования не являются частью этой статьи, я напишу для них отдельный блог.

Что такое анализ данных временных рядов?

Анализ данных временных рядов — это статистический метод, используемый для анализа точек данных, собранных с течением времени. Он включает в себя изучение закономерностей и тенденций данных для выявления любых основных закономерностей или связей между точками данных.

При анализе временных рядов данные обычно собираются через регулярные промежутки времени, например ежедневно, еженедельно или ежемесячно, а точки данных располагаются в хронологическом порядке. Анализ обычно включает в себя изучение статистических свойств данных, таких как среднее значение, дисперсия и корреляция, для выявления тенденций, закономерностей и взаимосвязей с течением времени.

Анализ временных рядов широко используется в различных областях, включая экономику, финансы, прогнозирование погоды, обработку сигналов и инженерию, и это лишь некоторые из них. Это помогает прогнозировать будущие тенденции и закономерности в данных и помогает принимать обоснованные решения на основе анализа.

Что такое исследовательский анализ данных (EDA)?

Исследовательский анализ данных (EDA) — это процесс анализа и обобщения наборов данных для получения аналитических сведений и выявления закономерностей, связей, аномалий, а также потенциальных проблем или выбросов. Это важный шаг в науке о данных, машинном обучении и статистическом анализе, который помогает лучше понять данные и сформулировать гипотезы для дальнейшего исследования.

Целью САПР является изучение данных, определение важных признаков и характеристик, а также создание визуализаций и статистических сводок, которые раскрывают значимую информацию. Этот процесс включает в себя различные методы, такие как гистограммы, диаграммы рассеяния, ящичковые диаграммы, тепловые карты, корреляционные матрицы и другие исследовательские графики для исследования данных.

EDA также используется для выявления проблем с качеством данных, таких как отсутствующие или несогласованные данные, выбросы и аномалии, и для разработки стратегий по их устранению. В целом, EDA помогает гарантировать, что данные подходят для предполагаемого анализа, а результаты являются точными и надежными.

После того, как вы завершили EDA, вам необходимо подготовить данные для дальнейшего анализа и прогнозирования. Подготовка данных является важным этапом анализа данных временных рядов, поскольку данные временных рядов часто являются зашумленными и несогласованными и могут содержать отсутствующие или ошибочные значения. Подготовка данных включает в себя очистку, преобразование и нормализацию данных, чтобы сделать их пригодными для анализа.

  1. Качество данных: Данные временных рядов могут содержать ошибки и пропущенные значения, что может повлиять на точность анализа. Подготовка данных помогает выявить и удалить такие ошибки и пропущенные значения из данных.
  2. Согласованность данных: Данные временных рядов собираются через разные промежутки времени и могут иметь несоответствия из-за изменений в измерительных приборах или методах. Подготовка данных обеспечивает согласованность данных и возможность их сравнения во времени.
  3. Выбор функции: Данные временных рядов могут содержать множество переменных или объектов, которые не имеют отношения к анализу. Подготовка данных помогает в выборе релевантных признаков для анализа.
  4. Нормализация: Данные временных рядов могут иметь разные масштабы и диапазоны для разных переменных, что затрудняет их сравнение. Подготовка данных нормализует данные, делая их сопоставимыми.

В целом, подготовка данных необходима для анализа данных временных рядов, чтобы гарантировать, что данные являются точными, согласованными и актуальными для анализа и прогнозирования.

Набор данных, используемый для САПР

Для анализа я использовал набор данных о ветроэнергетике от Kaggle. Этот набор данных содержит различные характеристики погоды, турбины и ротора. Данные были зафиксированы с января 2018 года по март 2020 года. Показания записывались с интервалом в 10 минут.

Определимся с целью данного анализа. Целью данного анализа является прогнозирование выработки ветровой энергии на еженедельной основе.#Import the necessary libraries
import pandas as pd
import numpy as np
import seaborn as sns
import matplotlib.pyplot as plt

#Load the dataset
data = pd.read_csv(«Turbine_Data.csv»)

В качестве первого шага выполните следующие команды для загруженных данных. Целью этого шага является ознакомление с набором данных. Внимательно изучите данные, проверьте столбцы и строки. Обратите внимание на пропущенные значения, типы данных и т.д.# show the first few rows of the dataset
data.head()
# show the last few rows of the dataset
data.tail()
# show information about the dataset such as data types and missing values
data.info()
# show statistical summaries of the dataset
data.describe()

Имя первого столбца, содержащего метку времени, пусто. Давайте переименуем его в «Date»# rename first column to «Date»
data.columns.values[0]=’Date’

Проверьте выбросы

Построим ящичковую диаграмму для значений выработки электроэнергии для всех месяцев в наборе данныхsns.set(rc={‘figure.figsize’:(14,10)})
sns.boxplot(x=data.index.month_name(), y=’ActivePower’, data=data, palette=’muted’)
plt.ylabel(‘ActivePower’)
plt.xlabel(‘Month’)
plt.title(«Power Generated»)

Ниже приведена сгенерированная блочная диаграмма.

Вы можете заметить, что на приведенной выше ящичковой диаграмме есть отрицательные значения активной мощности, что не должно быть возможным. Давайте проверим это дальше#Print ActivePower negative values
data[‘ActivePower’][data[‘ActivePower’]<0]

Ниже приведен выводDate
2018-01-01 00:00:00+00:00 -5.357727
2018-01-01 00:10:00+00:00 -5.822360
2018-01-01 00:20:00+00:00 -5.279409
2018-01-01 00:30:00+00:00 -4.648054
2018-01-01 00:40:00+00:00 -4.684632

2020-03-30 03:50:00+00:00 -7.005695
2020-03-30 04:00:00+00:00 -5.576951
2020-03-30 04:10:00+00:00 -4.945515
2020-03-30 04:20:00+00:00 -6.565684
2020-03-30 04:40:00+00:00 -9.814717
Name: ActivePower, Length: 15644, dtype: float64

Отрицательная генерация электроэнергии не имеет смысла. Отбросим строки с отрицательными значениями выработки электроэнергии#Filtering rows with postive power generation values
data = data[data[‘ActivePower’]>=0]

Проверка уникальных значений

Начнем с проверки уникальных значений в наборе данных.#Check unique values
data.nunique()

Ниже приведен вывод. Обратите внимание, что есть два столбца, которые имеют значение 1, что означает, что они имеют постоянное значение во всех строках.ActivePower 78449
AmbientTemperatue 77943
BearingShaftTemperature 52091
Blade1PitchAngle 33383
Blade2PitchAngle 33431
Blade3PitchAngle 33431
ControlBoxTemperature 1
GearboxBearingTemperature 52099
GearboxOilTemperature 52171
GeneratorRPM 51730
GeneratorWinding1Temperature 52188
GeneratorWinding2Temperature 52191
HubTemperature 31207
MainBoxTemperature 41376
NacellePosition 4990
ReactivePower 78413
RotorRPM 51641
TurbineStatus 121
WTG 1
WindDirection 4990
WindSpeed 78497

Обычно постоянные величины не имеют никакого значения в нашем анализе. Поэтому отбросим столбцы с постоянным значением, а именно — ‘WTG’ и ‘ControlBoxTemperature’# drop columns with constant values
data = data.drop([‘WTG’, ‘ControlBoxTemperature’], axis=1)
data.nunique()

Давайте быстро проверим вывод ниже после удаления двух столбцовActivePower 78449
AmbientTemperatue 77943
BearingShaftTemperature 52091
Blade1PitchAngle 33383
Blade2PitchAngle 33431
Blade3PitchAngle 33431
GearboxBearingTemperature 52099
GearboxOilTemperature 52171
GeneratorRPM 51730
GeneratorWinding1Temperature 52188
GeneratorWinding2Temperature 52191
HubTemperature 31207
MainBoxTemperature 41376
NacellePosition 4990
ReactivePower 78413
RotorRPM 51641
TurbineStatus 121
WindDirection 4990
WindSpeed 78497

Проверка пропущенных значений

Столбцы могут содержать нулевые значения, которые необходимо будет обработать. Давайте проверим нулевые значения в наборе данных.# count the number of missing values in each column
data.isnull().sum()

Ниже приведен вывод. Обратите внимание, что в каждом столбце очень много нулевых или отсутствующих значений.ActivePower 0
AmbientTemperatue 1032
BearingShaftTemperature 26932
Blade1PitchAngle 43397
Blade2PitchAngle 43480
Blade3PitchAngle 43480
GearboxBearingTemperature 26930
GearboxOilTemperature 26915
GeneratorRPM 26919
GeneratorWinding1Temperature 26901
GeneratorWinding2Temperature 26894
HubTemperature 27041
MainBoxTemperature 26952
NacellePosition 20429
ReactivePower 42
RotorRPM 26925
TurbineStatus 26577
WindDirection 20429
WindSpeed 308

Обработка пропущенных значений с помощью линейной интерполяции

Линейная интерполяция — это метод, используемый для оценки значения переменной в заданной точке на основе значений переменной в окружающих точках. Предполагается, что переменная изменяется с постоянной скоростью между соседними точками данных.

Применим линейную интерполяцию к отсутствующим точкам данных.## Linear Interpolation Method
data=data.interpolate(method =’linear’, limit_direction =’backward’)
data.isnull().sum()

Ниже приведен вывод после успешной линейной интерполяции пропущенных значений.Date 0
ActivePower 0
AmbientTemperatue 0
BearingShaftTemperature 0
Blade1PitchAngle 0
Blade2PitchAngle 0
Blade3PitchAngle 0
GearboxBearingTemperature 0
GearboxOilTemperature 0
GeneratorRPM 0
GeneratorWinding1Temperature 0
GeneratorWinding2Temperature 0
HubTemperature 0
MainBoxTemperature 0
NacellePosition 0
ReactivePower 0
RotorRPM 0
TurbineStatus 0
WindDirection 0
WindSpeed 0

Теперь, когда пропущенные значения обработаны, мы можем приступить к построению графика данных. Нас интересует переменная «Active Power», так как она содержит мощность, генерируемую с заданной меткой времени. Давайте быстро построим график зависимости генерируемой энергии от времени.

Построим график активной мощности# plot ‘ActivePower’ for the dates
data.reset_index().plot(x = ‘Date’, y = ‘ActivePower’)

Приведенный выше график неясен, так как в нем задействовано слишком много точек данных. Поскольку выборка данных выполняется каждые 10 минут, точек данных слишком много. Нас интересуют еженедельные данные. Там мы будем пересчитывать данные на ежедневной основе, свертывая их по «среднему» значению точек данных.# resample the data on daily frequency
data = data.resample(«D»).mean()
data.reset_index().plot(x = «Date», y = ‘ActivePower’)
plt.legend(loc=’best’)
plt.title(‘Turbine Data’)
plt.show(block=False)

Ниже приведен вывод ресемплинга. Довольно чисто, правда?

Корреляция переменных

Корреляция относится к степени, в которой две или более переменных связаны или связаны друг с другом. Корреляция — это статистическая мера, которая описывает силу и направление линейной связи между двумя переменными.# visualize the degree of correlation between variables
sns.heatmap(data.corr(), annot=True)

Ниже приведен результат, отображающий тепловую карту с коэффициентом корреляции Пирсона.

На приведенном выше графике коэффициент корреляции — это число в каждой ячейке. Давайте посмотрим на интересующую нас переменную «ActivePower». Учитывайте величину коэффициента корреляции.

Как правило, коэффициенты корреляции выше 0,5 или ниже -0,5 считаются сильными, в то время как коэффициенты от 0,3 до 0,5 или от -0,3 до -0,5 считаются умеренными.

Уберем переменные с низким коэффициентом, т.е. от 0.5 до -0.5# dropping columns with low coeffecient
data = data.drop([‘AmbientTemperatue’, ‘TurbineStatus’, ‘WindDirection’, ‘NacellePosition’, ‘MainBoxTemperature’,’HubTemperature’,’Blade3PitchAngle’,’Blade2PitchAngle’, ‘Blade1PitchAngle’], axis=1)

Теперь у нас есть подготовленные данные, которые мы можем использовать для прогнозирования. Мы можем нарисовать парную диаграмму в качестве последнего шага для понимания попарных отношений между переменными в наборе данных.

Многомерный анализ

Давайте создадим парную диаграмму, показывающую точечные диаграммы всех возможных пар переменных в DataFrame, с точками, раскрашенными в зависимости от значения столбца «ActivePower». Кроме того, диагональные подграфики будут отображать гистограммы каждой переменной в DataFrame.# creating a pairplot for ‘ActivePower’ values with remaining columns
sns.pairplot(data, hue=’ActivePower’, diag_kind=’hist’)
plt.show()

Из приведенного выше анализа можно сделать вывод, что выработка электроэнергии находится в прямой зависимости от скорости ветра и оборотов генератора.

На этом базовый набор данных САПР по ветроэнергетике завершен. Существует множество библиотек, которые при необходимости можно использовать для дальнейшего анализа.

Часть 2

В этой статье я проведу исследовательский анализ данных и анализ временных рядов с использованием набора данных, состоящего из продаж товаров в разных категориях из магазина в США в период с 2015 по 2018 год. Выполнив эти шаги, я проведу предварительную обработку данных и проведу многомерное прогнозное исследование.

Pixabay, Pexels

Здравствуйте, сегодня я в восторге от того, что наконец-то сяду перед своим компьютером и углубись в захватывающий проект по науке о данных, который я откладывал в течение довольно долгого времени. В этой статье я поделюсь первой частью своего двухэтапного проекта. С помощью исследовательского анализа данных и анализа временных рядов я отправлюсь в увлекательное путешествие в царство моего набора данных. В этой статье я раскрою секреты тенденций продаж в разных категориях товаров. Визуализируя данные, мы станем свидетелями колебаний показателей продаж, а затем углубимся в сложную динамику этих закономерностей с течением времени.

Я использовал библиотеку Plotly в качестве инструмента визуализации для получения ценной информации из набора данных. Plotly оказался чрезвычайно полезным при создании интерактивных графиков для визуализации данных. Я настоятельно рекомендую использовать эту библиотеку для визуализации данных. Учитывая, что мой набор данных является набором данных временных рядов, эта статья включает раздел исследовательского анализа данных, охватывающий различные аспекты анализа временных рядов, такие как стационарность, сезонность и цикличность.

Теперь я могу начать с импорта библиотек, которые я буду использовать в исследовании.import numpy as np
import pandas as pd
import seaborn as sns
from datetime import datetime as dt

import plotly.express as px
import plotly.graph_objects as go

import itertools

import statsmodels.api as sm
from statsmodels.tsa.seasonal import seasonal_decompose
import statsmodels.tsa.api as smt

import warnings
warnings.filterwarnings(‘ignore’)
#pd.set_option(‘display.max_columns’, None)
pd.set_option(‘display.max_rows’, None)
pd.set_option(‘display.width’, None)
pd.set_option(‘display.float_format’, lambda x: ‘%.f’ % x)
pd.set_option(‘display.expand_frame_repr’, False)

После импорта необходимых библиотек я перехожу к импорту набора данных, который буду использовать для этого проекта. Работа, представленная в этой статье, основана на моей записной книжке, опубликованной на Kaggle, где я поделился своими выводами и анализом. Используя возможности методов обработки и анализа данных, я стремлюсь предложить ценную информацию и проиллюстрировать пошаговый процесс, связанный с этим проектом. Набор данных служит ценным ресурсом для демонстрации практического применения различных методов анализа данных и прогнозирования.dataset = pd.read_csv(‘/kaggle/input/sales-forecasting/train.csv’)
df = dataset.copy()

Итак, после того, как я импортировал набор данных, я начинаю свою работу. На первом шаге я определяю функцию с именем ‘check_df’. Эта функция очень ценна для понимания общей структуры набора данных. Он предлагает доступ к различной информации, включая измерения набора данных, имена объектов, типы данных, а также возможность просмотра первой и последней 5 строк набора данных. Кроме того, я стремлюсь выявить ценные закономерности и тенденции в наборе данных, анализируя статистические свойства числовых переменныхdef check_df(dataset, head = 5):
print(‘#’*30 + ‘Shape of Dataset’ + ‘#’*30, end = ‘\n’*2)
print(dataset.shape, end = ‘\n’*2)
print(‘#’*30 + ‘General informations about to Dataset’ + ‘#’*30, end = ‘\n’*2)
print(dataset.info(), end = ‘\n’*2)
print(‘#’*30 + ‘First 5 Lines Of Dataset’ + ‘#’*30, end = ‘\n’*2)
print(dataset.head(head), end = ‘\n’*2)
print(‘#’*30 + ‘NaN values of Features’ + ‘#’*30, end = ‘\n’*2)
print(dataset.isnull().sum(), end = ‘\n’*2)
print(‘#’*30 + ‘Descriptive Statistics of Numerical Features’ + ‘#’*30, end = ‘\n’*2)
print(dataset.describe().T, end = ‘\n’*2)
print(‘#’*30 + ‘Quantiles of Numerical Features’ + ‘#’*30, end =’\n’*2)
print(dataset.describe([0,0.10, 0.25, 0.50,0.75,0.99]).T, end = ‘\n’*2)

check_df(df)

Когда я запускаю функцию с именем check_df, которую я создал, общий вид набора данных выглядит следующим образом:##############################Shape of Dataset##############################

(9800, 18)

##############################General informations about to Dataset##############################

<class ‘pandas.core.frame.DataFrame’>
RangeIndex: 9800 entries, 0 to 9799
Data columns (total 18 columns):
# Column Non-Null Count Dtype
— —— ————— ——
0 Row ID 9800 non-null int64
1 Order ID 9800 non-null object
2 Order Date 9800 non-null object
3 Ship Date 9800 non-null object
4 Ship Mode 9800 non-null object
5 Customer ID 9800 non-null object
6 Customer Name 9800 non-null object
7 Segment 9800 non-null object
8 Country 9800 non-null object
9 City 9800 non-null object
10 State 9800 non-null object
11 Postal Code 9789 non-null float64
12 Region 9800 non-null object
13 Product ID 9800 non-null object
14 Category 9800 non-null object
15 Sub-Category 9800 non-null object
16 Product Name 9800 non-null object
17 Sales 9800 non-null float64
dtypes: float64(2), int64(1), object(15)
memory usage: 1.3+ MB
None

##############################First 5 Lines Of Dataset##############################

Row ID Order ID Order Date Ship Date Ship Mode Customer ID Customer Name Segment Country City State Postal Code Region Product ID Category Sub-Category Product Name Sales
0 1 CA-2017-152156 08/11/2017 11/11/2017 Second Class CG-12520 Claire Gute Consumer United States Henderson Kentucky 42420 South FUR-BO-10001798 Furniture Bookcases Bush Somerset Collection Bookcase 262
1 2 CA-2017-152156 08/11/2017 11/11/2017 Second Class CG-12520 Claire Gute Consumer United States Henderson Kentucky 42420 South FUR-CH-10000454 Furniture Chairs Hon Deluxe Fabric Upholstered Stacking Chairs,… 732
2 3 CA-2017-138688 12/06/2017 16/06/2017 Second Class DV-13045 Darrin Van Huff Corporate United States Los Angeles California 90036 West OFF-LA-10000240 Office Supplies Labels Self-Adhesive Address Labels for Typewriters b… 15
3 4 US-2016-108966 11/10/2016 18/10/2016 Standard Class SO-20335 Sean O’Donnell Consumer United States Fort Lauderdale Florida 33311 South FUR-TA-10000577 Furniture Tables Bretford CR4500 Series Slim Rectangular Table 958
4 5 US-2016-108966 11/10/2016 18/10/2016 Standard Class SO-20335 Sean O’Donnell Consumer United States Fort Lauderdale Florida 33311 South OFF-ST-10000760 Office Supplies Storage Eldon Fold ‘N Roll Cart System 22

##############################NaN values of Features##############################

Row ID 0
Order ID 0
Order Date 0
Ship Date 0
Ship Mode 0
Customer ID 0
Customer Name 0
Segment 0
Country 0
City 0
State 0
Postal Code 11
Region 0
Product ID 0
Category 0
Sub-Category 0
Product Name 0
Sales 0
dtype: int64

##############################Descriptive Statistics of Numerical Features##############################

count mean std min 25% 50% 75% max
Row ID 9800 4900 2829 1 2451 4900 7350 9800
Postal Code 9789 55273 32041 1040 23223 58103 90008 99301
Sales 9800 231 627 0 17 54 211 22638

##############################Quantiles of Numerical Features##############################

count mean std min 0% 10% 25% 50% 75% 99% max
Row ID 9800 4900 2829 1 1 981 2451 4900 7350 9702 9800
Postal Code 9789 55273 32041 1040 1040 10024 23223 58103 90008 98115 99301
Sales 9800 231 627 0 0 8 17 54 211 2480 22638

Когда мы наблюдаем общую структуру датасета, мы видим, что мы работаем с датастом, содержащим 9800 строк наблюдений и 18 переменных. Мы также обнаруживаем отсутствующие наблюдения в переменной с именем «Почтовый индекс». Изучая недостающие наблюдения, мы стремимся разработать стратегии того, как восполнить эти ценности. Кроме того, в нашем наборе данных есть еще несколько переменных, которые, хотя и имеют числовые значения, на самом деле являются категориальными переменными в соответствии с содержащейся в них информацией. Например, столбцы «Почтовый индекс» и «Идентификатор строки» представлены в числовом виде, хотя на самом деле они имеют категориальный атрибут. Мы продолжим преобразование этих переменных в категориальные переменные.df.loc[(df[«Postal Code»].isnull()), :][[«City», «State»]].value_counts()

## Output ##

City State
Burlington Vermont 11
dtype: int64

Все 11 строк отсутствующих наблюдений в наборе данных относятся к одной и той же группе State-City. С помощью простого поиска в Google я нашел почтовый индекс для этого региона и заполнил эти значения для недостающих наблюдений. Для того, чтобы продолжить нашу работу, нам нужно изменить типы некоторых переменных. Столбцы ‘Дата отгрузки’ и ‘Дата заказа’ рассматриваются как ‘object’ как типы данных, но поскольку эти две переменные на самом деле содержат информацию о дате, я хочу преобразовать их в формат ‘datetime’.df.drop(columns = [‘Row ID’], inplace = True)

df.loc[(df[«Postal Code»].isnull()), «Postal Code»] = 5401

df[«Postal Code»] = df[«Postal Code»].astype(«str»)

def convert_datetime(dataframe):
date_cols = [col for col in dataframe.columns if ‘Date’ in col]

for col in date_cols:
dataframe[col] = pd.to_datetime(dataframe[col], errors = ‘coerce’)
return dataframe.info()

convert_datetime(df)

## Output ## :

<class ‘pandas.core.frame.DataFrame’>
RangeIndex: 9800 entries, 0 to 9799
Data columns (total 17 columns):
# Column Non-Null Count Dtype
— —— ————— ——
0 Order ID 9800 non-null object
1 Order Date 9800 non-null datetime64[ns]
2 Ship Date 9800 non-null datetime64[ns]
3 Ship Mode 9800 non-null object
4 Customer ID 9800 non-null object
5 Customer Name 9800 non-null object
6 Segment 9800 non-null object
7 Country 9800 non-null object
8 City 9800 non-null object
9 State 9800 non-null object
10 Postal Code 9800 non-null object
11 Region 9800 non-null object
12 Product ID 9800 non-null object
13 Category 9800 non-null object
14 Sub-Category 9800 non-null object
15 Product Name 9800 non-null object
16 Sales 9800 non-null float64
dtypes: datetime64[ns](2), float64(1), object(14)
memory usage: 1.3+ MB

Мы исправили отсутствующие наблюдения в наборе данных и отредактировали неправильный тип данных. После преобразования переменных Date в наборе данных из типа данных object в тип данных Datetime следующим шагом является категоризация переменных в наборе данных по различным типам, включая категориальные, числовые и т. д. Следуя этой категоризации, я буду использовать различные визуализации для изучения частот классов категориальных переменных.def grab_col_name(dataframe, cat_th = 5, car_th = 49):
cat_cols = [col for col in dataframe.columns if dataframe[col].dtypes == «O»]
num_but_cat = [col for col in dataframe.columns if dataframe[col].dtypes != «O» and dataframe[col].nunique() < cat_th]
cat_but_car = [col for col in dataframe.columns if dataframe[col].dtypes == «O» and dataframe[col].nunique() > car_th]

cat_cols = cat_cols + num_but_cat
cat_cols = [col for col in cat_cols if col not in cat_but_car]

num_cols = [col for col in dataframe.columns if dataframe[col].dtypes in [«int64», «float64»] and «ID» not in col.upper()]
num_cols= [col for col in num_cols if col not in num_but_cat]

date_cols = [col for col in dataframe.columns if ‘Date’ in col]

for col in date_cols:
dataframe[col] = pd.to_datetime(dataframe[col], errors = ‘coerce’)

return cat_cols, num_cols, cat_but_car,date_cols

grab_col_name(df)

## Output :

([‘Ship Mode’, ‘Segment’, ‘State’, ‘Region’, ‘Category’, ‘Sub-Category’],
[‘Sales’],
[‘Row ID’,
‘Order ID’,
‘Customer Name’,
‘City’,
‘Postal Code’,
‘Product ID’,
‘Product Name’],
[‘Order Date’, ‘Ship Date’])

##

cat_cols, num_cols, cat_but_car, date_cols = grab_col_name(df)

В настоящее время я классифицирую переменные в наборе данных на основе их типов. Далее я определю функцию для визуализации частоты классов категорий в категориальных данных.# Function for Categorical Variables:

def cat_summary(dataframe, plot=False, categorical_columns = cat_cols, threshold=10):
for col_name in categorical_columns:
print(«#» * 30 + » Unique Values Of » + col_name + » column » + «#» * 30)
print(col_name + «: » + str(dataframe[col_name].nunique()))

print(«#» * 30 + » Frequencies of Class Values » + «#» * 30)
value_counts = dataframe[col_name].value_counts()
if dataframe[col_name].nunique() > threshold:
top_n_values = value_counts.head(threshold)
#other_count = value_counts[threshold:].sum()
#top_n_values[‘Others’] = other_count
value_counts = top_n_values

print(pd.DataFrame({col_name: value_counts,
«Ratio»: value_counts / len(dataframe[col_name])}))

if plot:
rgb_values = sns.color_palette(«Set2», 6)
sns.set_theme(style=»darkgrid»)

if dataframe[col_name].nunique() > threshold:
plt.figure(figsize=(13,6))

ax = sns.barplot(x=value_counts.index, y=value_counts.values, palette=rgb_values)
ax.set_xticklabels(ax.get_xticklabels(), rotation=60)
for p in ax.patches:
ax.annotate(f’\n{p.get_height()}’, (p.get_x() + p.get_width() / 2., p.get_height()), ha=’center’, va=’top’, color=’white’, size=10)
plt.show()

# Function for Numerical Variables :

def num_summary(dataframe, plot = False):
quantiles = [0, 0.25, 0.50, 0.75, 1.0]

for col_name in num_cols:
print(«#»*30 + » Distribution of » + col_name+ » column» + «#»*30)
print(dataframe[col_name].describe(quantiles).T)

if plot:
sns.histplot(data = dataframe, x = col_name)
plt.xlabel(col_name)
plt.title(«The Distribution of «+ col_name)
plt.grid(True, alpha = 0.5)
plt.show(block = True)

# Function For Percent Format

def autopct_format(values):
def my_format(pct):
total = sum(values)
val = int(round(pct*total/100.0))
return ‘ ${v:d}’.format(v=val)
return my_formatcategorical_columns = [‘Ship Mode’, ‘Segment’, ‘State’,’Region’,
‘Category’, ‘Sub-Category’, ‘Customer Name’, ‘City’]

cat_summary(df,True,categorical_columns = categorical_columns,threshold = 10)############################## Unique Values Of Ship Mode column ##############################
Ship Mode: 4
############################## Frequencies of Class Values ##############################
Ship Mode Ratio
Standard Class 5859 0.60
Second Class 1902 0.19
First Class 1501 0.15
Same Day 538 0.05

Когда мы изучаем переменную Ship Mode, мы видим, что класс с именем Standard Class преимущественно находится в наборе данных. Заказы, как правило, отправляются стандартно.############################## Unique Values Of Segment column ##############################
Segment: 3
############################## Frequencies of Class Values ##############################
Segment Ratio
Consumer 5101 0.52
Corporate 2953 0.30
Home Office 1746 0.18

В переменной segment мы видим, что доминирует группа клиентов типа Потребитель. На втором месте стоит название Corporate.############################## Unique Values Of State column ##############################
State: 49
############################## Frequencies of Class Values ##############################
State Ratio
California 1946 0.20
New York 1097 0.11
Texas 973 0.10
Pennsylvania 582 0.06
Washington 504 0.05
Illinois 483 0.05
Ohio 454 0.05
Florida 373 0.04
Michigan 253 0.03
North Carolina 247 0.03

Когда мы смотрим на частоты штатов в наборе данных, мы видим, что Калифорния доминирует.############################## Unique Values Of Region column ##############################
Region: 4
############################## Frequencies of Class Values ##############################
Region Ratio
West 3140 0.32
East 2785 0.28
Central 2277 0.23
South 1598 0.16

Когда мы смотрим на переменную region, распределение кажется несколько однородным, а частоты классов кажутся близкими друг к другу. Но Запад виден более интенсивно.############################## Unique Values Of Category column ##############################
Category: 3
############################## Frequencies of Class Values ##############################
Category Ratio
Office Supplies 5909 0.60
Furniture 2078 0.21
Technology 1813 0.18

Мы видим, что в переменной Category, в которой хранятся категории продаваемых товаров, есть 3 класса. Я вижу, что класс канцелярских принадлежностей является доминирующим. Мебель я вижу на втором месте.############################## Unique Values Of Sub-Category column ##############################
Sub-Category: 17
############################## Frequencies of Class Values ##############################
Sub-Category Ratio
Binders 1492 0.15
Paper 1338 0.14
Furnishings 931 0.10
Phones 876 0.09
Storage 832 0.08
Art 785 0.08
Accessories 756 0.08
Chairs 607 0.06
Appliances 459 0.05
Labels 357 0.04

Когда я смотрю на распределение частот классов переменной Sub-categories, мы видим, что такие классы, как Binders и Papers, являются наиболее продаваемыми подкатегориями.############################## Unique Values Of City column ##############################
City: 529
############################## Frequencies of Class Values ##############################
City Ratio
New York City 891 0.09
Los Angeles 728 0.07
Philadelphia 532 0.05
San Francisco 500 0.05
Seattle 426 0.04
Houston 374 0.04
Chicago 308 0.03
Columbus 221 0.02
San Diego 170 0.02
Springfield 161 0.02

Когда мы смотрим на переменную City, мы видим, что класс с именем New York City является доминирующим.

Я использую как гистограммы, так и ящичковые диаграммы для визуализации распределения стоимости продаж. Когда мы исследуем статистические характеристики значений продаж с помощью функции check_df, мы обнаруживаем, что ряд имеет среднее значение 231, стандартное отклонение 627 и медианное значение 54. Кроме того, максимальное наблюдаемое значение составляет 22 638. Исходя из этих данных, очевидно, что наш набор данных демонстрирует неравномерное распределение. Учитывая, что медиана значений в наборе данных меньше среднего, мы ожидаем гистограмму с наклоном вправо. В таких случаях более уместно использовать медиану в качестве меры центральной тенденции для представления нашего набора данных.import plotly.graph_objects as go
import plotly.subplots as sp

fig_histogram = go.Figure()

fig_histogram.add_trace(go.Histogram(x=df[‘Sales’], name=’Sales’, histnorm=’probability’, opacity=0.7))

fig_histogram.update_layout(
title=’Sales Distribution (Histogram)’,
xaxis_title=’Sales’,
yaxis_title=’Probability’,
barmode=’overlay’,
bargap=0.1,
legend=dict(x=0.7, y=0.95),
autosize=False,
width=500,
height=400
)

fig_boxplot = go.Figure()

fig_boxplot.add_trace(go.Box(y=df[‘Sales’], name=’Sales’, boxpoints=’all’, jitter=0.3, pointpos=-1.8))

fig_boxplot.update_layout(
title=’Sales Distribution (Box Plot)’,
yaxis_title=’Sales’,
autosize=False,
width=500,
height=400
)

mean = df[‘Sales’].mean()
median = df[‘Sales’].median()
max_value = df[‘Sales’].max()

fig_histogram.add_shape(
type=’line’,
x0=mean,
x1=mean,
y0=0,
y1=1,
line=dict(color=’red’, width=3, dash=’dash’),
name=’Mean’
)

fig_histogram.add_shape(
type=’line’,
x0=median,
x1=median,
y0=0,
y1=1,
line=dict(color=’green’, width=3, dash=’dash’),
name=’Median’
)

fig_histogram.add_shape(
type=’line’,
x0=max_value,
x1=max_value,
y0=0,
y1=1,
line=dict(color=’blue’, width=3, dash=’dash’),
name=’Max’
)

fig_dashboard = sp.make_subplots(rows=1, cols=2, subplot_titles=(«Histogram», «Box Plot»))

fig_dashboard.add_trace(fig_histogram[‘data’][0], row=1, col=1)
fig_dashboard.add_trace(fig_boxplot[‘data’][0], row=1, col=2)

fig_dashboard.update_layout(height=600, width=1200, title_text=»Sales Distribution Dashboard»)

fig_dashboard.show()

Как мы и предсказывали выше, наш набор данных демонстрирует распределение с правым смещением. Еще одним примечательным моментом является то, что стоимость продаж находится на максимальном уровне. Это значение можно рассматривать как довольно экстремальное значение по сравнению с общим распределением набора данных и даже может считаться выбросом. Однако было бы неверно утверждать, что эта величина противоречит только общему распределению данных. Рассмотрим пример сценария; У нас есть Категория и Подкатегория, которые являются разными категориальными переменными в нашем наборе данных. Некоторые из продаж, которые происходят в категориях этих категориальных переменных, могут быть дорогостоящими продуктами, приобретенными в больших количествах. Это может привести к наблюдаемым значениям в нашем наборе данных. На этом этапе важно понять взаимосвязь между числовыми и категориальными переменными. Следующие этапы нашей работы мы будем проводить именно в этом направлении. Теперь я попробую исследовать категориальные переменные, которые взаимодействуют с целевой переменной Sales в наборе данных.from matplotlib import pyplot as plt
from matplotlib import figure
def desc_stats(dataframe):
desc = dataframe.describe(percentiles =[0, 0.25, 0.50, 0.75, 0.85, 0.95, 0.99, 1]).T
desc_df = pd.DataFrame(index= dataframe.columns,
columns= desc.columns,
data= desc)

f,ax = plt.subplots(figsize=(20,
desc_df.shape[0]*0.78))
sns.heatmap(desc_df,
annot=True,
cmap = «Wistia»,
fmt= ‘.2f’,
ax=ax,
linecolor=’white’,
linewidths = 1.3,
cbar = False,
annot_kws={«size»: 12})
plt.xticks(size = 18)
plt.yticks(size = 14,
rotation = 0)
plt.title(«Descriptive Statistics», size = 14)
plt.show()

desc_stats(df[[col for col in df.columns if df[col].dtype in [‘int64’, ‘float64’] ]])

Как мы заметили, максимальное значение переменной sales в наборе данных намного больше, чем другие квантильные значения. Это указывает на то, что объем продаж, как правило, низкий, но значительно возрастает после 85%. Например, значения до 75% показывают нормальное распределение, в то время как трехкратное увеличение наблюдается между 85% и 95%. Это может свидетельствовать о том, что высокая стоимость при перепродаже действительно встречается, хотя и редко, и что клиенты с более высокими доходами обычно совершают покупки в больших объемах. Изучив эту взаимосвязь более подробно, мы попытаемся понять факторы, влияющие на стоимость продаж.

Теперь я создам гистограмму, чтобы визуализировать 10 клиентов, приносящих наибольший доход. Эта линейчатая диаграмма будет использовать величину значений продаж для раскрашивания столбцов, представляющих каждого клиента. Это позволит нам более четко выделять клиентов с более высокой стоимостью продаж, что позволит нам легче идентифицировать наиболее важных клиентов с точки зрения дохода. Такая визуализация поможет нам определить основные источники дохода бизнеса и наиболее ценные сегменты клиентов.customer_sales = df.groupby(‘Customer Name’).agg({‘Sales’:»sum»}).reset_index().sort_values(by = «Sales», ascending = False)[:10]

fig = px.bar(customer_sales, x=’Customer Name’, y=’Sales’, hover_data=[‘Customer Name’, ‘Sales’],
labels={‘Customer Name’: ‘Customer Name’, ‘Sales’: ‘Sales’},
color = ‘Sales’)

fig.update_traces(hovertemplate=’Customer Name: %{x}<br>Sales: $%{y}<extra></extra>’)

fig.update_layout(title=’Top 10 Customer’,width=1200, height=600)

fig.update_traces(text=customer_sales[‘Sales’],
texttemplate=’$%{text:.2f}’,
textfont=dict(color=’gray’, size=14, family=’Arial’),
textangle=0,textposition=’outside’)

fig.show()

Я создаю Barplot для проверки переменной Sales относительно переменной State. В наборе данных 49 различных значений состояния, но я покажу 10 штатов с наибольшим количеством продаж.state_sales = df.groupby(‘State’).agg({‘Sales’:»sum», «Order ID»: lambda orderid: orderid.nunique()}).reset_index().sort_values(by = «Sales», ascending = False)[:10]

fig = px.bar(state_sales, x=’State’, y=’Sales’, hover_data=[‘State’, ‘Sales’],
labels={‘State’: ‘State’, ‘Sales’: ‘Sales’},
color = ‘Sales’)

fig.update_traces(hovertemplate=’State: %{x}<br>Sales: $%{y}<extra></extra>’)

fig.update_layout(title=’Total Sales by State’,width=1200, height=600)

fig.update_traces(text=state_sales[‘Sales’],
texttemplate=’$%{text:.2f}’,
textfont=dict(color=’gray’, size=14, family=’Arial’),
textangle=0,textposition=’outside’)

fig.show()

city_sales = df.groupby(‘City’).agg({‘Sales’:»sum», «Order ID»: lambda orderid: orderid.nunique()}).reset_index().sort_values(by = «Sales», ascending = False)[:10]

fig = px.bar(city_sales, x=’City’, y=’Sales’, hover_data=[‘City’, ‘Sales’],
labels={‘City’: ‘City’, ‘Sales’: ‘Sales’},
color = ‘Sales’)

fig.update_traces(hovertemplate=’City: %{x}<br>Sales: $%{y}<extra></extra>’)

fig.update_layout(title=’Total Sales by City’,width=1200, height=600)

fig.update_traces(text=city_sales[‘Sales’],
texttemplate=’$%{text:.2f}’,
textfont=dict(color=’gray’, size=14, family=’Arial’),
textangle=0,textposition=’outside’)

fig.show()

Когда я рассматриваю визуализацию, то вижу, что самые высокие продажи пришлись на Нью-Йорк — NYC. Тем не менее, города штата Калифорния Лос-Анджелес, Сан-Франциско и Сан-Диего также имеют значительные продажи. Одним из основных моментов на графике является то, что Калифорния опережает Нью-Йорк. Это показывает, что в то время как Нью-Йорк имеет наибольшее количество продаж на уровне города, Калифорния находится на вершине списка при сравнении на уровне штата. Это показывает, что общий объем продаж в Калифорнии высок, и ее различные города имеют значительную долю в общем объеме продаж.

Наряду со значениями продаж, я хотел бы проанализировать количество заказов, полученных магазином на товары, приобретенные из этих городов. Я стремлюсь определить, есть ли корреляция между продажной стоимостью и количеством заказов.

Мы изучили стоимость продаж по штатам и городам. Наконец, я создаю столбчатую диаграмму в переменной State, чтобы увидеть количество полученных заказов.city_sales = df.groupby(‘City’).agg({«Order ID»: lambda orderid: orderid.nunique()}).rename(columns = {‘Order ID’ : ‘Order Count’}).reset_index().sort_values(by = «Order Count», ascending = False)[:10]

fig = px.bar(city_sales, x=’City’, y=’Order Count’, hover_data=[‘City’, ‘Order Count’],
labels={‘City’: ‘City’, ‘Order Count’: ‘Order Count’},
color = ‘Order Count’)

fig.update_traces(hovertemplate=’City: %{x}<br>Order Count: $%{y}<extra></extra>’)

fig.update_layout(title=’Total Sales by City’,width=1200, height=600)

fig.update_traces(text=city_sales[‘Order Count’],
texttemplate=’%{text:.2f}’,
textfont=dict(color=’gray’, size=14, family=’Arial’),
textangle=0,textposition=’outside’)

fig.show()

Я создам кольцевую диаграмму, чтобы визуализировать общие суммы счетов-фактур, представляющих сегменты наших клиентов. В нашем наборе данных есть три основных сегмента клиентов: потребительские клиенты (Consumer), корпоративные клиенты (Corporate) и сегмент Home Office. Эта круговая диаграмма будет использоваться для отображения отношения суммы счета к общему доходу для каждого сегмента клиентов. Таким образом, мы сможем более четко видеть, какой клиентский сегмент генерирует основную часть операционного дохода. Выделяя разные сегменты цветами, мы сможем сделать визуальное сравнение сегментов клиентов через каждый срез на графике. Эта круговая диаграмма поможет нам понять вклад сегментов клиентов в доход и лучше сформировать наши бизнес-стратегии.segment_sales = df.groupby(‘Segment’).agg({«Sales»:»sum»}).sort_values(by = «Sales», ascending = False)
segment_sales = segment_sales[[‘Sales’]]
segment_sales .reset_index(inplace = True)
total_revenues = segment_sales[«Sales»].sum()
total_revenues = str(int(total_revenues))
total_revenues = ‘$’ + total_revenues

plt.rcParams[«figure.figsize»] = (10,6)
plt.rcParams[‘font.size’] = 11.0
plt.rcParams[‘font.weight’] = 6
colors = [‘#BC243C’,’#FE840E’,’#C62168′]
explode = (0.05,0.05,0.05)
fig1, ax1 = plt.subplots()
ax1.pie(segment_sales[‘Sales’], colors = colors, labels=segment_sales[‘Segment’], autopct= autopct_format(segment_sales[‘Sales’]),startangle=90,explode=explode)
centre_circle = plt.Circle((0,0),0.85,fc=’white’)
fig = plt.gcf()
fig.gca().add_artist(centre_circle)
ax1.axis(‘equal’)
label = ax1.annotate(‘Total Revenue \n’+str(total_revenues),color = ‘red’, xy=(0, 0), fontsize=12, ha=»center»)
plt.tight_layout()
plt.show()

В нашем наборе данных есть 3 сегмента: корпоративный, домашний офис и потребительский. Эти сегменты клиентов указывают на характер клиентов, которые совершают покупки в магазине. (Другой возможностью может быть указание типа счетов-фактур, выставленных магазином.)

Когда мы смотрим на распределение в наборе данных, мы видим, что большинство индивидуальных покупателей совершают покупки. За ними следуют потребители с долей в 688 494 доллара.

Ну интересно, в каких категориях и в каких подкатегориях проходили эти распродажи.topcat_subcat = df.groupby([‘Category’,’Sub-Category’]).agg({‘Sales’:’sum’}).sort_values(by = «Sales», ascending=False)[:10]
topcat_subcat = topcat_subcat[[«Sales»]].astype(int)
topcat_subcat = topcat_subcat.sort_values(«Category»)
topcat_subcat.reset_index(inplace=True)
topcat_subcat_1 = topcat_subcat.groupby([‘Category’]).agg({‘Sales’:’sum’})
topcat_subcat_1.reset_index(inplace=True)

plt.rcParams[«figure.figsize»] = (15,10)
fig, ax = plt.subplots()
ax.axis(‘equal’)
width = 0.1
outer_colors = [‘#FE840E’,’#009B77′,’#BC243C’]
inner_colors = [‘Orangered’,’tomato’,’coral’,»darkturquoise»,»mediumturquoise»,»paleturquoise»,»lightpink»,»pink»,»hotpink»,»deeppink»]
pie = ax.pie(topcat_subcat_1[‘Sales’], radius=1, labels=topcat_subcat_1[‘Category’],colors=outer_colors,wedgeprops=dict(edgecolor=’w’))
pie2 = ax.pie(topcat_subcat[‘Sales’], radius=1-width, labels=topcat_subcat[‘Sub-Category’],autopct= autopct_format(topcat_subcat[‘Sales’]),labeldistance=0.7,colors=inner_colors,wedgeprops=dict(edgecolor=’w’), pctdistance=0.53,rotatelabels =True)
fraction_text_list = pie2[2]
for text in fraction_text_list:
text.set_rotation(315)
centre_circle = plt.Circle((0,0),0.6,fc=’white’)
fig = plt.gcf()
fig.gca().add_artist(centre_circle)
ax1.axis(‘equal’)
plt.tight_layout()
plt.show()

Я создаю кольцевую диаграмму, чтобы посмотреть на распределение значений Sales в соответствии с переменной Ship mode в наборе данных.Top_shipping = df.groupby([«Ship Mode»]).agg({‘Sales’:’sum’}).sort_values(by = «Sales», ascending=False)
Top_shipping = Top_shipping[[«Sales»]]
Top_shipping.reset_index(inplace=True)
total_revenue_ship = Top_shipping[«Sales»].sum()
total_revenue_ship = str(int(total_revenue_ship))
total_revenue_ship = ‘$’ + total_revenue_ship

plt.rcParams[«figure.figsize»] = (13,5)
plt.rcParams[‘font.size’] = 12.0
plt.rcParams[‘font.weight’] = 6
colors = [‘#BC243C’,’#FE840E’,’#C62168′,»limegreen»]
fig1, ax1 = plt.subplots()
ax1.pie(Top_shipping[‘Sales’], colors = colors, labels=Top_shipping[‘Ship Mode’], autopct= autopct_format(Top_shipping[‘Sales’]), startangle=90)
centre_circle = plt.Circle((0,0),0.82,fc=’white’)
fig = plt.gcf()
fig.gca().add_artist(centre_circle)

ax1.axis(‘equal’)
label = ax1.annotate(‘Total Revenue \n’+str(total_revenue_ship),color = ‘red’, xy=(0, 0), fontsize=12, ha=»center»)
plt.tight_layout()
plt.show()

Когда мы изучаем приобретенные товары в соответствии с типом отгрузки, мы видим, что отгрузка была произведена максимум в категории стандартного класса. Наибольшие затраты на технологии были осуществлены в перевозках первого класса и в тот же день.

Мы рассмотрели распределение переменной Sales в соответствии с различными категориальными переменными, которые у нас есть. Теперь я хочу изучить переменную Product Name, которая содержит множество уникальных значений. Я хочу видеть товары с наибольшим количеством продаж и товары с наибольшим доходом соответственно, а затем я хочу видеть товары с наибольшими продажами с кольцевой диаграммой.df_pop_items = df[[‘Product Name’, ‘Sales’, ‘Order ID’]].groupby([‘Product Name’]).agg({‘Sales’ : ‘sum’,’Order ID’ : ‘count’}).rename(columns = {‘Order ID’ : ‘Order Count’, ‘Sales’ : ‘Total Sales’}).reset_index().sort_values(by = ‘Order Count’, ascending = False)[:10]
df_pop_items2 = df[[‘Product Name’, ‘Sales’, ‘Order ID’]].groupby([‘Product Name’]).agg({‘Sales’ : ‘sum’,’Order ID’ : ‘count’}).rename(columns = {‘Order ID’ : ‘Order Count’, ‘Sales’ : ‘Total Sales’}).reset_index().sort_values(by = ‘Total Sales’, ascending = False)[:10]

print(‘Best Selling Products (By Order Amount)’ , df_pop_items, sep = ‘\n’, end = ‘\n\n\n’)
print(‘Best Selling Products (By Sales Amount)’ , df_pop_items2, sep = ‘\n’, end = ‘\n\n’)

# Output :

Best Selling Products (By Order Amount)
Product Name Total Sales Order Count
1492 Staple envelope 1675.65 47
1498 Staples 755.47 46
537 Easy-staple paper 2414.16 44
259 Avery Non-Stick Binders 217.32 20
1495 Staple remover 263.09 18
1499 Staples in misc. colors 459.66 18
941 KI Adjustable-Height Table 4466.66 17
1510 Storex Dura Pro Binders 278.59 17
1496 Staple-based wall hangings 422.29 16
992 Logitech 910-002974 M325 Wireless Mouse for We… 1409.53 15

Best Selling Products (By Sales Amount)
Product Name Total Sales Order Count
404 Canon imageCLASS 2200 Advanced Copier 61599.82 5
649 Fellowes PB500 Electric Punch Plastic Comb Bin… 27453.38 10
444 Cisco TelePresence System EX90 Videoconferenci… 22638.48 1
785 HON 5400 Series Task Chairs for Big and Tall 21870.58 8
685 GBC DocuBind TL300 Electric Binding System 19823.48 11
687 GBC Ibimaster 500 Manual ProClick Binding System 19024.50 9
804 Hewlett Packard LaserJet 3310 Copier 18839.69 8
786 HP Designjet T520 Inkjet Large Format Printer … 18374.90 3
682 GBC DocuBind P400 Electric Binding System 17965.07 6
812 High Speed Automatic Electric Letter Opener 17030.31 3

В приведенном выше списке мы видим товары с наибольшими продажами с точки зрения дохода (Total Sales) и количества покупок (Item Count), предоставленных магазину. Наибольшая сумма продаж рассматривается как «Конверт со скобами» в зависимости от количества единиц.

Когда мы сортируем список таким образом, мы видим большее влияние на стоимость продаж некоторых продуктов. Например, цена, уплаченная за продукт под названием «Cisco TelePresence System EX90 Videoconferenci…» составляет $22638, даже если будет совершена одна покупка. Это довольно дорогая вещь.Top_products = df.groupby([«Product Name»]).agg({‘Sales’:’sum’}).sort_values(«Sales»,ascending=False).head(8)
Top_products = Top_products[[«Sales»]].round(2)
Top_products.reset_index(inplace=True)
total_revenue_products = ‘$’ +str(int(Top_products[«Sales»].sum()))

plt.rcParams[«figure.figsize»] = (13,7)
plt.rcParams[‘font.size’] = 12.0
colors = [‘#ff9999′,’#66b3ff’,’#99ff99′,’#ffcc99′,’#55B4B0′,’#E15D44′,’#009B77′,’#B565A7′] # colors are defined for the pie chart
explode = (0.05,0.05,0.05,0.05,0.05,0.05,0.05,0.05)
fig1, ax1 = plt.subplots()
ax1.pie(Top_products[‘Sales’], colors = colors, labels=Top_products[‘Product Name’], autopct= autopct_format(Top_products[‘Sales’]), startangle=90,explode=explode)
centre_circle = plt.Circle((0,0),0.80,fc=’white’)
fig = plt.gcf()
fig.gca().add_artist(centre_circle)

ax1.axis(‘equal’)
label = ax1.annotate(‘Total Revenue \n of these products \n’+str(total_revenue_products),color = ‘red’, xy=(0, 0), fontsize=12, ha=»center»)
plt.tight_layout()
plt.show()

Мы завершили исследование по поисковому анализу данных, основанному на категориальных переменных в наборе данных. Если вы помните, у нас было две переменные типа даты в наборе данных: Дата отгрузки и Дата заказа. На текущем этапе я хочу проанализировать зависящие от времени изменения нашей переменной Sales в соответствии с датой заказа между этими двумя переменными. Далее я снова воспользуюсь линейными диаграммами, чтобы визуализировать, как переменная Sales изменяется с течением времени для различных категориальных переменных.df_time = df.copy()

df_time[‘Order Date’] = pd.to_datetime(df[‘Order Date’], format=’%d/%m/%Y’)
df_time[‘Ship Date’] = pd.to_datetime(df[‘Ship Date’], format=’%d/%m/%Y’)

df_time[‘Month-Year’] = df_time[‘Order Date’].dt.strftime(‘%b-%Y’)

monthly_sales = df_time.groupby(‘Month-Year’).agg({‘Sales’: ‘sum’}).reset_index()
monthly_sales[‘Order’] = pd.to_datetime(monthly_sales[‘Month-Year’], format=’%b-%Y’)
monthly_sales = monthly_sales.sort_values(‘Order’)

fig = px.line(monthly_sales, x=’Month-Year’, y=’Sales’, title=’Monthly Sales’,
labels={‘Month-Year’: ‘Month-Year’, ‘Sales’: ‘Sales’})

fig.update_traces(hovertemplate=’Date: %{x}<br>Sales: %{y:,.2f}<extra></extra>’)

fig.update_xaxes(tickangle=90, tickmode=’array’, tickvals=monthly_sales[‘Month-Year’],
ticktext=monthly_sales[‘Month-Year’])

fig.update_layout(width=900, height=600)

fig.show()

Когда мы анализируем изменение стоимости продаж с течением времени, изначально кажется, что в ряду есть тенденция. Однако было бы преждевременно делать выводы, основываясь исключительно на внешнем виде. Чтобы обеспечить более точную интерпретацию, мы можем обратиться к тесту Расширенного Дики-Фуллера (ADF). Я продолжу анализ временных рядов.def time_breakdown(dataframe = dataset, date_col = ‘Order Date’,categorical_columns = cat_cols, threshold = 15):
dataframe[f’Month_Year_{date_col}’] = dataframe[date_col].dt.strftime(‘%b-%Y’)
for cat_col in categorical_columns:
if dataframe[cat_col].nunique() <= threshold:
monthly_sales = dataframe.groupby([f’Month_Year_{date_col}’,cat_col]).agg({‘Sales’: ‘sum’}).reset_index()
monthly_sales[‘Order’] = pd.to_datetime(monthly_sales[f’Month_Year_{date_col}’], format=’%b-%Y’)
monthly_sales = monthly_sales.sort_values(‘Order’)

fig = px.line(monthly_sales, x=f’Month_Year_{date_col}’, y=’Sales’, title=’Monthly Sales’,
labels={f’Month-Year_{date_col}’: f’Month_Year_{date_col}’, ‘Sales’: ‘Sales’}, color = cat_col)

fig.update_traces(hovertemplate=’Date: %{x}<br>Sales: %{y:,.2f}<extra></extra>’)

fig.update_xaxes(tickangle=90, tickmode=’array’, tickvals=monthly_sales[f’Month_Year_{date_col}’],
ticktext=monthly_sales[f’Month_Year_{date_col}’])

fig.update_layout(width=1000, height=600)

fig.show()time_breakdown(df_time, ‘Order Date’, categorical_columns = categorical_columns, threshold = 10)

Я создаю функцию для анализа изменения значений продаж с течением времени на основе категориальных переменных в нашем наборе данных. Эта функция эффективно представляет разбивку значений продаж по указанным категориальным столбцам в виде списка по отношению ко времени.

Теперь я генерирую новые переменные времени, используя столбец Order Date, который у нас есть. Далее я рассмотрю колебания продаж продукции с течением времени.df_time[‘WeekDay’] = df_time[‘Order Date’].dt.day_name()

start_date = df_time[‘Order Date’].min()
df_time[‘Week’] = ((((df_time[‘Order Date’] — start_date).dt.days)/7) +1).astype(int)

Ранее мы уже видели изменение продаж в соответствии с переменной Order Date. Но я не создавал этот график день за днем. Теперь пришло время воссоздать диаграмму.df_plot = df_time.groupby(‘Order Date’).agg({‘Sales’:’sum’}).rename(columns = {‘Sales’:’Total Sales’}).reset_index().sort_values(by = ‘Order Date’, ascending = True)

fig = px.line(df_plot, x=»Order Date», y=»Total Sales», title=’Total Sales by Date’)
fig.show()

Как видно на рисунке, в продажах нет видимой закономерности в отношении количества заказов. Однако при ближайшем рассмотрении можно выявить тенденцию сезонности. Когда мы строим линейный график на основе общей стоимости продаж, мы не видим четкой закономерности сезонности. Кроме того, в данных нет четкой дневной тенденции. Мы проведем тесты, чтобы убедиться в стационарности и сезонности данных. Моя следующая цель — выяснить, есть ли различия в продажах в разные дни недели. Другими словами, я планирую сгруппировать дни недели, чтобы определить, приходится ли пик продаж на конкретные дни.day_order = [‘Monday’, ‘Tuesday’, ‘Wednesday’, ‘Thursday’, ‘Friday’, ‘Saturday’, ‘Sunday’]

df_plot_dayofweek = df_time[[‘WeekDay’, ‘Sales’, ‘Order ID’]].groupby([‘WeekDay’]).agg({‘Sales’: ‘sum’, ‘Order ID’: ‘count’}).rename(columns = {‘Sales’: ‘Total Sales’, ‘Order ID’:’Order Count’}).reset_index()

df_plot_dayofweek = df_plot_dayofweek.set_index(‘WeekDay’).loc[day_order].reset_index()

fig_dayofweek = px.bar(df_plot_dayofweek, x = ‘WeekDay’, y = ‘Order Count’,hover_data=[‘WeekDay’, ‘Order Count’],
labels={‘WeekDay’: ‘WeekDay’, ‘Order Count’: ‘Order Count’},
color = ‘Order Count’)

fig_dayofweek.update_traces(hovertemplate=’WeekDay: %{x}<br>Order Count: %{y}<extra></extra>’)

fig_dayofweek.update_layout(title=’Total Sales by Day Of Week’,width=1200, height=600)

fig_dayofweek.update_traces(text=df_plot_dayofweek[‘Order Count’],
texttemplate=’%{text:.2f}’,
textfont=dict(color=’gray’, size=14, family=’Arial’),
textangle=0,textposition=’outside’)
fig_dayofweek.show()

Как видно на графике, пик продаж приходится на определенные вторники и постепенно снижается до четверга. В выходные дни они снова начинают подниматься.day_order = [‘Monday’, ‘Tuesday’, ‘Wednesday’, ‘Thursday’, ‘Friday’, ‘Saturday’, ‘Sunday’]

df_plot_dayofweek = df_time[[‘WeekDay’, ‘Sales’, ‘Order ID’]].groupby([‘WeekDay’]).agg({‘Sales’: ‘sum’, ‘Order ID’: ‘count’}).rename(columns = {‘Sales’: ‘Total Sales’, ‘Order ID’:’Order Count’}).reset_index()

df_plot_dayofweek = df_plot_dayofweek.set_index(‘WeekDay’).loc[day_order].reset_index()

fig_dayofweek = px.bar(df_plot_dayofweek, x = ‘WeekDay’, y = ‘Total Sales’,hover_data=[‘WeekDay’, ‘Total Sales’],
labels={‘WeekDay’: ‘WeekDay’, ‘Total Sales’: ‘Total Sales’},
color = ‘Total Sales’)

fig_dayofweek.update_traces(hovertemplate=’WeekDay: %{x}<br>Total Sales: $%{y}<extra></extra>’)

fig_dayofweek.update_layout(title=’Total Sales by Day Of Week’,width=1200, height=600)

fig_dayofweek.update_traces(text=df_plot_dayofweek[‘Total Sales’],
texttemplate=’$%{text:.2f}’,
textfont=dict(color=’gray’, size=14, family=’Arial’),
textangle=0,textposition=’outside’)
fig_dayofweek.show()

Когда мы изучаем стоимость продаж и количество продуктов, продаваемых отдельно по дням, мы замечаем, что самые высокие продажи продуктов происходят во вторник, за которым следуют вторник, суббота и понедельник. Однако с точки зрения выручки днями с самыми высокими продажами продукции являются суббота, вторник, понедельник и среда. Интересно, что Sunday, который занимает первое место по продажам продукции, опускается на более низкое место по выручке.

Основываясь на этом наблюдении, мы можем сделать вывод, что товары с более высокой денежной стоимостью, как правило, заказываются по субботам. Поскольку вторник возглавляет оба чарта, это означает, что товары, продаваемые по вторникам, не только многочисленны, но и ценны. По средам, как и по субботам, объем заказов меньше, но эти заказы, как правило, более ценны

Теперь я хочу свести свой анализ временных рядов к чуть меньшей линзе и посмотреть на ход продаж за год. Я хотел бы продолжить свою работу, создав переменную с именем single_year для изучения продаж, сделанных в 2017 году случайным образом из серии.SINGLE_YEAR = 2017

df_sales_single_year = df_time[df_time[‘Order Date’].dt.year == SINGLE_YEAR]
df_plot_single_year = df_sales_single_year[[‘Order Date’,’Sales’,’Order ID’]].groupby([‘Order Date’]).agg({‘Sales’: ‘sum’, ‘Order ID’:’count’}).rename(columns = {‘Sales’:’Total Sales’, ‘Order ID’: ‘Order Count’}).reset_index()
fig = px.bar(df_plot_single_year, x=’Order Date’, y=’Order Count’,facet_col_spacing=0, title=f’Total Sales by Day {SINGLE_YEAR}’)
fig.update_layout(bargap=0.0,bargroupgap=0.0, height = 600, width = 1200)
fig.show()

df_plot = df_sales_single_year[[‘Order Date’,’Category’, ‘Sales’,’Order ID’]].groupby([‘Order Date’,’Category’]).agg({‘Sales’:’sum’, ‘Order ID’:’count’}).rename(columns = {‘Sales’: ‘Total Sales’, ‘Order ID’: ‘Order Count’}).reset_index()
fig = px.bar(df_plot, x=’Order Date’, y=’Total Sales’,color=»Category», title=f’Item Sales by Date — {SINGLE_YEAR}’)
fig.update_layout(bargap=0.0,bargroupgap=0.0, height = 600, width = 1200, showlegend = True)
fig.show()

df_plot_weekly = df_sales_single_year[[‘Week’,’Category’,’Order ID’,’Sales’]].groupby([‘Week’,’Category’]).agg({‘Sales’:’sum’, ‘Order ID’: ‘count’}).rename(columns = {‘Sales’: ‘Total Sales’, ‘Order ID’ : ‘Order Count’}).reset_index()
fig = px.bar(df_plot_weekly, x=’Week’, y=’Total Sales’,color=»Category», title=f’Total Sales by Week — {SINGLE_YEAR}’)
fig.update_layout(bargap=0.0,bargroupgap=0.0, height = 600 , width = 1200, showlegend = True)
fig.show()

day_order = [«Monday», «Tuesday», «Wednesday», «Thursday», «Friday», «Saturday», «Sunday»]

df_plot_daily = df_time[[‘WeekDay’,’Category’,’Sales’]].groupby([‘WeekDay’,’Category’]).agg({‘Sales’:’sum’}).rename(columns = {‘Sales’:’Total Sales’}).reset_index()
df_plot_daily = df_plot_daily.set_index(‘WeekDay’).loc[day_order].reset_index()
fig = px.bar(df_plot_daily, x=’WeekDay’, y=’Total Sales’,color=»Category», title=f’Total Sales by Day of Week — {SINGLE_YEAR}’)
fig.update_layout(bargap=0.0,bargroupgap=0.0, height = 600 , width = 900, showlegend = True)
fig.show()

df_pop_cats = df_time[[‘Sub-Category’, ‘Sales’, ‘Order ID’]].groupby([‘Sub-Category’]).agg({‘Sales’ : ‘sum’,»Order ID»: lambda orderid: orderid.nunique()}).rename(columns = {‘Order ID’ : ‘Order Amount’, ‘Sales’ : ‘Total Sales’}).reset_index().sort_values(by = ‘Total Sales’, ascending = False)
df_pop_cats[:10]

df_plot_sc = df_sales_single_year[[‘Order Date’,’Sub-Category’, ‘Sales’]].groupby([‘Order Date’,’Sub-Category’]).agg({‘Sales’:’sum’}).rename(columns = {‘Sales’: ‘Total Sales’}).reset_index()
fig = px.bar(df_plot_sc, x=’Order Date’, y=’Total Sales’,color=»Sub-Category», title=f’Item Sales by Date — {SINGLE_YEAR}’)
fig.update_layout(bargap=0.0,bargroupgap=0.0, height = 600, width = 900)
fig.show()

df_plot_sc_weekly = df_sales_single_year[[‘Week’,’Sub-Category’,’Sales’]].groupby([‘Week’,’Sub-Category’]).agg({‘Sales’:’sum’}).rename(columns = {‘Sales’: ‘Total Sales’}).reset_index()
fig = px.bar(df_plot_sc_weekly, x=’Week’, y=’Total Sales’,color=»Sub-Category», title=f’Total Sales by Week — {SINGLE_YEAR}’)
fig.update_layout(bargap=0.0,bargroupgap=0.0, height = 600 , width = 900)
fig.show()

day_order = [«Monday», «Tuesday», «Wednesday», «Thursday», «Friday», «Saturday», «Sunday»]

df_plot_sc_daily = df_time[[‘WeekDay’,’Sub-Category’,’Sales’]].groupby([‘WeekDay’,’Sub-Category’]).agg({‘Sales’:’sum’}).rename(columns = {‘Sales’:’Total Sales’}).reset_index()
df_plot_sc_daily = df_plot_sc_daily.set_index(‘WeekDay’).loc[day_order].reset_index()
fig = px.bar(df_plot_sc_daily, x=’WeekDay’, y=’Total Sales’,color=»Sub-Category», title=f’Total Sales by Day of Week — {SINGLE_YEAR}’)
fig.update_layout(bargap=0.0,bargroupgap=0.0, height = 600 , width = 900)
fig.show()

При анализе изменений продаж в категории «Телефоны», которая может похвастаться наибольшей стоимостью продаж среди подкатегорий, мы замечаем, что пик продаж телефонов пришелся на понедельник и следующий вторник. В четверг наблюдались самые низкие продажи телефонов из-за общего снижения продаж в магазинах, но примечательно, что подкатегория телефонов по-прежнему представляла наибольшую долю ежедневных продаж по четвергам.

Еще одна отчетливая модель продаж проявляется в подкатегории «Машины», где продажи машин, по-видимому, достигают своего пика по средам и воскресеньям.

В разделе мебели подкатегория «Стулья» демонстрирует устойчивую тенденцию продаж. При ближайшем рассмотрении становится очевидно, что больше всего стульев продается по вторникам и субботам.

Прежде чем перейти к этапам моделирования и предварительной обработки, я намерен оценить стационарность и сезонность данных временных рядов в рамках анализа временных рядов. Моя цель состоит в том, чтобы определить, существует ли какая-либо сезонность или стационарность в наших данных.

Я выполнил свои исследования по анализу временных рядов на основе переменной Order Date. Теперь я хочу проверить стационарность временного ряда.df_time[‘Month-Year’] = df_time[‘Order Date’].dt.strftime(‘%b-%Y’)

monthly_sales = df_time.groupby(‘Month-Year’).agg({‘Sales’: ‘sum’}).reset_index()
monthly_sales[‘Order’] = pd.to_datetime(monthly_sales[‘Month-Year’], format=’%b-%Y’)
monthly_sales = monthly_sales.sort_values(‘Order’)

fig = px.line(monthly_sales, x=’Month-Year’, y=’Sales’, title=’Monthly Sales’,
labels={‘Month-Year’: ‘Month-Year’, ‘Sales’: ‘Sales’})

fig.update_traces(hovertemplate=’Date: %{x}<br>Sales: $%{y:,.2f}<extra></extra>’)

fig.update_xaxes(tickangle=90, tickmode=’array’, tickvals=monthly_sales[‘Month-Year’],
ticktext=monthly_sales[‘Month-Year’])

fig.update_layout(width=900, height=600)

fig.show()

На первый взгляд, мы можем легко заметить четкую тенденцию в данных временных рядов. Когда мы группируем данные о продажах по месяцам-годам, становится очевидно, что продажи постоянно растут с течением времени.

Другие тесты, которые я проведу, включают оценку сезонности, циклических закономерностей и стагнации. В первоначальном обзоре мы видим, что данные временных рядов проявляют сезонность, достигая пиковых точек через определенные промежутки времени. Понятие цикличности в чем-то схоже с сезонностью, хотя и может не оказывать существенного влияния на процессы прогнозирования. Понимание этого может быть ценным. Он относится к серии, повторяющей определенный узор через определенные промежутки времени. Например, мы ожидаем, что продажи мороженого вырастут в летние месяцы. Если бы мы смоделировали продажи мороженого, мы бы обнаружили, что в летние месяцы они неизменно следуют одной и той же схеме.plt.figure(figsize = (15,6))
plt.style.use(‘seaborn-colorblind’)
plt.grid(True, alpha = 0.6)

sns.kdeplot(df.sort_values(by= [‘Order Date’], ascending=True).loc[(df[‘Order Date’] == ‘2015-03-24’), ‘Sales’], label = ‘2015-March’)
sns.kdeplot(df.sort_values(by= [‘Order Date’], ascending=True).loc[(df[‘Order Date’] == ‘2016-03-24’), ‘Sales’], label = ‘2016-March’)
sns.kdeplot(df.sort_values(by= [‘Order Date’], ascending=True).loc[(df[‘Order Date’] == ‘2017-03-24’), ‘Sales’], label = ‘2017-March’)

plt.xlim(left = 0, right = 5000)
plt.show()

Как вы можете видеть, когда я смотрю на распределение данных о продажах в тот же день за 3 разных года, эти данные обычно не распределены. Вот почему я хочу использовать медианное значение, а не среднее в функции передискретизации, чтобы изучить данные о продажах в течение дня.df[‘Order Date’] = pd.to_datetime(df[‘Order Date’], format=’%d/%m/%Y’)
df_sorted = df.sort_values(by= [‘Order Date’], ascending=True)

df_sorted = df_sorted.set_index(«Order Date»)

df_time = pd.DataFrame(df_sorted[‘Sales’])

df_time = pd.DataFrame(df_time[‘Sales’].resample(‘D’).median())
df_time = df_time.interpolate(method=’linear’)from statsmodels.tsa.stattools import adfuller
adf = adfuller(df_time)

print(adf)
print(‘\nADF = ‘, str(adf[0]))
print(‘\np-value = ‘, str(adf[1]))
print(‘\nCritical Values: ‘)

for key, val in adf[4].items(): #for loop to print the p-value (1%, 5% and 10%) and their respective values
print(key,’:’,val)

if adf[0] < val:
print(‘Null Hypothesis Rejected. Time Series is Stationary’)
else:
print(‘Null Hypothesis Accepted. Time Series is not Stationary’)

Чтобы проверить, является ли ряд стационарным, мы используем тест Адфуллера в библиотеке statsmodels. По результатам этого теста мы получили результат, что наш набор данных является стационарным. Это интересно тем, что когда мы нарисовали линейный график, то увидели, что датасет находится в восходящем тренде. Давайте проведем другой наш тест.import matplotlib.pyplot as plt

def is_stationary(y):

# «HO: Non-stationary»
# «H1: Stationary»

p_value = sm.tsa.stattools.adfuller(y)[1]
if p_value < 0.05:
print(F»Result: Stationary (H0: non-stationary, p-value: \
{round(p_value, 3)})»)
else:
print(F»Result: Non-Stationary (H0: non-stationary, p-value: \
{round(p_value, 3)})»)

def ts_decompose(y, model=»additive», stationary=False):
result = sm.tsa.seasonal_decompose(y, model=model)
fig, axes = plt.subplots(4, 1, sharex=True, sharey=False)
fig.set_figheight(8)
fig.set_figwidth(10)

axes[0].set_title(«Decomposition for » + model + » model»)
axes[0].plot(y, ‘k’, label=’Original ‘ + model)
axes[0].legend(loc=’upper left’)

axes[1].plot(result.trend, label=’Trend’)
axes[1].legend(loc=’upper left’)

axes[2].plot(result.seasonal, ‘g’, label=’Seasonality & Mean: ‘ + str(round(result.seasonal.mean(), 4)))
axes[2].legend(loc=’upper left’)

axes[3].plot(result.resid, ‘r’, label=’Residuals & Mean: ‘ + str(round(result.resid.mean(), 4)))
axes[3].legend(loc=’upper left’)
plt.show(block=True)

if stationary:
is_stationary(y)

Теперь мы согласны с тем, что набор данных является стационарным. Проведенные нами тесты показали, что наша интерпретация фразы «Набор данных находится в восходящем тренде» не соответствует действительности.

Я хотел бы объяснить приведенный выше график просто. На этом рисунке верхняя диаграмма предоставляет нам информацию об уровне набора данных и позволяет нам сделать вывод о том, является ли ряд стационарным или нет. На втором графике представлена информация о тренде, на третьем — о сезонности. Наконец, график «Остатки» показывает остатки. Невязки показывают уровни отклонений относительно данных. В ожидаемом случае среднее значение остатков около нуля равно 0.

Скользящая средняя, экспоненциальное сглаживание и SARIMA

В предыдущей части мы рассмотрели различные методы, используемые для анализа и очистки данных, чтобы подготовить данные для прогнозирования. После завершения анализа и уточнения данных мы успешно достигаем нашей цели — получения первозданного набора данных, свободного от каких-либо пропущенных значений, выбросов или других аномалий.

Теперь мы продолжим работу с тем же набором данных, который был подготовлен в предыдущей статье, и рассмотрим дальнейшие шаги в процессе прогнозирования.

Обучение тестового сплита

Давайте разделим данные между тестом и обучением. 85% данных мы будем использовать в качестве обучающих данных, а остальные — в качестве тестовых.

Таким образом, мы успешно разделили данные на обучающие и тестовые подмножества. Мы будем использовать эти подмножества в методах прогнозирования, приведенных ниже.

Метод прогнозирования 1: Простая скользящая средняя

Простая скользящая средняя (SMA) является широко используемым методом прогнозирования временных рядов. Это статистический метод, который вычисляет среднее значение ряда за указанное количество периодов или временных интервалов. Метод SMA предполагает, что прошлые значения являются репрезентативными для будущих моделей, и использует это предположение для прогнозирования будущих значений

В приведенном выше коде я взял окно скользящей средней в 6 месяцев. Ниже приведен график тестового, обучающего и прогнозируемого наборов данных.

Теперь давайте проверим ошибку, чтобы оценить точность.

Точность модели составляет менее 50%, что означает, что она не подходит для этого набора данных. Важно отметить, что метод SMA является простым и интуитивно понятным подходом, но имеет некоторые ограничения. Он может не фиксировать сложные закономерности или изменения в базовых данных. Кроме того, он придает одинаковый вес всем точкам данных в пределах окна, что может не подходить для всех временных рядов. Другие продвинутые методы, такие как экспоненциальное сглаживание, ARIMA или алгоритмы машинного обучения, могут быть использованы для преодоления этих ограничений и предоставления более точных прогнозов.

Метод прогнозирования 2: Простое экспоненциальное сглаживание

Простое экспоненциальное сглаживание (SES) — это базовый метод прогнозирования временных рядов, который присваивает экспоненциально убывающие веса прошлым наблюдениям.

Рассмотрим код ниже. Здесь мы импортируем SimpleExpSmoothing из пакета statsmodel для создания модели

Выше я использовал уровень сглаживания 0.001 для получения наилучшей точности. Ниже приведен график тестового, обучающего и прогнозируемого наборов данных.

Теперь давайте проверим ошибку, чтобы оценить точность.

Точность модели составляет примерно 53,2%, что определенно лучше, чем раньше, но все еще не является хорошей моделью.

Метод прогнозирования 3: Сезонная авторегрессионная интегрированная скользящая средняя (SARIMA)

SARIMA, что расшифровывается как сезонная авторегрессионная интегрированная скользящая средняя, представляет собой модель прогнозирования временных рядов, которая сочетает в себе модель ARIMA с сезонностью. Она является расширением модели ARIMA и предназначена для обработки данных временных рядов, которые демонстрируют сезонные закономерности.

Модели ARIMA эффективны для прогнозирования данных стационарных временных рядов (данных с постоянным средним и дисперсией). Однако многие реальные временные ряды демонстрируют сезонность, когда закономерности повторяются через регулярные промежутки времени, например ежедневно, ежемесячно или ежегодно. SARIMA решает эту проблему, включая сезонные различия и сезонные компоненты в модель ARIMA.

Подтверждение стационарности

Для подтверждения стационарности данных временных рядов используются два наиболее распространенных метода: расширенный критерий Дики-Фуллера (ADF) и критерий Квятковски-Филлипса-Шмидта-Шина (KPSS). Эти тесты исследуют различные аспекты стационарности и могут использоваться для дополнения друг друга при оценке стационарности набора данных.

  1. Расширенный тест Дики-Фуллера (ADF). Тест ADF проверяет наличие единичного корня, указывающего на нестационарность. Если p-значение ниже выбранного уровня значимости, ряд считается стационарным.
Тест ADF

В приведенном выше коде dataFinal — это данные, которые мы подготовили в прошлой статье. Код дает p-значение 0,016, что меньше критического значения 0,05, что означает, что в соответствии с этой серией тестов является стационарным

2. Критерий Квятковски-Филлипса-Шмидта-Шина (KPSS): Тест KPSS фокусируется на нулевой гипотезе стационарности. Если p-значение выше уровня значимости, то ряд неподвижен.

Тест KPSS

Код дает p-значение 0.1, что больше критического значения 0.05, что означает, что нулевая гипотеза отвергнута и, соответственно, ряд неподвижен.

Теперь применим метод SARIMA для получения параметров модели, которые будут использоваться для аппроксимации

Ниже приведены выходные данные, в которых показано несколько шагов, выполненных для оценки модели, и предоставление параметров для нее.

Давайте проверим сводку модели.

Теперь, когда у нас есть параметры модели, давайте подогнаем модель и сделаем прогноз для тестовых данных

Теперь прогноз готов, давайте построим то же самое, чтобы проверить, как он выглядит с наборами данных для обучения и тестирования.

Вышеописанный выглядит довольно неплохо и лучше, чем предыдущие модели. Проверим погрешность и точность модели.

Модель показывает точность примерно 98,5%, что делает эту модель лучшей точностью среди всех моделей. Таким образом, это наша окончательная модель, которая может быть использована для прогнозирования ветрогенерации.

Заключение

Наш анализ набора данных о ветроэнергетике включал в себя тестирование трех различных моделей: Простая скользящая средняя, Простое экспоненциальное сглаживание и SARIMA. Результаты показали значительные различия в точности между моделями.

Метод простой скользящей средней, использующий скользящее окно в 6 месяцев, дал точность менее 50%. Этот подход может не подходить для учета сложных закономерностей и флуктуаций, присутствующих в наборе данных по ветроэнергетике.

С другой стороны, метод простого экспоненциального сглаживания показал себя немного лучше с точностью 46%. Несмотря на то, что это улучшение по сравнению с простой скользящей средней, она все еще не позволяет точно прогнозировать значения энергии ветра.

Наконец, модель SARIMA продемонстрировала высочайшую точность, достигнув впечатляющих 98,5%. Этот метод включает в себя компоненты сезонного, авторегрессионного и скользящего среднего, что позволяет ему фиксировать закономерности и зависимости, присущие данным временных рядов.

Важно отметить, что другие модели, такие как модели LSTM на основе искусственной нейронной сети (ANN) или рекуррентной нейронной сети (RNN), могут обеспечить еще более высокую точность в зависимости от конкретного варианта использования и характеристик набора данных. Эти модели известны своей способностью фиксировать сложные взаимосвязи и временные зависимости, что делает их особенно подходящими для анализа временных рядов.

В будущих анализах было бы полезно изучить потенциал моделей ANN или LSTM на наборе данных по ветроэнергетике, поскольку они могут обеспечить дальнейшее повышение точности. Тем не менее, важно тщательно оценить конкретные требования и характеристики данных, прежде чем выбрать наиболее подходящий подход к моделированию.

Прогнозирование временных рядов с помощью AutoGluon

Введение в мультимодальную библиотеку Amazon AutoML с задачей прогнозирования временной последовательности.

АвтоГлюон.

AutoGluon — это мультимодальная библиотека python с открытым исходным кодом для AutoML, запущенная Amazon. Библиотека, построенная на основе PyTorch, использует современные модели (SOTA) для достижения наилучшей производительности моделей для различных задач машинного обучения. Оснащенный моделями SOTA Deep Learning, AutoGluon предлагает решения для таких задач, как классификация изображений, обнаружение объектов, прогнозирование текста, сегментация изображений, прогнозирование временных рядов и многое другое, как на изображении обложки. Эта страница предназначена для того, чтобы показать читателям, как использовать AutoGluon для прогнозирования временных рядов, рассматривая хорошо известный Airlines Dataset с минимальным кодом.

Установите и импортируйте AutoGluon.

Рекомендуется создать новое окружение python для экспериментов Autogluon, чтобы избежать конфликтов с другими библиотеками. Дополнительные сведения см. в [1]. После создания нового окружения установите Autogluon, выполнив команду, приведенную ниже.

pip install autogluon
import autogluon

Импорт предиктора временных рядов.

AutoGluon ожидает данные в определенном формате с несколькими индексами. Один индекс должен быть меткой времени, а другой — уникальным идентификатором (это может быть любое значение). Обычный кадр данных pandas может быть преобразован в этот формат с помощью .TimeSeriesPredictorTimeSeriesDataFrame.from_data_frame()из autogluon.timeseries import TimeSeriesPredictor, TimeSeriesDataFrame

Импорт данных

import pandas

as PD data = pd.read_excel(‘airlinesData.xlsx’)data[‘Month’] = pd.to_datetime(data[‘Month’]) #convert to datetime
id = [‘airline’] * len(data)
#create список идентификаторов для добавления в dataframe data[‘id’] = id #add новый colum ‘id’ в dataframe

data.head()( )

png

Преобразование данных в кадр данных временных рядов.

data = TimeSeriesDataFrame.from_data_frame(
data,
id_column=»id»,
timestamp_column=»Month»
)
data.tail()

png

Разделение данных на обучение и тестирование.

data.shape
(96, 3)
data.columns
Index([‘Passengers’], dtype=’object’)
»’split data into train and test»’
train = data.head(77)
test = data.tail(19)

Инициализируйте и обучите предиктор.

Выбор метрики оценки зависит от того, нужен ли вам вероятностный прогноз или точечный прогноз. Ниже приведена таблица, содержащая различные метрики, но я настоятельно рекомендую вам обратиться к [2] для лучшего понимания.

В журнале обучения можно увидеть различные модели, обученные на данных обучения, выбранном параметре и метрике оценки. Библиотека использует статистические модели и модели SOTA Deep Learning для обучения. По умолчанию модель также будет опробована моделью. Его можно отключить, установив в методе.WeightedEnsembleenable_ensemble = False.fit()predictor = TimeSeriesPredictor(target=’Passengers’,
prediction_length=19,
eval_metric=»MASE»,).fit(train)No path specified. Models will be saved in: «AutogluonModels\ag-20240202_034220»
Beginning AutoGluon training…
AutoGluon will save models to ‘AutogluonModels\ag-20240202_034220’
=================== System Info ===================
AutoGluon Version: 1.0.0
Python Version: 3.8.18
Operating System: Windows
Platform Machine: AMD64
Platform Version: 10.0.22621
CPU Count: 12
GPU Count: 0
Memory Avail: 0.59 GB / 7.33 GB (8.0%)
Disk Space Avail: 262.60 GB / 476.08 GB (55.2%)
===================================================

Fitting with arguments:
{‘enable_ensemble’: True,
‘eval_metric’: MASE,
‘hyperparameters’: ‘default’,
‘known_covariates_names’: [],
‘num_val_windows’: 1,
‘prediction_length’: 19,
‘quantile_levels’: [0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9],
‘random_seed’: 123,
‘refit_every_n_windows’: 1,
‘refit_full’: False,
‘target’: ‘Passengers’,
‘verbosity’: 2}

Inferred time series frequency: ‘MS’
Provided train_data has 77 rows, 1 time series. Median time series length is 77 (min=77, max=77).

Provided dataset contains following columns:
target: ‘Passengers’

AutoGluon will gauge predictive performance using evaluation metric: ‘MASE’
This metric’s sign has been flipped to adhere to being higher_is_better. The metric score can be multiplied by -1 to get the metric value.
===================================================

Starting training. Start time is 2024-02-02 09:12:20
Models that will be trained: [‘SeasonalNaive’, ‘CrostonSBA’, ‘NPTS’, ‘AutoETS’, ‘DynamicOptimizedTheta’, ‘AutoARIMA’, ‘RecursiveTabular’, ‘DirectTabular’, ‘DeepAR’, ‘TemporalFusionTransformer’, ‘PatchTST’]
Training timeseries model SeasonalNaive.
-0.8728 = Validation score (-MASE)
0.02 s = Training runtime
6.75 s = Validation (prediction) runtime
Training timeseries model CrostonSBA.
-1.4745 = Validation score (-MASE)
0.00 s = Training runtime
15.24 s = Validation (prediction) runtime
Training timeseries model NPTS.
-2.3546 = Validation score (-MASE)
0.02 s = Training runtime
2.85 s = Validation (prediction) runtime
Training timeseries model AutoETS.
-0.8152 = Validation score (-MASE)
0.02 s = Training runtime
30.67 s = Validation (prediction) runtime
Training timeseries model DynamicOptimizedTheta.
-1.4460 = Validation score (-MASE)
0.02 s = Training runtime
29.97 s = Validation (prediction) runtime
Training timeseries model AutoARIMA.
-0.7476 = Validation score (-MASE)
0.02 s = Training runtime
21.59 s = Validation (prediction) runtime
Training timeseries model RecursiveTabular.
-0.5560 = Validation score (-MASE)
4.89 s = Training runtime
0.38 s = Validation (prediction) runtime
Training timeseries model DirectTabular.
-3.5581 = Validation score (-MASE)
0.65 s = Training runtime
0.08 s = Validation (prediction) runtime
Training timeseries model DeepAR.
-1.1489 = Validation score (-MASE)
55.81 s = Training runtime
0.08 s = Validation (prediction) runtime
Training timeseries model TemporalFusionTransformer.
-0.8554 = Validation score (-MASE)
144.14 s = Training runtime
0.02 s = Validation (prediction) runtime
Training timeseries model PatchTST.
-0.9520 = Validation score (-MASE)
32.19 s = Training runtime
0.03 s = Validation (prediction) runtime
Fitting simple weighted ensemble.
Ensemble weights: {‘AutoARIMA’: 0.09, ‘NPTS’: 0.07, ‘RecursiveTabular’: 0.04, ‘SeasonalNaive’: 0.35, ‘TemporalFusionTransformer’: 0.44}
-0.2127 = Validation score (-MASE)
1.67 s = Training runtime
31.59 s = Validation (prediction) runtime
Training complete. Models trained: [‘SeasonalNaive’, ‘CrostonSBA’, ‘NPTS’, ‘AutoETS’, ‘DynamicOptimizedTheta’, ‘AutoARIMA’, ‘RecursiveTabular’, ‘DirectTabular’, ‘DeepAR’, ‘TemporalFusionTransformer’, ‘PatchTST’, ‘WeightedEnsemble’]
Total runtime: 347.54 s
Best model: WeightedEnsemble
Best model score: -0.2127

To check the leaderboard

predictor.leaderboard()

png

Prediction.

Without specifying a model, the predictor will consider best model for prediction which is WeightedEnsemble in this case.predictions = predictor.predict(train)
predictions Model not specified in predict, will default to the model with the best validation score: WeightedEnsemble

png

By default, AutoGluon predicts quantile levels . To predict a different set of quantiles, you can use arguments as:[0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9]quantile_levels

predictor = TimeSeriesPredictor(eval_metric=”WQL”,quantile_levels=[0.1, 0.5, 0.75, 0.9]) [2].

Plot predictions.

import matplotlib.pyplot as plt

plt.figure(figsize=(20, 3))

item_id = «airline»
y_past = train.loc[item_id][«Passengers»]
y_pred = predictions.loc[item_id]
y_test = test.loc[item_id][«Passengers»]

plt.plot(y_past, label=»Past time series values»)
plt.plot(y_pred[«mean»], label=»Mean forecast»)
plt.plot(y_test, label=»Future time series values»)

plt.fill_between(
y_pred.index, y_pred[«0.1»], y_pred[«0.9″], color=»red», alpha=0.1, label=f»10%-90% confidence interval»
)
plt.legend();

png

To use a specific model from the trained models, choose from the leaderboard or the available models can be seen with following piece of code.available_models = predictor.get_model_names()
available_models[‘SeasonalNaive’,
‘CrostonSBA’,
‘NPTS’,
‘AutoETS’,
‘DynamicOptimizedTheta’,
‘AutoARIMA’,
‘RecursiveTabular’,
‘DirectTabular’,
‘DeepAR’,
‘TemporalFusionTransformer’,
‘PatchTST’,
‘WeightedEnsemble’]

В приведенном выше списке показаны модели, доступные с текущим доменом .TimeSeriesPredictor

Тренируйтесь с выбранной моделью из Model Zoo

Предпочтительная модель может быть выбрана из Model Zoo и использована. Полный список доступных моделей см. в [3]. Выбрав DeepAR в качестве интересующей модели из модели, обучение можно выполнить, как показано ниже.predictor = TimeSeriesPredictor(target=’Passengers’, prediction_length=19).fit(
train,
hyperparameters={
«DeepAR»: {},
},
)Beginning AutoGluon training…
AutoGluon will save models to ‘AutogluonModels\ag-20240201_193851’
=================== System Info ===================
AutoGluon Version: 1.0.0
Python Version: 3.8.18
Operating System: Windows
Platform Machine: AMD64
Platform Version: 10.0.22621
CPU Count: 12
GPU Count: 0
Memory Avail: 1.51 GB / 7.33 GB (20.6%)
Disk Space Avail: 262.48 GB / 476.08 GB (55.1%)
===================================================

Fitting with arguments:
{‘enable_ensemble’: True,
‘eval_metric’: WQL,
‘hyperparameters’: {‘DeepAR’: {}},
‘known_covariates_names’: [],
‘num_val_windows’: 1,
‘prediction_length’: 19,
‘quantile_levels’: [0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9],
‘random_seed’: 123,
‘refit_every_n_windows’: 1,
‘refit_full’: False,
‘target’: ‘Passengers’,
‘verbosity’: 2}

Inferred time series frequency: ‘MS’
Provided train_data has 77 rows, 1 time series. Median time series length is 77 (min=77, max=77).

Provided dataset contains following columns:
target: ‘Passengers’

AutoGluon will gauge predictive performance using evaluation metric: ‘WQL’
This metric’s sign has been flipped to adhere to being higher_is_better. The metric score can be multiplied by -1 to get the metric value.
===================================================

Starting training. Start time is 2024-02-02 01:08:51
Models that will be trained: [‘DeepAR’]
Training timeseries model DeepAR.
-0.1003 = Validation score (-WQL)
50.19 s = Training runtime
0.11 s = Validation (prediction) runtime
Not fitting ensemble as only 1 model was trained.
Training complete. Models trained: [‘DeepAR’]
Total runtime: 50.33 s
Best model: DeepAR
Best model score: -0.1003predictions = predictor1.predict(train)
predictionsModel not specified in predict, will default to the model with the best validation score: DeepAR

Прогнозирование с помощью DeepAR:predictions = predictor.predict(train)
predictions

png

Прогнозы сюжета:import matplotlib.pyplot as plt

plt.figure(figsize=(20, 3))

item_id = «airline»
y_past = train.loc[item_id][«Passengers»]
y_pred = predictions.loc[item_id]
y_test = test.loc[item_id][«Passengers»]

plt.plot(y_past, label=»Past time series values»)
plt.plot(y_pred[«mean»], label=»Mean forecast»)
plt.plot(y_test, label=»Future time series values»)

plt.fill_between(
y_pred.index, y_pred[«0.1»], y_pred[«0.9″], color=»red», alpha=0.1, label=f»10%-90% confidence interval»
)
plt.legend();

png

Преодоление ограничений древовидных моделей при прогнозировании временных рядов

Древовидные модели, такие как LightGBM, Random Forests и XGBoost, являются мощными инструментами для решения различных задач машинного обучения. Однако они сталкиваются с уникальной проблемой, связанной с прогнозированием временных рядов. Эти модели испытывают трудности с экстраполяцией тренда из-за присущей древовидным методам дизайна. Несмотря на то, что они отлично улавливают закономерности в данных, они ошибаются, когда дело доходит до проецирования этих закономерностей в будущее. В этой статье рассматриваются три стратегии решения этой проблемы и повышения предсказательной силы древовидных моделей для анализа временных рядов.

Основная проблема экстраполяции для древовидных моделей

Древовидные модели разбивают данные на ветви и листья на основе значений входных объектов. Эти разбиения выполняются для оптимизации некоторых критериев, таких как уменьшение примеси Джини для задач классификации или минимизация среднеквадратичной ошибки для задач регрессии. Затем каждый конечный узел делает прогноз на основе среднего целевого значения (для регрессии) или класса большинства (для классификации) обучающих выборок, попадающих в этот лист.

Ограничено диапазоном обучающих данных

Ограничение возникает из-за того, что эти модели по своей сути локальны; Их прогнозы основаны на конкретных экземплярах в обучающих данных, которые заполняют листья. Если новый экземпляр поступает со значениями признаков, выходящими за пределы диапазона, наблюдаемого во время обучения, модель применит прогнозирование ближайшего конечного узла. У него нет механизма для вывода трендов или закономерностей за пределами обучающих данных, что делает невозможным истинную экстраполяцию.

Точки принятия решений (узлы) в дереве определяются на основе диапазона значений, присутствующих в обучающем наборе данных. Для возраста разделение производится в 30 лет, а для дохода — в размере 50 000 долларов. Эти разбиения хорошо работают для данных в пределах наблюдаемых диапазонов, но создают проблемы, когда модель обнаруживает точки данных, выходящие за пределы этих диапазонов.

  • Проблема экстраполяции возраста: Если модель сталкивается с человеком в возрасте 31 года и старше, она больше не может принимать тонкие решения, основанные на возрасте. Выбор решения для людей старше 30 лет полностью зависит от их дохода, не принимая во внимание любые потенциальные возрастные факторы, которые не были учтены из-за ограниченного диапазона тренировочных данных.
  • Аналогичным образом, для людей старше 30 лет прогнозы модели зависят исключительно от того, выше или ниже их доход 50 000 долларов. Если человек имеет значительно более высокий доход (скажем, намного выше 50 000 долларов), модель относится к нему так же, как и к человеку с чуть более 50 000 долларов, потенциально упуская нюансы, которые более высокий доход может подсказать о его предпочтениях.

Нет механизма распознавания трендов за пределами обучающих данных

Структура дерева не позволяет распознавать или применять тренды или закономерности, которые могут возникнуть за пределами конкретных разбиений на основе обучающих данных. Например:

  • Если тенденция значительно увеличивает вероятность того, что продукт понравится людям в возрасте 40 лет и старше, дерево не может зафиксировать или применить эту тенденцию, поскольку оно не экстраполируется за пределы узла принятия решения «Возраст > 30 лет».
  • Если уровень дохода населения значительно возрастает, так что 50 000 долларов больше не являются значимым порогом для прогнозирования предпочтения продукта, дерево не может скорректировать свою границу принятия решения, чтобы отразить это изменение.

Проблема экстраполяции тренда

Для того, чтобы найти решение проблемы в случае прогнозирования временных рядов, важно сначала понять, почему древовидные модели испытывают трудности с экстраполяцией тренда. Эти модели полагаются на статистические характеристики признаков, присутствующих в обучающих данных, для составления прогнозов. Однако, когда дело доходит до прогнозирования данных временных рядов, этот подход может привести к сбою в прогнозировании будущих значений, которые следуют тренду, выходящему за рамки обучающего набора.

В следующих разделах будут показаны три возможных решения.

Стратегия 1: Устранение тренда во временных рядах

Концепция

Один из способов прогнозирования временного ряда заключается в удалении компонента тренда из исходного ряда, его оценке с помощью отдельной модели или метода, а затем добавлении обратно в прогнозы. Этот процесс называется детрендингом.

Плюсы и минусы

Плюсы: Этот метод прост и позволяет более четко анализировать циклические и сезонные закономерности данных без затеняющего эффекта тренда.

Минусы: Эффективность детрендинга в значительной степени зависит от точности метода оценки тренда. Неправильная оценка тренда может привести к потере информации и неточностям в прогнозе.

Две возможности для детренда

Устранение тренда с помощью линейной модели

В этом примере мы используем линейную модель для оценки, а затем удаляем тренд из временного ряда. После этого мы будем использовать модель для прогнозирования тренда на будущие периоды и добавим его обратно.import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from sklearn.linear_model import LinearRegression
from statsmodels.tsa.holtwinters import ExponentialSmoothing

# Загрузить набор
данных df = pd.read_csv(‘https://raw.githubusercontent.com/jbrownlee/Datasets/master/airline-passengers.csv’, parse_dates=[‘Месяц’], index_col=’Месяц’)# Линейная модель для времени тренда
= np.arange(len(df.index))model = LinearRegression()

model.fit(time.reshape(-1, 1), df[‘Пассажиры’])

# Оценка и удалить тренд тренд
= model.predict(time.reshape(-1, 1))detrended = df[‘Пассажиры’] — trend

# Прогнозирование будущего тренда еще на 12 месяцев
future_time = np.arange(len(df.index), len(df.index) + 12)future_trend = model.predict(future_time.reshape(-1, 1))

# Plot
plt.figure(figsize=(14, 7))plt.plot(df.index, df[‘Пассажиры’], label=’Original’)

plt.plot(df.index, trend, label=’Trend’)plt.plot(df.index, detrended, label=’Detrended’, linestyle=’—‘)

plt.plot(pd.date_range(df.index[-1], periods=13, closed=’ right’), future_trend, label=’Будущий тренд’, linestyle=’-.’)
plt.legend()plt.show()

Устранение тренда с помощью оценки лёсса и прогнозирования ETS

Для лессовой оценки (локально-взвешенная регрессия) можно использовать . Впоследствии мы применяем модель ETS (Error, Trend, Seasonality) без сезонной составляющей для прогнозирования будущих тенденций.statsmodelsfrom statsmodels.nonparametric.smoothers_lowess import lowess# Loess Smoothing loess_smoothed = lowess(df[‘Passengers’], time, frac=0.05)[:, 1]# Detrending detrended_loess = df[‘Passengers

‘]

— loess_smoothed# ETS Модель для прогнозирования тренда
model_ets = ExponentialSmoothing

(loess_smoothed

, trend=»additive», seasonal=None, initialization_method=»estimated»)fitted_model = model_ets.fit()# Прогнозирование тренда
future_trend_ets = fitted_model.forecast(12)# Plot
plt.figure(figsize=(14, 7))

plt.plot(df.index, df[‘Пассажиры’], label=’Оригинал’)

plt.plot(df.index, loess_smoothed, label=’Сглаженный лёсс’)plt.plot(df.index, detrended_loess, label=’Лёсс без тренда’, linestyle=’—‘)
plt.plot(pd.date_range(df.index[-1], periods=13, closed=’right’), future_trend_ets, label=’Future Trend ETS’, linestyle=’-.’)

plt.legend()plt.show()

Аргумент в лессовой оценке определяет степень сглаживания и может быть изменен в зависимости от длины временного ряда и желаемой гибкости тренда. Аналогичным образом можно настроить параметры модели ETS для достижения наилучших результатов для конкретного временного ряда.frac

Стратегия 2: Дифференциация временных рядов

Концепция

Дифференциация — это статистический метод, который преобразует временной ряд в стационарный. Он включает в себя вычисление разницы между последовательными наблюдениями во временном ряду. Цель этого процесса состоит в том, чтобы помочь стабилизировать среднее значение рядов и уменьшить тенденции и сезонность, которые могут затруднить анализ и моделирование данных. Взяв разницу между двумя последовательными наблюдениями, мы можем удалить трендовую составляющую ряда, что упростит выявление закономерностей и прогнозирование.

Плюсы и минусы

Плюсы: Это простой подход, позволяющий сделать серию более стационарной, что упрощает работу с древовидными моделями.

Недостатки: Дифференциация может привести к потере информации, а чрезмерная дифференциация может рандомизировать ряды, что затруднит моделирование.# Применить дифференциацию
df_diff = df[‘Passengers’].diff().fillna(0) # Использование fillna(0) для первого значения

NaN# Построить график разностного ряда
plt.figure(figsize=(14, 7))
plt.plot(df.index, df_diff, label=’Differenced Series’)plt.title(‘Разностная серия пассажиров’)plt.legend()plt.legend()plt.show()

Разностные временные ряды.

# Чтобы обратить разностность и восстановить исходную серию df_restored = df_diff.cumsum() + df[‘Passengers’].iloc[0]

# Постройте восстановленную серию
рядом с оригинальной
plt.figure(figsize=(14, 7))
plt.plot(df.index, df[‘Passengers’], label=’Original Series’)
plt.plot(df.index, df_restored, label=’Восстановленная серия’, linestyle=’—‘)
plt.title(‘Оригинальные и восстановленные серии по разнице’)plt.legend()
plt.show()

Стратегия 3: Использование метода линейного дерева

Концепция

Для бустинговых моделей интеграция линейного компонента в процесс построения дерева может позволить модели фиксировать как линейные, так и нелинейные закономерности, предлагая способ включения возможностей экстраполяции тренда.

Плюсы и минусы

Плюсы: Это позволяет модели экстраполировать тренды без прямого преобразования целевой переменной, охватывая более широкий спектр закономерностей в данных.

Недостатки: Этот подход характерен для определенных моделей и может увеличить сложность модели и требования к памяти.

Стандартное дерево регрессии

Убедитесь, что установлена последняя версия lightgbm.

В этом примере стандартная версия модели LightGBM использовалась для прогнозирования количества пассажиров за последние 12 месяцев в наборе данных AirPassengers. Модель была обучена на признаках, полученных из месяца каждого наблюдения, с целью зафиксировать сезонные закономерности в данных, не полагаясь напрямую на прошлые значения ряда. Этот подход предполагает, что месяц года существенно влияет на количество пассажиров, отражая сезонные тенденции путешествий.Импорт pandas как PD
Импорт numpy как NP
Импорт LightGBM как LGB
Импорт matplotlib.pyplot как PLT
из sklearn.metrics Импорт mean_squared_error
из SKLEARN.Preprocessing import OneHotEncoder

# Загрузка набора
данных url = ‘https://raw.githubusercontent.com/jbrownlee/Datasets/master/airline-passengers.csv’df
= pd.read_csv(url, parse_dates=[‘Месяц’], index_col=’Месяц’)# Извлечение месяца в виде функции
df[‘Month_Feature’] = df.index.month

# Однократное кодирование кодировщика признака
месяца = OneHotEncoder(sparse=False)

month_encoded = encoder.fit_transform(df[[‘Month_Feature’]])
month_encoded_df = pd. DataFrame(month_encoded, index=df.index, columns=[f’month_{i+1}’ for i in range(month_encoded.shape[1])])#

Соединяем закодированные объекты с исходным кадром
данных df = pd.concat([df, month_encoded_df], axis=1)# Определяем обучающий и тестовый наборы данных (последние 12 месяцев в качестве тестового набора)

train = df[:-12]test = df[-12:]# Подготовка наборов данных для LightGBM
features = [f’month_{i+1}’ for i in range(12)]X_train, y_train = train[features], train[‘Пассажиры’]
X_test, y_test = test[features], test[‘Пассажиры’]

# Настройка и обучение модели
LightGBM params = {
‘objective’: ‘регрессия’, ‘метрика’: ‘rmse’, ‘num_leaves’: 10, ‘learning_rate’: 0.05, ‘feature_fraction’: 0.9,

‘min_data_in_leaf’: 5
}
train_data = ЛГБ. Dataset(X_train, label=y_train)gbm = lgb.train(params, train_data, num_boost_round=200)# Прогноз на следующие 12 месяцев
y_pred = gbm.predict(X_test)# Вычисляем среднеквадратичное значение для прогнозов
rmse = np.sqrt(mean_squared_error(y_test, y_pred))
print(f’RMSE: {rmse}’)

# Построение фактических и прогнозируемых значений
plt.figure(figsize=(14, 7))plt.plot(train.index, y_train, label=’Обучающие данные’)plt.plot(test.index, y_test, label=’Фактические значения’, color=’green’)

plt.plot(test.index, y_pred, label=’Прогнозируемые значения’, linestyle=’—‘, color=’red’ )
plt.title(‘Прогноз авиапассажиров со стандартным LightGBM и функцией месяца’)plt.xlabel(‘Date’)plt.ylabel(‘Количество пассажиров’)plt.legend()

plt.show()

Активация опции линейного дерева

Этот вариант модели LightGBM включает параметр, который включает в себя линейные модели на листьях деревьев. Этот метод также был применен для прогнозирования последних 12 месяцев набора данных AirPassengers, используя месяц в качестве категориального признака. Метод линейного дерева направлен на то, чтобы объединить способность модели фиксировать как линейные, так и нелинейные закономерности, потенциально предлагая тонкое понимание сезонных эффектов в данных.linear_tree# Настройка и обучение модели LightGBM с помощью linear_tree
params = {
‘objective’: ‘регрессия’, ‘metric’: ‘rmse’, ‘linear_tree’: True, # Включить линейное дерево
‘num_leaves’: 31, ‘learning_rate’: 0.05, ‘feature_fraction’: 0.9,

}
train_data = ЛГБ. Dataset(X_train, label=y_train)gbm_linear = lgb.train(params, train_data, num_boost_round=200)# Прогноз на следующие 12 месяцев с использованием модели
линейного дерева y_pred_linear = gbm_linear.predict(X_test)# Вычислить RMSE для прогнозов
rmse_linear = np.sqrt(mean_squared_error(y_test, y_pred_linear)

)
print(f’RMSE (Linear Дерево): {rmse_linear}’)# Построение графика фактических и прогнозируемых значений для линейной древовидной модели
plt.figure(figsize=(14, 7))plt.plot(train.index, y_train, label=’Обучающие данные’)plt.plot(test.index, y_test, label=’Фактические значения’, color=’green’)

plt.plot(test.index, y_pred_linear, label=’Прогнозируемые значения ( Линейное дерево)’, linestyle=’—‘, color=’red’)plt.title(‘Прогноз авиапассажиров с линейным деревом LightGBM и функцией месяца’)plt.xlabel(‘Date’)plt.ylabel(‘Количество пассажиров’)plt.legend()plt.show()

Оба примера иллюстрируют применение моделей LightGBM для прогнозирования временных рядов, отличающихся в первую очередь использованием линейных деревьев. Ключевой вывод заключается в понимании того, как метод линейного дерева может предложить преимущества в обработке линейных трендов в данных при использовании lightgbm.

Заключение

Решение проблемы экстраполяции трендов в древовидных моделях требует вдумчивого рассмотрения данных и целей прогнозирования. Каждая из описанных выше стратегий предлагает способ повышения прогностической способности этих моделей, со своими преимуществами и ограничениями. Тщательно выбирая и применяя соответствующий метод, аналитики и специалисты по обработке и анализу данных могут улучшить свои прогнозы временных рядов, максимально используя имеющиеся в их распоряжении мощные инструменты моделирования.

Раскрытие возможностей VAR-моделирования для динамического прогнозирования временных рядов

1. Введение в VAR-моделирование

Векторная авторегрессия (VAR) — это ключевая статистическая модель в эконометрике и финансах, предназначенная для систем прогнозирования, в которых несколько переменных временных рядов взаимозависимы. Разработанные Кристофером Симсом в 1980 году, модели VAR отличаются от моделей с одним уравнением тем, что фиксируют динамическое взаимодействие между несколькими переменными.

Суть VAR заключается в том, что она позволяет данным раскрывать свою собственную историю, что делает ее особенно ценной для эмпирического анализа в экономике и финансах. Применяется для:

  • Прогнозирование экономических показателей: например, темпов роста ВВП и инфляции.
  • Анализ денежно-кредитной политики: понимание действий центрального банка.
  • : Изучение взаимодействия между ценами на акции и процентными ставками.

Для тех, кто хочет углубиться:

  • Книга Джеймса Д. Гамильтона «Анализ временных рядов» представляет собой обширное введение в систему VAR.
  • В статье Хельмута Люткепола «Новое введение в анализ множественных временных рядов» рассказывается о теории и применении этих моделей.
  • Раздел «Векторные модели авторегрессии и векторной коррекции ошибок» в книге Уильяма Х. Грина «Эконометрический анализ» предлагает доступный обзор.

В этой статье мы рассмотрим, как VAR можно применять в практических сценариях, таких как прогнозирование цен на электроэнергию с использованием связанных данных временных рядов.

2. Понимание данных

Успешное моделирование VAR начинается с глубокого понимания имеющихся данных. Этот шаг имеет решающее значение, поскольку качество и характеристики данных напрямую влияют на точность и надежность модели.

  1. Многомерная природа: Модели VAR основаны на взаимозависимости нескольких временных рядов. Например, при прогнозировании цен на электроэнергию учитывайте, как они могут быть связаны с погодными условиями и ценами на газ. Каждая переменная должна потенциально влиять на другие переменные или подвергаться влиянию других.
  2. Требование к стационарности: Фундаментальным аспектом моделирования VAR является стационарность. Проще говоря, статистические свойства ряда (например, среднее значение и дисперсия) должны быть постоянными с течением времени. Нестационарные данные могут приводить к недостоверным результатам, что делает необходимым тестирование и, при необходимости, преобразование данных для достижения стационарности.
  3. Частота и согласованность данных: Обеспечьте согласованное измерение данных временных рядов с течением времени. Независимо от того, идет ли речь о ежедневных, ежемесячных или ежеквартальных данных, единообразие частоты является ключом к поддержанию целостности вашей модели VAR.
  4. Длина временных рядов: Объем имеющихся у вас данных имеет значение. Более длительные временные ряды могут предоставить больше информации, но также создают такие проблемы, как повышенная сложность и возможность переобучения. Очень важно найти баланс.
  5. Качество данных: оценка данных на наличие отсутствующих значений, выбросов или ошибок. Высококачественные и чистые данные составляют основу любой надежной эконометрической модели.
  6. Исторический контекст: Понимание исторического контекста ваших данных. Экономические и финансовые временные ряды часто зависят от разовых событий, изменений в политике или экономических циклов. Их распознавание может помочь в более точной интерпретации выходных данных модели.
  7. Предварительный анализ: Прежде чем приступить к моделированию с помощью VAR, выполните исследовательский анализ данных. Визуализируйте свои серии, ищите тенденции, сезонные закономерности и потенциальные структурные разрывы. Этот шаг дает ценную информацию и направляет последующий процесс моделирования.

Таким образом, ваши данные — это не просто набор чисел; Это история, которая ждет, чтобы ее рассказали. Чем лучше вы понимаете эту историю, тем эффективнее вы можете использовать моделирование VAR для выявления сложных взаимосвязей внутри нее.

3. Подготовка данных

  • Подготовка данных является критически важным этапом в моделировании VAR. Для нашего проекта по прогнозированию цен на электроэнергию с использованием цен на газ и исторических данных о погоде мы провели детальный процесс очистки и трансформации.
  1. Импорт и упрощение данных: Мы начали с импорта наборов данных о ценах на электроэнергию, газ и погоду за прошлые периоды. Каждый набор данных был упрощен и теперь включает только релевантные переменные.
  2. Согласованность во временных рядах: мы преобразовали индексы всех наборов данных в объекты datetime для единообразия и пересчитывали их до ежедневной частоты. Это гарантирует, что наши данные временных рядов имеют согласованный ежедневный интервал, что имеет решающее значение для точного анализа VAR.
  3. Объединение источников данных: Затем наборы данных были объединены в один кадр данных. Этот процесс включал в себя:
  • Усреднение цен на электроэнергию и газ до суточных значений.
  • Выбор переменных погоды с коэффициентом вариации больше 0,5, гарантируя, что мы сосредоточимся на переменных со значительной изменчивостью.
  • Отбрасывание высоко коррелированных погодных переменных (корреляция > 0,8), чтобы избежать мультиколлинеарности.

4. Окончательный набор данных: Окончательный набор данных включает ежедневные значения цен на электроэнергию и газ, а также ключевые погодные индикаторы, такие как температура, дождь, снегопад и различные показатели облачности.

Сниппет обработанных данных выглядит следующим образом:electricity_price_per_mwh gas_price_per_mwh temperature rain snowfall cloudcover_total cloudcover_mid cloudcover_high shortwave_radiation
2021-09-01 110.47 45.78 13.35 0.0 0.0 24.5 12.5 0.0 60.5
2021-09-02 112.63 45.96 12.55 0.0 0.0 47.0 42.5 0.0 29.0

2023-05-28 44.45 31.10 13.70 0.0 0.0 17.0 12.0 1.0 201.0
2023-05-29 35.76 32.57 10.85 0.0 0.0 46.5 39.0 26.5 225.0
2023-05-30 55.55 31.50 8.70 0.0 0.0 35.0 56.0 0.0 24.0

Визуализация временных рядов

4. Проверка стационарности и обработка

4.1 Испытание на стационарность

В этом разделе мы углубимся в критический процесс тестирования на стационарность в данных временных рядов. Стационарность является основополагающим допущением при моделировании VAR. Это означает, что статистические свойства рядов, такие как среднее значение и дисперсия, остаются постоянными во времени. Обнаружение и исправление нестационарности имеет решающее значение, поскольку это может привести к ненадежным прогнозам и интерпретации моделей.

Для проверки стационарности мы использовали расширенный тест Дики-Фуллера (ADF). Этот тест помогает нам определить, содержит ли каждый временной ряд в нашем наборе данных (цены на электроэнергию, цены на газ, температуру и т. д.) единичный корень, указывающий на нестационарность.

Результаты испытаний ADF показывают, что в то время как некоторые ряды, такие как цены на электроэнергию и определенные погодные переменные, являются постоянными, другие, такие как цены на газ и температура, не являются постоянными. Нестационарные ряды демонстрируют тренды или другие характеристики, которые изменяются с течением времени, что делает их непригодными для непосредственного использования в моделировании VAR. Для этих рядов нам нужно будет применить такие методы, как дифференциация или декомпозиция, чтобы удалить тренды и стабилизировать их статистические свойства. После преобразования в стационарную форму эти серии могут быть повторно протестированы и, возможно, включены в модель VAR.

4.2 Создание стационарных временных рядов

В этом разделе мы рассмотрим задачу преобразования данных нестационарных временных рядов в форму, более подходящую для моделирования VAR. Мы фокусируемся на том, чтобы сделать серию статей о ценах на электроэнергию, газ, температуру и коротковолновое излучение стационарной. Для этого мы используем метод линейного тренда Холта, который является мощным инструментом для сглаживания данных временных рядов.import pandas as pd
import numpy as np
from statsmodels.tsa.holtwinters import Holt
import matplotlib.pyplot as plt

electricity_price = total_data[‘electricity_price_per_mwh’]
gas_price_series = total_data[‘gas_price_per_mwh’]
temperature_series = total_data[‘temperature’]
radiation_series = total_data[‘shortwave_radiation’]

electricity_price_model = Holt(total_data[‘electricity_price_per_mwh’])
electricity_price_results = electricity_price_model.fit(smoothing_level=0.1)
electricity_price_trend = electricity_price_results.fittedvalues
electricity_price_level = electricity_price_results.level
electricity_price_residual = total_data[‘electricity_price_per_mwh’] — electricity_price_trend

gas_price_model = Holt(total_data[‘gas_price_per_mwh’])
gas_price_results = gas_price_model.fit(smoothing_level=0.1)
gas_price_trend = gas_price_results.fittedvalues
gas_price_level = gas_price_results.level
gas_price_residual = total_data[‘gas_price_per_mwh’] — gas_price_trend

temperature_model = Holt(total_data[‘temperature’])
temperature_results = temperature_model.fit(smoothing_level=0.1)
temperature_trend = temperature_results.fittedvalues
temperature_level = temperature_results.level
temperature_residual = total_data[‘temperature’] — temperature_trend

radiation_model = Holt(total_data[‘shortwave_radiation’])
radiation_results = radiation_model.fit(smoothing_level=0.1)
radiation_trend = radiation_results.fittedvalues
radiation_level = radiation_results.level
radiation_residual = total_data[‘shortwave_radiation’] — radiation_trend

Мы выбираем smoothing_level 0,1, чтобы обеспечить значительное сглаживание уровневой составляющей нашего ряда. Этот выбор является осознанным, так как позволяет нам сосредоточиться на долгосрочных трендах, а не на краткосрочных колебаниях. Таким образом, мы смягчаем влияние краткосрочной волатильности и подчеркиваем основной тренд в каждом временном ряду.

Кроме того, мы рассчитываем остаток каждого ряда как разницу между фактическими значениями и соответствующими им оценочными трендами. Причина этого двоякая: во-первых, остатки часто имеют стационарные характеристики, что делает их более подходящими для моделирования VAR. Во-вторых, на будущих этапах мы намерены прогнозировать как трендовую, так и остаточную составляющие по отдельности, а затем объединять эти прогнозы для построения прогноза для основного временного ряда. Такой подход позволяет фиксировать как основной тренд, так и случайные флуктуации, присущие данным временных рядов.

На рисунке, представленном ниже, мы демонстрируем визуальное представление наших временных рядов вместе с их разложенными компонентами.

5. Определение порядка VAR

В моделировании VAR (векторная авторегрессия) определение правильного порядка (количества лагов) имеет решающее значение для точности и эффективности модели. Порядок модели VAR указывает количество предыдущих временных точек (лагов), используемых для прогнозирования текущего значения в ряду. Правильно подобранный порядок запаздывания отражает существенную динамику системы без переобучения данных.

Для определения оптимального порядка запаздывания мы часто полагаемся на информационные критерии, такие как информационный критерий Акаике (AIC), байесовский информационный критерий (BIC) или информационный критерий Ханнана-Куинна (HQIC). Эти критерии уравновешивают сложность модели (количество лагов) с хорошим соответствием, при этом более низкие значения указывают на лучшую модель.import pandas as pd
from statsmodels.tsa.api import VAR

data_for_var = pd.concat([total_data[[‘rain’,’snowfall’, ‘cloudcover_total’, ‘cloudcover_mid’, ‘cloudcover_high’]],
electricity_price_residual, gas_price_residual, temperature_residual, radiation_residual], axis=1, join=’inner’)

# Fitting a VAR model
model = VAR(data_for_var)

# Determining the optimal lag order
lag_order_results = model.select_order(maxlags=10)
print(lag_order_results.summary())VAR Order Selection (* highlights the minimums)
==================================================
AIC BIC FPE HQIC
—————————————————
0 30.26 30.32 1.387e+13 30.29
1 26.99 27.63*
5.258e+11 27.24*
2 26.94*
28.16 5.037e+11* 27.42
3 27.03 28.82 5.502e+11 27.73
4 27.12 29.48 5.997e+11 28.03
5 27.18 30.11 6.381e+11 28.32
6 27.24 30.75 6.797e+11 28.60
7 27.32 31.40 7.409e+11 28.91
8 27.38 32.04 7.869e+11 29.19
9 27.50 32.73 8.898e+11 29.53
10 27.59 33.39 9.763e+11 29.84
—————————————————

6. Оценка модели VAR

После определения оптимального порядка запаздывания для нашей модели VAR следующим шагом является оценка самой модели VAR. Это включает в себя подгонку модели к нашим данным временных рядов с выбранным порядком запаздывания. Оценка VAR-модели позволяет получить коэффициенты, описывающие взаимосвязи между каждой парой временных рядов в системе. Эти коэффициенты необходимы для понимания того, как изменения одной переменной влияют на другие, и для прогнозирования будущих значений временных рядов.# Assuming ‘data_for_var’ contains the relevant time series data
model = VAR(data_for_var)

# Fitting the model with the chosen lag order
model_fitted = model.fit(2, ic=’aic’, trend=’ct’)

# model_fitted.summary()Results for equation electricity_price_residual
================================================================================================
coefficient std. error t-stat prob
————————————————————————————————
const 5.566887 6.297416 0.884 0.377
trend 0.000486 0.011490 0.042 0.966
L1.rain 14.412011 40.236325 0.358 0.720
L1.snowfall 168.651593 73.688386 2.289 0.022
L1.cloudcover_total -0.110353 0.092959 -1.187 0.235
L1.cloudcover_
mid 0.063661 0.087542 0.727 0.467
L1.cloudcover_high -0.014581 0.063907 -0.228 0.820
L1.electricity_
price_residual 0.591264 0.041397 14.283 0.000
L1.gas_
price_residual 0.815831 0.301679 2.704 0.007
L1.temperature_
residual -2.201449 0.977021 -2.253 0.024
L1.radiation_residual 0.236744 0.076082 3.112 0.002
L2.rain -64.005815 40.195755 -1.592 0.111
L2.snowfall -8.597689 74.749992 -0.115 0.908
L2.cloudcover_
total -0.166000 0.092752 -1.790 0.073
L2.cloudcover_mid 0.092476 0.087292 1.059 0.289
L2.cloudcover_
high 0.091182 0.063286 1.441 0.150
L2.electricity_price_residual -0.121211 0.041038 -2.954 0.003
L2.gas_price_residual -0.178676 0.302302 -0.591 0.554
L2.temperature_residual 0.508086 0.949751 0.535 0.593
L2.radiation_
residual -0.204459 0.076367 -2.677 0.007
================================================================================================

7. Диагностическая проверка

На диагностическом этапе нашей модели VAR мы использовали статистику Дурбина-Уотсона для оценки автокорреляции в остатках. Автокорреляция, особенно в данных временных рядов, может быть признаком неправильной спецификации модели и может опровергнуть некоторые предположения, лежащие в основе стандартных статистических тестов.

Статистика Дурбина-Уотсона находится в диапазоне от 0 до 4, где:

  • Значение около 2 предполагает отсутствие автокорреляции.
  • Значения, значительно меньшие 2, указывают на положительную автокорреляцию.
  • Значения, значительно превышающие 2, предполагают отрицательную автокорреляцию.

Наши результаты:from statsmodels.stats.stattools import durbin_watson

# Compute Durbin-Watson statistics
dw_stats = durbin_watson(model_fitted.resid)

for col, stat in zip(data_for_var.columns, dw_stats):
print(f'{col}: {stat}’)Rain: 2.013
Snowfall: 2.004
Cloudcover Total: 2.033
Cloudcover Mid: 2.016
Cloudcover High: 2.000
Electricity Price Residual: 2.002
Gas Price Residual: 2.007
Temperature Residual: 1.975
Radiation Residual: 2.007

Все переменные показывают статистику Дурбина-Уотсона, очень близкую к 2, что указывает на то, что автокорреляция в невязках минимальна или отсутствует вовсе. Это позитивный знак, свидетельствующий о том, что наша модель VAR хорошо отражает динамику системы, не оставляя необъяснимых закономерностей в остатках.

8. Прогнозирование с помощью VAR

Оценив модель и убедившись, что она прошла диагностические проверки, мы можем использовать ее для прогнозирования будущих значений включенных переменных.import matplotlib.pyplot as plt

# Forecasting
forecast_steps = 10 # Number of steps to forecast ahead
forecast = model_fitted.forecast(model_fitted.endog, steps=forecast_steps)

forecast_df = pd.DataFrame(forecast, index=pd.date_range(start=data_for_var.index[-1] + pd.Timedelta(days=1), periods=forecast_steps, freq=’D’), columns=data_for_var.columns)

# Concatenate the predicted and last 50 observations of historical data for ‘electricity_price_residual’
concatenated_data = pd.concat([data_for_var[‘electricity_price_residual’].tail(50), forecast_df[‘electricity_price_residual’]])

# Plot the concatenated data
plt.figure(figsize=(12, 6))
plt.plot(data_for_var.index[-50:], data_for_var[‘electricity_price_residual’].tail(50), label=’Historical’, color=’blue’,marker = ‘o’)
plt.plot(forecast_df.index, forecast_df[‘electricity_price_residual’], label=’Forecast’, color=’orange’,marker = ‘*’)
plt.title(‘Historical and Predicted Electricity Price Residual’)
plt.xlabel(‘Date’)
plt.ylabel(‘Values’)
plt.legend()
plt.show()

9. Комбинирование компонента для прогнозирования первичных временных рядов

В нашем анализе мы разложили временные ряды цен на электроэнергию на трендовую и остаточную составляющие. Мы уже спрогнозировали остаточную составляющую с помощью модели VAR. Для построения комплексного прогноза цены на электроэнергию обратимся к прогнозированию трендовой составляющей. Для тренда мы будем использовать более простой подход: полиномиальную аппроксимацию. Этот метод эффективен для фиксации основного тренда во временном ряду, особенно на коротком горизонте. Комбинируя прогнозы на основе этих двух различных компонентов — тренда и остаточного — мы можем реконструировать более детальный и потенциально более точный прогноз исходного ценового ряда электроэнергии.electricity_price = total_data[‘electricity_price_per_mwh’]

# the fitted Holt model
electricity_price_model = Holt(total_data[‘electricity_price_per_mwh’])
electricity_price_results = electricity_price_model.fit(smoothing_level=0.1)
electricity_price_trend = electricity_price_results.fittedvalues
electricity_price_level = electricity_price_results.level
electricity_price_residual = total_data[‘electricity_price_per_mwh’] — electricity_price_trend

# Forecasting the trend component
forecast_steps = 10 # Number of steps to forecast ahead

# Forecasting future trend values
electricity_price_trend_forecast = electricity_price_results.forecast(steps=forecast_steps)

# Creating a date range for the forecast
electricity_price_trend_forecast_index = pd.date_range(start=electricity_price_trend.index[-1] + pd.Timedelta(days=1), periods=forecast_steps, freq=’D’)

# Creating a DataFrame for the forecasted trend
electricity_price_trend_forecast_df = pd.DataFrame(electricity_price_trend_forecast, index=electricity_price_trend_forecast_index, columns=[‘electricity_price_trend_forecast’])

# sum the predicted trend and predicted residual values to get the final forecast
forecast_df[‘electricity_price_residual_final’] = forecast_df[‘electricity_price_residual’] + electricity_price_trend_forecast_df[‘electricity_price_trend_forecast’]

# Plotting the final forecast of original electricity_price
plt.figure(figsize=(12, 6))
plt.plot(electricity_price[-100:], label=’Historical’, color=’blue’)
plt.plot(forecast_df.index, forecast_df[‘electricity_price_residual_final’], label=’Forecast’, color=’orange’, marker=’*’)
plt.title(‘Historical and Predicted Electricity Price’)
plt.xlabel(‘Date’)
plt.ylabel(‘Values’)
plt.legend()
plt.show()

10. Заключение и дополнительные ресурсы

В этой статье мы рассмотрели комплексный процесс построения и интерпретации модели векторной авторегрессии (VAR), специально разработанной для прогнозирования цен на электроэнергию с использованием соответствующих данных временных рядов.

Рекомендуемый ресурс для дальнейшего изучения:

  • «Прогнозирование: принципы и практика (3-е изд.)» Роба Хайндмана и Джорджа Атанасопулоса, OTexts. Эта книга является отличным ресурсом для всех, кто хочет углубить свое понимание методов прогнозирования, включая моделирование VAR. Он предлагает сочетание теоретических знаний и практического применения, что делает его подходящим как для начинающих, так и для опытных практиков в этой области.

Получите доступ к полному аналитическому путешествию

Полный код этого анализа временных рядов с использованием VAR-моделирования см. в моей записной книжке Jupyter на GitHub: VAR Modeling Repository.

Прогнозирование временных рядов с помощью трансформеров

Прогнозирование временных рядов является важнейшей задачей в различных областях, таких как финансы, энергетика, прогнозирование погоды и многое другое. Традиционные методы, такие как ARIMA, SARIMA и Prophet, широко используются для этой цели, но они могут испытывать трудности с улавливанием сложных временных зависимостей и долгосрочных закономерностей в данных.

С появлением глубокого обучения, в частности, архитектур трансформаторов, произошел сдвиг парадигмы в прогнозировании временных рядов, предлагающий многообещающие результаты в захвате сложных временных взаимосвязей.

Трансформеры, первоначально представленные для задач обработки естественного языка, показали значительный успех в моделировании последовательных данных благодаря своему механизму самовнимания и возможностям параллельной обработки.

Давайте посмотрим, как реализовать трансформаторы, чтобы сделать простую модель прогнозирования!

Прежде всего, давайте импортируем библиотеки,

import numpy as np
import pandas as pd
import torch
import torch.nn as nn
from torch.utils.data import DataLoader, TensorDataset
from sklearn.preprocessing import StandardScaler
from torch.optim.lr_scheduler import ReduceLROnPlateau

Теперь давайте импортируем набор данных,

names = ['year', 'month', 'day', 'dec_year', 'sn_value',
'sn_error', 'obs_num', 'unused1']
df = pd.read_csv(
"https://data.heatonresearch.com/data/t81-558/SN_d_tot_V2.0.csv",
sep=';', header=None, names=names,
na_values=['-1'], index_col=False)

Для этой статьи мы использовали набор данных солнечных пятен. Давайте взглянем на набор данных,


df.head()

Давайте выполним предварительную обработку данных,

# Data Preprocessing
start_id = max(df[df['obs_num'] == 0].index.tolist()) + 1
df = df[start_id:].copy()
df['sn_value'] = df['sn_value'].astype(float)
df_train = df[df['year'] < 2000]
df_test = df[df['year'] >= 2000]

spots_train = df_train['sn_value'].to_numpy().reshape(-1, 1)
spots_test = df_test['sn_value'].to_numpy().reshape(-1, 1)

scaler = StandardScaler()
spots_train = scaler.fit_transform(spots_train).flatten().tolist()
spots_test = scaler.transform(spots_test).flatten().tolist()

В приведенном выше коде мы разделяем данные для обучения и тестирования на основе столбца year. Мы изменяем их форму, а затем выполняем стандартизацию.

Теперь пришло время разбить временной ряд на последовательности,

# Sequence Data Preparation
SEQUENCE_SIZE = 10

def to_sequences(seq_size, obs):
x = []
y = []
for i in range(len(obs) - seq_size):
window = obs[i:(i + seq_size)]
after_window = obs[i + seq_size]
x.append(window)
y.append(after_window)
return torch.tensor(x, dtype=torch.float32).view(-1, seq_size, 1), torch.tensor(y, dtype=torch.float32).view(-1, 1)

Эта функция разделит данные на последовательность. Здесь размер последовательности равен 10.

Здесь мы используем первые 10 значений для прогнозирования 11-го значения. Это то, что делает описанная выше функция.

x_train, y_train = to_sequences(SEQUENCE_SIZE, spots_train)
x_test, y_test = to_sequences(SEQUENCE_SIZE, spots_test)

Давайте посмотрим, как устроены данные,

x_train[0]
y_train[0]

Теперь давайте посмотрим на следующий набор последовательностей.

x_train[1]

Мы видим, что предыдущий ответ является последним значением следующей последовательности. Последовательности формируются таким образом, чтобы модель могла обучаться.

Следующим шагом является определение наборов данных и создание загрузчиков данных,

# Setup data loaders for batch
train_dataset = TensorDataset(x_train, y_train)
train_loader = DataLoader(train_dataset, batch_size=32, shuffle=True)

test_dataset = TensorDataset(x_test, y_test)
test_loader = DataLoader(test_dataset, batch_size=32, shuffle=False)

Теперь все готово для обучения модели, кроме самой модели.

Создадим класс для модели.

В области архитектуры трансформаторов ключевым компонентом, обеспечивающим успех модели, является ее способность учитывать порядок последовательности. В отличие от традиционных RNN или LSTM, которые обрабатывают последовательности шаг за шагом и по своей сути соблюдают их порядок, преобразователи обрабатывают все токены в последовательности одновременно. В то время как эта параллельная обработка значительно повышает вычислительную эффективность и позволяет более эффективно захватывать дальние зависимости, это также означает, что преобразователи в своей исходной форме не обращают внимания на положение или порядок лексем в последовательности. Здесь в игру вступает концепция позиционного кодирования.

Позиционное кодирование — это механизм, предоставляющий преобразователю информацию о положении маркеров в последовательности. По сути, он вводит информацию о заказе в встраивания, не зависящие от позиции. При добавлении позиционных кодировок к вложениям маркеров перед их подачей в трансформатор положение каждой лексемы в последовательности становится различимым для модели.

Популярным методом генерации позиционных кодировок является использование синусоидальных функций.

Для заданной позиции «p» в последовательности и размерности «d» вложения позиционная кодировка вычисляется следующим образом:

# Positional Encoding for Transformer
class PositionalEncoding(nn.Module):
def __init__(self, d_model, dropout=0.1, max_len=5000):
super(PositionalEncoding, self).__init__()
self.dropout = nn.Dropout(p=dropout)

pe = torch.zeros(max_len, d_model)
position = torch.arange(0, max_len, dtype=torch.float).unsqueeze(1)
div_term = torch.exp(torch.arange(0, d_model, 2).float() * (-np.log(10000.0) / d_model))
pe[:, 0::2] = torch.sin(position * div_term)
pe[:, 1::2] = torch.cos(position * div_term)
pe = pe.unsqueeze(0).transpose(0, 1)
self.register_buffer('pe', pe)

def forward(self, x):
x = x + self.pe[:x.size(0), :]
return self.dropout(x)
  • d_model: размерность пространства встраивания модели.
  • dropout: Вероятность выпадения для регуляризации отсева (по умолчанию 0.1).
  • max_len: Максимальная длина последовательности для позиционного кодирования (по умолчанию 5000).

Теперь давайте определимся с моделью,

# Model definition using Transformer
class TransformerModel(nn.Module):
def __init__(self, input_dim=1, d_model=64, nhead=4, num_layers=2, dropout=0.2):
super(TransformerModel, self).__init__()

self.encoder = nn.Linear(input_dim, d_model)
self.pos_encoder = PositionalEncoding(d_model, dropout)
encoder_layers = nn.TransformerEncoderLayer(d_model, nhead)
self.transformer_encoder = nn.TransformerEncoder(encoder_layers, num_layers)
self.decoder = nn.Linear(d_model, 1)

def forward(self, x):
x = self.encoder(x)
x = self.pos_encoder(x)
x = self.transformer_encoder(x)
x = self.decoder(x[:, -1, :])
return x

device='cuda'
model = TransformerModel().to(device)
  • input_dim: Размерность входных данных, в данном случае мы используем только один вход — количество солнечных пятен.
  • d_model: Количество элементов во внутренних представлениях модели трансформатора (а также размер вложений). Этот параметр определяет, сколько данных модель может запоминать и обрабатывать.
  • nhead: Количество голов внимания в механизме самовнимания с несколькими головами.
  • num_layers: Количество слоев энкодера трансформатора. dropout: Вероятность отсева.

Далее необходимо определить оптимизатор, критерий и планировщик,

criterion = nn.MSELoss()
optimizer = torch.optim.Adam(model.parameters(), lr=0.001)
scheduler = ReduceLROnPlateau(optimizer, 'min', factor=0.5, patience=3, verbose=True)

Ладно, теперь все готово, приступим к тренировкам!

epochs = 15
early_stop_count = 0
min_val_loss = float('inf')

for epoch in range(epochs):
model.train()
for batch in train_loader:
x_batch, y_batch = batch
x_batch, y_batch = x_batch.to(device), y_batch.to(device)

optimizer.zero_grad()
outputs = model(x_batch)
loss = criterion(outputs, y_batch)
loss.backward()
optimizer.step()

# Validation
model.eval()
val_losses = []
with torch.no_grad():
for batch in test_loader:
x_batch, y_batch = batch
x_batch, y_batch = x_batch.to(device), y_batch.to(device)
outputs = model(x_batch)
loss = criterion(outputs, y_batch)
val_losses.append(loss.item())

val_loss = np.mean(val_losses)
scheduler.step(val_loss)

if val_loss < min_val_loss:
min_val_loss = val_loss
early_stop_count = 0
else:
early_stop_count += 1

if early_stop_count >= 5:
print("Early stopping!")
break
print(f"Epoch {epoch + 1}/{epochs}, Validation Loss: {val_loss:.4f}")

Оценка модели

# Evaluation
model.eval()
predictions = []
with torch.no_grad():
for batch in test_loader:
x_batch, y_batch = batch
x_batch = x_batch.to(device)
outputs = model(x_batch)
predictions.extend(outputs.squeeze().tolist())

rmse = np.sqrt(np.mean((scaler.inverse_transform(np.array(predictions).reshape(-1, 1)) - scaler.inverse_transform(y_test.numpy().reshape(-1, 1)))**2))
print(f"Score (RMSE): {rmse:.4f}")

Это дало бы среднеквадратичное значение прогнозов модели,

Источник

Источник

Источник

Источник

Источник

Источник

Источник

Источник

Источник

Источник

Источник

Источник

Источник

Источник