Combinando duas séries temporais calculando a média dos pontos de dados

Gostaria de combinar o previsto e o backcast (ou seja, os valores passados previstos) de um conjunto de dados de séries temporais em uma série temporal, minimizando o Erro de Previsão Quadrada Média.

Digamos que eu tenha séries temporais de 2001 a 2010 com uma lacuna para o ano de 2007. Consegui prever 2007 usando os dados de 2001-2007 (linha vermelha - chamada ) e fazer backcast usando os dados de 2008-2009 (azul claro linha - chame ). $Y_f$ $Y_b$

Gostaria de combinar os pontos de dados de e em um ponto de dados imputado Y_i para cada mês. Idealmente, eu gostaria de obter o peso modo a minimizar o Erro Médio de Previsão Quadrática (MSPE) de . Se isso não for possível, como eu encontraria a média entre os pontos de dados das duas séries temporais? $Y_f$ $Y_b$ $w$ $Y_i$

Y_{i} = w \cdot Y_{f} + (1 - w) \cdot Y_{b}

$Y_i = w\cdot Y_f + (1-w)\cdot Y_b$

Como um exemplo rápido:

tt_f <- ts(1:12, start = 2007, freq = 12)
tt_b <- ts(10:21, start=2007, freq=12)

tt_f
     Jan Feb Mar Apr May Jun Jul Aug Sep Oct Nov Dec
2007   1   2   3   4   5   6   7   8   9  10  11  12
tt_b
     Jan Feb Mar Apr May Jun Jul Aug Sep Oct Nov Dec
2007  10  11  12  13  14  15  16  17  18  19  20  21

Gostaria de obter (apenas mostrando a média ... Idealmente, minimizando o MSPE)

tt_i
     Jan Feb Mar Apr May Jun  Jul  Aug  Sep  Oct  Nov  Dec
2007 5.5 6.5 7.5 8.5 9.5 10.5 11.5 12.5 13.5 14.5 15.5 16.5

insira a descrição da imagem aqui

r time-series forecasting data-imputation OSlOlSO
fonte

Qual é o modelo de previsão (arima, ets, alguns outros)? (+1) para a sugestão de abordagem, pensei uma vez sobre esse assunto, mas permaneci dentro de Expectativa-Maximização após interpolação. Em princípio, o período de aprendizado pode ter importância, para dar um peso maior ao modelo com base em informações maiores (previsão vermelha na figura). Alguns critérios de precisão também podem ser potencialmente úteis para fazer pesos, para não serem tão deterministicamente vinculados aos comprimentos das séries temporais.

Dmitrij Celov 08/08

Desculpe por deixar de fora o modelo de previsão. O exemplo acima é simplesmente usando a predictfunção do pacote de previsão. No entanto, acho que vou usar o modelo de previsão da HoltWinters para prever e fazer backcast. Tenho séries temporais com pouco menos de 50 contagens e tentei a previsão de regressão de Poisson - mas, por algum motivo, com previsões muito fracas.

OSlOlSO 08/08

Os dados para contagens parecem ter uma quebra exatamente no local que você mostra, as previsões e retrocessos também ilustram a mesma coisa. Em Poisson, você fez a regressão de na tendência temporal ?

\log (c o u n t s)

$\log(counts)$

t

$t$

Dmitrij Celov 08/08

Você tem apenas contagens ou algumas séries temporais relacionadas adicionais sem NAvalores? Parece que tornar o MSPE no período de aprendizado pode ser enganoso, uma vez que os subperíodos são bem descritos por tendências lineares, mas no período perdido ocorre uma queda em algum lugar, e na verdade poderia ser qualquer ponto. Observe também que, como as previsões são de tendência colinear, sua média introduzirá duas quebras estruturais em vez de aparentemente uma.

Dmitrij Celov

Desculpe por voltar agora @Dmitij. Do que você está falando sobre esse 'intervalo'? Eu fiz o log (contagens) para a regressão GLM. E há um subconjunto dos dados de contagem que têm contagens menores que <6 que me forçarão a usá-los. Eu só tenho as contagens. Se você olhar para esta pergunta , terá uma idéia dos dados que tenho. As contagens acima são apenas para o grupo de idade '15up'. Se isso faz sentido?

