Блог

Я работаю с arima0() и co2 . Я хотел бы построить arima0() модель по моим данным. Я пробовал fitted() и curve() безуспешно.

Вот мой код:

 ###### Time Series

# format: time series
data(co2)

# format: matrix
dmn <- list(month.abb, unique(floor(time(co2))))
co2.m <- matrix(co2, 12, dimnames = dmn)

co2.dt <- pracma::detrend(co2.m, tt = 'linear')
co2.dt <- ts(as.numeric(co2.dt), start = c(1959,1), frequency=12)

# first diff
co2.dt.dif <- diff(co2.dt,lag = 12)

# Second diff
co2.dt.dif2 <- diff(co2.dt.dif,lag = 1)
  

Подготовив данные, я выполнил следующее arima0 :

 results <- arima0(co2.dt.dif2, order = c(2,0,0), method = "ML")
resultspredict <- predict(results, n.ahead = 36)
  

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

Сессия 1: Для начала…

Честно говоря, я в значительной степени обеспокоен вашим подходом к моделированию co2 временных рядов. Что-то не так произошло уже тогда, когда вы изменили тренд co2 . Зачем использовать tt = "linear" ? Вы устанавливаете линейный тренд в пределах каждого периода (т. Е. года) и берете остатки для дальнейшей проверки. Часто это не рекомендуется, поскольку это приводит к появлению искусственных эффектов в остаточных рядах. Я бы склонялся к тому, чтобы сделать tt = "constant" , т. е. просто отбросить среднегодовое значение. Это, по крайней мере, сохранит корреляцию с сезоном, как в исходных данных.

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

 data(co2)

## de-trend by dropping yearly average (no need to use `pracma::detrend`)
yearlymean <- ave(co2, gl(39, 12), FUN = mean)
co2dt <- co2 - yearlymean

## de-trend by dropping within season linear trend
co2.m <- matrix(co2, 12)
co2.dt <- pracma::detrend(co2.m, tt = "linear")
co2.dt <- ts(as.numeric(co2.dt), start = c(1959, 1), frequency = 12)

## compare time series and ACF
par(mfrow = c(2, 2))
ts.plot(co2dt); acf(co2dt)
ts.plot(co2.dt); acf(co2.dt)
  

Оба ряда с отклонением от тренда имеют сильный сезонный эффект, поэтому требуется дополнительное сезонное различие.

 ## seasonal differencing
co2dt.dif <- diff(co2dt, lag = 12)
co2.dt.dif <- diff(co2.dt, lag = 12)

## compare time series and ACF
par(mfrow = c(2, 2))
ts.plot(co2dt.dif); acf(co2dt.dif)
ts.plot(co2.dt.dif); acf(co2.dt.dif)
  

ACF для co2.dt.dif имеет более значительные отрицательные корреляции. Это признак чрезмерного отклонения от тренда. Поэтому мы предпочитаем co2dt . co2dt уже является стационарным, и больше никаких различий не требуется (в противном случае вы просто перегружаете его и вводите больше отрицательной автокорреляции).

Большой отрицательный всплеск с задержкой 1 для ACF co2dt.dif предполагает, что нам нужна сезонная средняя. Кроме того, положительный всплеск в зависимости от сезона подразумевает умеренный процесс AR в целом. Итак, рассмотрим:

 ## we exclude mean because we found estimation of mean is 0 if we include it
fit <- arima0(co2dt.dif, order = c(1,0,0), seasonal = c(0,0,1), include.mean = FALSE)
  

Работает ли эта модель хорошо, нам нужно проверить ACF остатков:

 acf(fit$residuals)
  

Похоже, что эта модель приличная (на самом деле довольно отличная).

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

 fit <- arima0(co2dt, order = c(1,0,0), seasonal = c(0,1,1), include.mean = FALSE)
  

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

 ## 3 years' ahead prediction (no prediction error; only mean)
predco2dt <- predict(fit, n.ahead = 36, se.fit = FALSE)
  

