読者です 読者をやめる 読者になる 読者になる

東京に棲む日々

データ分析、統計、ITを勉強中。未だ世に出ず。

状態空間モデル - R dlm (2)

前回より続く。

 

線形成長モデルは以下に定義される。

 

Y[t] = μ[t] + v[t],   v[t]~N(0, V)

μ[t] = μ[t-1] + β[t-1] + w1[t],   w1[t]~N(0, W1)

β[t] = β[t-1] + w2[t],   w2[t]~N(0, W2)

 

 

動的線形モデルで表現する場合、各要素は以下になる。

 

t=0,        θ[0]~N(m[0], C[0])

t>0,        Y[t] = Fθ[t] + v[t],  v[t]~N(0, V)                             観測値

              θ[t] = Gθ[t-1] + w[t],  w[t]~N(0, W)                     状態

 

f:id:High_School_Student:20141115132931j:plain

 

線形成長モデルは、ローカルレベルモデルにトレンドβが加わったもの。

 

W2=0と場合を考えてみると良い。そうするとβはtに依存しない定数となり、直線的に減少なり増加するトレンドを持つ時系列データを表現できる。初期値β[0]やW2の大きさで、線形の傾向を調整できる。

(実際のモデリングでW2=0とした場合、モデルは可制御でなくなる。可制御であるとはモデルが到達できない領域がなく、可制御でないとはモデルが到達できない領域がある、と解釈する。)

部分部分(ローカル)で線形な傾向を持つデータへ使用されるモデルだと考えることができる。

 

V = 10、W1 = 6、W2 = 3、μ[0] = 10、β[0] = 1 としてデータを作成してみる。 

W2が小さく、β[0]を大きく設定すると、グローバルに直線傾向が強いデータが作成される。

 

V <- 10
W_mu <- 6
W_trend <- 3

mu.0 <- 10
trend.0 <- 1 # このトレンドの初期値が大きければ、直線的なトレンドの傾向が強くなる

tm <- c(0)

mu <- c(mu.0)
trend <- c(trend.0)

y <- c(NA)

# 40サンプル
for(i in 1:40){
     tm <- append(tm, i)

     trend_next <- trend[i] + rnorm(1, mean=0, sd=sqrt(W_trend))
     trend <- append(trend, trend_next)

     mu_next <- mu[i] + trend_next + rnorm(1, mean=0, sd=sqrt(W_mu))
     mu <- append(mu, mu_next)

     y_next <- mu_next + rnorm(1, mean=0, sd=sqrt(V))
     y <- append(y, y_next)
}

df2 <- data.frame(time=tm, trend=trend, mu=mu, y=y)


minY <- min(na.omit(df2$y))-10
maxY <-max(na.omit(df2$y))+10

par(mfrow = c(2, 1))
plot(df2$time, df2$y, ylim=c(minY,maxY),type="o", ylab="y, mu", xlab="time") # 黒 : データ
par(new=TRUE)
plot(df2$time, df2$mu, ylim=c(minY,maxY), type="o", col="blue", ylab="", xlab="") # 青 : mu
plot(df2$time, df2$trend,type="o", ylab="trend", xlab="time") # 黒 : トレンド
abline(h=0, lty=2)

f:id:High_School_Student:20141115133257j:plain

黒 : データy

: μ

黒: トレンドβ (グラフは、トレンドは別に描いている。)

 

では、モデルを関数で定義し、dlmMLE(データ, パラメータ初期値ベクトル, モデル関数)でパラメータを推定してみる。

 

fn_mod2 <- function(parm) {
     V_estim <-parm[1]
     W_mu_estim <- parm[2]
     W_trend_estim <- parm[3]

     return( dlmModPoly( order=2, dV=exp(V_estim), dW=c(exp(W_mu_estim),exp(W_trend_estim)) ) )
}

 

parm_start <- c(10*runif(1), 10*runif(1), 10*runif(1))

 

fn_mod2(parm_start)

$FF
     [,1] [,2]
[1,]    1    0

$V
         [,1]
[1,] 111.6958

$GG
     [,1] [,2]
[1,]    1    1
[2,]    0    1

$W
         [,1]     [,2]
[1,] 1445.989    0.000
[2,]    0.000 1852.898

$m0
[1] 0 0

