理科系の勉強日記

Linux/Ubuntu/Mac/Emacs/Computer vision/Robotics

ポートフォリオ最適化でCVXOPTを学ぶ


リターンを確保しつつリスクを最小にするポートフォリオを構成するにはどの銘柄をどれだけ持てばいいのでしょう?

ポートフォリオに銘柄$i, (i = 1, \dots, N)$をどれだけの割合もつかということを$\boldsymbol{x} = [ x_1, \dots, x_N ]^{T}$、各銘柄のリターンの共分散行列を$C$とすると、ポートフォリオ全体の分散は$\boldsymbol{x}^{T}C\boldsymbol{x}$となります。各銘柄の期待リターンを$\boldsymbol{m} = [ m_1, \dots, m_N ]^{T}$とすると、ポートフォリオ全体のリターンは$\sum_{i=1}^{N}{m_i x_i}$ですね。ポートフォリオの分散を最小化しつつリターン$r$以上を求めるには、以下の条件で最適化を解くことになります。この最適化は非線形最適化のうち二次計画法で解ける形です。



\begin{aligned}
& \text{minimize } && \boldsymbol{x}^{T}C\boldsymbol{x} \\
& \text{subject to} && \sum_{i=1}^{N}{x_i} = 1 \\
& && \sum_{i=1}^{N}{m_i x_i} \geq r \\
& && 0 \leq x_i \leq 1 &&& (i = 1, \dots, N)
\end{aligned}


まずは米株のデータをQuandlで取得します。収益率は前日との変化率としました。

from datetime import datetime
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import seaborn as sns
import quandl

start = datetime(2010, 1, 1)
end = datetime(2018, 12, 31)
symbls = [ # DOW30種
"AAPL",
"AXP",
"BA",
"CAT",
"CSCO",
"CVX",
# 中略
"WMT",
"XOM"
]

# Quadlで株価取得
codes = ['WIKI/' + symbl for symbl in symbls]
data = quandl.get(codes, start_date=start, end_date=end)
types = [code + " - Adj. Close" for code in codes]
df = data[types]
df.columns = symbls
df_return = df.pct_change().dropna() # 収益率
N = len(df.columns) # 銘柄数

f:id:kenbell1988:20200209172954p:plainf:id:kenbell1988:20200209172951p:plainf:id:kenbell1988:20200209172948p:plain
APPLの株価と日々の収益率


ある銘柄の期待収益率は、ある期間での日々の収益率の平均値に日数を掛けたものとしました。各銘柄の収益率の共分散行列も求めます。こちらは日々の収益率の分散から計算しています。分散を最小化するポートフォリオには少なくとも15%程度の収益率を持たせることにしました。

# 期待(?)収益率 = 全区間の収益率の平均値に任意の期間を掛けたもの
duration = 252 # [days]
m = df_return.mean().values * duration

# 日々の収益率の共分散行列
cov = df_return.cov().values 

# 最適化するPFの収益率の最小値
r = 0.15 

f:id:kenbell1988:20200209173827p:plain
各銘柄の収益率の共分散行列

以上でポートフォリオの分散を最小にする銘柄選択のための数字が準備できました。これをCVXOPTで計算し、どの銘柄をどれだけの割合で持てばいいかを計算してみます。目的関数と制約条件をCVXOPTに読み込ませるため、多少の式変形を行います。CVXOPTが想定する2次計画問題は以下の形となっています。



\begin{aligned}
& \text{minimize } && \frac{1}{2}\boldsymbol{x}^{T}P\boldsymbol{x} + \boldsymbol{q} \boldsymbol{x}\\
& \text{subject to} && G\boldsymbol{x} \leq \boldsymbol{h} \\
& && A\boldsymbol{x} = \boldsymbol{b}
\end{aligned}

# 最適化ソルバへの入力を準備
g1 = np.eye(N)
g2 = -np.eye(N)
G = np.vstack((g1, g2))
G = np.vstack((G, -m))

# h = [1,1,...1, 0,0,...,0,-r]
h1 = np.ones(N)
h2 = np.zeros(N)
h = np.hstack((h1, h2))
h = np.hstack((h, -r))

最後にCVXOPTの最適化計算を行います。

import cvxopt
from cvxopt import matrix
cvxopt.solvers.options['show_progress'] = False

P = matrix(cov)
q = matrix(np.zeros(N))
A = matrix(np.ones(N).reshape(1,N))
b = matrix(1.0)
G = matrix(G)
h = matrix(h)

sol = cvxopt.solvers.qp(P,q, A=A, b=b,G=G, h=h)

