Введение для начинающих в анализ и прогнозирование временных рядов

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

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

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

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

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

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

Код, сопровождающий эту статью, можно найти на моем GitHub здесь.

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

Чтобы правильно работать с данными временных рядов в R, мы сначала должны определить данные как объект временных рядов, используя функцию ts( ).

testing.stationarity = read.table("testing.stationarity.txt")
head(testing.stationarity)
class(testing.stationarity)

testing.stationarity = ts(testing.stationarity)
head(testing.stationarity)
class(testing.stationarity)

Данные точно такие же, просто хранятся как разные классы в R.

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

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

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

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

1. Ожидаемое значение (среднее) постоянно во времени
2. Волатильность (дисперсия) временного ряда постоянна вокруг своего среднего значения с течением времени.
3. Ковариация временного ряда зависит только от временного лага

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

Тренд – это долгосрочная траектория временного ряда. Например, на фондовом рынке, таком как индекс S&P 500, за последние несколько десятилетий наблюдается общая тенденция роста.

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

Существует два метода проверки стационарности временного ряда:

1. Единичный корневой тест Филлипса-Перрона
2. Построение образца ACF

Единичный корневой тест Филлипса-Перрона, PP.test( ), используется для проверки нулевой гипотезы о том, что временной ряд интегрирован первого порядка, другими словами, временной ряд требует дальнейшего дифференцирования для достижения стационарности.

PP.test(testing.stationarity)

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

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

acf(testing.stationarity, main="Data: testing.stationarity", ylab="Sample ACF")

Как видно из приведенного выше графика АКФ, автокорреляция медленно затухает, что указывает на то, что временной ряд не является стационарным.

Удаление тренда

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

Есть несколько способов удалить тренды из временного ряда. Здесь я описал два самых простых метода:

1. Дифференциация
2. Удаление трендов методом наименьших квадратов

1. Отличие

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

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

# Difference data using lag=1 and differences=1
Xt = diff(testing.stationarity, lag=1, differences=1)

# Plot original data, differenced data, and their respective sample ACFs 
par(mfrow=c(2,2))
ts.plot(testing.stationarity, main="Data: testing.stationarity",
        ylab="Value")
ts.plot(Xt, main="Differenced data", ylab="Change in value")
acf(testing.stationarity, main="", ylab="Sample ACF")
acf(Xt, main="", ylab="Sample ACF")
par(mfrow=c(1,1))

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

Давайте выполним тест на единичный корень Филлипса-Перрона, чтобы подтвердить наше наблюдение.

# Perform unit root test on differenced data 
PP.test(Xt)

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

Дифференциация: бонус!

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

# Try difference values between 0-7 and store their respective variance as a data frame 
ts.var = var(testing.stationarity)
for(i in 1:7) {
  diff = diff(testing.stationarity, lag=1, differences = i)
  ts.var[i+1] = var(diff)
}
ts.var.df = data.frame(diff=0:7, var=ts.var)

# Plot variance against the number of times data is differenced 
plot(ts.var.df, type="l", ylab="Variance", xlab="d")

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

2. Удаление тренда наименьших квадратов

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

Здесь мы имеем другой временной ряд с восходящей положительной тенденцией с течением времени.

# Generate time series data 
set.seed(123)
n = 1000
sim = arima.sim(list(ar=0.9), n)
xt = 2000+cumsum(sim)

# Plot generated time series 
ts.plot(xt, col="blue", main="Time series with trend", ylab="Data")

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

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

# Fit a linear model on time series, extract fitted values and residuals 
time = time(xt)
fit = lm(xt ~ time)
yt = fit$fitted.values
zt = fit$residuals

# Plot time series with superimposed linear model and residuals 
par(mfrow=c(2,1))
ts.plot(xt, col="blue", main="Regression example", ylab="Data")
abline(fit, col="red")
plot(xt-yt, type="l", col="green", xlab="Time", ylab="Residuals")
par(mfrow=c(1, 1))

Удаление сезонности

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

Здесь я продемонстрирую три метода устранения сезонности из временного ряда:

1. Сезонная разница
2. Сезонные средства
3. Метод скользящих средних

Прежде чем мы начнем, давайте кратко обсудим набор данных ldeaths в R. Он представляет ежемесячные цифры смертности от бронхита, эмфиземы и астмы в Великобритании в период с 1974 по 1979 год как для мужчин, так и для женщин.

# Plot ldeaths 
plot(ldeaths, main="Monthly deaths from lung diseases in the UK", ylab="Deaths")
points(ldeaths, pch=20)

# Add red vertical line at the start of each year  
abline(v=1974:1980, col="red")

# Plot sample ACF of ldeaths 
acf(ldeaths, main="Sample ACF of ldeaths", ylab="Sample ACF", lag.max=36)

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

1. Сезонная разница

