第6回R勉強会@東京(Tokyo.R#06)で発表してきた
「Rで学ぶ金融工学入門〜現代ポートフォリオ理論偏〜」というタイトルで発表してきた。
以下、その資料
Tokyo.R(#06)
View more presentations from teramonagi.
以下、ソースコード
quadprog,RFinanceYJのインストールが必要。
library(quadprog) library(RFinanceYJ) #2005/4月〜2010/4月までの4つの投資信託のデータを月次で取得. #野村国内債券インデックスF(確定拠出年金)'01315027' #野村国内株式インデックスF(確定拠出年金)'0131202A' #野村外国債券インデックスF(確定拠出年金)'0131102B' #野村外国株式インデックスF(確定拠出年金)'01312022' codes <- c('01315027','0131202A','0131102B','01312022') assets.name <- c("domestic.bond","domestic.stock","foreign.bond","foreign.stock") NUM.ASSET <- length(codes) #基準価額取得 prices <- sapply(codes,function(x)quoteFundTsData(x,time.interval='monthly' ,since='2005-04-01',date.end='2010-04-30')$constant.value) head(prices) #リターンに変換(x月のリターン=x月末 x-1月末株価 - 1.) returns <- apply(prices,2,function(x)x[-1]/x[1:length(x)-1]-1) colnames(returns) <- assets.name head(returns) #相関係数のチェックと散布図の確認 cor(returns) pairs(returns) ############################################################################### #リターンの共分散行列 #月率のものを年率に換算するために12倍(分布がiidであることを仮定) returns.covmat.estimated <- 12*cov(returns) #それぞれの資産のリスク diag(returns.covmat.estimated)^0.5 risk.min <- min(diag(returns.covmat.estimated)^0.5) risk.max <- max(diag(returns.covmat.estimated)^0.5) #Risk:ポートフォリオウェイトとリターンの共分散行列からリスク(年率)を計算 #weight : ポートフォリオウェイトベクトル #returns.covariancematrix : リターンの共分散行列 Risk <- function(weight,returns.covmat){ weight <- matrix(weight,nrow=length(weight)) return(as.numeric(t(weight)%*%returns.covmat%*%weight)^0.5) } #分散効果のチェック weights <- c(0.0,0.7,0.0,0.3) Risk(weights,returns.covmat.estimated) ############################################################################### #リターンを設定(リスクプレミアム・ヒストリカル・GPIF値) returns.estimated <- 0.3*diag(returns.covmat.estimated)^0.5 returns.estimated <- 12*colMeans(returns) returns.estimated <- c(0.03,0.048,0.035,0.05) return.min <- min(returns.estimated) return.max <- max(returns.estimated) #http://www.gpif.go.jp/kanri/kanri03.html #Returnk:ポートフォリオウェイトとリターンから #ポートフォリオのリターンを計算 #weight : ポートフォリオウェイトベクトル #returns : リターンベクトル Return <- function(weights,returns){ return(sum(weights*returns)) } #外国債券・外国株式を50%ずつ保有↓場合のリターン Return(c(0,0,0.5,0.5),returns.estimated) ############################################################################### #合計が1になるようなnum個の数値の組み合わせを返す関数 combination.constantsum <- function(num,val.step,val.current=1) { if(num==1){ return(val.current) } result <- NULL for(x in seq(0,val.current,val.step)){ result <- rbind(result,cbind(x,combination.constantsum(num-1,val.step,val.current-x))) } return(result) } weights.combination <- combination.constantsum(NUM.ASSET,0.05) head(weights.combination) #投資可能領域のプロット lim.x <-c(risk.min*0.9 ,risk.max*1.1) lim.y <-c(return.min*0.9,return.max*1.1) plot(t(apply(weights.combination,1,function(x)c(Risk(x,returns.covmat.estimated) ,Return(x,returns.estimated)))),xlab="Risk",ylab="Return",xlim=lim.x,ylim=lim.y) par(new=TRUE) plot(diag(returns.covmat.estimated)^0.5,returns.estimated,xlab="Risk",ylab="Return" ,col=rainbow(NUM.ASSET),pch=19,cex=2,xlim=lim.x,ylim=lim.y) legend(lim.x[1],lim.y[2],assets.name,col=rainbow(NUM.ASSET),pch=19) ############################################################################### #最適化計算 NUM.OPTIMIZED <- 40 portfolios <- list() returns.target <- seq(return.min,return.max,length.out=NUM.OPTIMIZED) for(return.target in returns.target){ #制約条件行列(ウエイトの合計が1,期待リターン制約、空売りなし) A <- cbind(rep(1,NUM.ASSET),returns.estimated,diag(NUM.ASSET)) b <- c(1,return.target,rep(0,NUM.ASSET)) weight.opt <- solve.QP(returns.covmat.estimated,rep(0,NUM.ASSET),A,b,meq=2)$solution portfolios[[length(portfolios)+1]] <- list(weight=weight.opt,risk=Risk(weight.opt,returns.covmat.estimated) ,return=Return(weight.opt,returns.estimated)) } #最適化計算された効率的フロンティア上のポートフォリオをプロット plot(t(apply(weights.combination,1,function(x)c(Risk(x,returns.covmat.estimated) ,Return(x,returns.estimated)))),xlab="Risk",ylab="Return",xlim=lim.x,ylim=lim.y) par(new=TRUE) plot(t(sapply(portfolios,function(x)c(x$risk,x$return))) ,xlab="Risk",ylab="Return",col='red',pch=19,cex=2,xlim=lim.x,ylim=lim.y) legend(lim.x[1],lim.y[2],'effective frontier',col='red',pch=19) #効率的フロンティア上のウェイトをプロット portfolios.weight <- sapply(portfolios,function(x)x$weight) barplot(portfolios.weight,col=rainbow(NUM.ASSET),names=round(returns.target*100,1)) title("portfolio weight on effective frontier","estimated return") legend(NUM.OPTIMIZED*0.8,0.95,assets.name,col=rainbow(NUM.ASSET),pch=19,bg='white')