작성자: admin 작성일시: 2017-05-26 18:24:44 조회수: 1267 다운로드: 110
카테고리: 금융 공학 태그목록:

포트폴리오 최적화

In:
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt

주가 데이터 수집

In:
import pandas_datareader.data as wb

symbols = ['AAPL', 'MSFT', 'YHOO', 'DB', 'GLD']
noa = len(symbols)
data = pd.DataFrame()
for sym in symbols:
    data[sym] = wb.DataReader(sym, data_source='google',
                              end='2014-09-12')['Close']
data.columns = symbols
data.tail()
Out:
AAPL MSFT YHOO DB GLD
Date
2014-09-08 98.36 46.47 41.81 34.92 120.73
2014-09-09 97.99 46.76 40.78 34.71 120.87
2014-09-10 101.00 46.84 41.14 35.29 120.26
2014-09-11 101.43 47.00 41.26 35.25 119.47
2014-09-12 101.66 46.70 42.88 35.03 118.38
In:
(data / data.iloc[0] * 100).plot(figsize=(8, 5), grid=True)
plt.show()

수익률 계산

In:
rets = np.log(data / data.shift(1))
rets.mean() * 252
Out:
AAPL    0.256616
MSFT    0.087852
YHOO    0.196331
DB     -0.157479
GLD     0.016068
dtype: float64
In:
rets.cov() * 252
Out:
AAPL MSFT YHOO DB GLD
AAPL 0.072982 0.020341 0.023218 0.041231 0.005259
MSFT 0.020341 0.049669 0.024301 0.046324 0.002122
YHOO 0.023218 0.024301 0.093255 0.051310 -0.000810
DB 0.041231 0.046324 0.051310 0.178445 0.008726
GLD 0.005259 0.002122 -0.000810 0.008726 0.032505
In:
sns.heatmap(rets.cov() * 252, cmap=mpl.cm.bone_r)
plt.show()

포트폴리오 이론

  • 각 자산의 수익률이 $r_i$이고 자산 비중이 $w_i$이면
  • 포트폴리오 수익률은 $$ \begin{eqnarray} \mu_p &=& \text{E} [r_p] \\ &=& \text{E} \left[ \sum_i w_i r_i \right] \\ &=& \sum_i w_i \text{E} [r_i] \\ &=& \sum_i w_i \mu_i \\ &=& w^T \mu \end{eqnarray} $$ 이 식에서 $\mu$ 는 자산의 기대 수익률 벡터
  • 포트폴리오 분산은 $$ \begin{eqnarray} \sigma_p^2 &=& \text{E} \left[ (r_p - \mu_p)^2 \right] \\ &=& \text{E} \left[ (w^T r - w^T \mu)^2 \right] \\ &=& \text{E} \left[ w^T(r-\mu)(r-\mu)^Tw \right] \\ &=& w^T \text{E} \left[ (r-\mu)(r-\mu)^T \right] w \\ &=& w^T \Sigma w \end{eqnarray} $$ 이 식에서 $\Sigma$ 는 자산 수익률의 공분산 행렬
In:
np.random.seed(2)
weights = np.random.random(noa)
weights /= np.sum(weights)
weights
Out:
array([ 0.23349275,  0.01388454,  0.2943663 ,  0.2331326 ,  0.22512381])
In:
np.sum(rets.mean() * weights) * 252
Out:
0.08583483808126474
In:
np.dot(weights.T, np.dot(rets.cov() * 252, weights))
Out:
0.040143219666608443
In:
np.sqrt(np.dot(weights.T, np.dot(rets.cov() * 252, weights)))
Out:
0.20035772924099646

포트폴리오 수익률 시뮬레이션

In:
prets = []
pvols = []
for p in range (2500):
    weights = np.random.random(noa)
    weights /= np.sum(weights)
    prets.append(np.sum(rets.mean() * weights) * 252)
    pvols.append(np.sqrt(np.dot(weights.T, 
                        np.dot(rets.cov() * 252, weights))))
prets = np.array(prets)
pvols = np.array(pvols)
In:
plt.scatter(pvols, prets, c=prets/pvols, marker='o', cmap=mpl.cm.jet)
plt.grid(True)
plt.xlabel('expected volatility')
plt.ylabel('expected return')
plt.colorbar(label='Sharpe ratio')
plt.show()

포트폴리오 통계

  • 수익률, 변동성, 샤프 지수 계산
In:
def statistics(weights):
    ''' Return portfolio statistics.
    
    Parameters
    ==========
    weights : array-like
        포트폴리오 내의 증권 비중
    
    Returns
    =======
    pret : float
        포트폴리오 수익률의 기댓값
    pvol : float
        포트폴리오 변동성의 기댓값
    pret / pvol : float
        무위험 이자율이 0일 때의 포트폴리오 샤프 지수
    '''
    weights = np.array(weights)
    pret = np.sum(rets.mean() * weights) * 252
    pvol = np.sqrt(np.dot(weights.T, np.dot(rets.cov() * 252, weights)))
    return np.array([pret, pvol, pret / pvol])

샤프 지수 최대화

In:
import scipy.optimize as sco
In:
def min_func_sharpe(weights):
    return -statistics(weights)[2]
In:
cons = ({'type': 'eq', 'fun': lambda x:  np.sum(x) - 1})
In:
bnds = tuple((0, 1) for x in range(noa))
bnds
Out:
((0, 1), (0, 1), (0, 1), (0, 1), (0, 1))
In:
noa * [1. / noa,]
Out:
[0.2, 0.2, 0.2, 0.2, 0.2]
In:
%%time
opts = sco.minimize(min_func_sharpe, noa * [1. / noa,], method='SLSQP',
                       bounds=bnds, constraints=cons)
CPU times: user 140 ms, sys: 0 ns, total: 140 ms
Wall time: 147 ms
In:
opts
Out:
     fun: -1.0273846089134062
     jac: array([ -1.75386667e-05,   3.25929075e-02,   4.06056643e-05,
         1.53298907e+00,  -2.69711018e-06])
 message: 'Optimization terminated successfully.'
    nfev: 35
     nit: 5
    njev: 5
  status: 0
 success: True
       x: array([  6.93610646e-01,   2.00143721e-16,   3.00136168e-01,
         0.00000000e+00,   6.25318628e-03])
In:
opts['x'].round(3)
Out:
array([ 0.694,  0.   ,  0.3  ,  0.   ,  0.006])
In:
statistics(opts['x']).round(3)
Out:
array([ 0.237,  0.231,  1.027])
In:
plt.scatter(pvols, prets, c=prets/pvols, marker='o', cmap=mpl.cm.jet)
plt.grid(True)
plt.xlabel('expected volatility')
plt.ylabel('expected return')
plt.colorbar(label='Sharpe ratio')

pt_opts = statistics(opts['x']).round(3)
plt.scatter(pt_opts[1], pt_opts[0], marker="*", s=500, alpha=0.5)
plt.show()