In [86]:
import numpy as np
import pandas as pd
import numbers
# adopted from https://www.mathworks.com/matlabcentral/fileexchange/36487-weighted-generalized-hurst-exponent
def genhurstw(S, q = 1.0, delta = np.inf, maxT = 20, minT = 5, ddim = 0):
import numpy as np
import pandas as pd
import numbers
q = float(q) if isinstance(q, numbers.Integral) else q
delta = float(delta) if isinstance(delta, numbers.Integral) else delta
S = np.array(S) if isinstance(S, list) else (S.values if isinstance(S, (pd.DataFrame, pd.Series, pd.Panel, pd.Panel4D)) else S)
q = np.array(q) if isinstance(q, (list, numbers.Number)) else (q.values.copy() if isinstance(q, (pd.DataFrame, pd.Series, pd.Panel, pd.Panel4D)) else q)
delta = np.array(delta) if isinstance(delta, (list, numbers.Number)) else (delta.values.copy() if isinstance(delta, (pd.DataFrame, pd.Series, pd.Panel, pd.Panel4D)) else delta)
lpidxd = S.shape[:ddim]
rpidxd = S.shape[(ddim + 1):]
lpidxi = (slice(None, None, None),) * ddim
rpidxi = (slice(None, None, None),) * (S.ndim - ddim - 1)
rs1d = ((1,) * ddim) + (-1,) + ((1,) * (S.ndim - ddim - 1))
q = q.reshape(-1).reshape(rs1d)
delta = delta.reshape(-1).reshape(rs1d)
L = S.shape[ddim]
lq = q.shape[ddim]
deltas = delta.shape[ddim]
H = np.empty(lpidxd + (maxT - minT + 1, lq, deltas) + rpidxd, dtype = float)
cc = np.empty(lpidxd + (2,) + rpidxd, dtype = float)
k = 0
for Tmax in range(minT, maxT + 1):
x = np.arange(1, Tmax + 1, dtype = float).reshape(rs1d)
mcord = np.empty(lpidxd + (Tmax, lq, deltas) + rpidxd, dtype = float)
for tt in range(1, Tmax + 1):
dV = S[lpidxi + (slice(tt, L, tt),) + rpidxi] - S[lpidxi + (slice(0, L - tt, tt),) + rpidxi]
VV = S[lpidxi + (slice(0, L, tt),) + rpidxi]
N = float(dV.shape[ddim]) + 1.0
Np1 = N + 1.0
X = np.arange(1.0, Np1, dtype = float).reshape(rs1d)
Y = VV
mx = np.sum(X, axis = ddim) / N
SSxx = np.sum(X ** 2, axis = ddim) - N * (mx ** 2)
my = np.sum(Y, axis = ddim) / N
SSxy = np.sum(X * Y, axis = ddim) - N * mx * my
cc[lpidxi + (0,) + rpidxi] = SSxy/SSxx
cc[lpidxi + (1,) + rpidxi] = my - cc[lpidxi + (0,) + rpidxi] * mx
ddVd = dV - cc[lpidxi + (0,) + rpidxi]
VVVd = VV - cc[lpidxi + (0,) + rpidxi] * X - cc[lpidxi + (1,) + rpidxi] # X = 1:N
for d in range(0, deltas):
cdelta = delta[lpidxi + (d,) + rpidxi]
if np.isfinite(cdelta):
# extra coputation but this gives the same results as below for Inf values
aph = 1.0 / cdelta
r = np.arange(0, Np1, dtype = float).reshape(rs1d)
w = np.exp(-aph * r[lpidxi + (slice(0, -2),) + rpidxi]) * (1.0 - np.exp(-aph)) / (1.0 - np.exp(-aph * (N - 1)))
w1 = np.exp(-aph * r[lpidxi + (slice(0, -1),) + rpidxi]) * (1.0 - np.exp(-aph)) / (1.0 - np.exp(-aph * N))
for qq in range(0, lq):
mcord[lpidxi + (tt - 1, qq, d) + rpidxi] = np.sum(w * (np.abs(ddVd) ** q[lpidxi + (qq,) + rpidxi]), axis = ddim) / np.sum(w1 * (np.abs(VVVd) ** q[lpidxi + (qq,) + rpidxi]), axis = ddim)
else:
for qq in range(0, lq):
mcord[lpidxi + (tt - 1, qq, d) + rpidxi] = np.mean(np.abs(ddVd) ** q[lpidxi + (qq,) + rpidxi], axis = ddim) / np.mean(np.abs(VVVd) ** q[lpidxi + (qq,) + rpidxi], axis = ddim)
mx = np.mean(np.log10(x), axis = ddim)
SSxx = np.sum(np.log10(x) ** 2, axis = ddim) - Tmax * (mx ** 2)
for d in range(0, deltas):
for qq in range(0, lq):
my = np.mean(np.log10(mcord[lpidxi + (slice(None), qq, d) + rpidxi]), axis = ddim)
SSxy = np.sum(np.log10(x) * np.log10(mcord[lpidxi + (slice(None), qq, d) + rpidxi]), axis = ddim) - Tmax * mx * my
H[lpidxi + (k, qq, d) + rpidxi] = SSxy/SSxx
k = k + 1
mH = np.empty(lpidxd + (lq,) + rpidxd + (deltas,), dtype = float)
sH = np.empty(lpidxd + (lq,) + rpidxd + (deltas,), dtype = float)
for d in range(0, deltas):
for qq in range(0, lq):
mH[lpidxi + (qq,) + rpidxi + (d,)] = np.mean(H[lpidxi + (slice(None), qq, d) + rpidxi], axis = ddim) / q[lpidxi + (qq,) + rpidxi]
sH[lpidxi + (qq,) + rpidxi + (d,)] = np.std(H[lpidxi + (slice(None), qq, d) + rpidxi], axis = ddim) / q[lpidxi + (qq,) + rpidxi]
return mH, sH
In [87]:
rndn = np.cumsum(np.random.normal(0.0, 1.0, (10000, 20)), axis = 0)
import matplotlib.pyplot as plt
plt.plot(rndn)
plt.ylabel('random numbers')
plt.show()
In [88]:
H = genhurstw(rndn)
H
Out[88]:
In [89]:
H = genhurstw(rndn, np.array([1, 2]))
H
Out[89]:
In [90]:
H = genhurstw(rndn, np.array([1, 2]), np.exp(1000.0))
H
Out[90]:
In [91]:
H = genhurstw(rndn, np.array([1, 2]), 252.0)
H
Out[91]:
In [92]:
p = get_pricing(
symbols(['SPY', 'XLP', 'XLE', 'XLK', 'XLY', 'XLU', 'XLB', 'XLF', 'XLV','XLI', 'MDY', 'EWA', 'EWC', 'EWG', 'EWH', 'EWJ', 'EWW', 'EWS', 'EWU', 'EZU', 'EFA']),
fields = 'price',
start_date = '2002-01-04',
end_date = '2016-07-29',
frequency = 'daily',
)
p.fillna(method = 'ffill', inplace = True)
p = p / p.iloc[0, :]
p.plot()
Out[92]:
In [93]:
H = genhurstw(p, [1, 2], 252.0)
H
Out[93]:
In [94]:
((H[0][0, :, :] - 0.5) / H[1][0, :, :]).squeeze()
Out[94]:
In [95]:
((H[0][1, :, :] - 0.5) / H[1][1, :, :]).squeeze()
Out[95]:
In [96]:
H = genhurstw(p, [1, 2], [252.0, 126.0, 63.0, 21.0])
H
Out[96]:
In [97]:
((H[0][0, :, :] - 0.5) / H[1][0, :, :]).squeeze()
Out[97]:
In [98]:
((H[0][1, :, :] - 0.5) / H[1][1, :, :]).squeeze()
Out[98]:
In [ ]: