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

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

CVXOPTは凸最適化問題を解くPythonのフリーのライブラリです．

今回は，ポートフォリオ最適化を例にして，CVXOPTで２次計画問題を実装してみます．
比較として，SciPyを用いた例は，[Scipyで２次計画問題を解く〜ポートフォリオ最適化の例〜](https://qiita.com/ryoshi81/items/b323a363e5442a15db6d)をご覧ください．

## 準備

### CVXOPTをインストールする

[CVXOPT Installation instructions](http://cvxopt.org/install/index.html)を参考にしてください．


```python
pip install cvxopt
```
### データの準備

今回用いるCVXOPT以外のライブラリをインポートします．


```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]
```

期待リターン(収益率)を計算します．

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

リスク(分散共分散行列)を計算します．

```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/51b61aa4-0a8e-aed7-6c17-327a1580314e.png)


## CVXOPTで二次計画問題を解く

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

![finpy_study_171025_qp.png](https://qiita-image-store.s3.amazonaws.com/0/193486/5f7b4a81-994e-a1af-6c28-b5f64b2444aa.png)

### cvxopt.solvers.qp()を使う


```python
from cvxopt import matrix
from cvxopt.solvers import qp
```

[CVXOPT Quadratic Programming](http://cvxopt.org/userguide/coneprog.html#quadratic-programming)を確認すると，qpの引数はP, q, G, h, A, bがあります．よって，ポートフォリオ最適化問題を以下の形式に整理して，行列を求める必要があります．その分，SciPyより手間がかかります．

```math

\begin{align}
\min \quad & \frac{1}{2} {\bf x}^T P {\bf x} + {\bf q}^T {\bf x} \\
{\rm subject \: to} \quad & G {\bf x} \leq {\bf h} \\
& A {\bf x} = {\bf b} \\
\end{align}

```

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

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


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

### 目的関数

目的関数はポートフォリオの分散なので，Pは分散共分散行列になります．cvxopt.solvers.qpの引数はcvxopt.matrix型にする必要があります．


```python
P = matrix(np.array(cov))
```

目的関数はポートフォリオの分散なので，qはすべての要素が0の列ベクトルです．なお，「0」とするとmatrixのtc='i'の整数となってしまうため，「0.0」としてtc='d'となるようにします．


```python
q = matrix([0.0] * len(mean))
```

### 制約条件

Aとbは投資比率の和が1となる等式制約に当たります．Aは要素がすべて1.0の行ベクトルになります．


```python
A = matrix(1.0, (1, len(mean)))
```

投資比率の和が1なので，bは1.0です．


```python
b = matrix(1.0)
```

不等式制約はポートフォリオの期待リターンが目標リターン以上となる制約と下限制約にあたります．$G$は以下のようになります．

```math
\begin{eqnarray}
  G = \left(
\begin{array}{ccc}
1 & \cdots & 0 \\
\vdots & \ddots & \vdots \\
0 & \cdots & 1 \\
-r_1 & \cdots & -r_n
\end{array}
  \right)
\end{eqnarray}
```

```python
G = matrix(np.c_[-np.identity(len(mean)), np.array(-mean)].T)
```

hは下限0と目標リターン$r_G$の列ベクトルになります．

$$h=(0,\cdots,0,-r_G)^T$$


```python
h = matrix(np.r_[np.zeros(len(mean)), np.array(-tret, ndmin=1)])
```

まとめると，次のようになります．


```python
tret = mean.mean()
P = matrix(np.array(cov))
q = matrix([0.0] * len(mean))
A = matrix(1.0, (1, len(mean)))
b = matrix(1.0)
G = matrix(np.c_[-np.identity(len(mean)), np.array(-mean)].T)
h = matrix(np.r_[np.zeros(len(mean)), np.array(-tret, ndmin=1)])
```

### 最適化

最適化を実行します．


```python
%%time
sol = qp(P=P, q=q, G=G, h=h, A=A, b=b)
```

         pcost       dcost       gap    pres   dres
     0:  1.9211e+00  1.3849e+00  5e+01  6e+00  2e+01
     1:  2.2276e+00  1.8570e+00  7e+00  9e-01  2e+00
     2:  2.9679e+00  2.7243e+00  2e+00  2e-01  4e-01
     3:  3.2479e+00  3.1472e+00  2e-01  1e-02  3e-02
     4:  3.2302e+00  3.2133e+00  2e-02  3e-04  8e-04
     5:  3.2202e+00  3.2188e+00  1e-03  9e-06  2e-05
     6:  3.2192e+00  3.2191e+00  1e-04  5e-08  1e-07
     7:  3.2191e+00  3.2191e+00  9e-06  4e-11  1e-10
     8:  3.2191e+00  3.2191e+00  3e-07  2e-13  4e-13
    Optimal solution found.
    CPU times: user 4.7 ms, sys: 2.26 ms, total: 6.96 ms
    Wall time: 5.28 ms


最適投資比率を確認します．


```python
list(sol['x'])
```

最適値(ポートフォリオの分散)を確認します．


```python
2 * sol['primal objective']
```
    6.438220501299366

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

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


```python
mean_high = mean.max().round(3)
mean_low = round(mean.quantile(0.25), 3)
trets = np.linspace(mean_low, mean_high, 50)
```

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


```python
%time
tvols = []
for tret in trets:
    h = matrix(np.r_[np.zeros(len(mean)), np.array(-tret, ndmin=1)])
    sol = qp(P=P, q=q, G=G, h=h, A=A, b=b)
    tvols.append(np.sqrt(2 * sol['primal objective']))
tvols = np.array(tvols)
```

    CPU times: user 2 µs, sys: 1 µs, total: 3 µs
    Wall time: 4.77 µ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/4c7b3b79-a79e-31b4-c629-a98030c3dd10.png)

### SciPyとの比較

今回の例では決定変数30個だけですが，CVXOPTのほうがかなり速いです．

| 計測方法 |SciPy  |CVXOPT  |
|---|---|---|
|%%time |  65.8ms|  6.16ms|
|%%timeit |  44.6ms|  4.49ms|

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


## 参考

- [CVXOPT User's Guide](http://cvxopt.org/userguide/index.html)
- [Scipyで２次計画問題を解く〜ポートフォリオ最適化の例〜](https://qiita.com/ryoshi81/items/b323a363e5442a15db6d)