Markets

analysis
finance
r
Author

Jason Foster

Published

August 9, 2026

factors_r <- c("SP500", "DTWEXAFEGS") # "SP500" does not contain dividends; note: "DTWEXM" discontinued as of Jan 2020
factors_d <- c("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_xts <- na.omit(lag(exp(roll::roll_sum(returns_xts, width - order, min_obs = 1)) - 1, order))

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} \]

# volatility scale only
score_xts <- na.omit(momentum_xts / roll::roll_sd(momentum_xts, width, center = FALSE, min_obs = 1))
overall_xts <- xts::xts(rowMeans(score_xts), zoo::index(score_xts))
overall_xts <- overall_xts / roll::roll_sd(overall_xts, width, center = FALSE, min_obs = 1)
colnames(overall_xts) <- "Overall"
score_xts <- na.omit(merge(overall_xts, score_xts))

Outlier detection

Interquartile 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\):

outliers <- function(z) {
  
  n_cols <- ncol(z)
  result_ls <- list()
  
  for (j in 1:n_cols) {
    
    y <- z[ , j]
    
    if (n_cols == 1) {
      x <- 1:length(y)
    } else {
      x <- cbind(1:length(y), z[ , -j])
    }
    
    coef <- coef(lm(y ~ x))
    predict <- coef[1] + x %*% as.matrix(coef[-1])
    resid <- y - predict
    
    lower <- quantile(resid, prob = 0.25)
    upper <- quantile(resid, prob = 0.75)
    iqr <- upper - lower
    
    total <- y[(resid < lower - 1.5 * iqr) | (resid > upper + 1.5 * iqr)]
    
    result_ls <- append(result_ls, list(total))
    
  }
  
  result <- do.call(merge, result_ls)
  
  return(result)
  
}
outliers_xts <- outliers(score_xts[ , 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
ellipse <- function(n_sim, x, y, sigma) {

    data <- cbind(x, y)
    LV <- eigen(cov(data))
    L <- LV$values
    V <- LV$vectors

    c <- sqrt(qchisq(2 * pnorm(sigma) - 1, 2))
    t <- seq(0, 2 * pi, len = n_sim)
    phi <- atan2(V[2, 1], V[1, 1])
    a <- c * sqrt(L[1]) * cos(t)
    b <- c * sqrt(L[2]) * sin(t)
    R <- matrix(c(cos(phi), -sin(phi), sin(phi), cos(phi)), nrow = 2, ncol = 2)
    r <- t(rbind(a, b)) %*% R

    result <- sweep(r, 2, colMeans(data), "+") # 2D only

    return(result)

}
returns_x_xts <- na.omit(returns_xts)[ , factors] # extended history
ellipse_x_xts <- ellipse(1000, returns_x_xts[ , 1], returns_x_xts[ , 3], 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} \]

granger_test <- function(x, y, order) {

  # compute lagged observations
  lag_x <- lag(x, order)
  lag_y <- lag(y, 1:order)

  # collect series
  data <- merge(y, lag_y, lag_x)
  colnames(data) <- c("y", paste0("lag_y", 1:order), "lag_x")

  # fit full model
  fit <- lm(y ~ ., data = data)

  R <- matrix(c(rep(0, order + 1), 1), nrow = 1)
  coef <- fit$coefficients
  r <- 0 # technically a matrix (see Stack Exchange)

  wald <- t(R %*% coef - r) %*% solve(R %*% vcov(fit) %*% t(R)) %*% (R %*% coef - r)

  result <- 1 - pchisq(wald, 1)

  return(result)

}
roll_lead_lag <- function(x, y, width, order, sig_level) {
    
  n_rows <- nrow(x)
  x_name <- names(x)
  y_name <- names(y)
  x_y_ls <- list()
  y_x_ls <- list()
  
  for (i in width:n_rows) {
    
    idx <- max(i - width + 1, 1):i
    
    x_y <- granger_test(x[idx], y[idx], order)
    y_x <- granger_test(y[idx], x[idx], order)
    
    x_y_status <- (x_y < sig_level) && (y_x > sig_level)
    y_x_status <- (x_y > sig_level) && (y_x < sig_level)
    
    x_y_ls <- append(x_y_ls, list(x_y_status))
    y_x_ls <- append(y_x_ls, list(y_x_status))
    
  }
  
  result <- data.frame(do.call(c, x_y_ls), do.call(c, y_x_ls))
  result <- xts::xts(result, zoo::index(x)[width:n_rows])
  colnames(result) <- c(x_name, y_name)
  
  return(result)
    
}
sig_level <- 0.05
score_x_xts <- score_xts[ , "SP500"]
score_y_xts <- score_xts[ , "DGS10"]
lead_lag_xts <- roll_lead_lag(score_x_xts, score_y_xts, width, order, sig_level)