Повышение надежности регрессионной модели с помощью анализа временных рядов — Часть 1

Улучшение надежности регрессионной модели с использованием анализа временных рядов — Часть 1

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

Фото Мухаммеда Фаиза Зулкефли на Unsplash

Введение

Находясь в 1,5 часах езды от дома, Сингапур всегда меня восхищал. Окруженный куда более крупными соседями, маленькая страна одержала победу над трудностями. Из скромного начала после обретения независимости, он стал важным финансовым центром и международным портом. Несмотря на свой успех, Сингапур сталкивается с серьезным вызовом – нехваткой земли. Эта проблема привела к чрезвычайно высоким ценам на жилье.

Ожидая эту проблему, сингапурское правительство стремится обеспечить достаточное жилье для своих жителей, создавая Жилищную разработочную доску (HDB) в 1960 году. Через эту организацию правительство сумело предоставить качественное, санитарное и доступное общественное жилье для более 80% населения, согласно официальному веб-сайту HDB.

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

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

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

Инструменты

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

import pandas as pd  # манипуляции структурированными данными
import numpy as np  # научные и численные вычисления
import scipy  # научные и численные вычисления
import matplotlib.pyplot as plt  # визуализация данных
import seaborn as sns  # визуализация данных
import statsmodels.api as sm  # статистический анализ
import datetime as dt  # манипуляции с датой и временем

Извлечение, очистка и предварительная обработка данных

(Для полного кода от начала до конца процесса анализа, от извлечения данных до построения моделей регрессии, не стесняйтесь заглянуть в мои репозитории на GitHub)

Набор данных по ценам перепродажи ХДБ, общедоступный и бесплатный для использования, может быть загружен с этого веб-сайта. После получения файла CSV мы можем загрузить данные в DataFrame Pandas, результатом чего будет следующая таблица:

# Загрузка набора данных в dataframe # Здесь я загрузил набор данных: https://data.gov.sg/dataset/resale-flat-prices# Я выбрал временной промежуток между январем 2017 года и августом 2023 годаdata = pd.read_csv('Housing_Resale_Dataset.csv')data.head()
Сырой набор данных - изображение автора

Сырой набор данных состоит из 160 795 строк (количества транзакций) и 11 столбцов (количества переменных). Поскольку мы сосредоточены на прогнозировании цен перепродажи ХДБ, столбец ‘resale_price’ становится зависимой переменной, или целевой, а остальные становятся независимыми переменными, или признаками. Хотя набор данных был уже относительно аккуратным, мы выполнили некоторую очистку и предварительную обработку данных, включая переименование столбцов, извлечение новых признаков, удаление дубликатов и перестановку для получения конечного чистого набора данных, который выглядит следующим образом:

cleandata.head()
Очищенный набор данных - изображение автора

Изображение выше показывает, что мы можем извлечь информацию из столбца ‘resale_date’ (ранее назывался ‘month’ и не виден на изображении), чтобы создать две новые переменные, ‘resale_year’ и ‘resale_month’, что дает нам 13 столбцов в очищенном наборе данных. В то же время, количество строк незначительно уменьшилось до 160 454, что является результатом удаления дубликатов.

Анализ данных

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

Распределение цен

Во-первых, давайте посмотрим на диаграмму распределения цен на HDB-продажи:

Распределение цен - изображение автора
cleandata['resale_price'].describe()
Сводная статистика - изображение автора

С бинами размером 40 000 мы видим, что распределение ‘resale_price’ смещено вправо, среднее значение составляет S$486,519 и медиана S$455,000.

Средняя цена по городу

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

Средняя цена по городу - изображение автора

Ежемесячная средняя цена

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

Ежемесячная средняя цена - изображение автора

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

Ежемесячные транзакции

На изображении ниже показаны ежемесячные транзакции с января 2017 года по август 2023 года. График также показывает, что данные о транзакциях имеют сезонность: в декабре и январе ежемесячные транзакции обычно снижаются. Мы также видим аномалию, когда транзакции резко снизились в первом квартале 2020 года. Можно с уверенностью предположить, что эта резкая падеж была вызвана блокировкой из-за Covid-19.

