译文 | 量化投资教程：投资组合优化与R实践

require(fImport)
require(PerformanceAnalytics)
# 将股票数据加载到一个时间序列对象的函数
importSeries = function (symbol,from,to) {
# 从雅虎读取金融数据
input = yahooSeries(symbol,from=from,to=to)
# 列名
adjClose = paste(symbol,".Adj.Close",sep="")
inputReturn = paste(symbol,".Return",sep="")
CReturn = paste(symbol,".CReturn",sep="")
# 计算收益率并生成时间序列
input.Return = returns(input[,adjClose])
colnames(input.Return)[1] = inputReturn
input = merge(input,input.Return)
# 计算累积收益率并生成时间序列
input.first = input[,adjClose][1]
input.CReturn = fapply(input[,adjClose],FUN=function(x) log(x) - log(input.first))
colnames(input.CReturn)[1] = CReturn
input = merge(input,input.CReturn)
# 删掉一些没用的东西，如果你不知道就不用删除
rm(input.first,input.Return,input.CReturn,adjClose,inputReturn,CReturn)
# 返回时间序列
return(input)
}

# 获取短中期和长期政府债券的收益率序列

from = “2001-01-01″
to = “2011-12-16″
tlt = importSeries(“tlt”,from,to)
shy = importSeries(“shy”,from,to)
ief = importSeries(“ief”,from,to)
merged = merge(tlt,shy)
merged = merge(merged,ief)
vars = c(“tlt.Return”,“shy.Return”,“ief.Return”)
# 计算年回报率
t = table.AnnualizedReturns(merged[,vars],Rf=mean(merged[,“shy.Return”],na.rm=TRUE))
t

RP = WTLT* RTLT + WSHY * RSHY

Where: WTLT + WSHY= 1

Vp = WTLT 2*σ2TLT+ WSHY 2* σ2SHY + WTLT* WSHY * σTLT * σSHY *CorrTLT,SHY

Where:  WTLT+ WSHY = 1

