京橋のバイオインフォマティシャンの日常

南国のビーチパラソルの下で、Rプログラムを打ってる日常を求めて、、Daily Life of Bioinformatician in Kyobashi of Osaka

【R言語と株価プロット】S&P500インデックス (SPY: SPDR S&P500 ETF) と仮想通貨(ビットコイン/BTC、イーサリアム/ETH) の季節性変動(Seasonality)を考えてみた件

(この記事は、2021年3月に執筆した記事を再校正しています。)

はじめに

2021年3月のアメリカ株式市場は、初旬からテック株、小型株がかなりの暴落でした。まさか、ここまで下がるとは、、、*1

こうなると、常識的な相場のSeasonality(季節性)をちゃんと理解して、相場の乗り降りのタイミングを考えないといけないな、身を持って強く思えてきました(泣泣泣)。

この記事では、S&P500インデックスに連動する S&P500 ETF(ティッカーシンボル: SPY)に加えて、仮想通貨のビットコイン(BTC)やイーサリアム(ETH)の価格情報を取得して、Seasonality(季節変動)をグラフ化してみることにしました。 これらの実行には、無料の解析環境である、R言語を用いています。

また、このアイディアをもとに、 月ごとの上昇・下落を色分けしてプロット作成できる、seasonalityPlotというRパッケージも作成しました。

seasonalityPlot package (version 0.99.4) | seasonalityPlot

CRAN - Package seasonalityPlot

また、記事の後半部分に、Rでの実行コードと補足事項を記載しています。

seasonalityPlotパッケージを用いた季節変動性のプロット

関連パッケージのインストール

株価などの季節変動を可視化する、seasonalityPlotパッケージを用います。

インストール方法は選べる2タイプです。

#CRANパッケージのインストール(旧バージョン)
install.packages("seasonalityPlot", repos="http://cran.r-project.org")
library(seasonalityPlot)

##GitHubパッケージのインストール(新バージョン、開発版)
install.packages("devtools", repos="http://cran.r-project.org")
devtools::install_github("kumeS/seasonalityPlot")
library(seasonalityPlot)

seasonalityPlotパッケージを使うのに、その他のパッケージは、特段必要ないです。

SPYチャートの季節変動性プロット

まずは、SPYの季節変動(2000-2021年)をプロットします。

seasonPlot(Symbols="SPY",
           StartYear = 2000,
           EndYear = 2021)

続いて、SPYの季節変動(2010-2021年)をプロットします。

seasonPlot(Symbols="SPY",
           StartYear = 2010,
           EndYear = 2021)

ビットコイン・チャートの季節変動性プロット

ビットコイン(BTC-USD)の季節変動(2014-2021年)をプロットします。

seasonPlot(Symbols="BTC-USD",
           StartYear = 2014,
           EndYear = 2021)

イーサリアム・チャートの季節変動性プロット

イーサリアム(ETH-USD)の季節変動(2015-2021年)をプロットします。

seasonPlot(Symbols="ETH-USD",
           StartYear = 2015,
           EndYear = 2021)

R言語/quantmodを用いた、SPYチャートのプロット(過去バージョン)

1994年から2021年までのSPYチャート

まずは、1994年から2021年までのSPYの株価チャートを示します。圧巻の右上がり、、2021年3月の下がりもほぼノイズの域です、、

2021年1-3月のSPYチャート

こっちが、2021年1-3月のSPYの株価チャートです。

月中旬でピークで、月終わりに下がってくるパターンではあるが、あれ、上がってるけど、、、大型セクターは結構安定なことをよく物語っています。

SPYチャートのいろんなバリエーション

次に、3パターンの図を作成してみました*2

2000-2020年平均のパターンだと、2月中旬がピークで、3月3週目あたりが底となるみたいです。

2010-2020年平均のパターンも上記と形が似てるけど、3月初めからの下りが激しく、3月最終週あたりから上がりだすようです。

最後に、大統領選翌年平均のパターンだけど、最悪想定ではあるが、4月中旬まで底が続くパターンがあるようです。

まとめ

これくらいのパターンは想定して、春相場を見ていかないとダメと言うことなのね。

ネットを見てると、3月は ・ 債券のショート ・ 株 => 債券へのリバランス の傾向が強いらしく、2月までに株で儲かり過ぎたのだろう。

次は、8月あたりで売りつけれるかが山場になりそうだけど。。。

補足

strftime関数などにおけるformatについて

