# 多目標遺傳演算法優化 S&P 500 股票預測模型的特徵選擇
## 摘要
- 這篇研究主要基於多目標遺傳演算法，用於選擇最佳的股票特徵指標並應用於美國股市預測。它透過演化過程找到同時最佳化波動度和報酬率的股票特徵，並將這些特徵應用於預測股票價格的波動性和報酬率，主要通過多目標基因演算法和機器學習模型，並對資料進行了特徵選擇和預測。最終，根據波動度和報酬率將股票分為不同的群組，建立投資組合。並觀察總體績效和風險，我發現在2021年至2023年間，由基因演算法篩選的股票特徵能有效建構低波動率和高報酬率的投資組合，在所有投資組合中計算了投組報酬率、最大回落比和Calmar比率，大部分報酬率與風險皆優於sp500指數，總體來說，這個研究提供了一個基於多目標遺傳演算法的股票特徵選擇和投資組合建構方法，這對於投資者在股市中做出更明智的決策和配置資源非常有價值。
## 文獻回顧
- 股票市場的預測一直是投資者和研究人員關注的焦點之一。近年來，隨著技術分析和機器學習技術的發展，論文開始探索使用股票特徵選股的方法，以提高股票預測模型的準確性和效能，主要透過不同模型，或是最近流行的機器學習模型或深度學習，試圖對資本市場中的資產進行正確的預測，但往往在資本市場中，資產會有多種性質，因此如何選定特徵進行分析為主要課題，過去文獻很多都對價量資料進行討論，Tsantekidis et al. (2017)利用卷積神經網絡 (Convolutional Neural Networks, CNN) 從限價單簿 (limit order book) 預測股票價格。研究結果顯示，使用CNN模型能夠在股票價格預測方面取得良好的準確性，Guo et al. (2019)提出了一種基於多目標深度強化學習的股票選擇方法。通過在強化學習框架下設定多個目標，並建立投資組合。過往也有許多文獻對股票特徵進行研究，Zhang et al. (2019)探討了股票市場分析中的特徵選擇和分類方法。作者提出了一種基於綜合權重的特徵選擇方法，並結合支持向量機 (Support Vector Machine, SVM) 進行分類分析，這種方法能夠選擇具有關鍵影響力的特徵，Han et al. (2019)認為特徵選擇對於提高股票市場預測準確性具有重要作用，並且在不同市場條件下的表現存在差異。本研究也試圖以Sharpe ratio為基本概念建立衡量股票價值，並透過多個目標最佳化找出優秀的股票特徵因子，希望能找到在影響標準普爾500指數成分股的關鍵因子，並透過這些關鍵因子建立投資組合

- 股票市場的預測一直是投資者和研究人員關注的焦點之一。近年來，隨著技術分析和機器學習技術的發展，學術界開始探索使用股票特徵選股的方法，以提高股票預測模型的準確性和效能。主要的研究方法包括使用不同的模型，例如流行的機器學習模型或深度學習方法，以正確預測資本市場中的資產價格。然而，在資本市場中，不同資產具有多種性質，因此如何選擇適當的特徵進行分析成為主要的課題。過去的研究文獻中，Tsantekidis等人（2017）利用卷積神經網絡（CNN）從限價單簿（limit order book）預測股票價格。研究結果顯示，使用CNN模型能夠在股票價格預測方面取得良好的準確性。Guo等人（2019）提出了一種基於多目標深度強化學習的股票選擇方法。該方法在強化學習框架下設定多個目標並建立投資組合，以達到多個目標的最佳化。此外，Zhang等人（2019）探討了股票市場分析中的特徵選擇和分類方法。他們提出了一種基於綜合權重的特徵選擇方法，並結合支持向量機（SVM）進行分類分析。這種方法能夠選擇具有關鍵影響力的特徵。Han等人（2019）認為特徵選擇對於提高股票市場預測準確性具有重要作用，並且在不同市場條件下的表現存在差異。本研究試圖以Sharpe ratio為基本概念，建立衡量股票價值的方法，並透過多目標最佳化找出優秀的股票特徵因子。研究目標為找到影響標準普爾500指數成分股的關鍵因子，並利用這些關鍵因子建立一個優秀的投資組合。研究旨在提供投資者更具價值的股票選擇策略，並進一步優化投資組合的風險收益特性。這些研究結果對於投資者在股票市場中做出明智的選擇和決策提供了重要的參考依據。