require(fImport)
require(PerformanceAnalytics)
require(tseries)
require(stats)
options(scipen=100)
options(digits=4)
#Function to load stock data into a Time Series object
importSeries = function (symbol,from,to) {
#Read data from Yahoo! Finance
input = yahooSeries(symbol,from=from,to=to)
#Character Strings for Column Names
adjClose = paste(symbol,".Adj.Close",sep="")
inputReturn = paste(symbol,".Return",sep="")
CReturn = paste(symbol,".CReturn",sep="")
#Calculate the Returns and put it on the time series
input.Return = returns(input[,adjClose])
colnames(input.Return)[1] = inputReturn
input = merge(input,input.Return)
#Calculate the cumulative return and put it on the time series
input.first = input[,adjClose][1]
input.CReturn = fapply(input[,adjClose],FUN=function(x) log(x) - log(input.first))
colnames(input.CReturn)[1] = CReturn
input = merge(input,input.CReturn)
#Deleting things (not sure I need to do this, but I can't not delete things if
# given a way to...
rm(input.first,input.Return,input.CReturn,adjClose,inputReturn,CReturn)
#Return the timeseries
return(input)
}
#Get Return Series for Short Medium and Long Term Gov Bonds
from = "2001-01-01"
to = "2011-12-16"
tlt = importSeries("tlt",from,to)
shy = importSeries("shy",from,to)
ief = importSeries("ief",from,to)
merged = merge(tlt,shy)
merged = merge(merged,ief)
vars = c("tlt.Return","shy.Return","ief.Return")
#Calculate Annualized Returns
t = table.AnnualizedReturns(merged[,vars],Rf=mean(merged[,"shy.Return"],na.rm=TRUE))
t
#Get the annualized return and StdDev for each series from the table
rTLT = t['Annualized Return',"tlt.Return"]
rSHY = t['Annualized Return',"shy.Return"]
rIEF = t['Annualized Return','ief.Return']
sTLT = t['Annualized Std Dev',"tlt.Return"]
sSHY = t['Annualized Std Dev',"shy.Return"]
sIEF = t['Annualized Std Dev',"ief.Return"]
#Check the correlations
corr = cor(merged[,vars],use='complete.obs')
c = corr['tlt.Return','shy.Return']
#Assuming a barbell strategy of long and short holdings, what is the risk return profile
ws = NULL
wt = NULL
mu = NULL
sigma = NULL
#50 observations
n=50
#Loop through the weights of the barbell
for (i in 0:n){
wsi = i/n;
wti = 1-wsi;
mui = wsi * rSHY + wti * rTLT
sigmai = wsi*wsi*sSHY*sSHY + wti*wti*sTLT*sTLT + wsi*wti*sSHY*sTLT*c
ws = c(ws,wsi)
wt = c(wt,wti)
mu = c(mu,mui)
sigma = c(sigma,sigmai)
}
#Risk Return Profile Data Frame
rrProfile = data.frame(ws=ws,wt=wt,mu=mu,sigma=sigma)
#Plot the profile
plot(rrProfile$sigma,
 rrProfile$mu,
xlim=c(0,.022),
ylim=c(0,.08),
ylab="Expected Yearly Return",
xlab="Expected Yearly Variance",
main="Efficient Frontier for Government Bond Portfolios")
#Fit a quadratic function to the profile
fit = lm(rrProfile$sigma ~ rrProfile$mu + I(rrProfile$mu^2))
#How is the fit?
summary(fit)
#get the coefficients
coe = fit$coefficients
#Get predicted values risk values for each return
muf = NULL
sfit = NULL
for (i in seq(0,.08,by=.001)){
muf = c(muf,i)
s = coe[1] + coe[2]*i + coe[3]*i^2
sfit = c(sfit,s)
}
#plot the predicted frontier
lines(sfit,muf,col='red')
#now let's add the 3rd asset. Unless we want to do a grid search, we need
#to optimize the portfolio, minimizing the risk for each level of return
#portfolio.optim cannot have NA values in the time series, filter them out
m2 = removeNA(merged[,vars])
wSHY = NULL
wIEF = NULL
wTLT = NULL
er = NULL
eStd = NULL
#loop through finding the optimum portfolio for return levels between
#the minimum (rSHY) and the max(rTLT)
#
#portfolio.optim uses daily returns, so we have to adjust accordingly
for (i in seq((rSHY+.001),(rTLT-.001),length.out=100)){
pm = 1+i
pm = log(pm)/255
opt = portfolio.optim(m2,pm=pm)
er = c(er,exp(pm*255)-1)
eStd = c(eStd,opt$ps*sqrt(255))
wTLT = c(wTLT,opt$pw[1])
wSHY = c(wSHY,opt$pw[2])
wIEF = c(wIEF,opt$pw[3])
}
#Plot the efficient frontier for the 3 assets.
lines(eStd^2,er,col='blue')
legend(.014,0.015,c("'Barbell' Strategy","All Assets"),
col=c("red","blue"),
lty=c(1,1))
solution = data.frame(wTLT,wSHY,wIEF,er,eStd)