OSlOlSO

Respostas:

Supondo que você tenha os Erros de previsão ao quadrado da previsão e da backcast individualmente, eu recomendaria o seguinte: Seja w um vetor de comprimento 12, seja o mês em que você está interessado.

w=rep(NA,12);
for(w in 1:12){
w[m]=SPE_Backcast[m]/(SPE_Backcast[m]+SPE_Forecast[m]);
}

Agora w é o peso para a previsão e 1-w é o peso para o backcast.

Dennis Jaheruddin
fonte

Isso parece ponderar mais alto o valor mais baixo (a ponto de os números negativos poderem ter pesos> 1). Qual é o objetivo? Além disso, linha doiss/w/m/

naught101

Como você obteria erros de previsão ao quadrado negativos?

Owe Jessen

Seu objetivo é executar uma suavização de intervalo fixo (FI) da série temporal. O valor suavizado da observação no momento é definido como uma expectativa condicional onde a notação é para vector das observações de tempo para o tempo . Acima, supõe-se que o intervalo seja o intervalo que varia do tempo a , e é o comprimento de toda a série. O tempo está na lacuna e a expectativa pode ser escrita $t$

{\hat{Y}}_{t} := E (Y_{t} | Y_{1 : r}, Y_{s : n})

$\widehat{Y}_{t} := \mathbb{E}(Y_t|\mathbf{Y}_{1:r},\,\mathbf{Y}_{s:n})$

Y_{u : v} := [Y_{u}, Y_{u + 1}, \dots, Y_{v}]

$\mathbf{Y}_{u:v} := [Y_u,\,Y_{u+1}, \, \dots,\,Y_v]$

u

$u$

v

$v$

r + 1

$r+1$

s - 1

$s-1$

n

$n$

t

$t$

{\hat{Y}}_{t | 1 : r, s : n}

$\widehat{Y}_{t|1:r, s:n}$ para recordar sua natureza condicional.

O valor suavizado não possui a forma simples que você adivinha. Para uma série temporal estacionária gaussiana com estrutura de covariância conhecida, a estimativa de para na lacuna pode ser encontrada resolvendo um sistema linear. $\widehat{Y}_{t}$ $t$

Quando o modelo de série temporal pode ser colocado no formato State Space (SS), a suavização FI é uma operação padrão baseada na filtragem Kalman e pode ser feita, por exemplo, usando as funções R disponíveis. Você só precisa especificar que os valores na lacuna estão ausentes. O algoritmo de suavização estima o estado oculto que contém todas as informações relevantes sobre para na lacuna. Os modelos ARIMA podem ser colocados no formato SS. $\boldsymbol{\alpha}_t$ $Y_t$ $t$

Curiosamente, a suavização de FI pode ser escrita como uma combinação de dois filtros: um para a frente e outro para trás, levando a uma fórmula do tipo que você esperava, mas para a estimativa do estado oculto (previsão e backcast), mas não para a observação . Isso é conhecido como filtragem Rauch-Tung-Striebel . $\boldsymbol{\alpha}_t$ $Y_t$

Pelo menos nas versões multiplicativas, os procedimentos de previsão 'ad hoc', como Holt-Winters, se baseiam em modelos estocásticos, sem algoritmos FI simples, uma vez que não podem ser colocados na forma SS. A fórmula de suavização provavelmente pode ser aproximada usando o modelo SS, mas é muito mais simples usar modelos de séries temporais estruturais com transformações de log. As funções 'KalmanSmooth', 'tsSmooth' e 'StructTS' do pacote de estatísticas R podem fazer o trabalho. Você deve dar uma olhada nos livros de Harvey ou Durbin e Koopman citados nas páginas de ajuda da R. O algoritmo de suavização pode fornecer uma variação condicional para o estimado $Y_t$ e pode ser usado para criar intervalos de suavização, que geralmente tendem a ser maiores no meio do intervalo. Observe, no entanto, que a estimativa de modelos estruturais pode ser difícil.

AP <- log10(AirPassengers) 
## Fit a Basic Structural Model
fit <- StructTS(AP, type = "BSM")