%A: 曜日の名称(Sunday, Monday ... )
%a: 曜日の省略名(Sun, Mon ... )
%B: 月の名称(January, February ... )
%b: 月の省略名(Jan, Feb ... )
%j: 年中の通算日(001-366)
%U: 週を表す数。最初の日曜日が第1週の始まり(00-53)
%u: 月曜日を1とした、曜日の数値表現 (1..7)
%V: ISO 8601形式の暦週 (01..53)
%v: VMS形式の日付 (%e-%^b-%4Y)
%W: 週を表す数。最初の月曜日が第1週の始まり(00-53)
%w: 曜日を表す数。日曜日が0(0-6)

plot図XY軸の遊びをなくする引数

このオプションを用いることで、プロット領域の両端をxlim、およびylimで指定した範囲のプロットを生成します。

xaxs="i"
yaxs="i"

text()のちょっと高級版: 背景色とかが設定できる。

plotrix::boxed.labels(x,y=NULL,labels,
  bg=ifelse(match(par("bg"),"transparent",0),"white",par("bg")),
  border=TRUE,xpad=1.2,ypad=1.2,srt=0,cex=1,adj=0.5,xlog=FALSE,ylog=FALSE)

月の英語略語

month.abb
# [1] "Jan" "Feb" "Mar" "Apr" "May" "Jun" "Jul" "Aug" "Sep"
#[10] "Oct" "Nov" "Dec"

Rでの実行コード

#パッケージ・ロード
library(quantmod)
library(magrittr)
library(dygraphs)
library(plotrix)
library(htmltools)

#SPYの株価取得
NN01 <- 1994; NN02 <- 2021
Date <- c(paste0(NN01, "-01-01"), paste0(NN02, "-12-31"))
SPY <- quantmod::getSymbols("SPY", src = "yahoo", verbose = T, auto.assign=FALSE, from = Date[1], to=Date[2])

head(SPY)
#           SPY.Open SPY.High  SPY.Low SPY.Close SPY.Volume SPY.Adjusted
#1994-01-03 46.59375 46.65625 46.40625  46.46875     960900     28.06298
#1994-01-04 46.53125 46.65625 46.46875  46.65625     164300     28.17622
#1994-01-05 46.71875 46.78125 46.53125  46.75000     710900     28.23285
#1994-01-06 46.81250 46.84375 46.68750  46.75000     201000     28.23285
#1994-01-07 46.84375 47.06250 46.71875  47.03125     775500     28.40267
#1994-01-10 47.09375 47.59375 46.96875  47.59375     593700     28.74237

str(SPY)
#An ‘xts’ object on 1994-01-03/2021-03-25 containing:
#  Data: num [1:6856, 1:6] 46.6 46.5 46.7 46.8 46.8 ...
# - attr(*, "dimnames")=List of 2
#  ..$ : NULL
#  ..$ : chr [1:6] "SPY.Open" "SPY.High" "SPY.Low" "SPY.Close" ...
#  Indexed by objects of class: [Date] TZ: UTC
#  xts Attributes:  
#List of 2
# $ src    : chr "yahoo"
# $ updated: POSIXct[1:1], format: "2021-03-26 00:43:34"

#全体プロット
SPY[,4] %>%
  dygraphs::dygraph(main = "SPDR S&P 500 (SPY) Stock Price 26-years") %>% 
  dygraphs::dySeries("SPY.Close", label = "SPY") %>%
  dygraphs::dyRangeSelector(height = 40)  %>%
  htmltools::save_html(file="SPY_00.html")

#xtsオブジェクトをデータフレームに変換する
#データは終値だけにする
SPY00 <- data.frame(Date=index(SPY), SPY[,c(4)], row.names=1:nrow(SPY))
head(SPY00)
#        Date SPY.Close
#1 1994-01-03  46.46875
#2 1994-01-04  46.65625
#3 1994-01-05  46.75000
#4 1994-01-06  46.75000
#5 1994-01-07  47.03125
#6 1994-01-10  47.59375

