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

(1)

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

北川源四郎

時系列解析（１０）

(2)

2 季節調整とは

0.7

0.9

1.1

1.3

1.5

1

25

49

73

97

121

145 WHARD

東京⼤学北川源四郎数理⼿法VII （時系列解析）

月次データ

何らかの原因で特定の周期で繰り返す成分を

除去して本質的な現象を抽出する方法

(3)

季節調整モデル

y

_n

: 観測値

t

_n

: トレンド成分

s

_n

: 季節成分

w

_n

: 不規則成分（ノイズ）

n

t

s

w

y





非定常時系列をいくつかの成分に分解する方法

・

合成モデルによる成分分解

・本質的な変動を抽出する

観測モデル

(4)

4 トレンド成分モデル



k



n

k

n

k

n

t



0 

t



0 

t



v

1

2

1

2

3

1

2

3

3 n

n

k

t

v

k

t

v

k

t

v



















B





 1

多項式

局所的多項式

2 ,

~ (0,

)

k

n

t

v

N



 

状態空間表現

 

1

1 ,

1,

1

2

1

2 ,

,

1 0

1

0

0 n

n

k

x

t

F

G

H

t

k

x

F

G

H

t

_













 



_

_



_

_



_{ }







 





(5)

季節成分モデル

p

n

p

n

p

n

p

n

s

B s

s

B s









完全な周期

擬似的周期

2

2 (

)

2 n

n p

n

n p

n

p

n

s

u

s

u









2次の季節成分モデル

p:

周期

2.9 3 3.1 1 3 5 7 9 11131517192123252729313335

p

2 ,

~

(0,

)

(1

)

n

n p

n

p

n

S

u

N

B S

u













状態空間表現





1

0 ,

,

1

0 n

n

n p

S

x

F

G

S



 





_

_

_{ }





_

_

_{ }





_

_

_{ }



_

_







 





_

_

_{ }





_

_

_{ }









_

_

(6)

季節調整モデル

(1

)

(1

)

n

k

n

p

n

y

t

s

w

B t

v

B

s

u

 











実際には、このモデルはうまく機能しない

1－Bを共通の因子として含むため

1

1 

B

p

 

(1

B

)(1

 

B





B

p



)

(1



B e

)

_n

 　の任意の解

0 e

_n

=

c

n

t

c

s

c

  

  

y

n

  

t

n

s

n

w

n

  

t

n



s

n



w

n

₆

(7)

季節成分モデル（共通因子の除去）

1

1 (1

)(1

)

1

0

1

0 p

p

B







 



 













n

p

u

s

B









)

1 (



1 n

p

n

s

u

s





(

_

₁







_

_

₁

)



季節成分モデルの状態空間表現





1

0 ,

,

1 0

0

n n n

s

x



F

G

H











 









 









 



_

_







 





_

_

_{ }













(8)

季節調整モデル（基本型）

y

t

s

w

t

v

s

u

n

k

n

n p

n

  



 

_

 

_{ }





(

₁



₁

)

2

1

2

2 ~

(0,

)

~

(0,

)

~

(0,

)

n

w

N

v

N

u

N





₈

(9)

季節調整モデルの状態空間表現（合成モデル）

k=2, p=4の場合

x

F x

G v

y

H x

w

n







 1



x

_n



[

t

_n

,

t

_n

_

₁

|

s

_n

,

s

_n

_

₁

,

s

_n

_

₂

]

T



1

0

1

0

0 

0

1

0

1 ,

0

1

0

1

0

1

0

1

0

1

2 































H

G

F

(10)

季節調整法

Whard data

Trend

Seasonal

Noise

観測値＝トレンド＋季節成分＋偶然変動

₁₀

2

3

2

1

2

7

2 2,

12 ˆ

0.156 10

ˆ

_0.0248

ˆ

0.11 10

k

p















(11)

標準的方法 (Blsallfood data)

Trend

Noise

Seasonal

2

1

2

8

2 2,

12 ˆ

32.8907

ˆ

7.230 10

ˆ

1.100 10

1369.3

k

p

















AIC

標準的季節調整モデル

で問題が生じた例：

データの終盤でトレン

ドが波打っている

(12)

循環変動の抽出

p

_n



a p

₁

_n

_

₁







a p

_m

_n

_

_m



r

_n

循環変動成分

y

_n



t

_n



s

_n



p

_n



w

_n

₁₂

トレンド成分より短期的で定常変動成分

ARモデルで表現

(13)

季節調整モデル（循環変動成分つき）

n

m

n

m

n

p

n

k

v

p

a

p

a

p

v

s

v

t

3

1

2

1

1 )

