Microsoft PowerPoint - 時系列解析（１１）_講義用.pptx

東京⼤学 数理・情報教育研究センター

北川 源四郎

時系列解析（11）

1. 分散⾮定常モデル： 線形化・正規近似

2. 共分散⾮定常モデル：時変係数モデル

3. ⾮線形・⾮ガウス型状態空間モデル

分散・共分散⾮定常

⽇経225

地震波

⾮定常時系列のモデル

１．平均⾮定常 ・・・ トレンド，季節調整

２．分散⾮定常

３．共分散⾮定常

• 時変係数モデル

• 線形・ガウスモデル（カルマンフィルタ）で推定

するためには、近似が必要（線形化＋正規近似）

• ⾮線形・⾮ガウス型状態空間モデルを使うと直接

的なモデリングが可能

2

2

2

_log

_log

log

r

_n





_n



w

_n

n

r

2

log

r

_n

)

1

,

0

(

~

,

w

N

w

r

_n





_n

_n

_n

2

2

2

n

n

n

w

r





（線形モデルによる）時変分散の推定

⼆乗

対数

2

2

1

2

2

1

log

log

log

log

n

n

n

n

n

n









 













 



（線形モデルによる）時変分散の推定

2

2

1

2

2

2

log

log

log

log

log

n

n

n

n

n

n

r

w



















時変分散のモデル（例）

分散変化のモデル（例）

ランダムウォーク型

（定数付き）AR型

状態空間表現

2

2

1

2

2

2

log

log

log

log

log

n

n

n

n

n

n

r

w



















2

2

2

log

,

log ,

log

n

n

n

n

n

n

x





y



r





w

とおくと

1

n

n

n

n

n

n

x

x

y

x















• 状態空間表現

•



は正規分布には従わない

⼆種類のデータ変換

2

log

n n

z



r

n

r

2 n n

y



r

y

2

n



(

r

2 12n



r

22n

) / 2

2

_n

log 2

_n

z



y

の分布：⼆重指数分布

2

2

1

1

~

1

(

)

x



y



カイ⼆乗

分布

~

(0,1)

n

y

N

1

log

1

w



x





1

2

( )

h

h

x

h y

y

y

x









 





 

2 1 1 1 2 1 _{2 2} 1 1 1 2 ₂ 1 2

1

exp

2

2

2

1

( )

2

1

2

( )

( )

( )

( | 0,1)

2

2

exp

2

2

exp

dh

dx

x

y

x

h y

y

f y

h

x

x

x

dh

x

g x

x

x

dx

x









 _   _  















_

_



















 

_



_







1 1

log

( )

d

(

( ))

p w

g h

w

dw

 







1

2

 

2

( )

2

exp

x

g x





x





2

log

w

_n

⼆重指数分布（Gumbel分布）

2

2

2

2

2

1

2

2

1

χ

2

=

(

) / 2 ~

(

⾃由度２の 分布 指数分布

)

s



y



y



~

(0,1)

n

y

N

2

log

2

w



s





( )

exp

w

:

⼆重指数分布

p w



w



e

 

( )

exp

g s





s

 





1

( )

(

( ))

exp

exp

w

w

w

w

de

p w

g h

w

dw

e

e

w

e













1

1

( )

log

( )

w

h

h

w

h s

s

s

w

h

w

e













 Euler定数

時変分散・ボラティリティ

Kitagawa & Gersch (1985)

⼆重指数分布の正規分布による近似















2

2

exp

2

1

~

log

2

w

n

e

w

w



2

( ,

/ 2),

1.27036

N

 



 





2

log

~ exp

w

n

w

w



e

N

( ,

 

2

/ 6),



 

0.577216

log-ピリオドグラムの平滑化

Wahba (1980)

⼆重指数分布と正規近似

True

近似

1

2

2

( )

exp

w

w

e

g w











_

_









( )

exp

w

g w



w



e

2

log(

y

_n

)

2 2

log



2 2 2 1 2

log(

y

_m



y

_m

)

•

正規近似が良くなる

•

分散が1/3

•

データ数が1/2

近似

_True

2 1

log



分散変動(時変分散)モデル

k

m

m

m

m

m

t

v

z

t

w













2

6

( )

exp

w

(

,

)

h w



w



e



N











2

1

2

2

2

( )

exp

(

w

) ~

(

,

)

h w

w

e

N









１





2

2

-

型



2

1

-

型



時変分散の推定： MYE１F

型，トレンド次数＝2



2

_{= 6.6x10}

-6

_，



2

_{= 9.71ｘ10}

-1

log-lkhd = -2195.731

AIC=4399.46

-1 0 1 2 3 4 5 6 1 101 201 301 401 501 601 701 801 901 1001 1101 1201 0 200 400 600 800 1000 1200

2

2

-

# an earthquake wave data

data(MYE1F)

#

tvvar(MYE1F, 2, 6.6e-06, 1.0e-06)

tau2

6.60000e-06

sigma2

9.71228e-01

log-likelihood -2195.731

aic

4399.462

時変分散の推定とデータの等分散化

変換したデータ

時変分散

等分散化したデータ

Time 0 500 1000 1500 2000 2500 -4 0 -2 0 0 2 0 4 0

AR型の分散変動モデル

2

2

2

2

1

2

log

log

log

log

log

n

n

n

n

n

n

w

r

v









_







2

2

1

2

2

2

log

log

log

log

log

n

n

n

n

n

n

r

w



 



















ランダムウォークモデル

AR 型モデル

時変スペクトル

ARモデル

⾃⼰共分散

スペクトル

時変係数

⾮定常

時変

)