#年、週目などの情報取得
SPY00$Year <- substring(SPY00$Date, 1, 4)
SPY00$Week <- strftime(SPY00$Date, format = "%V")
SPY00$Year.Week <- paste0(SPY00$Year, SPY00$Week)
SPY00$Day <- sub("-", "", substring(SPY00$Date, 6, 10))
head(SPY00)
#        Date SPY.Close Year Week Year.Week  Day
#1 1994-01-03  46.46875 1994   01    199401 0103
#2 1994-01-04  46.65625 1994   01    199401 0104
#3 1994-01-05  46.75000 1994   01    199401 0105
#4 1994-01-06  46.75000 1994   01    199401 0106
#5 1994-01-07  47.03125 1994   01    199401 0107
#6 1994-01-10  47.59375 1994   02    199402 0110
tail(SPY00)
#           Date SPY.Close Year Week Year.Week  Day
#6851 2021-03-18    391.48 2021   11    202111 0318
#6852 2021-03-19    389.48 2021   11    202111 0319
#6853 2021-03-22    392.59 2021   12    202112 0322
#6854 2021-03-23    389.50 2021   12    202112 0323
#6855 2021-03-24    387.52 2021   12    202112 0324
#6856 2021-03-25    387.90 2021   12    202112 0325

#いったん、1年372日分作成
Days <- c()
for(n in 1:12){
Days <- c(Days, paste0(formatC(n, width=2, flag = "0"), formatC(seq(1, 31, 1), width=2, flag = "0")))
}
SPY01 <- data.frame(Date=Days)
Years <- paste0("Y", seq(NN01, NN02, 1))
SPY01[,Years] <- NA
head(SPY01)
#  Date Y1994 Y1995 Y1996 Y1997 Y1998 Y1999 Y2000 Y2001 Y2002 Y2003 Y2004 Y2005 Y2006 Y2007 Y2008 Y2009 Y2010
#1 0101    NA    NA    NA    NA    NA    NA    NA    NA    NA    NA    NA    NA    NA    NA    NA    NA    NA
#2 0102    NA    NA    NA    NA    NA    NA    NA    NA    NA    NA    NA    NA    NA    NA    NA    NA    NA
#3 0103    NA    NA    NA    NA    NA    NA    NA    NA    NA    NA    NA    NA    NA    NA    NA    NA    NA
#  Y2011 Y2012 Y2013 Y2014 Y2015 Y2016 Y2017 Y2018 Y2019 Y2020 Y2021
#1    NA    NA    NA    NA    NA    NA    NA    NA    NA    NA    NA
#2    NA    NA    NA    NA    NA    NA    NA    NA    NA    NA    NA
#3    NA    NA    NA    NA    NA    NA    NA    NA    NA    NA    NA

#日付、年ごとにデータ整理
for(m in 1:length(Years)){
a <- SPY00[SPY00$Year == sub("Y", "", Years[m]),]
SPY01[c(SPY01$Date %in% a$Day),m+1] <- a$SPY.Close
}

#ここは、特に実行しなくて良い
#欠損数のカウント
#SPY01$MissingNum <- sapply(data.frame(t(SPY01[,-1])), function(x) sum(is.na(x)))
#株価データが無い日を消す。
#table(SPY01$MissingNum == length(Years))
#SPY01 <- SPY01[SPY01$MissingNum != length(Years),]
#head(SPY01)

#それぞれの年初日を「0%」とする
SPY02 <- SPY01
for(m in 1:length(Years)){
b <- SPY02[,m+1]
SPY02[,m+1] <- SPY02[,m+1]/b[!is.na(b)][1]*100 - 100
}
head(SPY02)
#  Date     Y1994     Y1995      Y1996     Y1997      Y1998    Y1999
#1 0101        NA        NA         NA        NA         NA       NA
#2 0102        NA        NA  0.0000000 0.0000000  0.0000000       NA
#3 0103 0.0000000 0.0000000  0.2765904 1.4352047         NA       NA
#4 0104 0.4034970 0.4778157 -0.6789037        NA         NA 0.000000
#5 0105 0.6052455 0.4778157 -0.8800603        NA  0.2242152 1.143002
#6 0106 0.6052455 0.5802048         NA 0.5487547 -1.3773222 3.581407

#株価データが無い日を削除して、株価データを繋げる
SPY03 <- data.frame(Date=1:nrow(SPY01))
SPY03[,Years] <- NA
for(m in 1:length(Years)){
d <- SPY02[,m+1][!is.na(SPY02[,m+1])]
SPY03[,m+1] <- c(d, rep(NA, nrow(SPY01)-length(d)))
}

#株価データが無い(欠損数2以上)行を削除する
SPY03$MissingNum <- sapply(data.frame(t(SPY03[,-1])), function(x) sum(is.na(x)))
SPY03 <- SPY03[c(SPY03$MissingNum < 2),]

