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

東京に棲む日々

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

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

今回は季節要素モデル。

 

例えば、データが四半期(周期=4)だった場合を考える。

 

Y[t-1] = α1 + v[t-1]

Y[t] = α2 + v[t]

Y[t+1] = α3 + v[t+1]

Y[t+2] = α4 + v[t+2]

 

α1+α2+α3+α4 = w[t],  w[t] ~ N(0, w_seas)

 

α1,2,3,4が各四半期成分となり、和が平均0で分散w_seasの正規分布に従うとする。

 

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

 

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:20141115141140j:plain

ちなみに、

  f:id:High_School_Student:20141115141203j:plain

 

dlmModSeas()によって定義することができる。(dlmModTrig()を使うともっと柔軟に周期成分を捉えることができるようであるが、文献をさっと見ただけではよく理解できなかったので、必要があれば後々調べて行く。)

 

例えば、V=2、w_seas=4とした場合。 

dlmModSeas(frequency=4, dV=2, dW=c(4,0,0))

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

$V
     [,1]
[1,]    2

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

$W
     [,1] [,2] [,3]
[1,]    4    0    0
[2,]    0    0    0
[3,]    0    0    0

$m0
[1] 0 0 0

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

 

今回はこの季節要素モデル単体ではなく、ローカルレベルモデルと組み合わせてみる。ここら辺が状態空間モデルの柔軟なところである。

 

例えば、ローカルレベルモデルでV=3、W=6とし、季節要素モデル(四半期)でV=2、v_seas=4とした場合、以下のようにdlmModPoly()とdlmModSeas()を組み合わせる。

 

dlmModPoly( order=1, dV=3, dW=6 ) + dlmModSeas(frequency=4, dV=2, dW=c(4,0,0))

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

$V
     [,1]
[1,]    5

$GG
     [,1] [,2] [,3] [,4]
[1,]    1    0    0    0
[2,]    0   -1   -1   -1
[3,]    0    1    0    0
[4,]    0    0    1    0

$W
     [,1] [,2] [,3] [,4]
[1,]    6    0    0    0
[2,]    0    4    0    0
[3,]    0    0    0    0
[4,]    0    0    0    0

$m0
[1] 0 0 0 0

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

 

この場合、モデルは以下に定義される。

 

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

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

  α1+α2+α1+α2 = w2[t],         w2[t] ~ N(0, w_seas)

 

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

 

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:20141115141852j:plain

 ちなみに、

  f:id:High_School_Student:20141115141923j:plain

 

これまで見てきた各モデルの要素のベクトルや行列を、ブロックとして左や上から詰めたり、行列に対角に挿入したりの形式で表現される。

Vは、ローカルレベルと季節要素モデルそれぞれの観測誤差の和(3+2=5)になっている。

 

では、各パラメータを指定し、サンプルデータ(40行)を出力してみる。

 

alpha1 <- 10
alpha2 <- 20
alpha3 <- -5
alpha4 <- -alpha1 -alpha2 -alpha3 # -25
seas <- c(alpha1, alpha2, alpha3, alpha4)

V <- 3+2
w_loc <- 6
w_seas <- 4


tm <- c(0)
mu <- c(100) # μの初期値は100
alpha <- c(NA)
y <- c(NA)

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

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

     alpha_next <- seas[i%%4+1] + rnorm(1, mean=0, sd=sqrt(w_seas))
     alpha <- append(alpha, alpha_next)

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

df3 <- data.frame(time=tm, mu=mu, alpha=alpha, y=y)

   time        mu      alpha         y