(





















n

t

s

p

w

y









1

2

1

0 ,

,

1 0

0

1

0 n

m

n

p

a

p

x

F

G

H

p











 









 









 











 









 





 













(14)

状態空間表現

x

F x

G v

y

H x

w

n









 1



1

0

1

0

1

0 

0

1

0

1

0

1 ,

0

1

0

1

0

1

0

1

0

1

0

1

2 ]

,

|

,

|

,

[

2

1

10

1

1 













































_

H

G

a

F

p

s

t

x

_n

₁₄

(15)

標準的方法

循環変動成分を加えた方法

Trend

Cycle

Noise

Seasonal

(16)

長期予測の精度

循環変動成分を含む方法

Gersch & Kitagawa (1983)

標準的方法

(17)

長期予測

標準的方法

_{循環変動成分を含む方法}

季節成分の予測精度は同じでもトレンド

の予測精度が大きく異なる

(18)

18 data(Blsallfood)

season(Blsallfood, 2, 1)

$tau2

[1] 7.229607e-01 1.099998e-08

$sigma2

[1] 32.89071

$lkhood

[1] -669.6522

$aic

[1]

1369.304

season(Blsallfood, 2, 1, 2)

$tau2

[1] 7.229536e-03 1.099771e-08 9.999101e-01

$sigma2

[1] 27.02489

$lkhood

[1] -648.2725

$aic

[1]

1336.545

$arcoef

[1] 1.3579988 -0.5400455

(19)

data(Blsallfood)

(20)

season(Blsallfood, 2, 1, filter=c(1,132))

長期予測（n=133,156)

(21)

season(Blsallfood, 2, 1, 2, filter=c(1,132))

(22)

22 season(Whard, 2, 1, log=TRUE)

$tau2

[1] 2.593525e-02 1.099988e-08

$sigma2

[1] 0.0001501694

$lkhood

[1] 342.5324

$aic

[1]

-655.0648

season(Whard, 2, 1, 0, 6, year=1967, log=TRUE)

$tau2

[1] 1.269665e-01 1.100012e-08

$sigma2

[1] 5.774088e-05

$lkhood

[1] 358.2322

$aic

[1]

-674.4645

(23)

曜日効果の検出



₁









₇



0 td

d

n



 



 





 





1

6

1

6

7

1

7

6

7 ,

,

(

)

(

)

(

)



td

_n





₁

d

_n

_,

₁









₇

d

_n

_,

₇

曜日効果項

d

_{n j}

_,

_{第n月中の}

_j

_{番目の曜日の数}

y

_n



t

_n



s

_n



t d

_n



w

_n

日月火水木金土

J

4

5

4

4 F

4

4 M

4

5

5 A

5

4

4 M

4

5

4

4 d

_{n j}_,

(24)

季節調整モデル

_(一般形)

n

t

s

p

td

r

w

y





n

s

t

n

td

p

n

w

r

トレンド

季節成分

定常変動

曜日効果

外生変数

ノイズ

x

F x

G v

y

H x

w

n









1

状態空間表現

1

2

1

3

1

7

7 (

)

k

n

n p

n

m

n m

n

t

v

s

v

p

a p

v

td



d



d



 







 













成分モデル

Kitagawa & Gersch (1984)

(25)

曜日効果を含む季節調整モデルの状態空間表現

k=2, p=4,ar=0,trade=6の場合

1

2

1

6 *

*

1

6

7

2

1 1 0

1

0 0 0

1

1 0 1

1

0

0 0 0

,

0

1

0 0 0

1

0 0 0

0 0

0

1 0 0

1 0 1 0 0

,

n

j

t

S

F

G

x

S

H

d





_

_









































































_

_



_

_

_{ }

_















































































_

_











(26)

26 曜日調整なし

0.7

0.9

1.1

1.3

1.5

1

25

49

73

97

121

145

0.7 0.9 1.1 1.3 1.5 1 25 49 73 97 121 145

-0.1

0

0.1

1

25

49

73

97

121

145 -0.1

0

0.1

1

25

49

73

97

121

145 トレンド

ノイズ

季節成分

観測データ

(27)

曜日調整つき

0.8 0.9 1 1.1 1.2 1.3 1.4 1.5 1 13 25 37 49 61 73 85 97 109 121 133 145 -0.1 0 0.1 1 13 25 37 49 61 73 85 97 109 121 133 145 -0.1 0 0.1 1 13 25 37 49 61 73 85 97 109 121 133 145 0.8 0.9 1 1.1 1.2 1.3 1.4 1.5 1 13 25 37 49 61 73 85 97 109 121 133 145 -0.008 -0.004 0 0.004 0.008 1 2 3 4 5 6 7