Сезонное дифференцирование означает вычитание каждой точки данных на предыдущую точку данных с фиксированным запаздыванием.

# Difference ldeaths using lag=12 i.e. January 1975 minus January 1974, February 1975 minus February 1974, and so on 
sdiff.ldeaths = diff(ldeaths, lag=12, differences=1)

# Plot original data, differenced data, and their respective sample ACFs   
par(mfrow=c(2,2))
ts.plot(ldeaths, main="Data: ldeaths", ylab="Number of deaths")
acf(ldeaths, main="Sample ACF of ldeaths", ylab="Sample ACF")
ts.plot(sdiff.ldeaths, main="Data: sdiff.ldeaths", ylab="Difference in number of deaths")
acf(sdiff.ldeaths, main="Sample ACF of sdiff.ldeaths", ylab="Sample ACF")
par(mfrow=c(1,1))

2. Сезонные средства

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

# Generate ldeaths as dataframe 
ldeaths.df = data.frame(year=rep(1974:1979, each=12),
                        month=rep(1:12, 6),
                        value=ldeaths)
head(ldeaths.df, 12)

# Monthly averages of ldeaths dataset 
xbars = aggregate(value ~ month, data = ldeaths.df, mean)
xbars

# Subtract each month in ldeaths by their respective means
yt = ldeaths - xbars$value

# Plot ldeaths after subtracting seasonal means 
par(mfrow=c(2, 1))
plot(yt, main="Monthly deaths from lung diseases in the UK", ylab="Deaths")
points(yt, pch=20)
acf(yt, main="Sample ACF of the series ldeaths less seasonal means", ylab="Sample ACF", lag.max=36)
par(mfrow=c(1, 1))

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

3. Метод скользящих средних

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

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

# Decompose ldeaths into its trend, seasonal and random components 
plot(decompose(ldeaths))

# Store trend, seasonal and random as individuals variables  
decomp = decompose(ldeaths) 
trend = decomp$trend
seasonal = decomp$seasonal
random = decomp$random

# Plot data, trend and seasonal + trend  
ts.plot(ldeaths, ylab="", main="Components of time series: ldeaths", col="grey")
lines(trend, col="red")
lines(seasonal+trend, col="blue")
legend("topright", legend=c("Data", "Trend", "Seasonal + trend"), col=c("grey", "red", "blue"), lty=1)

Примерка модели

Модель ARIMA(p, d, q) является одной из наиболее часто используемых моделей временных рядов. Он состоит из трех гиперпараметров:

• p = Авторегрессия (AR)
• d = разность (I)
• q = скользящая средняя (MA)

Хотя в пакете forecast существуют такие функции, как auto.arima( ), которые могут помочь нам определить эти гиперпараметры, в этом разделе мы сосредоточимся на некоторых примерах, чтобы научиться вручную выбирать значения для p, d и q.

В частности, мы будем исследовать графики ACF и Partial ACF двух отдельных процессов временных рядов, чтобы определить гиперпараметры для наших моделей ARIMA.

Выбор правильной модели ARIMA: Пример 1

# Read time series data 
data = read.csv("fittingmodelEg1.csv", header=F)
data = ts(data[, 1])

# Plot data, ACF and partial ACF 
m = matrix(c(1, 1, 2, 3), 2, 2, byrow=TRUE)
layout(m)
ts.plot(data, ylab="")
acf(data,main="")
pacf(data,main="")
par(mfrow=c(1,1))

Ранее мы рассмотрели, что такое автокорреляционная функция (ACF), которая сообщает нам о корреляции ряда с его запаздывающими значениями. С другой стороны, функция частичной автокорреляции (PACF) измеряет корреляцию остатков следующего значения запаздывания, которые не учитывались при предыдущем запаздывании.

Глядя на приведенный выше временной ряд, мы видим, что очевидной тенденции нет, а дисперсия также кажется достаточно постоянной с течением времени.

АКФ не уменьшается медленно, что говорит о том, что нет необходимости в дальнейшем дифференцировании процесса, другими словами, установить d=0. Более того, АКФ также резко падает внутри доверительного интервала после лага 3. Это говорит о том, что мы должны установить гиперпараметр скользящего среднего равным 3, то есть q=3.

PACF, с другой стороны, затухает более постепенно.

Учитывая приведенную выше информацию, мы могли бы попробовать подобрать к этим данным модель MA (3) или, что то же самое, ARIMA (p = 0, d = 0, q = 3). На практике мы исследовали несколько моделей, прежде чем принять окончательное решение.

Выбор правильной модели ARIMA: Пример 2

# Read time series data 
data2 = read.csv("fittingmodelEg2.csv", header=F)
data2 = ts(data2[, 1])

