Klasifikace, identifikace a statistická analýza nestacionárních náhodných procesů

Proceedings of International Scientific Conference of FME Session 4: Automation Control and Applied Informatics

Paper 26

Klasifikace, identifikace a statistická analýza nestacionárních náhodných procesů MORÁVKA, Jan1 1

Ing., Ph.D., TŘINECKÉ ŽELEZÁRNY inženýring, a.s., Středisko elektro, MaR a ASŘ, 739 70 Třinec, [email protected], http://www.trz.cz

Abstrakt: V ekonomické i technologické praxi se často setkáváme s nestacionárními náhodnými procesy. Pro korektní analýzu dat těchto procesů je třeba znát moderní výsledky a metody, které umožňují správnou klasifikaci, identifikaci a statistickou analýzu. V opačném případě dochází k různým problémům a chybám, mezi něž zejména patří jev zdánlivé korelace. Mezi moderní výsledky patří klasifikace procesů na typy DS/TS, testování tzv. jednotkových kořenů, kointegrace a kotrendová analýza nestacionárních časových řad. Klíčová slova: nestacionární náhodné procesy typu DS/TS, zdánlivá korelace, klasifikace, identifikace, statistická analýza

1 Úvod Nestacionární náhodné procesy (časové řady) se poměrně hojně vyskytují jak v ekonomických systémech [CIPRA, T. 1986], [ARLT, J. 1999], tak i v systémech technologických. V technologii tyto procesy vznikají hlavně jako výstupy integračních soustav, tj. zásobníků, nádrží, homogenizačních hromad, skládek apod. [MORÁVKA, J. 1999]. Pro korektní statistickou analýzu nestacionárních náhodných procesů (klasifikace, identifikace, zjišťování závislostí, regresní analýza) je třeba znát moderní metody [ARLT, J. 1999]. Regresní analýza s nestacionárními veličinami naráží na známý jev (problém) tzv. zdánlivé korelace [SWOBODA, H. 1977], [SEGER, J., HINDLS, R. & HRONOVÁ, S. 1998]. Pro systémy obsahující nestacionární procesy byl pro odlišení zdánlivých vazeb zaveden pojem kointegrace [JOHANSEN, S. 1997], [BIERENS, H. 1999a], [ARLT, J. 1999]. Pro definování nelineárních vazeb mezi nestacionárními veličinami byla rozpracována tzv. kotrendová analýza [BIERENS, H. 1999b].

2 Klasifikace nestacionárních signálů V literatuře jsou nestacionární procesy definovány jako procesy u kterých se s časem mění jejich charakteristiky 1. a 2.řádu, tj. střední hodnota, rozptyl, autokovariance i autokorelace (a tím i spektrum). V teorii i praxi se nejčastěji uvažují dvě skupiny nestacionárních procesů s proměnlivou střední hodnotou a s proměnlivým rozptylem, autokovariancí a autokorelací. Podle terminologie Box-Jenkinsovy metodologie jsou nestacionární procesy označovány jako ARIMA – AutoRegressive Integrated Moving Average [CIPRA, T. 1986], [LJUNG, L. 1987], [TŮMA, J. 1998], [ARLT, J. 1999]. Z hlediska hlubší statistické analýzy je však nutné nestacionární procesy rozdělit na procesy DS – diferenčně stacionární a TS – trendově stacionární [ARLT, J. 1999]. V tab.1 jsou přehledně uvedeny jejich vlastnosti, označení a rozdíly. DS proces lze stacionarizovat diferencováním, TS proces očištěním od deterministického trendu. V běžné literatuře přitom není definováno kritérium volby způsobu stacionarizace [LJUNG, L. 1987], [TŮMA, J. 1998], což nakonec vede na nejednoznačnou a heuristickou volbu.

Tab. 1. Základní klasifikace nestacionárních procesů Integrovaný Nestacionarita Proces Box-Jenkins proces I(d) ARIMA(p,d,q) v rozptylu, AR(1) : ϕ 1 = 1 autokovarianci a DS I(1) autokorelaci ARIMA(0,1,0) ve střední TS I(1)+MA(1) ARIMA(0,1,1) hodnotě

Jiné označení

Běžné označení

náhodná procházka

drift trend

U procesů ARIMA(p,d,q) ~ I(d), kde d znamená řád diferencování procesu. Odhad tohoto řádu lze získat buď testováním postupně diferencované řady na bílý šum pomocí ACF/PACF [ARLT, J. 1999], nebo hledáním minimálního rozptylu této řady [CIPRA, T. 1986]. Náhodná procházka je mezním případem procesu AR(1) s koeficientem autoregrese ϕ 1 = 1. Pro nulové počáteční podmínky v čase t = 0 je také označována jako Brownův, či Wienerův proces [TŮMA, J. 1998]. Integrované procesy typu I(0) tvoří procesy AR(p), MA(q), ARMA(p,q) a proces bílého šumu (BŠ), tj. AR(0) = ARMA(0,0). Modely (diferenční rovnice) nestacionárních procesů DS / TS včetně diferenčních rovnic jejich prvních zpětných diferencí a rovnic procesů s odstraněným lineárním trendem (Tt) mají následující tvar: • DS procesy ! náhodná procházka s konstantou ~ I(1) (1) X t = c + X t −1 + ε t , X t = 0 = X 0 , (2) X t = c + ϕ 1 X t −1 + ε t , X t = 0 = X 0 ~ AR(1) : ϕ 1 = 1 , t (3) X t = X 0 + ct + ∑ ε i , i =1

∆X t = X t − X t −1 = c + ε t

~ bílý šum N (c, σ ε ) , t

X t − Tt = X t − ( X 0 + ct ) = ∑ ε i

~ tzv. stochastický trend .

(4) (5)

i =1

! proces ARIMA(p,d,q) ~ I(d) X t = δ + X t −1 + ψ ( B)ε t , X t = 0 = X 0 , t

X t = X 0 + δt + ψ ( B ) ∑ ε i ,

(6) (7)

i =1

∆X t = δ + ψ ( B)ε t

~ ARMA( p, q) , t

X t − Tt = X t − ( X 0 + δt ) = ψ ( B)∑ ε i ~ I (d ) .

(8) (9)

i =1

• TS procesy TS proces s lineárním trendem ~ I(1) + MA(1) X t = a + bt + ε t , X t = 0 = X 0 = a , !

∆X t = X t − X t −1 = b + ε t − ε t −1 = b + ∆ε t ~ MA(1) ⇒ X t = b + X t −1 + ε t − ε t −1 ~ I (1) + MA(1) / ARIMA(0,1,1) , X t − Tt = X t − (a + bt ) = ε t ~ bílý šum N (0, σ ε ) .

(10) (11) (12) (13)

3 Identifikace nestacionárních signálů Identifikaci nestacionárních signálů je možné provést jak pomocí klasické analýzy časových řad, tak jejími moderními metodami prostřednictvím testů tzv. jednotkových kořenů. 3.1 Identifikace pomocí klasické analýzy časových řad Ze srovnání rovnic diferencí procesů a rovnic procesů s odstraněným lineárním trendem je zřejmé, že pro identifikaci jejich typů lze použít Box-Jenkinsovu metodologii (ACF, PACF a identifikační body [CIPRA, T. 1986]). Pro názornost byly v programu Matlab vygenerovány procesy xds (typu DS / I(1) – náhodná procházka) a xts (typu TS) s parametry: n=100, seed=12345, σ ε =1, pro xds platí: xds0=5, c=0.5 a pro xts je a=5, b=0.6 – obr. 1.

xds

60

xts

50 40 30 20 10 0

10

20

30

40

50

60

70

80

90

100

Obr. 1. Průběhy nestacionárních procesů xds(t) a xts(t) (čárkovaně) Jenom vizuálním posouzením zřejmě nebude možné rozlišit oba typy procesů. Po provedení Box-Jenkinsovy identifikace původních hodnot, jejich diferencí (funkce diff) a trendově očištěných hodnot procesů (funkce detrend) dostaneme výsledky (tab.2-Statgraphics) a algoritmus (obr.2) korespondující s dříve odvozenými vztahy. Tab. 2. Základní klasifikace generovaných nestacionárních procesů xds(t) a xts(t) Procesy Hodnoty Model Závěr xds AR(1) / I(1) procesy nelze rozlišit! původní xts AR(1) / I(1) bílý šum (BŠ) diff xds diferencované xds je DS / I(1) proces diff xts MA(1) AR(1) / I(1) očištěné od lineárního detrend xds xts je TS proces trendu (Tt) bílý šum (BŠ) detrend xts