,

0

(

~

,

2

1



N

w

w

y

a

y

_n

_j

_n

_n

m

j

jn

n











時変係数 AR モデル

時変係数AR モデルの推定

1

,

は時間とともに変動

m

n

jn

n

j

n

nj

j

y

a y

_

w

a









2

,

~ (0, )

k

jn

jn

jn

a

v

v

N







1

n

n

n

n

n

n

n

x

Fx

Gv

y

H x

w











状態空間表現

係数変化のモデル

~

(0, )

~

(0, )

n

n

v

N

Q

w

N

R

, , ,

, ,

n

n

x F G H

Q R

を定める

時変係数ARモデル



Kronecker 積

1

)

1

(

)

1

(

)

1

(

_

_

_

H

G

F



1

0



,

0

1

,

0

1

1

2

)

1

(

)

1

(

)

2

(

































H

G

F





































































B

a

B

a

B

a

B

a

b

b

b

b

a

a

a

a

B

A

m m pq p q m m    































1 1 11 1 1 11 1 1 11

( )

( )

( )

1

2

2

1

1

,

(

,

,

)

,

(

,

,

)

( , ,

,

)

k

k

m

m

k

n

n

n m

m

T

k

n

n

mn

F

F

I

G

G

I

H

H

y

y

Q

I

R

x

a

a

I B

B







































状態空間表現（k = 1 の場合）





1

1,

1

1,

,

1

,

1

1

2

2

2

1

1

1

1

,

,

,

n

n

n

mn

m n

m n

n

n

n

n m

n

mn

a

a

v

a

a

v

a

y

y

y

w

a

Q

R

























 

















 

_



_















 

















 









 



_

_





_

_











_

_





















_

_































1

1

1

1,

1

1

1

1

2

,

1

2

1

1

1

1

2

1

1

2

1

1

1

1

1

,

,

, 0,

, 0

n

n

n

m n

m n

n

n

m n

m n

m n

n

mn

n

n

n m

n

a

a

v

a

a

a

a

v

a

a

a

a

y

y

y

a



































































_{ }

_







_











_

_



















_{ }

_















_{ }

_

















































_

_





_

_



















_



























n

w

























状態空間表現（k = 2 の場合

）

2

2

2

,

Q

R

















_

_













システムノイズの等分散仮定について

Q = diag{



2

,…,



2

}: の仮定は妥当か？

1 2 1 2

2

1

2

2

2

1

1

1

2

2

2

2

1

2

1

1

2

2

( , ) 1

,

AR

( )

( , )

( , )

( , )

|

( , ) |

(

)

,

~ (0, )

係数のフーリエ変換

の時間変化の滑らかさに関する評価

周波数領域-

で⼆乗積分

各次数を同じ割合で加算しているので

m

ijf

nj

j

n

m

k

k

ijf

nj

j

m

k

k

nj

j

k

n

nj

nj

nj

A f n

a e

f

p

f

A f n

A f n

A f n

a e

f

A f n

df

a

a

v

v

N























 

 









 



















ARモデルを⽩⾊化フィルタと

考えたときの周波数応答関数

時変係数ARモデルのAIC

m

k=1

k=2

m

k=1

k=2

1

6492.5

6520.4

6

4831.9

4873.8

2

5527.7

5643.2

7

4821.6

4878.7

3

5070.0

5134.5

8

4805.1

4866.9

4

4820.0

4853.0

9

4813.4

4884.9

5

4846.0

4886.0

10

4827.1

4911.9

data(MYE1F) # an earthquake wave data

z <-

tvar(MYE1F, trend.order = 2, ar.order = 8,

span = 20, tau2.ini = 6.6e-06, delta = 1.0e-06)

z

tau2

1.60000e-06

sigma2

1.43071e+01

log-likelihood -7284.520

aic

14589.041

Rによる時変係数ARモデルの推定

時変スペクトル

)

