Lineární a adaptivní zpracování dat 9. Modely časových řad II.
Daniel Schwarz
Investice do rozvoje vzdělávání
Opakování K čemu je dobré vytvářet modely procesů generující časové řady? Dekompozice časový řad: jaké přístupy znáte? Vyjmenujte složky, do kterých časové řady rozkládáme. Které z nich tvoří užitečnou informaci o procesu? Je nějaký rozdíl mezi sezónní a cyklickou složkou? Jaký? Definujte stacionární stochastický proces. Jaké máme možnosti pro očištění časových řad od trendu? Jaké máme možnosti pro očištění časových řad od sezónních vlivů?
Bi0440
© Institute of Biostatistics and Analyses
Modely časových řad Jakoukoli stacionární časovou řadu či signál s náhodnou složkou generuje stochastický proces, kterému lze přiřadit jeden z těchto modelů: • • • •
čistě rekursivní model nerekursivní model s klouzavým průměrem kombinovaný model bílý šum
Bi0440
AR – autoregressive MA – moving average ARMA ν
© Institute of Biostatistics and Analyses
Modely časových řad Jakoukoli stacionární časovou řadu či signál s náhodnou složkou generuje stochastický proces, kterému lze přiřadit jeden z těchto modelů: • • • •
čistě rekursivní model nerekursivní model s klouzavým průměrem kombinovaný model bílý šum
AR – autoregressive MA – moving average ARMA ν
V případě nestacionarity procesu, který časovou řadu generuje, modelujeme nejdříve časovou řadu oproštěnou od trendových a sezónních (cyklických) vlivů a ve výsledném modelu pak tyto vlivy opět uplatníme. Hovoří se pak i o modelech ARIMA a SARIMA. Bi0440
© Institute of Biostatistics and Analyses
Bílý šum Náhodný proces označujeme za bílý šum, pokud jeho střední hodnota a autokorelační funkce (ACF) splňují tyto podmínky: Diracova distribuce
μν = Ε{ν n } = 0 , Rνν (n1 , n2 ) = Ε{ν (n1 )ν (n2 )} =
N0 δ (n1 − n2 ). 2
1000
0.6 0.4
800
0.2
600
Rww(n1,n2)
w(n)
Empirická ACF
0 -0.2 -0.4
400 200 0
-0.6 0
20
40
60 n
Bi0440
80
100
-200 -100
-50
0 50 n1-n2 © Institute of Biostatistics and Analyses
100
Autoregresní (AR) model
xn = a1xn−1 + a2 xn−2 + ...+ ap xn− p +ν n xn…………… časová řada / signál ν n……………bílý šum p …………… řád AR modelu
ai…………… parametry modelu Odhad parametrů
xn
νn
+
z‐1
z‐1 a1
z‐1 a2
xn Bi0440
© Institute of Biostatistics and Analyses
ap
Autoregresní (AR) model
xn = a1xn−1 + a2 xn−2 + ...+ ap xn− p +ν n • AR model je lineární regrese aktuální hodnoty řady proti jedné a více předcházejícím hodnotám. • Aktuální hodnota časové řady je dána lineární kombinací jejich předchozích hodnot. • AR modely mají přímou a jasnou interpretaci. • Rekurzivní systém ‐ buzení IIR filtru bílým šumem. • Stanovení parametrů modelu: metoda nejmenších čtverců neboli minimalizace rozptylu ν. xn
νn
+
z‐1
z‐1 a1
Bi0440
z‐1 a2
© Institute of Biostatistics and Analyses
ap
(MA) model s klouzavým průměrem
xn =ν n + c1ν n−1 + c2ν n−2 + ...+ cqν n−q xn…………… časová řada / signál ν n……………bílý šum q …………… řád MA modelu μ …………… střední hodnota náhodného procesu
ci…………… parametry modelu νn
355
koncentrace CO2
350
Odhad parametrů
345 340 335 330 325 1974
Bi0440
1976
1978
1980 1982 čas
1984
z‐1
z‐1 c1
z‐1 c2
z‐1 cq‐1
xn 1986
cq +
1988
© Institute of Biostatistics and Analyses
xn
(MA) model s klouzavým průměrem
xn =ν n + c1ν n−1 + c2ν n−2 + ...+ cqν n−q • MA model je lineární regrese aktuální hodnoty řady ovšem tentokrát proti vzorkům bílého šumu. Komplikace: v naměřených datech obvykle tyto šumové hodnoty nejsou k dispozici • MA modely mají horší interpretaci než AR modely. • Nerekurzivní systém ‐ buzení FIR filtru bílým šumem. • Stanovení parametrů modelu: nelineární iterační techniky.
νn
z‐1
z‐1 c1
z‐1 c2
z‐1 cq‐1
cq +
Bi0440
© Institute of Biostatistics and Analyses
xn
ARMA model ARMA(p, q) kombinuje AR(p) a MA(q) modely.
xn = a1xn−1 + a2 xn−2 + ...+ ap xn− p +ν n − − c1ν n−1 − c2ν n−2 − ...− cqν n−q Boxova‐Jenkinsova metodologie zahrnuje: • identifikaci modelu, • odhad modelu, • validaci modelu.
Určení řádů p, q Výpočet parametrů ai, ci
Kontrola rozložení residuí Bi0440
© Institute of Biostatistics and Analyses
ARMA model: identifikace 1. 2.
Je časová řada / signál stacionární? Vykazuje časová řada / signál sezónnost? ANO
NE
zjištění periody T, zahrnutí členu AR(T) nebo MA(T) do modelu, případně sezónní diference
ARIMA (autoregressive integrated moving average model) Identifikace, odhad, validace stacionárního ARMA modelu na diferencovaných datech a následná úprava modelu na nestacionární model ARIMA. Př. model AR(2): yn = -0.406yn-1-0.146yn-2-0.00521 = xn-xn-1 xn-xn-1=-0406(xn-1-xn-2)-0.164(xn-2-xn-3)-0.00521 xn=0.594xn-1+0.242xn-2+0.164xn-3-0.00521 Bi0440
© Institute of Biostatistics and Analyses
ARMA model: identifikace Určení řádů p a q a) na základě zkušenosti a experimentování b) spektrum: každé výrazné maximum v rozsahu <0,fvz/2> vyžaduje jeden pár pólů, což zvyšuje řád o 2. c) kritéria na základě autokorelační funkce (ACF) a parciální autokorelační funkce (PACF)
Srovnávání teoretických průběhů ACF, PACF procesů známých řádů s ACF, PACF naměřených časových řad
Bi0440
© Institute of Biostatistics and Analyses
Parciální autokorelační funkce - PACF
PACF představuje korelaci mezi rezidui regresí dvou proměnných vůči třetí proměnné.
Pro zpoždění k vyjadřuje PACF vztah mezi dvěma zpožděnými vzorky časové řady xn a xn+k tak, že nezapočítává lineární vliv vzorků ležících mezi nimi. PACF(k): jedná se o autokorelační koeficient pro zpoždění k očištěný od vlivu xn+1, xn+2,…, xn+k‐1
Bi0440
© Institute of Biostatistics and Analyses
ARMA model: identifikace xn
Příklad: AR(1) proces
νn
xn + axn −1 = ν n
+
z‐1 ‐a1
xn = ν n − axn −1 =
= ν n − a(ν n −1 − axn − 2 ) =
= ν n − aν n −1 + (− a ) ν n − 2 + ... + (− a ) 2
Bi0440
n −1
(ν n−1 − ax0 )
© Institute of Biostatistics and Analyses
ARMA model: identifikace Příklad: AR(1) proces
xn + axn −1 = ν n xn = ν n − axn −1 =
= ν n − a(ν n −1 − axn − 2 ) =
= ν n − aν n −1 + (− a ) ν n − 2 + ... + (− a )
n −1
2
(ν n−1 − ax0 )
|a|<1:
1 − (− a ) E{xn } = μ = μν , 1+ a 1 ≈ μν 1+ a n
Bi0440
{ }
D{xn } = E xn2 − μ 2 =
{
}
σ ν2 1− a
,
2
Rxx ( k ) = E xn xn − k = (− a ) . k
© Institute of Biostatistics and Analyses
ARMA model: identifikace Příklad: AR(1) proces
xn + axn −1 = ν n
Bi0440
{
}
Rxx ( k ) = E xn xn − k = (− a ) . k
© Institute of Biostatistics and Analyses
ARMA model: identifikace Příklad: AR(1) proces
xn + axn −1 = ν n
{
}
Rxx ( k ) = E xn xn − k = (− a ) . k
V literatuře lze nalézt teoretické vlastnosti základních stacionárních procesů AR(1), AR(2), MA(1), MA(2) atd. Tyto teoretické vlastnosti (např. průběhy teoretické autokorelační funkce) se pak srovnávají s empirickými vlastnostmi získanými z naměřených časových řad a výsledky porovnání slouží pro identifikaci p, q.
Bi0440
© Institute of Biostatistics and Analyses
ARMA model: identifikace
Tvar ACF Exponenciála klesající k nule.
Model
Změny kladných a záporných hodnot, postupný pokles k nule.
AR(p) model. Pro určení p se vychází z PACF, která je pro AR(p) nulová od zpoždění p+1.
Jeden nebo několik vrcholů, zbytek zanedbatelný, nulový.
MA model. Řád odpovídá hodnotě zpoždění, od kterého je ACF nulová.
Průběh klesající až po několika zpožděních
Kombinovaný model ARMA.
Vše zanedbatelné, nulové
Data jsou náhodná, generuje je bílý šum.
Vysoké hodnoty ve stejných intervalech
Zahrnout AR člen s řádem odpovídajícím periodě.
Neklesá k nule
Nejedná se o stacionární řadu / signál.
Bi0440
© Institute of Biostatistics and Analyses
ARMA model: identifikace
očištění od trendu
Detekce sezónních komponent
Perioda 12 vzorků
Bi0440
© Institute of Biostatistics and Analyses
ARMA model: identifikace Zbývá identifikovat ještě nesezónní komponenty signálu / řady.
Sezónní diferencování:
Bi0440
© Institute of Biostatistics and Analyses
ARMA model: identifikace Zbývá identifikovat ještě nesezónní komponenty signálu / řady.
Bi0440
© Institute of Biostatistics and Analyses
ARMA model: identifikace Zbývá identifikovat ještě nesezónní komponenty signálu / řady. Hledáme tzv. bod useknutí, tj. zpoždění, od kterého je PACF nulová/zanedbatelná.
Bi0440
© Institute of Biostatistics and Analyses
ARMA model ARMA(p, q) kombinuje AR(p) a MA(q) modely.
xn = a1xn−1 + a2 xn−2 + ...+ ap xn− p +ν n − − c1ν n−1 − c2ν n−2 − ...− cqν n−q Odhad parametrů modelu ‐ iterační algoritmy: ‐ nelineární metoda nejmenších čtverců ‐ odhad na základě maximální věrohodnosti (MLE) Vhodnější tvar rovnice (pro SW nástroje):
xn + a1xn−1 + a2 xn−2 + ...+ ap xn− p =ν n + c1ν n−1 + c2ν n−2 + ...+ cqν n−q
Bi0440
© Institute of Biostatistics and Analyses
ARMA model ARMA(p, q) kombinuje AR(p) a MA(q) modely.
xn = a1xn−1 + a2 xn−2 + ...+ ap xn− p +ν n − − c1ν n−1 − c2ν n−2 − ...− cqν n−q Odhad parametrů modelu ‐ iterační algoritmy: ‐ nelineární metoda nejmenších čtverců ‐ odhad na základě maximální věrohodnosti (MLE) Výsledný AR(2) model (pro data očištěná od lineárního trendu a sezónních vlivů)
xn +1.745⋅ xn−1 + 0.8745⋅ xn−2 =ν n Bi0440
© Institute of Biostatistics and Analyses
ARMA modely Validace modelu
xn +1.745⋅ xn−1 + 0.8745⋅ xn−2 =ν n
‐Zpětné ověření předpokladů kladených na náhodné chyby, tj. analýza residuí RESIDUA = CHYBY PREDIKCE ‐ Residua by měla představovat bílý šum.
Bi0440
© Institute of Biostatistics and Analyses
ARMA modely Validace modelu AR(2)
xn +1.745⋅ xn−1 + 0.8745⋅ xn−2 =ν n
Validace modelu AR(4)
xn + 2.943⋅ xn−1 + 3.765⋅ xn−2 + 2.494⋅ xn−3 + 0.728⋅ xn−4 =ν n
Bi0440
© Institute of Biostatistics and Analyses
ARMA modely Validace modelu AR(2)
xn +1.745⋅ xn−1 + 0.8745⋅ xn−2 =ν n
Validace modelu ARMA(4,4)
xn + 2.345⋅ xn−1 + 2.581⋅ xn−2 +1.645⋅ xn−3 + 0.5163⋅ xn−4 = =ν n − 3.789⋅ν n−1 + 5.397⋅ν n−2 − 3.425⋅ν n−3 + 0.8168⋅ν n−4
Bi0440
© Institute of Biostatistics and Analyses
ARMA modely Validace modelu AR(2)
xn +1.745⋅ xn−1 + 0.8745⋅ xn−2 =ν n
Ani postupným navyšováním složitosti (řádu) modelu se nedospělo k uspokojivým výsledkům Validace modelu ARMA(4,4) (ACF residuí neodpovídá ACF bílého šumu).
xn + 2.345⋅ xn−1 + 2.581⋅ xn−2 +1.645⋅ xn−3 + 0.5163⋅ xn−4 = V takovém případě je nejlepší návrat k co =ν n − 3.789⋅ν n−1 + 5.397⋅ν n−2 − 3.425nejjednoduššímu modelu: ⋅ν n−3 + 0.8168⋅ν n−4 v ukázaném příkladu se jedná o AR(2).
Bi0440
© Institute of Biostatistics and Analyses
ARMA modely Posouzení kvality předpovídání ‐ aplikace modelu na řadu zkrácenou o m pozorování, ‐ předpověď hodnot xn‐m+1, xn‐m+2,…,xn ‐ porovnání
xn +1.745⋅ xn−1 + 0.8745⋅ xn−2 =ν n
Model: AR(2) Horizont predikce: 1 Shoda: 82.1 %
Bi0440
© Institute of Biostatistics and Analyses
ARMA modely Posouzení kvality předpovídání ‐ aplikace modelu na řadu zkrácenou o m pozorování, ‐ předpověď hodnot xn‐m+1, xn‐m+2,…,xn ‐ porovnání
xn +1.745⋅ xn−1 + 0.8745⋅ xn−2 =ν n
Model: AR(2) Horizont predikce: 5 Shoda: < 1 %
Bi0440
© Institute of Biostatistics and Analyses
ARMA modely Posouzení kvality předpovídání ‐ aplikace modelu na řadu zkrácenou o m pozorování, ‐ předpověď hodnot xn‐m+1, xn‐m+2,…,xn ‐ porovnání
xn + 2.943⋅ xn−1 + 3.765⋅ xn−2 + 2.494⋅ xn−3 + 0.728⋅ xn−4 =ν n
Model: AR(4) Horizont predikce: 1 Shoda: 92.0 %
Bi0440
© Institute of Biostatistics and Analyses
ARMA modely Posouzení kvality předpovídání ‐ aplikace modelu na řadu zkrácenou o m pozorování, ‐ předpověď hodnot xn‐m+1, xn‐m+2,…,xn ‐ porovnání
xn + 2.943⋅ xn−1 + 3.765⋅ xn−2 + 2.494⋅ xn−3 + 0.728⋅ xn−4 =ν n
Model: AR(4) Horizont predikce: 5 Shoda: 3.3 %
Bi0440
© Institute of Biostatistics and Analyses
ARMA modely Posouzení kvality předpovídání ‐ aplikace modelu na řadu zkrácenou o m pozorování, ‐ předpověď hodnot xn‐m+1, xn‐m+2,…,xn ‐ porovnání x n + 2 . 345 ⋅ x n − 1 + 2 . 581 ⋅ x n − 2 + 1 . 645 ⋅ x n − 3 + 0 . 5163 ⋅ x n − 4 = =ν
n
− 3 . 789 ⋅ ν
n −1
+ 5 . 397 ⋅ ν
n−2
− 3 . 425 ⋅ ν
n−3
+ 0 . 8168 ⋅ ν
n−4
Model: ARMA(4,4) Horizont predikce: 1 Shoda: 98.6 %
Bi0440
© Institute of Biostatistics and Analyses
ARMA modely Posouzení kvality předpovídání ‐ aplikace modelu na řadu zkrácenou o m pozorování, ‐ předpověď hodnot xn‐m+1, xn‐m+2,…,xn ‐ porovnání x n + 2 . 345 ⋅ x n − 1 + 2 . 581 ⋅ x n − 2 + 1 . 645 ⋅ x n − 3 + 0 . 5163 ⋅ x n − 4 = =ν
n
− 3 . 789 ⋅ ν
n −1
+ 5 . 397 ⋅ ν
n−2
− 3 . 425 ⋅ ν
n−3
+ 0 . 8168 ⋅ ν
n−4
Model: ARMA(4,4) Horizont predikce: 5 Shoda: 23.9 %
Bi0440
© Institute of Biostatistics and Analyses
LTI systém a jeho popis yn = xn −1.745⋅ yn−1 − 0.8745⋅ yn−2 yn xn
+
z‐1
z‐1 ‐1.7
Bi0440
‐0.9
© Institute of Biostatistics and Analyses
LTI systém a jeho popis yn = xn −1.745⋅ yn−1 − 0.8745⋅ yn−2 yn xn
+
z‐1
z‐1 ‐1.7
‐0.9
Nesystematickou složku předložené časové řady (měření koncentrace CO2 v ovzduší) očištěné od trendu a sezónních vlivů modelujeme jako buzení IIR filtru druhého řádu.
Bi0440
© Institute of Biostatistics and Analyses
Shrnutí Popis a identifikace systémů a procesů
z‐1
z‐1 c1
z‐1 c2
z‐1 cq‐1
cq +
Analýza, Simulace, Predikce, Monitoring, Diagnostika, Řízení
Bi0440
© Institute of Biostatistics and Analyses
9. cvičení 1) Z předložených dat co2.csv představujících monitoring koncentrace CO2 v ovzduší identifikujte a odhadněte model procesu, který tato data vygeneroval. Využijte aditivní přístup k dekompozici časových řad a dále Box‐Jenkisovu metodiku pro identifikaci struktury modelu. Pro odhad parametrů modelu využijte funkce z System Identification Toolbox Matlabu. Nezapomeňte svůj model validovat pomocí analýzy residuí a po té se pokuste předpovědět koncentraci ovzduší v následujících letech.
Bi0440
© Institute of Biostatistics and Analyses
9. cvičení
c= 1.4519*time-2537.2
Bi0440
© Institute of Biostatistics and Analyses
9. cvičení
Bi0440
© Institute of Biostatistics and Analyses
9. cvičení
Bi0440
© Institute of Biostatistics and Analyses
9. cvičení
Bi0440
© Institute of Biostatistics and Analyses
9. cvičení
Bi0440
© Institute of Biostatistics and Analyses
9. cvičení
Bi0440
© Institute of Biostatistics and Analyses
9. cvičení xdd n = 3.74 xdd n −1+6.721xdd n − 2 + 7.47 xdd n −3 + 5.429 xd n − 4 + 2.44 xdd n −5 + 0.5302 xdd n −6 +ν n xd n − xd n −12 = 3.74( xd n −1 − xd n −13 ) + 6.721(xd n − 2 − xd n −14 ) + 7.47( xd n −3 − xd n −15 ) + 5.429( xd n − 4 − xd n −16 ) + 2.44( xd n −5 + xd n −17 ) + 0.5302( xd n −6 − xd n −18 ) + ν n xd n = 3.74 xd n −1+6.721xd n − 2 + 7.47 xd n −3 + 5.429 xd n − 4 + 2.44 xd n −5 + 0.5302 xd n −6 + + xd n −12 − 3.74 xd n −13 − 6.721xd n −14 + 7.47 xd n −15 − 5.429 xd n −16 − 2.44 xd n −17 − 0.5302 xd n −18 + ν n xn = 1.4519n − 2537.2 + 3.74 x n −1 +6.721xn − 2 + 7.47 xn −3 + 5.429 xn − 4 + 2.44 xn −5 + 0.5302 xn −6 + + xn −12 − 3.74 xn −13 − 6.721xn −14 − 7.47 xn −15 − 5.429 xn −16 − 2.44 xn −17 − 0.5302 xn −18 + ν n
Bi0440
© Institute of Biostatistics and Analyses
9. cvičení
Bi0440
© Institute of Biostatistics and Analyses
9. cvičení
Bi0440
© Institute of Biostatistics and Analyses
9. cvičení
Bi0440
© Institute of Biostatistics and Analyses
9. cvičení
Bi0440
© Institute of Biostatistics and Analyses
;
ffgf
Otázky ?
[email protected]
50 Bi0440
© Institute of Biostatistics and Analyses