Feb 14 2015

東大物価指数・USDJPY・日経平均株価・WTI原油価格

zoo,xts,TTR,gdata,quantmod パッケージ

東大日次物価指数プロジェクトというサイトがあります。
そこでは、東大日次物価指数のグラフとデータが公開されています。
(詳しくはサイトでご覧ください。）

ここでは、Rを使って、東大物価指数・USDJPY・日経平均株価・WTI原油価格のデータの取得とグラフの作成をしてみます。
2015年１月10日現在、以下のコードでうまくいきました。
（サイト消滅、移転、データ公開取りやめなどの理由によりNGになることもあり）

データは
東大物価指数は上記。
USDJPY : 日銀のHP
日経平均株価 : yahooファイナンス

日経新聞社のサイトのデータは正確だけど注意書きが・・・
WTI原油価格 : Energy Information Administration
より入手。感謝です。

データを入手。必要な部分だけ取り出してxtsクラスに変換。

library(zoo)
library(xts)
library(TTR)
library(gdata)
#東大日次物価指数プロジェクト
prices <- read.csv("http://www.cmdlab.co.jp/price_u-tokyo/download/HistoricalDailyData.csv", header=TRUE,na.strings="NA")
prices<-prices[,1:2]
#####
#USDJPY:日銀のHP
USDJPY<-read.csv("http://www.stat-search.boj.or.jp/ssi/mtshtml/csv/d.csv",skip=6,header=F, na.strings="NA")
USDJPY<-data.frame(date=USDJPY[,1],dy=USDJPY[,5])
prices.xts <- as.xts(zoo(prices[,-1]), as.POSIXct(prices$date))
USDJPY.xts <- as.xts(zoo(USDJPY[,-1]), as.POSIXct(USDJPY$date))
#####
#日経平均株価:yahooファイナンス
market<-read.csv("http://real-chart.finance.yahoo.com/table.csv?s=%5EN225&d=0&e=8&f=2015&g=d&a=0&b=4&c=1984&ignore=.csv",skip=0,header=T, na.strings="NA")
nikkei<-market[,c(1,5)]
nikkei.xts <- as.xts(zoo(nikkei[,-1]), as.POSIXct(nikkei$Date))
#グラフにするとおかしいところがある
#調べてみる
subset(nikkei.xts, subset=nikkei.xts<=5000)
#2014/12/1のデータがおかしい。
#nikkei.xts["2014-12-01"]
#2014-12-01 175.9
#正しくは、17590.10
#訂正する
nikkei.xts["2014-12-01"]<-17590.10
#####
#WTI原油価格 Energy Information Administration 
#zorin OS なので fileEncoding="CP932" をつけて、
d<-read.xls("http://www.eia.gov/dnav/pet/hist_xls/RWTCd.xls",sheet=2,skip=2,fileEncoding="CP932" )
names(d)<-c("date","WTI")
wti<-d
#head(wti)
#「%d」「%m」「%Y」以外のフォーマットを使う場合，
#日本語環境では数値の代わりに日本語の文字列が入ることにより NA になってしまう．
#事前に関数 Sys.setlocale("LC_TIME","C") を実行する必要がある．
Sys.setlocale("LC_TIME", "C")
wti$date<-strptime(d$date,format ="%b %d, %Y")
wti.xts <- as.xts(zoo(wti[,-1]), as.POSIXct(wti[,1]))
#いつまでのデータが入手できたのか確認
tail(prices.xts);tail(USDJPY.xts);tail(nikkei.xts);tail(wti.xts)

2015-01-02 -0.001627
2015-01-03 -0.004865
2015-01-04 -0.007219
2015-01-05 0.000507
2015-01-06 -0.007007
2015-01-07 -0.000212

2015-01-02 NA
2015-01-03 NA
2015-01-04 NA
2015-01-05 120.29
2015-01-06 119.30
2015-01-07 118.95

2014-12-31 17450.77
2015-01-02 17450.77
2015-01-05 17408.71
2015-01-06 16883.19
2015-01-07 16885.33
2015-01-08 17167.10

2014-12-26 54.59
2014-12-29 53.46
2014-12-30 54.14
2014-12-31 53.45
2015-01-02 52.72
2015-01-05 50.05

ここでは最終データはそろえずにグラフ化します。(目立たないだけでグラフのx軸はそろってはいません）

グラフにする最初のデータは2012-01-01ですが、NAになっているデータもあります。

重要イベント（ここでは４つ）に目印をつける準備

#自公連立政権:2012年12月26日
ab<-"2012-12-26"
#消費税８％
zei<-"2014-04-01"
#日銀の追加緩和:黒田バズーカ第一弾 ２０１３年４月４日
#日銀の追加緩和:黒田バズーカ第二弾 ２０１４年１０月３１日
bazooka<-c("2013-04-04","2014-10-31")

グラフ１（東大物価指数・WTI原油価格円換算・WTI原油価格ドル）

#png("prices01.png", width=1000, height=800)
par(mfrow=c(3,1),mar = c(3.2, 5, 1.8, 3))
plot(prices.xts["2012-01-01::"],main="東大物価指数(2012-01-01～)")
abline(h = 0, col = "blue", lty =2)
abline(v = as.POSIXct(zei),col = "red", lty = 1)
abline(v = as.POSIXct(ab),col = "red", lty = 1)
abline(v = as.POSIXct(bazooka),col = "red", lty = 2)
#7日移動平均線(欠損値なし）
lines(SMA(prices.xts, 7),col="red")
text(as.POSIXct("2012-12-26"),0.02,"自公政権",pos=4,col="red",cex=1,font = 2)
text(as.POSIXct(zei),0.02,"消費税：８％",pos=4,col="red",cex=1,font = 2)
text(as.POSIXct("2013-04-04"),0.02,"バズーカ:'13-04-04",pos=4,col="red",cex=1,font = 2)
text(as.POSIXct("2014-10-31"),0.02,"バズーカ:'14-10-31",pos=4,col="red",cex=1,font = 2)
#WTI原油価格にUSDJPYをかける（円換算）
plot(wti.xts["2012-01-01::"]*USDJPY.xts["2012-01-01::"],main="WTI原油価格 円換算(2012-01-01～)")
abline(v = as.POSIXct(zei),col = "red", lty = 1)
abline(v = as.POSIXct(ab),col = "red", lty = 1)
abline(v = as.POSIXct(bazooka),col = "red", lty = 2)
#1バレル9000円に横線をいれてみる
abline(h = 9000, col = "blue", lty =2)
plot(wti.xts["2012-01-01::"],main="WTI原油価格 ドル(2012-01-01～)")
abline(v = as.POSIXct(zei),col = "red", lty = 1)
abline(v = as.POSIXct(ab),col = "red", lty = 1)
abline(v = as.POSIXct(bazooka),col = "red", lty = 2)
#1バレル79ドルに横線をいれてみる
abline(h = 79, col = "blue", lty =2)
par(mfrow=c(1,1))
#dev.off()

一番上のグラフの横線は上昇率０％
真ん中のグラフの横線は1バレル9000円。偶然にもバズーカは2回ともこの水準。
一番下のグラフの横線は1バレル79ドル

グラフ2（東大物価指数・アメリカドル / 日本円・日経平均株価・日経平均株価ドル換算）

#png("prices02.png", width=1000, height=800)
par(mfrow=c(4,1),mar = c(3.2, 5, 1.8, 3))
plot(prices.xts["2012-01-01::"],main="東大物価指数(2012-01-01～)")
abline(h = 0, col = "blue", lty =2)
abline(v = as.POSIXct(zei),col = "red", lty = 1)
abline(v = as.POSIXct(ab),col = "red", lty = 1)
abline(v = as.POSIXct(bazooka),col = "red", lty = 2)
#7日移動平均線(欠損値なし）
lines(SMA(prices.xts["2012-01-01::"], 7),col="red")
text(as.POSIXct("2012-12-26"),0.018,"自公政権",pos=4,col="red",cex=0.8,font = 2)
text(as.POSIXct(zei),0.018,"消費税：８％",pos=4,col="red",cex=0.8,font = 2)
text(as.POSIXct("2013-04-04"),0.018,"バズーカ:'13-04-04",pos=4,col="red",cex=0.8,font = 2)
text(as.POSIXct("2014-10-31"),0.018,"バズーカ:'14-10-31",pos=4,col="red",cex=0.8,font = 2)
plot(USDJPY.xts["2012-01-01::"],main="アメリカ ドル / 日本 円(2012-01-01～)")
abline(v = as.POSIXct(zei),col = "red", lty = 1)
abline(v = as.POSIXct(ab),col = "red", lty = 1)
abline(v = as.POSIXct(bazooka),col = "red", lty = 2)
plot(nikkei.xts["2012-01-01::"],main="日経平均株価(2012-01-01～)")
abline(v = as.POSIXct(zei),col = "red", lty = 1)
abline(v = as.POSIXct(ab),col = "red", lty = 1)
abline(v = as.POSIXct(bazooka),col = "red", lty = 2)
plot(nikkei.xts["2012-01-01::"]/USDJPY.xts["2012-01-01::"],main="日経平均株価 ドル換算(2012-01-01～)")
abline(v = as.POSIXct(zei),col = "red", lty = 1)
abline(v = as.POSIXct(ab),col = "red", lty = 1)
abline(v = as.POSIXct(bazooka),col = "red", lty = 2)
par(mfrow=c(1,1))
#dev.off()

日経平均株価 (ドル換算)はバズーカ１では上昇しているが、バズーカ２ではいまのところ上昇はしていない。

グラフ３：週平均のグラフを作成
グラフ４：移動平均を加える

#各週ごとの平均値
#png("prices03.png", width=1000, height=800)
par(mfrow=c(4,1),mar = c(3.2, 5, 1.8, 3))
plot(apply.weekly(prices.xts["2012-01-01::"],mean,na.rm=TRUE),main="東大物価指数 週平均 (2012-01-01～)")
abline(h = 0, col = "blue", lty =2)
abline(v = as.POSIXct(zei),col = "red", lty = 1)
abline(v = as.POSIXct(ab),col = "red", lty = 1)
abline(v = as.POSIXct(bazooka),col = "red", lty = 2)
text(as.POSIXct("2012-12-26"),0.005,"自公政権",pos=4,col="red",cex=1,font = 2)
text(as.POSIXct(zei),0.005,"消費税：８％",pos=4,col="red",cex=1,font = 2)
text(as.POSIXct("2013-04-04"),0.005,"バズーカ:'13-04-04",pos=4,col="red",cex=1,font = 2)
text(as.POSIXct("2014-10-31"),0.005,"バズーカ:'14-10-31",pos=4,col="red",cex=1,font = 2)
plot(apply.weekly(USDJPY.xts["2012-01-01::"],mean,na.rm=TRUE),main="アメリカ ドル / 日本 円 週平均 (2012-01-01～)")
abline(v = as.POSIXct(zei),col = "red", lty = 1)
abline(v = as.POSIXct(ab),col = "red", lty = 1)
abline(v = as.POSIXct(bazooka),col = "red", lty = 2)
plot(apply.weekly(nikkei.xts["2012-01-01::"],mean,na.rm=TRUE),main="日経平均株価 週平均 (2012-01-01～)")
abline(v = as.POSIXct(zei),col = "red", lty = 1)
abline(v = as.POSIXct(ab),col = "red", lty = 1)
abline(v = as.POSIXct(bazooka),col = "red", lty = 2)
plot(apply.weekly(wti.xts["2012-01-01::"]*USDJPY.xts["2012-01-01::"],mean,na.rm=TRUE),main="WTI原油価格  円換算 週平均 (2012-01-01～)")
abline(v = as.POSIXct(zei),col = "red", lty = 1)
abline(v = as.POSIXct(ab),col = "red", lty = 1)
abline(v = as.POSIXct(bazooka),col = "red", lty = 2)
#1バレル9000円に横線をいれてみる
abline(h = 9000, col = "blue", lty =2)
par(mfrow=c(1,1))
#dev.off()
###
###
#移動平均（注：休場のときは線形補間を行います）
#na.fill(USDJPY.xts,"extend"):先頭、終端であれば、NA ではない値をそのまま延長、中間であれば線形補間を行う
#png("prices04.png", width=1000, height=800)
par(mfrow=c(3,1),mar = c(3.2, 5, 1.8, 3))
plot(USDJPY.xts["2012-01-01::"],main="アメリカ ドル / 日本 円(2012-01-01～)")
abline(v = as.POSIXct(zei),col = "red", lty = 1)
abline(v = as.POSIXct(ab),col = "red", lty = 1)
abline(v = as.POSIXct(bazooka),col = "red", lty = 2)
#25日移動平均線(休場（欠損値）があるため工夫が必要)
lines(SMA(na.fill(USDJPY.xts["2012-01-01::"],"extend"),25),col="red")
plot(nikkei.xts["2012-01-01::"],main="日経平均株価(2012-01-01～)")
abline(v = as.POSIXct(zei),col = "red", lty = 1)
abline(v = as.POSIXct(ab),col = "red", lty = 1)
abline(v = as.POSIXct(bazooka),col = "red", lty = 2)
#25日移動平均線
lines(SMA(na.fill(nikkei.xts["2012-01-01::"],"extend"),25),col="red")
#日経平均株価をUSDJPYで割る（ドル換算）
plot(nikkei.xts["2012-01-01::"]/USDJPY.xts["2012-01-01::"],main="日経平均株価 ドル換算(2012-01-01～)")
abline(v = as.POSIXct(zei),col = "red", lty = 1)
abline(v = as.POSIXct(ab),col = "red", lty = 1)
abline(v = as.POSIXct(bazooka),col = "red", lty = 2)
#25日移動平均線
lines(SMA(na.fill(nikkei.xts["2012-01-01::"],"extend")/na.fill(USDJPY.xts["2012-01-01::"],"extend"),25),col="red")
par(mfrow=c(1,1))
#dev.off()

備考

Fetching Public Data With R
パブリックデータの取り込みを助けるパッケージを紹介している。

quantmodパッケージを使うと上の方法より簡単でした。
（東大日次物価指数の取り込みには使えない）

library(quantmod)
#このパッケージを使うと読み込んだ時点で、xts classになっている
getSymbols("^N225",from="2012-01-01")
head(N225);tail(N225)
plot(N225)
#グラフにすると明らかにおかしいところがある
#調べてみる
subset(N225, subset=N225$N225.Close<=5000)
#N225.Open N225.High N225.Low N225.Close N225.Volume N225.Adjusted
#2014-12-01    174.75    176.49   174.74      175.9           0         175.9
#正しくは
#日付	始値	高値	安値	終値
#2014/12/01	17,475.10	17,649.02	17,474.27	17,590.10
#訂正する(小数点以下を無視してもいいなら１００倍する)
N225["2014-12-01",c(1:6)]<-c(17475.10,17649.02,17474.27,17590.10,0,17590.10)
#chartSeriesという関数できれいなチャートが描ける
#chartSeries(N225)
#We can retrieve foreign exchanges or precious metal data from oanda. 
#Beware that there is a limit of 500 days per request.
#分けて取り込み、結合する。
getFX("USD/JPY",from="2012-01-01", to = '2012-12-31')
d1<-USDJPY
getFX("USD/JPY",from="2013-01-01", to = '2013-12-31')
d2<-USDJPY
getFX("USD/JPY",from="2014-01-01")
d3<-USDJPY
USDJPY<-rbind(d1,d2,d3)
tail(USDJPY)
#FREDは更新が遅い
#getSymbols("EXJPUS",src="FRED")
#tail(USDJPY);tail(EXJPUS)
#chartSeries(USDJPY)
#WTI crude prices(WTI原油価格)
getSymbols('DCOILWTICO',src='FRED')
tail(DCOILWTICO) 
#chartSeries(DCOILWTICO["2012-01-01::"])
#次のようにしてもよいが
#最終データの古さ、変換の手間を考えるとおすすめできない
#library(Quandl)
#WTI<-Quandl("DOE/RWTC")
#class(WTI) #xts classに変換する必要有り

2014年4月１日から２０１５年１月11日までのグラフを描いてみます。

xlimオプションを使ってx軸を揃えます

start<-"2014-04-01" ; end<-"2015-01-11"
r<-as.POSIXct(c(start,end))
t<-paste(start,"::",end,sep="")
#png("qm01.png", width=1000, height=800)
par(mfrow=c(3,1),mar = c(3.2, 5, 1.8, 3))
#N225は終値
plot(N225$N225.Close[t],main=paste("日経平均株価(",start,"～",end,")"),xlim=c(r[1],r[2]))
plot(USDJPY[t],main=paste("アメリカ ドル / 日本 円(",start,"～",end,")"),xlim=c(r[1],r[2]))
plot(DCOILWTICO[t]*USDJPY["2012-01-01::"],main=paste("WTI原油価格 円換算(",start,"～",end,")"),xlim=c(r[1],r[2]))
par(mfrow=c(1,1))
#dev.off()