多目標遺傳演算法NSGA-Ⅱ與其Python實現多目標投資組合優化問題

對於單目標優化問題，一般的遺傳演算法可以較為簡單的得到較好的結果。但是，當問題擴充套件到多目標時，原先的遺傳演算法便不再適用了。因為目標之間通常有著較深的相互關係，一個目標的優化通常會影響到其餘的目標，很難能夠得到所有目標都達到最優的解。這時候，如何尋找合適的適應度函式便成解決多目標遺傳演算法的關鍵。如今，相關的演算法已經有很多種了。包括妥協演算法（compromise approach），GWASF-GA，SPEA2，NSGA-Ⅱ。其中NSGA-Ⅱ的使用非常廣泛。

NSGA-Ⅱ

NSGA-Ⅱ的優點

1.NSGA-Ⅱ提出了快速的非支配(non-dominated)排序，很好的降低了演算法的複雜度。一般的多目標演算法複雜度為$O(MN^3)$，而NSGA-Ⅱ可以做到$O(MN^2)$
2.NSGA-Ⅱ改進了原先NSGA演算法為保留解多樣性而採用的共享函式。提出了擁擠比較運算元（crowded-comparison operator)，從而避免了人為輸入引數的不確定性。

快速非支配排序

快速非支配排序的核心思想主要是通過計算比較得到種群中每個個體p的被支配度$N_p$,通過支配度的大小得到多層非支配曲面。具體來說過程如下：
對於種群中的每一個個體p，我們計算兩個實體。第一個是其被支配度$N_p$，即P個體被其餘個體所支配的數量。支配的定義為如果個體p中所有目標均不優於個體q中對應目標，則稱個體p被個體q所支配。第二個實體是個體p的支配集合$S_p$。這一步所需要的計算複雜度為$O(MN^2)$,因為最壞的情況下，數目為N的種群中每一個個體都要與其餘個體比較，這一步為$N^2$，那麼對於$M$個目標則為$MN^2$。
接下來可以開始尋找非支配曲面了。對於最優非支配曲面(Pareto-optimal front)，其中的個體$N_p$為0。接著，對於最優非支配曲面中的每一個個體，尋找其相應的$S_p$，對於其中所有的個體q，將其$N_p$減1。對於此時$N_p$為0的個體，我們將其歸入集合Q，Q便是第二非支配曲面。按照相同的步驟，我們可以得到所有的非支配曲面。由於對於所有個體p，其$N_p$最大為N-1，因此，考慮最壞的情況，即最大$N_p$為N-1，其非支配曲面的個數為N，則為了得到N個曲面所需要的最大遍歷的次數為$N^2$。這樣，對於M個目標，其計算複雜度為$O(MN^2)$。快速非支配排序用演算法表達如下：

擁擠比較運算元

一般來說，在尋求多目標遺傳演算法的解時，我們希望解集合能夠儲存多樣性。在NSGA-Ⅱ中，通過運用擁擠比較運算元來實現此目的。簡單來說，擁擠比較運算元的原理如下：
首先定義密度衡量準則（density estimation metric），密度衡量準則定義為一個個體在其每個目標上上下相鄰兩點的距離。即某個個體與其附近的個體相距越遠，密度越小。定義擁擠距離$i_{distance}$定量計算每個個體的密度。擁擠距離的計算需要首先對同一條非支配曲面中每個個體在每個目標上進行升序排序，其中最大值與最小值的個體在此目標上的擁擠距離設為無窮。個體$i$的擁擠距離即等於它在每個目標上$i-1$與$i+1$個體的函式值之差的絕對值之和。在計算時需要標準化每一個目標。兩目標的密度表示圖如下：

擁擠距離的計算用演算法流程表示如下：

有了擁擠距離，便可以定義擁擠比較運算元。簡單來說，此運算元的計算流程如下。對於越優的非支配曲面中的個體，其適應度函式值越大。而對於同一非支配曲面中的個體，根據擁擠距離確定適應度函式，擁擠距離越大，密度越小，適應度函式值越高。論文中描述如下：

