Identifying a changed spending pattern

Question

I have a list of numbers and I am trying to detect a drastic point difference from the previous numbers and see if the pattern has changed. For instance,

x = [5000,5500,6250,4800,3950,7200,5500,800,1200,900,500,400,300,200]

Above, there is high spending until 800 and then it seems that there is high spending before the 800 and low spending after the 800. All in all after 800, the spending has decreased a good bit. I want to try and divide the list based on this drastic point and then check if there is a different pattern (i.e. there is high spending before and low spending after or if there is low spending before and low spending after or low spending before high spending after or high spending before and high spending after). Essentially, I am looking for this sort of inflection point and trying to detect if there are two different classes of numbers. I know I could set a threshold and check each number or assume a normal distribution and check standard deviations. Is there a better way to approach this problem from a more statistical point of view? Note that the numbers could be of larger scale, purchases could be tens of thousands of dollars and then drop to thousands of dollars and this should indicate a change. I'm wondering if there is some statistical based method that is better than using standard deviations and means.

For a set of numbers like x above the following does not work well

def reject_outliers(data, m = 2.):
    d = np.abs(data - np.median(data))
    mdev = np.median(d)
    s = d/mdev if mdev else 0.
    return data[s<m]

Answer 1

You describe a circumstance in which spending might vary around a constant level (under the null hypothesis) or change suddenly to vary around a second level at some point during the series. One way to measure this is to compare the variance in spending under the two models. In the first case it would be the mean squared difference between spending and the overall average. In the second case, the differences would be relative to the two levels of spending for the best possible model: that is, among all possible points where the series could break, choose the one where the magnitude of the effect is largest. Let's call the resulting mean squared difference the "two-level variance."

A permutation test can do a fine job of assessing the data. Its null hypothesis is that the sequence of the data doesn't matter. The sampling distribution of the two-level variance therefore is obtained by examining all possible reorderings of the data and computing its value for each one. In practice there are so many reorderings that we sample from them randomly, thereby approximating the sampling distribution. The data will show significant evidence of a break when it is rare for the reordered data to exhibit a two-level variance that is as low (or lower) than the two-level variance of the data themselves.

Here, in the histogram at left, is an approximate permutation distribution of the two-level variance for the data given in the question. The vertical dotted red line marks the value of the two-level variance for the actual data. The p-value is very low because no permuted versions of the data exhibited a smaller two-level variance.

The plot at the right shows the (obvious) break where the symbols change between indexes 7 and 8.

Here is R code to reproduce (or extend) these calculations. It requires about one second to compute ten thousand values of the permutation distribution--more than enough for practical work. (The computation of the test statistic can be optimized and made to scale directly with the amount of data.)

#
# The two-level variance.
#
stat <- function(y) {
  r2 <- function(i) {
    r1 <- function(z) crossprod(z - mean(z))
    (r1(y[1:i]) + r1(y[-(1:i)])) / length(y)
  }
  r <- sapply(1:length(y), r2)
  i <- which.min(r)
  return(list(Value=r[i], Index=i))
}
#
# The data.
#
y <- c(5000,5500,6250,4800,3950,7200,5500,800,1200,900,500,400,300,200)
#
# Approximate the permutation distribution of the two-level variance.
#
set.seed(17)
sim <- replicate(1e4, stat(sample(y))$Value)
#
# Compute the two-level variance of the data, then compute its p-value.
#
stat.data <- stat(y)
sim <- c(sim, stat.data$Value)
p.value <- mean(sim <= stat.data$Value)
#
# Plot the results.
#
par(mfrow=c(1,2))
hist(sim, xlab="Statistic", main="Permutation Distribution")
abline(v=stat.data$Value, col="red", lwd=2, lty=3)
plot(y, main="Estimated Break", sub=paste("p =", signif(p.value, 3)))
points(y[1:stat.data$Index], pch=16)
par(mfrow=c(1,1))

Answer 2

You could use a sequential t-test as described at the link below and implemented in Visual Basic for Application (Excel).

http://www.beringclimate.noaa.gov/regimes/

Answer 3

You could consider k-means clustering.

# Data
x = c(5000,5500,6250,4800,3950,7200,5500,800,1200,900,500,400,300,200)

# Calculate within-group sum of squares for different number of clusters
n <- length(x)
wss <- rep(0, 6)
wss[1] <- (n - 1)*var(x)
for (i in 2:6) {
  wss[i] <- sum(kmeans(x, centers = i)$withinss)
}

# Plot WGSS vs number of groups
plot(1:6,
     wss,
     type = "b",
     xlab = "Number of groups",
     ylab = "Within groups sum of squares")

# Perform k-means clustering
clusters <- kmeans(x, centers = 2, iter.max = 1e4, nstart = 1e2)

# Create results data frame
xClustered <- data_frame(i = 1:length(x), x = x) %>% 
  mutate(cluster = factor(clusters$cluster))

# Plot clusters
ggplot(data = xClustered, aes(x = i, y = x, color = cluster, shape = cluster)) +
  geom_point() +
  theme(axis.text.x = element_text(angle = 90, hjust = 1)) +
  labs(x = "index", y = "x", color = "cluster", fill = "cluster")

Answer 4

Look to the work of Tsay and Balke on the idea of a "level shift". It is a deterministic variable that looks like this 0,0,0,0,0,0,0,1,1,1,1,1,1,etc.

Here are results using Autobox (I am affiliated with this company). An outlier was detected at period 6 and a level shift down 4,635 at period 8.

Y(T) = 5250.0 test
+[X1(T)][(- 4635.7 )] :LEVEL SHIFT 8 +[X2(T)][(+ 1950.0 )] :PULSE 6 + + [A(T)]

When running the tso() package it also confirms these results.