SPY04 <- SPY03
colnames(SPY04)
# [1] "Date"       "Y1994"      "Y1995"      "Y1996"      "Y1997"      "Y1998"      "Y1999"      "Y2000"     
# [9] "Y2001"      "Y2002"      "Y2003"      "Y2004"      "Y2005"      "Y2006"      "Y2007"      "Y2008"     
#[17] "Y2009"      "Y2010"      "Y2011"      "Y2012"      "Y2013"      "Y2014"      "Y2015"      "Y2016"     
#[25] "Y2017"      "Y2018"      "Y2019"      "Y2020"      "Y2021"      "MissingNum"
dim(SPY04)
#[1] 248  30

#図 1. 2021年 SPY のSeasonality
Main <- "図 1. 2021年 SPY のSeasonality"
Y <- c(29)
SPY04$Mean <- NA
SPY04$Mean <- apply(data.frame(SPY04[,c(Y)]), 1, function(x) mean(x[!is.na(x)]))

COL <- "blue1"
St <- SPY04$Mean; St00 <- St[!is.na(St)]
par(family="HiraKakuProN-W3", lwd=1, xpd=F, cex=1, mgp=c(0.5, 1, 0), mai=c(0.5, 0.75, 0.5, 0.5))
plot(St, type="n", axes = F, xlab="© 2021 京橋のバイオインフォマティシャンの日常 by skume", 
     ylab="", cex.lab=0.75, xlim=c(0,nrow(SPY04)), 
     ylim=c(min(St00) - (max(St00)- min(St00))*0.1, max(St00) + (max(St00)- min(St00))*0.15), 
     xaxs="i", yaxs="i", main=Main)
axis(side=2, labels=paste0(seq(-30, 30, by=1), "%"), at=seq(-30, 30, by=1), las=2)
abline(v=seq(0, nrow(SPY04), length.out = 13), col="grey", lty=3, lwd=0.5)
abline(h=seq(-30, 30, by=1), col="black", lty=1, lwd=0.3)
lines(St, col=COL, lwd=1.2)
plotrix::boxed.labels(seq(0, nrow(SPY04), length.out = 13)[1:12]+10,  min(St00) - (max(St00)- min(St00))*0.05, 
                      month.abb, cex=0.7, bg = "grey20", xpad = 1.2, ypad = 1.2)
legend("topleft", legend="SPDR S&P500 ETF", col=COL, lwd=1, cex=0.7)
quartz.save(file = paste0("./SPY_01.png"), type = "png", dpi = 300); dev.off()

#図 2. 2000-2020年 SPY のSeasonality
Main <- "図 2. 2000-2020年 SPY のSeasonality"
Y <- c(8:28)
#head(SPY04[,c(Y)])
SPY04$Mean <- NA
SPY04$Mean <- apply(data.frame(SPY04[,c(Y)]), 1, function(x) mean(x[!is.na(x)]))

COL <- "blue1"
St <- SPY04$Mean; St00 <- St[!is.na(St)]
par(family="HiraKakuProN-W3", lwd=1, xpd=F, cex=1, mgp=c(0.5, 1, 0), mai=c(0.5, 0.75, 0.5, 0.5))
plot(St, type="n", axes = F, xlab="© 2021 京橋のバイオインフォマティシャンの日常 by skume", 
     ylab="", cex.lab=0.75, xlim=c(0,nrow(SPY04)), 
     ylim=c(min(St00) - (max(St00)- min(St00))*0.1, max(St00) + (max(St00)- min(St00))*0.15), 
     xaxs="i", yaxs="i", main=Main)
axis(side=2, labels=paste0(seq(-30, 30, by=1), "%"), at=seq(-30, 30, by=1), las=2)
abline(v=seq(0, nrow(SPY04), length.out = 13), col="grey", lty=3, lwd=0.5)
abline(h=seq(-30, 30, by=1), col="black", lty=1, lwd=0.3)
lines(St, col=COL, lwd=1.2)
plotrix::boxed.labels(seq(0, nrow(SPY04), length.out = 13)[1:12]+10,  min(St00) - (max(St00)- min(St00))*0.05, 
                      month.abb, cex=0.7, bg = "grey20", xpad = 1.2, ypad = 1.2)
legend("topleft", legend="SPDR S&P500 ETF", col=COL, lwd=1, cex=0.7)
quartz.save(file = paste0("./SPY_02.png"), type = "png", dpi = 300); dev.off()

