おすすめの投資先は金（ゴールド）！プログラミングで投資分析をした話

#import libraries
from pandas_datareader import data
import pandas as pd
from pylab import mpl, plt
import numpy as np
import datetime
plt.style.use('seaborn')
mpl.rcParams['font.family'] = 'serif'
%matplotlib inline

#date
start = '2000-01-01'
end = datetime.date.today()

#get date from yahoo! finance
def price_df(item, start, end):
    file = data.DataReader(item, 'yahoo', start, end)
    date = file.index
    price = file['Close']
    return date, price

#S&P 500 price
SP500_date, SP500_price = price_df('^GSPC', start, end)

#Gold price
gold_date, gold_price = price_df('GC=F', start, end)

#show a price chart
plt.figure(figsize=(10, 6))
plt.plot(gold_date, gold_price, label = 'Gold', color='g')
plt.plot(SP500_date, SP500_price, label = 'S&P 500', color='b')
plt.legend(loc=2, fontsize=18)
plt.xlabel('year', fontsize=18)
plt.ylabel('price(USD)', fontsize=18)
plt.tick_params(labelsize=18)
plt.savefig("SP500_gold_1.png")

#make a DataFrame
df = pd.DataFrame()
df['gold'] = gold_price
df['SP500'] = SP500_price

#calculate log returns
rets = np.log(df / df.shift(1))
rets.dropna(inplace=True)

return_date = rets.index
gold_return = rets['gold'].cumsum().apply(np.exp)
SP500_return = rets['SP500'].cumsum().apply(np.exp)

#show a return chart
plt.figure(figsize=(10, 6))
plt.plot(return_date, gold_return, label = 'Gold', color='g')
plt.plot(return_date, SP500_return, label = 'S&P 500', color='b')
plt.legend(loc=2, fontsize=18)
plt.xlabel('year', fontsize=18)
plt.ylabel('Return', fontsize=18)
plt.tick_params(labelsize=18)
plt.savefig("SP500_gold_2.png")

#OLS analysis
reg = np.polyfit(rets['gold'], rets['SP500'], deg=1)

#scatter plot
ax = rets.plot(kind='scatter', x='gold', y='SP500', figsize=(10,6), fontsize=18)
ax.plot(rets['gold'], np.polyval(reg, rets['gold']), 'r', lw=2)
plt.xlabel('Gold', fontsize=18)
plt.ylabel('S&P500', fontsize=18)
plt.savefig("SP500_gold_3.png")

#correlation
rets.corr()

#import libraries
import math
import scipy.stats as scs
import statsmodels.api as sm

#number of item to invest
noa = 2

#calculate return
def port_ret(weights):
    return np.sum(rets.mean() * weights) * 252

#calculate volatility
def port_vol(weights):
    return np.sqrt(np.dot(weights.T, np.dot(rets.cov() * 252, weights)))

#make plot data
prets = []
pvols = []
for p in range(2500):
    weights = np.random.random(noa)
    weights /= np.sum(weights)
    prets.append(port_ret(weights))
    pvols.append(port_vol(weights))
prets = np.array(prets)
pvols = np.array(pvols)

#import libraries
import scipy.optimize as sco

#calculate sharp ratio
def min_func_sharpe(weights):
    return -port_ret(weights) / port_vol(weights)

#restriction
cons = ({'type':'eq', 'fun':lambda x: np.sum(x) - 1})
bnds = tuple((0, 1) for x in range(noa))

#initial raio
eweights = np.array(noa * [1. / noa,])

#calculate the best value of sharp ratio
%%time
opts = sco.minimize(min_func_sharpe, eweights,
                   method='SLSQP', bounds=bnds,
                   constraints=cons)

#the result of calculation
opts['x'].round(3)
port_ret(opts['x']).round(3)
port_vol(opts['x']).round(3)
port_ret(opts['x']) / port_vol(opts['x'])

#calculate minimum volatility
optv = sco.minimize(port_vol, eweights,
                   method='SLSQP', bounds=bnds,
                   constraints=cons)

#the result of calculation
optv['x'].round(3)
port_vol(optv['x']).round(3)
port_ret(optv['x']).round(3)
port_ret(optv['x']) / port_vol(optv['x'])

#show a efficient frontier
plt.figure(figsize=(10, 6))
plt.scatter(pvols, prets, c=prets / pvols,
           marker='.', alpha=0.8, cmap='coolwarm')
plt.plot(port_vol(opts['x']), port_ret(opts['x']),
        'y*', markersize=15.0)
plt.plot(port_vol(optv['x']), port_ret(optv['x']),
        'r*', markersize=15.0)
plt.xlabel('expected volatility')
plt.ylabel('expected return')
plt.colorbar(label='Sharpe ratio')
plt.savefig("SP500_gold_5.png");