粒子フィルタ使って実装したい (2)
はじめに
今回は、先生に見せるように書いたJupyter Notebookをブログに書き起こしました。
少しメモっぽい書き方になっていますが、ご了承ください。
Markowitz Portfolioの最適化
資産
4つの資産を持つと仮定し、リターン時系列の長さを1000とするとします.
まず最初は、numpy.random.randn
を用いて、正規分布に従うサンプルリターンを使います.
# number of assets n_assets = 4 # number of observations n_obs = 1000 returns = np.random.randn(n_assets, n_obs)
これらのリターンによって、異なるリターンとリスクのポートフォリオの(選択)範囲を作り出すことができます。すべての資産を投資するので、重みベクトルの総和が必要になります。
def rand_weights(n): ''' Produces n random weights that sum to 1''' k = np.random.randn(n) return k / sum(k)
例えば、以下のようになります。
[ 0.551635 -0.94880205 0.60015727 0.79700978] [ 1.44902231 0.1502679 -0.67985275 0.08056255]
期待リターン
平均リターン
Portfolioのそれぞれの重み
とする。例えば今回は4つの資産に投資すると仮定したので、となり、はスカラです。
リスク
リターンの共分散行列
とします。前回の中間発表でもしたように最小分散ポートフォリオ問題では、リターンとリスクさえわかっていればポートフォリオを作り出すことができます。
plotしたときにわかりやすくするため、リスクが1以上なものは除外します。
def random_portfolio(returns): ''' Returns the mean and standard deviation of returns for a random portfolio ''' p = np.asmatrix(np.mean(returns, axis=1)) w = np.asmatrix(rand_weights(returns.shape[0])) C = np.asmatrix(np.cov(returns)) mu = p@w.T sigma = np.sqrt(w * C * w.T) # This recursion reduces outliers to keep plots pretty if sigma > 1: return random_portfolio(returns) return mu, sigma n_portfolios = 500 print(random_portfolio(returns)) means, stds = np.column_stack([ random_portfolio(returns) for _ in range(n_portfolios) ]) # print(means)
結果としては
(matrix([[-0.01047962]]), matrix([[ 0.68538163]]))
Frontier
フロンティアは、最小分散ポートフォリオのリターンがプラスで、リスクが最小になるようなものを選択するためのものです。
重みの総和がになるような制約のもとでの最適化問題なので、scipy.optimise.minimize
かcvxopt
パッケージを用いられます。
def optimal_portfolio(returns): ret_len = len(returns) returns = np.asmatrix(returns) N = 100 mus = [10**(5. * t/N - 1.) for t in range(N)] # Convert to cvxopt matrices S = opt.matrix(np.cov(returns)) pbar = opt.matrix(np.mean(returns, axis=1)) # Create constraint matrices G = -opt.matrix(np.eye(ret_len)) # negative n x n identity matrix h = opt.matrix(0., (ret_len,1)) A = opt.matrix(1., (1,ret_len)) b = opt.matrix(1.) # Calculate efficient frontier weights using quadratic programming portfolios = [solvers.qp(mu*S, -pbar, G, h, A, b)['x'] for mu in mus] # calculate risks and returns for frontiner returns = [blas.dot(pbar, x) for x in portfolios] risks = [np.sqrt(blas.dot(x, S*x)) for x in portfolios] # calculate the 2nd degree polynomial of the frontier curve m1 = np.polyfit(returns, risks, 2) x1 = np.sqrt(m1[2] / m1[0]) # calculate the optimal portfolio wt = solvers.qp(opt.matrix(x1 * S), -pbar, G, h, A, b)['x'] return np.asarray(wt), returns, risks returns = np.random.randn(n_assets, n_obs) weights, returns, risks = optimal_portfolio(returns)
重みは以下のようになりました。
[[ 2.09381842e-08] [ 5.37940204e-02] [ 9.46205311e-01] [ 6.48181669e-07]]
ソースコード
さいごに
これ、適当にメモなのでそこまで深くかけてません。(2.5)でアルゴリズムをもう少し詳しく書こうと思います。