[譯]投資組合優化與R實踐

概述

最近，我在研究投資組合優化的問題，主要針對的是股票持倉的組合優化，我們會在這個分析過程中發現一些有意思的現象，并一步一步優化、檢驗我們的風控模型。本文將有四個部分分別闡述具體步驟。

• 在 第一部分（原文） 中，我將解釋什么是杠鈴策略，并初步建立風控模型，比較持倉策略和風險收益的關系。

• 在 第二部分（原文） 中，我將解釋什么是無風險利率假定，討論多項式擬合的情形。

• 在 第三部分（原文） 中，我將解釋如何通過放松約束最優化求解過程以避免非凹的情形，并做了實例演示。

• 在 第四部分（原文） 中，我將對比大盤策略、等權策略以及之前的優化策略之間的優劣。

請注意,本文不應該被作為投資建議。本文數據是基于之前觀察到的收益來模擬的，和歷史上的數據并不太一致。這些技術可以幫助了解如何更好地分配一個投資組合。它不應該被用作唯一的投資決策，如果你正在尋找的建議應該找到一個合格的專業機構。

第一部分

數字特征計算

當看到三個政府 ETF 債券(TLT、IEF、SHY)調整后的股息回報率,我注意到中間到期債券(IEF)風險收益情況比長期債券(TLT)更好。我以表格形式顯示結果。在本文中,我們將重新分析和圖形化展示我們的結果：

首先，用如下函數來獲取ETF的回報序列

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

結果如下：

杠鈴策略

如果你經常看娛樂投資電視臺,你最終會聽到“杠鈴策略”這個術語。這是指一個極端的投資組合分配方案。所有的權重都是極端情況,你可以想象這是一個杠鈴。在政府債券的投資組合,這將意味著購買期限長或短而不是持有中間。那么什么樣的風險收益情況下你會采用這個策略?

首先,我們將風險定義為投資組合的方差。有各種各樣的理由不使用方差,但它是從最古老的50年代開始這種類型的分析都是全新的。我們將定義收益為預期收益。在上面的表中,年回報率表示持有資產的預期收益為1年,標準差的平方表示風險。

假設投資組合只包括持有長期和短期債券,我們便需要計算投資組合的預期收益和風險。收益的計算是很容易的,這是兩種持倉的加權平均收益,權重就是每個資產的投入資本百分比。

RP = WTLT* RTLT + WSHY * RSHY
Where: WTLT + WSHY= 1

顯然這兩種資產具有相關性(在馬科維茨于1952年的博士論文發表之前,投資經理不了解相關性并且默認假設為1 -馬科維茨因此獲得了諾貝爾獎)。假設回報是正態分布的，那么投資組合方差將是：

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

風控模型

我們可以根據這兩個部分的知識改變持倉權重并為我們的杠鈴策略建立風險收益模型。

# 檢查相關性
corr = cor(merged[,vars],use=‘complete.obs’)
c = corr[‘tlt.Return’,‘shy.Return’]
# 假設一個杠鈴策略是持有長期和短期資產
# 定義風險、收益
ws = NULL
wt = NULL
mu = NULL
sigma = NULL
# 50個觀察
n=50
# 遍歷杠鈴策略的權重
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)
}
# 風險收益的數據集
rrProfile = data.frame(ws=ws,wt=wt,mu=mu,sigma=sigma)
# 繪圖
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”)

注意,上面的方程是二次的。我們可以配合我們剛剛創建的點畫出拋物線。注意,而習慣上把風險放在X軸上,而把擬合方差(風險)作為因變量放在Y軸。

# 為模型擬合一個二次函數
fit = lm(rrProfile$sigma ~ rrProfile$mu + I(rrProfile$mu^2))

接下來,我們圖上添加擬合線。

# 得到回歸系數
coe = fit$coefficients
# 得到每個回歸預測的風險值
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)
}
# 畫出預測邊值
lines(sfit,muf,col=‘red’)

tseries 包中的 portfolio.optim 也許是一個更好的選擇。我們只需要輸入預期收益率,它會直接返回出來最優組合權重。我們可以在最低預期回報率(比如 100% 持有 SHY)到最高預期回報率(比如 100% 持有 TLT)之間修改輸入的收益。注意, portfolio.optim 使用日回報率做計算,因此代碼將不得不做一些處理。我們假設一年255個交易日。

