銀投資はおすすめ？Pythonのプログラムで分析した話

#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 price data 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

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

#silver
silver_date, silver_price = price_df('SI=F', start, end)

#make a DataFrame
df = pd.DataFrame()
df['gold'] = gold_price
df['silver'] = silver_price

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

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

#make a scatter
ax = rets.plot(kind='scatter', x='silver', y='gold', figsize=(10,6), fontsize=18)
ax.plot(rets['silver'], np.polyval(reg, rets['silver']), 'r', lw=2)
plt.xlabel('Silver', fontsize=18)
plt.ylabel('Gold', fontsize=18)
plt.savefig("silver_gold_1.png")

#correlation
rets.corr()

#get DGS10 from FRED
tx = data.DataReader('DGS10', 'fred', start, end).dropna()

#get T10YIE from FRED
ty = data.DataReader('T10YIE', 'fred', start, end).dropna()

# calculate real interest rate
test = tx['DGS10'] - ty['T10YIE']
test = test.dropna()
test_date = test.index

#show a graph
fig, ax1 = plt.subplots(figsize=(10, 6))
plt.plot(test_date, test, label = 'DGS10-T10YIE', color='g')
plt.legend(loc=2, fontsize=18)
plt.xlabel('year', fontsize=18)
plt.ylabel('percent(%)', fontsize=18)
plt.tick_params(labelsize=18)
ax2 = ax1.twinx()
plt.plot(silver_date, silver_price, label = 'Silver', color='b')
plt.legend(loc=6, fontsize=18)
plt.ylabel('Silver price(USD)', fontsize=18)
plt.tick_params(labelsize=18)
plt.savefig("silver2.png")

#DGS10
DGS10 = tx['DGS10']
DGS10_date = tx.index 

#T10YIE
T10YIE = ty['T10YIE']
T10YIE_date = ty.index

#monetary basa total
tz = data.DataReader('BOGMBASE', 'fred', start, end).dropna()
BOGMBASE_date = tz.index
BOGMBASE = tz['BOGMBASE']

#show a economic analysis
fig, ax1 = plt.subplots(figsize=(10, 6))
plt.plot(DGS10_date, DGS10, label = 'DGS10', color='g')
plt.plot(T10YIE_date, T10YIE, label = 'T10YIE', color='b')
plt.legend(loc=2, fontsize=18)
plt.xlabel('year', fontsize=18)
plt.ylabel('percent(%)', fontsize=18)
plt.tick_params(labelsize=18)
ax2 = ax1.twinx()
plt.plot(BOGMBASE_date, BOGMBASE, label = 'MonetaryBase', color='black')
plt.legend(loc=6, fontsize=18)
plt.ylabel('Millions of Dollars', fontsize=18)
plt.tick_params(labelsize=18)
plt.savefig("silver3.png")

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

#show  a price chart
fig, ax1 = plt.subplots(figsize=(10, 6))
plt.plot(SP500_date, SP500_price, label = 'S&P 500', color='g')
plt.plot(gold_date, gold_price, label = 'Gold', color='y')
plt.plot(silver_date, silver_price, label = 'Silver', 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("silver4.png")

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

#calculate returns
rets = np.log(df / df.shift(1))
date = df.index
gold_return = rets['gold'].cumsum().apply(np.exp)
silver_return = rets['silver'].cumsum().apply(np.exp)
SP500_return = rets['SP500'].cumsum().apply(np.exp)

#show a return graph
plt.figure(figsize=(10, 6))
plt.plot(date, SP500_return, label = 'S&P 500', color='g')
plt.plot(date, gold_return, label = 'Gold', color='y')
plt.plot(date, silver_return, label = 'Silver', 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("silver5.png")

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

#number of symbols
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)))

#calculate 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)

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

#start point
eweights = np.array(noa * [1. / noa,])

#calculate the best 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)

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

#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("silver7.png");