- 神經網路本文回顧了幾篇相關的論文，並結合我自身的研究，探討了股票特徵選股在預測模型中的應用。Abdullah et al. (2014)介紹了結合技術分析和機器學習技術的方法，用於預測股票價格。該研究評估了不同的技術指標的效能，並發現這種結合方法能夠提高預測準確性。Perera et al. (2018)提供了對金融時間序列預測中機器學習技術的綜述，並探討了不同機器學習技術在股票預測中的應用。Karakaya and Cevikcan (2018)介紹了不同的特徵選擇技術，並討論了其在股票預測中的應用，過去的文獻回顧大多致力於透過不同面向的分析，總體、個體等等的關係，探討如何透過股票特徵去最佳化預測準確性，本研究更進一步透過sp500中250檔股票，試圖通過多目標演算法建構有效股票特徵因子，並透過這些因子進行投資組合的建構。

## 資料與方法
- 這邊使用的股價資料為yifinance股價資料庫，總共抓取了250檔股票，從2019/01/01至2023/04/30資料，資料分兩部分，2019/01/01至2019/10/30資料為訓練模型與多目標演算法部分，剩下的部分為測試投資組合績效期。每一檔股票，總共250檔，每一檔股票都會計算180種股票特徵，建構股票特徵資料，股票特徵主要分為10種類別，分別為價量資料、週期指標、數值變換、動量指標、圖表形態識別、價格變換、統計函數、波動性指標、交易量指標、重疊研究指標，其中有些指標因為數值過大，例如指數函數，故刪除，最後共選出了174個股票特徵。
- 我主要透過多目標基因演算法，用於解決特徵選擇的最佳化問題。該問題的目標是從一組股票特徵中選擇最佳的子集，以用於股票價格預測。以下是該程式的主要組件和操作的詳細介紹，我將此演算法適應度函數設定為隨機森林模型測試集的的均方誤差（Mean Squared Error），並將最佳化設定為極小化兩個隨機森林模型測試集的均方誤差，分別為訓練波動度和訓練收益率的模型。最後，基因演算法的基因設定為由0和1組成的二進制串列，代表是否選擇此股票特徵作為預測因子，基因數量設定為20，迭代次數為50，交叉機率為0.5、突變機率為0.2
- 通過多目標基因演算法篩選出最佳股票特徵後，透過這些特徵作為預測股票波動度和報酬率的因子，建構未來預測波動度和報酬率，並從原本250檔股票中篩選出排名前30%且波動率為後30%的股票，建立低波動高報酬率投資組合，下一章會進行實證研究。

### 備註
- 刪除cos、exp、ad、adosc、obv、sinh