marketPortfolio = function(merged,rf,returnNames, weightNames,graph=FALSE){
      # 为投资组合创建空数据框以初始化
      weights = data.frame(t(rep(NA,length(weightNames))))
      colnames(weights) = weightNames
      weights = weights[–1,]
      # 计算年化收益
      t = table.AnnualizedReturns(merged[,returnNames])
      # 优化范围
      maxRet = max(t[‘Annualized Return’,]) – .005
      minRet = min(t[‘Annualized Return’,]) + .005
      #portfolio.optim 没有 NA 值，进行过滤
      m2 = removeNA(merged[,returnNames])
      er = NULL
      eStd = NULL
      # 在每个收益水平上循环搜索最优组合
      # portfolio.optim 是日收益，做出相应调整
      for (i inseq(minRet,maxRet,length.out=500)){
            pm = 1+i
            pm = log(pm)/255
            opt = portfolio.optim(m2,pm=pm)
            er = c(er,exp(pm*255)–1)
            eStd = c(eStd,opt$ps*sqrt(255))
            w = t(opt$pw)
            colnames(w) = weightNames
            weights = rbind(weights,w)
      }
      solution = weights
      solution$er = er
      solution$eStd = eStd
      #找到最小 Std 和最大 Er 的下标
      minIdx = which(solution$eStd == min(solution$eStd))
      maxIdx = which(solution$er == max(solution$er))
      # 获取结果子集
      subset = solution[minIdx:maxIdx,c(“er”,“eStd”)]
      subset$nAbove = NA
      #对于子集中的每一个值, 计算点的总数，在下面的点和RF资产之间画线
      for (i in seq(1,maxIdx–minIdx+1)){
            toFit = data.frame(er=rf,eStd=0)
            toFit = rbind(toFit,subset[i,c(“er”,“eStd”)])
            fit = lm(toFit$er ~ toFit$eStd)
            poly = polynomial(coef = fit$coefficients)
            toPred = subset
            colnames(toPred) = c(“actEr”,“eStd”)
            toPred$er = predict(poly,toPred[,“eStd”])
            toPred$diff = toPred$er – toPred$actEr
            subset[i,“nAbove”] = nrow(toPred[which(toPred$diff > 0),])
      }
      # 得到切点
      # 线以下是最大化
      max = max(subset$nAbove)
      er = subset[which(subset$nAbove == max),“er”]
      eStd = subset[which(subset$nAbove == max),“eStd”]
      # 市场投资组合的下标
      idx = which(solution$er == er & solution$eStd == eStd)
      # 画线
      if (graph){
            maxStd = max(solution$eStd) + .02
            maxRetg = max(solution$er) + .02
            plot(solution$eStd,
                        solution$er,
                        xlim=c(0,maxStd),
                        ylim=c(0,maxRetg),
                        ylab=“Expected Yearly Return”,
                        xlab=“Expected Yearly Std Dev”,
                        main=“Efficient Frontier”,
                        col=“red”,
                        type=“l”,
                        lwd=2)
            abline(v=c(0), col=“black”, lty=“dotted”)
            abline(h=c(0), col =“black”, lty=“dotted”)
            toFit = data.frame(er=rf,eStd=0)
            toFit = rbind(toFit,solution[idx,c(“er”,“eStd”)])
            fit = lm(toFit$er ~ toFit$eStd)
            abline(coef=fit$coefficients,col=“blue”,lwd=2)
      }
      # 返回投资组合权重、eStd 和 eR
      out = solution[idx,]
      return (out)
}

require(polynom)
require(fImport)
require(PerformanceAnalytics)
require(tseries)
require(stats)
from = “2003-01-01″
to = “2011-12-16″
xom = importSeries(“xom”,from,to)
ibm = importSeries(“ibm”,from,to)
ief = importSeries(“ief”,from,to)
merged = merge(xom,ibm)
merged = merge(merged,ief)
vars = c(“xom.Return”,“ibm.Return”,“ief.Return”)
vars2 = c(“xom”,“ibm”,“ief”)
mp = 
marketPortfolio(merged,.02,vars,vars2,graph=TRUE)
mp