曜日効果

0.8 0.9 1 1.1 1.2 1.3 1.4 1.5 1 13 25 37 49 61 73 85 97 109 121 133 145

--- T + S + TD

--- 観測データ

(28)

28 曜⽇効果パターン

-0.01 -0.006 -0.002 0.002 0.006 1 2 3 4 5 6 7 2.6 3 3.4 3.8 1 25 49 73 97 121 145 169 193 217 241 265 -0.2 0 0.2 1 25 49 73 97 121 145 169 193 217 241 265 2.6 3 3.4 3.8 1 25 49 73 97 121 145 169 193 217 241 265 東京⼤学北川源四郎数理⼿法VII （時系列解析）

(29)

曜日調整つき季節調整

2.6

3

3.4

3.8

1

25

49

73

97

121

145

169

193

217

241

265 --- T + S + TD

--- 観測データ

(30)

(31)

data(Whard) # Wholesale hardware data

(32)

32 曜日効果の時間変化の検出

-0.012

-0.006

0

0.006

1 23 45 67 89 111 133 155 177 199 221 243 265

Sunday

Monday

Tuesday

Wednesday

Thirsday

Friday

Saturday

-0.012 -0.009 -0.006 -0.003 0 0.003 0.006 0.009 1 2 3 4 5 6 7 -0.012 -0.009 -0.006 -0.003 0 0.003 0.006 0.009 1 2 3 4 5 6 7

jn

n

j

jn









_,

_1



（曜日効果を時変係数にする）

(33)

制約モデルの導入

土日と平日だけのモデル

1

7 ,

2

3

4

5

6 





















1

2

1

2

5

2 



5 



0 



 



（自由パラメータは個）

1

(34)

景気動向の推定

GDP系列の季節調整

Cycleの階差

GDP系列のCycle

経企庁（福田）

₃₄

(35)

低S/Nデータからの信号抽出

観測モデル

n

r

s

w

y





n

w

s

r

成分モデル

常微動

地震波

)

,

0 (

~

)

,

0 (

~

)

,

0 (

~

2

1

1 



N

w

N

v

N

u

v

s

b

s

u

r

a

r

n

j

n

j

n

j

n

m

j

n

j

n

j

n





















常微動

地震波

観測ノイズ

x

F x

G w

y

H x

n









 1



地震波常微動＋観測雑音

(36)

微小地震の抽出

n

r

s

w

y





n

r

s

w

Observed

常微動（背景雑音）

地震波

観測ノイズ

₃₆

(37)

信号分離のためのモデル

観測モデル

n

r

s

w

y





n

r

s

w

常微動

地震波

観測ノイズ

)

,

0 (

~

)

,

0 (

~

)

,

0 (

~

2

1

1 



N

w

N

v

N

u

v

s

b

s

u

r

a

r

n

j

n

l

j

n

j

n

m

j

n



















成分モデル

(38)

状態空間モデル

2

1

2

1

2

1

1 ,

0

0 ,

]

0

1

0

1 [

0

1

0

1 ,

0

1

0

1

1 ,







































































n

l

m

l

n

m

n

R

Q

H

G

b

a

F

s

r

x















x

F x

G w

y

H x

n









 1



₃₈

(39)

微小信号の抽出

観測時系列

地震波

(40)

微小地震波の抽出

観測データ

地震波

常微動

0 400 800 1200 1600 2000 2400 2800

15

0 -15

15

0 -15

15

0 -15

4

2

0 -2

-4

-6

時変分散

(in log

₁₀

)

_{10,000 times}

実用化には

自己組織型

状態空

間

モデルを

利用

₄₀

(41)

観測値

S-波

P-波

背景雑音

(42)

42 地下水位データ（地震の影響の抽出）

1981

1985

1984

1983

1982

1990

1989

1988

1987

1986

産総研観測井

気圧，地球潮汐，降雨

などの影響を受ける

地震の影響の検出が困難

(43)

欠測値と異常値の処理

状態空間モデル

)

,

(

)

,

0 (

)

1 (

~

)

(

Mixture

)

(

)

(

Cauchy

)

,

0 (

~

)