Ежемесячные транзакции - изображение автора

Анализ временных рядов

Фото Jake Hills на Unsplash

Как упоминалось во вступлении, наш акцент в этой истории на понимании ценообразования на HDB-квартиры. Таким образом, в качестве предварительного условия для анализа временных рядов мы сначала создадим наш набор данных временных рядов, создав объект Pandas с ежемесячными средними ценами и использовав «resale_date» в качестве индекса.

price_time = cleandata.groupby('resale_date')['resale_price'].mean()price_time.head()
Ежемесячное среднее значение цены - изображение автора

Декомпозиция

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

import statsmodels.api as sm# Мы используем "аддитивную" модель, так как сезонность, кажется, не усиливается со временем# Мы добавили 6 периодов экстраполяции, чтобы тенденция распространилась на первые и последние 6 месяцевprice_decomp = sm.tsa.seasonal_decompose(price_time, model = 'additive', extrapolate_trend=6) price_decomp.plot();
Декомпозиция временных рядов - изображение автора

С помощью простого скользящего среднего мы уже увидели тенденцию наших временных рядов. Между тем, выполняя декомпозицию, мы также можем видеть, что у данных есть сезонность в рамках 12-месячного цикла. Из графика выше видно, что сезонность незначительна по сравнению с средней ценой и находится в диапазоне от -4000 до 2500. Чтобы увидеть это лучше, мы выделим график сезонности и изменим линейный график на столбчатую диаграмму:

price_season = price_decomp.seasonal # сезонность price_season_month = price_season.reset_index()price_season_month['resale_month'] = price_season_month['resale_date'].dt.strftime('%b') # извлечь информацию о месяцахmonths_order = ['Янв', 'Фев', 'Мар', 'Апр', 'Май', 'Июн', 'Июл', 'Авг', 'Сен', 'Окт', 'Ноя', 'Дек']price_season_month['resale_month'] = pd.Categorical(price_season_month['resale_month'], categories= months_order, ordered = True) # создать порядок месяцев в данныхpm = price_season_month.groupby('resale_month')['seasonal'].mean()pm.plot(kind = 'bar')
Сезонность средней месячной цены - изображение автора

Разделение данных

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

# мы используем логарифмическое преобразование, потому что модель работает лучше таким образом
transform = np.logprice_time = price_time.apply(transform)# Разделение на обучающую и тестовую выборки
price_test = price_time['Сен 2022':] # тестовый набор данных (последние 12 месяцев)
price_train = price_time[:'Авг 2022'] # обучающий набор данных (предшествующие данные)

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

Моделирование временных рядов

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

Существует множество различных методов моделирования временных рядов. Классические методы включают скользящую среднюю, экспоненциальное сглаживание и ARIMA (авторегрессионное интегрированное скользящее среднее), в то время как более продвинутые техники включают Facebook Prophet и LSTM-сети. В данном контексте классические методы работают довольно хорошо, и среди них лучший результат показывает SARIMA (сезонное авторегрессионное интегрированное скользящее среднее), которое является расширением ARIMA.

Для построения модели SARIMA нам необходимо определить параметры: p, d и q модели ARIMA, а также P, D, Q и s параметры сезонности. Объяснение параметров можно прочитать здесь и здесь. Хотя существуют более строгие способы определения параметров, в данном случае мы используем прагматический подход с использованием сеточного поиска и выбираем модель на основе наименьшего значения критерия информационной Акаике (AIC).

from statsmodels.tsa.statespace.sarimax import SARIMAX # импорт функции SARIMA из библиотеки statsmodelss = 12 # сезонность состоит из цикла продолжительностью 12 месяцевfor p in range(1, 4):    for d in range(0, 3):        for q in range(0, 4):            for P in range(1, 4):                for D in range(0, 3):                    for Q in range(0, 4):                        model = SARIMAX(price_train, order=(p, d, q), seasonal_order=(P, D, Q, s))                        results = model.fit()                        print(f’ARIMA{p},{d},{q}-SARIMA{P},{D},{Q}-{s} - AIC: {results.aic}’)# обратите внимание, что выполнение сеточного поиска может потребовать значительных вычислительных ресурсов и занять некоторое время# поэтому выключите код, когда получите нужные параметры