+ DS / I(1)

TS

TS

Z

Z

∆Xt = Xt - Xt-1

eXt = Xt - T t

ACF, PACF

ACF, PACF

∆Xt = BŠ ?

-

MA(1)

AR(1) / I(1)

TS

eXt = BŠ ?

+

TS

K

K

Obr. 2. Algoritmy určení typů procesů pomocí analýzy časových řad Pozn.: Rozlišení modelů AR(1) a I(1) umožňují až tzv. testy jednotkových kořenů (testují hypotézu H0: ϕ 1 = 1). Závěr: Na základě výsledků a uvedených vztahů lze tedy jednoznačně uskutečnit základní identifikaci nestacionárních procesů na procesy DS a TS včetně jejich zařazení do podskupin integrovaných procesů I(0), I(1) a I(d). 3.2 Identifikace pomocí testů tzv. jednotkových kořenů Identifikaci a rozdělení nestacionárních procesů lze provést pomocí testů tzv. jednotkových kořenů (unit-roots) [ARLT, J. 1999]. Testy vycházejí z diferenční rovnice modelu AR(1) a testují H0, zda autoregresní koeficient 1. řádu ϕ1 = 1. V případě nezamítnutí této nulové hypotézy jde o proces typu DS / I(1), při jejím zamítnutí je proces TS – viz obr. 3. TS Z Testy jedn. kořenů + DS

H0 ?

-

TS

K

Obr. 3. Algoritmus určení procesů DS/TS pomocí testů jednotkových kořenů 3.2.1 Nestacionární procesy typu DS / I(1) Tyto procesy jsou nejčastěji chápány jako integrované procesy 1. řádu I(1) a jsou popsány modelem, či procesem tzv. náhodné procházky (random walk process) s konstantou: (14) X t = c + X t −1 + ε t , X t = 0 = X 0 . Jako příklad pro testování jednotkových kořenů použijeme 2 procesy x1(t) a x2(t) – viz obr. 4, které byly vygenerovány v programu Matlab s parametry: n=100, seed=12345, σ ε =1, x10=x20= 5, c=0 pro x1(t) a c=0.5 pro x2(t).

*+,-./0123425678&96:41/+,;6<&8-9=6<"8-9 )! $! (! #! '! "! &! !

!

"!

#!

$!

%!

&!!

Obr. 4. Průběhy nestacionárních I(1)/DS procesů x1(t) a x2(t) (horní graf) Z obrázku lze vidět, že proces x2(t) má na pohled spíše charakter TS procesu, zatímco proces x1(t) má téměř stacionární průběh. Objektivní charakterizaci procesů lze provést pouze pomocí testů tzv. jednotkových kořenů a to testů nestacionarity i stacionarity – tab. 3. Tab. 3. Výsledky testů nestacionarity a stacionarity DS procesů (EasyReg) Testy nestacionarity Testy stacionarity Proces ADF BG KPSS PP x1(t) -1.956 > -3.45 -11.946 > -21.5 3.829 < 12.71 0.135 < 0.146 x2(t) -1.956 > -3.45 -11.946 > -21.5 3.829 < 12.71 0.135 < 0.146 Pozn.: ADF ... Rozšířený (Augmented) Dickey-Fullerův test, PP ... Phillips-Perronův neparametrický test, BG .. test autorů Bierens-Guo, KPSS ... test autorů Kwiatkowski, Phillips, Schmidt a Shin. U testů jsou uvedeny testační statistiky a kritické hodnoty pro hladinu významnosti α = 0.05. Pro testy nestacionarity je nulová hypotéza H0: jednotkový kořen s driftem, alternativní hypotéza H1: lineární TS proces – pro testy stacionarity jsou hypotézy obrácené. Závěr: Z tabulky lze vidět, že oba testy stacionarity dávají na uvedené hladině významnosti nesprávné závěry – neodhalily DS procesy a považují je za TS. Na hladině významnosti α =0.10 však už správně odmítají nulovou hypotézu. 3.2.2 Nestacionární procesy typu TS Tyto procesy jsou popsány jednoduchým modelem s explicitně vyjádřeným lineárním trendem: (15) X t = a + bt + ε t , X t = 0 = X 0 = a . Je logické, že při regresi s procesy tohoto typu bude regresorem i proměnná t - čas, nebo pořadí (jako relativní, normovaný čas). Na obr. 5 je uveden příklad TS procesů x1(t) a x2(t)

vygenerovaných v programu Matlab s parametry: n=100, seed=12345, σ ε =1, a=5, b=0.0 pro x1(t) a b=0.6 pro x2(t). *+,-./01234256>?6:41/+,;6<&8-9=6<"8-9 )! $! (! #! '! "! &! !

!

"!

#!

$!

%!

&!!

Obr. 5. Průběhy nestacionárních TS procesů x1(t) a x2(t) (horní graf) Identifikaci procesů lze uskutečnit pomocí testů jednotkových kořenů (testů nestacionarity a stacionarity) – tab. 4. Tab. 4. Výsledky testů nestacionarity a stacionarity TS procesů (EasyReg) Testy nestacionarity Testy stacionarity Proces ADF BG KPSS PP -1.909 > -3.45 -100.7 < -21.5 0.152 < 12.71 0.096 < 0.146 x1(t) -1.912 > -3.45 -101.8 < -21.5 0.129 < 12.71 0.095 < 0.146 x2(t) Závěr: Z tabulky lze vidět, že pro dané procesy test ADF dává nesprávné závěry – neodhalil TS procesy a považuje je za DS. Závěr: Z obou tabulek testů DS a TS procesů vyplývá, že nejspolehlivějším testem z uvedených je neparametrický Phillips-Perronův test. 3.3 Důsledky záměny typů procesů DS/TS Při záměně, nerozlišování, nebo nesprávné identifikaci typů nestacionárních procesů DS/TS dochází k určitým nesnázím [ARLT, J. 1999]: A. Použitím modelu TS při odstraňování trendu z DS procesu dochází k následujícím vážným problémům a chybám: 1. Index determinace R2 pro DS proces s c = 0 se pohybuje okolo hodnoty 0.44 bez ohledu na délku časové řady. Pro c≠0 se hodnota R2 zvyšuje s rostoucí délkou časové řady a dosahuje limitní hodnoty 1 bez ohledu na hodnotu konstanty, či variability časové řady

2. Rezidua získaná z tohoto modelu mají rozptyl v průměru ve výši 14% skutečného 2 rozptylu σ ε a to bez ohledu na hodnotu konstanty c, či délku časové řady 3. Průměrné hodnoty ACF reziduí jsou funkcí délky časové řady T a hodnoty oscilují s periodou asi (2/3)T, tzn. že tato funkce indikuje zdánlivý dlouhý cyklus. Tento výsledek nezávisí na tvaru ACF procesu ε t 4. Pro testování parametru u časové proměnné v modelu TS je t-test velmi slabý – v 87% případů dojde k zamítnutí správné nulové hypotézy b = 0 a v 73% dojde k zamítnutí při použití náhodné složky typu AR(p). B. Použití modelu DS na časovou řadu tvořenou procesem TS: V tomto případě je situace podstatně lepší. Odhad parametru c MNČ je nestranný a má přibližně normální (nebo Studentovo t) rozdělení. Snížení vydatnosti tohoto odhadu není vážný problém, pokud je současně modelována autokorelace nesystematické (náhodné) složky. Závěr: Z uvedeného vyplývá, že očištění obou typů časových řad diferencováním je univerzálnější a méně „nebezpečné”, než odstranění trendu. Z tohoto důvodu lze diferencování obecně použít na nestacionární procesy, i když neznáme jejich typ.

