Scipyで２次計画問題を解く〜ポートフォリオ最適化の例〜

# SciPyで２次計画問題を解く〜ポートフォリオ最適化の例〜

２次計画問題はSciPyのオプティマイザーのメソッド'SLSQP'で解くことができます．

今回は，ポートフォリオ最適化を例にして，Scipyで２次計画問題を実装し，効率的フロンティアを描いてみます．

[SciPyの公式ページ](https://docs.scipy.org/doc/scipy-0.19.1/reference/generated/scipy.optimize.minimize.html)

## 準備

はじめに今回用いるライブラリをインポートしておきます．


```python
import numpy as np
import pandas as pd
import pandas_datareader.data as pdr
import matplotlib.pyplot as plt
%matplotlib inline
```

今回の例では，Kenneth R. French教授のData Libraryの30業種ポートフォリオのデータを用います．

[Kenneth R. Frenchデータライブラリ](http://mba.tuck.dartmouth.edu/pages/faculty/ken.french/data_library.html)

```python
name = '30_Industry_Portfolios'
df_dict = pdr.DataReader(name=name, data_source='famafrench')
df = df_dict[0]
df.head(10)
```
<div>
<table border="1" class="dataframe">
  <thead>
    <tr style="text-align: right;">
      <th></th>
      <th>Food</th>
      <th>Beer</th>
      <th>Smoke</th>
      <th>Games</th>
      <th>Books</th>
      <th>Hshld</th>
      <th>Clths</th>
      <th>Hlth</th>
      <th>Chems</th>
      <th>Txtls</th>
      <th>...</th>
      <th>Telcm</th>
      <th>Servs</th>
      <th>BusEq</th>
      <th>Paper</th>
      <th>Trans</th>
      <th>Whlsl</th>
      <th>Rtail</th>
      <th>Meals</th>
      <th>Fin</th>
      <th>Other</th>
    </tr>
    <tr>
      <th>Date</th>
      <th></th>
      <th></th>
      <th></th>
      <th></th>
      <th></th>
      <th></th>
      <th></th>
      <th></th>
      <th></th>
      <th></th>
      <th></th>
      <th></th>
      <th></th>
      <th></th>
      <th></th>
      <th></th>
      <th></th>
      <th></th>
      <th></th>
      <th></th>
      <th></th>
    </tr>
  </thead>
  <tbody>
    <tr>
      <th>2010-01</th>
      <td>-0.92</td>
      <td>-3.70</td>
      <td>-3.15</td>
      <td>1.12</td>
      <td>1.23</td>
      <td>0.85</td>
      <td>-3.11</td>
      <td>0.00</td>
      <td>-4.04</td>
      <td>-8.66</td>
      <td>...</td>
      <td>-6.84</td>
      <td>-7.10</td>
      <td>-7.88</td>
      <td>-5.49</td>
      <td>-4.04</td>
      <td>-2.55</td>
      <td>-1.92</td>
      <td>-0.98</td>
      <td>-1.18</td>
      <td>2.24</td>
    </tr>
    <tr>
      <th>2010-02</th>
      <td>2.92</td>
      <td>0.25</td>
      <td>5.00</td>
      <td>6.02</td>
      <td>-1.78</td>
      <td>3.69</td>
      <td>6.07</td>
      <td>0.37</td>
      <td>5.46</td>
      <td>15.68</td>
      <td>...</td>
      <td>2.92</td>
      <td>2.84</td>
      <td>6.30</td>
      <td>3.05</td>
      <td>7.25</td>
      <td>5.89</td>
      <td>4.32</td>
      <td>4.50</td>
      <td>2.71</td>
      <td>1.87</td>
    </tr>
    <tr>
      <th>2010-03</th>
      <td>4.44</td>
      <td>6.11</td>
      <td>6.36</td>
      <td>12.06</td>
      <td>5.53</td>
      <td>2.13</td>
      <td>11.39</td>
      <td>3.61</td>
      <td>7.34</td>
      <td>8.42</td>
      <td>...</td>
      <td>7.88</td>
      <td>5.22</td>
      <td>8.13</td>
      <td>5.86</td>
      <td>8.11</td>
      <td>5.99</td>
      <td>6.14</td>
      <td>8.27</td>
      <td>8.14</td>
      <td>9.07</td>
    </tr>
    <tr>
      <th>2010-04</th>
      <td>-1.46</td>
      <td>-1.76</td>
      <td>-3.03</td>
      <td>10.96</td>
      <td>2.32</td>
      <td>1.21</td>
      <td>4.23</td>
      <td>-2.22</td>
      <td>3.08</td>
      <td>11.43</td>
      <td>...</td>
      <td>3.82</td>
      <td>0.46</td>
      <td>3.89</td>
      <td>4.21</td>
      <td>3.82</td>
      <td>3.20</td>
      <td>1.85</td>
      <td>7.31</td>
      <td>0.95</td>
      <td>3.86</td>
    </tr>
    <tr>
      <th>2010-05</th>
      <td>-5.37</td>
      <td>-3.88</td>
      <td>-7.71</td>
      <td>-7.05</td>
      <td>-10.12</td>
      <td>-4.35</td>
      <td>-5.67</td>
      <td>-8.01</td>
      <td>-9.99</td>
      <td>-9.87</td>
      <td>...</td>
      <td>-5.91</td>
      <td>-8.31</td>
      <td>-7.07</td>
      <td>-8.89</td>
      <td>-5.98</td>
      <td>-4.01</td>
      <td>-5.83</td>
      <td>-4.76</td>
      <td>-9.12</td>
      <td>-10.05</td>
    </tr>
    <tr>
      <th>2010-06</th>
      <td>-1.97</td>
      <td>-1.63</td>
      <td>3.54</td>
      <td>-8.80</td>
      <td>-8.70</td>
      <td>-3.24</td>
      <td>-9.31</td>
      <td>-1.63</td>
      <td>-7.79</td>
      <td>-14.55</td>
      <td>...</td>
      <td>-4.33</td>
      <td>-6.78</td>
      <td>-5.99</td>
      <td>-2.30</td>
      <td>-7.66</td>
      <td>-5.20</td>
      <td>-10.56</td>
      <td>-5.21</td>
      <td>-7.31</td>
      <td>-7.58</td>
    </tr>
    <tr>
      <th>2010-07</th>
      <td>5.04</td>
      <td>8.37</td>
      <td>11.12</td>
      <td>7.85</td>
      <td>5.29</td>
      <td>4.22</td>
      <td>8.34</td>
      <td>2.15</td>
      <td>14.24</td>
      <td>6.81</td>
      <td>...</td>
      <td>9.74</td>
      <td>8.23</td>
      <td>6.82</td>
      <td>7.64</td>
      <td>10.84</td>
      <td>3.20</td>
      <td>4.65</td>
      <td>6.90</td>
      <td>6.62</td>
      <td>8.03</td>
    </tr>
    <tr>
      <th>2010-08</th>
      <td>-0.46</td>
      <td>-0.01</td>
      <td>-0.05</td>
      <td>-1.32</td>
      <td>-10.04</td>
      <td>-3.96</td>
      <td>-7.41</td>
      <td>-1.66</td>
      <td>-0.11</td>
      <td>-5.50</td>
      <td>...</td>
      <td>-2.23</td>
      <td>-5.26</td>
      <td>-8.21</td>
      <td>-6.06</td>
      <td>-4.65</td>
      <td>-7.55</td>
      <td>-4.18</td>
      <td>0.46</td>
      <td>-8.38</td>
      <td>-7.12</td>
    </tr>
    <tr>
      <th>2010-09</th>
      <td>1.57</td>
      <td>4.95</td>
      <td>10.04</td>
      <td>14.16</td>
      <td>12.91</td>
      <td>3.98</td>
      <td>17.54</td>
      <td>9.00</td>
      <td>9.96</td>
      <td>13.39</td>
      <td>...</td>
      <td>7.87</td>
      <td>11.62</td>
      <td>13.48</td>
      <td>9.16</td>
      <td>9.87</td>
      <td>9.90</td>
      <td>13.59</td>
      <td>7.06</td>
      <td>7.46</td>
      <td>10.01</td>
    </tr>
    <tr>
      <th>2010-10</th>
      <td>5.86</td>
      <td>1.93</td>
      <td>5.27</td>
      <td>12.57</td>
      <td>7.13</td>
      <td>4.51</td>
      <td>3.34</td>
      <td>2.00</td>
      <td>8.18</td>
      <td>3.28</td>
      <td>...</td>
      <td>4.80</td>
      <td>8.01</td>
      <td>4.14</td>
      <td>3.49</td>
      <td>5.35</td>
      <td>4.95</td>
      <td>1.32</td>
      <td>6.04</td>
      <td>1.97</td>
      <td>-0.11</td>
    </tr>
  </tbody>
</table>
<p>10 rows × 30 columns</p>
</div>

### 期待リターン(収益率)を計算する

![finpy_study_171025_リターン.png](https://qiita-image-store.s3.amazonaws.com/0/193486/85b194e1-6e87-b548-1f58-bb48cc2e7dc3.png)

```python
mean = df.mean()
```

平均は１次元の形式（pandas.Seriesや１次元のnumpy.ndarray）にしておきます．

### リスクを計算する

![finpy_study_171025_リスク.png](https://qiita-image-store.s3.amazonaws.com/0/193486/60c49045-bed8-72b9-5fb2-98056c89c12a.png)

```python
cov = df.cov()
std = np.sqrt(np.diag(cov))
```

### プロットする

各資産を標準偏差と期待リターンでプロットし，シャープレシオ(リスクフリーレートは考慮しない)で色付けします．


```python
plt.figure(figsize=(10, 6))
plt.scatter(std, mean, c=mean/std, marker='o')
plt.grid(True)
plt.xlabel('expected volatility')
plt.ylabel('expected return')
```

![output_16_1.png](https://qiita-image-store.s3.amazonaws.com/0/193486/0f1c7414-33f1-c7d0-cd32-7136e1e8bd66.png)


これらの資産への最適な投資比率を決定するのが目的です．

## SciPy Optimization で二次計画問題を解く

### H.M.Markowitzの平均・分散モデル

H.M.Markowitz(1952)は“Portfolio Selection”で，ポートフォリオの期待リターンが目標リターン以上になるという制約条件のもと，ポートフォリオの分散を最小化する平均・分散モデルを提案しています．

![finpy_study_171025_qp.png](https://qiita-image-store.s3.amazonaws.com/0/193486/730e0255-5d01-61e3-632e-7fc03a54a62f.png)

### scipy.optimize.minimize() を使う


```python
import scipy.optimize as sco
```

scipy.optimize.minimize()の主な引数は以下のものがあります．



|parameters|内容|type|
|:--|:--|:--|
|fun|目的関数|関数|
|x0|初期解|ndarray などシーケンス|
|method|ソルバー|'SLSQP':逐次最小二乗法|
|bounds|上下限制約|tupleのシーケンス|
|constraints|制約条件|dict または dictのシーケンス|


一つずつ説明していきます．

### 目的関数：最小化する関数を定義する

平均・分散モデルではポートフォリオの分散を目的関数とします．
投資比率を引数として，以下の関数を定義します．


```python
def min_func_var(weights):
    return np.dot(weights.T, np.dot(cov, weights))
```

### 初期解

初期解を設定します．今回は資産数$n$とし，各資産への投資比率が$\frac{1}{n}$となるように資産数分の初期解をリストとして生成します．


```python
x0 = [1. / len(mean)] * len(mean)
```

### 目標リターン：自分で決める

目標リターンを設定します．
目標リターンは，投資家が恣意的に決めますが，今回は個々の資産の期待リターンの水準を参考に決めます．


```python
tret = mean.mean()  # e.g.: tret = 0.5
```

### 制約条件

制約条件は等式制約と不等式制約があります．：
- 'type'に'eq'を指定すると等式制約で，'fun'に左辺=0となる関数を定義します．
- 'type'に'ineq'を指定すると不等式制約で，'fun'は左辺>=0となる関数を定義します
<br>

平均・分散モデルでは等式制約として投資比率の和が１となる制約，不等式制約としてポートフォリオのリターンが目標リターン以上となる制約を加えます．


```python
cons = [{'type': 'eq', 'fun': lambda x: np.sum(x) - 1}, 
        {'type': 'ineq', 'fun': lambda x: np.sum(mean * x) - tret}]
```

### 上下限制約

上下限制約：各決定変数の(最小，最大)の組を決定変数の数分並べたシーケンスとして定義します．

通常，空売り禁止制約として下限０を設定します．


```python
bnds = [(0, None)] * len(mean)
```

### 最適化

最適化を実行します．


```python
%%time
opts = sco.minimize(fun=min_func_var, x0=x0, method='SLSQP', bounds=bnds, constraints=cons)
```

    CPU times: user 64.7 ms, sys: 1.99 ms, total: 66.7 ms
    Wall time: 65.8 ms

最適化の結果を確認します．

```python
opts
```


         fun: 6.4382200570155153
         jac: array([ 12.86967444,  12.88338792,  17.07558435,  16.22200376,
            16.85564399,  12.86223888,  12.91390479,  12.88969362,
            17.00302589,  18.76555657,  17.04112363,  18.64207184,
            17.00152731,  16.5060941 ,  16.78500116,  14.17968714,
            13.70201218,  16.07571757,  14.00374794,  12.86757839,
            14.81756985,  14.18152452,  15.26657546,  16.06373262,
            14.38778698,  13.2118395 ,  13.79812872,  12.89780128,
            13.23550874,  13.95237762,   0.        ])
     message: 'Optimization terminated successfully.'
        nfev: 293
         nit: 9
        njev: 9
      status: 0
     success: True
           x: array([  8.09082690e-02,   1.81076080e-01,   0.00000000e+00,
             0.00000000e+00,   0.00000000e+00,   2.23625879e-01,
             5.29605656e-14,   4.31335643e-02,   0.00000000e+00,
             0.00000000e+00,   0.00000000e+00,   0.00000000e+00,
             0.00000000e+00,   0.00000000e+00,   0.00000000e+00,
             1.87447714e-14,   6.68065863e-14,   4.68453401e-14,
             3.45044089e-14,   2.70424480e-01,   0.00000000e+00,
             1.36019668e-14,   0.00000000e+00,   0.00000000e+00,
             1.52070197e-14,   4.56967363e-14,   1.35369146e-14,
             2.00831727e-01,   5.60252799e-14,   1.37926779e-14])

最適解と最適値をそれぞれ取り出します．

|Returns|内容|
|:--|:--|
|x|最適解|
|fun|最適値|
|success|最適化に成功したか|

```python
opts['x']
```

    array([  8.09082690e-02,   1.81076080e-01,   0.00000000e+00,
             0.00000000e+00,   0.00000000e+00,   2.23625879e-01,
             5.29605656e-14,   4.31335643e-02,   0.00000000e+00,
             0.00000000e+00,   0.00000000e+00,   0.00000000e+00,
             0.00000000e+00,   0.00000000e+00,   0.00000000e+00,
             1.87447714e-14,   6.68065863e-14,   4.68453401e-14,
             3.45044089e-14,   2.70424480e-01,   0.00000000e+00,
             1.36019668e-14,   0.00000000e+00,   0.00000000e+00,
             1.52070197e-14,   4.56967363e-14,   1.35369146e-14,
             2.00831727e-01,   5.60252799e-14,   1.37926779e-14])


```python
opts['fun']
```

    6.4382200570155153


```python
opts['success']
```

    True

### 効率的フロンティアを描く

![finpy_study_171025_frontier.png](https://qiita-image-store.s3.amazonaws.com/0/193486/3328bb2b-68c8-2bf1-0751-0fe2611e4a00.png)

目標リターンを複数パターン作ります．


```python
mean_high = mean.max().round(3)
mean_low = round(mean.quantile(0.25), 3)
```


```python
trets = np.linspace(mean_low, mean_high, 50)
```

各目標リターンについて，最適化を実行します．


```python
%time
tvols = []
for tret in trets:
    cons = [{'type': 'eq', 'fun': lambda x: np.sum(x) - 1}, 
            {'type': 'ineq', 'fun': lambda x: np.sum(mean * x) - tret}]
    res = sco.minimize(fun=min_func_var, x0=x0, method='SLSQP', bounds=bnds, constraints=cons)
    tvols.append(np.sqrt(res['fun']))
tvols = np.array(tvols)
```

    CPU times: user 2 µs, sys: 1e+03 ns, total: 3 µs
    Wall time: 6.2 µs


効率的フロンティアを描きます．


```python
plt.figure(figsize=(10, 6))
plt.scatter(tvols, trets, c=trets/tvols, marker='x')
plt.grid(True)
plt.xlabel('expected volatility')
plt.ylabel('expected return')
```

![output_60_1.png](https://qiita-image-store.s3.amazonaws.com/0/193486/9bfac5c0-aaff-3c4c-f02e-91c8617db25f.png)

ハイリターンを狙う場合は右上のリスク・リターンのプロットを実現する最適投資比率になるように資産配分し，ローリスクを狙う場合には左下のプロットを実現する最適投資比率になるように資産配分します．

## 免責事項
- 本内容はポートフォリオ理論や数理計画問題の一般論であり，投資を勧めるものでは一切ありません．
- 本投稿で計算した数値は仮定のものであり，その数値に一切責任は負いません．
- 本内容は個人の見解であり，所属組織とは一切関係ありません．

## 参考

- [SciPyの公式ページ](https://docs.scipy.org/doc/scipy-0.19.1/reference/generated/scipy.optimize.minimize.html)
- [Kenneth R. Frenchデータライブラリ](http://mba.tuck.dartmouth.edu/pages/faculty/ken.french/data_library.html)
- [Python-Highchartsのドリルダウンでポートフォリオを可視化する](https://qiita.com/ryoshi81/items/ea2c6113fe312fb1b52f)
- [CVXOPTで２次計画問題を解く〜ポートフォリオ最適化の例〜](https://qiita.com/ryoshi81/items/8b0c6add3e367f94c828)