1     0 100.00000         NA        NA
2     1 101.03818  17.955605 117.69493
3     2 101.25233  -4.091898  98.77220
4     3 100.06346 -23.062593  72.51658
5     4  96.04669   8.901576 102.28330
6     5  95.93837  20.853973 114.93580
7     6  95.13444  -4.282241  88.50010
8     7  97.25849 -24.034810  71.62968
9     8  91.51833   8.773162 103.92921
10    9  89.47757  23.573238 111.84888
11   10  88.85872  -3.415835  85.99318
12   11  85.27012 -25.549430  60.39846
13   12  83.69630  10.338072  95.27553
14   13  81.79563  19.841001  96.37555
15   14  84.31110  -7.703048  76.44482
16   15  82.14169 -24.087354  62.47521
17   16  83.97272   8.371363  94.25735
18   17  86.51272  22.397388 109.34330
19   18  84.96344  -6.344667  78.87386
20   19  88.88974 -29.737400  56.65816
21   20  86.93938   8.035606  93.87915
22   21  85.26500  17.822354 106.99937
23   22  86.98058  -5.784289  79.34397
24   23  90.86066 -25.089460  64.02854
25   24  88.00473   7.444889 100.68429
26   25  92.70302  16.964471 110.13245
27   26  98.23367  -4.378634  91.80236
28   27 103.50861 -25.343409  75.63039
29   28 101.71510   7.622101 105.03880
30   29  99.40863  20.895265 118.87026
31   30  98.33315  -5.755254  97.06423
32   31  98.25762 -24.370844  72.00850
33   32  96.20491  10.091563 103.16654
34   33  92.37165  23.834555 117.60553
35   34  90.52986  -2.675492  84.79213
36   35  87.80711 -23.419353  62.24392
37   36  87.24032   7.758830  91.15328
38   37  84.09756  20.816570 105.28377
39   38  85.87142  -7.620894  76.79782
40   39  83.18068 -26.800049  57.26111
41   40  79.01976  11.497422  87.86101

 

プロット。
minY <- min(na.omit(df3$y))-10
maxY <-max(na.omit(df3$y))+10
par(mfrow = c(2, 1))
plot(df3$time, df3$y, ylim=c(minY,maxY),type="o", ylab="y, mu", xlab="time") # 黒 : データ
par(new=TRUE)
plot(df3$time, df3$mu, ylim=c(minY,maxY), type="o", col="blue", ylab="", xlab="") # 青 : mu
abline(h=100, lty=2)
plot(df3$time, df3$alpha, type="o", ylab="seasonal", xlab="time") # 黒 : alpha
abline(h=0, lty=2)

f:id:High_School_Student:20141115142348j:plain

黒 : データy

: ローカルレベルモデル部分μ

黒: 季節要素部分α (グラフは、季節要素は別に描いている。)

 

では、関数定義と推定。

fn_mod3 <- function(parm) {
     V <- parm[1]
     w_loc <- parm[2]
     w_seas <- parm[3]

     return( dlmModPoly( order=1, dV=exp(V), dW=exp(w_loc) ) + dlmModSeas(frequency=4, dV=0, dW=c(exp(w_seas),0,0)) )
}

 

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

res_fn_mod3 <- dlmMLE(y=df3$y[-1], parm=parm_start, build=fn_mod3)


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

3.613708

w_loc_est <- exp(res_fn_mod3$par[2])

11.18024

w_seas_est <- exp(res_fn_mod3$par[3])

 

0.03253725

 

V= 3.613、w_loc= 11.180、w_seas= 0.032となる。これまでの推定結果で見てきた通り、真の値V = 5、w_loc=6、w_seas=4と異なる結果。

 

カルマンフィルタの実行。

mod3_est <- fn_mod3( res_fn_mod3$par )

fit_mod3_est <- dlmFilter(df3$y[-1], mod3_est)

str(fit_mod3_est, 1)


fit_mod3_est$m # 41×4の行列。必要なのは最初の2列。

           [,1]       [,2]       [,3]        [,4]
 [1,]   0.00000   0.000000   0.000000   0.0000000
 [2,]  29.42375  88.271165 -29.423722 -29.4237217
 [3,]  78.80987  19.962313  38.885103 -29.4237080
 [4,]  86.52002 -14.003451  12.252181  31.1749761
 [5,]  97.81681   4.466474 -25.300197   0.9554599
 [6,]  96.28570  18.841810   5.404486 -24.9552884
 [7,]  91.08020  -1.571212  20.276052   6.3781730
 [8,]  94.54692 -23.586630  -2.546678  19.6680823
 [9,]  96.33134   7.245955 -24.060232  -2.8045851