4 Problém zdánlivé korelace a regrese Definice: O zdánlivé korelaci mezi dvěma proměnnými (veličinami, časovými řadami, náhodnými procesy) mluvíme v případě, kdy mezi proměnnými skutečná závislost neexistuje, ale zdá se, že existuje. (Tento jev je duálním jevem ke zdánlivé nonkorelaci). Ve statistické literatuře se kromě názvu zdánlivá korelace (spurious correlation) používají ještě názvy klamná, pseudo, nepravá, iluzorní, nesmyslná (nonsense), náhodná a podmíněná (conditional) korelace. Vznik a rozčlenění: Podle příčin vzniku lze zdánlivou korelaci rozdělit v podstatě na tři skupiny, které se mohou i vzájemně kombinovat (kuriózní a zábavné příklady jsou z knihy [SWOBODA, H. 1977]): 1. Nesmyslná / náhodná / nevysvětlená korelace: žádná závislost nemůže z principu existovat (nebo o ní nevíme), ale jenom náhodně mají obě proměnné podobný průběh. Například: „čím kratší sukně, tím vyšší burzovní kurzy”, „čím větší dovoz pomerančů, tím větší počet úmrtí na rakovinu”, „lidé trpící ledvinovými kameny mají oblíbenou barvu zelenou” 2. Podmíněná korelace: zdánlivá závislost je způsobena vlivem společné třetí proměnné, společného faktoru (který na obě sledované proměnné nějak působil). Nejčastěji je společnou (třetí) proměnnou: čas, pořadí, trend (u časových řad a náhodných procesů), teplota, tlak, referenční/vztažná veličina, akční veličina, funkce svazující obě veličiny. Například: „s přibývajícím věkem ženy kladou chodidla více k vnější straně” (vliv výchovy na přelomu 19. a 20. století), „mezi podílem kuřáků a nekuřáků existuje výrazná záporná korelace” (procentní podíly se navzájem doplňují na 100 %), „čím mají děti vyšší tělesnou váhu, tím jsou zručnější” (korelace obou znaků s věkem: starší děti jsou šikovnější i těžší v poměru k mladším). 3. Nesprávně interpretovaná korelace: zdánlivá korelace pochází buď z nesprávného výkladu skutečně existující korelace, nebo je zaměněna orientace kauzality (obrácení příčiny a účinku). Například: „mezi tělesnou váhou a četnými nemocemi existuje kladná korelace. Obézní budou proto zdravější, budou-li užívat odtučňovací tablety”, „nepoměrně více lidí umírá v posteli než na ulici, na moři nebo ve vykřičené krčmě. Postel je proto prokazatelně nejnebezpečnější místo pobytu!”.

Nebezpečí zdánlivé korelace spočívá v tom, že vede na ukvapené a nesprávné závěry (včetně obrácení orientace kauzality), nebo podporuje manipulované interpretace (v politice, reklamě apod.). Detekce: Indikaci zdánlivé korelace lze uskutečňovat: ! pro průřezová data pomocí srovnání znamínek a velikostí hodnot korelační a parciální korelační matice [HEBÁK, P. & HUSTOPECKÝ, J. 1990], [MELOUN, M. & MILITKÝ, J. 1994] ! u časových řad (náhodných procesů) pomocí srovnání výsledků regrese původních a upravených (stacionarizovaných / očištěných, nebo doplněných o společný faktor) dat. Pokud je lineární regresní člen pro původní data statistický významný a pro upravená data statisticky nevýznamný, pak jde o zdánlivou korelaci a regresi [SEGER, J., HINDLS, R. & HRONOVÁ, S. 1998] ! u nestacionárních integrovaných DS / I(d) procesů byl zaveden pojem kointegrace, která slouží k odlišení jejich zdánlivé a skutečné korelace. Řešení zdánlivé korelace obecně spočívá ve: ! vyloučení společného faktoru (pomocí stacionarizace: odstranění trendu, cykličnosti, sezónnosti / odstranění driftu - diferencováním) [BAKYTOVÁ, H. AJ. 1979], [CYHELSKÝ, L. AJ. 1986], [SEGER, J., HINDLS, R. & HRONOVÁ, S. 1998] a následné testování korelace mezi nepravidelnými (náhodnými) složkami analyzovaných dat, nebo ! vyjádření společného faktoru (jeho zavedením do regresní rovnice) [KAŇOKOVÁ, J. 1989], [BIERENS, H. 1999c], [HEBÁK, P. & HUSTOPECKÝ, J. 1990], [MELOUN, M. & MILITKÝ, J. 1994]. Pozn.: V předchozí kapitole bylo konstatováno, že diferencování je poměrně obecnější způsob stacionarizace nestacionárních procesů, než odstranění deterministického trendu. Tento způsob má však omezení: u systému nestacionárních procesů existuje mnohem univerzálnější a účinnější způsob testování i řešení vzájemných vazeb pomocí kointegrace (nebo pomocí nelineární kotrendové analýzy). Příklad zdánlivé korelace dvou TS procesů je uveden např. v [SEGER, J., HINDLS, R. & HRONOVÁ, S. 1998]: „V luxusní francouzské restauraci jednoho pražského hotelu byla v deseti dnech za sebou sledována tržba za poledního provozu (obědy Xt – tis. Kč) a za večerního provozu (večeře Yt – tis. Kč). Je třeba zjistit, jak je silná závislost mezi vývojem tržeb v poledním a večerním provozu, nebo jinak řečeno – zda by se dalo z vývoje poledních tržeb usuzovat na očekávanou výši tržeb večerních”. Časové řady jsou zobrazeny na obr. 6, odkud se nám jeví, že spolu poměrně těsně souvisí – mají podobný trend. Graf průběhu tržeb 190

46

188 44 186 184

42

182 40

180 POLEDNE (L) VEČER (R)

38

178 176

36

0

2

4

6

8

DEN

Obr. 6. Průběhy tržeb v poledne a večer (v tis. Kč)

10

174

A. Průřezová (prostorová) data – korelace: I když je známo, že obě veličiny jsou časové řady, přesto kvůli demonstraci bude zajímavé pohlížet na ně jako na průřezová data a otestovat jejich vzájemné korelační vztahy (nakonec i samotní renomovaní autoři [SEGER, J., HINDLS, R. & HRONOVÁ, S. 1998] zvolili právě tento přístup). Pro otestování přítomnosti zdánlivé korelace je třeba vypočíst matici výběrových (párových) korelačních koeficientů a matici výběrových parciálních korelačních koeficientů [HEBÁK, P. & HUSTOPECKÝ, J. 1990], [MELOUN, M. & MILITKÝ, J. 1994]. Parciální korelační koeficient je korelační koeficient počítaný mezi dvěma veličinami s vyloučením vlivu ostatních proměnných [ANDĚL, J. 1985]. Umožňuje tedy odhalit souvislosti, které by jinak zůstaly skryté v souvislostech „všeho se vším”, anebo naopak se ukáže, že nějaká zdánlivě výrazná souvislost (zdánlivá korelace), indikována korelačním koeficientem, je vlastně jen zprostředkována jinými proměnnými (společnou třetí proměnnou) [KOSCHIN, F. AJ. 1992]. Porovnání hodnot obou typů korelačních koeficientů lze také využít k určení parazitních a významných proměnných při vícenásobné regresi [MELOUN, M. & MILITKÝ, J. 1994]. Pro dané veličiny – tržby v poledne (P), večer (V) a čas, dny (D) - dostaneme následující matice párových (horní hodnoty) a parciálních (dolní hodnoty) korelačních koeficientů (tab. 5 – QC Expert): Tab. 5. Matice korelačních koeficientů veličiny P V D P 1 +0.73 V 1 (-0.06) +0.79 +0.94 D 1 (+0.45) +0.86 Pozn.: Hodnoty koeficientů byly zaokrouhleny na 2 desetinná místa, hodnoty v závorkách jsou statisticky nevýznamné na hladině významnosti α =0.05. Kritická hodnota párových korelačních koeficientů r(n=10, α =0.05) = 0.63 a parciálních korelačních koeficientů r’(9, 0.05) = 0.67 [ANDĚL, J. 1985]. Nevýznamnost parciálního korelačního koeficientu r’ proměnných P, D je pravděpodobně způsobena malým rozsahem výběru (n=10). Z tabulky je zřejmé, že statisticky nevýznamná hodnota i jiné znaménko (než u párového) parciálního korelačního koeficientu veličin P,V indikuje zdánlivou (podmíněnou, zprostředkovanou) korelaci. Stejnou hodnotu lze dostat po výpočtu párového korelačního koeficientu pro veličiny očištěné od lineárního trendu (korelační matice pouze proměnných P,V) [SEGER, J., HINDLS, R. & HRONOVÁ, S. 1998]. Situaci při zdánlivé korelaci analyzovaných veličin znázorňuje obr. 7. P

r = +0.73

V

akce v čase

společný faktor

V

P

neuvažován (=> zdánlivá korelace)

uvažován

r' = (-0.06)

Obr. 7. Schéma zdánlivé a skutečné korelace mezi veličinami P a V