# Plot data, ACF and partial ACF 
m = matrix(c(1, 1, 2, 3), 2, 2, byrow=TRUE)
layout(m)
ts.plot(data2, ylab="")
acf(data2,main="")
pacf(data2,main="")
par(mfrow=c(1,1))

Как и в примере 1, в самих данных нет очевидной тенденции, и переменная кажется достаточно постоянной во времени.

АКФ не показывает неуклонно снижающейся тенденции, поэтому данные кажутся стационарными. В результате мы можем установить d=0.

С другой стороны, PACF внезапно падает внутри доверительного интервала после задержки 2, что предполагает, что мы должны установить гиперпараметр авторегрессии равным 2, другими словами, p = 2.

Основываясь на приведенной выше информации, мы могли бы попытаться подогнать к этим данным модель AR (2) или, что то же самое, модель ARIMA (p = 2, d = 0, q = 0).

Проверка соответствия модели

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

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

1. График остатков
2. Тест Юнга-Бокса
3. Информационный критерий Акаике (AIC)

1. График остатков

График остатков, как следует из названия, отображает остатки данных временных рядов, которые не учитываются моделью.

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

См. ниже график остатков после подгонки модели MA(3) к примеру 1, который мы видели в предыдущем упражнении.

# Fit MA(3) model to data and extract residuals 
ma3 = arima(data, order=c(0, 0, 3))
residuals = ma3$residuals

# Plot residuals and ACF of residuals 
par(mfrow=c(2, 1))
ts.plot(residuals, main="MA(3) residuals", ylab="Residuals", col="blue")
acf(residuals, main="", ylab="ACF")
par(mfrow=c(1, 1))

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

В результате мы также можем заключить, что модель MA(3) обеспечивает хорошее соответствие данным.

2. Тест Льюнга-Бокса

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

Box.test(residuals, lag=5, type="Ljung", fitdf=3)

Из приведенного выше вывода, поскольку значение p больше 0,05, мы не отвергаем нулевую гипотезу и делаем вывод, что модель хорошо соответствует данным.

3. Информационный критерий Акаике (AIC)

И последнее, но не менее важное: AIC — это мера компромисса между качеством подгонки модели и количеством параметров, используемых в модели.

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

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

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

# Read time series data 
data3 = read.csv("fittingmodelEg3.csv", header=F)
data3 = ts(data3[, 1])

# Plot data without differencing and differenced data 
m = matrix(c(1, 1, 4, 4, 2, 3, 5, 6), 2, 4, byrow=TRUE)
layout(m)
ts.plot(data3, main="Data without differencing", ylab="")
acf(data3, main="", ylab="Sample ACF")
pacf(data3, main="", ylab="Sample PACF")
d = diff(data3)
ts.plot(d, main="Differenced data", ylab="")
acf(d, main="", ylab="Sample ACF")
pacf(d, main = "", ylab="Sample PACF")
par(mfrow=c(1, 1))

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

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

До сих пор мы знали, что наша модель ARIMA должна иметь ARIMA(p, d=1, q). Чтобы помочь нам определить значения p и q, мы развернем AIC. В частности, мы выберем комбинацию значений для p и q на основе той, которая дает нам самое низкое значение AIC.

# Try values 0-2 for both p and q, record their respective AIC and put them into a data sframe 
aic.result = numeric(3)
for (p in 0:2) {
  for (q in 0:2) {
    aic = arima(d, order=c(p, 0, q))$aic
    aic.result = rbind(aic.result, c(p, q, aic))
  }
}
aic.result = aic.result[-1, ]
colnames(aic.result) = c("p", "q", "AIC")
aic.result

Мы видим, что p = 2 и q = 2 дают самый низкий AIC, поэтому мы должны подогнать модель ARIMA (2, 1, 2) к этим конкретным данным.

Прогнозирование с использованием модели ARIMA

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

Для простоты давайте спрогнозируем на 100 шагов вперед от конечной точки данных в нашем исходном процессе временных рядов.

# Fit ARIMA(2, 0, 2) to differenced data, since data has already been differenced, we can set d=0 
fit = arima(d, order=c(2, 0, 2))
# Predict 100 steps ahead using ARIMA(2, 1, 2) model 
predictions = predict(fit, n.ahead=100)
predictions = predictions$pred

# Aggregate predictions with the final point of past data 
predictions.with.trend = tail(data3, 1) + cumsum(predictions)
predictions.with.trend = ts(predictions.with.trend, start=501, frequency=1)

# Plot past data and forecasts 
xlim = c(0, 600)
ylim = c(floor(min(data3, predictions.with.trend)), ceiling(max(data3, predictions.with.trend)))
ts.plot(data3, xlim=xlim, ylim=ylim, col="blue", main="Past data and forecasts", ylab="")
lines(predictions.with.trend, col="red")

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

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

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

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

Не знаете, что читать дальше? Вот некоторые предложения.