[10,]  93.62298  18.768130   7.954189 -23.7808371
[11,]  90.42758  -3.736562  19.320389   8.2521456
[12,]  86.20605 -24.885758  -3.006170  19.7127197
[13,]  86.81212   8.330082 -24.986893  -3.0549997
[14,]  79.70472  18.238856   9.520411 -24.5146983
[15,]  79.69382  -3.246505  18.240248   9.5210856
[16,]  84.88081 -23.592474  -3.907571  17.9172445
[17,]  84.73424   9.556773 -23.574268  -3.8994040
[18,]  89.55919  18.674913   8.953554 -23.8090888
[19,]  84.56592  -4.521445  19.188505   9.1835815
[20,]  81.61569 -24.265807  -4.218537  19.3252328
[21,]  83.87268   9.476095 -24.493025  -4.3140614
[22,]  86.65627  19.688730   9.193996 -24.5991631
[23,]  84.42293  -4.547578  19.884117   9.2753542
[24,]  87.53267 -24.244067  -4.819311  19.7703504
[25,]  90.37753   9.628814 -24.488764  -4.9157097
[26,]  90.36216  19.773962   9.630145 -24.4882853
[27,]  95.09226  -4.426344  19.408591   9.4886199
[28,]  98.81982 -24.085032  -4.714118  19.2968316
[29,]  96.37653   9.249946 -23.898850  -4.6455587
[30,]  98.76652  19.528847   9.066798 -23.9612860
[31,] 100.96262  -4.429610  19.374948   9.0114099
[32,]  97.21961 -24.305715  -4.167304  19.4693449
[33,]  94.93014   8.790651 -24.147025  -4.1128300
[34,]  97.33823  19.684315   8.622898 -24.2009000
[35,]  90.98728  -4.651370  20.096420   8.7605263
[36,]  87.57253 -24.498573  -4.429692  20.1702737
[37,]  83.67656   8.424338 -24.248019  -4.3496794
[38,]  84.75778  20.263011   8.354492 -24.2690601
[39,]  82.03062  -4.567471  20.429461   8.4060369
[40,]  81.65193 -24.298429  -4.544341  20.4365929
[41,]  79.99277   8.273214 -24.197904  -4.5145095

 

df3_1 <- data.frame( df3, m_loc=fit_mod3_est$m[,1], m_seas=fit_mod3_est$m[,2],
a_loc=c(NA,fit_mod3_est$a[,1]), a_seas=c(NA,fit_mod3_est$a[,2]), f=c(NA,fit_mod3_est$f) )

   time        mu      alpha         y     m_loc     m_seas     a_loc      a_seas            f