(

Gauss

2 2 2 2 2

















N

w

r

w

r

N

w

r









観測ノイズのモデル

Noise

n

Observatio

:

Noise

System

:

Signal

:

n

Observatio

:

n n n n

w

v

t

y

n

w

t

y

v

t









_1

0 0.2 0.4 0.6 0.8 1 41 81 121 161 201 241 281 321 361 401 Gauss Cauchy 0.2 0.4 0.6 0.8 Mixture 0 500 1000 1500 2000 1 21

異常値

欠測値

Model

AIC

Gauss

-8741

Cauchy -8655

Mixture -8936

(44)

観測データ

1981

₁₉₈₂

1983

1984

1985

1986

1987

1988

1989

₁₉₉₀

トレンド周りに強い変動

成分が存在

地震の影響の検出が困難

₄₄

(45)

地下水位データ（産総研・榛原）

単位m，観測間隔1時間

地下水位データとトレンド

(46)

46

(47)





0

1

0 ,

,

1

0

1 n

n

m

n

n m

t

a

x

F

G

a

H

p

_

 





 

 





 





 



 





 





 





 











気圧効果モデル

n

w

Hx

y

Gv

Fx

x









_1

n

y

 

t

P





n n n n

y

t

P



観測値

トレンド

気圧効果

観測ノイズ

成分モデル

1

0 n

n

m

n

i

n i

i

t

w

P

a p











状態空間モデル

モデル

次数

AIC



2

トレンド

1 -21311

0.181×10

-2

(48)

48 気圧効果項

原データと気圧補正値

(49)

気圧補正後の残差のヒストグラム

気圧補正後の残差の自己相関関数

(50)

50 

1 

,

0

1 ,

1

1 ,

0 0  









 













































n n m n n m n n

et

p

H

G

F

b

a

t

x

地球潮汐モデル

n

w

Hx

y

Gv

Fx

x









_1

n

t

P

E

y







n n n n n

E

P

t

y



観測値

トレンド

気圧効果

地球潮汐効果

観測ノイズ

成分モデル

i n l i i n i n m i i n n n k

et

b

E

p

a

P

w

t

   







0 0

状態空間モデル

(51)

地球潮汐補正後の地下水位，潮汐効果項，残差

地球潮汐補正後の残差の自己相関関数

モデル

AIC



2 トレンド

-21311

0.181×10

-2

気圧効果

-57836

0.120×10

-5

潮汐効果

-59580

0.887×10

-6

推定したモデル

(52)

地震の影響の検出（降雨モデル）

成分モデル

y

_n



t

_n



P

_n



E

_n



R

_n





_n

n

P

t

y

n

R

E



観測値

トレンド

気圧効果

潮汐効果

降雨効果

観測ノイズ

n

i

n

k

i

n

k

i

n

i

n

l

i

n

i

n

m

i

n

k

v

r

d

R

c

R

et

b

E

p

a

P

w

t

























1

0

0 ,

,

₅₂

(53)

モデル

AIC



2 _残差分布

_{自己相関関数}

トレンド

-21311

0.181×10

-2

負の偏り

3日程度

気圧効果

-57836

0.120×10

-5

正の偏り

12.5時間周期

潮汐効果

-59580

0.887×10

-6

1日程度まで正

降雨効果

-61810

0.130×10

-6

ほぼ対象

2時間周期

降雨補正後の残差の自己相関関数

推定したモデルと残差の特徴

(54)

AIC Values

₅₄

61809

6 61810

5 61800

4 61803

3 61734

2 61675

1 AIC



k

25, 1

m



l



l

0

1

2

3 m

…

23 -57819

-59368

-59498

-59488

24 -57815

-59374

-59536

-59526

25 -57830

-59393

-59580

-59569

26 -57836

-59386

-59575

-59563

27 -57832

-59379

-59566

-59554

(55)

降雨補正後の地下水位，残差，降雨効果項

(56)

56 地下水位の観測値

気圧効果

潮汐効果

降水効果

min AIC model

m=25, l=2, k=5

地震の影響の抽出

地下水位の補正値

_{M=4.8, D=48km}

(57)

移動平均 vs. 統計的モデリング

トレンド

（移動平均）

気圧・地球潮汐

気圧・地球

潮汐・降雨

機械的統計処理

統計的モデリング

原データ

(58)

58 1981

1982

1983

1984

1985

1986

1987

₁₉₈₈

1989

₁₉₉₀

M=7.0 D=375km M=4.8 D=48km M=5.7 D=66km M=7.7 D=622km M=6.0 D=113km M=6.2 D=150km M=5.0 D=57km M=7.9 D=742km M=6.8 D=128km M=6.0 D=126km M=6.7 D=226km M=5.7 D=122km M=6.5 D=96km

1981

₁₉₈₂

1983

1984

1985

1986

1987

1988

1989

₁₉₉₀

地震の影響の検出

(59)

地震の影響

地震

水位

降雨

水位

> 16cm > 4cm > 1cm

1 log

45 .

2 



D

M

0 log

45 .

2 



D

M

マグニチュード

水位変化

距離

(60)

得られた知⾒

 M > 2.6 log D + 0.2 の場合、影響が検出される

 変化量は M － 2.62 log D の関数で決まる

 地震がない場合およそ6cm/年の割合で水位が上昇

この地域のひずみの増加に対応か？

₆₀

(61)

OBS (海底地震計) Data

Bottom

OBS

Source

Sea Surface

－地下構造の探索－

(62)

62 100 1000 10000 100000 351 451 551 651 751 851 951

Time

Location

Amplitude

-10000 0 10000 20000 30000 40000 50000 1 501 1001 1501 2001 2501 3001 3501 4001 4501 5001 5501

Location

観測時系列（調整済，982ch）

(63)

時空間モデルによる分解

j

n

j

n

j

n

j

n

r

s

w

y

_,



_,



_,



_,

s

j

n

j

h

n

j

n

r

j

n

j

k

n

j

n

u

s

u

r

j

,

1 ,

,

1 ,

,











X

0

0 X

0 W

of

time

Arrival

:

)

