Como ajustar uma função de etapa robusta a uma série temporal?

7

Tenho uma série temporal um tanto barulhenta que paira em diferentes níveis.

Por exemplo, os seguintes dados:

insira a descrição da imagem aqui

Eu tenho os dados de linha sólida disponíveis e gostaria de obter uma estimativa para a linha tracejada. Deve ser constante por partes.

Quais algoritmos são apropriados para testar aqui?

Até agora, minhas idéias pairam em torno de splines P de 0 grau (mas como descobrir onde colocar os nós?) Ou em modelos de quebra estrutural. Uma árvore de regressão é a melhor ideia que tenho atualmente, mas, idealmente, eu estaria procurando um método que levasse em conta o fato de que os dois níveis em y = 250 estão com valores y iguais. Se bem entendi, uma árvore de regressão dividiria esses dois intervalos em dois grupos diferentes, cada um com uma média diferente.

O código R que o gerou é este:

set.seed(20181118)
true_fct = stepfun(c(100, 200, 250), c(200, 250, 300, 250))
x = 1:400
y = true_fct(x) + rt(length(x), df=1)
plot(x, y, type="l")
lines(x, true_fct(x), lty=2, lwd=3)
Alexander Engelhardt
fonte
2
Se seus dados realmente se parecerem com os simulados, dificilmente você poderá fazer melhor do que calcular uma mediana em janela com uma janela muito pequena: isso detectaria de maneira confiável todos os saltos. Estime os níveis usando as medianas das respostas dentro de cada intervalo detectado. Você poderia, portanto, indicar se as suposições implícitas da simulação - grandes saltos, medianas constantes e erros de Student t - são precisamente as suposições que devemos fazer?
whuber
11
Obrigado por seu comentário! Tenho duas observações: (1) Como obteria os intervalos da mediana da janela? (2) As suposições são medianas constantes e saltos perceptíveis, mas não sei nada sobre a distribuição de erros, além do fato de que grandes valores discrepantes podem acontecer.
Alexander Engelhardt
Às vezes, métodos não paramétricos simples funcionam quando o problema é simples. Gostaria que você simulasse um conjunto de dados mais desafiador / realista, onde há uma estrutura de arima incorporada e talvez um pulso sazonal ou dois. Abordagens abrangentes para problemas como esse precisam considerar e isolar estruturas e anomalias autoregressivas durante o processamento. Você pode postar outra pergunta e incluir o conjunto de dados um pouco mais realista.
IrishStat
Gostaria também de acrescentar que o nível / mudanças de passo são tão grandes vis-a-vis o processo de erro métodos não paramétricos podem desempenhar um papel útil e menos de modo que a relação fica menor
IrishStat

Respostas:

7

Um método simples e robusto para lidar com esse ruído é calcular medianas.

Uma mediana contínua sobre uma janela curta detectará todos os saltos, exceto os menores, enquanto as medianas da resposta dentro de intervalos entre os saltos detectados estimarão de forma robusta seus níveis. (Você pode substituir esta última estimativa por qualquer estimativa robusta que não seja afetada pelos valores discrepantes.)

Você deve ajustar essa abordagem com dados reais ou simulados para obter taxas de erro aceitáveis. Por exemplo, para a simulação na pergunta, achei bom usar o segundo e o 98º percentil para definir limites para detectar os saltos. Em outras circunstâncias - como quando muitos saltos podem ocorrer - mais percentis centrais funcionariam melhor.

Aqui está o resultado mostrando (a) os três saltos como pontos vermelhos e (b) os quatro níveis estimados como linhas azuis claras.

Figura

Estima-se que os saltos ocorram nos índices 100, 200, 250 (que é exatamente onde a simulação os faz ocorrer) e os níveis resultantes são estimados em 199,6, 249,8, 300,0 e 250,2: tudo dentro de 0,4 dos valores subjacentes verdadeiros.

Esse excelente comportamento persiste com repetidas simulações (removendo o set.seed comando no início).

Aqui está o Rcódigo.

#
# Rolling medians.
#
rollmed <- function(x, k=3) {
  n <- length(x)
  x.med <- sapply(1:(n-k+10), function(i) median(x[i + 0:(k-1)]))
  l <- floor(k/2)
  c(rep(NA, l), x.med, rep(NA, k-l))
}
y.med <- rollmed(y, k=5)
#
# Changepoint analysis.
#
dy <- diff(y.med)
fourths <- quantile(dy, c(1,49)/50, na.rm=TRUE)
thresholds <- fourths + diff(fourths)*2.5*c(-1,1)
jumps <- which(dy < thresholds[1] | dy > thresholds[2]) + 1

points(jumps, y.med[jumps], pch=21, bg="Red")
#
# Plotting.
#
limits <- c(1, jumps, length(y)+1)
y.hat <- rep(NA, length(jumps)+1)
for (i in 1:(length(jumps)+1)) {
  j0 <- limits[i]
  j1 <- limits[i+1]-1
  y.hat[i] <- median(y[j0:j1])
  lines(x[j0:j1], rep(y.hat[i], j1-j0+1), col="skyblue", lwd=2)
}
whuber
fonte
+1, mas a parte "análise do ponto de mudança" do código pode não estar totalmente clara para alguns usuários. Talvez você possa comentar o que está acontecendo lá?
Tim
@ Tim Obrigado pela sugestão. O objetivo do primeiro parágrafo é explicar esse algoritmo. Gostaria de subestimar os detalhes de sua implementação porque eles não são importantes: basta aplicar qualquer método robusto de triagem externa aos resíduos.
whuber
Você pode considerar zoo::rollmedianuma função semelhante para simplificar seu código.
usεr11852
@ usεr11852 Obrigado. Estou ciente,zoo mas eleito para não usá-lo, porque sou preguiçoso! Era mais rápido e fácil escrever rollmeddo que revisar as chamadas de argumento para qualquer função que já estivesse disponível. Além disso, gosto de rollmedilustrar com clareza o que estou fazendo, em vez de ocultar os detalhes atrás de uma caixa preta.
whuber
Sem problemas. :) (eu estava certo de que sabia de zoo, eu era incerto se você não usá-lo por opção ou por acidente Boa resposta, em qualquer caso +1.)
usεr11852
3

