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

R dlm をつかってみるだけ

  • やったのは以下。コード参照。
    • フィルタリング(ランダムウォーク・プラス・ノイズモデル)
    • 平滑化(ランダムウォーク・プラス・ノイズモデル)
    • フィルタリング(線型成長モデル)
    • 平滑化(線型成長モデル)
  • そもそもフィルタリングと平滑化とは何だったか(関連: 状態空間モデルのまとめ絵 - クッキーの日記)。
    • フィルタリング: 直接観測できない「状態」を、その時点の観測値までから推定。
    • 平滑化: 直接観測できない「状態」を、すべての観測値を手に入れた後から推定。
  • 以下の日足(直近250日)を拾ってきてつかった。R 内蔵のデータでなくわざわざこれをつかう意味はない。
  • このデータに対してこれらのモデルが(そもそもDLMが)適当ということではない。実際にはちゃんと考えてモデルを決める。
フィルタリング(ランダムウォーク・プラス・ノイズモデル)

   f:id:cookie-box:20160109164730p:plain:w600

平滑化(ランダムウォーク・プラス・ノイズモデル)

   f:id:cookie-box:20160109165142p:plain:w600

フィルタリング(線型成長モデル)

   f:id:cookie-box:20160109165633p:plain:w600

平滑化(線型成長モデル)

   f:id:cookie-box:20160109165754p:plain:w600

グラフの点線はフィルタ後/平滑化後の状態分布の95%点。
線型成長モデルについては、状態の平均には不確実性はなく(分散0)、状態の変化(成長)に不確実性がある(こちらに分散がある)と考えたモデルにしているので、グラフには分散をプロットしていない。

ランダムウォーク・プラス・ノイズモデルを描画するコードは以下(フィルタリングは dlmFilter の1行だけ)。

library(dlm)

# データ読み込み
nsa <- read.table("I101.csv", header=F, skip=2, sep=",", stringsAsFactors=FALSE)
nsa <- nsa[rev(1:nrow(nsa)),] # 日付逆順のため
colnames(nsa) <- c("date", "open", "high", "low", "close")
nsa <- nsa[as.Date(nsa$date) >= as.Date("2015-03-01"),] # 絞り込み
nsa <- nsa[as.Date(nsa$date) <= as.Date("2015-07-31"),] # 絞り込み

# グラフ枠描画
plot(1:nrow(nsa), nsa$close, type="n", main="Nikkei Stock Average", 
     xlab="", ylab="closing price", xaxt="n")
indices <- seq(1, nrow(nsa), by=5)
text(indices, par("usr")[3], labels=nsa$date[indices], srt=90, pos=2, offset=0, xpd=TRUE)
abline(v=indices, h=seq(18000, 22000, by=500), col="dimgray", lty=3)

# 凡例を溜めておく用
note <- data.frame(name="REAL", col="black", stringsAsFactors=FALSE)

# ランダムウォーク・プラス・ノイズモデルのテンプレート
rw.template <- dlm(FF=1, V=0, GG=1, W=0,
                   m0=18797.94, C0=500*500) # 事前平均は2月終値

# 描画するグラフの切り替えフラグ
flag.process <- "smoothing" # "filtering", "smoothing"

# 描画用の便利関数(ランダムウォーク・プラス・ノイズモデル)
add.plot.rw <- function(V, W, col) {
  rw <- rw.template 
  V(rw) <- V
  W(rw) <- W
  nsaFilt <- dlmFilter(nsa$close, rw) # フィルタリング
  if (flag.process == "filtering") { # フィルタリングならフィルタ化分布を描画
    lines(0:nrow(nsa), nsaFilt$m, lwd=2, col=col)
    sigma <- sqrt(unlist(dlmSvd2var(nsaFilt$U.C, nsaFilt$D.C))) # 各点のフィルタ後分散も描画
    lines(0:nrow(nsa), nsaFilt$m + qnorm(0.025) * sigma, lwd=2, lty=3, col=col)
    lines(0:nrow(nsa), nsaFilt$m + qnorm(0.975) * sigma, lwd=2, lty=3, col=col)
  } else if (flag.process == "smoothing") { # 平滑化なら、さらに平滑化分布を描画
    nsaSmooth <- dlmSmooth(nsaFilt) # 平滑化
    lines(0:nrow(nsa), nsaSmooth$s, lwd=2, col=col)
    sigma <- sqrt(unlist(dlmSvd2var(nsaSmooth$U.S, nsaSmooth$D.S))) # 各点の平滑化後分散も描画
    lines(0:nrow(nsa), nsaSmooth$s + qnorm(0.025) * sigma, lwd=2, lty=3, col=col)
    lines(0:nrow(nsa), nsaSmooth$s + qnorm(0.975) * sigma, lwd=2, lty=3, col=col)
  }
  note <- rbind(note, c(paste0("V=", V, ", W=", W), col))
  return(note)
}

# 分散を色々変えて描く
note <- add.plot.rw(V=40000, W=0,    col="blue")  
note <- add.plot.rw(V=40000, W=100,  col="forestgreen") 
note <- add.plot.rw(V=40000, W=2500, col="orange")   

# 元データと凡例
lines(1:nrow(nsa), nsa$close, lwd=2, col="black")
legend("bottomright", legend=note$name, lwd=2, col=note$col, box.lwd=NA)

線型成長モデルの場合は該当箇所を以下に置き換える。

# 線型成長モデルのテンプレート
lg.template <- dlm(FF=matrix(c(1, 0),       nrow=1), V=0, 
                   GG=matrix(c(1, 0, 1, 1), nrow=2), W=0*diag(2),
                   m0=c(18797.94, 0), C0=500*500*diag(2)) # 事前平均は2月終値

# 描画用の便利関数(線型成長モデル)
add.plot.lg <- function(V, W1, W2, col) {
  lg <- lg.template 
  V(lg) <- V
  W(lg)[1, 1] <- W1
  W(lg)[2, 2] <- W2
  nsaFilt <- dlmFilter(nsa$close, lg) # フィルタリング
  if (flag.process == "filtering") { # フィルタリングならフィルタ化分布を描画
    lines(0:nrow(nsa), nsaFilt$m[,1], lwd=2, col=col)
  } else if (flag.process == "smoothing") { # 平滑化なら、さらに平滑化分布を描画
    nsaSmooth <- dlmSmooth(nsaFilt) # 平滑化
    lines(0:nrow(nsa), nsaSmooth$s[,1], lwd=2, col=col)
  }
  note <- rbind(note, c(paste0("V=", V, ", W1=", W1, ", W2=", W2), col))
  return(note)
}

# 分散を色々変えて描く
note <- add.plot.lg(V=40000, W1=0, W2=0,   col="blue")        
note <- add.plot.lg(V=40000, W1=0, W2=1,   col="forestgreen") 
note <- add.plot.lg(V=40000, W1=0, W2=100, col="orange")