| 股票特徵種類 | 該類別數量 |
| ----------------------- | ---------- |
| Price Volume Indicators | 6 |
| Cycle Indicators | 5 |
| Math Operators | 11 |
| Math Transform | 15 |
| Momentum Indicators | 30 |
| Pattern Recognition | 61 |
| Price Transform | 4 |
| Statistic Functions | 9 |
| Volatility Indicators | 3 |
| Volume Indicators | 3 |
| Overlap Studies | 17 |
- Cycle Indicators:週期性因素的分析(HT_TRENDLINE、KAMA、TRIX、DPO)
- Math Transform:測量資產價格的速度和幅度變化。它可以顯示股票或其他資產價格在一段時間內的變化趨勢(正弦（Sin）和餘弦（Cos）函數、自然對數（Ln）、算術平方根（Sqrt）)
- Momentum Indicators:測量資產價格的速度和幅度變化。它可以顯示股票或其他資產價格在一段時間內的變化趨勢(RSI、ROC、CMF、CCI、ADX)
- Pattern Recognition:提供了多種圖表形態的技術指標，可以用來偵測股票價格走勢中常見的模式
- Price Transform:Weighted Close Price（WCLPRICE）、Median Price（MEDPRICE）
- Statistic Functions:CORREL、LINEARREG、STDEV、VAR
- Volatility Indicators:衡量股票或市場波動性(Average True Range (ATR)、Bollinger Bands、Chaikin's Volatility (CHV))
- Volume Indicators:AD、CMF、OBV、VWAP、PVI
- Overlap Studies:用過去一段時間的股價資訊，計算出一個或多個股價的平均值或變動率，以及相應的上下限。這些指標可以用來幫助分析股價的走勢、支撐位和阻力位(EMA、MACD、SMA)

## 實證結果

### 股票特徵篩選
- 資料期間:2019/01/01~2019/10/30
- 資料類別:250檔sp500上市股票
- 花費時間:684分鐘43.1秒
- 篩選出最適股票特徵，用於預測未來波動度和報酬率的因子
## 建立投資組合
### 第一組投組
1. 訓練+測試集:2019/10/31 ~ 2020/12/31，並計算每支股票平均波動度、平均報酬率
2. 建構2021/01/01~2021/04/30投組
3. 計算投組報酬率、最大回落比、calmar ratio
4. 並與持有大盤做比較
5. 持有投組:共25檔

| 持有標的 | | | | |
| --- | -------- | ---- | ---- | --- |
| MRK| JNJ | PG | WMT | T |
| CVX| VZ | NKE| LMT| MO |
| GILD| SPGI | CME| CCI | D |
| WFC | BAX | FOX| BXP | WBA |
| MMC| DLR | MSCI | RCL | NDAQ|

7. 投組績效:
8. sp500績效

| portfolio return & risk | portfolio | sp500 |
| -------- | -------- | -------- |
| return| :apple:6.6% | 6.53% |
| mdd | :apple:3.71% | 6.02% |
| calmar ratio| :apple:1.78 | 1.08 |

### 第二組投組
1. 訓練+測試集:2020/02/25 ~ 2021/4/30，並計算每支股票平均波動度、平均報酬率
2. 建構2021/05/01~2021/08/31投組
3. 計算投組報酬率、最大回落比、calmar ratio
4. 並與持有大盤做比較
5. 持有投組:共17檔

| 持有標的 | | | | |
| --- | -------- | ---- | ---- | --- |
| AAPL| ADBE | CRM | AMD | AMGN |
| HON| APD | EW| DD| AEP |
| REGN| BSX | PGR| CNC | ANET |
| VRSN | AME | | | |
7. 投組績效:

8. sp500績效

| portfolio return & risk | portfolio | sp500 |
| -------- | -------- | -------- |
| return| :apple:10.54% | 6.54% |
| mdd | :apple:2.03% | 2.99% |
| calmar ratio| :apple:5.2 | 2.19 |

### 第三組投組
1. 訓練+測試集:2020/06/25 ~ 2021/8/31，並計算每支股票平均波動度、平均報酬率
2. 建構2021/09/01~2021/12/31投組
3. 計算投組報酬率、最大回落比、calmar ratio
4. 並與持有大盤做比較
5. 持有投組:共16檔

| 持有標的 | | | | |
| --- | -------- | ---- | ---- | --- |
| DIS| VZ | T | INTC | CVS |
| VRTX| BAX | FOX| ECL| ANSS |
| STZ| KMB | CDW| AIG | MCHP |
| PENN | | | | |
6. 投組績效:
7. sp500績效:

| portfolio return & risk | portfolio | sp500 |
| -------- | -------- | -------- |
| return| 5.99% | :apple:9.62% |
| mdd | 5.37% | :apple:3.06% |
| calmar ratio| 1.12 | :apple:3.15 |

### 第四組投組
1. 訓練+測試集:2020/10/25 ~ 2021/12/31，並計算每支股票平均波動度、平均報酬率
2. 建構2022/01/01~2022/04/15投組
3. 計算投組報酬率、最大回落比、calmar ratio
4. 並與持有大盤做比較
5. 持有投組:共21檔

| 持有標的 | | | | |
| -------- | --- | ---- | ---- | ---- |
| TSM | MDT | C | LMT | GILD |
| IBM | BA | BDX | SYK | D |
| DUK | AEP | FOX | SWKS | FISV |
| ALL | HIG | PEAK | PSX | RCL |
| PENN | | | | |

6. 投組績效:
7. sp500績效:

| portfolio return & risk | portfolio | sp500 |
| -------- | -------- | -------- |
| return| :apple:5.12% | -2.62% |
| mdd | :apple:5.0% | 10.74% |
| calmar ratio| :apple:1.02 | -0.24 |

### 第五組投組
1. 訓練+測試集:2021/2/25 ~ 2022/4/15，並計算每支股票平均波動度、平均報酬率
2. 建構2022/4/16~2022/08/16投組
3. 計算投組報酬率、最大回落比、calmar ratio
4. 並與持有大盤做比較
5. 持有投組:共16檔

| 持有標的 | | | | |
| -------- | --- | ---- | ---- | ---- |
| C | BA | GM | SPG | EBAY |
| FISV | WBA | VFC | DLR | ALB |
| WY | CTAS | KMB | COF | IR |
| DLB | | | | |

6. 投組績效:
7. sp500績效:

| portfolio return & risk | portfolio | sp500 |
| -------- | -------- | -------- |
| return| :apple:7.49% | 3.83% |
| mdd | :apple:11.91% | 13.7% |
| calmar ratio| :apple:0.63 | 0.28 |

### 第六組投組
1. 訓練+測試集:2021/6/25 ~ 2022/8/16，並計算每支股票平均波動度、平均報酬率
2. 建構2022/8/17~2022/12/31投組
3. 計算投組報酬率、最大回落比、calmar ratio
4. 並與持有大盤做比較
5. 持有投組:共19檔

| 持有標的 | | | | |
| -------- | ---- | ---- | ---- | ---- |
| MA | CSCO | AMD | QCOM | HON |
| AXP | GM | DD | LUV | ROST |
| BKNG | ECL | ADI | HLT | DXCM |
| DOW | AIG | MCHP | STE | |
6. 投組績效:
7. sp500績效:

| portfolio return & risk | portfolio | sp500 |
| -------- | -------- | -------- |
| return| :apple:-1.08% | -8.86% |
| mdd | 18.35% | :apple:17.15% |
| calmar ratio| :apple:-0.06 | -0.52 |

### 第七組投組
1. 訓練+測試集:2021/10/25 ~ 2022/12/31，並計算每支股票平均波動度、平均報酬率
2. 建構2023/1/1~2023/4/30投組
3. 計算投組報酬率、最大回落比、calmar ratio
4. 並與持有大盤做比較
5. 持有投組:共19檔

| 持有標的 | | | | |
| -------- | ---- | --- | --- | ---- |
| GOOGL | GOOG | AMD | PM | FIS |
| AXP | SPGI | PLD | CME | EW |
| EBAY | IQV | EXC | BXP | ANSS |
| KEY | PSA | DOW | ESS | PKI |
| HAS | ZBRA | | | |
6. 投組績效:

7. sp500績效:

| portfolio return & risk | portfolio | sp500 |
| -------- | -------- | -------- |
| return| :apple:9.85% | 8.97% |
| mdd | 8.66% | :apple:6.18% |
| calmar ratio| 1.14 | :apple:1.45 |

## conclusion
### 回測期間:2021/01/01 ~ 2023/4/30
### 每四個月調整投資組合持股，共調整6次，7個投資組合，績效如下
:apple:**:代表優於sp500指數績效**
**投組一代表2021/01/01持有投組，依此類推**

| 投資組合 | return | mdd | calmar ratio |
| -------- | ------ | --- | ------------ |
| 投組1 | :apple:6.6% | :apple:3.71% | :apple:1.78 |
| 投組2 | :apple:10.54% | :apple:2.03% | :apple:5.2 |
| 投組3 | 5.99% | 5.37% | 1.12 |
| 投組4 | :apple:5.12% | :apple:5.0% | :apple:1.02 |
| 投組5 | :apple:7.49% | :apple:11.91% | :apple:0.63 |
| 投組6 | :apple:-1.08% | 18.35% | :apple:-0.06 |
| 投組7 | :apple:9.85% | 8.66% | 1.14 |

:crown:**:代表優於投資組合績效**
| 投資組合 | return | mdd | calmar ratio |
| -------- | ------ | --- | ------------ |
| sp500 | 6.53% | 6.02% | 1.08 |
| sp500 | 6.54% | 2.99% | 2.19 |
| sp500 | :crown: 9.62% | :crown: 3.06% | :crown: 3.15 |
| sp500 | -2.62% | 10.74% | -0.24 |
| sp500 | 3.83% | 13.7%| 0.28 |
| sp500 | -8.86% | :crown: 17.15% | -0.52 |
| sp500 | 8.97% | :crown: 6.18% | :crown: 1.45 |

- 根據上述績效比較結果，可以發現在大多數測試期間（其中七次調整投資組合中的五次），投資組合在報酬率和風險（Calmar比率）衡量方面均優於標準普爾500指數。這意味著通過多目標最佳化策略，我成功地篩選出了最佳的股票特徵，並構建了一系列優秀的投資組合。這些投資組合能夠有效地過濾低報酬和高風險的股票，從而實現了在大部分實證期間擊敗指數的目標。這一成果的達成，使整體投資組合相較於指數表現出更高的報酬率和更低的風險水平。觀察這七個投資組合的走勢，可以發現它們與標準普爾500指數之間存在高度相關性。這意味著即使在追求優異績效的同時，我依然保持了與市場整體走勢的一定一致性。這種相關性的存在，使得我的投資策略在市場波動期間能夠更好地應對風險，並且在市場上表現出較為穩定的投資回報。

#抓美股資料
import pandas_datareader.data as pdr
import datetime
import pandas as pd

#技術指標
import talib
all_ta_label = talib.get_functions()
len(all_ta_label) #有158類
#分類
all_ta_groups = talib.get_function_groups()
all_ta_groups.keys()
all_ta_groups['Cycle Indicators']
table = pd.DataFrame({'技術指標類別名稱:':list(all_ta_groups.keys()),
                      '該類別指標總數:':list(map(lambda x: len(x), all_ta_groups.values()))})

#Abstract API
from talib import abstract
#打包成函數
def get_stock_data(stock_code, start_date, end_date):
    #取得股票資料
    df = pdr.get_data_yahoo([stock_code], start_date, end_date)
    df.columns = [col.lower() for col in df.columns]
    ta_list = talib.get_functions()
    for x in ta_list:
        try:
            # x 為技術指標的代碼，透過迴圈填入，再透過 eval 計算出 output
            output = eval('abstract.'+x+'(df)')
            # 如果輸出是一維資料，幫這個指標取名為 x 本身；多維資料則不需命名
            output.name = x.lower() if type(output) == pd.core.series.Series else None
            # 透過 merge 把輸出結果併入 df DataFrame
            df = pd.merge(df, pd.DataFrame(output), left_on = df.index, right_on = output.index)
            df = df.set_index('key_0')
        except:
            print(x)
    # 將股票價格和報酬率向後滯後一期
    df['return'] = df['adj close'].pct_change()
    df['change'] = (df['adj close'] > df['open']).astype(int)
    # 刪除包含缺失值的行
    data = df
    data = data.drop(['acos', 'asin'], axis=1)
    data = data.dropna()
    data = data.reset_index(drop=True)
    data = data.astype('float')

    return data

#寫機器學習函數
from sklearn.ensemble import RandomForestClassifier
from sklearn.metrics import roc_auc_score
from sklearn.model_selection import train_test_split
from sklearn.preprocessing import StandardScaler
from sklearn.ensemble import RandomForestRegressor
from sklearn.metrics import mean_squared_error
import numpy as np
import copy
def predict_stock_price_ver1(data, keep_cols):
    # 根据 keep_cols 的值删除需要删除的列
    # 用除了漲跌的欄位作為自變數，去預測股價漲跌
    df = copy.deepcopy(data)
    df.fillna(0, inplace=True)
    df['stddev_lagging'] = df['stddev'].shift(-1)
    df = df.dropna()
    scaler = StandardScaler()
    X = df.drop(['stddev_lagging'], axis=1)
    drop_cols = [col for i, col in enumerate(X.columns) if not keep_cols[i]]
    X = X.drop(drop_cols, axis=1)
    X = X.dropna()
    X = scaler.fit_transform(X)
    y = df['stddev_lagging'].values.reshape(-1,1)
    y = scaler.fit_transform(y)
    # 分割資料集為訓練集和測試集
    X_train, y_train = X, y
    # 創建隨機森林模型
    rf = RandomForestRegressor(random_state=10)
    rf.fit(X_train, y_train)
    y_pred = rf.predict(X_train)
    # 計算測試集的 AUC 分數
    mse_score = mean_squared_error(y_train, y_pred)
    return mse_score
def predict_stock_price_ver2(data, keep_cols):
    # 根据 keep_cols 的值删除需要删除的列
    # 用除了漲跌的欄位作為自變數，去預測股價漲跌
    df = copy.deepcopy(data)
    df.fillna(0, inplace=True)
    df['return_lagging'] = df['return'].shift(-1)
    df = df.dropna()
    scaler = StandardScaler()
    X = df.drop(['return_lagging'], axis=1)
    drop_cols = [col for i, col in enumerate(X.columns) if not keep_cols[i]]
    X = X.drop(drop_cols, axis=1)
    X = X.dropna()
    X = scaler.fit_transform(X)
    y = df['return_lagging'].values.reshape(-1,1)
    y = scaler.fit_transform(y)
    # 分割資料集為訓練集和測試集
    X_train, y_train = X, y
    # 創建隨機森林模型
    rf = RandomForestRegressor(random_state=10)
    rf.fit(X_train, y_train)
    y_pred = rf.predict(X_train)
    # 計算測試集的 AUC 分數
    mse_score = mean_squared_error(y_train, y_pred)
    return mse_score
def predict_sp500_ver1(sp500_top_30, keep_cols):
    data_combined = None
    for stock_code in sp500_top_30:
        try:
            #取得股票資料
            start_date = datetime.datetime(2020, 1, 1)
            end_date = datetime.datetime(2022, 12, 31)
            data = get_stock_data(stock_code, start_date, end_date)
            if data_combined is None:
                data_combined = data
            else:
                data_combined = data_combined.append(data)
            print(f'股票代碼{stock_code}資料抓取成功')
        except Exception as e:
            print(f"股票代碼 {stock_code} 資料抓取失敗: {e}")
            continue
    if data_combined is not None:
        avg_mse_score_std = predict_stock_price_ver1(data_combined, keep_cols)
        print("預測完成")
        return avg_mse_score_std
    else:
        print("沒有可用的股票資料")
        return None
#算平均mse_ret分數
def predict_sp500_ver2(sp500_top_30, keep_cols):
    data_combined = None
    for stock_code in sp500_top_30:
        try:
            #取得股票資料
            start_date = datetime.datetime(2020, 1, 1)
            end_date = datetime.datetime(2022, 12, 31)
            data = get_stock_data(stock_code, start_date, end_date)
            if data_combined is None:
                data_combined = data
            else:
                data_combined = data_combined.append(data)
            print(f'股票代碼{stock_code}資料抓取成功')
        except Exception as e:
            print(f"股票代碼 {stock_code} 資料抓取失敗: {e}")
            continue
    if data_combined is not None:
        avg_mse_score_ret = predict_stock_price_ver2(data_combined, keep_cols)
        print("預測完成")
        return avg_mse_score_ret
    else:
        print("沒有可用的股票資料")
        return None
    
#寫基因演算法
import deap
from deap import algorithms, base, creator, tools
import random
random.seed(42)
#設定基因和適應度函數
creator.create('FitnessMin', base.Fitness, weights=(-1.0, -1.0))
creator.create('Individual', list, fitness=creator.FitnessMin)

#設定工具箱
toolbox = base.Toolbox()
toolbox.register('attr_bool', random.randint, 0, 1)
toolbox.register('individual', tools.initRepeat, creator.Individual, toolbox.attr_bool, n=180)
toolbox.register('population', tools.initRepeat, list, toolbox.individual)
toolbox.register('mate', tools.cxTwoPoint)
toolbox.register('mutate', tools.mutFlipBit, indpb=0.05)

#設定適應度函數
def evaluate_individual(individual):
    keep_cols = individual
    sp500_stock_code = ['AAPL', 'MSFT']
    #計算mse平均值
    avg_mse_score_std = predict_sp500_ver1(sp500_stock_code, keep_cols)
    avg_mse_score_ret = predict_sp500_ver2(sp500_stock_code, keep_cols)

    # 返回一個元組，包含兩個目標
    return (avg_mse_score_std, avg_mse_score_ret)

toolbox.register('evaluate', evaluate_individual)
toolbox.register('select', tools.selNSGA2)
#設定演化參數
population_size = 100
number_of_generations = 50
probability_of_crossover = 0.5
probability_of_mutation = 0.2

#初始化族群
population = toolbox.population(n=population_size)

#執行演化
for generation in range(number_of_generations):
    #交配變異
    offspring = algorithms.varAnd(population, toolbox, cxpb=probability_of_crossover, mutpb=probability_of_mutation)
    #計算每個個體的適應度
    fitnesses = map(toolbox.evaluate, offspring)
    #跟新適應度
    for ind, fit in zip(offspring, fitnesses):
        ind.fitness.values = fit
    #跟新種群
    population = toolbox.select(offspring, k=population_size)

#獲得最佳解
best_individual = tools.selBest(population, k=1)[0]
best_fitness = best_individual.fitness.values

print('Best individual:', best_individual)
print('Best fitness:', best_fitness)

#最佳化之後的特徵選擇
def predict_std(data, keep_cols):
    # 根据 keep_cols 的值删除需要删除的列
    # 用除了漲跌的欄位作為自變數，去預測股價漲跌
    df = copy.deepcopy(data)
    df.fillna(0, inplace=True)
    df['stddev_lagging'] = df['stddev'].shift(-1)
    df = df.dropna()
    scaler = StandardScaler()
    X = df.drop(['stddev_lagging'], axis=1)
    drop_cols = [col for i, col in enumerate(X.columns) if not keep_cols[i]]
    X = X.drop(drop_cols, axis=1)
    X = X.dropna()
    X = scaler.fit_transform(X)
    y = df['stddev_lagging'].values.reshape(-1,1)
    y = scaler.fit_transform(y)
    # 分割資料集為訓練集和測試集
    X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.3, shuffle=False)
    # 創建隨機森林模型
    rf = RandomForestRegressor(random_state=10, n_jobs=-1)
    rf.fit(X_train, y_train)
    y_pred = rf.predict(X_test)
    # 計算平均標準差（stddev）分數
    avg_stddev = np.mean(y_pred)
    return avg_stddev