W

(

W

of

time

Arrival

:

)

W

(

)

W

(

),

W

(

X

j

T

h

T

k









s

n

m

n

m

n

r

n

v

s

b

s

b

s

v

r

a

r

a

r

,

1

1 ,

1

1 













_

wave

Reflection

e

Direct wav

,

j

n

j

n

s

r

Time series model

Spatial model

Basic observation model

J-1

J

(64)

経路モデルと到着時刻

23 13 03 3 2 2

2

2 ,

/

v

d

D

d

h

v

d

_ij



_i _i _j



_i





3 1 3 2 23 2 2 1 2 2 13 2 1 1 1 2 03 2 0 1 0 12 02 1 2 2 12 2 1 1 1 2 02 2 0 1 0 12 02 1 2 2 12 2 1 1 1 2 02 2 0 1 0 01 1 1 2 01 2 0 1 0 01 1 1 2 01 2 0 1 0 2 2 0 1 0 2 2 0 1 0 2 2 0 1 0

2

2 1)

Wave(01232

)

2

3 (

2

3 1)

Wave(00012

)

2 (

2 Wave(0121)

)

3 (

3 Wave(0001)

)

(

Wave(01)

25 )

Wave(00000

9 Wave(000)

Wave(0)

d

v

d

h

v

d

h

v

d

h

v

d

D

v

d

h

v

d

h

v

d

D

v

d

h

v

d

h

v

d

D

v

d

h

v

d

D

v

d

h

v

D

h

v

D

h

v

D

h

v

                



















経路モデル

到着時刻

0

h

0

d

1

h

2

h

D

1

d

Water

OBS

Source

2

d

Parallel Layer Structure

Width

h

₀

,

h

₁

,

h

₂

,

h

₃



km

Velocity

v

₀

,

v

₁

,

v

₂

,

v

₃



km/sec

(65)

経路モデルと到着時刻（OBS4）

0 1 2 3 4 5 6 7 8 9 10 11 12 13 1 11 21 31 41 51 61 71 81 91 101

W0

W000

W00000

W0000000

W01

W011

W0121

W000121

W01221

W012321

W00012321

W0123321

W0111221

W0122221

W01111

W0001

W00011

W0001221

0 10 20 30 40 0 5 10 15 20 W0 W000 W0^5 W0^7 W0^9 W01 W011 W0121 W01221 W012321

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

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

北川 源四郎

時系列解析（１０）

2

季節調整とは

0.7

0.9

1.1

1.3

1.5

1

25

49

73

97

121

145

WHARD

月次データ

何らかの原因で特定の周期で繰り返す成分を

除去して本質的な現象を抽出する方法

季節調整モデル

y

n

: 観測値

t

n

: トレンド成分

s

n

: 季節成分

w

n

: 不規則成分（ノイズ）

n

n

n

n

t

s

w

y







非定常時系列をいくつかの成分に分解する方法

・

合成モデルによる成分分解

・本質的な変動を抽出する

観測モデル

4

トレンド成分モデル



k





n

k

n

k

n

n

t



0



t



0



t



v

1

1

2

1

2

3

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

北川源四郎

_n

_n

_n

_n