marketPortfolio <- function(merged,
                            rf,
                            returnNames,
                            weightNames,
                            graph=FALSE,
                            points=500,
                            maxWeight=.334,
                            Debug=FALSE){
# 初始化组合权重数据帧
weights = data.frame(t(rep(NA,length(weightNames))))
colnames(weights) = weightNames
weights = weights[-1,]
# 年化收益
t = table.AnnualizedReturns(merged[,returnNames])
# 优化范围
maxRet = max(t['Annualized Return',]) - .005
minRet = min(t['Annualized Return',]) + .005
# 设置回报上下限
# 在.005%到50%之间
maxRet = min(.5,maxRet)
minRet = max(0.005,minRet)
# 如果都是预期回报均为负数，
# 那么只返回收益最大的投资组合
if (maxRet < 0){
  minRet = maxRet
  points = 1
}
# 调试打印
if (Debug){
  print("Max Return")
  print(maxRet)
  print("Min Return")
  print(minRet)
}
# portfolio.optim 不能处理含有 NA 值的时间序列，这里过滤掉 NA 值
m2 = removeNA(merged[,returnNames])
er = NULL
eStd = NULL

# 在回报水平之间循环，寻找最优的投资组合

# portfolio.optim 使用日收益,因此我们必须作出相应调整

ok = FALSE
for (i in seq(minRet,maxRet,length.out=points)){
pm = 1+i
pm = log(pm)/255
# 调试
if (Debug){
   print("Finding Optimum for")
   print("ER")
   print(pm)
   print(exp(pm*255)-1)
}
# 最优化 各组权重 <= 最大权重
# tryCatch 捕获异常值
opt = tryCatch(portfolio.optim(m2,
                               pm = pm,
                               reshigh = rep(maxWeight,length(weightNames)),
                               shorts = FALSE),
                               error = function(err) return(NULL))
# 打印最优结果
if (Debug) print(opt)
# 检查是否存在可行解
if (!is.null(opt)){
    er = c(er, exp(pm*255) - 1)
    eStd = c(eStd, opt$ps*sqrt(255))
    w = t(opt$pw)
    colnames(w) = weightNames
    weights = rbind(weights,w)
    # 更新变量
    ok = (ok | TRUE)
} else {
    print("ERROR IN THIS TRY")
    # 更新变量
    ok = (ok | FALSE)
  }
}
# 如果在边界没有可行解，则返回NULL
if (!ok){
  return (NULL)
}
solution = weights
solution$er = e
solution$eStd = eStd
# 找到最小标准差和最大平均收益点的索引值
minIdx = which(solution$eStd == min(solution$eStd))
maxIdx = which(solution$er == max(solution$er))
if (Debug){
  print(minIdx)
  print(maxIdx)
}
# 切分结果
subset = solution[minIdx:maxIdx,c("er","eStd")]
subset$nAbove = NA

# 对每个分片的值,计算在划线下面的点的数量

for (i in seq(1,maxIdx-minIdx+1)){
    toFit = data.frame(er=rf,eStd=0)
    toFit = rbind(toFit,subset[i,c("er","eStd")])
    fit = lm(toFit$er ~ toFit$eStd)
    poly = polynomial(coef = fit$coefficients)
    toPred = subset
    colnames(toPred) = c("actEr","eStd")
    toPred$er = predict(poly,toPred[,"eStd"])
    toPred$diff = toPred$er - toPred$actE
    subset[i,"nAbove"] = nrow(toPred[which(toPred$diff > 0),])
}

# 得到基准线以下的数量最大化的切点
max = max(subset$nAbove)
er = subset[which(subset$nAbove == max),"er"]
eStd = subset[which(subset$nAbove == max),"eStd"]
# 如果找到多个组合,返回第一个
if (length(er) > 1){
   er = er[1]
   eStd = eStd[1]
}
# 市场投资组合的索引
idx = which(solution$er == er & solution$eStd == eStd)
if (Debug){
    print("solution")
    print(er)
    print(eStd)
    print(solution[idx,])
}
# 如果需要可以加一条基准线
if (graph){
    maxStd = max(solution$eStd) + .02
    maxRetg = max(solution$er) + .02
    plot(solution$eStd,
    solution$er,
    xlim=c(0,maxStd),
    ylim=c(0,maxRetg),
    ylab="Expected Yearly Return",
    xlab="Expected Yearly Std Dev",
    main="Efficient Frontier",
    col="red",
    type="l",
    lwd=2)
    abline(v=c(0), col="black", lty="dotted")
    abline(h=c(0), col ="black", lty="dotted")
    toFit = data.frame(er=rf,eStd=0)
    toFit = rbind(toFit,solution[idx,c("er","eStd")])
    fit = lm(toFit$er ~ toFit$eStd)
    abline(coef=fit$coefficients,col="blue",lwd=2)
}
# 返回市场投资组合权重、eStd 、eR
out = solution[idx,]
return (out)
}

