R

原油価格推移

gdata、xts、lattice、tseries、vars、MSBVAR、urca、lmtest パッケージ

IMFのデータ(月次)
1
2
3
4
5
6
7
8
9
10
11
12
library(gdata)
#IMFのデータ
d<-read.xls("http://www.imf.org/external/np/res/commod/External_Data.xls",skip=6)
#原油価格の項目を抽出
oil<-d[,c(1,36:39)]
names(oil)<-c("year_month","average","Brent","Dubai","WTI")
oil.ts<-ts(oil[,-1], start=c(1980,1), frequency=12)
#head(oil.ts);tail(oil.ts)
#png("imf01.png",width=1000,height=800)
plot(oil.ts[,2:4],plot.type="single",col=1:3,las=1,ylab="原油価格（ドル）",main="原油価格(Brent,Dubai,WTI)の推移 単位：ドル")
legend("topleft",colnames(oil.ts[,2:4]),lty=c(1,1),lwd=c(2.5,2.5),col=1:3) 
#dev.off()

2010年以降のグラフ

1
2
3
4
png("imf02.png",width=1000,height=800)
plot(window(oil.ts[,2:4],c(2010,1),c(2014,12)),plot.type="single",col=1:3,las=1,ylab="原油価格（ドル）",main="原油価格(Brent,Dubai,WTI)の推移 単位：ドル")
legend("topleft",colnames(oil.ts[,2:4]),lty=c(1,1),lwd=c(2.5,2.5),col=1:3) 
dev.off()

EIAのデータ(日次) WTI,Brent

EIA

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
library(xts)
library(gdata)
library(lattice)
d1<-read.xls("http://www.eia.gov/dnav/pet/hist_xls/RCLC1d.xls",sheet=2,skip=3,header=F)
d2<-read.xls("http://www.eia.gov/dnav/pet/hist_xls/RBRTEd.xls",sheet=2,skip=3,header=F)
wti<-d1;brent<-d2
#head(wti);head(brent)
names(wti)<-c("date","wti")
names(brent)<-c("date","brent")
Sys.setlocale("LC_TIME", "C")
wti$date<-strptime(wti$date,format ="%b %d%Y")
brent$date<-strptime(brent$date,format ="%b %d%Y")
#mergeするall=Tをつけること
oil <- merge(wti,brent, all=T)
oil.xts <- as.xts(zoo(oil[,-1]), as.POSIXct(oil[,1]))
#install.packages("xtsExtra", repos="http://R-Forge.R-project.org")
#library(xtsExtra)
#plot.xts(oil.xts, screens = factor(1, 1), auto.legend = TRUE)
#png("oil01.png",width=1000,height=800)
xyplot(oil.xts, superpose=TRUE, xlab="Year", ylab="USD",main="Crude Oil Prices")
#dev.off()

1
2
3
4
5
#月平均
#apply.monthly(oil.xts,mean,na.rm=T)
#png("oil02.png",width=1000,height=800)
xyplot(apply.monthly(oil.xts,mean,na.rm=T), superpose=TRUE, xlab="Year", ylab="USD",main="Crude Oil Prices")
#dev.off()

1
2
3
4
#2011年以降(日次)
#png("oil03.png",width=1000,height=800)
xyplot(oil.xts['2011-01-01::'], superpose=TRUE, xlab="Year", ylab="USD",main="Crude Oil Prices")
#dev.off()

2011年以降(日次)についてGrangerの因果検定を行ってみる。
（注意）正しい方法かどうかはわかりません。

1
2
3
4
5
6
7
8
#欠損値を削除。
x<-na.omit(data.frame(coredata(oil.xts['2011-01-01::'])))
#ccf(x[,1],x[,2])
library(vars)
library(tseries)
# (1) 個々の原系列に対して単位根検定を実施
adf.test(x$wti);adf.test(x$brent)
adf.test(diff(x$wti,lag=1));adf.test(diff(x$brent,lag=1))

原系列はどちらも単位根過程。
１階差分をとれば単位根なし

1
2
# (2)共和分関係を推定  
po.test(x)

Phillips-Ouliaris Cointegration Test

data: x
Phillips-Ouliaris demeaned = -14.6339, Truncation lag parameter = 10, p-value = 0.15

1
2
oil.vecm<-ca.jo(x,ecdet="none",type="eigen",K=2,spec="longrun")
summary(oil.vecm)

Values of teststatistic and critical values of test:

test 10pct 5pct 1pct
r <= 1 | 0.84 6.50 8.18 11.65
r = 0 | 7.45 12.91 14.90 19.19

  • どちらの検定でも共和分なし
1
2
3
4
5
6
7
8
# (3) 差分系列に対してVARモデルを推定して因果性分析
x.d<-data.frame(wti=diff(x$wti,lag=1),brent=diff(x$brent,lag=1))
VARselect(x.d)
#VAR Model, lag=1
V.1<-VAR(x.d,p=1,type="both")
#VAR-Model, lag=3
V.3<-VAR(x.d,p=3,type="both")
serial.test(V.1);serial.test(V.3)
  • lag=3を採用
1
2
3
# (4) Grangerの因果検定
library(MSBVAR)
granger.test(x.d,p=3)

F-statistic p-value
brent -> wti 0.8288626 4.780541e-01
wti -> brent 19.6839451 2.105205e-12

  • wti -> brent のGrangerの因果あり。
1
2
3
4
5
6
7
library(lmtest)
#wti -> brent
grangertest(x.d$wti,x.d$brent, order = 3)
grangertest(brent ~ wti, order = 3, data = x.d)
#brent -> wti
grangertest(x.d$brent,x.d$wti, order = 3)
grangertest(wti ~ brent, order = 3, data = x.d)
  • wti -> brent のGrangerの因果あり。
1
2
3
4
5
6
7
8
9
#FIAR
library(Matrix)
library(np)
#library(lavaan) 
library(devtools)
#FIAR:FIARパッケージはcranにはもうないが、githubに関数があったので読み込んで使う。
#source_url('https://raw.githubusercontent.com/cran/FIAR/master/R/ARorder.R')
source_url('https://raw.githubusercontent.com/cran/FIAR/master/R/partGranger.R')
#ARorder(x.d)

partGranger関数は偏Granger因果検定のための関数のようですが、Granger因果検定でも使えます。

1
2
3
4
#wti -> brent
partGranger(x.d[,c("wti","brent")],nx=1,ny=1,order=3)
#brent -> wti
partGranger(x.d[,c("brent","wti")],nx=1,ny=1,order=3)

partGranger(x.d[,c(“wti”,”brent”)],nx=1,ny=1,order=3)

$orig
0.05726271

$prob
2.105237e-12

partGranger(x.d[,c(“brent”,”wti”)],nx=1,ny=1,order=3)

$orig
0.00247855

$prob
0.4780541

  • wti -> brent のGrangerの因果あり。