B. Časové řady – regrese: Výsledky lineární regrese pro původní data jsou uvedeny v tab.6, přičemž lineární statický regresní model byl uvažován ve tvaru: (16) Yt = b0 + b1 X t + ε t , t = 1, 2, ... 10 . Tab.6. Výsledky lineární regrese pro původní neupravená data (EasyReg) b1 F R2 [%] DW JB BP Pozn. b0 136.7 1.093 9.01 53 1.99 1.044 0.044 vše OK Pozn.: F ... Fisherova statistika modelu, DW ... statistika Durbin-Watsonova testu autokorelace reziduí, JB ... statistika Jarque-Beraova testu normality reziduí, BP ... statistika Breusch-Paganova testu homoskedasticity reziduí. Hladina významnosti α = 0.05. Model je po všech stránkách naprosto korektní i významný a tak by se zdálo, že hypotéza o závislosti večerních a poledních tržeb je prokázána. Podezření na zdánlivou korelaci však vyplývá z několika skutečností: - věcná úvaha o reálnosti výsledku: není důvod očekávat silnější vztah mezi polední a večerní tržbou, čili asi jde o nesmyslnou korelaci, - zřejmě působil společný trendový faktor např. začátek sezóny, vliv reklamy apod., takže jde také o podmíněnou korelaci, - fakt, že skutečná korelace se projevuje mezi nesystematickými (náhodnými) složkami časových řad a ne mezi složkami systematickými (trend, cykličnost, sezónnost). Existuje více možností úprav dat: v případě, že jde o procesy DS (po otestování jednotkových kořenů) použijeme pro regresi diferencované řady. Pokud však umíme otestovat kointegraci, tak ji u těchto řad napřed použijeme a rozhodneme o dalším postupu. Pokud jde o procesy TS, pak buď provedeme očištění od trendu, nebo trendovou proměnnou (datum, čas, pořadí) zařadíme do modelu. Pokud mají řady smíšený DS/TS charakter, nebo charakter neznáme, pak použijeme diferencované řady. Za předpokladu nepoužití (neznalosti) testů jednotkových kořenů a kointegrace, provedeme regresní analýzu s daty očištěnými od lineárního trendu, se zavedením relativního času (pořadí, dny) a také s diferencovanými hodnotami. Výsledky jsou uvedeny v tab.7,8,9 (hodnoty v závorkách jsou statisticky nevýznamné na hladině α = 0.05). Tab. 7. Výsledky regrese s daty očištěnými od lineárního trendu (Statgraphics, EasyReg) b0 b1 F R2 [%] DW JB BP (0.000) (-0.050) (0.027) (0.3) 2.015 0.117 0.311 Model s trendově očištěnými veličinami (odchylkami od trendů) měl tvar: (17) e y = b0 + b1e x + ε t ⇒ (Yt − T y ) = b0 + b1 ( X t − Tx ) + ε t . Lineární trendové funkce pro obě veličiny byly nalezeny ve formě: (18) T y = 173.933 + 1.394t , (19) Tx = 36.800 + 0.782t . Výsledky signalizují, že model je korektní, ale statisticky nevýznamný včetně svých parametrů. Tato skutečnost znamená, že mezi veličinami je pouze zdánlivá (nesmyslná) korelace. Tab. 8. Výsledky lineární regrese se zavedenou další proměnnou – časem (EasyReg) b0 b1 bt F R2 [%] DW JB BP 1.433 (-.050) 175.8 25.49 88 2.015 0.117 0.996 Při této regresi byl uvažován model: (20) Yt = b0 + b1 X t + bt t + ε t .

Model je naprosto korektní, ale lineární koeficient b1 vyjadřující závislost na poledních tržbách je statisticky nevýznamný. Srovnáním s výsledky modelu pro původní data lze tedy konstatovat, že jde o zdánlivou (podmíněnou) korelaci. Tab. 9. Výsledky lineární regrese pro diferencované časové řady (EasyReg) b0 b1 F R2 [%] DW JB BP (1.625) (-0.325) (1.16) (14) 2.695 0.826 2.520 Byl uvažován model s diferencemi veličin: (21) ∆Yt = b0 + b1 ∆X t + ε t , t = 1, 2, ... 9 , přičemž z výsledků vyplývá, že model je korektní (kritická hodnota Durbin-Watsonovy statistiky pro negativní autokorelaci je 4-dL( α =0.05, n=9, k=1) = 3.176 [CYHELSKÝ, L. AJ. 1986]), ale přitom statisticky nevýznamný i se svými koeficienty. Závislost mezi řadami je tedy pouze zdánlivá.

5 Kointegrace Pojem kointegrace jako první teoreticky popsal Granger (1981), který se zajímal o zdánlivou závislost mezi nestacionárními procesy. Princip kointegrace je jednoduchý: i když se procesy chovají individuálně (jsou nestacionární), přesto mezi nimi může existovat vztah, který se chová ustáleně (je stacionární). Jinak řečeno: proměnné se nazývají integrovanými, pokud jejich trend může být odstraněn diferencováním. Integrované proměnné jsou kointegrované, když existuje jejich lineární kombinace, která nemá stochastický trend [BENKWITZ, A. AJ. 1999]. Odborně řečeno: i když všechny komponenty vektorového náhodného procesu zt mají jednotkové kořeny (tj. když zt je vícenásobný integrovaný proces typu I(1)), přesto mohou existovat lineární kombinace ξTzt bez jednotkového kořene (tj. jsou stacionárními I(0) procesy). Tyto lineární kombinace mohou být interpretovány jako dlouhodobé vztahy mezi komponentami zt, nebo podle ekonometrické terminologie jako statická ekvilibria [BIERENS, H. 1999a]. Princip kointegrace je ústřední myšlenkou modelování integrovaných časových řad, protože: 1. Střední hodnotu stacionární lineární kombinace integrovaných časových řad I(d) je možné chápat jako (dlouhodobé) ekvilibrium (rovnovážný vztah), které spojuje uvažované časové řady. Analýza vztahů mezi procesy I(d) má smysl pouze tehdy, jsou-li tyto kointegrované, tj. 2. jsou-li spjaté společným stochastickým stacionárním trendem. Není-li tomu tak, každý proces má jiný směr vývoje. Potom při zkoumání vztahů mezi nimi pomocí regresní analýzy vzniká stav tzv. zdánlivé regrese. Test kointegrace procesů je tedy současně metodou pro rozlišení mezi pravou a zdánlivou regresí (i korelací). 3. Skupinu kointegrovaných časových řad lze popsat modelem korekce chyby (ECM – Error Correction Model), jehož prostřednictvím lze odlišit dlouhodobé a krátkodobé vztahy mezi časovými řadami. Může být prostředkem řešení rozporu mezi statistickým a ekonometrickým přístupem modelování nestacionárních časových řad. Oba přístupy použité izolovaně jsou problematické. ECM umožňuje spojit přístup statistický (zkoumající vlastnosti stacionarizovaných diferencovaných časových řad, ale zbavující se důležitých informací obsažených v původních nestacionarizovaných časových řadách) a přístup ekonometrický (kladoucí důraz na ekvilibrium časových řad, ale přehlížející problém zdánlivé regrese). Systém kointegrovaných procesů se nejčastěji popisuje pomocí k-rozměrného vektorového autoregresního modelu řádu p - VAR(p) typu I(1), který lze vyjádřit v diferenční formě tzv. modelu korekce chyb – ECM (Error Correction Model) [ARLT, J. 1999].

