« 米国版電子ブロック | Home | RによるPIC24HJ12GP202専用ICSPソフト 暫定公開 »
Median polishで株価変動のトレンド除去
By kmgs | 12 月 29, 2008
Tukeyにより1977年に発表されたmedian polishは、トレンドの除去手法の一つとして知られており、ググってみると"中央値分散分析"などと訳されることもあるようです(マイクロアレイ屋さんにとっては、RMA法のコアとなるメソッドとして有名かもしれません)。アルゴリズムとしては、行列データに対して(i) 行ごとに中央値を差し引く、(ii) 列ごとに中央値を差し引く、という処理を収束するまで繰り返すものです。これにより、行列の各要素は(行ごとのトレンド)+(列ごとのトレンド)+(要素固有の誤差)+(オフセット値)という足し算に分解されます(加法分解)。
たとえば、株価の入った行列(行:日付、列:銘柄)に適用した場合、median polishにより(平均株価の推移ベクトル)+(銘柄の平均株価ベクトル)+(銘柄固有のパターン行列)という意味合いを持つ成分に分解されます。Rには組み込み関数としてmedpolish()がありますので、こいつを使って株価のトレンド除去を行ってみます。赤が正、緑が負の変動を示し、銘柄は意図的に除去してあります。
いつものようにRソースです。これを"kabu_trend.r"として保存してください。
なお、"core30.txt"については、2008年12月25日のエントリをご参照ください。
# kabu_trend.r
parse_html <- function(u, noVol=F) {
Sys.sleep(1)
a <- read.delim(u, as.is=T)
# 株価データらしき部分を抽出して整形、yahoo側の仕様変更があると使えなくなるかも?
sel <- a[grep("<t[dh]><small>(<b>)?[0-9]+(</b>)?", a[,1]),1]
sel <- gsub("<[^>]+>", "", sel)
sel <- gsub("[年月]", "-", sel)
sel <- gsub("[日,]", "", sel)
# 行列に整形、適宜データ型を揃える
if (noVol) {
# 日付/始値/高値/安値/終値
mat <- matrix(sel, length(sel)/5, 5, by=T)
nums <- matrix(as.numeric(mat[,-1]), nrow(mat), ncol(mat)-1)
dat <- data.frame(as.Date(mat[,1]), nums)
colnames(dat) <- c("date", "open", "high", "low", "close")
} else {
mat <- matrix(sel, length(sel)/7, 7, by=T)
nums <- matrix(as.numeric(mat[,-1]), nrow(mat), ncol(mat)-1)
dat <- data.frame(as.Date(mat[,1]), nums)
colnames(dat) <- c("date", "open", "high", "low", "close", "volume", "adj_close")
}
# 日付順をひっくり返してリターン
dat[nrow(dat):1,]
}
# 例:
# x <- kabu("7203.t")
kabu <- function(code, noVol=F) {
# HTMLの取得
u <- url(paste("http://table.yahoo.co.jp/t?s=", code, "&g=d", sep=""), enc="EUC-JP")
parse_html(u, noVol)
}
kabu_old <- function(code, days=50, noVol=F) {
# HTMLの取得
d <- strsplit(as.character(Sys.Date()), "-")[[1]]
u <- url(paste("http://table.yahoo.co.jp/t?c=2001&a=1&b=1&f=", d[1], "&d=", d[2], "&e=", d[3], "&g=d&s=", code, "&y=", days, sep=""), enc="EUC-JP")
parse_html(u, noVol)
}
core30 <- function() {
a <- read.delim("core30.txt", as.is=T)
l <- list()
str <- paste(a[,2], a[,3])
for (i in 1:nrow(a)) {
print(str[i])
l[[str[i]]] <- rbind(kabu_old(a[i,2], 100), kabu_old(a[i,2], 50), kabu(a[i,2]))
}
l
}
index <- function() {
l <- list()
code <- 998407
l[["NIKKEI"]] <- rbind(kabu_old(code, 100, T), kabu_old(code, 50, T), kabu(code, T))
code <- 998405
l[["TOPIX"]] <- rbind(kabu_old(code, 100, T), kabu_old(code, 50, T), kabu(code, T))
code <- "usdjpy=x"
l[["USDJPY"]] <- rbind(kabu_old(code, 100, T), kabu_old(code, 50, T), kabu(code, T))
code <- "eurjpy=x"
l[["EURJPY"]] <- rbind(kabu_old(code, 100, T), kabu_old(code, 50, T), kabu(code, T))
l
}
plot_list <- function(l) {
# 2009/01/11 終値から調整後終値に
mat <- log2(sapply(l, function(X) X[,7]))
mat <- sweep(mat, 2, colMeans(mat))
rownames(mat) <- as.character(l[[1]][,1])
mm <- max(abs(mat))
heatmap(mat, scale="none", mar=c(15,4), Rowv=NA, zlim=c(-mm, mm), col=c(rgb(0,100:0/100,0), rgb(1:100/100,0,0)))
}
medpol <- function(l, ind) {
p <- par(ask=T)
on.exit(par(p))
# 株価
# 2009/01/11 終値から調整後終値に
mat <- log2(sapply(l, function(X) X[,7]))
rownames(mat) <- as.character(l[[1]][,1])
mres <- medpolish(mat)
# 株価指標等
mat2 <- log2(sapply(ind, function(X) X[,5]))
rownames(mat) <- as.character(ind[[1]][,1])
print(colnames(mat2))
# 行トレンドと日経平均およびTOPIXとの相関
df <- scale(data.frame(core30_trend=mres[[2]], nikkei=mat2[,"NIKKEI"], topix=mat2[,"TOPIX"]))
matplot(df, type="l")
legend(100, 1, colnames(df), col=1:3, lty=1:3)
# 列トレンドと銘柄ごとの平均株価(幾何平均)の相関
plot(colMeans(mat), mres[[3]], xlab="Geometric means", ylab="Column effect")
med <- mres[[4]]
mm <- max(abs(med))
heatmap(med, scale="none", mar=c(15,4), Rowv=NA, zlim=c(-mm, mm), col=c(rgb(0,100:0/100,0), rgb(1:100/100,0,0)))
}以下のコマンドで実行します。株価取得は2分くらい掛かります。
> source("kabu_trend.r")
> l <- core30() # 株価の取得 (過去150取引日)
> ind <- index() # 株価指数等の取得 (過去150取引日)
> medpol(l, ind) # 作図1枚目の図は、行トレンドと日経平均およびTOPIXをそれぞれz変換したプロットです。100日目前後からのブレが目立ちますが、Core30の特性かもしれませんので、まあいいでしょう。
2枚目の図は、列トレンドと銘柄ごとの平均株価(幾何平均)の相関です。これはほぼ直線的な関係があります。
3枚目の図が誤差マトリックスのヒートマップであり、冒頭に示したものです(再掲)。
真ん中のクラスターは相場と相関する動きを示す銘柄(製造業中心)と思われます。対して、左側のクラスターは相場よりも落ち込みが激しいもの(金融系)、右側のクラスターは暴落相場のなかで耐えているもの(すなわち相対的な価値が上がっているもの、インフラや製薬等)になると考えられます。
同じ業種でも、一部銘柄では異なる動きが見えます。
また、今回はmedian polishに先立ち対数変換を行っているので、意味合い的には加法分解ではなく、乗法分解となります。
[改訂履歴]
2009/01/11: 使用する株価を、終値から調整後終値に変更 (株式分割の影響回避)
Topics: R |