Juliaでポートフォリオ最適化

Juliaの凸最適化のパッケージ[Convex](https://jump.dev/Convex.jl/stable/)とソルバーのラッパーのパッケージ[SCS](https://github.com/jump-dev/SCS.jl)でポートフォリオ最適化を試してみます。

ポートフォリオ最適化の説明やPythonとの比較は以前書いた記事を参考に。
- [Scipyで２次計画問題を解く〜ポートフォリオ最適化の例〜](https://qiita.com/ryoshi81/items/b323a363e5442a15db6d)
- [CVXOPTで２次計画問題を解く〜ポートフォリオ最適化の例〜](https://qiita.com/ryoshi81/items/8b0c6add3e367f94c828)
- [PyPortfolioOptでポートフォリオ最適化](https://qiita.com/ryoshi81/items/c2445a5705b495002c1d)

この記事は以下を参考にしています。
https://jump.dev/Convex.jl/v0.12/examples/portfolio_optimization/portfolio_optimization2/

# 準備

## パッケージをインストールする

パッケージの追加モードで以下のコマンドを実行します。

```
add Convex
add SCS
```

## リスクとリターン

パッケージをインポートします。

```julia
using Convex, SCS
using Plots
```

リスクと期待リターンは、今回は実際のデータからではなく、仮の値を用います。

期待リターンμを

```julia
μ = [11.5; 9.5; 6] / 100
```

とします。なお、ギリシャ文字μは「\mu」と入力後「Tab」で変換できます。

次に、分散共分散行列Σを

```julia
Σ = [166 34 58
 34 64 4
 58 4 100]　/　100^2
```

とします。投資対象資産数は

```julia
n = length(μ)
```

です。

# ポートフォリオ最適化

## 二次計画問題を解く

目的関数

$$\lambda \times risk − (1 − \lambda ) \times return$$

を最小化する問題を解きます。
効率的フロンティアを作成するため、λを複数パターン用意します。

```julia
N = 101
λ_vals = range(0.01, stop=0.99, length=N)
```

ConvexのVariableで決定変数を作成します。

```julia
w = Variable(n)
```

ポートフォリオのリターンは

```julia
ret  = dot(w, μ)
```

です。ポートフォリオのリスクは

```julia
risk = quadform(w, Σ)
```

です。 数理計画問題を、minimizeの1つ目の引数に目的関数、2つ目の引数に制約条件を入れて、定義します。

```julia
λ = 0.5
problem = minimize(λ * risk - (1 - λ) * ret, sum(w) == 1)
```

以下のように返ります。

```
minimize
└─ + (convex; real)
   ├─ * (convex; positive)
   │  ├─ 0.5
   │  └─ * (convex; positive)
   │     ├─ …
   │     └─ …
   └─ - (affine; real)
      └─ * (affine; real)
         ├─ …
         └─ …
subject to
└─ == constraint (affine)
   ├─ sum (affine; real)
   │  └─ 3-element real variable (id: 144…983)
   └─ 1

status: `solve!` not called yet
```

ソルバーにSCS.Optimizerを指定して実行します。

```julia
solve!(problem, SCS.Optimizer)
```

実行結果が表示されます。

```
----------------------------------------------------------------------------
	SCS v2.1.4 - Splitting Conic Solver
	(c) Brendan O'Donoghue, Stanford University, 2012
----------------------------------------------------------------------------
Lin-sys: sparse-indirect, nnz in A = 21, CG tol ~ 1/iter^(2.00)
eps = 1.00e-05, alpha = 1.50, max_iters = 5000, normalize = 1, scale = 1.00
acceleration_lookback = 10, rho_x = 1.00e-03
Variables n = 6, constraints m = 10
Cones:	primal zero / dual free vars: 2
	linear vars: 1
	soc vars: 7, soc blks: 2
Setup time: 9.30e-05s
----------------------------------------------------------------------------
 Iter | pri res | dua res | rel gap | pri obj | dua obj | kap/tau | time (s)
----------------------------------------------------------------------------
     0| 1.83e+19  1.94e+19  1.00e+00 -3.24e+19  2.99e+18  3.81e+19  3.30e-05 
    79| 3.50e-11  5.84e-11  3.16e-11 -6.43e-02 -6.43e-02  2.26e-16  4.48e-04 
----------------------------------------------------------------------------
Status: Solved
Timing: Solve time: 4.49e-04s
	Lin-sys: avg # CG iterations: 3.30, avg solve time: 6.40e-07s
	Cones: avg projection time: 7.48e-08s
	Acceleration: avg step time: 3.97e-06s
----------------------------------------------------------------------------
Error metrics:
dist(s, K) = 1.1102e-16, dist(y, K*) = 2.7756e-17, s'y/|s||y| = -3.3910e-18
primal res: |Ax + s - b|_2 / (1 + |b|_2) = 3.5044e-11
dual res:   |A'y + c|_2 / (1 + |c|_2) = 5.8361e-11
rel gap:    |c'x + b'y| / (1 + |c'x| + |b'y|) = 3.1561e-11
----------------------------------------------------------------------------
c'x = -0.0643, -b'y = -0.0643
============================================================================
```

最適解は

```julia
evaluate(w)
```

```
3-element Vector{Float64}:
  1.6436100090749952
  0.978590254896493
 -1.6222002640539048
```

で、そのときの期待リターンは

```julia
evaluate(ret)
```

```
0.18464920941555704
```

リスクは

```julia
evaluate(risk)
```

```
0.05602697621995912
```

です。

## 効率的フロンティア

効率的フロンティアをプロットするため、リスクとリターンの組み合わせを格納する箱を用意します。

```julia
MeanVarA = zeros(N, 2)
```

λを複数パターンで実行します。

```julia
for i = 1:N
    λ = λ_vals[i]
    problem = minimize(λ * risk - (1 - λ) * ret, sum(w) == 1)
    solve!(problem, SCS.Optimizer)
    MeanVarA[i,:]= [evaluate(ret),evaluate(risk)[1]]
end
```

決定変数に上下限制約を付ける場合です。

```julia
w_lower = 0
w_upper = 1
λ = 0.5
```

制約条件は第2引数以降で並べて書きます。

```julia
problem = minimize(λ * risk - (1 - λ) * ret,
    sum(w) == 1,
    w_lower <= w,
    w <= w_upper)
```

```julia
solve!(problem, SCS.Optimizer)
```

最適解は

```julia
evaluate(w)
```

効率的フロンティアを描くため、複数のλのパターンで実行します。

```julia
MeanVarB = zeros(N, 2)
for i = 1:N
    λ = λ_vals[i]
    problem = minimize(λ * risk - (1 - λ) * ret,
        sum(w) == 1,
        w_lower <= w,
        w <= w_upper)
    solve!(problem, SCS.Optimizer)
    MeanVarB[i,:]= [evaluate(ret), evaluate(risk)[1]]
end
```

効率的フロンティアをプロットします。

```julia
plot(sqrt.(MeanVarA[:,2]),
    MeanVarA[:,1],
    xlim = (0,0.25),
    ylim = (0,0.15),
    title = "Markowitz Efficient Frontier",
    xlabel = "Standard deviation",
    ylabel = "Expected return",
    label = "no bounds on w")
plot!(sqrt.(MeanVarB[:,2]),
    MeanVarB[:,1],
    label = "with 0<w<1")
scatter!(sqrt.(diag(Σ)),μ,color=:red,label = "assets")
```

![スクリーンショット 2021-10-30 11.22.35.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/193486/f968986b-9257-b86a-bb3b-ba44380bf370.png)

期待リターンが目標リターン以上となる制約条件を設定する場合のコードです。

まず目標リターンを複数パターン設定します。

```julia
retG_vals = range(minimum(μ),stop=maximum(μ),length=N)
```

　実行します。

```julia
MeanVarC = zeros(N, 2)
for i = 1:N
    retG = retG_vals[i]
    problem = minimize(risk,
        sum(w) == 1,
        retG <= ret,
        w_lower <= w,
        w <= w_upper)
    solve!(problem, SCS.Optimizer)
    MeanVarC[i,:]= [evaluate(ret), evaluate(risk)[1]]
end
```

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

```julia
plot(sqrt.(MeanVarC[:,2]),
    MeanVarC[:,1],
    xlim = (0,0.25),
    ylim = (0,0.15),
    title = "Markowitz Efficient Frontier",
    xlabel = "Standard deviation",
    ylabel = "Expected return",
    label = "Expected return >= Goal return")
scatter!(sqrt.(diag(Σ)),μ,color=:red,label = "assets")
```

![スクリーンショット 2021-10-30 11.24.58.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/193486/910e4323-1423-59a0-bb78-39068fce6b43.png)

# おわり

PythonのScipyやCVXOPTと比べても、Juliaだと数式のように直感的に書けます。