## Fit with a gap
AP.gap <- AP
AP.gap[73:96] <- NA
fit.gap <- StructTS(AP.gap, type = "BSM", optim.control = list(trace = TRUE))

# plot in orginal (non-logged) scale
plot(AirPassengers, col = "black", ylab = "AirPass")
AP.missing <- ts(AirPassengers[73:96], start=1955, , freq=12)
lines(AP.missing, col = "grey", lwd = 1)

## smooth and sum 'level' and 'sea' to retrieve series
sm <- tsSmooth(fit.gap)
fill <- apply(as.matrix(sm[ , c(1,3)]), 1, sum)
AP.fill <- ts(fill[73:96], start=1955, , freq=12)
lines(10^AP.fill, col = "red", lwd = 1)

Preenchimento suavizado

Yves
fonte

Acho interessante a sua abordagem sugerida, de usar os meios dos modelos dianteiro e traseiro.

Vale ressaltar que, em qualquer sistema que exibe estrutura caótica, é provável que as previsões sejam mais precisas em períodos mais curtos. Esse não é o caso para todos os sistemas, por exemplo, um pêndulo amortecido pode ser modelado por uma função com o período errado; nesse caso, todas as previsões de médio prazo provavelmente estarão erradas, enquanto as de longo prazo serão todas muito preciso, pois o sistema converge para zero. Mas parece-me, pelo gráfico da pergunta, que essa pode ser uma suposição razoável a ser feita aqui.

Isso implica que é melhor confiarmos mais nos dados de previsão da parte anterior do período ausente e mais nos dados de conversão para a última parte. A maneira mais simples de fazer isso seria usar um peso linearmente decrescente para a previsão e o oposto para o back-cast:

> n <- [number of missing datapoints] 
> w <- seq(1, 0, by = -1/(n+1))[2:(n+1)]

Isso dá um pouco de peso ao backcast no primeiro elemento. Você também pode usar n-1, sem os subscritos no final, se quiser usar apenas o valor da previsão no primeiro ponto interpolado.

> w
 [1] 0.92307692 0.84615385 0.76923077 0.69230769 0.61538462 0.53846154
 [7] 0.46153846 0.38461538 0.30769231 0.23076923 0.15384615 0.07692308

Não tenho seus dados, então vamos tentar isso no conjunto de dados AirPassenger em R. Vou remover um período de dois anos perto do centro:

> APearly <- ts(AirPassengers[1:72], start=1949, freq=12)
> APlate <- ts(AirPassengers[97:144], start=1957, freq=12)
> APmissing <- ts(AirPassengers[73:96], start=1955, freq=12)
> plot(AirPassengers)
# plot the "missing data" for comparison
> lines(APmissing, col="#eeeeee")
# use the HoltWinters algorithm to predict the mean:
> APforecast <- hw(APearly)[2]$mean
> lines(APforecast, col="red")
# HoltWinters doesn't appear to do backcasting, so reverse the ts, forecast, 
# and reverse again (feel free to edit if there's a better process)
> backwards <- ts(rev(APlate), freq=12)
> backcast <- hw(backwards)[2]$mean
> APbackcast <- ts(rev(backcast), start=1955, freq=12)
> lines(APbackcast, col='blue')
# now the magic: 
> n <- 24 
> w <- seq(1, 0, by=-1/(n+1))[2:(n+1)]
> interpolation = APforecast * w + (1 - w) * APbackcast
> lines(interpolation, col='purple', lwd=2)

E aí está a sua interpolação.

saída gráfica

Claro, não é perfeito. Eu acho que isso é resultado dos padrões na parte anterior dos dados serem diferentes dos da última parte (o pico de julho a agosto não é tão forte nos anos anteriores). Mas, como você pode ver na imagem, é claramente melhor do que apenas a previsão ou a transmissão inversa. Imagino que seus dados possam obter resultados um pouco menos confiáveis, pois não há uma variação sazonal tão forte.

Meu palpite seria que você também poderia tentar isso, incluindo os intervalos de confiança, mas não tenho certeza da validade de fazê-lo da maneira mais simples possível.

naught101
fonte