Použité package

depmixS4 quantmod

Cieľ

Predikcia ceny kukurice na ďalší ďeň.

Data

Denné ceny futurov kukurice 2005-2016

Z nich som zobral Close cenu a spravil 2 features: returns tj. \(r=\frac{y_{i}-y_{i-1}}{y_{i-1}}\) a 5 dňový moving avarage z returns.

Data som rozdelil na trenovaciu a validizačnú množinu. 70% trenovaciu a 30% validizačnú.

HMM model

Pomocou depmixS4 som postavil 11 stavový HMM model s pozorovaniami, ktorých emission probalitities sa riadia 2-rozmenrým normálnym rozdelením. Následne bol model natrénovaný na trénovacích dátach pomocou EM algoritmu.

Viterbiho cesta na základe fitovaného modelu.

Clustering

Následne som rozdelil trenovaciu množinu “returns” na 11 časti a podľa odhadov posteriornych pravdepodobnosti v akom stave sa systém nachádza.

Zhluknuté clustre dát podľa stavu v ktorom stave sa nachádza od 1 až po 11

Trenovanie ESN

Na každej z 11 clustrov dát som na trénoval rekurentnú neurónovú sieť. Presnejšie echo state network (ESN). tj. na konci trénovania som mal 11 špecialistov, každého pre jeden druh dát. Každý z clustrov bol ďalej rozdelení na trénovaciu a validizačnú množinu pre optimalizáciu parametrov ESN.

Parametre ESN

Kvôly veľkosti jednotlivých clustrov som sa obmedzil na rezervoir veľkosti 50 neurónov. Všeobecné pravidlo je čím viac tým lepšie ale má to obmedzenie veľkosti trénovacej množiny neurónovej siete.

Ďalšími parametrami sú leaking rate a regularizacia rezervoirovej matice na vhodný spektrálny polomer. To som robil pomocou funkcie optim() a metódou L-BFGS-B (A limited memory algorithm for bound constrained optimization) kde som minimalizoval štvorec chýb predikcií na valizačnej množine daného clustru.

Na záver som ešte nastavil regularizačný parameter pri ridge regresii, ktorý sa používa pri trenovaní výstupnej matice ESN.

Predikcie

Na základe natrénovanej HMM som robil predikcie pravdepodobnosti všetkých stavov na ďalší deň pomocou forward algoritmu tj. \(P(X_{t+1} \mid X_{t} , Y_{1:t}) = \sum_{X} P(X_{t+1} \mid X_{t}) P(X_{t} \mid Y_{1:t})\) kde \(Y_{1:t})\) sú dáta returns a 5 dňový SMA returns z predchádzajúcich 20 dní.

Súčasne s tým bežalo 11 ESN expertov ktorí robili každý predikcie hodnoty returns na ďalší ďeň. Výsledná hodnota predikcie \(returns_{t+1}\) bola \(returns_{t+1}=\sum_{i} w_{i} y_{i}\) kde \(w_{i}=P(X_{t+1}=i \mid X_{t} , Y_{t-20:t})\) a \(y_{i}\) je predikcia i-teho ESN experta.

Obrázok ukazuje porovnanie mojich predikcií a skutočnej hodnoty returns. Červená je skutočná, čierna sú moje predikcie.

Ako vizuálne vidno nie je to moc dobré.

Evaluácia

Ako benchmarkový model som použil tzv. Wienerov proces. Čo prakticky znamená, že predikcie returns vyzerajú nasledovne: \(returns_{t+1} \thicksim \mathcal{N}(\mu \frac{1}{dt} , \sigma^{2} \frac{1}{dt})\) kde \(\mu\) je výberový priemer z returns za posledných 20 dní, \(\sigma\) je výberová disperzia za posledných 20 dní a \(dt = 20\).

Následne som zrátal MSE pre moj HMM model a pre Wienerov proces.

# MSE pre HMM 
mse
## [1] 0.0002855381
# MSE pre Wienerov proces 
mse2
## [1] 0.0002685905

Tiež som chcel nejak štatisticky otestovať oba modely. Zrátal som si štvorce chýb oboch modelov a chcel som testovať hypotézu H0: mse (HMM) - MSE2(Wiener)=0. Bohužiaľ na párový t-test neboli spolnene predpoklady tj. štvorce chýb oboch modelov pochádzajú z normálneho rozdelenia. To som otestoval pomocou Kolmogorovho-Smirnovho testu.

ks.test(seHMM,pnorm)
## 
##  One-sample Kolmogorov-Smirnov test
## 
## data:  seHMM
## D = 0.5, p-value < 2.2e-16
## alternative hypothesis: two-sided
ks.test(seNaive,pnorm)
## 
##  One-sample Kolmogorov-Smirnov test
## 
## data:  seNaive
## D = 0.5, p-value < 2.2e-16
## alternative hypothesis: two-sided

Tak som použil neparametricku alternativu Wilcoxonov párový test

wilcox.test(seHMM, seNaive, paired = TRUE,conf.int = TRUE)
## 
##  Wilcoxon signed rank test with continuity correction
## 
## data:  seHMM and seNaive
## V = 202220, p-value = 0.0001605
## alternative hypothesis: true location shift is not equal to 0
## 95 percent confidence interval:
##  8.359904e-06 3.043495e-05
## sample estimates:
## (pseudo)median 
##  -5.816121e-06

Ktorý hypotézu H0 zamieta.

Pre zaujímavosť párový t-test dáva rovnaký výsledok

t.test(seHMM, seNaive, paired = TRUE,conf.int = TRUE)
## 
##  Paired t-test
## 
## data:  seHMM and seNaive
## t = 3.7842, df = 837, p-value = 0.0001652
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  8.157212e-06 2.573794e-05
## sample estimates:
## mean of the differences 
##            1.694758e-05

Záver

Oba párové testy ukazujú ,že chyby môjho HMM modelu sú štatisticky horšie ako Wienerov proces (vidno aj z intervalov spoľahlivosti), čo je v podstate predikcia na základe pseudonáhodného generátora z normálneho rozdelenia.

Toto som robil pre HMM modely s 2 až 12 stavmi. HMM s 11 stavmi vychádzala najlepšie ale stále nedostatočne.

Plus ako prvý nástrel bola 2 stavová HMM kde som nerobil lineárnu kombinaciu expertov ale som iba medzi nimi prepinal podľa odhadu stavu v ktorom sa HMM bude nachádzať v t+1