# 現在讓我們添加第三個標的。
# 除非我們想做一個格點搜索,否則我們需要對每個級別的回報減少風險來優化投資組合。
# portfolio.optim 在時間序列中不能有 NA 值。

m2 = removeNA(merged[,vars])
wSHY = NULL
wIEF = NULL
wTLT = NULL
er = NULL
eStd = NULL

# 在回報水平之間不斷循環搜索找到最優的投資組合，包括最小值(rSHY)和最大值(rTLT)
# portfolio.optim 使用日回報數據，因此我們不得不做出相應的調整

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])
}

# 畫出三個標的的有效邊界。

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)

結論

如下圖:

總資產組合中有效邊界的藍線表示其優于杠鈴策略。對于每個風險水平,預期回報都是更高的。從圖表上看,這表明添加 IEF 到組合將優化組合。進一步,我們看到杠鈴策略回報的逼近最大值,用三個標的組合的組合策略比之前的風險少了一半。

相關代碼

第二部分

在前面的文章中,我們構建了一個投資組合的有效邊界的債券，下一步，我們要找到超級有效的(或市場)的投資組合。如果您有不熟悉的概念,第二部分可以在維基百科上參考一些資料。

無風險利率假定

如果你不愿意看維基百科,我也會解釋相關概念的。如果你有一個保底回報率(無風險利率),那么資產位于圖表的y軸。在邊界的切點處畫一條切線，切點代表著非常有效的投資組合。你可以混合持有一定權重的組合標的和無風險資產,實現比邊界曲線更好的風險回報比。

明白了嗎?非常棒！

所以我們需要找到線和切點。首先,讓我們假定一個無風險利率。有些人會使用3個月的國債收益率。為了和數據匹配，我們需要將它處理成一年期的。我的銀行給我一個2%的年保底收益率，所以我將用2%。

多項式擬合

我們如何找到切點?當我們有兩個標的時,我們知道我們有一個二階多項式。當我們有三個標的時有一些存在缺陷的面（非凸時求極值較困難），在這種情況下我們停止投資 SHY，轉向投資 TLT。我們可以擬合高階多項式,但我們不能確保我們有一個凹面。或者我可以說,我們不能保證我們的切點總是高于邊值。同樣地，我們也可以想象一下二次的情形或許有切點存在負值。

作為一個例子,這里雖然六階多項式的擬合符合缺陷,但我們的切線點不是有用的。

只有一個實根,其余的都是虛根，我們需要另一種方法。

我們可以為第一部分里的邊值擬合一個多項式;此時在持倉組合中只有 SHY 和 IEF。雖然這樣也行得通，但是這不太通用。我想找到一個可以不管是什么邊值形狀都適用的通用解決方案。下個部分，我們會繼續討論這個問題。

第三部分

上一節，我們討論了用擬合曲線尋找有效邊值來建立投資組合所存在的問題。由于邊值存在的缺陷，我們不能保證你和曲線在投資組合的解空間內是凹的。我們需要其他方法來解決這個問題。

放松約束

本文所用方法是在無風險利率和每個邊值之間都畫一條線來計算這條線和邊值的差值是多少。資本市場線應該是不超過所有邊值的。

CMLi <= EFi

我們所找到的這個邊值的尺度意味著我們也許不能找到準確地市場投資組合。為避開這一點，我放松了上述約束：

Portfolioj that max Count(CMLi < EFi)

我整理以下R函數。注意,我已經轉向使用標準差作為風險度量尺度，這是更傳統的選擇。

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)
}

例子

讓我們使用埃克森美孚(XOM),IBM(IBM),中期政府債券ETF(IEF)這個組合做測試。這里假定你有importSeries()函數的定義。

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

日志的輸出是:

創建的圖:

結論

這個投資組合優化的給了我們一個發現更低邊界例子。這是一個不正常的現象。

這就是為什么我們結果子集只包括部分邊界的頂點(min(StdDev))和最大的回報。因為我們發現邊界最小收益到最大的收益,我們保證序列是有序的,所以只考慮了上部邊界。

一個更精確的方法是找到的區域包含市場組合的邊值然后用網格搜索尋找最優投資組合。上節我們討論了在一個范圍中擬合曲線的方法。如果有需求,我也可以用上面的方法再做一次。出于演示目的,我想我們應該足夠了。

第四部分

這節將對投資組合優化系列做一個總結,我們將基于組合優化和測試結果對CAPM市場投資組合構建一個交易策略。