$C0
      [,1]  [,2]
[1,] 1e+07 0e+00
[2,] 0e+00 1e+07

状態の初期値m[0]とC[0]を指定しなければ、それらの要素は0と1e+07にデフォルトで指定される。

 

res_fn_mod2 <- dlmMLE(y=df2$y[-1], parm=parm_start, build=fn_mod2)

 

dlmMLE()でのパラメータ推定結果は以下。

V_est <- exp(res_fn_mod2$par[1])

9.692269

W_mu_est <- exp(res_fn_mod2$par[2])

3.757845

W_trend_est <- exp(res_fn_mod2$par[3])

7.397736

 

V= 9.692、W1= 3.757、W2= 7.397となる。真の値V = 10、W1 = 6、W2 = 3と異なる。

関数に渡す初期値を大きくしたり小さくしたりで推定をシミュレーションしてみると、渡す初期値の大きさで結果が大きく変わる。

 

今回のデータ数は40。データ数を数百程度に増やしてシミュレーションしてみても、データによって結果は近かったり遠かったりとまちまち。

そもそも一系列のデータyから、複数の分散などを推定する訳である。一意に値が定まる訳でもないであろうし、推定自体簡単なことではないことだと想像がつく。推定に関しては今後も調査継続の必要性あり。

 

観測値1期先予測(f)と、状態のフィルタリング(m)をプロットしてみる。

f:id:High_School_Student:20141115134236j:plain

黒 (○): データy

(○): μ

黒(×) : 観測値1期先予測

(×): フィルタリングμ

 

黒(○): トレンドβ   (グラフは、トレンドは別に描いている。)

黒(×): フィルタリングβ

 

yやμのフィルタリングはある程度上手く行っているように見えるが、トレンドβの推定に関しては、ちょっと微妙な気がしないでもない。

真の値(V = 10、W1 = 6、W2 = 3)をモデルに指定しての推定も試ししてみたが同様の結果で、上手く行くとは限らない。

継続調査の必要あり。

 

分散を見てみる。

 

フィルタリング

dlmSvd2var(fit_mod2_est$U.C, fit_mod2_est$D.C)

これの最後のリスト。

41
         [,1]      [,2]
[1,] 7.495514  4.031254
[2,] 4.031254 13.754984

 

1期先観測値予測

dlmSvd2var(fit_mod2_est$U.R, fit_mod2_est$D.R)

これの最後のリスト。

40
         [,1]     [,2]
[1,] 33.07085 17.78624
[2,] 17.78624 21.15272

 

推定したフィルタや予測値の分散共分散行列は対角行列とならない。

モデルの定義では、状態のイノベーション項の分散共分散行列(mod2_est$W)は対角行列とするので、最初、これらの計算結果の出力に違和感があり戸惑ったが、線形成長モデルのμとβ間に相関があることは自明である。

 

k=10期先予測。

 

k <- 10 # k=10

mod2_pred <- dlmForecast(fit_mod2_est, nAhead=k)

str(mod2_pred, 1)

List of 4
$ a: num [1:10, 1:2] 94.2 95.3 96.4 97.5 98.6 ...
$ R:List of 10
$ f: num [1:10, 1] 94.2 95.3 96.4 97.5 98.6 ...
$ Q:List of 10

aは状態なので、μとβの2列の行列となる。

Rは状態の分散共分散行列(2×2)。

f、Qは長さ10のベクトルと解釈できる。

 

a_mu_pred <- unlist(mod2_pred$a)[,1]
a_trend_pred <- unlist(mod2_pred$a)[,2]
f_pred <- unlist(mod2_pred$f)


tm <- c( df2_1$time, length(df2_1$time):(length(df2_1$time)+length(f_pred)-1) )
tr <- c(df2_1$trend, rep(NA,length(f_pred)))
mu <- c(df2_1$mu, rep(NA,length(f_pred)))
y <- c(df2_1$y, rep(NA,length(f_pred)))
m_tr <- c(df2_1$m_trend, rep(NA,length(f_pred)))
m_mu <- c(df2_1$m_mu, rep(NA,length(f_pred)))
a_tr <- c(df2_1$a_trend, a_trend_pred)
a_mu <- c(df2_1$a_mu, a_mu_pred)
f <- c(df2_1$f, f_pred)