,

0

(

~

,

2

1



N

w

w

y

a

y

_n

_j

_n

_n

m

j

jn

n











2 1

2

2

1

( )

m j

n

n

ijf

jn

p

f

e

a















時変係数 AR モデル

時変スペクトル

Rによる時変スペクトルの計算

z <- tvar(MYE1F, trend.order = 2, ar.order = 8,

span = 20, tau2.ini = 6.6e-06, delta = 1.0e-06)

# 時変スペクトル

spec <- tvspc(z$arcoef, z$sigma2)

係数の急激な変化について

• トレンドモデルによる変化はゆっくりした変化を仮定している

• 地震波などでは突然別のモデルに変化することがある。

• 対応１

 変化点既知の場合

k=1

の場合： その時点で



2

を⼤きくすればよい

k=2

の場合： それだけでは屈折点となる

• x

_n|n-1

とV

_n|n-1

を初期化するかV

_n|n-1

の対⾓成分に

⼤きな値を⼊れる

• 対応２： ⾮ガウス型モデルを利⽤する（変化点未知でよい）

z <- tvar(MYE1F, trend.order = 2, ar.order = 8,

span = 40, outlier = c(630, 1026), tau2.ini = 6.6e-06,

delta = 1.0e-06)

Rによる時変係数ARモデルの推定（構造変化を仮定）

# n=630と1026の2か所で構造変化があったと仮定

#

z <- tvar(MYE1F, trend.order = 2, ar.order = 8, span = 40,

outlier = c(630, 1026), tau2.ini = 6.6e-06, delta = 1.0e-06)

#

# 時変係数ARモデルから時変スペクトルを計算

# 時変スペクトルを3次元表⽰

#

spec <- tvspc(z$arcoef, z$sigma2)

plot(spec, theta = 30, phi = 40, expand = 0.5)

時変係数と時変スペクトル

時変係数AR

局所定常構造

時変係数AR

本と同様の3次元プロット（未公開）

########################################## # 時変スペクトルの3次元表⽰ ########################################## # seismic data data(MYE1F)

z <- tvar(MYE1F, trend.order = 2, ar.order = 8, span = 20, outlier = c(630, 1026), tau2.ini = 6.6e-06, delta = 1.0e-06) spec <-tvspc(z$arcoef, z$sigma2) ######################################### # 最初のスペクトル(n=0) nf <- 201 dt <- 2 dy <- 0.2 nf1 <- nf-dt t <- 1:nf tt <- 1:nf1 plot(t,spec$z[,1],type="l", xlim=c(1,501),ylim=c(-2,18)) y <- spec$z[,1] z <- y # ⿃瞰図(n=1,80) nrep <- 1:80 for (i in nrep){ par(new=T) t <- t+dt for (j in 1:nf) z[j] <- spec$z[j,i]+dy*(i-1) for (j in tt){ # z[j] <- max(spec$z[j,i]+dy*(i-1),y[j+dt]) z[j] <- max(z[j],y[j+dt]) } y <- z plot(t,y,type="l", xlim=c(1,501),ylim=c(-2,18),xaxt='n',yaxt="n")

スペクトルの滑らかさの制約

1 2 1 2

2

2

2

2

1

( , )

(2 )

k k

m

k

k

k

nj

j

A f n

f

R

df



j a

















2

4 2

1

0

2

,

~ (0,(

) )

jn

jn

jn

a



u

u

N





j





( )

k

_,

( )

k

_{~ (0,(}

k

_{) )}

2

jn

jn

jn

k

a



u

u

N

j





1

1, 1

1

, 1

1

1

n

n

n

mn

m n

mn

a

a

u

a

a

u















 













 

_



_{  }







 













 





 

_{ }

_{ }

_





  













2 1 2

0

~

,

0

1

0

n mn

v

N

v

w







_

_



 

 



_

_



 

 



_

_



 

_

 

_





 

 

 

_

_{  }

_

_









 

 







東京⼤学 数理・情報教育研究センター

北川 源四郎

時系列解析（１２）

• 拡張カルマンフィルタ

• ガウス和フィルタ

• ⾮ガウス型フィルタ

• 粒⼦フィルタ

• アンサンブルカルマンフィルタ

⾮線形・⾮ガウス型フィルタ

⾮ガウス型モデリングの必要性

構造変化

異常値（外れ値）

離散過程

⾮線形性

n

n

n

f

x

v

x



(

_

₁

)



0 0.1 0.2 0.3 0.4

⾮対称分布

状態空間モデルの拡張

⾮線形・⾮ガウス型

線形・ガウス型

n

n

n

n

n

n

w

Hx

y

Gv

Fx

x









_1

x

f x

v

y

h x

w

n

n

n

n

n

n







(

,

)

(

,

)

1

関数：⾮線形

分布：⾮ガウス型

⾮線形・⾮ガウス型状態空間モデル

1

(

,

)

(

)

n

n

n

n

n

n

x

f x

v

y

h x

w









:

:

:

n

n

n

n

x

y

v

w

: 状 態

時 系 列

シ ス テ ム ノ イ ズ

観 測 ノ イ ズ

~ ( )

~ ( )

n

n

v

q v

w

r w

(

,

)

n

n

n

y



h x

w

観測モデルの拡張

ただし，h(x,w) の逆関数 g(y,x) が存在して

(

,

)

n

n

n

w

g y

x

y

g

y



 





で微分可能（

/ が存在 ）

n n

x

n

n

x

n

n

y

e

w

w

e

y

g











_

(例)

⾮線形・⾮ガウス型状態空間モデル

⾮ガウス型分布の例

⾮線形関数の例

1

(

,

)

(

,

)

n

n

n

n

n

n

x

f x

v

y

h x

w







f(x)

⾮線形変換

f(x,v)

積型

h(x,w)

積型

（e

x

w

など

）

y=h(x,w)

から

w=k(y,x)

と

書ける必要がある。

h(x)+w

の⽅が簡単

・

_{コーシー分布}

0 0.7 1 201 2 2

1

(

)

( )

x

p x















・

_{ピアソン分布族}



2 2



( )

(

)

b

p x



C

x











  



2 1 1 1 2 2

( )

b

b

C

b





_









・

⼆重指数分布

・混合分布

1 1 2 2

( )

( , ) (1

) ( , )

p x





x



 

 

x



0.2 0.4 0.6 0.8

状態推定

⾮線形

予測

フィルタ

)

