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)