GLPK/投資組合最佳化

本示例取自投資理論，屬於金融經濟學的分支。其目標是在已知歷史表現的一組潛在投資中分配資本，以在收益和投資者的風險承受能力之間取得最佳平衡。

概要

本教程示例展示瞭如何使用 GLPK/GMPL 來最佳化投資組合，其中投資風險使用平均絕對偏差 (MAD) 標準來衡量。除了投資組合最佳化外，本示例還重點介紹了在其他 GMPL 應用中可能實用的技術，包括

Cholesky 分解正定矩陣。
從多元正態分佈中生成樣本，用於模擬歷史收益。
從 GMPL 中為 GLPK 的隨機數生成器設定種子。注意：新的randseed選專案前正在評估中。

相同作者的其他示例可在運籌學入門處找到。

平均絕對偏差 (MAD) 投資組合最佳化

投資組合最佳化指的是在金融資產集之間分配資本的過程，以實現收益和風險之間所需的平衡。該問題的經典馬科維茨方法使用投資組合收益的預期方差來衡量風險。該標準為最佳投資組合中資產的相對權重提供了一個二次規劃問題。

1991 年，Konno 和 Yamazaki ^[1] 提出了一個線性規劃模型，用於投資組合最佳化，其中風險由平均絕對偏差 (MAD)（相對於預期收益的）來衡量。使用 MAD 作為風險度量標準產生的投資組合具有馬科維茨投資組合不具備的幾個理想特性，包括二階隨機支配。

如 Konno 和 Yamazaki 最初所提出的那樣，首先需要一個收益歷史記錄 $R_{i}(t_{n})$ ，用於資產集中的每項資產 $S$ 。時間 $t_{n}$ 的收益由資產價格的變化決定， $R_{i}(t_{n})={(P_{i}(t_{n})-P_{i}(t_{n-1}))}/{P_{i}(t_{n-1})}$ 。對於每項資產，預期收益透過以下公式估算：

{\bar {R}}_{i}\approx {\frac {1}{N}}\sum _{n=1}^{N}R_{i}(t_{n})

.

投資者指定一個最低要求收益 $R_{p}$ 。投資組合最佳化問題是確定分配給每項資產的總投資的比例，即 $w_{i}$ ，其最小化與平均值的平均絕對偏差

\min _{w_{i}}{\frac {1}{N}}\sum _{n=1}^{N}|\sum _{i\in S}w_{i}(R_{i}(t_{n})-{\bar {R}}_{i})|

受所需回報和固定總投資約束

{\begin{array}{rcl}\sum _{i\in S}w_{i}{\bar {R}}_{i}&\geq &R_{p}\\\quad \\\sum _{i\in S}w_{i}&=&1\end{array}}

最小所需回報的值， $R_{p}$ 表達了投資者的風險承受能力。 $R_{p}$ 的值越小，則可行解空間越大，導致投資組合的 MAD 風險指標值越低。增加 $R_{p}$ 會導致投資組合具有更高的風險。風險和回報之間的關係是投資理論的基本原理。

此公式不會對權重 $w_{i}$ 設定個體邊界。特別是， $w_{i}<0$ 對應於相關資產的賣空。如果投資者對賣空、單一資產的最大投資額或投資組合多元化的其他要求設定限制，則可以新增約束條件。根據可用投資機會的集合，額外的約束條件可能會導致不可行的投資問題。

實現細節

本示例的目的是展示在 GMPL 中對投資組合最佳化問題的建模。主要任務是將平均絕對偏差目標轉化為線性形式，並模擬給定一組潛在投資的均值和協方差的的歷史回報。

MAD 目標的重新表述

目標函式的以下公式改編自 Gerald Curnuejols 和 Reha Tutuncu (2007) ^[2] 的 "金融最佳化方法"，該公式反過來遵循 Feinstein 和 Thapa (1993) ^[3]。透過引入決策變數 $y_{n}\geq 0$ 和 $z_{n}\geq 0$ ， $n=1,\ldots ,N$ ，簡化模型使目標變為

\min _{w_{i},y_{n},z_{n}}{\frac {1}{N}}\sum _{n=1}^{N}(y_{n}+z_{n})

受制於約束

{\begin{array}{rcl}\sum _{i\in S}w_{i}{\bar {R}}_{i}&\geq &R_{p}\\\\\sum _{i\in S}w_{i}&=&1\\\\y_{n}-z_{n}&=&\sum _{i\in S}w_{i}(R_{i}(t_{n})-{\bar {R}}_{i})\quad {\mbox{for}}\quad n=1,\ldots ,N\end{array}}

