Quandmodで金融データをとってStanで同時相関の変化を見る
中盤はテクニカルな話題なので、統計や経済に馴染みのない人は読み飛ばして最後だけ読んでください
金融関係のレポートを見ていると、よく相関係数の変化を示すグラフが登場します。
(zero hedge http://www.zerohedge.com/news/2014-12-30/spot-difference)
（日銀「最近の株価と為替の同時相関関係の強まりについて」）
だいたいの場合において、こうした相関係数はその地点の前100日〜500日くらいの期間をとってヒストリカルな相関（sample correlation）をローリングして可視化するだけといった簡単なものになっています。
統計的推定の文脈でいうと、これはfiltering（過去のデータのみを用いてその地点の真の値を推定する）で得られる値に近く、smoothing（後のデータも加味して後だしじゃんけんで推定する）ではありません。
予測の文脈では過去のデータのみしか使えないのでfilteringが大事なのですが、だいたいこういう相関係数は「いつどの指数同士がどう相関していたか」という振り返りに使われているので、smoothingのほうが適切です。もちろん「smoothingっぽくするには前後のx日をとればいい」ということにもなりますが、期間xのとりかた如何で相関係数の変化の具合も変わるので、もっと洗練された方法はないのか、という話になります。
ということで、相関係数ρ自体がランダムウォークで時間変化する非線形状態空間モデルをStanで推定してしまいましょう。
モデルとしてはこんな感じです。
t期の各系列の対数リターン(y_t1,y_t2)が正規分布に従うという典型的なランダムウォークですが、ノイズの部分が相関を持つというものです。こうした構造は統計用語ではSeemingly Unrelated(SUR model, 見かけ上無相関のモデル)という用語で呼ばれます。
さらにこの相関係数ρをランダムウォークさせてしまおうというのが今回の目論みです。ここまでモデルが複雑になると、一般に太刀打ちが難しい（世の中に存在するパッケージソフトやライブラリでは到底推定不可能）のですが、StanはMCMC(マルコフ連鎖モンテカルロ)という汎用的なソリューションを用いているので、モデルを記述するだけで解決してくれます。
ちなみにこれは、2変量確率的ボラティリティ変動モデル(Bivariate Stochastic Volatility Model)の特殊ケースにあたります。
相関係数は-1〜1という制約があるので、ランダムウォークしてもそこに収まるようにうまくロジスティック関数で変換しています。
（各系列の分散の変化も考える場合には、コレスキー分解済みの三角行列の各要素をランダムウォークさせて、それの積をとって共分散行列をつくったほうが楽かもしれません。）
【もととなるデータの取得】
Rのライブラリquantmodが楽です。データソースの指定に関しては、株はYahooやGoogle,FXはOandaなど、ほかのソースからもとってくることが出来ますが、FREDhttp://research.stlouisfed.org/fred2/(FRBが管理している経済データ集)がおすすめです。
library(vars)
library(quantmod)
library(dplyr)
library(car)
library(ggplot2)
library(reshape2)
library(PerformanceAnalytics)
data(Canada)master_quote=read.csv("./master_quotes.csv",stringsAsFactors=F)
dat_all=list()
daily_return=list()
com_indx=c()
daily_return_bundle=list()
close_bundle=list()get_data=function(quote,src){
if(src=="yahooj"){
tmp_dat=data.frame(na.omit(rfj(quote,start="2000-1-1")))
ret=tmp_dat["Close"]
rownames(ret)=as.Date(gsub(" JST","",tmp_dat$Date))
colnames(ret)=quote
r=as.xts(ret)
index(r)=as.Date(gsub(" JST","",tmp_dat$Date))
r
}else{
na.omit(getSymbols(quote,src=src,auto.assign=F))
}
}for(i in 1:nrow(master_quote)){
quote=master_quote[i,"quote"]
src=master_quote[i,"src"]
adj=master_quote[i,"adj"]
dat_allquote=get_data(quote,src)
index(dat_allquote)=index(dat_allquote)+adj
daily_returnquote=dailyReturn(dat_allquote,type="log")
if(i==1){
com_indx=index(daily_returnquote)
}else{
com_indx=intersect(com_indx,index(daily_returnquote))
}
}
for(i in 1:nrow(master_quote)){
quote=master_quote[i,"quote"]
src=master_quote[i,"src"]
daily_return_bundlequote=daily_returnquote[as.Date(com_indx)]close_bundlequote=if(src=="oanda"||src=="FRED"||src=="yahooj"){dat_allquote[index(dat_allquote)%in%as.Date(com_indx)]}else{Cl(dat_allquote)[as.Date(com_indx)]}
colnames(daily_return_bundlequote)=quote
colnames(close_bundlequote)=quote
}
drb.df=data.frame(daily_return_bundle)
close.df=data.frame(close_bundle)drb.ts=ts(drb.df)
close.ts=ts(close.df)roll.cor = function(x) {
cor(x)[1,2]
}pdf("test.pdf",family="Japan1")
chart.Correlation(drb.df)
cnames=colnames(drb.df)
for(cnum1 in 1:length(cnames)){
for(cnum2 in 1:length(cnames)){
if(cnum1<cnum2){
cname1=cnames[cnum1]
cname2=cnames[cnum2]
drb.df_filter=drb.df[c(cname1,cname2)]
roll.cor.vals = rollapply(as.zoo(drb.df_filter), width=300,
by.column=FALSE,
FUN=roll.cor,
align="right")
roll.cor.vals.df=data.frame(roll.cor.vals)
correlation=roll.cor.vals.df[rownames(close.df),]
close.df_selected=close.df[c(cname1,cname2)]
colnames(close.df_selected)=c(master_quote[master_quote$quote==cname1,"name"],master_quote[master_quote$quote==cname2,"name"])
d=cbind(correlation,close.df_selected)
d$date=rownames(d)
d_plot=melt(d,id="date")
d_plot$date=as.Date(d_plot$date)
g=ggplot(d_plot,aes(date, value, col=variable, group=1))+
geom_line()+
facet_grid(variable~., scale='free_y')
plot(g)
}
}
}dev.off()
write.csv(drb.df,file="drb.df.csv")
write.csv(close.df,file="close.df.csv")
参照しているcsvデータは以下になります。
quote,src,adj,name
NIKKEI225,FRED,0,日経平均
NASDAQCOM,FRED,0,ナスダック
SP500,FRED,0,S&P500
DCOILWTICO,FRED,0,WTI原油
DEXJPUS,FRED,0,ドル円(FED)
VIXCLS,FRED,0,VIX
TEDRATE,FRED,0,TEDスプレッド
USD6MTD156N,FRED,0,6ヵ月LIBOR
USEPUINDXD,FRED,0,政策の不確実性
GOLDAMGBD228NLBM,FRED,0,金(LONDON)
DTWEXM,FRED,0,ドルインデックス
DEXUSAL,FRED,0,米ドル豪ドル(FRED)
実行すると、ローリングでの相関係数の変化がまとめてpdfファイルに出力されます。
【Stanによる推定】
ただ、これだと物足りないのでStanで相関係数ρの変化をモデリングします。
ここでは日経平均とドル円の変化率（対数収益率）の同時相関を上記のモデルで捉えます。データの数が多いと収束が怪しいので、10営業日ごとにリターンをまとめたものを用います。
Stanのモデルは以下になります。
data {
int<lower=1> T;
vector[2] y[T];
}parameters {
real<lower=0.0001> sigma[2];
real rho_logistic[T];
real mu_rho;
real<lower=-0.999,upper=0.999> phi_rho;
real<lower=0.0001> eta_rho;
}
transformed parameters {
cov_matrix[2] Sigma[T];
real rho[T];
vector[2] mu;
for(t in 1:T){
rho[t] <- 2/(1+exp((-1)*rho_logistic[t]))-1;
Sigma[t][1,1]<-sigma[1]*sigma[1];
Sigma[t][1,2]<-sigma[1]*sigma[2]*rho[t];
Sigma[t][2,1]<-sigma[1]*sigma[2]*rho[t];
Sigma[t][2,2]<-sigma[2]*sigma[2];
}
mu[1]<-0;
mu[2]<-0;
}
model {
y[1]~multi_normal(mu',Sigma[1]);
rho_logistic[1]~normal(mu_rho,eta_rho*eta_rho/(1-phi_rho*phi_rho));
for (t in 2:T){
rho_logistic[t]~normal(mu_rho*(1-phi_rho)+phi_rho*rho_logistic[t-1],eta_rho*eta_rho);
y[t]~multi_normal(mu',Sigma[t]);
}
}
Rでキックするコードは以下になります。
library(rstan)
library(reshape2)
library(ggplot2)
library(dplyr)model.text=scan("sur.rstan",what = character(),sep="@" , blank.lines.skip = F)
model.text=paste(model.text,collapse="\n")
cnames=c("NIKKEI225","DEXJPUS")
by=10dat=read.csv("drb.df.csv",stringsAsFactor=F)
close.df=read.csv("close.df.csv",stringsAsFactors=F)
dat=cbind(aggregate(.~floor((1:nrow(dat)/by)), dat[cnames], FUN=sum),aggregate(.~floor((1:nrow(dat)/by)), dat[c("X")], FUN=last))d=list()
d[["T"]]=nrow(dat)
d[["y"]]=dat[cnames]*100fitchan <- stan(model_code = model.text,
data = d,
iter = 3000, warmup=500,
chain = 1)
traceplot (fitchan,par=c("sigma"),inc_warmup=F)fit=extract(fitchan,permuted=T)
fit.posterior=get_posterior_mean(fitchan)# From Small Data Scientist Memorandum
# (http://heartruptcy.blog.fc2.com/blog-entry-163.html)
plot_timecourse <- function(la_m, time_scale){
qua <- apply(la_m, 2, quantile, prob = c(0.025, 0.5, 0.975))
est <- data.frame(t=time_scale, t(qua))
colnames(est) <- c('t', 'p2.5', 'p50', 'p97.5')
est_med <- data.frame(t=time_scale, value=apply(la_m, 2, mean))
p <- ggplot()
p <- p + geom_ribbon(data=est, aes(x=t, ymin=p2.5, ymax=p97.5), fill="darkorange3", alpha=I(1/3))
p <- p + geom_line(data=est_med, aes(x=t, y=value), color="red")
p <- p + labs(x="time", y="value")
p
}pdf("time-series_correlation_SUR_model.pdf",font="Japan1")
plot(plot_timecourse(data.frame(fit$rho),as.Date(dat$X)))
d=data.frame(correlation=apply(fit$rho,2,mean),close.df[close.df$X%in%dat$X,cnames])
d$date=dat$X
d_plot=melt(d,id="date")
d_plot$date=as.Date(d_plot$date)g=ggplot(d_plot,aes(date, value, col=variable, group=1))+
geom_line()+
facet_grid(variable~., scale='free_y')
plot(g)
ggdat=melt(fit["sigma"])
graph <- ggplot (ggdat, aes(x = value))+geom_histogram(aes(y=..density..),colour="blue", fill="grey",alpha=.1)
graph <- graph + geom_density(alpha=.5,fill="blue") +facet_wrap(Var2~L1,ncol=1,scale="free")
plot(graph)dev.off()
結果は以下になります。
相関係数の変化のグラフは次のようになります。
90%信用区間が薄い色がついたところ、事後平均が赤い線です。
原系列とならべると次のようになります。
ポートフォリオ理論だと相関係数は一定と仮定することが多く、普通は期間をきめてヒストリカルをとることが多いですが、時期による相関係数の変化は無視できないほど大きく、しかも急に変動することがあるということがわかりました。
一方で、90%信頼区間ベースではかなり推定値が広がってしまう、つまりきっちり推定するのが難しいということもわかりました。
グラフからは、長期的に見て負相関（円高株高）から正相関（円安株高）へと変化してきていて、例えば2009年では一時的にその傾向が弱まったということが言えそうです。
円安株高の同時相関の高まりに関しては、日本株に投資する海外投資家がヘッジのために立てている円ポジションのリバランスが要因なのではないかという説が有力です（上述日銀ペーパーなど）。つまり、株高になるとそれに応じてヘッジのための追加の円shortのポジションが発生するために、円売りの圧力がかかるということです。これについては海外投資家比率の高まりと、アルゴリズムトレードによって統計的裁定機会の解消速度があがったこと（＝より短期スパンで見たときにも強い同時相関が観測される）が指摘されています。
そして巷に言われるように、円安は（長期的に）輸出企業の業績を上げることを通じて株高の呼び水となり、正のフィードバックが働いているようです。