Se você ainda estiver interessado em suavizar as penalidades L0, daria uma olhada na seguinte referência: "Visualização de alterações genômicas por suavização segmentada usando uma penalidade L0" - DOI: 10.1371 / journal.pone.0038230 (uma boa introdução à O Whittaker mais suave pode ser encontrado no artigo P. Eilers "Um perfeito mais suave" - ​​DOI: 10.1021 / ac034173t). Obviamente, para atingir seu objetivo, você precisa trabalhar um pouco em torno do método.

Em princípio, você precisa de 3 ingredientes:

  1. O mais suave - eu usaria o mais suave de Whittaker. Além disso, usarei o aumento da matriz (ver Eilers e Marx, 1996 - "Suavização flexível com splines B e penalidades", p.101).
  2. Regressão quantílica - usarei o pacote R quantreg (rho = 0,5) para preguiça :-)
  3. Penalidade L0 - Seguirei a mencionada "Visualização de alterações genômicas por suavização segmentada usando uma penalidade L0" - DOI: 10.1371 / journal.pone.0038230

Obviamente, você precisaria também de uma maneira de selecionar a quantidade ideal de suavização. Isso é feito pelos meus olhos de carpinteiro neste exemplo. Você pode usar os critérios em DOI: 10.1371 / journal.pone.0038230 (pág. 5, mas eu não tentei no seu exemplo).

Você encontrará um pequeno código abaixo. Deixei alguns comentários como guia.

# Cross Validated example
rm(list = ls()); graphics.off(); cat("\014")

library(splines)
library(Matrix)
library(quantreg)

# The data
set.seed(20181118)
n = 400
x = 1:n
true_fct = stepfun(c(100, 200, 250), c(200, 250, 300, 250))
y = true_fct(x) + rt(length(x), df = 1)

# Prepare bases - Identity matrix (Whittaker)
# Can be changed for B-splines
B = diag(1, n, n)

# Prepare penalty - lambda parameter fix
nb = ncol(B)
D = diff(diag(1, nb, nb), diff = 1)
lambda = 1e2

# Solve standard Whittaker - for initial values
a = solve(t(B) %*% B + crossprod(D), t(B) %*% y, tol = 1e-50)    

# est. loop with L0-Diff penalty as in DOI: 10.1371/journal.pone.0038230
p = 1e-6
nit = 100
beta = 1e-5

for (it in 1:nit) {
  ao = a

  # Penalty weights
  w = (c(D %*% a) ^ 2  + beta ^ 2) ^ ((p - 2)/2)
  W = diag(c(w))

  # Matrix augmentation
  cD = lambda * sqrt(W) %*% D
  Bp = rbind(B, cD)
  yp =  c(y, 1:nrow(cD)*0)

  # Update coefficients - rq.fit from quantreg
  a = rq.fit(Bp, yp, tau = 0.5)$coef

  # Check convergence and update
  da = max(abs((a - ao)/ao))
  cat(it, da, '\n')
  if (da < 1e-6) break
}

# Fit 
v = B %*% a

# Show results
plot(x, y, pch = 16, cex = 0.5)
lines(x, y, col = 8, lwd = 0.5)
lines(x, v, col = 'blue', lwd = 2)
lines(x, true_fct(x), col = 'red', lty = 2, lwd = 2)
legend("topright", legend = c("True Signal", "Smoothed signal"), 
       col = c("red", "blue"), lty = c(2, 1))

insira a descrição da imagem aqui PS. Esta é minha primeira resposta no Cross Validated. Espero que seja útil e claro o suficiente :-)

Gi_F.
fonte
1

Eu consideraria o uso de Outliers de papel de Ruey Tsay , mudanças de nível e mudanças de variação no modelo de diferenciação de séries temporais com os outliers AR1 e 21.

insira a descrição da imagem aqui

insira a descrição da imagem aqui

Desativamos a diferença e as mudanças de nível são especificamente mencionadas.

insira a descrição da imagem aqui

Tom Reilly
fonte
11
Gostaria de saber se você ignorou a ênfase em "robusto" no título da pergunta, porque qualquer método que identifique 18 parâmetros espúrios (correspondentes aos valores discrepantes introduzidos na simulação) além dos 3 saltos reais dificilmente pode ser considerado robusto (ou parcimonioso, por esse assunto).
whuber
Essa é uma solução robusta. Não sei por que você é contra a identificação e o ajuste de valores discrepantes, mas há um mundo de pesquisas apoiando isso e, é claro, nossas experiências. Essas outras variáveis ​​são discrepantes. Eu adicionei um gráfico que mostra os dados históricos e uma versão limpa para contrastar a diferença.
Tom Reilly
11
Você poderia ser explícito sobre qual é a sua estimativa da função step?
whuber
11
Há um sinalizador no período 100 (x3), 200 (x2), 250 (x4) que mostra a etapa. O operador de diferenciação torna um pouco mais difícil de ver, mas o efeito é o mesmo. Eu adicionei um modelo sem diferenciar.
Tom Reilly