正如 Feinstein 和 Thapa 所討論的，這個版本將問題簡化為 $N+2$ 個約束，在 $2N+{\mbox{card}}(S)$ 個決策變數中，注意到 card 函式。

為 GLPK 偽隨機數生成器設定種子

不幸的是，GMPL 沒有提供一種方法來設定 GLPK 中的偽隨機數生成器的種子。相反，以下 GMPL 程式碼片段使用 gmtime() 函式來查詢自 1970 年以來的秒數。丟棄前導數字可以避免後續的溢位錯誤，而平方可以返回一個每秒變化很大的數字。提取最低位數字可以產生一個介於 0 到 100,000 之間的數字，該數字決定在後續使用之前對 GLPK 的偽隨機數生成器進行取樣的次數。這種技巧沒有保證其統計特性。

 /* A workaround for the lack of a way to seed the PRNG in GMPL */
param utc := prod {1..2} (gmtime()-1000000000);
param seed := utc - 100000*floor(utc/100000);
check sum{1..seed} Uniform01() > 0;

歷史回報的模擬

對於這個實現，歷史回報是透過假設了解資產回報的均值和協方差來模擬的。我們從一個平均回報向量開始 ${\bar {R}}$ 和協方差矩陣 $\Sigma$ ，由下式估計

\Sigma _{ij}\approx {\frac {1}{N-1}}\sum _{n=1}^{N}(R_{i}(t_{n})-{\bar {R}}_{i})(R_{j}(t_{n})-{\bar {R}}_{j})

歷史回報的模擬需要從多變數正態分佈中生成樣本。為此，我們計算 Cholesky 分解，其中，對於一個正半定 $\Sigma$ ， $\Sigma =CC^{T}$ ，並且 $C$ 是一個下三角矩陣。以下 GMPL 程式碼片段實現了 Cholesky-Crout 演算法。

/* Cholesky Lower Triangular Decomposition of the Covariance Matrix */
param C{i in S, j in S : i >= j} :=
    if i = j then
        sqrt(Sigma[i,i]-(sum {k in S : k < i} (C[i,k]*C[i,k])))
    else
        (Sigma[i,j]-sum{k in S : k < j} C[i,k]*C[j,k])/C[j,j];

如果沒有錯誤檢查，除非 $\Sigma$ 是正定的，否則該程式碼片段將失敗。對於真實世界的資料，協方差矩陣通常是正定的，所以這通常不是問題。但是，如果嘗試在投資資產集中包含無風險資產（如政府債券），這就會成為一個問題。

一旦計算出 Cholesky 因子 $C$ ，模擬收益向量 $R(t_{n})$ 由 $R(t_{n})={\bar {R}}+CZ(t_{n})$ 給出，其中 $Z(t_{n})$ 的元素是從均值為零、方差為一的正態分佈中獨立抽取的樣本。

/* Simulated returns */
param N default 5000;
set T := 1..N;
param R{i in S, t in T} := Rbar[i] + sum {j in S : j <= i} C[i,j]*Normal(0,1);

PortfolioMAD.mod

將此模型儲存為PortfolioMAD.mod.

/* Portfolio Optimization using Mean Absolute Deviation

Jeff Kantor 
December 4, 2009
Revised: December 6, 2009 to fix problem with random variate generation.
Revised: December 7, 2009 to add a 'seeding' of the PRNG
Revised: July 8, 2010 reformatted for GLPK Wikibook */

/* Stock Data */

set S;                                              # Set of stocks
param Rbar{S};                                      # Means of projected returns
param Sigma{S,S};                                   # Covariance of projected returns
param Rp  default (1/card(S))*sum{i in S} Rbar[i];  # Lower bound on portfolio return

/* Generate sample data */

/* Cholesky Lower Triangular Decomposition of the Covariance Matrix */
param C{i in S, j in S : i >= j} :=
    if i = j then
        sqrt(Sigma[i,i]-(sum {k in S : k < i} (C[i,k]*C[i,k])))
    else
        (Sigma[i,j]-sum{k in S : k < j} C[i,k]*C[j,k])/C[j,j];

 /* A workaround for the lack of a way to seed the PRNG in GMPL */
param utc := prod {1..2} (gmtime()-1000000000);
param seed := utc - 100000*floor(utc/100000);
check sum{1..seed} Uniform01() > 0;

/* Simulated returns */
param N default 5000;
set T := 1..N;
param R{i in S, t in T} := Rbar[i] + sum {j in S : j <= i} C[i,j]*Normal(0,1);

/* MAD Optimization */

var w{S};                     # Portfolio Weights with Bounds
var y{T} >= 0;                # Positive deviations (non-negative)
var z{T} >= 0;                # Negative deviations (non-negative)