print("Optimized Portfolio")
for i, e in enumerate(sol["x"]):
    if e > 10e-3:
        print(f"{df_return.columns[i]}\t: {e:.4f}[%]")

risk_e = 2 * sol["primal objective"]
x_e = np.array(list(sol["x"]))
r_e = m.dot(x_e)

print("-------------------")
print(f"PF variance \t: {risk_e:.8f}")
print(f"PF std \t\t: {np.sqrt(risk_e):.8f}")
print(f"PF return \t: {r_e:.8f}")

# ===> Result
# Optimized Portfolio
# BA	: 4.46[%]
# CVX	: 2.35[%]
# DOW	: 32.37[%]
# KO	: 11.50[%]
# MCD	: 4.10[%]
# MMM	: 30.72[%]
# NKE	: 2.94[%]
# PG	: 2.70[%]
# TRV	: 3.79[%]
# WMT	: 3.69[%]
# -------------------
# PF variance : 0.00001859
# PF std 		: 0.00431213
# PF return 	: 0.15000158

収益率を15%以上に保ちながら分散を最小にするには、以上のように銘柄を選択すればいいようです。この計算を信じる限りは。

。。。


問題は各銘柄の期待収益率、収益率の分散や共分散行列がどう求めるかでしょうね。
やってから検索してQiitaでよく似たことをやっている記事がありました。先に見なくてよかった。。。

pytorchハマったところ備忘録

Cuda assert fails: device-side assert triggered

BCELossにnanを入力していました。BCELossの入力は[0,1]なのでエラーが出ますね。BCEWithLogitsを使えば[0,1]以外も扱えますが、nanをだしている方をなんとかしました。nanが出ていた理由は学習率が高すぎたというものでした。schedulerでburn in期間を設定して、最初は小さなものから徐々に大きくすることにしました。

logをコピーしすぎてGPUメモリが足りない

for batch in dataloader:
    (中略)
    running_loss += loss # loss.item()でないと勾配情報をどんどんコピーしていまい、メモリリーク

Accumulated gradients

大きなバッチサイズを確保したいがGPUのメモリが足りないので、勾配を貯めることにしました。 BNはbatch1つ分のままになってしまうようですので注意が必要です。

accumulation_steps = 4 # batch_size * accumulation_stepsだけ蓄積
for itr, batch in enumerate(dataloader):
    imgs, targets = batch
    loss = model(imgs, targets)
    loss = loss / accumulation_steps
    loss.backward()
    if (itr + 1) % accumulation_steps ==0:
        optimizer.step()
        optimizer.zero_grad()

````

B spline曲線と曲面

B spline曲線について勉強したかったんですけどずっと放置していたので、ついに実装して遊んでみました。 basic functionとknot vectorが曲線にどう関係してくるのか実装するまでわかっていませんでした。


S(t) = \sum_{i=0}^n{N_{i,p}(t)P_i}

Basic function