Model ECM obsahuje jak krátkodobé vztahy mezi procesy (tj. vztahy mezi stacionarizovanými, diferencovanými procesy), tak i vztahy dlouhodobé (tj. vztahy mezi nediferencovanými procesy). Informace o těchto vztazích jsou obsaženy v parametrické matici Π. Konstrukce modelu EC umožňuje oddělit oba druhy vztahů a zkoumat je samostatně, přičemž mohou nastat tři situace: 1. h(Π) = r = k, tj. matice Π má plnou hodnost (počet kointegračních vektorů r je roven počtu procesů k). To znamená, že k-rozměrná časová řada je generována stacionárním vektorovým procesem {Xt}(systém obsahuje pouze procesy typu I(0) a ne I(1)). 2. 0 < h(Π) = r < k, v tomto případě nezmizí nediferencovaný člen modelu EC, ale současně nelze vektorový proces {Xt} považovat za stacionární. Protože matice Π je nenulová, lze najít mezi časovými řadami dlouhodobý vztah a stacionarizaci jejich individuálním diferencováním bez ztráty informace nelze provést. Prakticky to znamená, že některé časové řady lze stacionarizovat jejich diferencováním, neboť nejsou obsaženy v žádném dlouhodobém vztahu s jinými. Některé časové řady však nemohou být stacionarizovány diferencováním, neboť jejich lineární kombinace s jinými časovými řadami již stacionární jsou – tyto řady jsou kointegrované. Touto situací se detailně zabývá Grangerova věta (1987), která dokazuje, že kointegrovaný systém časových řad může být vyjádřen ve třech formách: ve formě modelu VAR, EC a VMA. 3. h(Π) = r = 0, tzn. že matice Π je nulová a model EC neobsahuje nediferencovaný člen. Krozměrná časová řada je generována nestacionárním vektorovým procesem {Xt} a její stacionarizaci lze provést individuálním diferencováním jednotlivých časových řad. Diferencováním nedochází ke ztrátě informace o dlouhodobém vztahu mezi časovými řadami, neboť žádný neexistuje. Kointegraci lze definovat taktéž v jednorovnicových modelech. Princip kointegrace byl detekován a aplikován v oblasti ekonometrie (makroekonomicko-finanční aplikace), jeho rozpracování a aplikace v oblasti technometrie (jakožto průniku matematiky, matematické statistiky a technologie) je teprve v počáteční fázi. 5.1 Ekonometrický pohled V příspěvku [JOHANSEN, S. 1997] je uveden jednoduchý a přehledný systém kointegrovaných procesů Yt, Zt ve tvaru: (22) Yt = (1 − a )Yt −1 + aZ t −1 + ε 1t ,  Z t = aYt −1 + (1 − a ) Z t −1 + ε 2 t , přičemž pro koeficient vazby platí, že a ∈<0, 1>. Pro a = 0 jsou oba procesy typu náhodná procházka (integrované procesy řádu I(1)) a nejsou vzájemně kointegrované. Pro tento případ však může nastat jev zdánlivé korelace (zvláště pro větší počet hodnot) [BIERENS, H. 1999c]. Pro a = 1 jsou procesy prakticky totožné (proces je roven druhému zpožděnému procesu s přídavným šumem). Schéma vazeb vektorového systému je znázorněna na obr. 8.

Y0

ε1t

Σ

a

Yt

1-a

z-1

Z0 a

ε2t

Σ

Zt

z-1

1-a

Obr. 8. Blokové schéma kointegrovaných procesů Yt, Zt (Johansen) V programu Matlab lze vytvořit M-funkci, např. Coin2_J se vstupními parametry: {n, seed, µ ε 1 , µ ε 2 , σ ε 1 , σ ε 2 , a, Y0, Z0} a s výstupy {Yt, Zt} včetně jejich grafického průběhu a uložení hodnot do ASCII souboru. Pro parametry n = 100, seed = 13579463, µ ε 1 = µ ε 2 = 0 , σ ε 1 = σ ε 2 = 1, a = 0 (tj. procesy jsou nezávislé, nekointegrované), Y0 = 15 a Z0 = 5 byly vygenerovány integrované procesy a jejich průběhy jsou znázorněny na obr.9. 4

e1t

3

e2t

2 1 0 -1 -2 0

10

20

Yt

30

40

50

60

70

80

90

100

40

50

60

70

80

90

100

Zt

20 10

0

10

20

30

Obr. 9. Průběhy integrovaných procesů Yt (horní průběh), Zt při zdánlivé korelaci (a = 0) Při vizuálním „expertním” posouzení průběhů se zdá, že mezi procesy existuje poměrně těsná záporná korelace, což nakonec celkem „potvrdí” i korelační diagram (obr. 10).

12

Zt x Yt

10

8

6

4

2 13

14

15

16

17

18

19

20

21

22

23

24

25

26

Obr. 10. Korelační diagram procesů Yt a Zt (při zdánlivé korelaci) Pokud uskutečníme tzv. naivní regresní analýzu (tj. bez předchozího ověření předpokladů), dostaneme pro model Zt = a + b.Yt výsledky (tab. 10), které nám „potvrdí” zdánlivou korelaci jako korelaci skutečnou - model i koeficienty jsou statisticky významné, tj. mezi veličinami existuje lineární závislost. Tab. 10. Základní výsledky jednoduché lineární regrese (Statgraphics) Parametr / statistika Hodnota Poznámka / závěr absolutní člen a 18.108 statisticky významný lineární člen b -0.545 statisticky významný Fisherův F-test modelu 38.895 statisticky významný koeficient determinace R2 0.284 statisticky významný Pozn.: Uvažovaná hladina významnosti α = 0.05, celkový F-test modelu je v tomto případě ekvivalentní testu koeficientu determinace, výsledky byly zaokrouhleny na 3 desetinná místa. Teprve při analýze reziduí (tab.11) vyplyne nekorektnost modelu způsobena výraznou pozitivní autokorelaci reziduí (kritická hodnota Durbin-Watsonovy statistiky dL(α=0.05, n=100, k=1) je 1.65). Tab. 11. Analýza reziduí jednoduché lineární regrese (EasyReg) Vlastnost Test Hodnota Poznámka / závěr 0.193 pozitivní autokorelace autokorelace Durbin – Watson normalita Jarque – Berra 5.368 akceptována homoskedasticita Breusch – Pagan 0.596 akceptována Autokorelace reziduí signalizuje nesprávný model, který však nemá význam hledat prostředky vhodnými pro statistickou analýzu stacionárních časových řad. Při analýze nestacionárních procesů, nejprve otestujeme tzv. jednotkové kořeny (pro odlišení TS/DS procesů) a v případě integrovaných procesů (DS, I(1)-náhodné procházky) testujeme jejich kointegraci (tab.12,13). Tab. 12. Výsledky testů jednotkových kořenů (EasyReg) Test Proces Poznámka / závěr ADF PP Yt -2.177 > -3.45 -8.852 > -21.5 nestacionární proces ~ I(1) Zt -2.061 > -3.45 -6.895 > -21.5 nestacionární proces ~ I(1)

Pozn.: ADF ... Rozšířený (Augmented) Dickey-Fullerův test, PP ... Phillips-Perronův neparametrický test. U testů jsou uvedeny testační statistiky a kritické hodnoty pro α = 0.05. Nulová hypotéza H0: jednotkový kořen s driftem, alternativní hypotéza H1: lineární TS proces. Tab. 13. Výsledky testů kointegrace (EasyReg) Řád VAR systému Ekvilibrium Poznámka / závěr p=1 p=2 konstanta r=0 r=0 r = 0 => procesy nejsou konstanta + kointegrované r=0 r=0 lineární trend Pozn.: VAR ... Vector AutoRegressive. Počet kointegrovaných vektorů r byl stanoven podle Johansenova přístupu. Závěr: Oba procesy jsou nestacionární typu I(1) (náhodná procházka) a nejsou kointegrované, tj. neexistuje mezi nimi dlouhodobý vztah (ekvilibrium). Testování kointegrace lze využít pro rozlišení mezi skutečnou a zdánlivou korelací mezi nestacionárními procesy typu DS / I(1). Nekointegrované procesy lze jednotlivě stacionarizovat diferencováním. Mezi diferencemi procesů Yt a Zt, tj. mezi procesy ε1t a ε2t se už žádným způsobem (vizuálně, ani klasickou regresní analýzou) nepodaří zjistit korelaci – viz horní graf obr. 9. 5.2 Technometrický pohled Nevýhodou ekonometrického přístupu pro techniky jsou těžko interpretovatelné, umělé a nevhodné modely pro signály dynamických soustav technologických procesů. V následující části budou definovány diskrétní modely soustav Si1 a Sp1, jako i způsob testování jejich rozlišení (identifikace) na základě vstupních a výstupních signálů. 5.2.1 Diskrétní modely soustav Si1 a Sp1 Základním zdrojem nestacionarity technologických procesů jsou integrační (astatické) soustavy – např. různé zásobníky, nádrže, hromady apod. V dalším budeme uvažovat základní integrační soustavu s astatismem 1. řádu (Si1) se spojitým přenosem (23) 1 GS ( s) = , TI s kde TI je integrační časová konstanta soustavy. Ze spojitých přenosů lze získat tzv. celkové diskrétní přenosy - invariantní vzhledem k přechodové funkci, nebo jinak přenosy s uvažovaným A-Č převodníkem na vstupu (vzorkovač) a Č-A převodníkem na výstupu (vzorkovač a tvarovač 0. řádu) - pomocí vztahu uvedeného např. ve [VÍTEČEK, A. 1988]: (24) G (s) G SC ( z ) = (1 − z −1 ) Z {L−1 { S } t = kT } , s kde t = k.T je diskrétní čas v násobcích k = {0, 1, 2, ...} periody vzorkování T. Pro uvažovanou integrační soustavu tak dostaneme celkový diskrétní přenos (model) (25) T z −1 G SC ( z ) = TI 1 − z −1 s odpovídající diferenční rovnicí: (26) T y t = a1 y t −1 + b1u t −1 = y t −1 + u t −1 , TI odkud je zřejmé, že pro koeficienty platí: (27) T . a1 ≡ 1, b1 = TI

