Dekompozičné metódy kladú dôraz na prácu so systematickými zložkami – trendom, sezónnou a cyklickou. Na ich vyjadrenie pomocou regresnej analýzy slúži statický model, v ktorom je odozva \(y_t\) vysvetľovaná hodnotami regresorov v rovnakom čase, \(x_{jt}\).
Po odstránení systematických zložiek v časovom rade zostane iba reziduálna zložka, ktorá by mala byť stacionárna. Popis jej autokorelačnej štruktúry si už vyžaduje zostavenie dynamického modelu, v ktorom je časová závislosť vyjadrená prostredníctvom minulých hodnôt stochastického procesu.
Každý stochastický proces, ktorý neobsahuje žiadnu (dlhodobú) systematickú zložku, môže byť vyjadrený ako lineárna kombinácia postupnosti nekorelovaných, rovnako rozdelených náhodných premenných \(\varepsilon_t\). Nazýva sa lineárny proces, formálne \[ y_t-\mu=\varepsilon_t+\psi_1 \, \varepsilon_{t-1} + \psi_2 \, \varepsilon_{t-2} + ... = \sum_{i=0}^\infty \psi_i \, \varepsilon_{t-i}, \] pričom \(E(\varepsilon_t)=0\), \(var(\varepsilon_t)=\sigma_\varepsilon^2\), \(\psi_0=1\) a \(\sum_{i=0}^\infty |\psi_i|<\infty\). Lineárny proces je kauzálny, pretože \(y_t\) nezávisí od budúcich hodnôt (na rozdiel napr. od procesu klasických kĺzavých priemerov).
Lineárny proces je stacionárny, ak navyše platí \(\sum_{i=0}^\infty \psi_i^2<\infty\), ako to bude vidieť z vyjadrenia jeho rozptylu.
Pre lepšie porozumenie vlastnostiam lineárneho procesu si v nasledujúcom odvodíme jeho základné charakteristiky:
Ešte skôr než popíšeme modely lineárneho procesu, musíme spomenúť charakteristiku, ktorá úzko súvisí s ACF a používa sa na identifikáciu lineárneho procesu: parciálnu autokorelačnú funkciu (PACF).
Korelácia medzi dvoma náhodnými premennými je často spôsobená tým, že obidve tieto premenné sú korelované s treťou premennou. Veľká časť korelácie medzi premennými \(y_t\) a \(y_{t-k}\) môže byť spôsobená ich koreláciou s medziľahlými premennými \(y_{t-1}, y_{t-2},..., y_{t-k+1}\).
Parciálna autokorelačná funkcia je teda (v základnom tvare s \(\mu=0\)) definovaná ako podmienená stredná hodnota súčinu, \(\rho^*(k)=E(y_ty_{t-k}|y_{t-1},\ldots,y_{t-k+1})\) a podáva informáciu o korelácii premenných \(y_t\) a \(y_{t-k}\) očistené o vplyv \((k-1)\) premenných ležiacich v čase medzi nimi.
Hodnotu PACF pre posunutie \(k\) vyjadruje regresný koeficient \(\phi_{kk}\) v autoregresii \(k\)-teho stupňa: \[ y_t = \phi_{k1}y_{t-1} + \phi_{k2}y_{t-2} + \ldots + \phi_{kk}y_{t-k} + e_t, \] kde \(e_t\) je reziduálny člen nekorelovaný s premennými \(y_{t-j}\), pre \(j\geq1\). Postupnými úpravami \[ \begin{split} y_ty_{t-j} & = \phi_{k1}y_{t-1}y_{t-j} + \ldots + \phi_{kk}y_{t-k}y_{t-j} + e_ty_{t-j}\\ E(y_ty_{t-j}) & = \phi_{k1}E(y_{t-1}y_{t-j}) + \ldots + \phi_{kk}E(y_{t-k}y_{t-j}) + E(e_ty_{t-j})\\ \gamma(j) & = \phi_{k1}\gamma(j-1) + \ldots + \phi_{kk}\gamma(j-k) + 0 \quad \textrm{ak } \mu=0\\ \rho(j) & = \phi_{k1}\rho(j-1) + \ldots + \phi_{kk}\rho(j-k). \end{split} \] pre \(j=1,\ldots,k\) potom dostávame systém rovníc, z ktorého napr. použitím Cramerovho pravidla postupne pre \(k=1,2,\ldots\) získame hodnoty parciálnej autokorelačnej funkcie (PACF), \(\rho^*(k)=\phi_{kk}\).
Vo všeobecnosti môžeme ľubovoľný časový rad \(y_t\) zapísať (dynamicky) ako súčet dvoch častí – toho, čo môže a toho, čo nemôže byť vysvetlené s použitím vedomostí z minulosti (sústredených v množine \(\Omega_{t-1}\)) – teda
\[
y_t=E(y_t|\Omega_{t-1}) + e_t,
\] kde člen \(e_t\) sa nazýva chybový (angl. error term), pričom \(E(e_t|\Omega_{t-1}) = 0\).
Jeden – v praxi často používaný – model predpokladá, že
čiže \[ y_t = \varphi_1 y_{t-1} + ... + \varphi_p y_{t-p} + \varepsilon_t, \] kde \(\varphi_i\) sú reálne konštanty – parametre modelu, ktoré treba odhadnúť. Tento jednoduchý, no užitočný predpovedný model sa nazýva autoregresný model stupňa (stupňa) \(p\), skrátene AR(\(p\)).
Predpoklad modelu AR(\(p\)) o tom, že \(y_t\) závisí od \(p\) minulých hodnôt, môže znieť zavádzajúco, veď aj \(y_{t-p}\) závisí od \(y_{t-p-1},\ldots,y_{t-2p}\), takže v skutočnosti \(y_t\) závisí od všetkých minulých hodnôt procesu. Ich priamy vplyv na aktuálnu hodnotu však s časom odznieva, preto model používa iba tie relevantné, z nedávnej minulosti. Napr. ak postupne rozpíšeme AR(1) proces \[ \begin{split} y_1 &= \varphi_1 y_0 + \varepsilon_1 \\ y_2 &= \varphi_1 y_1 + \varepsilon_2 = \varphi_1\varphi_1 y_0 + \varphi_1\varepsilon_1 + \varepsilon_2 \\ \vdots \\ y_t &= (\varphi_1)^t y_0 + (\varphi_1)^{t-1}\varepsilon_1 + (\varphi_1)^{t-2}\varepsilon_2 + \ldots + \varepsilon_t \end{split} \] vidíme, že vplyv minulých pozorovaní je v konečnom dôsledku vplyvom náhodných šokov \(\varepsilon_t\) (teda skutočne ide o lineárny stochastický proces), ktoré odznievajú exponenciálnym tempom, podľa veľkosti \(\varphi_1\).
Ak je stupeň uvažovaného AR modelu veľký (takže treba odhadnúť veľa, povedzme \(m\) parametrov), vhodnou náhradou môže byť zmiešaný, autoregresný model s kĺzavými súčtami stupňa \(p,q\) (angl. Autoregressive Moving Average model of order \(p,q\)), čiže ARMA(\(p,q\)), \[ y_t = \varphi_1 y_{t-1} + ... + \varphi_p y_{t-p} + \varepsilon_t + \theta_1 \varepsilon_{t-1} + ... + \theta_q \varepsilon_{t-q}, \] v ktorom \(p+q<m\). Pri niektorých časových radoch je dokonca vhodnejší jeho druhý špeciálny prípad, model kĺzavých súčtov MA(\(q\)), \[ y_t = \varepsilon_t + \theta_1 \varepsilon_{t-1} + ... + \theta_q \varepsilon_{t-q}, \] s parametrami \(\theta_1,\ldots,\theta_q\).
V nasledujúcej časti charakterizujeme dve dôležité vlastnosti zmiešaného lineárneho modelu.
Ak definujeme operátor spätného posunu \(B^j\) (angl. backshift), s ktorým platí \(B^j y_t = y_{t-j}\), potom ARMA(\(p,q\)) možno prepísať do jednoduchého tvaru \[ \varphi(B) y_t = \theta(B) \varepsilon_t, \] kde \(\varphi(x) = 1 - \varphi_1 x - \varphi_2 x^2 - ... - \varphi_p x^p\) je autoregresný (AR) polynóm, a \(\theta(x) = 1 + \theta_1 x + \theta_2 x^2 + ... + \theta_q x^q\) zas polynóm kĺzavých súčtov (skrátene MA polynóm).
V príklade z predošlej časti sme videli, že AR(1) predstavuje jednoduchý lineárny stochastický proces, pre ktorý z Woldovej reprezentácie vyplýva \(\psi_i=(\varphi_1)^i\). Aby bol AR(1) proces stacionárny, musí byť účinok minulých šokov \(\varepsilon_{t=i}\) tlmený, konkrétne musí platiť \(|\varphi_1|<1\).
Pre AR (a zmiešaný) proces vyššieho stupňa, je podmienka stacionarity iba trochu zložitejšia: proces ARMA(\(p, q\)) je stacionárny vtedy a len vtedy, ak všetky nulové body AR polynómu \(\varphi(x)\) ležia v komplexnej rovine mimo jednotkového kruhu. Teda pre všetky korene \(1 - \varphi_1 x - \dots - \varphi_p x^p=0\) platí \(|x|>1\). Až potom ku \(\varphi(B)\) existuje \(\varphi(B)^{-1}\) a \[ \begin{split} \varphi(B) y_t &= \theta(B) \varepsilon_t \quad /\varphi(B)^{-1} \\ y_t &= \varphi(B)^{-1}\theta(B) \varepsilon_t = \sum_{j=0}^\infty\psi_j \varepsilon_{t-j}, \end{split} \] takže stacionárny proces ARMA(\(p, q\)) môžeme formálne vyjadriť ako MA(\(\infty\)).
Duálnou vlastnosťou ku stacionarite je tzv. invertibilita: proces ARMA(\(p, q\)) je invertibilný (čiže \(\varepsilon_t\) možno vyjadriť ako lineárnu kombináciu súčasnej a minulých hodnôt \(y_t\)) práve vtedy, ak všetky nulové body MA polynómu \(\theta(x)\) ležia v komplexnej rovine mimo jednotkového kruhu. Vďaka tomu invertibilný proces ARMA(\(p, q\)) môžeme formálne vyjadriť ako AR(\(\infty\)).
Teraz podrobnejšie rozoberieme vlastnosti AR a MA modelov, aby sme zistili, ako sa dajú identifikovať v stochastickom procese.
Autoregresný proces prvého stupňa, AR(1) v tvare \(y_t = \varphi_1 y_{t - 1} + \varepsilon_t\) alebo \((1 - \varphi_1 B) y_t = \varepsilon_t\), je stacionárny, ak pre koreň rovnice \(\varphi(X)=0\), čiže rovnice \(1-\varphi_1 x=0\), platí \(|x|>1\), teda \(|\frac1\varphi_1|>1\) a tak \(|\varphi_1|<1\).
Ako sme už ukázali, stacionárny AR(1) sa dá vyjadriť v tvare stacionárneho lineárneho procesu \(y_t=\varepsilon_t+\varphi_1 \varepsilon_{t-1}+\varphi_1^2 \varepsilon_{t-2}+ ... = \sum_{k=0}^\infty \varphi_1^k \varepsilon_{t-k}\). V nasledujúcom sú odvodené jeho základné charakteristiky.
Pre \(0<\varphi_1 <1\) funkcia \(\rho(k)\) exponenciálne klesá k nule. Pre \(-1<\varphi_1 <0\) tiež klesá k nule, ale oscilačne (Obr. 3.1).
# helper function for subsequent plots
ggARMAacf <- function(x, ylab = c("ACF", "PACF")) {
x |> rbind(value = _) |> t() |> as.data.frame() |>
tibble::rownames_to_column("lag") |> dplyr::mutate(lag = as.integer(lag)) |>
ggplot() + aes(x = lag, y = value) + geom_line() + geom_point() +
geom_hline(yintercept = 0, linetype = 3) +
scale_x_continuous(breaks = scales::breaks_pretty()) + # ensure integers
ylab(ylab[1])
}ARMAacf(ar=c(0.8), lag.max = 20) |> ggARMAacf()
ARMAacf(ar=c(-0.8), lag.max = 20) |> ggARMAacf()
Autoregresný proces stupňa \(p\), \[ y_t = \varphi_1 y_{t - 1} + \varphi_2 y_{t - 2} +...+\varphi_p y_{t - p} + \varepsilon_t \] je vždy invertibilný, ale stacionárny je len vtedy, ak všetky korene \(\varphi(x)=0\) ležia mimo jednotkového kruhu.
Charakteristiky sú dané nasledujúcimi vzťahmi:
Pre reálne korene rovnice \(\varphi(x) = 0\) je ACF zložená z exponenciálne klesajúcich kriviek, pre komplexné korene zas z exponenciálne klesajúcich sínusových kriviek.
Proces kĺzavých súčtov stupňa \(q\) \[
y_t = \varepsilon_t + \theta_1 \, \varepsilon_{t-1}+\theta_2 \, \varepsilon_{t-2}+...+\theta_q \, \varepsilon_{t-q}
\]
je vždy stacionárny. A invertibilný je práve vtedy, ak všetky korene \(\theta(x)=0\) ležia mimo jednotkového kruhu. Pre charakteristiky platí:
Priebeh parciálnej autokorelačnej funkcie na príklade MA(1) ilustruje Obr. 3.2.
ARMAacf(ma = c(0.8), lag.max=20, pacf=T) |> ggARMAacf()
ARMAacf(ma = c(-0.8), lag.max=10, pacf=T) |> ggARMAacf()
V prípade AR(\(p\)) platí, že PACF od \(k>p\) klesne k nule, \[ \begin{split} \rho^*(k) &= 0 \qquad \text{pre } k > p, \\ \rho(k) &\neq 0,\qquad \forall k. \end{split} \] Pre MA(\(q\)) je to naopak, ku nule klesne ACF, \[ \begin{split} \rho(k) &= 0 \qquad \text{pre } k > q, \\ \rho^*(k) &\neq 0,\qquad \forall k. \end{split} \] Bod \(k\), za ktorým sú hodnoty ACF resp. PACF nulové, označíme súhrnne \(k=k_0\) (niekedy sa nazýva bodom useknutia).
V nasledujúcich príkladoch si všimnime bod useknutia \(k_0=1\) pri prvých 4 procesoch, a \(k_0=2\) pri ostatných.
kmax <- 10
p1 <- ARMAacf(ar=c(0.7), lag.max = kmax) |> ggARMAacf()
p2 <- ARMAacf(ar=c(0.7), lag.max = kmax, pacf = TRUE) |> ggARMAacf(ylab = "PACF")
p1 + p2
p1 <- ARMAacf(ar=c(-0.7), lag.max = kmax) |> ggARMAacf()
p2 <- ARMAacf(ar=c(-0.7), lag.max = kmax, pacf = TRUE) |> ggARMAacf(ylab = "PACF")
p1 + p2
p1 <- ARMAacf(ma=c(0.7), lag.max = kmax) |> ggARMAacf()
p2 <- ARMAacf(ma=c(0.7), lag.max = kmax, pacf = TRUE) |> ggARMAacf(ylab = "PACF")
p1 + p2
p1 <- ARMAacf(ma=c(-0.7), lag.max = kmax) |> ggARMAacf()
p2 <- ARMAacf(ma=c(-0.7), lag.max = kmax, pacf = TRUE) |> ggARMAacf(ylab = "PACF")
p1 + p2
p1 <- ARMAacf(ar=c(0.9, -0.8), lag.max = kmax) |> ggARMAacf()
p2 <- ARMAacf(ar=c(0.9, -0.8), lag.max = kmax, pacf = TRUE) |> ggARMAacf(ylab = "PACF")
p1 + p2
p1 <- ARMAacf(ma=c(0.2, -1.2), lag.max = kmax) |> ggARMAacf()
p2 <- ARMAacf(ma=c(0.2, -1.2), lag.max = 50, pacf = TRUE) |> ggARMAacf(ylab = "PACF")
p1 + p2
Parciálna ACF na Obr. 3.8 konverguje ku nule veľmi pomaly, pretože proces MA(2) nie je invertibilný, čo vidieť aj na koreňoch MA plynómu – jeden z nich leží vo vnútri jednotkového kruhu.
c(1, 0.2, -1.2) |> polyroot() |> round(2)[1] 1.00+0i -0.83+0iUkázali sme, že ak reziduálnu zložku možno považovať za lineárnu kombináciu postupnosti nekorelovaných rovnako rozdelených náhodných premenných \(\varepsilon_t\), vieme ju popísať triedou modelov ARMA, ktoré určíme pomocou tzv. Box-Jenkinsovej metodológie. Je to relatívne starý matematický aparát, prvý raz bol popísaný v r. 1976, no dodnes sa často uplatňuje pri modelovaní ČR. Ako sme videli, za základ analýzy berie reziduálnu zložku tvorenú korelovanými (alebo závislými) náhodnými premennými.
Medzi výhody tejto metodológie patria nasledujúce:
Ako nevýhody sa uvádza napr. to, že
Praktický postup určenia modelu daného časového radu v rámci Box-Jenkinsovej metodológie zahŕňa tri kroky:
Identifikovať model znamená určiť vonkajšie parametre (nazývané aj hyper-parametre), teda tie, ktoré sa nedajú odhadnúť klasickou optimalizačnou metódou spolu s vnútornými parametrami. V kontexte modelov ARMA ide o stupne AR a MA polynómov (\(p,q\)). Používajú sa dve metódy:
V neparametrickom prístupe využijeme informácie z predošlej podkapitoly a na základe priebehu výberovej ACF \(\hat\rho(k)\) a PACF \(\hat\rho^*(k)\) určíme typ a stupeň vhodného modelu, pre zhrnutie pozri Tab. 3.1.
Bod useknutia \(k_0\) tu predstavuje najväčšiu hodnotu časového posunu, pre ktorý je ACF, resp. PACF ešte nenulová. Na určenie \(k_0\) v teste nulovosti hodnôt ACF sa môže použiť presnejšia, Bartlettova aproximácia rozptylu výberovej autokorelácie, ktorá (namiesto bieleho šumu) predpokladá proces MA(\(k_0\)). Potom \[ |\hat\rho(k)|<2\sqrt{\frac1n(1+2\sum_{i=1}^{k_0}\hat\rho(j)^2)}\quad\text{pre } k>k_0 \] na hladine významnosti 5%.
Pre test nulovosti PACF je k dispozícii iba doteraz používaná, tzv. Quenouilleova aproximácia, ktorá v nulovej hypotéze predpokladá biely šum. Takže ak \(H_0\) platí, potom asymptoticky \[ \hat\rho^*(k) \sim N\left(0, \frac1n\right) \quad\text{pre } k>k_0. \]
| funkcia | AR(p) | MA(q) | ARMA(p,q) |
|---|---|---|---|
| ACF \(\rho(k)\) | neexistuje \(k_0\), tvar U krivky | \(k_0=q\) | neexistuje \(k_0\), po \(\max(0,q-p)\) tvar U krivky |
| PACF \(\rho^*(k)\) | \(k_0=p\) | neexistuje \(k_0\), ohraničená U krivkou | neexistuje \(k_0\), po \(\max(0,p-q)\) ohraničená U krivkou |
Prirodzene, ak je ACF všade blízke 0, časový rad je pravdepodobne realizáciou procesu bieleho šumu. Naopak, ak je pokles k nule veľmi pomalý, proces je nestacionárny.
V parametrickom prístupe sa najprv určia horné hranice \(p_{max}\) a \(q_{max}\) pre správny stupeň modelu ARMA(\(p,q\)). Ak by bol odhad hyperparametrov \(p,q\) založený iba na minimalizácii odhadu rozptylu \(\hat{\sigma}_{p,q}^2\) chybovej zložky procesu ARMA(\(p,q\)), preferoval by neúmerne vysoké hodnoty. Namiesto toho je lepšie minimalizovať funkciu \[ A_{p,q} = (1 + w_n(p, q) )*\hat{\sigma}_{p,q}^2, \qquad p = 0, 1, ..., p_{max}, \quad q = 0, 1, ..., q_{max}, \] kde \(w_n(p, q)\) je tzv. penalizačná funkcia, ktorá znevýhodňuje voľbu príliš veľkého stupňa modelu a musí mať tieto dve vlastnosti:
Voľbou \(w_n(p, q)\) a zlogaritmovaním \(A_{p,q}\) vznikne konkrétne kritérium, najčastejšie používané sú:
Ak sa z nejakého dôvodu obmedzíme iba na AR, stupeň \(p\) možno určiť aj subjektívne, z grafu reziduálneho rozptylu odhadnutých modelov AR(1), AR(2) …, kedy disperzia od kritického bodu už prudko neklesá.
(TO DO: urobiť príklady na základe simulovaných procesov)
Po identifikácii modelu (t.j. určení triedy modelu a jeho hyperparametrov) prichádza k slovu odhad jeho parametrov. Pripomeňme požiadavku, aby bol časový rad stacionárny a s nulovou strednou hodnotou. V opačnom prípade môžeme získať nekorektný model, prípadne výpočet môže byť numericky nestabilný.
Na odhad parametrov AR modelov slúži napr.
parametre všeobecnejších modelov ARMA sa odhadujú pomocou
Používa sa na odhad parametrov modelu AR(\(p\)) pre dané \(p\). Vychádza z rovnice AR procesu, ktorá postupnými úpravami
\[ \begin{split} y_t - \varphi_1 y_{t-1} - \ldots - \varphi_p y_{t-p} & = \varepsilon_t \\ y_t y_{t-k} - \varphi_1 y_{t-1} y_{t-k} - \ldots - \varphi_p y_{t-p} y_{t-k} & = \varepsilon_t y_{t-k} \\ E(y_t y_{t-k}) - \varphi_1 E(y_{t-1} y_{t-k}) - \ldots - \varphi_p E(y_{t-p} y_{t-k}) & = E(\varepsilon_t y_{t-k}), \end{split} \qquad \begin{split} / & y_{t-k} \\ / & \text{E} \\ k &= 0,1,...,p \end{split} \] vedie na systém (\(p+1\)) rovníc o (\(p+1\)) neznámych1 \[ \begin{split} \gamma(k) - \varphi_1 \gamma(k-1) -\varphi_2 \gamma(k-2) -...- \varphi_p \gamma(k-p) & = 0, \\ \gamma(0) - \varphi_1 \gamma(1) -\varphi_2 \gamma(2) -...- \varphi_p \gamma(p) & = \sigma_{\varepsilon}^2, \end{split} \qquad \begin{split} k & = 1, 2, ..., p, \\ k & = 0, \end{split} \] ktorý sa rieši finitnou metódou pomocou inverznej matice. V maticovom zápise \[ \boldsymbol{\Gamma\varphi}=\boldsymbol\gamma,\quad \sigma_{\varepsilon}^2=\gamma(0)-\boldsymbol{\varphi'\gamma}, \] kde \(\boldsymbol\Gamma = \{\gamma(k-j)\}_{j,k=1}^p\) je \(p\times p\) matica, \(\boldsymbol\varphi=(\varphi_1,\ldots,\varphi_p)'\) a \(\boldsymbol\gamma=(\gamma(1),\ldots,\gamma(p))'\). Potom metódou momentov a úpravami, kedy namiesto teoretickej použijeme výberovú autokovariančnú funkciu \(\hat\gamma(k)\), dostávame odhady \(\hat{\boldsymbol\varphi}=\hat\Gamma^{-1}\hat{\boldsymbol\gamma}\) a \(\sigma_{\varepsilon}^2=\hat\gamma(0)-\hat{\boldsymbol\gamma}'\hat\Gamma^{-1}\hat{\boldsymbol\gamma}\), alebo analogicky pomocou ACF \[ \hat{\boldsymbol\varphi}=\hat{\mathbf{R}}^{-1}\hat{\boldsymbol\rho}, \qquad \hat{\sigma}_{\varepsilon}^2=\hat\gamma(0)\left[1-\hat{\boldsymbol\rho}'\hat{\mathbf{R}}^{-1}\hat{\boldsymbol\rho}\right], \] kde \(\mathbf{R}=\boldsymbol{\Gamma}/\gamma(0)\).
Je založený na rekurzívnom riešení Yule – Walkerových rovníc. Počíta vhodné modely AR(1), AR(2), …, AR(\(p_{max}\)), kde \(p_{max}\) volíme, pričom skutočný stupeň \(p \leq p_{max}\).
Výhodou Levinson–Durbinovho algoritmu je, že okrem priameho riešenia YW rovníc (bez použitia inverznej matice) získame aj odhady PACF – výstupom je \(p_{max}\) modelov AR(\(p\)), \(p=1,...,p_{max}\), pričom posledný koeficient každého z nich predstavuje hodntu PACF \(\hat\rho^*(p)\). To umožní vybrať vhodný stupeň AR priamo v procese odhadu.
Používa sa na odhad koeficientov modelu ARMA(\(p,q\)), ale aj AR(\(p\)) (pre \(q = 0\)), resp. MA(\(q\)) (pre \(p = 0\)). Využíva fakt, že invertibilný model ARMA(\(p,q\)) môžeme aproximovať modelom AR(\(k\)) pre dostatočne veľké \(k\).
Najprv odhadneme model AR(\(k\)) pre dostatočne veľké \(k\): \[ y_t=\phi_1 \, y_{t-1}+\phi_2 \, y_{t-2}+...+\phi_k \, y_{t-k}+z_t \tag{3.1}\] Odhady koeficientov \(\hat{\phi}_1,\ldots,\hat{\phi}_k\) získame klasickou metódou najmenších štvorcov. Tieto odhady dosadíme do (3.1), čím získame odhady rezíduí \[ \hat{z}'_t=y_t-\hat{\phi}_1 \, y_{t-1}-\ldots-\hat{\phi}_k \, y_{t-k}, \quad t=k+1,...,n. \] Po dosadení do modelu ARMA(\(p,q\)), \[ y_t-\varphi_1 \, y_{t-1}-\varphi_2 \, y_{t-2}-\ldots-\varphi_p \, y_{t-p}=\varepsilon_t+\theta_1 \, \hat{z}'_{t-1}+\ldots+\theta_q \, \hat{z}'_{t-q} \] klasickou metódou najmenších štvorcov odhadneme parametre \(\varphi_1,\ldots,\theta_q\) modelu ARMA(\(p,q\)). Rozptyl chybovej zložky bude \[ \hat{\sigma}_{\varepsilon}^2=\sum_{t=t_1+1}^n \frac{\hat{\varepsilon}_t^2}{n-t_1}, \quad\text{ kde } t_1=\max (q+k,p). \]
Rozširuje Long AR metódu o určenie stupňov \(k_{max},p,q\) pomocou informačných kritérií. Pozostáva z nasledujúcich krokov:
Predpokladáme, že biely šum \(\varepsilon_t\) stacionárneho ARMA(\(p,q\)) procesu s nulovou strednou hodnotou je gaussovský, teda platí \(\varepsilon_t \sim N(0, \sigma_\varepsilon^2 )\). Združená hustota pravdepodobnosti náhodného vektora \(\boldsymbol{\varepsilon}=(\varepsilon_1,..., \varepsilon_n)\) je \[ f(\boldsymbol{\varepsilon})=(2 \pi \sigma_\varepsilon^2)^{-\frac{n}{2}} \; e^{-\frac{1}{2\sigma_\varepsilon^2} \sum_{t=1}^n \varepsilon_t^2}. \] Dosadením \(\varepsilon_t = y_t-\varphi_1y_{t-1}-\ldots-\varphi_py_{t-p}-\theta_1z_{t-1}-\ldots-\theta_qz_{t-q}\) dostávame hustotu ako funkciu parametrov \(\varphi_1,\ldots,\theta_q\) (a \(\sigma_\varepsilon^2\)). Zlogaritmovaním získame tzv. funkciu vierohodnosti (angl. likelihood function) \[ L(\boldsymbol{\varphi},\boldsymbol{\theta},\sigma_\varepsilon^2) = -\frac{n}{2} \ln (2 \pi) -\frac{n}{2} \ln (\sigma_\varepsilon^2)-\frac{1}{2\sigma_\varepsilon^2} \sum_{t=1}^n \varepsilon_t^2 \] ktorej maximalizáciou (ekvivalentne minimalizáciou reziduálneho súčtu štvorcov \(S(\boldsymbol{\varphi},\boldsymbol{\theta})=\sum_{t=1}^n \varepsilon_t^2\)) dostaneme odhady parametrov \(\hat\varphi_1,\ldots,\hat\theta_q\) a \(\hat\sigma_\varepsilon^2=\frac{S(\boldsymbol{\hat\varphi},\boldsymbol{\hat\theta})}{n-p-q}\). Pretože hodnoty \(y_0, y_{-1}, ...\) a \(\varepsilon_0, \varepsilon_{-1},...\) nie sú známe, začneme od \(t=p+1\) a v prípade \(q>p\), neznáme hodnoty \(\{\varepsilon_t\}_{t\leq0}\) nahradíme strednou hodnotou, čiže 0 (hovoríme o podmienenej metóde max.vierohodnosti), inak by sme museli spätne extrapolovať.
Metódou maximálnej vierohodnosti získame nevychýlené efektívne odhady parametrov.
Napríklad pre ARMA(1,1), čiže \(\varepsilon_t=y_t-\varphi_1y_{t-1}-\theta_1\varepsilon_{t-1}\) bude vierohodnostná funkcia v tvare \[ L(\boldsymbol{\varphi},\boldsymbol{\theta},\sigma_\varepsilon^2) = -\frac n2 \ln(2\pi) - \frac n2 \ln(\sigma_\varepsilon^2) - \frac1{2\sigma_\varepsilon^2}\sum_{t=1}^n(y_t-\varphi_1y_{t-1}-\theta_1\varepsilon_{t-1})^2 \] a štandardným postupom hľadania lokálnych extrémov \[ \begin{split} \frac{\partial L(\boldsymbol{\varphi},\boldsymbol{\theta},\sigma_\varepsilon^2)}{\partial\varphi_1} = -\frac1{\sigma_\varepsilon^2} \sum_{t=1}^n(y_t-\varphi_1y_{t-1}-\theta_1\varepsilon_{t-1})(-y_{t-1}) & \\ \sum_{t=1}^n(y_t-\hat\varphi_1y_{t-1}-\hat\theta_1\varepsilon_{t-1})(-y_{t-1}) &= 0 \\ \hat\varphi_1\sum_{t=1}^ny_{t-1}^2 + \hat\theta_1\sum_{t=1}^n\varepsilon_{t-1}y_{t-1} &= \sum_{t=1}^ny_ty_{t-1} \\ \frac{\partial L(\boldsymbol{\varphi},\boldsymbol{\theta},\sigma_\varepsilon^2)}{\partial\theta_1} = 0 \quad\Rightarrow \quad \hat\varphi_1\sum_{t=1}^ny_{t-1}\varepsilon_{t-1} + \hat\theta_1\sum_{t=1}^n\varepsilon_{t-1}^2 & = \sum_{t=1}^ny_t\varepsilon_{t-1} \end{split} \] dostávame dve rovnice o dvoch neznámych \(\hat\varphi_1, \hat\theta_1\).
Na záver ešte spomenieme asymptotické vlastnosti maximálne vierohodných odhadov.
Nech \(\boldsymbol{\beta} = (\varphi_1, ..., \varphi_p, \theta_1, ..., \theta_q)^T\) sú skutočné parametre stacionárneho a invertibilného procesu ARMA(, ) a \(\hat{\boldsymbol{\beta}}\) je odhad vektora \(\boldsymbol{\beta}\), získaný metódou maximálnej vierohodnosti alebo podmienenou metódou maximálnej vierohodnosti. Pre \(n \rightarrow \infty\) platí \(\sqrt{n}\; \left(\hat{\boldsymbol{\beta}} - \boldsymbol{\beta} \right)\sim N(0, \mathbf{V}(\beta)),\) kde \(\mathbf{V}(\beta)\) je asymptotická kovariančná matica odhadov parametrov typu \((p + q) \times (p + q)\).
V praxi nie sú skutočné hodnoty parametrov \(\boldsymbol{\beta}\) známe. Preto sa používa aproximácia asymptotickej kovariančnej matice pomocou informačnej matice.
Informačná matica je stredná hodnota matice 2. parciálnych derivácií logaritmu vierohodnostnej funkcie podľa jej parametrov. Inverzná matica k informačnej matici je kovariančná matica odhadov parametrov. Aproximácia spočíva v tom, že namiesto stredných hodnôt sa používajú výberové kvantily. Štandardné chyby odhadov parametrov sú druhé odmocniny diagonálnych prvkov kovariančnej matice vydelené \(\sqrt{n}\).
Posledným bodom postupu je overenie adekvátnosti odhadnutého modelu, či už vyšetrením vlastností modelu alebo vlastností rezíduí, napr. nasledujúcimi metódami:
Dobrý (lineárny) model časových radov určite spĺňa podmienku nekorelovanosti rezíduí (ktoré samozrejme zároveň majú strednú hodnotu nulovú). V opačnom prípade predpovedný model nevyužíva všetky informácie dostupné v časovom rade. Vlastnosť normality a konštantného rozptylu rezíduí nie je nevyhnutná, skôr je užitočná pri konštrukcii intervalových odhadov/predpovedí (niekedy sa dá zabezpečiť vhodnou transformáciou časového radu).
Okrem testovania hodnoty ACF pre každý posun osobitne a to známym testom, ktorým sme už vyšetrovali náhodnosť rezíduí po modeli systematických zložiek, môžeme použiť testy z triedy Portmanteau testov, kde pravdivosť štatistických hypotéz formulovaných \[ H_0\colon \rho(k)=0, \; k = 1, ..., m, \qquad H_1\colon \exists k: 0 < k \leq m, \; \rho(k)\neq 0 \] skúmame na základe výberovej autokorelačnej funkcie rezíduí stacionárneho procesu ARMA(\(p,q\)) pre prvých \(m\) hodnôt spolu (\(m \leq \frac{n}{4}\)), pričom testovacia štatistika, najčastejšie
má pri platnosti \(H_0\) rozdelenie \(\chi^2(m - p - q)\). Test je jednostranný, teda ak \(Q > \chi_{1-\alpha}^2(m - p - q)\), model ARMA(\(p,q\)) nepovažujeme za vhodný (na hladine významnosti \(\alpha\)).
Jarque-Bera test normality rezíduí, teda hypotézy \(H_0\), že rezíduá \(\hat\varepsilon_t\) sú generované normálnym rozdelením, je založený na súčasnom testovaní šikmosti a špicatosti.
Za predpokladu platnosti \(H_0\) a nulovej strednej hodnoty rezíduí sú výberová šikmosť \(SK\) a špicatosť \(SP\) dané vzťahmi a rozdelením (pre dostatočne veľké \(n\)): \[ SK=\frac{\frac{1}{n}\sum_{t=1}^n \hat{\varepsilon}_t^3}{\left[\frac{1}{n}\sum_{t=1}^n \hat{\varepsilon}_t^2 \right]^{\frac{3}{2}}}\sim N\left(0,\frac6n\right), \qquad SP=\frac{\frac{1}{n}\sum_{t=1}^n \hat{\varepsilon}_t^4}{\left[\frac{1}{n}\sum_{t=1}^n \hat{\varepsilon}_t^2 \right]^{2}}\sim N\left(3,\frac{24}n\right) \] Testovacia štatistika Jarque-Bera testu, \[ JB = n \left( \frac{SK^2}{6}+ \frac{(SP-3)^2}{24}\right) \] má potom rozdelenie \(\chi^2(2)\).
Testujeme hypotézu \(H_0\), že časový rad rezíduí \(\hat{\varepsilon}_t\) je rad s podmienenou homoskedasticitou, formálne \[ H_0\colon var(\varepsilon_t|y_{t-1},y_{t-2},...) = \mbox{ konšt.} \] Engleho test je založený na Lagrangeových multiplikátoroch.2 Postup je nasledovný:
Alternatívne sa YW rovnice môžu použiť na odhad ACF pre známy model AR.↩︎
Názov Lagrangeov multiplikátor vyplýva z toho, že test meria, ako veľmi by sa zlepšila vierohodnosť (fit) modelu, ak by sme uvoľnili umelé obmedzenie o konštantnom rozptyle; ide teda o štatistickú analógiu tieňovej ceny (\(\lambda\)) z matematickej optimalizácie.↩︎