df2_2 <- data.frame( time=tm, trend=tr, mu=mu, y=y, m_trend=m_tr, m_mu=m_mu, a_trend=a_tr, a_mu=a_mu, f=f )

 

   time       trend          mu          y     m_trend       m_mu     a_trend         a_mu            f
1     0  1.00000000  10.0000000         NA  0.00000000   0.000000          NA           NA           NA
2     1 -0.39189751  11.3792097  15.271752  7.63587107  15.271745  0.00000000   0.00000000   0.00000000
3     2 -0.27076775  11.2985444   7.616363 -7.65529646   7.616393  7.63587107  22.90761608  22.90761608
4     3 -0.86313367   6.4755666   4.722802 -5.03339589   4.091164 -7.65529646  -0.03890386  -0.03890386
5     4 -1.62576596   5.9510299  11.275962  0.30797550   8.824934 -5.03339589  -0.94223215  -0.94223215
6     5 -0.35411010   9.2943873   8.234765 -0.06646885   8.434325  0.30797550   9.13290905   9.13290905
7     6  0.41779153   8.7413314   3.370188 -2.14387961   4.501004 -0.06646885   8.36785609   8.36785609
8     7 -0.98943300   5.4297462   5.423797 -0.86825547   4.728817 -2.14387961   2.35712451   2.35712451
9     8 -0.77556706   7.6265669   8.362854  1.00459543   7.342526 -0.86825547   3.86056160   3.86056160
10    9  0.11205921  10.8966641  13.466506  3.13395497  12.306263  1.00459543   8.34712184   8.34712184
11   10 -0.28963512   8.6545542   9.477838  0.65404066  10.829192  3.13395497  15.44021832  15.44021832
12   11 -1.13876371   9.0680897  16.081112  2.56641211  15.039004  0.65404066  11.48323282  11.48323282
13   12  0.15358498  12.3559981  15.920635  1.86567026  16.302491  2.56641211  17.60541614  17.60541614
14   13 -0.04898526  10.3677927  16.102531  1.00652330  16.570706  1.86567026  18.16816149  18.16816149
15   14  1.44431649  10.4584454   8.183906 -2.90039216  10.312904  1.00652330  17.57722945  17.57722945
16   15  0.37100105  14.2978069  15.276352  0.37037322  13.494011 -2.90039216   7.41251195   7.41251195
17   16 -1.31472734   9.0689054  15.552997  1.07270896  15.170273  0.37037322  13.86438459  13.86438459
18   17 -1.66637109   4.1789927   2.907222 -4.47396319   5.929774  1.07270896  16.24298190  16.24298190
19   18 -5.26607314  -1.1893775   2.527121 -4.02837868   2.284309 -4.47396319   1.45581074   1.45581074
20   19 -4.45189524  -8.2568128 -15.576649 -9.78168978 -12.441492 -4.02837868  -1.74407013  -1.74407013
21   20 -1.25306575  -9.2654350 -10.873195 -5.06094989 -13.445672 -9.78168978 -22.22318215 -22.22318215
22   21  0.11424844  -9.2279975 -15.507593 -3.81357973 -16.187323 -5.06094989 -18.50662166 -18.50662166
23   22  3.82678457 -11.4033086 -13.048461 -0.92188742 -14.624233 -3.81357973 -20.00090297 -20.00090297
24   23  4.10467212  -3.5449498  -3.869699  3.93462463  -6.516162 -0.92188742 -15.54612067 -15.54612067
25   24  3.95499814   0.1072197   4.519471  6.88810939   2.910026  3.93462463  -2.58153764  -2.58153764
26   25  5.71548099   9.2240111   7.336158  5.86411239   7.894166  6.88810939   9.79813530   9.79813530
27   26  6.87021215  20.3360060  14.913178  6.34446366  14.651420  5.86411239  13.75827829  13.75827829
28   27  8.72990299  32.6152238  29.464687  9.86684844  27.545231  6.34446366  20.99588336  20.99588336
29   28  7.33977724  41.4001828  42.137907 11.83243660  41.066797  9.86684844  37.41207992  37.41207992
30   29  6.89292503  48.5144668  52.967099 11.86066321  52.951717 11.83243660  52.89923396  52.89923396
31   30  7.21971539  56.3389234  57.208671  8.69809251  58.932053 11.86066321  64.81238035  64.81238035
32   31  3.81404530  61.1853035  64.224510  7.28160456  64.996398  8.69809251  67.63014524  67.63014524
33   32  5.05958442  60.8899078  60.959617  2.57400845  63.524931  7.28160456  72.27800214  72.27800214
34   33  3.45025926  72.4182725  71.442502  4.79652835  70.231383  2.57400845  66.09893913  66.09893913
35   34  2.57374664  75.5788671  79.736036  6.75475372  78.668939  4.79652835  75.02791110  75.02791110
36   35  1.90856435  80.3943913  80.145944  4.55960767  81.342146  6.75475372  85.42369226  85.42369226
37   36  2.23365459  87.6591308  88.498695  5.63973962  87.910098  4.55960767  85.90175389  85.90175389
38   37  3.05362375  90.0239314  91.289013  4.69940675  91.801429  5.63973962  93.54983746  93.54983746
39   38  1.06568152  90.6666095  89.647751  1.84903939  91.201004  4.69940675  96.50083577  96.50083577
40   39  0.46012447  89.9896574  91.953965  1.39315362  92.202392  1.84903939  93.05004300  93.05004300
41   40  2.34790032  93.2382139  92.920977  1.11258411  93.073868  1.39315362  93.59554516  93.59554516
42   41          NA          NA         NA          NA         NA  1.11258411  94.18645220  94.18645220
43   42          NA          NA         NA          NA         NA  1.11258411  95.29903632  95.29903632
44   43          NA          NA         NA          NA         NA  1.11258411  96.41162043  96.41162043
45   44          NA          NA         NA          NA         NA  1.11258411  97.52420454  97.52420454
46   45          NA          NA         NA          NA         NA  1.11258411  98.63678866  98.63678866
47   46          NA          NA         NA          NA         NA  1.11258411  99.74937277  99.74937277
48   47          NA          NA         NA          NA         NA  1.11258411 100.86195688 100.86195688
49   48          NA          NA         NA          NA         NA  1.11258411 101.97454100 101.97454100
50   49          NA          NA         NA          NA         NA  1.11258411 103.08712511 103.08712511
51   50          NA          NA         NA          NA         NA  1.11258411 104.19970922 104.19970922

 

