factors_r = ["SP500", "DTWEXAFEGS"] # "SP500" does not contain dividends; note: "DTWEXM" discontinued as of Jan 2020
factors_d = ["DGS10", "BAMLH0A0HYM2"]Price momentum
One month reversal and 2-12 month momentum are two ends of the spectrum. The general trend indicates that positive acceleration leads to reversals and negative acceleration leads to rebounds, i.e. acceleration is not sustainable. An unsustainable acceleration leading to reversal can reconcile the one-month reversal and 2-12 month momentum.
order = 20# "Momentum, Acceleration, and Reversal"
momentum_df = (np.exp(returns_df.shift(order).rolling(width - order, min_periods = 1).sum()) - 1).dropna()
Time-series score
Suppose we are looking at \(n\) independent and identically distributed random variables, \(X_{1},X_{2},\ldots,X_{n}\). Since they are iid, each random variable \(X_{i}\) has to have the same mean, which we will call \(\mu\), and variance, which we will call \(\sigma^{2}\):
\[ \begin{aligned} \mathrm{E}\left(X_{i}\right)&=\mu\\ \mathrm{Var}\left(X_{i}\right)&=\sigma^{2} \end{aligned} \]
Let’s suppose we want to look at the average value of our \(n\) random variables:
\[ \begin{aligned} \bar{X}=\frac{X_{1}+X_{2}+\cdots+X_{n}}{n}=\left(\frac{1}{n}\right)\left(X_{1}+X_{2}+\cdots+X_{n}\right) \end{aligned} \]
We want to find the expected value and variance of the average, \(\mathrm{E}\left(\bar{X}\right)\) and \(\mathrm{Var}\left(\bar{X}\right)\).
Expected value
\[ \begin{aligned} \mathrm{E}\left(\bar{X}\right)&=\mathrm{E}\left[\left(\frac{1}{n}\right)\left(X_{1}+X_{2}+\cdots+X_{n}\right)\right]\\ &=\left(\frac{1}{n}\right)\mathrm{E}\left(X_{1}+X_{2}+\cdots+X_{n}\right)\\ &=\left(\frac{1}{n}\right)\left(n\mu\right)\\ &=\mu \end{aligned} \]
Variance
\[ \begin{aligned} \mathrm{Var}\left(\bar{X}\right)&=\mathrm{Var}\left[\left(\frac{1}{n}\right)\left(X_{1}+X_{2}+\cdots+X_{n}\right)\right]\\ &=\left(\frac{1}{n}\right)^{2}\mathrm{Var}\left(X_{1}+X_{2}+\cdots+X_{n}\right)\\ &=\left(\frac{1}{n}\right)^{2}\left(n\sigma^{2}\right)\\ &=\frac{\sigma^{2}}{n} \end{aligned} \]
def sd(x):
n_rows = sum(~np.isnan(x))
if (n_rows > 1):
result = np.sqrt(np.nansum(x ** 2) / (n_rows - 1))
else:
result = np.nan
return result# volatility scale only
score_df = (momentum_df / momentum_df.rolling(width, min_periods = 1).apply(sd, raw = False)).dropna()overall_df = score_df.mean(axis = 1)
overall_df = overall_df / overall_df.rolling(width, min_periods = 1).apply(sd, raw = False)score_df.insert(loc = 0, column = "Overall", value = overall_df)
score_df = score_df.dropna()
Outlier detection
import statsmodels.api as smInterquartile range
Outliers are defined as the regression residuals that fall below \(Q_{1}-1.5\times IQR\) or above \(Q_{3}+1.5\times IQR\):
- https://stats.stackexchange.com/a/1153
- https://stats.stackexchange.com/a/108951
- https://robjhyndman.com/hyndsight/tsoutliers/
def outliers(z):
n_cols = z.shape[1]
result_ls = []
for j in range(n_cols):
y = z.iloc[:, j]
trend = pd.DataFrame({"trend": range(len(y))}, index = y.index)
if (n_cols == 1):
x = sm.add_constant(trend)
else:
x = sm.add_constant(pd.concat([trend, z.drop(z.columns[j], axis = 1)], axis = 1))
coef = sm.OLS(y, x).fit().params
predict = coef.iloc[0] + np.dot(x.iloc[:, 1:], coef.iloc[1:])
resid = y - predict
lower = resid.quantile(0.25)
upper = resid.quantile(0.75)
iqr = upper - lower
total = y[(resid < lower - 1.5 * iqr) | (resid > upper + 1.5 * iqr)]
total = pd.DataFrame({"date": total.index, "symbol": total.name, "values": total})
result_ls.append(total)
result = pd.concat(result_ls, ignore_index = True)
result = result.pivot_table(index = "date", columns = "symbol", values = "values")
return resultoutliers_df = outliers(score_df[factors])
Contour ellipsoid
The contours of a multivariate normal (MVN) distribution are ellipsoids centered at the mean. The directions of the axes are given by the eigenvectors of the covariance matrix and squared lengths are given by the eigenvalues:
\[ \begin{aligned} ({\mathbf{x}}-{\boldsymbol{\mu}})^{\mathrm{T}}{\boldsymbol{\Sigma}}^{-1}({\mathbf{x}}-{\boldsymbol{\mu}})=c^{2} \end{aligned} \]
Or, in general parametric form:
\[ \begin{aligned} X(t)&=X_{c}+a\,\cos t\,\cos \varphi -b\,\sin t\,\sin \varphi\\ Y(t)&=Y_{c}+a\,\cos t\,\sin \varphi +b\,\sin t\,\cos \varphi \end{aligned} \] where \(t\) varies from \(0,\ldots,2\pi\). Here \((X_{c},Y_{c})\) is the center of the ellipse and \(\varphi\) is the angle between the x-axis and the major axis of the ellipse.
Specifically:
\[ \begin{aligned} &\text{Center: }\boldsymbol{\mu}=(X_{c},Y_{c})\\ &\text{Radius: }c^{2}=\chi_{\alpha}^{2}(df)\\ &\text{Length: }a=c\sqrt{\lambda_{k}}\\ &\text{Angle of rotation: }\varphi=\text{atan2}\left(V_{k}(2),V_{k}(1)\right) \end{aligned} \]
Note: the radius converts a normal quantile into a chi-squared quantile, i.e. \(c^{2}=\chi_{\alpha}^{2}(2)\) with \(\alpha=2\Phi(\sigma)-1\), so a one-sigma ellipse covers \(2\Phi(1)-1\approx68\%\) of the distribution, i.e. the familiar two-sided normal coverage.
# https://www.visiondummy.com/2014/04/draw-error-ellipse-representing-covariance-matrix/
# https://maitra.public.iastate.edu/stat501/lectures/MultivariateNormalDistribution-I.pdf
# https://en.wikipedia.org/wiki/Multivariate_normal_distribution
# https://en.wikipedia.org/wiki/Ellipse#General_parametric_form
from scipy.stats import chi2, norm
def ellipse(n_sim, x, y, sigma):
data = np.concatenate((x, y), axis = 1)
L, V = np.linalg.eig(np.cov(data.T, ddof = 1))
idx = L.argsort()[::-1]
L = L[idx]
V = V[:, idx]
c = np.sqrt(chi2.ppf(2 * norm.cdf(sigma) - 1, 2))
t = np.linspace(0, 2 * np.pi, n_sim)
phi = np.arctan2(V[1, 0], V[0, 0])
a = c * np.sqrt(L[0]) * np.cos(t)
b = c * np.sqrt(L[1]) * np.sin(t)
R = np.array([[np.cos(phi), np.sin(phi)], [-np.sin(phi), np.cos(phi)]])
r = np.matmul(np.column_stack([a, b]), R)
result = np.add(r, np.mean(data, axis = 0)) # 2D only
return resultreturns_x_df = returns_df.dropna()[factors] # extended history
returns_x_mat = np.asarray(returns_x_df)
ellipse_x_mat = ellipse(1000, returns_x_mat[:, [0]], returns_x_mat[:, [2]], 1)
Granger causality
Granger causality tests whether lagged values of one series improve predictions of another series. Regress each series on lagged values of itself and the other series, then use a Wald test on the coefficient of the other series, i.e. the restriction \(R\hat{\beta}=r\) (note: condition on all own lags up to the momentum order but test a single cross lag at the momentum order):
\[ \begin{aligned} \left(R\hat{\beta}-r\right)^\mathrm{T}\left(R\hat{V}R^\mathrm{T}\right)^{-1}\left(R\hat{\beta}-r\right)\xrightarrow{d}\chi_{Q}^{2} \end{aligned} \]
- https://github.com/cran/lmtest/blob/master/R/waldtest.R
- https://en.wikipedia.org/wiki/Wald_test#Test(s)_on_multiple_parameters
- https://math.stackexchange.com/a/1591946
def granger_test(x, y, order):
# compute lagged observations
lag_x = x.shift(order).rename("lag_x")
lag_y = pd.concat([y.shift(i) for i in range(1, order + 1)], axis = 1)
lag_y.columns = ["lag_y" + str(i) for i in range(1, order + 1)]
# collect series
df = pd.concat([y.rename("y"), lag_y, lag_x], axis = 1)
z = sm.add_constant(df.drop("y", axis = 1))
# fit full model
fit = sm.OLS(df["y"], z, missing = "drop").fit()
R = np.append(np.zeros(order + 1), 1)
coef = fit.params
r = 0 # technically a matrix (see Stack Exchange)
matmul = np.dot(R, coef) - r
matmul_mid = np.linalg.inv(np.atleast_2d(np.dot(R, np.dot(fit.cov_params(), R.T))))
wald = np.dot(matmul.T, np.dot(matmul_mid, matmul))
result = 1 - chi2.cdf(wald, 1)
return resultdef roll_lead_lag(x, y, width, order, sig_level):
n_rows = len(x)
x_name = x.name
y_name = y.name
x_y_ls = []
y_x_ls = []
for i in range(width - 1, n_rows):
idx = range(max(i - width + 1, 0), i + 1)
x_y = granger_test(x.iloc[idx], y.iloc[idx], order)
y_x = granger_test(y.iloc[idx], x.iloc[idx], order)
x_y_status = (x_y < sig_level) and (y_x > sig_level)
y_x_status = (x_y > sig_level) and (y_x < sig_level)
x_y_ls.append(x_y_status)
y_x_ls.append(y_x_status)
result = pd.DataFrame({x_name: x_y_ls, y_name: y_x_ls}, index = x.index[(width - 1):])
return resultsig_level = 0.05score_x_df = score_df.loc[:, "SP500"]
score_y_df = score_df.loc[:, "DGS10"]lead_lag_df = roll_lead_lag(score_x_df, score_y_df, width, order, sig_level)