NSGA-Ⅱ演算法流程

在解釋完快速非支配排序演算法和擁擠比較運算元後，我們便可以設計整個演算法流程了。首先定義引數：定義$P_i$為第$i$代種群，而$P_0$為初始種群。種群中個體數為$N$。而$Q_i$為通過遺傳演算法產生的第$i$子代，數目也為$N$。演算法描述如下：

1.隨機產生初始代種群$P_0$，個體數為$N$
2.通過交叉運算元與遺傳運算元得到子代$Q_0$
3.合併種群$P_0$與$Q_0$，得到個數為$2N$的種群。
4.通過快速非支配排序演算法求得種群中每個個體的適應度函式。
5.通過擁擠比較運算元執行自然選擇過程。同時，在選擇過程中引入精英機制（elitism），保留排在前面的非支配曲面。具體來說，首先保留最優非支配曲面的所有解，如果個數小於N，則繼續保留第二非支配曲面的解，以此類推，直至無法保留此非支配曲面的所有解，即保留總數大於N，則運用擁擠比較運算元選取適應度函式值更優的解，直至總數達到N，完成選擇過程。
6.得到個體數為N的下一代種群$P_1$
7.按照模型的規定進化次數重複執行2-6步，直至完成演算法。

到此為止，NSGA-Ⅱ演算法就闡述完了，如果想要更深入的瞭解此演算法，可以自行搜尋論文$閱讀^{[1]}$。接下來，我們嘗試用python實現演算法並運用到實際中解決特定的多目標優化問題。

Python實現NSGA-Ⅱ演算法並解決多目標投資組合優化問題

模型描述

均值-方差理論自從馬科維茨在1959年提出後，逐漸成為了現代投資理論的基石。在均值-方差理論中，馬科維茨將資產的收益率定義為其收益率分佈的均值，而將其收益率分佈的方差定義為資產的風險，而資產組合的優化問題被定義為尋找收益率均值高同時方差小的投資組合。近幾十年，均值-方差理論得到了多方面的發展，越來越多的目標被考慮進來，比如收益率分佈的偏度與峰度，資產的流動性，熵，VaR，換手率等；另一方面，基於可信度與可能性理論的投資組合優化理論也得到充分的研究與發展。
本篇文章中， 我們具體考慮均值-方差-換手率（expectation-variance-turnover rate）投資組合模型。定義投資組合中資產$i$的比例為$x_i$，投資組合中總資產數為$n$。同時，為了模擬更為真實的投資環境，我們禁止買空操作，同時，加入Yager熵用於分散投資，並且設定投資組合中每個資產的投資上限與下限為$u$和$l$。關於這部分內容，由於不是本文重點，就不做具體闡述，有興趣的讀者可以自行查閱相關資料。這樣，我們的投資組合優化模型可以表示為：

很明顯，這是一個三目標的優化問題，一般需要運用啟發式演算法求解，這裡我們採用上文闡述的NSGA-Ⅱ進行求解。

資料選取

這裡我們選用深交所與上交所上市的12只股票在2014年至2016年的日交易資料，資料來源為同花順。股票選取如下：

原始資料格式如下：

資料處理

接下來通過dataProcessing.py對資料進行處理：

"""
data process
"""
import pandas as pd

#define function to output processed file
def outputFile(data,source):
    data.to_excel(source+'.xlsx',encoding='utf-8')

name=['000002','600690','002001','600009','000001','002008','002236','002384','002304','600885','000046','000858']
suffix='_price'

returnData=pd.DataFrame()
turnoverData=pd.DataFrame()
for i in range(len(name)):
    data=pd.read_excel('../rowData/'+name[i]+suffix+'.xlsx')
    data.index=pd.to_datetime(data['時間'])
    returnData[name[i]]=data['漲幅']
    turnoverData[name[i]]=data['換手%']