信頼区間は描かず、予測平均のみを示す。

par(mfrow = c(2, 1))
minY <- min(na.omit(df2_2$y))-20
maxY <-max(na.omit(df2_2$y))+30
plot(df2_2$time, df2_2$y, ylim=c(minY,maxY),type="o", ylab="y, mu, f, m_mu", xlab="time") # 黒 : データ
par(new=TRUE)
plot(df2_2$time, df2_2$mu, ylim=c(minY,maxY), type="o", col="red", ylab="", xlab="") # 赤 : mu
par(new=TRUE)
plot(df2_2$time, df2_2$f, ylim=c(minY,maxY), type="o", pch=4, lty=3, ylab="", xlab="") # 黒× : f
par(new=TRUE)
plot(df2_2$time, df2_2$a_mu, ylim=c(minY,maxY), type="o", pch=4, lty=3, col="red", ylab="", xlab="") # 赤× : a_mu
minTY <- min(na.omit(df2_2$trend))-5
maxTY <- max(na.omit(df2_2$trend))+5
plot(df2_2$time, df2_2$trend, ylim=c(minTY,maxTY), type="o", ylab="trend", xlab="time") # 黒 : トレンド
par(new=TRUE)
plot(df2_2$time, df2_2$a_trend, ylim=c(minTY,maxTY), type="o", pch=4, lty=3, ylab="trend", xlab="time") # ? : a_trend
abline(h=0, lty=2)

 

f:id:High_School_Student:20141115135558j:plain

黒 (○): データy

(○): μ

黒(×) : 観測値1期、k期先予測 f

(×): 1期、k期先予測 μ

 

黒(○): トレンドβ   (グラフは、トレンドは別に描いている。)

黒(×): 1期、k期先予測トレンドβ

 

トレンドβは最後に1期先予測を行った後、定数となる。μはその時点からベータの累積和となり、直線的に増加または減少する。μと観測値予測fは同じになり、図ではfはμに上書きされている。