N_{i,0}(u) =  \left\{ \begin{array}{ll}
    1 & (u_i \leq u < u_{i+1}) \\
    0 & (otherwise)
  \end{array} \right. \\
N_{i,p}(u) = \frac{u-u_i}{u_{i+p}-u_i}N_{i, p-1}(u) + \frac{u_{i+p+1}-u}{u_{i+p+1}-u_{i+1}}N_{i+1,p-1}(u)

f:id:kenbell1988:20190304093622p:plainf:id:kenbell1988:20190304093638p:plain
B spline curve

f:id:kenbell1988:20190304093646p:plainf:id:kenbell1988:20190304093700p:plain
Knot vector, basic function

ついでにB spline曲面も実装しました。u方向の制御点m個とv方向の制御点n個をつかって、control netなる m \times nの行列 Pをつかって以下のように表現できます[1]。実装はgithubのjupyter notebook参照。


S(u, v) = \sum_{i=0}^n\sum_{j=0}^m{N_{i,p}(u)N_{j,q}(v)P_{i,j}}

f:id:kenbell1988:20190304093926p:plain
B spline surface

github.com

参考

[1]https://graphics.stanford.edu/courses/cs348a-09-fall/Handouts/surfaces8.pdf

ガウス過程回帰を実装して遊ぶ

ガウス過程回帰を実装して遊んだ

ガウス過程の本が発売されるということで復習として簡単にメモしようとおもったけど、はてなTeX記法にうまく馴染めなかった1ので 数式は載せずに結果だけ簡単に貼り付けて日記としておきます。

結果

ガウス過程回帰、ハイパーパラメタをうまく決めないとゴミみたな結果が得られるようです。 適当な回数イテレーションを回してハイパーパラメタを決めましょう。やり方はいろいろあると思いますが、MCMCで決めるのがいいかなと考えています。 勾配法はうまくいったためしがないので。。。

f:id:kenbell1988:20190211143121p:plainf:id:kenbell1988:20190211143126p:plainf:id:kenbell1988:20190211143205p:plain
ガウス過程回帰の結果.(左)入力、(中)適当なパラメタ、(右)よさそうなパラメタ


  1. インラインで数式がうまく反映されない。中身を変えずにブロックにすると正しく表示される。やる気失うタイプのトラブル。

続きを読む

Jupyter notebookのthemeを変えたら出力セルが見切れる件

最高なjupyter notebookをさらに最高にするためのjupyter notebook themesをインストールしたが、
どうも表示の調子がよくない。端が切れている。
f:id:kenbell1988:20180929172403p:plain
OS:Mac OS High Sierra
ブラウザ:Google Chrome v64
(Safariでも同様)

このままでは辛いので、もとのcssを直接編集して対応することにした。
~/.jupyter/custom/custom.cssにあるだろうと予想をつけて、
適当にoutputで検索しているとそれらしいものを発見した。
山勘でpaddingを増やしたところうまくいってしまった。

div.output_subarea {
 overflow-x: auto;
 padding: 1.8em !important;  /* もともとコレ padding: 0.8em !important;*/
 -webkit-box-flex: 1;
 -moz-box-flex: 1;
 box-flex: 1;
 flex: 1;
}

f:id:kenbell1988:20180929173018p:plain
1.8が妥当かはまったく検討なし。

以上。

Juliaで主成分分析

周りの人と読んでいる本(ベイズ推論による機械学習)の影響でJuliaをさわってみました。
まだあまりむずかしいことができないので、操作に慣れるために主成分分析をしてみました。
(主成分分析を一発で実行する関数もあるかもしれませんが、練習なので教科書通りやってみます)

using PyCall, PyPlot
using Distributions
# データの作成
n = 300
μ  = [1.0, 3.0]
Σ  = [10 3
        3 2]
X = rand(MvNormal( μ, Σ), n);

MvNormal関数で2次元のガウス分布に従うデータを作成しました。 全角の文字も変数にできるのでギリシャ文字を変数にできて楽しい。

# データから平均を計算して、データから引く
mu = [mean(X[1, :])
            mean(X[2, :])]

X_m = copy(X)
X_m[1, :] = X[1, :] - mu[1];
X_m[2, :] = X[2, :] - mu[2];

2x300のデータから2x1の形のベクトルをうまく引き算する方法がわからなかったのでダサい方法で。。

# 共分散行列の計算
Σ = 1.0 / n *  X_m * X_m'

2×2 Array{Float64,2}: 10.5301 3.24329 3.24329 2.1232

# 固有値分解
λ, S = eig(Σ)
println(λ)
println(S)

# 線形変換
y = S' * X_m;

[1.01741, 11.6359] [0.322704 -0.9465; -0.9465 -0.322704]

# 可視化用の第一成分、第二成分
p_ax1 = [mu[1] mu[1] + S[1,1]*4
                mu[2] mu[2] + S[2,1] *4]
p_ax2 = [mu[1] mu[1] + S[1,2]*4
                mu[2] mu[2] + S[2,2]*4 ]

2×2 Array{Float64,2}: 1.15528 -2.63072 3.10447 1.81366

figure("result")
clf()
subplot(121)
title("PCA result")
plot(X[1,:], X[2, :], "+")
plot(p_ax1[1, :], p_ax1[2, :])
plot(p_ax2[1, :], p_ax2[2, :])

axis("equal")

subplot(122)
title("Linear Transformed")
plot(y[1,:], y[2, :], "+")
axis("equal")

show()

プロットもPyPlotを使えば(pythonになじみがあれば)あまり迷うことなくできました。 次は画像を使って見ようかな、データ用意するの面倒だけど。

f:id:kenbell1988:20180512142948p:plain

セルオートマトンによる渋滞シミュレーション

はじめに

年末年始に渋滞学という本を読んだ。
車の渋滞だけではなく、緊急時の避難や蟻の行列、通信についても書かれていて勉強になった。

渋滞学 (新潮選書)

渋滞学 (新潮選書)

渋滞とセルオートマトン

交通をモデル化して解析することを交通流解析と呼ぶ。交通流モデルは連続モデルと離散モデルに分けられ、連続モデルでは最適速度過程がよく用いられる。離散モデルでは、セルオートマトンがよく用いられており、「渋滞学」でもこれについて説明されていた。

続きを読む