# 读取已存的收益值load(“D:\Workspace\PortfolioOptimization\returns.rda”)
returns = results
rm(results)
stocks = colnames(returns)# 获取大盘权重stockWgt = read.table(“D:\Workspace\PortfolioOptimization\stocks.csv”,
                          header=TRUE,sep=“,”)[,“Weight”]# 计算大盘权重投资组合的回报results = as.matrix(returns) %*% stockWgt
colnames(results) = c(“CapWeight Portfolio”)
results = as.timeSeries(results)# 计算等权重投资组合的回报ret = t(as.matrix(rep(1/length(stocks),length(stocks))))
resEqual = as.matrix(returns) %*% t(ret)
colnames(resEqual) = “EqualWeight Portfolio”# 汇总结果到一个时间序列对象中results = cbind(results,resEqual)# 获取收益值的日期序列dates = time(returns)
dates = dates@Data# 从日期计算季度qtrs = quarters(dates)
qtrs = cbind(dates,qtrs)
keep = NULL
lastQtr = “n”#遍历日期和季度序列,只保留每个季度第一天的日期for (i in seq(1,nrow(qtrs))){      if (qtrs[i,2] == lastQtr){            if (i == nrow(qtrs)){
                  keep = c(keep,1)
            } else {
                  keep = c(keep,0)
            }
      }else {
            keep = c(keep,1)
      }
 
      lastQtr = qtrs[i,2]
}
qtrs = cbind(qtrs,keep)# 获取每个季度第一天的下标indx = which(qtrs[,3] == 1)# 对每个周期的第一个季度,使用等权策略res = as.matrix(returns[indx[1]:(indx[2]–1),]) %*% t(ret)#对每个周期基于上一季度的数据进行循环计算大盘组合的表现for (i in seq(2,length(indx)–1)){      print(“Running “)      print(i)     
      # 得到上季度股票回报的子集 
      subset = returns[indx[i–1]:(indx[i]–1),]      s = start(subset)
      e = end(subset)      print(“Fitting for:”)      print(s)      print(e)     
      # 计算大盘投资组合
      mp = marketPortfolio(subset,.02,stocks,stocks,
                           graph=TRUE,
                           points=500,
                           Debug=FALSE)     
      #如果优化失败,使用等权策略
      if (is.null(mp)){
            ret = t(as.matrix(rep(1/length(stocks),length(stocks))))
      } else {
            ret = as.matrix(mp[,stocks])
      }     
      # 本季度的子集
      subRes = returns[indx[i]:(indx[i+1]–1),]      s = start(subRes)
      e = end(subRes)      print(“Calculating Returns for:”)      print(s)      print(e)     
      # 计算当前季度的大盘策略的收益并追加在收益序列后面
      subRes = as.matrix(subRes) %*% t(ret)
      res = rbind(res,subRes)
}# 循环计算时，序列的最后一天不计算收益subRes = returns[nrow(returns),]
subRes = as.matrix(subRes) %*% t(ret)
res = rbind(res,subRes)# 添加组合优化策略colnames(res) = “Portfolio Optimization”
res = as.timeSeries(res)
results = cbind(results,res)#计算年化收益统计特征table.AnnualizedReturns(results)#计算并绘制相关性png(“d:\Workspace\PortfolioOptimization\mpCorr.png”)
chart.Correlation(results,histogram=TRUE,pch=“+”)
dev.off();##计算并绘制累积收益png(“d:\Workspace\PortfolioOptimization\mpReturns.png”)
chart.CumReturns(results,
            main=“Total Returns CapWeight vs PortOpt”,
            legend.loc=“topleft”)
dev.off()