Pro diferenci výstupního signálu platí vztah odpovídající rovnici přímky procházející počátkem (28) T ∆y t = y t − y t −1 = b1u t −1 = u t −1 , TI ze kterého je zřejmá jeho použitelnost (ve smyslu výsledků lineárního regresního modelu) pro parametrickou identifikaci soustavy Si1, tj. určení její integrační časové konstanty (perioda vzorkování T je obecně známá) (29) T TI ≈ . b1 Další soustavou, která může být aproximací integrační soustavy 1.řádu, je proporcionální soustava se setrvačností 1.řádu (Sp1) se spojitým přenosem (30) k1 G S (s) = , T1 s + 1 kde k1 je koeficient přenosu a T1 časová konstanta soustavy. Tato soustava má celkový diskrétní přenos: T (31) − k1 (1 − c) z −1 T1 G SC ( z ) = c e , ( 0 , 1 ) = ∈ 1 − cz −1 a odpovídající diferenční rovnici: (32) y t = a1 y t −1 + b1u t −1 = cyt −1 + k1 (1 − c)u t −1 , odkud je zřejmé, že pro koeficienty platí: T (33) − T1 a1 = c = e ∈ (0, 1), b1 = k1 (1 − c) > 0 . Odhady parametrů soustavy z odhadů regresních koeficientů lze získat pomocí vztahů: ~ ~ (34) b T k1 = 1 , T1 = − . 1 − a1 ln(a1 ) Pro tzv. kvazidiferenci (viz např. [GARAJ, V. & ŠUJAN, I. 1980]) výstupního signálu platí také vztah odpovídající rovnici přímky procházející počátkem ~ y = y − cy = k (1 − c)u , (35) ∆ 1 t t t −1 t −1 odkud lze vidět, že kvazidiference pro c = 1 přechází na diferenci procesu. Ze srovnání diferenčních rovnic je zřejmé, že obě soustavy mají velice podobné diskrétní modely (obr.11). ut

b1 z-1

k1(1-c)

y0

Σ

b1 z-1

T/TI

y0

Σ

yt

a1

a1 c

yt

ut

z-1

a1 = c = e-T/T 1 ε (0, 1)

1

z-1

a1 = 1

Obr. 11. Diskrétní modely soustav Sp1 (vlevo) a Si1 (vpravo) Z hlediska klasifikace procesů podle moderních a klasických metod analýzy časových řad bude mít výstup yt soustavy Sp1 charakter stacionárního procesu I(0) / AR(1) a výstup soustavy Si1 charakter nestacionárního procesu I(1) / ARIMA(0,1,0) – za předpokladu, že na vstupy soustav bude působit signál ut typu náhodný Gaussovský proces (bílý šum).

Pokud zobrazíme závislost koeficientů diferenční rovnice a1, b1 diskrétního modelu soustavy Sp1 na parametrech soustavy k1, T1 a modelu T, dostaneme pro normovaný poměr časových parametrů T1/T následující průběhy (obr. 12): 1

1.5 b1( 0.5, T1, T )

1

b1( 1 , T1, T )

a1( T1, T ) 0.5

b1( 1.5, T1, T ) 0.5 0

0.1

1

10

100

0

0.1

1

10

T1

T1

T

T

100

Obr. 12. Závislosti koeficientů a1(T1,T) a b1(k1,T1,T) na poměru T1/T Z obrázku je jasné, že pro poměr časových parametrů T1/T > 10 (čili pro časovou konstantu T1 řádově větší než perioda vzorkování T) se hodnota koeficientu a1 → 1 a soustava Sp1 se přibližuje svým chováním soustavě Si1. Pro uvedený velký poměr časových parametrů budou pravděpodobně selhávat testy typů procesů DS/TS (vycházející z klasické analýzy časových řad) a testy jednotkových kořenů, které zřejmě nerozliší stacionární a nestacionární procesy výstupních signálů těchto soustav. Podobnost výstupních procesů diskrétních modelů soustav DSi1 a DSp1 se v případě a1 = c → 1, TI >> T projeví i na jejich diferencích (36) T ∆y t = u t −1 → 0 , TI (37) ∆y t = k1 (1 − c)u t −1 → 0 . 5.2.2 Identifikace soustav Si1 a Sp1 Identifikace uvedených soustav bude uskutečněna pomocí testování nestacionarity a typů procesů jejich výstupních signálů. Pro další analýzu budou vygenerovány odezvy soustav Si1 a Sp1 na stejný vstupní stacionární signál typu Gaussovský náhodný proces (bílý šum) s nenulovou střední hodnotou. Budou prozkoumány možnosti určení typů výstupních nestacionárních procesů (DS/TS) a určení typů soustav (Si1/Sp1), ze kterých procesy pocházejí (vycházejí). V programu Matlab lze vytvořit M-funkce pro diskrétní modely soustav Si1 a Sp1, např. DSi1 se vstupními parametry: {n, u, TI, T, y0} s výstupem yI, DSp1 se vstupními parametry: {n, u, k1, T1, T, y0} s výstupem yP, včetně jejich grafického průběhu a uložení hodnot do ASCII souboru. Pro parametry simulace: n=100, T=1, parametry vstupního stacionárního náhodného procesu s nenulovou střední hodnotou ut: seed=135791, µ u = 5, σ u = 1, byly vygenerovány: výstupní proces yI soustavy Si1 s parametry: TI=500, yI0=12 a výstupní proces yP soustavy Sp1 s parametry: k1=2.3, T1=100, yP0=10. Konfigurace diskrétních modelů soustav je znázorněna na obr.13. yIt DSi1 ut yPt DSp1 Obr. 13. Konfigurace modelů soustav DSi1 a DSp1

Průběhy vstupního stacionárního a výstupních nestacionárních procesů lze sledovat na obr.14. 9

ut

8 7 6 5 4 3 0 13

10 yI

20

30

40

50

60

70

80

90

100

30

40

50

60

70

80

90

100

yP

12 11 10

0

10

20

Obr. 14. Průběhu vstupního (ut) a výstupních procesů (yI a yP) Výstupní procesy soustav lze testovat a identifikovat pomocí: - klasické analýzy časových řad, - moderní metody testování jednotkových kořenů - regresní a korelační analýzy (tj. pomocí regresních rovnic a korelačních grafů). A. Klasická analýza časových řad: Při testování procesů pomocí klasické analýzy časových řad testujeme původní, diferencované a trendově očištěné procesy – podle dříve uvedené metodiky ověřované na generovaných DS a TS procesech. Výsledky analýzy jsou uvedeny v tab. 14 (Statgraphics). Tab. 14. Základní klasifikace nestacionárních procesů yI(t) a yP(t) Procesy Hodnoty Model Závěr yI I(1) procesy nelze rozlišit! původní yP I(1) bílý šum (BŠ) diferencované procesy diff yI diferencované nelze rozlišit! bílý šum (BŠ) diff yP AR(1) / I(1) trendově očištěné očištěné od lineárního detrend yI procesy nelze rozlišit! trendu (Tt) detrend yP AR(1) / I(1) Závěr: Jak je z tabulky zřejmé, pomocí klasické analýzy časových řad nelze rozlišit typy výstupních procesů DS/TS soustav Si1 a Sp1! B. Moderní analýza časových řad: Pokud použijeme testy jednotkových kořenů (a to testy nestacionarity i stacionarity) pro původní, diferencované a trendově očištěné procesy, dostaneme výsledky uvedené v tab. 15.