|

(

x

_n

Y

_n

_

₁

p

)

|

(

x

_n

Y

_n

p

1

|

1

|

n



,

n

n



n

V

x

n

n

n

n

V

x

_|

,

_|

カルマンフィルタ

初期値

予測

フィルタ

y

_n

n

 1

⾮ガウス型予測の導出

1

1

1

1

)

(

,

|

)

|

(



_

_

_











n

n

n

n

n

n

Y

p

x

x

Y

dx

x

p

1

1

1

1

1

1

1

1

1

( ,

|

)

( |

,

) (

|

)

( |

) (

|

)

n

n

n

n

n

n

n

n

n

n

n

n

p x x

Y

p x

x

Y

p x

Y

p x

x

p x

Y























dy

y

x

p

x

p

(

)





(

,

)







)

|

(

)

(

)

,

(

x

y

p

y

p

x

y

p



)

|

(

)

,

|

(

x

_n

x

_n

_

₁

Y

_n

_

₁



p

x

_n

x

_n

_

₁

p

1

1

1

1

1

( |

_n

_n

)

( |

_n

_n

) (

_n

|

_n

)

_n

p x Y

_



p x

x

_

p x

_

Y

_

dx

_







⾮ガウス型フィルタの導出

1

1

{ ,

,

} {

,

}

n

n

n

n

Y



y



y



Y

_

y

)

|

(

)

,

|