def predict_ret(data, keep_cols):
    # 根据 keep_cols 的值删除需要删除的列
    # 用除了漲跌的欄位作為自變數，去預測股價漲跌
    df = copy.deepcopy(data)
    df.fillna(0, inplace=True)
    df['return_lagging'] = df['return'].shift(-1)
    df = df.dropna()
    scaler = StandardScaler()
    X = df.drop(['return_lagging'], axis=1)
    drop_cols = [col for i, col in enumerate(X.columns) if not keep_cols[i]]
    X = X.drop(drop_cols, axis=1)
    X = X.dropna()
    X = scaler.fit_transform(X)
    y = df['return_lagging'].values.reshape(-1,1)
    y = scaler.fit_transform(y)
    # 分割資料集為訓練集和測試集
    X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.3, shuffle=False)
    # 創建隨機森林模型
    rf = RandomForestRegressor(random_state=10, n_jobs=-1)
    rf.fit(X_train, y_train)
    y_pred = rf.predict(X_test)
    avg_return = np.mean(y_pred)
    return avg_return

def predict_sp500(sp500_top_250, best_individual):
    data_dict = {}
    for stock_code in sp500_top_250:
        try:
            #取得股票資料
            start_date = datetime.datetime(2018, 10, 31)
            end_date = datetime.datetime(2020, 12, 31)
            data = get_stock_data(stock_code, start_date, end_date)
            # 根據最佳化結果進行特徵選擇
            keep_cols = best_individual
            #進行預測
            avg_std_score = predict_std(data, keep_cols)
            avg_ret_score = predict_ret(data, keep_cols)
            
            print(f"股票代碼 {stock_code} 預測成功")
        except Exception as e:
            print(f"股票代碼 {stock_code} 預測失敗: {e}")
            continue

        #預測結果添加到字典中
        data_dict[stock_code] = {'隱含波動度':avg_std_score, '隱含報酬率':avg_ret_score}
    df = pd.DataFrame.from_dict(data_dict, orient="index")
    

    return df