Tab. 15. Výsledky testů nestacionarity a stacionarity výstupních procesů (EasyReg) Testy nestacionarity Testy stacionarity Proces Poznámka ADF PP BG KPSS TS yI DS DS DS procesy nelze rozlišit! TS yP DS DS DS TS diff yI DS TS TS diferencované procesy nelze rozlišit! TS diff yP DS TS TS TS detrend yI DS DS DS trendově očištěné procesy nelze rozlišit! TS detrend yP DS DS DS Pozn.: Testy PP a KPSS jsou silnější než testy ADF a BG. Pomocí testů PP a KPSS bylo možné odlišit pouze diferencované procesy (které se jevily jako procesy TS) od původních a trendově očištěných procesů (které měly charakter procesů DS, tj. náhodných procházek typu I(1)). Závěr: Z výsledků v tabulce je možné učinit závěr, že pomocí testů jednotkových kořenů také nelze rozlišit výstupní procesy uvedených soustav Si1 a Sp1! C. Regresní a korelační analýza časových řad: Pro otestování (identifikaci) typů soustav (Si1/Sp1) zkusíme použít regresní a korelační analýzu (tj. regresní rovnice a korelační grafy). Můžeme přitom využít buď: ! výše uvedené vztahy pro (kvazi)diference výstupních signálů soustav ve formě výpočtů regresních rovnic odpovídajících diskrétnímu modelu soustavy Si1 ! nebo regresní rovnice odpovídající diskrétnímu modelu soustavy Sp1. 1) Regresní rovnice typu Si1 a) Pro testování použijeme nejprve příslušný lineární regresní model bez absolutního členu (který však přesně odpovídá pouze modelu soustavy DSi1) (38) ∆y t = b1u t −1 + ε t , kde ε t označuje rezidua. Výsledky regrese jsou uvedeny v tab. 16. Tab. 16. Výsledky regrese pro model bez absolutního členu (Statgraphics, EasyReg) R2 [%] DW JB Homoskedasticita Poznámka Procesy b1 2.823 1.662 + perfektní model 100 ∆yI, ut-1 0.0020 33 2.078 2.594 + korektní model ∆yP, ut-1 0.0028 Závěr: Oba modely jsou statisticky významné i korektní (správné), přičemž model pro výstupní signál ze soustavy DSi1 je perfektní (R2 = 100 %). Z této skutečnosti můžeme usoudit, že signál yI pochází z integrační soustavy. b) Pokud uskutečníme výpočet pro umělý (neodpovídající modelům dynamických soustav) lineární regresní model s absolutním členem, dostaneme výsledky uvedené v tab.17. Regresní model má v tomto případě tvar (39) ∆y t = a 0 + b1u t −1 + ε t . Pro soustavu DSi1 dostaneme rezidua ε t odečtením rovnice diskrétního modelu soustavy od uvedené regresní rovnice (40) ε t = −a0 , odkud je vidět, že rezidua jsou nezávislá na vstupní (vysvětlující) proměnné ut-1 i na výstupní (vysvětlované) proměnné ∆yI , mají charakter Gaussovského náhodného procesu se střední

hodnotou blízkou nule. To dále znamená, že absolutní člen a0 by měl být statisticky nevýznamný. Pro soustavu DSp1 upravíme nejprve rovnici diskrétního modelu na tvar s diferencí výstupního signálu na levé straně (41) ∆y t = y t − y t −1 = (c − 1) y t −1 + k1 (1 − c)u t −1 = (c − 1) y t −1 + b1u t −1 , odkud rezidua ε t dostaneme odečtením této rovnice od uvažované regresní rovnice s absolutním členem (42) ε t = (c − 1) y t −1 − a0 . V tomto případě je zřejmé, že rezidua jsou závislá na zpožděné (o jeden krok) výstupní (vysvětlované) proměnné yIt-1, z čehož vyplývá jejich autokorelace 1.řádu. Jejich obecně nenulová hodnota je závislá ne velikosti poměru časových parametrů T1/T. Absolutní člen a0 bude nenulový a zřejmě i statisticky významný. Tab. 17. Výsledky regrese pro model s absolutním členem (Statgraphics, EasyReg) b1 R2 [%] DW JB BP Poznámka Procesy a0 absolutní člen je 2.823 1.669 0.078 100 ∆yI, ut-1 (3.10-9) 0.002 nevýznamný autokorelace, 98.6 (0.016) (8.760) 0.245 ∆yP, ut-1 -0.107 0.023 nenormalita! Závěr: Model pro výstupní signál yP ze soustavy DSp1 je statisticky významný, ale není korektní (správný) – rezidua vykazují významnou pozitivní autokorelaci a nejsou normálně rozdělena. Koeficient b1 pro tento model je asi o jeden řád větší, než v modelu bez absolutního členu. Tato skutečnost také podporuje závěr, že s modelem pro yP není něco v pořádku. Model pro výstupní signál yI ze soustavy DSi1 je perfektní (R2 = 100 %) a korektní (správný). Z této skutečnosti můžeme znovu usoudit, že signál yI pochází z integrační soustavy. Uvedené skutečnosti můžeme podchytit v algoritmu identifikace typu soustav / jejich procesů (obr. 15). TS

Z ∆yt = (a0) + b1ut-1 + ε

Si1

+

perfektní model ?

Sp1

K

Obr. 15. Algoritmus testování typů soustav (jejich výstupních procesů) Pro vizuální posouzení lze zobrazit korelační diagramy (bodové X-Y grafy) pro odpovídající proměnné (signály) – viz obr. 16.

Obr. 16. Korelační diagramy diferencí výstupů a zpožděného vstupu soustav Závěr: Z grafů a předchozích vztahů je zřejmé, že přesná lineární závislost mezi diferencí výstupního signálu a o jeden krok zpožděným vstupním signálem odpovídá pouze integrační soustavě Si1. Tato skutečnost je tedy kritériem pro rozlišení typů soustav a tím i typů signálů. Výstupní signál yI z integrační (astatické) soustavy Si1 je pravým nestacionárním signálem, zatímco výstupní signál yP z proporcionální (statické) soustavy Sp1 je kvazinestacionárním signálem, který se po uplynutí přechodového děje stává signálem stacionárním. Je samozřejmé, že v praxi budou výstupy soustav zatíženy poruchami a výsledky nebudou tak jednoznačné. Přesto lze předpokládat, že uvedeným způsobem půjde odlišit výstupní signály a tím i soustavy Si1 a Sp1. 2) Regresní rovnice typu Sp1 a) Pro testování použijeme nejprve příslušný lineární regresní model bez absolutního členu, který přesně odpovídá oběma diskrétním modelům soustav DSi1 (a1=1) i DSp1 (a1<1) (43) y t = a1 y t −1 + b1u t −1 + ε t , kde ε t označuje rezidua. Výsledky regrese jsou uvedeny v tab. 18. Tab. 18. Výsledky regrese pro model bez absolutního členu (Statgraphics, EasyReg) Procesy a1 b1 R2 [%] DW JB Homoskedasticita Poznámka yI, ut-1 1.000 0.002 100 2.823 1.67 + perfektní model yP, ut-1 0.990 0.023 100 2.815 2.83 + perfektní model Závěr: Oba modely jsou statisticky významné, korektní (správné) a perfektní (R2 = 100 %). Podle velikosti koeficientu a1 lze jednoznačně posoudit příslušnost signálů k typům soustav (pro a1 = 1 jde o výstup z integrační soustavy). Tato úvaha výrazně koresponduje s myšlenkou testování tzv. jednotkových kořenů procesů DS/TS používanými v moderní analýze časových řad (v ekonometrii).

b) Pokud uskutečníme výpočet pro umělý (neodpovídající modelům dynamických soustav) lineární regresní model s absolutním členem, dostaneme výsledky uvedené v tab.19. Regresní model má v tomto případě tvar (44) y t = a 0 + a1 y t −1 + b1u t −1 + ε t . Tab. 19. Výsledky regrese pro model s absolutním členem (Statgraphics, EasyReg) a1 b1 R2 [%] DW Dh JB BP Poznámka Procesy a0 perfektní model yI, ut-1 (0.000) 1.000 0.002 100 2.823 (-4.10) 1.68 0.60 perfektní model yP, ut-1 (0.000) 0.990 0.023 100 2.813 (-4.02) 3.05 1.60 Pozn.: Oba modely jsou statisticky významné, korektní (správné) a perfektní (R2 = 100 %). Durbinova h statistika (Dh) sice vykazuje statisticky významnou hodnotu (vzhledem k N(0,1)) a signalizuje zápornou autokorelaci reziduí, avšak vzhledem ke statisticky nevýznamné hodnotě DW a ke korektně definovanému modelu je tento údaj nevhodný a matoucí. Absolutní člen regresní rovnice a0 je v obou případech nulový a statisticky nevýznamný. Závěr: Podle velikosti koeficientu a1 lze opět jednoznačně posoudit příslušnost signálů k typům soustav. Uvedené skutečnosti můžeme jednoduše podchytit (pro uvažovaný teoretický případ) v algoritmu identifikace typu soustav a jejich výstupních procesů (obr. 17). TS

Z yt = (a0) + a1yt-1 + b1ut-1 + ε

Si1

-

a1 < 1 ?

+

Sp1

K