(

y

_n

x

_n

Y

_n

₁

p

y

_n

x

_n

p

_



( , )

( | )

( )

( | ) ( )

( )

p x y

p x y

p y

p y x p x

p y





1

1

1

1

1

1

1

1

(

|

)

(

|

,

)

(

,

|

)

(

|

)

(

|

,

) (

|

)

(

|

)

(

|

) (

|

)

(

|

)

n

n

n

n

n

n

n

n

n

n

n

n

n

n

n

n

n

n

n

n

n

n

n

p x

Y

p x

Y

y

p y

x

Y

p y

Y

p y

x Y

p x

Y

p y

Y

p y

x

p x

Y

p y

Y

























⾮ガウス型平滑化の導出

1

1

1

1

1

1

1

1

1

1

1

1

1

1

( ,

|

)

(

|

) ( |

,

)

(

|

) ( |

, )

(

,

| )

(

|

)

(

| )

(

| , ) ( | )

(

|

)

(

| )

(

| ) ( | )

(

|

)

(

| )

n

n

N

n

N

n

n

N

n

N

n

n

n

n

n

n

n

N

n

n

n

n

n

n

n

n

N

n

n

n

n

n

n

n

N

n

n

p x x

Y

p x

Y

p x

x

Y

p x

Y

p x

x

Y

p x

x Y

p x

Y

p x

Y

p x

x Y p x Y

p x

Y

p x

Y

p x

x p x Y

p x

Y

p x

Y







































( , )

( ) ( | )

p x y



p y p x y

1

1

(

_n

|

_n

,

_N

)

(

_n

|

_n

, )

_n

p x

x

_

Y



p x

x

_

Y

( , )

( | )

( )

( | ) ( )

( )

p z x

p x z

p z

p z x p x

p z





1

1

(

_n

|

_n

, )

_n

(

_n

|

_n

)

p x

_

x Y



p x

_

x

1

1

1

1

1

1

(

|

)

( ,

|

)

(

|

) (

|

)

(

|

)

(

|

)

n

N

n

n

N

n

n

n

n

N

n

n

n

n

n

p x

Y

p x

x

Y

dx

p x

x

p x

Y

p x

Y

dx

p x

Y











_

_















⾮ガウス型フィルタ・平滑化

⼀期先予測

p x

_n

x p x

_n

_n

Y

_N

( | )



( | )





(



1

| ) (



1

| )

p x Y

p y

x

p x Y

p y Y

n

n

n

n

n

n

n

n

(

|

)

(

|

) (

|

)

(

|

)







1

1

p x Y

( |

_n

_n

_

)

p x x

( |

_n

_n

) (

p x

_n

|

Y

_n

)

dx

_n

















1

1

1

1

1

フィルタ

平滑化

分布の近似

True

正規近似

.

区分線形近似

階段関数

混合正規

0. 線形・正規モデル近似

カルマンフィルタ・平滑化

1. 正規分布近似

拡張カルマンフィルタ・平滑化

2. 区分線形（階段）近似

⾮ガウス型フィルタ・平滑化

3. 混合正規分布近似

ガウス和フィルタ・平滑化

3. 粒子近似

遂次モンテカルロフィルタ・平滑

化

拡張カルマンフィルタ

予 測

フィルタ

| | n n n n n n

n

n

n

n

n

n

x

x

x

x

f

F

x

g

G

x











 

_













 

_







| 1 n n n

n

n

n

_x

_x

h

H

x











 

_







| 1

1| 1

| 1

1| 1

(

)

n n

n

n

n

T

T

n n

n n

n

n

n

n

n

x

f

x

V

F V

F

G Q G



 



 







1

|

1

|

1

|

|

1

|

1

|

|

1

(

)

(

(

))

(

)

T

T

n

n n

n

n

n n

n

n

n n

n n

n

n

n

n n

n n

n

n

n n

K

V

H

H V

H

R

x

x

K

y

h x

V

I

K H V



























密度関数の数値近似

密度関数

近似

記号

p(x

_n

|Y

_n-1

) {d;t

₀

,…,t

_d

;p

₁

,…,p

_d

}

p(t)

p(x

_n

|Y

_n

) {d;t

₀

,…,t

_d

;f

₁

,…,f

_d

}

f(t)

p(x

_n

|Y

_N

) {d;t

₀

,…,t

_d

;s

₁

,…,s

_d

}

s(t)

q(v

_n

)

{2d+1;t

_-d

,…,t

_d

;q

_-d

,…,p

_d

}

q(t)

1次トレンドモデルの場合

⼀期先予測

n

n

n

n

n

n

w

t

y

v

t

t









_1

数値積分による実現

0 1

1

( )

(

) ( )

(

) ( )

d j j

t

i

i

i

t

d

_t

i

t

j

p

p t

q t

s f s ds

q t

s f s ds

















 

(

) ( )

( )

(

)

n

i

i

i

i

n

i

i

r y

t

p t

f

f t

C

r y

t

p

C











0 1

1

1

(

) ( )

(

) ( )

(

)

d j j

t

n

t

d

_t

n

t

j

d

n

j

j

C

r y

t p t dt

r y

t p t dt

t

r y

t p















 





 



フィルタ

数値積分による実現(平滑化）

平滑化

0 1

1

1

(

) ( )

( )

( )

( )

(

) ( )

( )

( )

( )

d j j

t

i

i

i

i

t

d

_t

i

i

t

j

d

i

j

j

i

j

j

q t

u s u

s

s t

f t

du

p u

q t

u s u

f t

du

p u

q

s

t f t

p



















  



 



各ステップの後で密度関数の全積分が１になるように規格化する