Для этих данных наши оптимальные параметры – 3, 1 и 3 для параметров ARIMA и 1, 1, 3 и 12 для параметров сезонности. В блоке кода ниже мы можем увидеть, как подогнать обучающий набор данных под модель с использованием библиотеки Statsmodels.

p, d, q = 3, 1, 3P, D, Q, s = 1, 1, 3, 12sarima_model = SARIMAX(price_train, order=(p, d, q), seasonal_order=(P, D, Q, s))sarima_result = sarima_model.fit()

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

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

sarima_result.resid.plot()
График остатков модели — Изображение автора

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

sarima_result.resid[13:].plot()
Остаточный график модели без выбросов - изображение автора

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

Фактические данные против прогноза, анализ временных рядов - изображение автора

Оценка модели

Из приведенных выше графиков мы видим значительные сходства между прогнозом, полученным как на обучающем, так и на тестовом наборе данных, и фактическими значениями. Прогноз также тесно следует фактическим значениям в терминах паттерна и тренда. Следовательно, на основе визуального анализа, мы можем заключить, что наша модель SARIMA хорошо может прогнозировать данные на обучающем наборе данных. Более того, эта точность также распространяется на тестовый набор данных. Это означает, что модель может достаточно хорошо обобщать данные, которые не были видимы при обучении, по крайней мере, в короткий период времени (в нашем случае – 12-месячный период).

Помимо визуального анализа, мы также можем использовать определенные метрики для оценки точности нашей модели временных рядов. Одной из этих метрик является средняя абсолютная ошибка (MAE). Мы вычисляем MAE, взяв абсолютное значение средней разницы между предсказанными значениями и фактическими значениями в нашем наборе данных. В сущности, MAE показывает, насколько близко, в среднем, наша модель соответствует фактическим значениям. Более низкое значение MAE означает, что прогнозы нашей модели близки к фактическим данным, а более высокое значение MAE означает, что наша модель имеет большие ошибки. Для вычисления MAE мы можем легко использовать функцию, предоставленную библиотекой Scikit-learn.

from sklearn.metrics import mean_absolute_error as MAE# не забудьте перевернуть наше логарифмическое преобразование MAE_train = MAE(price_train[13:].apply(reverse_transform),train_predict[13:].apply(reverse_transform))MAE_test = MAE(price_test.apply(reverse_transform),test_predict.apply(reverse_transform))print(f'MAE обучающего набора данных = {MAE_train:.2f}')print(f'MAE тестового набора данных = {MAE_test:.2f}')
Средние абсолютные ошибки - изображение автора

Из вычислений выше мы получили значения MAE для обучающего набора данных и тестового набора данных в размере S$5748 и S$7549 соответственно. Чтобы поставить значения MAE в перспективу, мы можем сравнить их со средним значением наших данных временных рядов, равным S$481,070. Это означает, что значения MAE довольно малы по сравнению с среднемесячной ценой. Эти значения указывают на то, что наша модель имеет только незначительные отклонения от фактических значений.

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

Продолжение следует…

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

Дополнительная информация:

[1] Housing Development Board, HDB History and Towns, 2023. [Online]. Available: https://www.hdb.gov.sg/cs/infoweb/about-us/history [Accessed: October 30th, 2023]

[2] Джей Браунли, Как создать модель ARIMA для прогнозирования временных рядов с помощью Python, 2020 . [Online]. Доступно: https://machinelearningmastery.com/arima-for-time-series-forecasting-with-python/ [Просмотрено: 30 октября 2023]

[3] Джей Браунли, Представление SARIMA для прогнозирования временных рядов с помощью Python, 2019 . [Online]. Доступно: https://machinelearningmastery.com/sarima-for-time-series-forecasting-in-python/ [Просмотрено: 30 октября 2023]