Modele ARMA w R

Najbardziej ogólnym modelem opisującym stacjonarny proces stochastyczny jest model ARMA (autoregressive moving average).

Modele $ARMA$ posiadają dwa parametry $p$ i $q$ dlatego oznaczamy je przez $ARMA(p,q)$.

Parametr $p$ jest rzędem autoregresji a $q$ rzędem średniej ruchomej.

Model ARMA(p,q) dla procesu $Y_{t}$ ma postać:

$Y_{t}=\alpha_1Y_{t-1}+\ldots+\alpha_pY_{t-p}+\epsilon_{t}-\beta_1\epsilon_{t-1}-\ldots-\beta_q\epsilon_{t-q}$,

gdzie $\epsilon_t$ to proces resztowy (biały szum).

Szczególnymi przypadkami modeli ARMA są:

• model autoregresyjny AR(p) - równocześnie ARMA(p,0)

• model średniej ruchomej MA(q) - równocześnie ARMA(0,q).

Proces ARMA często zapisuje operatorowo jako:

$A(u)Y_t=B(u)\epsilon_t$

gdzie

     $A(u)=1-\alpha_1u-\alpha_2u^2-\ldots-\alpha_pu^p$

oraz

     $B(u)==1-\beta_1u-\beta_2u^2-\ldots-\beta_qu^q$

Proces $ARMA(p,q)$ jest stacjonarny, jeżeli wszystkie pierwiastk zespolone wielomianu $A(u)$, co do modułu są większe od 1. Na stacjonarność nie ma wpływu średnia ruchoma, która zawsze jest stacjonarna.

W środowisku R do analizy modeli $ARMA$ służy biblioteka fArma. Przed jej pierwszym użyciem trzeba ją zainstalować (tylko raz)

install.packages("fArma")

a następnie załadować

library("fArma")

Przyjmijmy, że chcemy zasymulować proces $ARMA(2,3)$

$Y_t=0.12Y_{t-1}-0.4Y_{t-2}+\epsilon_t-0.9\epsilon_{t-1}-0.8\epsilon_{t-2}+0.1\epsilon_{t-3}$

W tym celu najpierw opiszemy model

model <- list( ar = c(0.12, -0.4), ma = c(0.9, 0.8, -0.1) )

Wektor ar określa parametry części autoregresyjnej a ma średniej ruchomej.

Określamy jak długi ma być szereg casowy

n <- 200

Szereg czasowy otrzymamy za pomocą polecenia

y <- armaSim(model=model, n=n)

Albo za pomocą jednego polecenia

y <- armaSim(model=list( ar = c(0.12, -0.4), ma = c(0.9, 0.8, -0.1) ), n=200)

Możemy również użyć polecenia ts aby ustalić inny zakres dat oraz częstotliwość danych

y_t <- ts(y, start=2000, frequency = 1)

Dzieki temu otrzymamy szereg danych rocznych (frequency=1) dla którego pierwsza obserwacja pochodzi z 2000 roku (start=2000).

Jeżeli chcemy utworzyć szereg czasowy dla danych kwartalnych (frequency=4) dla którego pierwsza obserwacja pochodzi z II kwartału 2001 roku  (start=c(2001 ,2)) należy użyć polecenia

y_t <- ts(y, start=c(2000,2), frequency = 4)

Do estymacji modeli $ARMA$ używamy polecenia armaFit

fit <- armaFit(formula = ~ arma(2, 3), data = y_t, include.mean=F)

Jeżeli parametr include.mean przyjmie wartość TRUE to szacowana będzie również stała.

Używając dodatkowo polecenia summary otrzymujemy 4 wykresy diagnostyczne wraz z pełniejszym opisem oszacowanego modelu

par(mfrow = c(2, 2))

summary(fit)

Jeżeli nie chcemy otrzymać wykresów

summary(fit, doplot=F)

Za pomocą polecenia coef możemy dotrzeć do poszczególnych parametrów modelu

coef(fit)[[3]]

Do prognozwoania (w przykładzi na 5 okresów) używamy polecenia predict

predict(fit, 5)