-
Notifications
You must be signed in to change notification settings - Fork 14
/
quant_research.py
163 lines (151 loc) · 5.73 KB
/
quant_research.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
import os
import math
import numpy as np
import pandas as pd
from scipy.stats import norm
import matplotlib.pyplot as plt
from backend.config import currencies
file = 'history/AUD_JPY_D.csv'
folder = 'history/'
def monteCarloSimulation(file, n_days, n_iters):
"""
file: csv file containing at least a 'Close' column with asset prices inside
n_days: the number of days to simulate retruns for
n_iters: the number of different simulations to run
returns the price_matrix of the simulation
"""
# import data
df = pd.read_csv(file)
prices = df['Close']
# set formula's variables
log_returns = np.log(1 + prices.pct_change())
u = log_returns.mean()
var = log_returns.var()
drift = u - (0.5 * var)
std = log_returns.std()
drift = pd.Series(drift)
std = pd.Series(std)
z = norm.ppf(np.random.rand(n_days, n_iters))
# generate future daily returns
daily_returns = np.exp(drift.values + std.values * z)
p0 = prices.iloc[-1] # last price in our original data is the first price in our simulation
price_matrix = np.zeros_like(daily_returns)
price_matrix[0] = p0
# start the simulation
for t in range(1, n_days):
price_matrix[t] = price_matrix[t-1] * daily_returns[t]
plt.plot(price_matrix)
plt.show()
return price_matrix
# Brownian motion
def Brownian(seed, N):
np.random.seed(seed)
dt = 1/N # timestep
b = np.random.normal(0., 1., int(N)) * np.sqrt(dt) # Brownian normal increments
W = np.cumsum(b) #Brownian path
return W, b
# Geometric Brownian Motion
def GBM(p0, mu, sigma, W, T, N):
t = np.linspace(0., 1., N+1)
P = list()
P.append(p0)
for i in range(1, int(N + 1)):
drift = (mu - 0.5 * sigma ** 2) * t[i]
diffusion = sigma * W[i-1]
p_temp = p0 * np.exp(drift + diffusion)
P.append(p_temp)
return P, t
# Euler Maruyama Approximation
def eulerApprox(p0, mu, sigma, b, T, N, M):
dt = M * (1/N) # step size
L = N / M
wi = [p0]
for i in range(0, int(L)):
Winc = np.sum(b[(M * (i - 1) + M) : (M * i + M)])
w_i_new = wi[i] + mu * wi[i] * dt + sigma * wi[i] * Winc
wi.append(w_i_new)
return wi, dt
# plot these stochastic compared to a given asset
def plot_stochastics(file, seed1, seed2, seed3):
df = pd.read_csv(file)
n = len(df) - 1
returns = df['Close'].pct_change()
p0 = df['Close'][n]
N = 2.0 ** 6
W1 = Brownian(seed1, N)[0]
W2 = Brownian(seed2, N)[0]
W3 = Brownian(seed3, N)[0]
T = 1.0
mu = np.mean(returns) * 252.0
sigma = np.std(returns) * np.sqrt(252.0)
gbm1 = GBM(p0, mu, sigma, W1, T, N)[0]
gbm2 = GBM(p0, mu, sigma, W2, T, N)[0]
gbm3 = GBM(p0, mu, sigma, W3, T, N)[0]
t = GBM(p0, mu, sigma, W1, T, N)[1]
plt.plot(t, gbm1, label='GMB1', ls='--')
plt.plot(t, gbm2, label='GBM2', ls='--')
plt.plot(t, gbm3, label='GBM3', ls='--')
plt.plot(t, df['Close'][-65:], label='Actual')
plt.ylabel('Asset Price, $')
plt.title('Geometric Brownian Motion')
plt.legend(loc='upper left')
plt.show()
# Markowitz Portfolio optimization, efficient frontier, MPT, etc.
def optimalPortfolio(asset_list, n_portfolios=2500):
"""
asset_list: list of asset symbols to extract
type: CSV or GET; CSV from local file or GET from yahoo
"""
data = pd.DataFrame()
# get asset pricing from local data
if type == 'csv':
for asset in asset_list:
folder = 'history/'
name_window = '%s_D' % asset
for files in os.walk(folder):
for f in files[2]:
if name_window in f:
f = os.path.join(folder, f)
df = pd.read_csv(f)
data[asset] = df['Close']
elif type == 'GET':
for asset in asset_list:
s = stock(asset).chart_table(range='1y')
data[asset] = s['close']
# convert prices to returns and calculate daily mean returns and covariance
returns = data.pct_change()
n = returns.shape[1]
mean_daily_returns = returns.mean()
cov_matrix = returns.cov()
# set array to hold results
results = np.zeros((3+n, n_portfolios))
# calculate the n portfolios
for i in range(n_portfolios):
# set random weights; rebalance weights to sum 1
weights = np.random.random(n)
weights /= np.sum(weights)
# calculate portfolio return and volatility
port_return = np.sum(mean_daily_returns * weights) * 252
port_std = np.sqrt(np.dot(weights.T, np.dot(cov_matrix, weights))) * np.sqrt(252)
# store results and calculate/store Sharpe ratio
results[0,i] = port_return
results[1,i] = port_std
results[2,i] = results[0,i] / results[1,i]
# iterate through weight vector and add weights to the results
for j in range(len(weights)):
results[j+3,i] = weights[j]
# plot results
cols = ['Returns', 'Std', 'Sharpe'] + asset_list
results_frame = pd.DataFrame(results.T, columns=cols)
max_sharpe_port = results_frame.iloc[results_frame['Sharpe'].idxmax()]
min_vol_port = results_frame.iloc[results_frame['Std'].idxmin()]
plt.scatter(results_frame.Std, results_frame.Returns, c=results_frame.Sharpe, cmap='RdYlBu')
plt.colorbar()
plt.ylabel('Returns')
plt.xlabel('Volatility')
plt.title('Portfolio optimization test for %d portfolios with %d assets' % (n_portfolios, n))
plt.scatter(max_sharpe_port[1], max_sharpe_port[0], marker=(5,1,0), color='r', s=1000)
plt.scatter(min_vol_port[1], min_vol_port[0], marker=(5,1,0), color='g', s=1000)
plt.show()
print("Portfolio weights with the best returns: ", max_sharpe_port)
print("Portfolio weights with the lowest volatility: ", min_vol_port)