Получение значений AR (p) из модели ARMA в R

Question

Мне интересно получить значения части AR (p) из любой модели ARIMA в R.

Например, предположим, что я оцениваю следующую модель ARIMA

Мне интересно, чтобы эти значения были определены как:

Что я делаю, это:

  ar_model <- auto.arima(mydata) # fit an ARIMA model as decided by the auto.arima
  fitted_ar_model <- fitted(ar_model)
  resid_ar_model <- residuals(ar_model)
  ar_p <- fitted_ar_model - resid_ar_model

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

Answer 1

Я попытаюсь ответить здесь. К сожалению, написание формул в SO действительно громоздко. Итак, во-первых, я думаю, что вам нужно провести различие между «эпсилоном» вашего DGP (процесс генерации данных), также называемым шоками или ошибками (которые ненаблюдаемы), и моделью невязки . Шоки не видны, потому что вы не знаете истинных значений параметров модели (в данном случае параметров AR и MA) - вы можете только оценить их. Таким образом, существует разница (в простой модели AR (1) без константы для простоты написания с использованием нотации «hat» для оценок) между y_t - phi_1 * y_ {t-1} и y_t - \ hat {phi_1} * y_ {t-1} - первый является «шоком», а второй - «остаточным». Обычно вы также должны иметь обозначение \ hat {y_t} = \ hat {phi_1} * y_ {t-1} и, следовательно, остаточное значение, определенное как e_t = y_t - \ hat {y_t}.

Теперь так, учитывая я думаю, что вы хотели бы получить y_t - e_t - \ hat {phi} * e_ {t-1}, \ hat {phi} - это предполагаемый параметр MA1. Ниже приведен пример, где я явно подгонял ARMA (2,1) без константы к временному ряду LakeHuron. Я рассмотрел остатки модели, используя residuals (resid в df), а также рассчитал вручную (resid2 в df) - они немного отличаются из-за аппроксимации начальных условий, что метод residuals делает то, что я действительно не копался. Вы увидите, что значения становятся все ближе и ближе, когда t увеличивается, а влияние начальных условий уменьшается (это связано с тем, что у нас есть стационарная модель ARMA). Таким же образом я вычислил «AR_fit», как я назвал его один раз, как ar_1 * y_ {t-1} + ar_2 * y_ {t-2} (ar_fit1 in df) и второй раз как y_t - e_t - ma_1 * e_ {t-1} (ar_fit2 in df). Опять же, вы увидите, что значения сходятся при увеличении t (из-за тех же начальных условий, которые использовались для вычисления остатков с использованием функции residuals).

library(forecast)
data <-  LakeHuron - mean(LakeHuron)

nobs <- length(data)

ar_model <- arima(data, order = c(2, 0, 1), include.mean = FALSE)
fitted <-  fitted(ar_model)
resid <- residuals(ar_model)

ar1 <- ar_model$coef[1]
ar2 <- ar_model$coef[2]
ma1 <- ar_model$coef[3]

df <- data.frame(data, fitted, resid)

resid2 = data[3:nobs] - ar1*data[2:(nobs-1)] - ar2*data[1:(nobs-2)] - ma1*resid[2:(nobs-1)]
resid2 <- c(NA, NA, resid2)
df$resid2 <- resid2

ar_fit1 <- ar1*data[2:(nobs-1)] + ar2*data[1:(nobs-2)]
ar_fit1<- c(NA, NA, ar_fit1)
df$ar_fit1 <- ar_fit1   

ar_fit2 <- data[3:nobs] - resid[3:nobs] - ma1*resid[2:(nobs-1)]
ar_fit2<- c(NA, NA, ar_fit2)
df$ar_fit2 <- ar_fit2