程序代码：

# -*- coding: utf-8 -*-

#   NumPyの読み込み
import numpy as np
#   Pandasの読み込み
import pandas as pd
#   Cvxoptの読み込み
import cvxopt as co
#   MatplotlibのPyplotモジュールの読み込み
import matplotlib.pyplot as plt
#   日本語フォントの設定


#   最小分散ポートフォリオの計算
def MV_Portfolio(ReturnData, TargetReturn):
    ##        ReturnData: 収益率データ
    ##      TargetReturn: 目標収益率
    ##            Output: 最適ポートフォリオ
    co.solvers.options['show_progress'] = False
    MeanVector = np.mean(ReturnData, axis=0)
    Deviation = ReturnData - MeanVector
    NumberOfPeriods = ReturnData.shape[0]
    NumberOfAssets = ReturnData.shape[1]
    NumberOfElements = NumberOfAssets + NumberOfPeriods
    P = co.spmatrix(1.0/NumberOfPeriods, \
                    range(NumberOfAssets, NumberOfElements), \
                    range(NumberOfAssets, NumberOfElements))
    q = co.matrix(0.0, (NumberOfElements, 1))
    G = co.spmatrix(-1.0, range(NumberOfAssets), range(NumberOfAssets), \
                    (NumberOfAssets, NumberOfElements))
    h = co.matrix(0.0, (NumberOfAssets, 1))
    A = co.sparse([ \
        co.sparse([ \
            [co.matrix(Deviation)], \
            [co.spmatrix(-1.0, range(NumberOfPeriods), range(NumberOfPeriods))]\
        ]), \
        co.sparse([ \
            [co.matrix(np.c_[MeanVector, np.ones(NumberOfAssets)]).trans()], \
            [co.spmatrix([], [], [], (2, NumberOfPeriods))] \
        ]) \
        ])
    b = co.matrix(np.r_[np.zeros(NumberOfPeriods), TargetReturn, 1.0])
    return co.solvers.qp(P, q, G, h, A, b)

#   最小絶対偏差ポートフォリオの計算
def AD_Portfolio(ReturnData, TargetReturn):
    ##        ReturnData: 収益率データ
    ##      TargetReturn: 目標収益率
    ##            Output: 最適ポートフォリオ
    co.solvers.options['show_progress'] = False
    MeanVector = np.mean(ReturnData, axis=0)
    Deviation = ReturnData - MeanVector
    NumberOfPeriods = ReturnData.shape[0]
    NumberOfAssets = ReturnData.shape[1]
    NumberOfElements = NumberOfAssets + NumberOfPeriods
    c = co.matrix([co.matrix(0.0, (NumberOfAssets, 1)), \
                   co.matrix(1.0/NumberOfPeriods, (NumberOfPeriods, 1))])
    G = co.sparse([ \
        co.spmatrix(-1.0, range(NumberOfPeriods), \
                    range(NumberOfAssets, NumberOfElements)), \
        co.sparse([ \
            [co.matrix(-Deviation)], \
            [co.spmatrix(-1.0, range(NumberOfPeriods), range(NumberOfPeriods))]\
        ]), \
        co.spmatrix(-1.0, range(NumberOfAssets), range(NumberOfAssets), \
                     (NumberOfAssets, NumberOfElements)) \
        ])
    h = co.matrix(0.0, (2*NumberOfPeriods + NumberOfAssets, 1))
    A = co.sparse([ \
            [co.matrix(np.c_[MeanVector, np.ones(NumberOfAssets)]).trans()], \
            [co.spmatrix([], [], [], (2, NumberOfPeriods))] \
        ])
    b = co.matrix([TargetReturn, 1.0])
    return co.solvers.lp(c, G, h, A, b)

#   最適期待ショートフォール・ポートフォリオの計算
def ES_Portfolio(ReturnData, TargetReturn, Alpha):
    ##        ReturnData: 収益率データ
    ##      TargetReturn: 目標収益率
    ##             Alpha: 下方確率
    ##            Output: 最適ポートフォリオ
    co.solvers.options['show_progress'] = False
    MeanVector = np.mean(ReturnData, axis=0)
    NumberOfPeriods = ReturnData.shape[0]
    NumberOfAssets = ReturnData.shape[1]
    NumberOfElements = NumberOfAssets + NumberOfPeriods
    c = co.matrix( \
        [co.matrix(0.0, (NumberOfAssets, 1)), \
         co.matrix(1.0/(Alpha*NumberOfPeriods), (NumberOfPeriods, 1)), -1.0])
    G = co.sparse([ \
        co.spmatrix(-1.0, range(NumberOfPeriods), \
                    range(NumberOfAssets, NumberOfElements), \
                    (NumberOfPeriods, NumberOfElements + 1)), \
        co.sparse([ \
        [co.matrix(-ReturnData)], \
        [co.spmatrix(-1.0, range(NumberOfPeriods), range(NumberOfPeriods))], \
        [co.matrix(1.0, (NumberOfPeriods, 1))] \
        ]), \
        co.spmatrix(-1.0, range(NumberOfAssets), range(NumberOfAssets), \
                     (NumberOfAssets, NumberOfElements + 1)) \
        ])
    h = co.matrix(0.0, (2*NumberOfPeriods + NumberOfAssets, 1))
    A = co.sparse([ \
            [co.matrix(np.c_[MeanVector, np.ones(NumberOfAssets)]).trans()], \
            [co.spmatrix([], [], [], (2, NumberOfPeriods + 1))] \
        ])
    b = co.matrix([TargetReturn, 1.0])
    return co.solvers.lp(c, G, h, A, b)

#   株価データの読み込み
stockvalue = pd.read_csv('stocks.csv', index_col=0)
#   最小分散フロンティアの計算
R = (stockvalue.diff()/stockvalue.shift(1))[1:] * 100.0
T = R.shape[0]
mu = R.mean().values
Sigma = (R.cov().values)*((T-1.0)/T)
V_Target = np.linspace(mu.min(), mu.max(), num=50)
V_SD = np.array( \
        [np.sqrt(2.0*MV_Portfolio(R.values, Target)['primal objective']) \
         for Target in V_Target])
V_AD = np.array( \
        [2.0*AD_Portfolio(R.values, Target)['primal objective'] \
         for Target in V_Target])
V_ES = np.array( \
        [ES_Portfolio(R.values, Target, 0.1)['primal objective'] \
         for Target in V_Target])
#   フロンティアのグラフの作成
()
fig1 = plt.figure(1, facecolor='w')
plt.subplot(1, 3, 1)
plt.plot(V_SD, V_Target, 'b-')
plt.plot(np.sqrt(np.diagonal(Sigma)), mu, 'rx')
plt.xlabel(u'standard deviation')
plt.ylabel(u'Expected rate of return')
plt.subplot(1, 3, 2)
plt.plot(V_AD, V_Target, 'b-')
plt.plot((R - R.mean()).abs().mean().values, mu, 'rx')
plt.xlabel(u'Absolute deviation')
plt.subplot(1, 3, 3)
plt.plot(V_ES, V_Target, 'b-')
plt.plot((-R[R <= R.quantile(0.1)]).mean().values, mu, 'rx')
plt.xlabel(u'Expected shortfall')

[此贴子已经被作者于2017-2-5 20:11编辑过]