returnData=returnData.fillna(0)
turnoverData=turnoverData.fillna(0)
outputFile(returnData,'../processedData/returns')
outputFile(turnoverData,'../processedData/turnovers')

整合所選股票日收益率與換手率資料，對於資料缺失的部分，一般為股票停牌，我們用0代替。

NSGA-Ⅱ

在NSGA2Selection.py檔案中，我們首先定義3個目標的計算函式：

"""
NSGA2 selection
"""
import pandas as pd
import numpy as np

#compute the expectation of the returns of the portfolio
def expectation(returns,weights):
    eArray=np.array(returns.mean())
    return np.dot(eArray,weights.T)

#compute the variance of the returns of the portfolio
def variance(returns,weights):    
    return np.dot(weights.T, np.dot(returns.cov()*len(returns.index),weights))
    
#compute the skewness of the returns of the portfolio
def turnover(turnovers,weights):
    tArray=np.array(turnovers.mean())
    return np.dot(tArray,weights.T)

接著，定義函式用於將計算結果儲存：

#generate value set for every objective
def objectiveValueSet(population,populationSize,fun,data):
    valueList=[0 for i in range(populationSize)]
    for i in range(populationSize):
        valueList[i]=fun(data,np.array(population[i][0]))
    return valueList

做完準備工作後，便可以編寫快速非支配排序演算法了：

#Fast Nondominated Sorting Approach
def NSGA(population,populationSize,returns,turnovers):
    #build list and compute return,variance,turnover rate for every entity in the population
    returnList=objectiveValueSet(population,populationSize,expectation,returns)
    varianceList=objectiveValueSet(population,populationSize,variance,returns)
    turnoverList=objectiveValueSet(population,populationSize,turnover,turnovers)
    
    #define the dominate set Sp
    dominateList=[set() for i in range(populationSize)]
    #define the dominated set
    dominatedList=[set() for i in range(populationSize)]
    #compute the dominate and dominated entity for every entity in the population
    for i in range(populationSize):
        for j in range(populationSize):
            if returnList[i]> returnList[j] and varianceList[i]<varianceList[j] and turnoverList[i]>turnoverList[j]:
                dominateList[i].add(j)
                
            elif returnList[i]< returnList[j] and varianceList[i]>varianceList[j] and turnoverList[i]<turnoverList[j]:
                dominatedList[i].add(j)
    #compute dominated degree Np
    for i in range(len(dominatedList)):
        dominatedList[i]=len(dominatedList[i])
    #create list to save the non-dominated front information
    NDFSet=[]
    #compute non-dominated front
    while max(dominatedList)>=0:
        front=[]
        for i in range(len(dominatedList)):
            if dominatedList[i]==0:
                front.append(i)
        NDFSet.append(front)
        for i in range(len(dominatedList)):
            dominatedList[i]=dominatedList[i]-1                                
    return NDFSet

其中，population是種群儲存當前種群資訊的list，大小為populationSize×2×12，其中二階的兩個長度為12的list分別用於儲存12只股票的投資比例和是否投資的資訊，具體定義方法將在後面提到。
dominatedList與dominateList分別用於儲存種群中每個個體的支配集合與被支配度。使用set()的原因是為了避免重複新增。
最後得到的NDFSet儲存著非支配曲面的資訊。注意到其中會存在部分空list，這是由於支配度相差大於1所導致的，但是不影響後續建模。


接著編寫擁擠距離函式：

