4 Predpovedanie
Predpovedaním zvyčajne rozumieme odhad budúcich hodnôt na základe prítomných a minulých hodnôt časového radu. V kontexte regresného modelu, ktorý môže okrem endogénnych hodnôt (teda patriacich tomu istému stochastickému procesu) obsahovať aj hodnoty exogénnych premenných (náhodných aj deterministických), rozumieme pod pojmom predpoveď (angl. prediction) odhad rozdelenia vysvetľovanej premennej (odozvy) \(y_t\) pre dané hodnoty vysvetľujúcich premenných (prediktorov, regresorov) \(x_{1,t}, x_{2,t}, \ldots, x_{k,t}\) 1.
Z hľadiska času konštrukcie rozlišujeme predpovede
- pred udalosťou (ex ante) – hodnoty prediktorov je potrebné určiť (napr. ako predpovede z iných modelov) ak nie sú zrejmé (ako pri deterministických funkciách času, napr. \(x_{1,t}=t\), \(x_{2,t}=cos(\frac{2\pi t}{L})\) a podobne),
- po udalosti (ex post) – vtedy sú hodnoty prediktorov známe, pretože mohli byť priamo pozorované.
Ex ante predpovede teda slúžia na podporu rozhodovania, zatiaľčo ex post na vyhodnotenie predpovedných modelov.
Táto kapitola čerpá najmä z publikácií (Cipra 2008, kap. 3.6.5, 8.2.3 a 10.6) a (Hyndman a Athanasopoulos 2018, kap. 3.4, 3.5 a 8.8).
4.1 Bodové a intervalové
Keďže predpovede sú konštruované ako štatistické odhady, klasifikujeme ich rovnako:
- Bodová predpoveď je odhad podmienenej strednej hodnoty,
\[ \hat y_t=E(y_t|x_{1,t},\ldots,x_{k,t}) = \hat\beta_0 + \hat\beta_1 x_{1,t} + \ldots + \hat\beta_k x_{k,t}=\hat{\boldsymbol{\beta}}'\mathbf{x}_t, \] ktorý je nevychýlený, ak sa model chová rovnako i mimo výberový súbor.2 Chyba bodovej predpovede \(e_t=y_t-\hat y_t\) má smerodajnú odchýlku \[ \hat\sigma_{e_t}=\hat{\sigma}_\varepsilon\sqrt{1+\mathbf{x}_t'(\mathbf{X}'\mathbf{X})^{-1}\mathbf{x}_t}, \] kde \(\sigma_\varepsilon\) je smerodajná odchýlka rezíduí modelu a \(\mathbf{X}\) regresná matica obsahujúca bázové funkcie lineárneho modelu (použitá na jeho odhad). Vzťah sa dá interpretovať tak, že chyba predpovede je tvorená súčtom dvoch nekorelovaných chýb:- odchýlky, ktorou sa \(y_t\) líši od svojej strednej hodnoty (t.j. chyba reziduálnej zložky),
- odchýlky v predpovedi tejto strednej hodnoty (vplyvom neznalosti skutočných regresných parametrov).
- Intervalová predpoveď je interval, v ktorom sa s pravdepodobnosťou \((1-\alpha)\) nachádza skutočná hodnota \(y_t\). Konštruuje sa pomocou kvantilov hypotetického rozdelenia pravdepodobnosti, prípadne kvantilov získaných simuláciou. Za predpokladu normality je tento interval symetrický okolo stredu v bodovej predpovedi, s hranicami \[ \hat{y}_t \pm t_{1-\alpha/2}(n-k-1)\hat{\sigma}_{e_t}, \] kde \(t_{p}(\mathit{df})\) je \(p\)-ty kvantil Studentovho \(t\)-rozdelenia s \(\mathit{df}\) stupňami voľnosti.
Konštrukcia ex ante predpovedí pomocou regresných modelov systematických zložiek je jednoduchá, pretože prediktory sú tvorené deterministickými premennými. To znamená, že časová premenná \(t\) a jej transformácie – spojité (goniometrické funkcie) aj nespojité (indikačné funkcie) – sú známe v ktoromkoľvek okamihu predpovede. Potom vektor bázových funkcií (riadok matice plánu \(\mathbf{X}\) vytvorenej na konštrukciu predpovedí) v čase \(t\), napríklad pre model s lineárnym trendom a jednoduchou sezónnou zložkou s periódou \(L=12\), bude mať tvar \(\mathbf{x}_t = (1, t, \cos\frac{2\pi t}{12}, \sin\frac{2\pi t}{12})'\).
Na druhú stranu, konštrukcia predpovedí pomocou ARMA modelu podmieňuje strednú hodnotu minulými hodnotami stochastického procesu. Tie treba v jednotlivých krokoch odhadovať, čo predpoveď komplikuje. Aj tak je však v rámci Box-Jenkins metodológie výpočet pomerne jednoduchý. Predpokladajme, že časový rad \(\{y_t\}\) je realizáciou stochastického procesu \[ y_t - \varphi_1 y_{t-1} - \ldots - \varphi_p y_{t-p} = \varepsilon_t + \theta_1 \varepsilon_{t-1} + ... + \theta_q \varepsilon_{t-q}, \] potom bodová predpoveď hodnoty \(y_{t+h}\) konštruovaná v čase \(t\) je tzv. \(h\)-kroková predpoveď \(\hat{y}_{t+h|t}=E(y_{t+h}|y_t,y_{t-1},...)\), a hodnotu \(h>0\) sa nazýva horizont predpovede. Prakticky sa vypočíta pomocou vzťahu \[ \hat{y}_{t+h|t}=\hat\varphi_1 y_{t+h-1|t} + \ldots + \hat\varphi_p y_{t+h-p|t} + \hat\theta_1 \varepsilon_{t+h-1|t} + \ldots + \hat\theta_q \varepsilon_{t+h-q|t}, \] pričom, ak hodnoty prediktorov nie sú v čase \(t\) známe, aproximujú sa podmienenou strednou hodnotou, čiže \[ \begin{split} & \left. \begin{split} y_{t+j|t} &= y_{t+j} \\ \varepsilon_{t+j|t} &\approx y_{t+j}-\hat{y}_{t+j|t+j-1} \end{split} \right\} \quad \forall j\leq0 \\ &\left. \begin{split} y_{t+j|t} &\approx \hat{y}_{t+j|t} \\ \varepsilon_{t+j|t} &\approx 0 \end{split} \right\} \quad \forall j\colon\ 0 < j < h \end{split} \]
To znamená, že kým môžeme, dosádzame pozorované hodnoty, resp. rezíduá modelu. Následne sa dajú dosadiť už iba hodnoty predpovedané v predošlých krokoch a chybovú zložku aproximujeme jej strednou hodnotou, teda nulou.
Výpočet intervalovej predpovede z ARMA modelu je náročnejší, no za predpokladu normality a stacionarity (model sa dá rozšíriť ako MA(\(\infty\)) s koeficientami \(\psi_1,\psi_2\ldots\)) je ju možné pri požadovanej úrovni spoľahlivosti (\(1-\alpha\)) aproximovať intervalom \(\hat{y}_{t+h|t} \mp q_{\alpha/2}\hat\sigma_h\), kde \(q_{p}\) je \(p\)-kvantil N(0,1) rozdelenia a rozptyl predpovede \(\hat{\sigma}_h^2=\hat{\sigma}_\varepsilon^2(\sum_{j=0}^{h-1}\hat{\psi}_j^2)\) pričom \(\psi_0=1\). Z toho je vidieť, že výpočet intervalu jednokrokovej predpovede pomocou ARMA modelu alebo viackrokovej predpovede pomocou MA modelu je najjednoduchší. Treba ešte podotknúť, že takýto intervalový odhad je oproti skutočnému užší, pretože počíta iba s variabilitou v chybovej zložke \(\varepsilon_t\), zanedbáva tak neistotu v určení parametrov a stupňov modelu.
Alternatívnym spôsobom výpočtu predpovedných intervalov (odporúčaný napríklad pri porušení predpokladu normality predpovedných chýb) je použitie simulácie možnej trajektórie procesu. Vtedy sa vygeneruje \(N\) trajektórií, \(\{(\tilde{y}_{t+1|t}^{(i)},\ldots,\tilde{y}_{t+h|t}^{(i)})\}_{i=1}^N\), podľa vzťahu \[ \tilde{y}_{t+h|t}^{(i)}=\hat\varphi_1 y_{t+h-1|t}^{(i)} + \ldots + \hat\varphi_p y_{t+h-p|t}^{(i)} + \hat\theta_1 \varepsilon_{t+h-1|t}^{(i)} + \ldots + \hat\theta_q \varepsilon_{t+h-q|t}^{(i)} + \tilde{\varepsilon}_{t+h|t}^{(i)} \]
pričom hodnoty prediktorov \(i\)-tej trajektórie sa v čase \(t\) určia podľa kľúča \[
\begin{split}
& \left.
\begin{split}
y_{t+j|t}^{(i)} &= y_{t+j} \\
\varepsilon_{t+j|t}^{(i)} &= y_{t+j}-\hat{y}_{t+j|t+j-1}
\end{split}
\right\} \quad \forall j\leq0
\\
&\left.
\begin{split}
y_{t+j|t}^{(i)} &= \tilde{y}_{t+j|t}^{(i)} \\
\varepsilon_{t+j|t}^{(i)} &= \tilde{\varepsilon}_{t+j|t}^{(i)}
\end{split}
\right\} \quad \forall j\colon\ 0 < j < h
\end{split}
\] Hodnoty \(\tilde{\varepsilon}_{t+j|t}^{(i)}\) pre \(0<j\leq h\) sa vygenerujú pre každú trajektóriu zvlášť, a to
- buď z predpokladaného rozdelenia (metóda Monte Carlo)
- alebo náhodným výberom z odhadu minulých realizácií \(\{\hat\varepsilon_{t+j}\}_{j<0}\) – čiže z rezíduí modelu (metóda bootstrap).
Potom z veľkého počtu takýchto simulácií v každom okamihu \(t+h\) bude \((1-\alpha)\%\)-ný predpovedný interval zdola ohraničený výberovým \((\alpha/2)\)-kvantilom a zhora \((1-\alpha/2)\)-kvantilom súboru \((\tilde{y}_{t+h|t}^{(1)},\ldots,\tilde{y}_{t+h|t}^{(N)})\). Prirodzene, tento súbor hodnôt trajektórií v danom časovom okamihu sa môže použiť aj na výpočet bodovej predpovede (vo forme strednej hodnoty alebo mediánu).
Pre úplnosť témy spomenieme, že najjednoduchšou metódou, aká sa na predpovedanie používa, je aproximácia budúcich hodnôt stochastického procesu konštantnou úrovňou, či už vo forme
- priemernej hodnoty pozorovaného časového radu, pričom \(\hat\sigma_h=\hat\sigma_{\varepsilon}^2\sqrt{1+1/n}\), alebo
- ako tzv. naivnú predpoveď, tzn. poslednú pozorovanú hodnotu, prípadne poslednú hodnotu pozorovanú v tej istej sezóne. Potom \(\hat\sigma_h=\hat\sigma_{\varepsilon}^2\sqrt{h}\) a v prípade sezónnej naivnej \(\hat\sigma_h=\hat\sigma_{\varepsilon}^2\sqrt{k+1}\), kde \(k=\lfloor(h-1)/L\rfloor\) a \(L\) je perióda sezonality. Pre odvodenie pozri napr. blog.
Pri stacionárnych časových radoch priemerná hodnota vedie k triezvym predpovediam, ktoré presnosťou môžu predčiť mnohé sofistikované (ale nesprávne) modely. V prípade nestacionárnych časových radov má väčšiu šancu na úspech skôr naivná predpoveď, teda posledná pozorovaná hodnota. Obe sa však považujú za elementárne predpovedacie metódy, s ktorými sa porovnávajú iné, vhodnejšie modely.
4.2 Miery presnosti
Modely časových radov sa posudzujú nielen z hľadiska ich popisnej schopnosti, ale aj predpovednej, to znamená, ako dobre popisuje realitu v čase, v ktorom nebol trénovaný. Model, ktorý dáta dobre popisuje, nemusí dobre predpovedať. Problém preurčenia (vytvorenia príliš komplikovaného) modelu je rovnako závažný ako neschopnosť zachytiť signál (systematický vzor) v dátach.
Preto sa na začiatku pred výstavbou modelov časový rad rozdelí na trénovaciu vzorku \(y_1,\ldots,y_n\) a validačnú (testovaciu) vzorku \(y_{n+1},\ldots,y_{n+m}\), ktorú zvyčajne tvorí posledných cca 20% pozorovaní. Zatiaľčo identifikácia a odhad parametrov modelu prebehne za pomoci trénovacej vzorky, validačná vzorka je použitá na výpočet predpovedných chýb \(\hat e_t=y_t-\hat{y}_{t}\), kde \(t=n+1,\ldots,n+m\). V tejto súvislosti predpovede \(\hat{y}_{t}\) rozlišujeme
- jednokrokové \(\{\hat{y}_{t|t-1}\}_{t=n+1}^{n+m}\) – na ich výpočet sa použijú pozorované hodnoty z validačnej vzorky,
- viackrokové \(\{\hat{y}_{n+j|n}\}_{j=1}^m\),
v zásade podľa toho, akým spôsobom (pre aký horizont) chceme model ďalej využívať na predpovede.
Na vyhodnotenie predpovednej schopnosti je potrebné zhrnúť (okamžité) chyby predpovedí \(\hat e_t\) do jednej číselnej charakteristiky – strednej predpovednej chyby. Takých sa používa niekoľko druhov:
- chyby, ktoré závisia na mierke \(y_t\) (scale-dependent errors), hodia sa na porovnanie modelov toho istého časového radu, napr.
- stredná kvadratická chyba (mean squared error) \[
MSE=\frac1m\sum_{t=n+1}^{n+m}\hat e_t^2
\] alebo jej odmocnina \(RMSE=\sqrt{MSE}\), ktoré sú analógiou mier používaných na ohodnotenie popisných schopností modelu (reziduálny rozptyl \(\sigma_\varepsilon^2\) resp. št.odchýlka \(\sigma_\varepsilon\)). Výhoda \(RMSE\) spočíva v rovnakých jednotkách ako časový rad.
- stredná absolútna chyba (mean absolute error) \[
MAE=\frac1m\sum_{t=n+1}^{n+m}|\hat e_t|
\] je oproti \((R)MSE\) miernejšia voči veľkým chybám, ktoré vzniknú predovšetkým prítomnosťou odľahlých hodnôt v časovom rade.
- stredná kvadratická chyba (mean squared error) \[
MSE=\frac1m\sum_{t=n+1}^{n+m}\hat e_t^2
\] alebo jej odmocnina \(RMSE=\sqrt{MSE}\), ktoré sú analógiou mier používaných na ohodnotenie popisných schopností modelu (reziduálny rozptyl \(\sigma_\varepsilon^2\) resp. št.odchýlka \(\sigma_\varepsilon\)). Výhoda \(RMSE\) spočíva v rovnakých jednotkách ako časový rad.
- chyby v relatívnej mierke (percentage errors), ktoré nezávisia na mierke \(y_t\), preto sa hodia na porovnanie predpovednej schopnosti medzi rôznymi časovými radmi, nadobúdajú zvyčajne hodnoty od 0% do 100% a nehodia sa pre časové rady s hodnotami blízkymi či dokonca rovnými nule, alebo takými, v ktorých nula nemá zásadný význam (napr. teplota v \(^\circ\)C a \(^\circ\)F oproti Kelvinom). Najznámejšie sú
- stredná absolútna percentuálna chyba (mean absolute percentage error) \[ MAPE=\frac{100}{m}\sum_{t=n+1}^{n+m}\left|\frac{\hat e_t}{y_t}\right| \]
- symetrická MAPE \[ sMAPE=\frac{100}{m}\sum_{t=n+1}^{n+m}\left|\frac{\hat e_t}{(y_t+\hat{y}_t)/2}\right| \] ktorá dáva rovnaký výsledok aj po zámene skutočných a predpovedaných hodnôt.
- stredná absolútna chyba preškálovaná pomocou MAE naivných predpovedí na trénovacej vzorke (mean absolute scaled error) \[ MASE=\frac1m\sum_{t=n+1}^{n+m}\left| \frac{\hat e_t}{\frac1{n-1}\sum_{t=2}^n|y_t-y_{t-1}|} \right| \] ktorá nadobúda hodnoty menšie ako 1, ak dáva model v priemere lepšie predpovede ako naivná predpoveď počas trénovacieho obdobia. Hodí sa na porovnanie schopnosti predpovednej metódy (triedy modelu) medzi rôznymi časovými radmi.
4.3 Diebold-Mariano test
Porovnanie mier presnosti bodových predpovedí nám nič nepovie o tom, či je takto zostavené poradie modelov platné pre celý základný súbor alebo je len zhodou okolností stojacich za náhodným výberom. Na detekciu významných rozdielov medzi predpoveďami potrebujeme formálny štatistický test, akým je napr. Diebold-Mariano (pozri (Diebold 2015) pre kritický pohľad jedného z autorov). Nulová hypotéza hovorí, že predpovede dvoch modelov sú rovnako presné.
Označme jednokrokovú predpovednú chybu \(i\)-teho modelu ako \(\hat e_{i,t}\) a zvolme tzv. stratovú (angl. loss) funkciu \(g\), najčastejšie \(g(x)=x^2\). Ak platí \(H_0\), rozdiel \(d_t = g(\hat e_{1,t}) - g(\hat e_{2,t})\) bude mať strednú hodnotu \(E(d_t) = 0\) a za predpokladu kovariančnej stacionarity časového radu \(\{d_t\}\) je asymptotické rozdelenie výberového priemeru \(\bar{d}=\frac{1}{m}\sum_{t=n+1}^{n+m} d_t\) dané vzťahom \(\sqrt{m} (\bar{d}-E(d_t)) \sim N(0,\sigma_{\bar{d}}^2)\).
Testovacia štatistika \[ DM=\frac{\bar{d}}{\sqrt{\hat\sigma_{\bar{d}}^2}} \] kde \(\hat\sigma_{\bar{d}}^2=\frac{1}{m}\sum_{t=n+1}^{n+m} \left( d_t-\bar{d}\right)^2\) je konzistentný odhad rozptylu \(\hat\sigma_{\bar{d}}^2\), má za predpokladu platnosti \(H_0\colon E(d_t) = 0\) asymptoticky normálne rozdelenie \(N(0, 1)\). Vlastnosti testovacej štatistiky všeobecne pre \(h\)-krokové predpovede sumarizuje napr. (Franses, Van Dijk, a Opschoor 2014).
Vyhodnotenie predpovednej schopnosti viacerých modelov potom môže byť urobené tak, že výsledky testov medzi dvojicami modelov sa zapisujú do tabuľky (matice) \(A\), konkrétne \[ A_{ij}=\begin{cases} 1 & \text{ak model $i$ predpovedá lepšie ako model $j$} \\ -1 & \text{ak model $i$ predpovedá horšie ako model $j$} \\ 0 & \text{inde} \end{cases} \] potom poradie modelov sa určí zostupným usporiadaním riadkových súčtov tabuľky \(A\).