#図 3. 2010-2020年 SPY のSeasonality
Main <- "図 3. 2010-2020年 SPY のSeasonality"
Y <- c(18:28)
#head(SPY04[,c(Y)])
SPY04$Mean <- NA
SPY04$Mean <- apply(data.frame(SPY04[,c(Y)]), 1, function(x) mean(x[!is.na(x)]))

COL <- "red1"
St <- SPY04$Mean; St00 <- St[!is.na(St)]
par(family="HiraKakuProN-W3", lwd=1, xpd=F, cex=1, mgp=c(0.5, 1, 0), mai=c(0.5, 0.75, 0.5, 0.5))
plot(St, type="n", axes = F, xlab="© 2021 京橋のバイオインフォマティシャンの日常 by skume", 
     ylab="", cex.lab=0.75, xlim=c(0,nrow(SPY04)), 
     ylim=c(min(St00) - (max(St00)- min(St00))*0.1, max(St00) + (max(St00)- min(St00))*0.15), 
     xaxs="i", yaxs="i", main=Main)
axis(side=2, labels=paste0(seq(-30, 30, by=1), "%"), at=seq(-30, 30, by=1), las=2)
abline(v=seq(0, nrow(SPY04), length.out = 13), col="grey", lty=3, lwd=0.5)
abline(h=seq(-30, 30, by=1), col="black", lty=1, lwd=0.3)
lines(St, col=COL, lwd=1.2)
plotrix::boxed.labels(seq(0, nrow(SPY04), length.out = 13)[1:12]+10,  min(St00) - (max(St00)- min(St00))*0.05, 
                      month.abb, cex=0.7, bg = "grey20", xpad = 1.2, ypad = 1.2)
legend("topleft", legend="SPDR S&P500 ETF", col=COL, lwd=1, cex=0.7)
quartz.save(file = paste0("./SPY_03.png"), type = "png", dpi = 300); dev.off()

#図 4. 大統領選翌年 SPYのSeasonality
#'97 '01 '05 '09 '13 '17
Main <- "図 4. 大統領選翌年 SPYのSeasonality"
Y <- c(5, 9, 13, 17, 21, 25)
head(SPY04[,c(Y)])
#      Y1997     Y2001     Y2005      Y2009       Y2013     Y2017
#1 0.0000000 0.0000000  0.000000  0.0000000  0.00000000 0.0000000
#2 1.4352047 4.8034934 -1.221946 -0.1183315 -0.22593592 0.5949196
#3 0.5487547 3.6754003 -1.903575  0.5486252  0.21223949 0.5150013

SPY04$Mean <- NA
SPY04$Mean <- apply(data.frame(SPY04[,c(Y)]), 1, function(x) mean(x[!is.na(x)]))
colnames(SPY04)

COL <- "green1"
St <- SPY04$Mean; St00 <- St[!is.na(St)]
par(family="HiraKakuProN-W3", lwd=1, xpd=F, cex=1, mgp=c(0.5, 1, 0), mai=c(0.5, 0.75, 0.5, 0.5))
plot(St, type="n", axes = F, xlab="© 2021 京橋のバイオインフォマティシャンの日常 by skume", 
     ylab="", cex.lab=0.75, xlim=c(0,nrow(SPY04)), 
     ylim=c(min(St00) - (max(St00)- min(St00))*0.1, max(St00) + (max(St00)- min(St00))*0.15), 
     xaxs="i", yaxs="i", main=Main)
axis(side=2, labels=paste0(seq(-30, 30, by=1), "%"), at=seq(-30, 30, by=1), las=2)
abline(v=seq(0, nrow(SPY04), length.out = 13), col="grey", lty=3, lwd=0.5)
abline(h=seq(-30, 30, by=1), col="black", lty=1, lwd=0.3)
lines(St, col=COL, lwd=1.2)
plotrix::boxed.labels(seq(0, nrow(SPY04), length.out = 13)[1:12]+10,  min(St00) - (max(St00)- min(St00))*0.05, 
                      month.abb, cex=0.7, bg = "grey20", xpad = 1.2, ypad = 1.2)
legend("topleft", legend="SPDR S&P500 ETF", col=COL, lwd=1, cex=0.7)
quartz.save(file = paste0("./SPY_04.png"), type = "png", dpi = 300); dev.off()

参考資料

stackoverflow.com

rstudio.github.io

docs.ruby-lang.org

stackoverflow.com

*1:じっちゃま曰く、阿波踊りして、降りなかったヤツが悪いのですがね

*2:自由に転載可です