第6回R勉強会@東京(Tokyo.R#06)で発表してきた

「Rで学ぶ金融工学入門〜現代ポートフォリオ理論偏〜」というタイトルで発表してきた。
以下、その資料

以下、ソースコード
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')