Obr. 17. Algoritmus testování typů soustav (a jejich výstupních procesů) Je logické, že v praktických případech mohou jednoznačnost identifikace ztížit jak aditivní šum na vstupech/výstupech soustav, tak mimořádně velké časové konstanty T1 soustav Sp1. Závěr: Přehledné výsledky možnosti a úspěšnosti identifikace soustav Si1 a Sp1 na základě jejich vstupně/výstupních signálů (proměnných, veličin, procesů) jsou uvedeny v tab. 20. Tab. 20. Možnosti identifikace soustav Si1 a Sp1 na základě jejich vstupů a výstupů Metoda Data / model Rovnice Závěry KAČR původní, diferencované, nelze identifikovat, rozlišit! trendově očištěné MAČR původní, diferencované, nelze identifikovat, rozlišit! trendově očištěné ∆y t = b1u t −1 + ε t OK: yI, R2 = 100% DSi1 ∆y t = a 0 + b1u t −1 + ε t OK: yI, a0 = 0, R2 = 100% RA y t = a1 y t −1 + b1u t −1 + ε t OK: yI, a1 = 1 Dsp1 yt = a 0 + a1 yt −1 + b1u t −1 + ε t OK: yI, a1 = 1

Pozn.: KAČR ... klasická analýza časových řad, MAČR ... moderní analýza časových řad, RA ... regresní analýza. Shrnutí: Z předchozích dílčích závěrů shrnutých do tabulky lze konstatovat, že nejspolehlivější identifikační metodou soustav Si1 a Sp1 je statistická identifikace pomocí regresní analýzy s vhodnými modely. Identifikace pomocí klasické a moderní analýzy časových řad se ukázala jako nevhodná a nepoužitelná, což je v rozporu s předchozími výsledky, kdy byl tento způsob identifikace nestacionárních procesů typu DS/TS naopak velice vhodný. Toto konstatování jenom dále dotvrzuje závěr, že pro technometrickou analýzu signálů dynamických soustav technologických procesů nelze automaticky použít všechny výsledky teoretických přístupů ekonometrie a moderní i klasické analýzy časových řad.

6 Závěr Náplní výzkumu v oblasti řízení technologických procesů v hutnictví je [MORÁVKA, J. 1999a], [MORÁVKA, J. 1999b], [MORÁVKA, J. 2000]: - nejprve analýza (deterministická i statistická) struktury, kauzality a závislostí (vazeb, interakcí) mezi veličinami, - pak jejich (sub)optimální řízení. Cíle optimalizace technologických procesů jsou výrobně-ekonomické, týkají se procesu výroby i samotných výrobků a patří mezi ně zejména: - minimalizace poruch, nestacionarit technologického procesu (např. zastavení aglomeračního pásu, průvalů na ZPO apod.) - maximalizace životnosti technologických agregátů / minimalizace jejich opotřebení (např. Cu vložek krystalizátorů na ZPO apod.) - maximalizace (optimalizace) kvality výrobků (např. aglomerátu, předlitků apod.) - minimalizace vad (reklamací) výrobků - maximalizace výrobnosti zařízení - minimalizace nákladů na výrobky (zvýšení konkurenční schopnosti produktů na trhu). Technologické procesy (např. v hutnictví) jsou přitom obecně dynamické systémy vícerozměrné, nelineární, stochastické, spojitě-diskrétní s rozloženými parametry – typu MIMO (Multiple Input Multiple Output). Pokud obsahují integrační subsystémy a proporcionální setrvačné subsystémy s velkými časovými konstantami (což je v hutnických systémech velice běžné), pak se v systémech vyskytují nestacionární vstupní, stavové a výstupní procesy (typu DS, ale i TS). Klasifikace, identifikace a statistická analýza jednotlivých (jednorozměrových) nestacionárních procesů je už poměrně dobře zvládnutá. Analýza systémů nestacionárních procesů pomocí P, V a H kanonických struktur víceparametrových systémů [ŠULC, B. 1999] a aplikace těchto poznatků je však teprve ve stadiu zkoumání a rozpracování. Je tedy celkem logické, že uvedené náročné úkoly výzkumu řízení technologie se neobejdou bez znalostí a použití moderních metod statistické analýzy (i) nestacionárních náhodných procesů v mnohorozměrových dynamických systémech technologických procesů.

7 Literatura ANDĚL, J. 1985. Matematická statistika. 2.vyd. Praha : SNTL/ALFA, 1985. 352 s. ARLT, J. 1999. Moderní metody modelování ekonomických časových řad. 1.vyd. Praha: Grada Publishing, s.r.o., 1999. 312 s. ISBN 80-7169-539-4. BAKYTOVÁ, H. AJ. 1979. Základy štatistiky. 2.vyd. Bratislava : ALFA, 1979. 392 s.

BENKWITZ, A. AJ. 1999. Multiple Time Series Analysis, Co-Integration. [online]. Humboldth University, Berlin, Germany, 1999. Dostupné z 6 s. BIERENS, H. 1999a. Cointegration Analysis [online]. Pensylvania State University, PA, 1999. Dostupné z 29 s. BIERENS, H. 1999b. Nonparametric Nonlinear Co-Trending Analysis, With an Application to Interest and Inflation in the U.S. [online]. Pensylvania State University, PA, 1999. Dostupné z 39 s. BIERENS, H. 1999c. Free Econometrics Software for Easy Regression Analysis [online]. Pensylvania State University, PA, 1999. Dostupné z . CIPRA, T. 1986. Analýza časových řad s aplikacemi v ekonomii. 1.vyd. Praha : SNTL/ALFA, 1986. 248 s. CYHELSKÝ, L. AJ. 1986. Teorie statistiky. 2.upravené vyd. Praha : SNTL/ALFA, 1986. 344 s. GARAJ, V. & ŠUJAN, I. 1980. Ekonometria. 1.vyd. Bratislava : ALFA/SNTL, 1980. 288 s. HEBÁK, P. & HUSTOPECKÝ, J. 1990. Průvodce moderními statistickými metodami. 1.vyd. Praha : SNTL, 1990. 296 s. ISBN 80-03-00534-5. JOHANSEN, S. 1997. Mathematical and Statistical Modelling of Cointegration. In Sborník EUI (European University Institute) Florence, Italy : Working Paper ECO No. 97/14, 1997, 16 s. KAŇOKOVÁ, J. 1989. Teorie statistiky pro řízení a plánování. 1.vyd. Praha/Bratislava: SNTL/ALFA, 1989. 398 s. KOSCHIN, F. AJ. 1992. STATGRAPHICS aneb statistika pro každého. 1.vyd. Praha : Grada, 1992. 349 s. ISBN 80-85424-70-3. LJUNG, L. 1987. System Identification : Theory for the User. 1. vyd. Englewood Cliffs New Jersey : PTR Prentice Hall, 1987. 519 s. ISBN 0-13-881640-9. MELOUN, M. & MILITKÝ, J. 1994. Statistické zpracování experimentálních dat. 1.vyd. Praha : PLUS, 1994. 839 s. ISBN 80-85297-56-6. MORÁVKA, J. 1999a. Ekologická optimalizace provozu A2 SP4. (Analytická studie projektu č. 5098030). Třinec : TŽi a.s. SA - Středisko automatizace, červen-říjen 1999. 40 s. MORÁVKA, J. 1999b. Hierarchický distribuovaný systém řízení aglomeračního procesu. Disertační práce. Ostrava : KATŘ FS VŠB-TU Ostrava, září 1999. 170 s. MORÁVKA, J. 2000. Základní rozbor možností statistického zpracování technologických dat ZPO 1. Úvodní studie 1.etapy projektu č.5100004 “Statistické zpracování technologických dat ZPO 1”. Třinec : SPPČ TŽi, a.s., březen 2000. 18 s. SEGER, J., HINDLS, R. & HRONOVÁ, S. 1998. Statistika v hospodářství. 1.vyd. Praha : ETC, 1998. 636 s. ISBN 80-86006-56-5. SWOBODA, H. 1977. Moderní statistika. I.vyd. Praha : Svoboda, 1977. 352 s. ŠULC, B. 1999. Teorie automatického řízení s počítačovou podporou. 1.vyd. Praha : skripta FS ČVUT Praha, 1999. 154 s. ISBN 80-01-01974-8. TŮMA, J. 1998. Složité systémy řízení I. – Regulace soustav s náhodnými poruchami. 1.vyd. Ostrava : skripta FS VŠB-TU Ostrava, 1998. 158 s. VÍTEČEK, A. 1988. Matematické metody automatického řízení (Transformace L a Z). 1.vyd. Ostrava : skripta FSE VŠB Ostrava, 1988. 156 s.