1     0 100.00000         NA        NA   0.00000   0.000000        NA          NA           NA
2     1 101.03818  17.955605 117.69493  29.42375  88.271165   0.00000   0.0000000 0.000000e+00
3     2 101.25233  -4.091898  98.77220  78.80987  19.962313  29.42375 -29.4237219 3.275159e-05
4     3 100.06346 -23.062593  72.51658  86.52002 -14.003451  78.80987 -29.4237081 4.938616e+01
5     4  96.04669   8.901576 102.28330  97.81681   4.466474  86.52002 -29.4237059 5.709632e+01
6     5  95.93837  20.853973 114.93580  96.28570  18.841810  97.81681  19.8782635 1.176951e+02
7     6  95.13444  -4.282241  88.50010  91.08020  -1.571212  96.28570   0.7089925 9.699469e+01
8     7  97.25849 -24.034810  71.62968  94.54692 -23.586630  91.08020 -25.0830134 6.599718e+01
9     8  91.51833   8.773162 103.92921  96.33134   7.245955  94.54692   6.4652266 1.010121e+02
10    9  89.47757  23.573238 111.84888  93.62298  18.768130  96.33134  19.6188627 1.159502e+02
11   10  88.85872  -3.415835  85.99318  90.42758  -3.736562  93.62298  -2.9414819 9.068150e+01
12   11  85.27012 -25.549430  60.39846  86.20605 -24.885758  90.42758 -23.8359724 6.659161e+01
13   12  83.69630  10.338072  95.27553  86.81212   8.330082  86.20605   8.1792081 9.438526e+01
14   13  81.79563  19.841001  96.37555  79.70472  18.238856  86.81212  19.7118112 1.065239e+02
15   14  84.31110  -7.703048  76.44482  79.69382  -3.246505  79.70472  -3.2445687 7.646015e+01
16   15  82.14169 -24.087354  62.47521  84.88081 -23.592474  79.69382 -24.5148283 5.517899e+01
17   16  83.97272   8.371363  94.25735  84.73424   9.556773  84.88081   9.5828005 9.446361e+01
18   17  86.51272  22.397388 109.34330  89.55919  18.674913  84.73424  17.9168992 1.026511e+02
19   18  84.96344  -6.344667  78.87386  84.56592  -4.521445  89.55919  -3.8193778 8.573981e+01
20   19  88.88974 -29.737400  56.65816  81.61569 -24.265807  84.56592 -23.8506414 6.071528e+01
21   20  86.93938   8.035606  93.87915  83.87268   9.476095  81.61569   9.1591109 9.077480e+01
22   21  85.26500  17.822354 106.99937  86.65627  19.688730  83.87268  19.3309907 1.032037e+02
23   22  86.98058  -5.784289  79.34397  84.42293  -4.547578  86.65627  -4.2835635 8.237271e+01
24   23  90.86066 -25.089460  64.02854  87.53267 -24.244067  84.42293 -24.6118931 5.981104e+01
25   24  88.00473   7.444889 100.68429  90.37753   9.628814  87.53267   9.2930274 9.682569e+01
26   25  92.70302  16.964471 110.13245  90.36216  19.773962  90.37753  19.7756597 1.101532e+02
27   26  98.23367  -4.378634  91.80236  95.09226  -4.426344  90.36216  -4.9158219 8.544634e+01
28   27 103.50861 -25.343409  75.63039  98.81982 -24.085032  95.09226 -24.4708669 7.062139e+01
29   28 101.71510   7.622101 105.03880  96.37653   9.249946  98.81982   9.5023191 1.083221e+02
30   29  99.40863  20.895265 118.87026  98.76652  19.528847  96.37653  19.2944622 1.156710e+02
31   30  98.33315  -5.755254  97.06423 100.96262  -4.429610  98.76652  -4.6343585 9.413216e+01
32   31  98.25762 -24.370844  72.00850  97.21961 -24.305715 100.96262 -23.9567479 7.700587e+01
33   32  96.20491  10.091563 103.16654  94.93014   8.790651  97.21961   9.0036746 1.062233e+02
34   33  92.37165  23.834555 117.60553  97.33823  19.684315  94.93014  19.4692037 1.143993e+02
35   34  90.52986  -2.675492  84.79213  90.98728  -4.651370  97.33823  -4.1063123 9.323192e+01
36   35  87.80711 -23.419353  62.24392  87.57253 -24.498573  90.98728 -24.2055755 6.678170e+01
37   36  87.24032   7.758830  91.15328  83.67656   8.424338  87.57253   8.7579920 9.633052e+01
38   37  84.09756  20.816570 105.28377  84.75778  20.263011  83.67656  20.1733599 1.038499e+02
39   38  85.87142  -7.620894  76.79782  82.03062  -4.567471  84.75778  -4.3484429 8.040934e+01
40   39  83.18068 -26.800049  57.26111  81.65193 -24.298429  82.03062 -24.2680267 5.776259e+01
41   40  79.01976  11.497422  87.86101  79.99277   8.273214  81.65193   8.4061765 9.005810e+01

 

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

minY <- min(na.omit(df3$y))-20
maxY <-max(na.omit(df3$y))+20
par(mfrow = c(2, 1))
plot(df3$time, df3$y, ylim=c(minY,maxY),type="o", ylab="y, mu", xlab="time") # 黒 : データ
par(new=TRUE)
plot(df3$time, df3$mu, ylim=c(minY,maxY), type="o", col="blue", ylab="", xlab="") # 青 : mu
abline(h=100, lty=2)
par(new=TRUE)
plot(df3_1$time, df3_1$f, ylim=c(minY,maxY), type="o", pch=4, lty=3, ylab="", xlab="") # ? : f
par(new=TRUE)
plot(df3_1$time, df3_1$m_loc, ylim=c(minY,maxY), type="o", pch=4, lty=3, col="blue", ylab="", xlab="") # ? : m_mu

minTY <- min(na.omit(df3_1$alpha))-10
maxTY <- max(na.omit(df3_1$alpha))+10
plot(df3$time, df3$alpha, ylim=c(minTY,maxTY), type="o", ylab="seasonal", xlab="time") # 黒 : alpha
abline(h=0, lty=2)
par(new=TRUE)
plot(df3_1$time, df3_1$m_seas, ylim=c(minTY,maxTY), type="o", pch=4, lty=3, ylab="seasonal", xlab="time") # ? : m_trend
abline(h=0, lty=2)

 

f:id:High_School_Student:20141115143450j:plain

黒 (○): データy

(○): ローカルレベルモデル部分μ

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

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

 

黒(○): 季節要素α   (グラフは、季節要素は別に描いている。)

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

 

前回の線形成長モデルのトレンド成分より、今回の季節要素の方が予測しやすいようである。5順くらいしたところからフィルタリング値が適応していく様子が分かる。