読者です 読者をやめる 読者になる 読者になる

主に、強化学習

情報系の大学2年生が確率に関連したことを多めに書いてるブログ

粒子フィルタ使って実装したい (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)

f:id:reonreon3reon:20161010103428p:plain

これらのリターンによって、異なるリターンとリスクのポートフォリオの(選択)範囲を作り出すことができます。すべての資産を投資するので、重みベクトルの総和が必要になります。

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]

 R = \mathbf{p^{T} w}

 \hspace{12pt} R \equiv \, \, 期待リターン
 \hspace{12pt} \mathbf{p^{T}} \equiv \, \, 平均リターン
 \hspace{12pt} \mathbf{w} \equiv \, \, Portfolioのそれぞれの重み

とする。例えば今回は4つの資産に投資すると仮定したので、 p, w \in R^{4}となり、Rはスカラです。

 \sigma = \sqrt{\mathbf{w^{T} C w}}

 \hspace{12pt} \sigma \equiv \, \, リスク
 \hspace{12pt} \mathbf{C} \equiv \, \, リターンの共分散行列

とします。前回の中間発表でもしたように最小分散ポートフォリオ問題では、リターンとリスクさえわかっていればポートフォリオを作り出すことができます。

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

f:id:reonreon3reon:20161010103957p:plain

Frontier

フロンティアは、最小分散ポートフォリオのリターンがプラスで、リスクが最小になるようなものを選択するためのものです。

重み \mathbf{w}の総和が1になるような制約のもとでの最適化問題なので、scipy.optimise.minimizecvxoptパッケージを用いられます。

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)

f:id:reonreon3reon:20161010104205p:plain

重みは以下のようになりました。

[[  2.09381842e-08]
 [  5.37940204e-02]
 [  9.46205311e-01]
 [  6.48181669e-07]]

ソースコード

さいごに

これ、適当にメモなのでそこまで深くかけてません。(2.5)でアルゴリズムをもう少し詳しく書こうと思います。