Давайте построим co2dt подобранную модель и прогноз вместе:

 fittedco2dt <- co2dt - fit$residuals
ts.plot(co2dt, fittedco2dt, predco2dt, col = 1:3)
  

Результат выглядит очень многообещающим!

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

 fittedco2 <- fittedco2dt   yearlymean
  

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

 ts.plot(co2, fittedco2, col = 1:2)
  

Сессия 2: Лучшая идея для моделирования временных рядов

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

Сезонная структура co2 действительно сильна, поэтому нам в любом случае нужны сезонные различия:

 data(co2)
co2dt <- diff(co2, lag = 12)
par(mfrow = c(1,2)); ts.plot(co2dt); acf(co2dt)
  

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

 co2dt.dif <- diff(co2dt)
par(mfrow = c(1,2)); ts.plot(co2dt.dif); acf(co2dt.dif)
  

Отрицательные всплески внутри сезона и между сезонами предполагают, что процесс MA необходим для обоих. Я не буду работать с co2dt.dif ; мы можем работать с co2 напрямую:

 fit <- arima0(co2, order = c(0,1,1), seasonal = c(0,1,1))
acf(fit$residuals)
  

Теперь остатки совершенно некоррелированы! Итак, у нас есть ARIMA(0,1,1)(0,1,1)[12] модель для co2 рядов.

Как обычно, подогнанные значения получаются путем вычитания остатков из данных:

 co2fitted <- co2 - fit$residuals
  

Прогнозы выполняются одним вызовом predict :

 co2pred <- predict(fit, n.ahead = 36, se.fit = FALSE)
  

Давайте построим их вместе:

 ts.plot(co2, co2fitted, co2pred, col = 1:3)
  

О, это просто великолепно!

Сессия 3: Выбор модели

К настоящему времени история должна была закончиться; но я хотел бы провести сравнение с auto.arima from forecast , который может автоматически выбрать «лучшую» модель.

 library(forecast)
autofit <- auto.arima(co2)

#Series: co2 
#ARIMA(1,1,1)(1,1,2)[12]                    
#
#Coefficients:
#         ar1      ma1     sar1     sma1     sma2
#      0.2569  -0.5847  -0.5489  -0.2620  -0.5123
#s.e.  0.1406   0.1204   0.5880   0.5701   0.4819
#
#sigma^2 estimated as 0.08576:  log likelihood=-84.39
#AIC=180.78   AICc=180.97   BIC=205.5
  

auto.arima выбрал ARIMA(1,1,1)(1,1,2)[12] , который намного сложнее, поскольку включает как сезонные, так и несезонные различия.

Наша модель, основанная на пошаговом исследовании, предполагает ARIMA(0,1,1)(0,1,1)[12] :

 fit <- arima0(co2, order = c(0,1,1), seasonal = c(0,1,1))

#Call:
#arima0(x = co2, order = c(0, 1, 1), seasonal = c(0, 1, 1))
#
#Coefficients:
#      ma1     sma1
#  -0.3495  -0.8515
#s.e.   0.0497   0.0254
#
#sigma^2 estimated as 0.08262:  log likelihood = -85.98,  aic = 177.96
  

Значения AIC лучше отражают нашу модель. То же самое делает BIC:

 BIC = -2 * loglik   log(n) * p
  

У нас есть n <- length(co2) данные и p <- length(fit$coef) 1 параметры (дополнительный для sigma2 ), поэтому наша модель имеет BIC

 -2 * fit$loglik   log(n) * p
# [1] 196.5503
  

Итак, auto.arima данные перегружены.

На самом деле, как только мы видим ARIMA(1,1,1)(1,1,2)[12] , у нас возникают серьезные подозрения в его чрезмерной подгонке. Потому что разные эффекты «компенсируют» друг друга. Это происходит с дополнительными сезонными MA и несезонными AR, введенными auto.arima , поскольку AR вводит положительную автокорреляцию, в то время как MA вводит отрицательную.