minimize MAD: (1/card(T))*sum {t in T} (y[t] + z[t]);

s.t. C1: sum {s in S} w[s]*Rbar[s] >= Rp;
s.t. C2: sum {s in S} w[s] = 1;
s.t. C3 {t in T}: (y[t] - z[t]) = sum{s in S} (R[s,t]-Rbar[s])*w[s];

solve;

/* Report */

/* Input Data */
printf "\n\nStock Data\n\n";
printf "         Return   Variance\n";
printf {i in S} "%5s   %7.4f   %7.4f\n", i, Rbar[i], Sigma[i,i];

printf "\nCovariance Matrix\n\n";
printf "     ";
printf {j in S} " %7s ", j;
printf "\n";
for {i in S} {
    printf "%5s  " ,i;
    printf {j in S} " %7.4f ", Sigma[i,j];
    printf "\n";
}

/* MAD Optimal Portfolio */
printf "\n\nMinimum Absolute Deviation (MAD) Portfolio\n\n";
printf "  Return   = %7.4f\n",Rp;
printf "  Variance = %7.4f\n\n", sum {i in S, j in S} w[i]*w[j]*Sigma[i,j];
printf "         Weight\n";
printf {s in S} "%5s   %7.4f\n", s, w[s];
printf "\n";

data;

/* Data for monthly returns on four selected stocks for a three
year period ending December 4, 2009 */

# Simulation Horizon
param N := 5000;

# Minimum acceptable investment return
param Rp := 0.01;

# Historical returns on assets
param : S : Rbar :=
    AAPL    0.0308
    GE     -0.0120
    GS      0.0027
    XOM     0.0018 ;
    
# Covariance on asset returns
param   Sigma :
            AAPL    GE      GS      XOM  :=
    AAPL    0.0158  0.0062  0.0088  0.0022
    GE      0.0062  0.0136  0.0064  0.0011
    GS      0.0088  0.0064  0.0135  0.0008
    XOM     0.0022  0.0011  0.0008  0.0022 ;

end;

結果

典型的輸出如下。結果表明，與單個資產的風險回報表現相比，精心設計的投資組合可以表現出明顯改善的風險回報表現。

Stock Data

         Return   Variance
 AAPL    0.0308    0.0158
   GE   -0.0120    0.0136
   GS    0.0027    0.0135
  XOM    0.0018    0.0022

Covariance Matrix

         AAPL       GE       GS      XOM 
 AAPL    0.0158   0.0062   0.0088   0.0022
   GE    0.0062   0.0136   0.0064   0.0011
   GS    0.0088   0.0064   0.0135   0.0008
  XOM    0.0022   0.0011   0.0008   0.0022

Minimum Absolute Deviation (MAD) Portfolio

  Return   =  0.0100
  Variance =  0.0036

         Weight
 AAPL    0.2794
   GE    0.0002
   GS    0.1120
  XOM    0.6084

可能的擴充套件

此示例有許多擴充套件，可以使其在現實世界中更加有用

新增 Cholesky 分解的錯誤檢查。
新增權重係數的上限和下限（例如，沒有賣空）約束。
新增為投資組合指定無風險資產的能力。
新增約束以實施多元化要求。
新增風險回報權衡的引數分析（需要外部指令碼）。
開發一種從 GMPL 內部對 GLPK 的隨機數生成器進行播種的更好方法。
新增一個ODBC 資料庫介面，用於經驗（而不是模擬）歷史收益資料。

參考

↑ Konno, Hiroshi; Yamazaki, Hiroaki (1991). "Mean-absolute deviation portfolio optimization model and its applications to Tokyo stock market". Management Science. 37: 519–531.
↑ Curnuejols, Gerald; Tutuncu, Reha (2007). Optimization Methods in Finance. Cambridge University Press.
↑ Feinstein, Charles D.; Thapa, Mukund N. (1993). "A Reformulation of a Mean-Absolute Deviation Portfolio Optimization Model". Management Science. 39: 1552–1553.

[1] Konno, Hiroshi; Yamazaki, Hiroaki (1991). "Mean-absolute deviation portfolio optimization model and its applications to Tokyo stock market". Management Science. 37: 519–531.

[2] Curnuejols, Gerald; Tutuncu, Reha (2007). Optimization Methods in Finance. Cambridge University Press.

[3] Feinstein, Charles D.; Thapa, Mukund N. (1993). "A Reformulation of a Mean-Absolute Deviation Portfolio Optimization Model". Management Science. 39: 1552–1553.

[1]

[2]

[3]