#crowded distance
def crowdedDistance(population,Front,minmax,returns,turnovers):
    #create distance list to save the information of crowded for every entity in front
    distance=pd.Series([float(0) for i in range(len(Front))], index=Front)
    #save information of weight for every entity in front
    FrontSet=[]
    for i in Front:
        FrontSet.append(population[i])

    #compute and save objective value for every entity in front
    returnsList=objectiveValueSet(FrontSet,len(FrontSet),expectation,returns)
    varianceList=objectiveValueSet(FrontSet,len(FrontSet),variance,returns)
    turnoverList=objectiveValueSet(FrontSet,len(FrontSet),turnover,turnovers)   
    returnsSer=pd.Series(returnsList,index=Front)
    varianceSer=pd.Series(varianceList,index=Front)
    turnoverSer=pd.Series(turnoverList,index=Front)
    #sort value
    returnsSer.sort_values(ascending=False,inplace=True)
    varianceSer.sort_values(ascending=False,inplace=True)
    turnoverSer.sort_values(ascending=False,inplace=True)
    #set the distance for the entities which have the min and max value in every objective 
    distance[returnsSer.index[0]]=1000
    distance[returnsSer.index[-1]]=1000
    distance[varianceSer.index[0]]=1000
    distance[varianceSer.index[-1]]=1000
    distance[turnoverSer.index[0]]=1000
    distance[turnoverSer.index[-1]]=1000
    for i in range(1,len(Front)-1):
        distance[returnsSer.index[i]]=distance[returnsSer.index[i]]+(returnsSer[returnsSer.index[i-1]]-returnsSer[returnsSer.index[i+1]])/(minmax.iloc[0,1]-minmax.iloc[0,0])
        distance[varianceSer.index[i]]+=(varianceSer[varianceSer.index[i-1]]-varianceSer[varianceSer.index[i+1]])/(minmax.iloc[1,1]-minmax.iloc[1,0])
        distance[turnoverSer.index[i]]+=(turnoverSer[turnoverSer.index[i-1]]-turnoverSer[turnoverSer.index[i+1]])/(minmax.iloc[2,1]-minmax.iloc[2,0])
    distance.sort_values(ascending=False,inplace=True)
    return distance

這裡由於使用無窮大不方便計算，使用1000作為無窮大的代替。minmax裡儲存著每個目標的最大值與最小值，用於標準化每個目標值，對於一階矩，可直接進行計算，對於二階矩方差使用單目標遺傳演算法得到近似結果，具體程式碼便不贅述了，具體可仿照本人之前的文章編寫。


有了計算擁擠距離的函式，接下來編寫擁擠比較運算元：

def crowdedCompareOperator(population,populationSize,NDFSet,minmax,returns,turnovers):
    newPopulation=[]
    #save the information of the entity the new population have
    count=0
    #save the information of the succession of the front 
    number=0
    while count<populationSize:            
        if count + len(NDFSet[number])<=populationSize:
            if number==0:
                #save the information of the first non-dominated front
                firstFront=[i for i in range(len(NDFSet[number]))]
            for i in NDFSet[number]:
                newPopulation.append(population[i])
            count+=len(NDFSet[number])
            number+=1
        else:
            if number==0:
                firstFront=[i for i in range(populationSize)]
            n=populationSize-count
            distance=crowdedDistance(population,NDFSet[number],minmax,returns,turnovers)
            for i in range(n):
                newPopulation.append(population[distance.index[i]])
            number+=1
            count+=n

    return newPopulation,firstFront

這裡firstFront用於儲存自然選擇後種群中處於最優非支配曲面上的個體，因為最後進行分析時一般只需要用到最優分支配曲面。


接下來編寫選擇主函式：

def NSGA2Selection(population,populationSize,minmax,returns,turnovers):
    NDFSet=NSGA(population,populationSize*2,returns,turnovers)
    newPopulation,firstFront=crowdedCompareOperator(population,populationSize,NDFSet,minmax,returns,turnovers)
    return newPopulation,firstFront,NDFSet

有了NSGA2選擇運算元，我們接著編寫遺傳演算法主函式檔案。不過首先我們在tools.py定義工具函式

"""
function tools
"""


#initiate stock list signal
def stocksignal(population,populationSize):
    for i in range(populationSize):
        stocklist=[]
        for j in range(12):
            if population[i][0][j]>0:
                stocklist.append(1)
            else:
                stocklist.append(0)
        if len(population[i])==1:
            population[i].append(stocklist)
        else:
            population[i][1]=stocklist
    return population

