前記事で過去のデータを使って投資ポートフォリオのリターン/リスク比を検討する方法を提案しました。本文中でも触れましたが、過去のデータはあくまで過去のデータで、短い期間からの推計や、政治経済の変化を見落とすと高値づかみになるリスクがあります。twitterでも過去のデータを使うリスクに関して同様の指摘をいただきました。過去のデータがあるからこそ未来予測ができるので過去データを全て否定するのは慎みたいところですが、未来予測についてはJP morganが各アセットの今後の期待値を出しているので、それを使ってシミュレーションをしてみました。前回に引き続き、R (RStudio)とplot.lyがメインの武器になります。
まずはデータの取得です。今回はRで自動取得する便利技は使えないので、JPMAMのWebサイトにアクセスして手動でデータを取得します。各アセットの期待リターン、期待リスク（ボラティリティ）、それぞれのアセットクラス間の相関係数が掲載されています。この三種の値があれば理論上はポートフォリオのリターンとリスクが計算できます。机上の空論、取らぬ狸の皮算用という言葉が頭をよぎりますが気にせずにファイルを作ります。テーブルにまとめるとこんな感じの気が滅入るテーブルができました。残りの作業はR上でやります。
Rに読み込んだテーブルをリターン成分、リスク成分、相関係数成分の3つに分けます。
まずはテストのポートフォリオを設定するために、新興国株、米国株、日本国債のシンプルな３つのポートフォリオでテストします。表の中でこれらはそれぞれ28, 38, 3行目にあるのでその数字から表を必要な成分だけに集約することは簡単です。
以下のサイトを参考に計算しました。
Expected Return, Variance And Standard Deviation Of A Portfolio
Portfolio Variance
Portfolio Calculations
ポートフォリオの期待リターンはウェイトと期待リターンの加重平均、期待リスクはウェイトと期待リスクのテンソル積（ちょっと難しいですがこれ）に相関係数をかけたものと表現できます。数式さえ立てれば計算は全てRがやってくれます。リターンとリスクがわかれば前回と同様にプロットして期待シャープレシオも計算できます。
ポートフォリオのウェイトを前回と同様に発生させて計算を進めていけばplot.lyのグラフが描画できます。
こんな感じになりました。（※PC/Mac推奨）3種類だとだいぶきれいなグラフになりますね。効率的フロンティアのラインもきれいに見れます。後は任意のリスク許容度からリターンを最大にする配分を選べばOKです。
応用編で、シャープレシオを最大化するアセットアロケーションも掲載します。シャープレシオを最大化するウェイトを計算する関数が公開されていたので、それで今回の3アセットクラスで計算すると、なんと米国株 0%, 新興国株 75%, 日本国債 25%という結果になりました。追加のプロットなどしたい方は前回記事を参考にしてみてください。
読書の記録をやろうやろうと思いながらなぜかあらぬ方向にブログ記事が向かっておりどんどんとマニアックな世界に入り込んでしまっていますが、どこかで修正しないと。
こんな感じのコードを書いています。
####ライブラリ読み込み#### library(dygraphs) library(dplyr) library(foreach) library(ggplot2) library(stringr) library(plotly) library(tseries) library(PortfolioOptim) library(readr)
JPMAMのデータを読み込んでリターン・リスクのデータと、相関係数のデータに分ける。
####JPMAMのデータを読み込む####
jpmam <- read_csv(file.path("data", "JPMAM.csv"), quote = ",", col_names = TRUE,
locale = locale(encoding = "UTF-8"))
asset_return <- jpmam[,c(1:5)]
asset_cor <- jpmam[,c(6:ncol(jpmam))]
ここで好きなポートフォリオを設定してみてください。1行目を実行すると各アセットクラスと対応する番号が分かります。
colnames(asset_cor) #アセットクラスと対応する番号が表示されます。今回は米国株、新興国株、日本国債の28, 38, 3を使います test <- c(28, 38, 3) #好きなアセットクラスに対応する番号を記入 num <- 1000 #プロットする数を設定 int <- 2.5 #何%刻みで測定するか1 or 2.5 or 5 or 10
以下で選んだポートフォリオからさらにテーブルを絞り込みます。
あとは前記事と同じ感じでウェイトとリスク・リターンの表を作ります。
今回はPortfolioAnalysisのパッケージは使いません。
####以下でデータを整理していきます####
return_table <- t(as.matrix(asset_return[test,3]))
var_table <- t(as.matrix(asset_return[test,4]))
cov_table <- as.matrix(asset_cor[test, test])
####全ての組み合わせのポートフォリオを再現する####
random <- combn(100/int+length(test)-1,length(test)-1)
random <- rbind(0, random , 100/int+length(test))
w <- (t(apply(random, 2, diff)) - 1) * int / 100
####ウェイトとデータからリスク（標準偏差）、リターンを求める####
temp <- NULL
sr <- t(apply(w, 1, function(i) {
temp[1] <- sqrt(sum(t(var_table * i) %x% (var_table * i) * cov_table)) * 100
temp[2] <- sum(return_table * i) * 100
sharpe_ratio_manual <- round(temp[2]/temp[1], 4)
temp[3] <- sharpe_ratio_manual[1]
temp
}))
データからplot.lyでグラフを作成します。
####データを使ってplot.lyでグラフを作成####
plotdata <- as.data.frame(sr)
rownames(plotdata) <- NULL
colnames(plotdata) <- c("StdDev", "Return", "SharpeRatio")
plot_weight <- as.data.frame(w)
rownames(plot_weight) <- NULL
colnames(plot_weight) <- colnames(cov_table)
plot_weight <- apply(plot_weight, 1, str_c, collapse = ", ")
plotdata_p <- cbind(plotdata, plot_weight)
pf_tickers <- str_c(colnames(cov_table), collapse = ", ")
p <- plot_ly(
plotdata_p, x = ~StdDev, y = ~Return,
# Hover text:
text = ~paste(pf_tickers, "\n" ,plot_weight)) %>%
layout(title = "Estimated Return and Risk from JPMAM")
p
Sys.setenv("plotly_username"="ユーザー名")
Sys.setenv("plotly_api_key"="APIキー")
api_create(p, filename = "タイトル")
ポートフォリオを最適化するウェイトはこのWeb上のPDFを参考にして関数を設定すると自動的に求められます。
####以下の2つの関数でポートフォリオを最適化するウェイトが求められます####
effFrontier = function (averet, rcov, nports = 20, shorts=T, wmax=1)
{
mxret = max(abs(averet))
mnret = -mxret
n.assets = ncol(averet)
reshigh = rep(wmax,n.assets)
if( shorts )
{
reslow = rep(-wmax,n.assets)
} else {
reslow = rep(0,n.assets)
}
min.rets = seq(mnret, mxret, len = nports)
vol = rep(NA, nports)
ret = rep(NA, nports)
for (k in 1:nports)
{
port.sol = NULL
try(port.sol <- portfolio.optim(x=averet, pm=min.rets[k], covmat=rcov,
reshigh=reshigh, reslow=reslow,shorts=shorts),silent=T)
if ( !is.null(port.sol) )
{
vol[k] = sqrt(as.vector(port.sol$pw %*% rcov %*% port.sol$pw))
ret[k] = averet %*% port.sol$pw
}
}
return(list(vol = vol, ret = ret))
}
maxSharpe = function (averet, rcov, shorts=T, wmax = 1)
{
optim.callback = function(param,averet,rcov,reshigh,reslow,shorts) {
port.sol = NULL
try(port.sol <- portfolio.optim(x=averet, pm=param, covmat=rcov,
reshigh=reshigh, reslow=reslow, shorts=shorts), silent = T)
if (is.null(port.sol)) {
ratio = 10^9
} else {
m.return = averet %*% port.sol$pw
m.risk = sqrt(as.vector(port.sol$pw %*% rcov %*% port.sol$pw))
ratio = -m.return/m.risk
assign("w",port.sol$pw,inherits=T)
}
return(ratio)
}
ef = effFrontier(averet=averet, rcov=rcov, shorts=shorts, wmax=wmax, nports = 100)
n = ncol(averet)
reshigh = rep(wmax,n)
if( shorts ) {
reslow = -reshigh
} else {
reslow = rep(0,n)
}
max.sh = which.max(ef$ret/ef$vol)
w = rep(0,ncol(averet))
xmin = optimize(f=optim.callback, interval=c(ef$ret[max.sh-1], upper=ef$ret[max.sh+1]),
averet=averet,rcov=rcov,reshigh=reshigh,reslow=reslow,shorts=shorts)
return(w)
}
ウェイトを計算します。
####シャープレシオを最大化するウェイトを計算#### weights <- maxSharpe(return_table,cov_table,shorts=F) names(weights) <- colnames(cov_table) round(weights * 100, 0)
以上です。グラフを描画すると結構感動するんですが、実際の投資の責任は取れませんので各自でご検討ください。
コメント