添付資料 - TOKYO TECH OCW

データ解析
第十五回「時系列解析」
鈴木 大慈
理学部情報科学科
西八号館 W707 号室
[email protected]
1 / 37
今日の講義内容
時系列解析
TOPIX30 データで時系列解析
2 / 37
構成
1
時系列解析の基本事項
2
TOPIX30 データで時系列解析
3 / 37
2000
3000
x
4000
5000
6000
時系列
1992
1993
1994
1995
1996
1997
1998
Time
4 / 37
時系列
確率過程は時刻 t でインデックス付けされた確率変数の列 {Xt }t ．
時間的な系列がランダムに生成されていると考える．
i.i.d. 系列との違いは異なる時刻で依存性を持つ点．
t は連続値でも離散値でも良い．今回は離散値のみを考える．
5 / 37
ガウス過程
Definition (ガウス過程)
確率過程 {Xt }t がガウス過程である.
⇔ ∀t1 , . . . , tN に対して Xt1 , . . . , XtN の同時分布が多変量ガウス分
布に従う．
6 / 37
定常性
Definition (強定常)
確率過程 {Xt }t が強定常である．
⇔ ∀t1 , . . . , tN , ∀h に対して Xt1 , . . . , XtN の同時分布が
Xt1 +h , . . . , XtN +h の同時分布と等しい．
「時間をずらしても分布が変わらない」と言っている．
→ 過去のデータから未来の挙動が予測できる．
7 / 37
−0.3
−0.2
−0.1
z
0.0
0.1
0.2
0.3
定常性
0
50
100
150
Index
定常 (Xt = 0.7Xt + 0.1ϵt , AR(1) 過程)
8 / 37
z
0.0
−20
−10
0.5
0
x
10
1.0
20
30
1.5
定常性
0
50
100
150
200
0
Index
50
100
150
Index
定常でない
(なぜか？)
8 / 37
平均と自己共分散
平均:
µt := E[Xt ]
(∀t)
定常ならば µt は t によらない．µ = µt と書く．
ある時刻に平均から大きく離れてもまた戻ってくる．
自己共分散:
γ(t, s) := Cov(Xt , Xs ) = E[(Xt − µt )(Xs − µs )]
定常ならば γ(t, s) は時間差 t − s のみに依存する．
γ(h) = γ(t, t + h) と書く．
9 / 37
弱定常
強定常はやや強い条件である．平均と共分散 (相関) だけでも時間
に依存しなければ十分では？
Definition (弱定常)
確率過程 {Xt }t が弱定常である．
⇔ 平均 µt が時間によらず，自己共分散 γ(t, s) が時間差 t − s に
のみ依存．
10 / 37
ガウス過程の定常性
ガウス過程の場合
弱定常
⇔
強定常
ガウス過程は時系列解析で最も基本的な確率過程である．
11 / 37
非定常から定常へ: 季節成分の除去
非定常時系列は季節成分を除去することで定常になったりする．
m <- decompose(co2)
# co2 は timeseries データ
plot(m) # 原データ、傾向成分、季節成分、不規則成分
340
0.0
−0.5
observed
seasonal
random
0.5
−3 −2 −1 0
1
2
3 320
330
trend
350
360
320 330 340 350 360
Decomposition of additive time series
1960
1970
1980
1990
Time
12 / 37
非定常から定常へ：季節成分の除去 (stl)
関数 stl は decompose より高度．局所多項式近似を使って平滑化．
stllc <- stl(co2, "periodic")
plot(stllc)
340
−0.5
0.0
remainder
0.5
320
330
trend
350
360
−3 −2 −1
0
seasonal
1
2
3
data
320 330 340 350 360
1960
1970
1980
1990
time
13 / 37
非定常から定常へ：平滑化
核関数で平滑化
xsmooth <- kernapply(x,kernel("daniell", 10))
# Daniell 核関数で平滑化
2000
3000
x
4000
5000
6000
1992
1993
1994
1995
1996
1997
1998
Time
より詳しくはスクリプト参照
14 / 37
自己相関
自己相関:
γ(t, s)
.
ρ(t, s) := √
γ(t, t)γ(s, s)
弱定常ならば時間差 t − s のみに依存．
acf(stllc$time.series[,"remainder"])
先の CO2 データから季節成分・トレンドを取り除き自己相関を計算．
stllc$time.series[, "remainder"]
−0.2
0.0
0.2
ACF
0.4
0.6
0.8
1.0
Series
0.0
0.5
1.0
Lag
1.5
2.0
15 / 37
自己共分散
自己共分散もオプションを指定すれば計算できる．
acf(stllc$time.series[,"remainder"],
type = "covariance")
stllc$time.series[, "remainder"]
−0.02
0.00
ACF (cov)
0.02
0.04
0.06
Series
0.0
0.5
1.0
1.5
2.0
Lag
16 / 37
AR 過程
ちょっと前までの実現値に依存して今の値が決まる確率過程
例：株価が 3 日連続で上がり続けたら，次は下りやすい．
17 / 37
AR 過程
ちょっと前までの実現値に依存して今の値が決まる確率過程
例：株価が 3 日連続で上がり続けたら，次は下りやすい．
AR 過程：
Xt =
p
∑
ϕi Xt−i + ϵt .
i=1
ただし，ϵt ∼ N(0, σ 2 ) (i.i.d.).
何時刻前まで遡るかの値 p を AR 次数と呼ぶ．
17 / 37
R で AR 過程の推定
ar(x, aic = TRUE, order.max = NULL,
method = c("yule-walker", "burg", "ols", "mle",
"yw"),
デフォルトで AR 次数は AIC で決定．
order.max で AR 次数を探索する最大の次数を指定．
method で推定方法を指定．デフォルトは”yule-walker”．
18 / 37
Yule-Walker 法
時系列は定常であると仮定する．自己回帰モデルの定義より．
Xt = Xt−1 ϕ1 + · · · + Xt−p ϕp + ϵt
⇒Xt−h Xt = Xt−h Xt−1 ϕ1 + · · · + Xt−h Xt−p ϕp + Xt−h ϵt
⇒γ(h) = γ(h − 1)ϕ1 + · · · + γ(h − p)ϕp . (両辺期待値を取った)
これより次の Yule-Walker 方程式を得る:

 
 
γ(1)
γ(0)
γ(1)
· · · γ(p − 1)
ϕ1
γ(2)  γ(1)


γ(0)
· · · γ(p − 2) ϕ2 

 

  .. 
 ..  = 
..
..
..
..

 .  

.
.
.
.
.
γ(p)
γ(p − 1) γ(p − 2) · · ·
γ(0)
ϕp
| {z } |
{z
} | {z }
γ
Γ
ϕ
γ と Γ をサンプルからの推定量に置き換え，方程式を解く:
b =b
b.
ϕ
Γ−1 γ
19 / 37
0
50
ar.co2$aic
100
150
AR 次数の決定: AIC
0
5
10
15
20
25
Index
> ar.co2 <- ar(stllc$time.series[,"remainder"])
> ar.co2$aic
で AIC の値を確認できる．
20 / 37
AR(1) 過程の非定常性：単位根
Xt = ϕ1 Xt−1 + ϵt
AR(1) 係数 ϕ1 が |ϕ1 ≥ 1| であると分散が発散してゆき，AR 過程は
非定常になる．
21 / 37
R で単位根の検定
Dickey-Fuller の単位根検定
> adf.test(UKgas)
Augmented Dickey-Fuller Test
data: UKgas
Dickey-Fuller = -1.6079, Lag order = 4, p-value = 0.7393
alternative hypothesis: stationary
22 / 37
ベクトル自己回帰 (VAR)
Xt ∈ Rd がベクトルであるとする．
例：複数の株価
ベクトル AR 過程 (VAR):
Xt = A1 Xt−1 + · · · + Ap Xt−p + ϵt .
なお，Ai ∈ Rd×d は行列．
23 / 37
R でベクトル自己回帰
VAR の AR 次数も推定しつつ VAR を推定：
library(vars)
varsel <- VARselect(tsx,lag.max=5) #モデル選択
var.topix <- VAR(tsx,p=varsel$selection[1]) #AIC でモ
デル選択した次数を採用
24 / 37
構成
1
時系列解析の基本事項
2
TOPIX30 データで時系列解析
25 / 37
TOPIX CORE 30 銘柄時系列
TOPIX CORE 30 指標に含まれる 30 銘柄．
過去 250 日分 (2013/6/28–2014/7/4) の始値，安値，高値，終
値 (日足) が格納されている．
26 / 37
380
4400
420
......
460
4800
5000
..................
4600
0
50
100
Time
150
200
250
......
3800
...............
3400
5500
500 600 650 700 750 800 850 2000
3000
6500
2200
2400
...............
6000
..................
2600
280
300
320
340
4200
...............
3800
3600
.....................
3400
.....................
3200
1000 1100 1200 1300
3000
.........&.................................
70003400
例
ts(topix30[, 1:10])
0
50
100
Time
150
200
250
27 / 37
収益率の計算
前日と比べて上がったか下がったか．
対数収益率:
Rt = log(Xt /Xt−1 ).
こちらを分析対象とする．
28 / 37
自己相関
Series ts(logrt[, 1])
0.4
ACF
0.2
0.4
0.0
0.2
−0.2
0.0
ACF
0.6
0.6
0.8
0.8
1.0
1.0
Series ts(logrt[, 20])
0
5
10
15
Lag
東京三菱 UFJ
20
0
5
10
15
20
Lag
日本たばこ産業
29 / 37
AR モデルへの当てはめ
東京三菱 UFJ を当てはめ．
ufjts <- ts(logrt[,20]) #時系列データ構造に変換
ar.ufj <- ar(ufjts) #AR
30 / 37
15
10
5
0
ar.ufj$aic
20
25
30
AIC
0
5
10
15
20
0:(length(ar.ufj$aic) − 1)
UFJ データの AIC
31 / 37
推定された AR モデル
> ar.ufj$ar #AR 係数
[1] 0.10819647 0.05554658 -0.10994533 -0.11853056
32 / 37
残差の正規性のチェック
ちゃんとデータが AR 過程に従っているか．残差の正規性をチェッ
ク．線形回帰の場合と同様の考え方．
0.00
−0.02
−0.04
Sample Quantiles
0.02
0.04
Normal Q−Q Plot
−3
−2
−1
0
1
2
3
Theoretical Quantiles
33 / 37
残差の正規性検定
> shapiro.test(ar.ufj$resid[5:230])
Shapiro-Wilk normality test
data: ar.ufj$resid[5:230]
W = 0.9907, p-value = 0.1603
正規性の検定，棄却はされない．
34 / 37
おまけ：対数収益率の正負を判別
次の８つの銘柄を用いて，次に三菱ＵＦＪの株価が上がるか下が
るかを判別．(Rt > 0 か Rt < 0 かを判別)
三井物産
三菱商事
三菱ＵＦＪフィナンシャル・グループ
三井住友フィナンシャルグループ
みずほフィナンシャルグループ
野村ホールディングス
東京海上ホールディングス
三菱地所
説明変数としては，これらの銘柄の過去 3 日間の対数収益率を用
いた．
推定手法は L1 正則化付きロジスティック回帰．
35 / 37
0.0
0.2
0.4
tedata$y
0.6
0.8
1.0
結果
0
10
20
30
40
Index
150 サンプルで推定，残りの 96 サンプルで判別誤差を計算．
判別誤差： 0.6304348
常に上がると答えた場合：0.5
常に下がると答えた場合：0.5
ランダムな解答より良い正答率．
36 / 37
講義情報ページ
http://www.is.titech.ac.jp/~s-taiji/lecture/dataanalysis/dataanalysis.html
スクリプトやデータはすべてこのページ載せてあります．
37 / 37