def yagersEntropy(porprotion,n):
    H=0
    for i in range(n):
        H+=abs(porprotion[i]-1/n)
        
    return H

其中stocksignal函式用於計算是否投資某個具體股票，投資則為1，為投資則為0，此資訊在交叉與遺傳運算元中將會用到。yagersEntropy函式用於計算投資組合的雅各熵。


編寫主函式檔案multiObjective.py,首先匯入各運算元模組並定義引數：

import pandas as pd
import tools
import random
from crossover import crossover
from mutation import mutation
from repair import repair
from NSGA2Selection import NSGA2Selection

returns=pd.read_excel('../processedData/returns.xlsx')
turnovers=pd.read_excel('../processedData/turnovers.xlsx')
minmax=pd.read_excel('../processedData/minmaxValue.xlsx')

populationSize=50
generation=50
Pm=0.5
Pc=0.9
l=0.05
u=0.3

crossover與mutation模組為交叉與變異運算元，但由於研究問題不公佈在上面，有興趣的讀者可以自行查閱文獻按照相關規則設計。repair模組為修復運算元，用於在交叉與變異後調整比例使得投資組合重新滿足各項限制（如最大最小值，雅閣熵），由於研究問題同樣就不公佈了，設計起來其實也不是很難，相關文獻也有很多。
Pm為變異概率，為了提高種群多樣性與避免區域性最小值，將變異概率調高至0.5；Pc為交叉概率，l與u分別為每隻股票所能投資的下限與上限。


接著初始化種群：

#initiate population    
population=[[] for i in range(populationSize)]
selectedPopulation=[[] for i in range(populationSize)]
for i in range(populationSize):
    entity=[0 for i in range(12)]
    n=1
    label=[i for i in range(12)]
    for j in range(12):
        if n>u:
            t=random.choice(label)
            entity[t]=random.uniform(l,u)
            n-=entity[t]
            label.remove(t)
        elif n<l:            
            entity[t]+=n
            n=0
        else:
            t=random.choice(label)
            entity[t]=random.uniform(l,n)
            n-=entity[t]         
            label.remove(t)              
    population[i].append(entity)
    selectedPopulation[i].append(entity)
newPopulation=tools.stocksignal(population,populationSize)

執行迴圈：

#GA loop
for i in range(generation):
    offspring=crossover(newPopulation,populationSize,Pc,l,u)
    offspring=mutation(offspring,populationSize,Pm)
    offspring=repair(offspring,populationSize,l,u)
    selectPopulation=selectedPopulation+offspring
    selectedPopulation,firstFront,NDFSet=NSGA2Selection(selectPopulation,populationSize,minmax,returns,turnovers)
    newPopulation=tools.stocksignal(selectedPopulation,populationSize)
    print (i)

如果執行迴圈可以發現，最優非支配面上的個體將逐漸增加，證明通過遺傳演算法種群在不斷進行優化。


完成迴圈後，儲存最優非支配面上個體的資訊並儲存以做後續分析：

name=['000002','600690','002001','600009','000001','002008','002236','002384','002304','600885','000046','000858']
result=pd.DataFrame([[0 for j in range(len(name))] for i in range(len(firstFront))])
for i in range(len(firstFront)):
    for j in range(len(name)):
        result.iloc[i,j]=newPopulation[firstFront[i]][0][j]

result.columns=name
result.to_excel('../resultData/result2.xlsx',encoding='utf-8')

這樣，我們便使用NSGA-Ⅱ演算法解決了多目標投資組合優化問題，對於得到的結果，我們可以通過運用歷史資料進行回測評價效果，具體便不贅述。
相關程式碼與資料已經上傳至Github，由於編寫倉促，程式碼有許多不完善的地方，歡迎大家進行討論！

[1]Deb, Kalyanmoy, Pratap, Amrit, Agarwal, Sameer, et al. A fast and elitist multiobjective genetic algorithm: NSGA-II[J]. IEEE Transactions on Evolutionary Computation, 2002, 6(2):182-197.