def filter_stocks(df):
    volatility_quantiles = df['隱含波動度'].quantile([0.3, 0.7])
    df['波動度分類'] = pd.cut(df['隱含波動度'], bins=[-np.inf, volatility_quantiles.iloc[0],
                                            volatility_quantiles.iloc[1],
                                            np.inf],labels=['低波動', '中等波動', '高波動'])
    returns_quantiles = df['隱含報酬率'].quantile([0.3, 0.7])
    df['報酬率分類'] = pd.cut(df['隱含報酬率'], bins=[-np.inf, returns_quantiles.iloc[0], returns_quantiles.iloc[1], np.inf], labels=['低報酬', '中等報酬', '高報酬'])

    # 選取低波動且高報酬的股票代碼
    filtered_stocks = df[(df['波動度分類'] == '低波動') & (df['報酬率分類'] == '高報酬')].index.tolist()

    return filtered_stocks


#視覺化投組資料
def get_stock_data_only(stock_code, start_date, end_date):
    #取得股票資料
    df = pdr.get_data_yahoo([stock_code], start_date, end_date)
    df.columns = [col.lower() for col in df.columns]
    df['return'] = df['adj close'].pct_change()
    df = df.dropna()
    return df

def plot_portfolio_returns(filtered_stocks, start_date, end_date):
    portfolio_returns = []
    portfolio_returns_sp_500=[]
    stock_code_sp_500 = 'SPY'
    dates = pd.date_range(start=start_date, end=end_date, freq='B')
    for date in dates:
        # 取得當天的股票資料
        stock = []
        for stock_code in filtered_stocks:
            stock_data = get_stock_data_only(stock_code, date - timedelta(days=1), date + timedelta(days=1))
            if not stock_data.empty:
                stock.append(stock_data.iloc[0]['return'])
            else:
                stock.append(0)     
        avg_return = sum(stock) / len(stock)
        portfolio_returns.append(avg_return)

        stock_data_sp_500 = get_stock_data_only(stock_code_sp_500, date - timedelta(days=1), date + timedelta(days=1))
        if not stock_data_sp_500.empty:
            return_sp_500 = stock_data_sp_500.iloc[0]['return']
        else:
            return_sp_500 = 0
        portfolio_returns_sp_500.append(return_sp_500)
   
    df = pd.DataFrame({'profit':np.cumsum(portfolio_returns), 'sp 500 return':np.cumsum(portfolio_returns_sp_500)}, index=dates)
    df.plot(grid=True, figsize=(16, 6))
    plt.legend()
    plt.ylabel('Profit')
    plt.xlabel('date')
    plt.title('Portfolio Returns')
    plt.show()
    df = df.reset_index(drop=True)

    # 計算總利潤、股權、回撤百分比和回撤
    df['equity'] = df['profit']+1
    df['drawdown_percent'] = (df['equity'] / df['equity'].cummax())-1
    df['drawdown'] = df['equity'] - df['equity'].cummax()
    fig, ax = plt.subplots(figsize = (16,6))
    high_index = df[df['profit'].cummax() == df['profit']].index
    df['profit'].plot(label = 'Total Profit', ax = ax, c = 'r', grid=True)
    plt.fill_between(df['drawdown'].index, df['drawdown'], 0, facecolor = 'r', label = 'Drawdown', alpha = 0.5)
    plt.scatter(high_index, df['profit'].loc[high_index], c = '#02ff0f', label = 'High')
    plt.legend()
    plt.ylabel('Accumulated Profit Return')
    plt.xlabel('Time')
    plt.title('Portfolio Profit & Drawdown',fontsize  = 16)
    plt.show()
    profit = df['profit'].iloc[-1]
    mdd = abs(df['drawdown_percent'].min())
    calmarratio = profit/mdd
    print("portfolio return & risk ")
    print(f'return: ${np.round(profit, 4)*100}%')
    print(f'mdd: {np.round(mdd, 4)*100}%')
    print(f'calmar ratio: {np.round(calmarratio, 2)}')

    #計算持有大盤
    df_sp500 = pd.DataFrame({'profit':np.cumsum(portfolio_returns_sp_500)})
    df_sp500['equity'] = df_sp500['profit']+1
    df_sp500['drawdown_percent'] = (df_sp500['equity'] / df_sp500['equity'].cummax())-1
    df_sp500['drawdown'] = df_sp500['equity'] - df_sp500['equity'].cummax()
    fig, ax = plt.subplots(figsize = (16,6))
    high_index = df_sp500[df_sp500['profit'].cummax() == df_sp500['profit']].index
    df_sp500['profit'].plot(label = 'Total Profit', ax = ax, c = 'r', grid=True)
    plt.fill_between(df_sp500['drawdown'].index, df_sp500['drawdown'], 0, facecolor = 'r', label = 'Drawdown', alpha = 0.5)
    plt.scatter(high_index, df_sp500['profit'].loc[high_index], c = '#02ff0f', label = 'High')
    plt.legend()
    plt.ylabel('Accumulated Profit Return')
    plt.xlabel('Time')
    plt.title('SP500 Profit & Drawdown',fontsize  = 16) 
    plt.ylim(top=df['profit'].max())
    plt.show()
    profit = df_sp500['profit'].iloc[-1]
    mdd = abs(df_sp500['drawdown_percent'].min())
    calmarratio = profit/mdd
    print("sp 500 return & risk ")
    print(f'return: ${np.round(profit, 4)*100}%')
    print(f'mdd: {np.round(mdd, 4)*100}%')
    print(f'calmar ratio: {np.round(calmarratio, 2)}')
    
    return portfolio_returns, portfolio_returns_sp_500