© Institut biostatistiky a analýz

SPEKTRÁLNÍ ANALÝZA ČASOVÝCH ŘAD

SPEKTRÁLNÍ ANALÝZA ČASOVÝCH ŘAD

prof. Ing. Jiří Holčík, CSc.prof. Ing. Jiří Holčík, CSc.

holcikholcik@iba.muni.cz@iba.muni.cz

© Institut biostatistiky a analýz

V.PARAMETRICKÉ METODY ODHADU VÝKONOVÉHO

SPEKTRApokračování

V.PARAMETRICKÉ METODY ODHADU VÝKONOVÉHO

SPEKTRApokračování

MA A ARMA MODELYMA A ARMA MODELY

© Institut biostatistiky a analýz

MA MODELYMA MODELY

opakování

0m)mT(

qm0bb

qm0

)mT(

yy

mq

0kmkk

2xyy

jindy0

0m}xx{E0}x{E

xby

2x

nmnn

q

0mmnmn

© Institut biostatistiky a analýz

MA MODELYMA MODELY

(b0z0+b1z1+…+bqzq).(b0z0+b1z-1+…+bqz-q) =

= b02+b0b1z+…+b0bqzq+b0b1z-1+b1

2+b2b1z+…+bqb1zq-1+…+

+b0bqz-q+b1bqz-q+1+b2bqz-q+2+…+bq2 =

= (b0b0+b1b1+…+bqbq).z0 +(b0b1+b1b2+…+bq-1bq).z-1 +

d0 d-1

+(b1b0+b2b1+…+bqbq-1).z1 +…+(…).zq

d-1

q

qm

mm

1 zd)z(D)z(B).z(B

mq

ékmkkm qmprobbd

© Institut biostatistiky a analýz

MA MODELYMA MODELY

Tedy

a výkonové spektrum

odhad spektra:

qmd

qm0)mT(

m2x

yy metoda momentůmetoda momentů

q

qm

fmT2jyy

MA e).mT()f(yy

q

qm

fmT2jyy

MA e).mT(r)f(Pyy

!!! A TO JE KLASIKA !!!!!! A TO JE KLASIKA !!!

© Institut biostatistiky a analýz

MA MODELYMA MODELY

ALTERNATIVA:stanovení {bk}založené na aproximaci MA procesu AR

procesem vysokého řádu,tj. p » q

v tom případě platí, že B(z) = 1/A(z), resp. B(z).A(z)=1 a proto

q

1kknkn 0n0

0n1a.ba

© Institut biostatistiky a analýz

MA MODELYMA MODELY

ALTERNATIVA:protože p » q, můžeme odhad {bk} zpřesnit pomocí metody nejmenších

čtvercůurčíme kvadratickou chybu

a tu minimalizujeme výběrem parametrů {bk}

kde;q...,,2,1j,ipro,a.a)ji(R

jip

0njinnaa

0kpro0a;1a,1a.bae k0

p

0n

2q

0kknkn

aa1aa.rRb

.q...,,2,1ipro,a.a)i(rip

0ninnaa

vymyslel to Durbin a ukázalo se, že je to přibližně odhad s maximální pravděpodobností, je-li proces normální

vymyslel to Durbin a ukázalo se, že je to přibližně odhad s maximální pravděpodobností, je-li proces normální

© Institut biostatistiky a analýz

URČENÍ ŘÁDU MA MODELUURČENÍ ŘÁDU MA MODELU

pro vyjádření úzkých pásem je třeba hodně vysoký řád;

intuitivně:

hodnoty odhadu autokorelační funkce musí rychle klesat k nule, protože yy(mT)=0 pro m>q

test, zda ryy(mT)0 je založený na srovnávání ryy(qT) s rozptylem hodnot ryy(qT) pro m<q

© Institut biostatistiky a analýz

jinak:založeno na testu „bělosti“ posloupnosti, která je

výsledkem působení soustavy inverzní k odhadnutému MA modelu na analyzovanou posloupnost

Akaikovo informační kritérium

URČENÍ ŘÁDU MA MODELUURČENÍ ŘÁDU MA MODELU

N

q2lnAIC 2

wq

© Institut biostatistiky a analýz

ARMA MODELYARMA MODELY

opět si stručně zopakujeme:

0k0

0k)kT(a1b

)T)kn((xb)T)kn((ya)nT(y

2x

xx0

q

0kk

p

1kk

qm0)mTkT(hb)kTmT(a

qm)kTmT(a)mT( q

mkk

2x

p

1kyyk

p

1kyyk

yy

2p

1kk

2q

1kk

2x

ARMAyy

)fkT2jexp(a1

)fkT2jexp(b1.T

)f(P

odhadvýkonovéhospektra

odhadvýkonovéhospektra

© Institut biostatistiky a analýz

ARMA MODELYARMA MODELY

ARMA(1,1):

u(kT) je jednotkový skok

mnoho teoretických postupů bez praktického významu, protože: jsou výpočetně příliš náročné (maticové operace, iterační

optimalizační postupy); iterační optimalizace nezaručuje konvergenci řešení nebo

konvergenci ke skutečně optimálnímu řešení;

proto suboptimální postupy : používají metody nejmenších čtverců řešení soustavy lineárních

rovnic; oddělený výpočet AR a MA parametrů

2m),TmT(a)mT(

b)0(a)T(

)bab1()T(a)0(

)1k(u)a(b)kT(u.)a()kT(h

yy1yy

12xyy1yy

1121

2xyy1yy

1k11

k1

© Institut biostatistiky a analýz

ARMA MODELYARMA MODELYMETODA PRVNÍMETODA PRVNÍ

pro m>q je

dosadíme yy a řešíme p lineárních rovnic (pro modely vyšších řádů horší výsledky díky slabším odhadům yy(mT) pro větší m);

přeurčená soustava lineárních rovnic (pro m>q)

její řešení opět optimalizací metodou nejmenších čtverců

)a(f.resp),a(f)mT( kkyy

© Institut biostatistiky a analýz

ARMA MODELYARMA MODELYMETODA PRVNÍMETODA PRVNÍ

předpokládejme, že známe odhady autokorelační funkce až do zpoždění M > p+q

protože Ryy je rozměru (M-q) x p a M-q > p, lze použít metodu nejmenších ;

výsledkem minimalizace je

(může být použito i váhování k potlačení méně spolehlivých odhadů AK funkcí)

yyyy

yy

yy

yy

p

2

1

yyyyyy

yyyyyy

yyyyyy

)M(r

)2q(r

)1q(r

a

a

a

.

)pM(r)2M(r)1M(r

)2pq(r)q(r)1q(r

)1pq(r)q(r)q(r

raR

yyTyy

1yy

Tyy ..).( rRRRa

© Institut biostatistiky a analýz

ARMA MODELYARMA MODELYMETODA PRVNÍMETODA PRVNÍ

analyzovanou sekvenci vyfiltrujeme FIR filtrem a dostaneme

kaskádní zapojení ARMA(p,q) [HARMA(z)=B(z)/A(z)] s AR(p) reprezentuje přibližně MA(q) proces s HMA(z)=B(z)

p

1k

kk z.a1)z(A

)z(A

;1N...,,1,0n,)kTnT(y.a)nT(y)nT(vp

1kk

)z(A)z(HAR

© Institut biostatistiky a analýz

z filtrované sekvence v(nT) pro pnN-1 je rvv(mT), odhad výkonového spektra potom je

(opět můžeme „woknovat“ k potlačení méně spolehlivých odhadů autokorelačních funkcí nebo filtrace AR(q) oběma směry dvě různé posloupnosti autokorelační funkce)

ARMA MODELYARMA MODELYMETODA PRVNÍMETODA PRVNÍ

q

qmvv

MAyy )fmT2jexp()mT(r)f(P

2p

1kk

MAvvARMA

yy

)fkT2jexp(a1

)f(P)f(P

© Institut biostatistiky a analýz

vstupní/výstupní identifikace metodou nejmenších

opakování:

nelinearita vztahu pro výpočet yy(mT) plyne z toho, že neznáme yy(mT-kT), protože neznáme x(nT)

ARMA MODELYARMA MODELYMETODA DRUHÁMETODA DRUHÁ

q

mkxyk

p

1kyykyy )kTmT(b)kTmT(a)mT(

q

0kk

p

1kk )kTnT(xb)kTnT(ya)nT(y

© Institut biostatistiky a analýz

y = H.θ +

y = [y(0) y(T) … y(NT-T)]T

= [x(0) x(T) … x(NT-T)]T

θ = [-a1 –a2 … -ap b1 b2 … bq ]T

matice vstupních/výstupních dat o rozměru N x (p+q)

ARMA MODELYARMA MODELYMETODA DRUHÁMETODA DRUHÁ

1N...,,1,0n),nT(x)kTnT(xb)kTnT(ya)nT(yq

1kk

p

1kk

1qN3N2N1pN3N2N

1q101p10

q21p21

xxxyyy

xxxyyy

xxxyyy

H

© Institut biostatistiky a analýz

minimalizace nejmenšími

počáteční podmínky y-p, …, y-1, x-q, …, x-1 buď specifikovat nebo nulové

odhad hodnot vstupního šumového signálu z analyzované posloupnosti AR modelem vysokého řádu

ARMA MODELYARMA MODELYMETODA DRUHÁMETODA DRUHÁ

yHHHΘ ..).( T1T

© Institut biostatistiky a analýz

(už tu částečně byla)

koeficienty {ck} se určí nějakým AR algoritmem a z nich {ak} a {bk}

je-li p>q, pak

nebo v časové oblasti

ARMA MODELYARMA MODELYMETODA TŘETÍMETODA TŘETÍ

1k

kkzc1)z(Ckde,

)z(C

1

)z(A

)z(B

qpM,zc1)z(CM

1k

kk

,)z(C

1

)z(A

)z(B

q

0k0nknk 1b,...;2,1n,ac.b

© Institut biostatistiky a analýz

protože by mělo platit an = 0 pro n>p, je

pro určení se spočítají z

ARMA MODELYARMA MODELYMETODA TŘETÍMETODA TŘETÍ

qp,...,2p,1pnpro,0c.b kn

q

0kk

qp

2p

1p

q

2

1

p2qp1qp

q2pp1p

q1p1pp

c

c

c

b

b

b

.

ccc

ccc

ccc

}b,...,b,b{ q21

}a,...,a,a{ p21

1c;p,...,2,1npro,ac.b 0nkn

q

0kk

© Institut biostatistiky a analýz

co je maticově

ARMA MODELYARMA MODELYMETODA TŘETÍMETODA TŘETÍ

q

1

qp2p1pp

12

1

p

2

1

b

b

1

c

0

0

ccc

1cc

01c

a

a

a

.

© Institut biostatistiky a analýz

URČENÍ ŘÁDŮ ARMA MODELUURČENÍ ŘÁDŮ ARMA MODELU

dodatečné kritérium:posouzení bělosti posloupnosti po inverzní filtraci analyzované

posloupnosti navrženou ARMA soustavou

odhad p

det(R’yy)=0, pokud je rozměr modelu větší než řád analyzovaného procesu

je odhad rozptylu vstupní chybové posloupnosti

wpqwpq ,N

)qp(2ln

q)AIC(p,

)q(r)2pq(r)1pq(r

)2pq(r)q(r)1q(r

)1pq(r)1q(r)q(r

'R

yyyyyy

yyyyyy

yyyyyy

yy

rozšířené modifikované Y.-W. rovnice