值得重申的是:我所說不應該被當做投資建議。這些結果是基于之前觀察到的收益并且是一些模擬值。這些技術可以幫助了解如何更好地分配一個投資組合。它不應該被當作是唯一的投資決策。如果你正在尋找的建議，還是找一個合格的專家比較好。

在馬科維茨的工作的基礎上,特雷諾,夏普等人開發了資本資產定價模型(CAPM)。在1990年，他們因為這項工作與馬科維茨共同獲得了諾貝爾獎。CAPM是一個一般均衡模型。模型假定市場價格反映了所有可獲得的信息并且反映一個標的的“公平”價值。在此基礎上,市場投資組合可以證明是市值加權組合。市值(或市值)被定義為股價乘以流通股的數量，也就是公司的股本總額。公司的權重是該公司市值除以所有證券的總市值。

資本加權指標和指數基金已成為標準。標準普爾是大多數人考慮的標準“市場投資組合”。我們將參考一個市值加權策略對我們的投資組合優化策略進行測試。

現在的CAPM還存在諸多漏洞，有很多方法都能發現這些問題。一種方法是說現在的價格不是公允價值,而是將均值當作公允價值。在這種情況下,當價格高于公允價值,市值加權組合將對過定價過高的證券資產給予過大的權重。當它用均值取代后, 投資組合的表現將由于權重超標而變差。

這個理論是著名的 羅伯特?阿諾特 提出的。我強烈推薦這本書,《基本面指數:一種更好的投資方式》。他認為,任何打破用價格打破相關性的投資組合策略隨著時間的推移都將跑贏資本化指數。他在書中提到他創造了一個新的指數,他簡單地假定每個標的都是等權重的(標準普爾發布了這個指數)。正因為如此,我們還將在標的同等權重條件下測試我們的策略。

組合優化策略

這是我們的投資組合優化策略:

1. 每個季度初,用上一季度收益計算市場投資組合。

2. 對當前季度使用當前組合。

3. 下個季度的開始,循環回到第一步

4. 在我們的投資組合中至少需要3個股票。

5. 沒有做空。

6. 用2%作為無風險利率。

7. 每次分析的第一個季度如果優化失敗就使用同等權重的投資組合。

當價格走勢是按季度選取的這種策略往往會跑贏大盤。比如如果上個季度的收益和波動性可以準確預測本季度的值的情況就是這樣。此外，我們也不考慮交易成本。而且，2%的無風險利率是靜態的,嚴格的說，我們應該在每個季度開始時使用3個月國債的利率。這就是為什么這只是一個例子，我們假定了很多美好的假設。

首先,我們需要修改我們以前創建的marketPortfolio()函數。你可以在這里找到它。新函數：

marketPortfolio = function(merged, rf, returnNames, weightNames,graph=FALSE, points=500, maxWeight=.334, Debug=FALSE)

改進之處

我們的改進之處:

1. 如果我們需要可以添加一個調試選項打印輸出

2. 增加容錯功能。有時直接得到一個可行解是不可能的,這個函數需要相應的檢測并且處理錯誤。

3. 資本的最大回報我們限定在50%以下，這個值太大會導致其他奇怪的行為。

4. 同樣,把最小收益的下界定在.005%。

5. 如果最大收益是< 0,那么簡單地找到最小方差投資組合。

6. 添加一個maxWeight選項,讓我們限制每個證券標的的權重。

我們考慮的股票池現在是28個道瓊斯成分股（因為某些原因雅虎金融值提供了28只而不是30只）。我預先下載并存儲這些收益值為一個 .rda 文件。 你可以在這里得到它 。我計算每個股票的初始市值權重并存儲為 .csv 文件。 你可以在這里得到它 。

# 讀取已存的收益值
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()

我們的年回報率表:

結論

我們的投資組合優化策略優于大盤權重策略,但跑輸了等權重策略。如果你支持阿諾特的話就覺得這沒什么奇怪的了,這只是因為我們沒有打破價格的相關性罷了。

這是相關性的圖表:

我們已經創建了一個和大盤權重策略非常相關的策略，但是還是不如等權策略。等權策略和大盤權重策略的關聯度是非常有趣的。

這是收益繪制的時間序列:

有趣的是,可以看到在圖中綠色的部分顯示我們的投資組合在2009年3月份的市場底部開始有一個快速反彈。這大大跑贏了大盤權重組合。