Combinando dos series de tiempo promediando los puntos de datos

Me gustaría combinar el pronóstico y el pronóstico (es decir, los valores pasados pronosticados) de un conjunto de datos de series de tiempo en una serie de tiempo minimizando el Error de predicción cuadrático medio.

Digamos que tengo series temporales de 2001-2010 con una brecha para el año 2007. He podido pronosticar 2007 usando los datos de 2001-2007 (línea roja - llamada ) y hacer una usando los datos de 2008-2009 (azul claro línea: ). $Y_f$ $Y_b$

Me gustaría combinar los puntos de datos de y en un punto de datos imputado Y_i para cada mes. Idealmente, me gustaría obtener el peso manera que minimice el error medio de predicción cuadrado (MSPE) de . Si esto no es posible, ¿cómo podría encontrar el promedio entre los dos puntos de datos de las series temporales? $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 un ejemplo 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

Me gustaría obtener (solo mostrando el promedio ... Idealmente minimizando el 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

ingrese la descripción de la imagen aquí

r time-series forecasting data-imputation OSlOlSO
fuente

¿Cuál es el modelo de pronóstico (arima, ets, algún otro)? (+1) para la sugerencia de enfoque, una vez pensé de esa manera, pero me mantuve dentro de Expectation-Maximization después de la interpolación. En principio, el período de aprendizaje podría ser importante, para dar mayor peso al modelo basado en información más grande (pronóstico rojo en la imagen). Algunos criterios de precisión también podrían ser potencialmente útiles para hacer pesas, ya que no están vinculados de manera tan determinista a las longitudes de series de tiempo.

Dmitrij Celov

Perdón por dejar de lado el modelo de pronóstico. El anterior es simplemente usar la predictfunción del paquete de pronóstico. Sin embargo, creo que voy a utilizar el modelo de predicción de HoltWinters para predecir y retroceder. Tengo series de tiempo con poco <50 recuentos, y probé el pronóstico de regresión de Poisson, pero por alguna razón a predicciones muy débiles.

OSlOlSO

Los datos para los recuentos parecen tener un descanso exactamente en el lugar que muestra, las previsiones y los retrocesos también ilustran lo mismo. En Poisson, ¿realizó la regresión de en la tendencia de tiempo ?

\log (c o u n t s)

$\log(counts)$

t

$t$

Dmitrij Celov

¿Tiene recuentos o algunas series temporales relacionadas adicionales sin NAvalores? Parece que hacer que el período de aprendizaje MSPE pueda ser engañoso ya que los subperíodos están bien descritos por tendencias lineales, pero en el período perdido se produce un desplegable en algún lugar, y en realidad podría ser cualquier punto. Tenga en cuenta también que, dado que las previsiones tienen una tendencia colineal, su promedio introducirá dos interrupciones estructurales en lugar de aparentemente una.

Dmitrij Celov

Lo siento por volver ahora @Dmitij. ¿De qué está hablando este "descanso"? Hice el registro (conteos) para la regresión GLM. Y hay un subconjunto de los datos de recuento que tienen recuentos inferiores a <6, lo que me obligará a usarlo. Solo tengo los recuentos. Si observa esta pregunta , tendrá una idea de los datos que tengo. Los recuentos anteriores son solo para el grupo de edad '15up'. Si esto tiene sentido?

OSlOlSO

Respuestas:

Suponiendo que tiene los Errores de predicción al cuadrado para el pronóstico y la retransmisión individualmente, recomendaría esto: Sea w un vector de longitud 12, sea el mes que le interese.

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

Ahora w es el peso para el pronóstico y 1-w es el peso para el backcast.

Dennis Jaheruddin
fuente

Esto parece ponderar más el valor más bajo (hasta el punto de que los números negativos pueden terminar teniendo pesos> 1). ¿Cuál es el punto de? Además, línea doss/w/m/

nada101

¿Cómo obtendrías errores de predicción al cuadrado negativo?

Le debo a Jessen el

Su propósito es realizar un suavizado de intervalo fijo (FI) de las series de tiempo. El valor suavizado de la observación en el tiempo se define como una expectativa condicional donde la notación es para vector de las observaciones de vez en cuando . Arriba, se supone que la brecha es el intervalo que va desde el tiempo hasta , es la longitud de toda la serie. El tiempo está en la brecha y la expectativa podría escribirse $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}$ recordar su naturaleza condicional.

El valor suavizado no tiene la forma simple que supones. Para una serie temporal estacionaria gaussiana con estructura de covarianza conocida, el estimado para en la brecha se puede encontrar resolviendo un sistema lineal. $\widehat{Y}_{t}$ $t$

Cuando el modelo de serie temporal se puede poner en forma de espacio de estado (SS), el suavizado FI es una operación estándar basada en el filtrado de Kalman y se puede hacer, por ejemplo, utilizando las funciones R disponibles. Simplemente necesita especificar que faltan los valores en el espacio. El algoritmo de suavizado estima el estado oculto que contiene toda la información relevante sobre para en la brecha. Los modelos ARIMA se pueden poner en forma SS. $\boldsymbol{\alpha}_t$ $Y_t$ $t$

Curiosamente, el suavizado de FI se puede escribir como una combinación de dos filtros: uno hacia adelante y otro hacia atrás, lo que lleva a una fórmula del tipo que esperaba, pero para la estimación de estado oculto (pronóstico y retroceso), pero no para la observación . Esto se conoce como filtrado Rauch-Tung-Striebel . $\boldsymbol{\alpha}_t$ $Y_t$

Al menos en las versiones multiplicativas, los procedimientos de pronóstico 'ad hoc' como Holt-Winters se basan en modelos estocásticos sin algoritmos FI simples, ya que no se pueden poner en forma SS. La fórmula de suavizado probablemente se puede aproximar usando el modelo SS, pero es mucho más simple usar modelos de series de tiempo estructurales con transformaciones de registros. Las funciones 'KalmanSmooth', 'tsSmooth' y 'StructTS' del paquete de estadísticas R pueden hacer el trabajo. Debería echar un vistazo a los libros de Harvey o de Durbin y Koopman citados en las páginas de ayuda de R. El algoritmo de suavizado puede proporcionar una varianza condicional para el estimado $Y_t$ y puede usarse para construir intervalos de suavizado, que generalmente tienden a ser más grandes en el medio del espacio. Sin embargo, tenga en cuenta que la estimación de los modelos estructurales puede 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)

Relleno alisado

Yves
fuente

Encuentro interesante su enfoque sugerido, de tomar los medios de los lanzamientos anteriores y posteriores.

Una cosa que vale la pena señalar es que en cualquier sistema que exhiba una estructura caótica, es probable que los pronósticos sean más precisos en períodos más cortos. Ese no es el caso para todos los sistemas, por ejemplo, un péndulo amortiguado podría ser modelado por una función con un período incorrecto, en cuyo caso es probable que todos los pronósticos a mediano plazo sean incorrectos, mientras que los a largo plazo serán todos muy preciso, ya que el sistema converge a cero. Pero me parece, a partir del gráfico en la pregunta, que esta podría ser una suposición razonable para hacer aquí.

Eso implica que podríamos estar mejor confiando más en los datos de pronóstico para la primera parte del período faltante, y más en los datos retrospectivos para la última parte. La forma más sencilla de hacer esto sería utilizar un peso que disminuya linealmente para el pronóstico, y lo contrario para el reverso:

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

Esto le da un poco de peso al backcast en el primer elemento También puede usar n-1, sin los subíndices al final, si desea usar solo el valor de pronóstico en el primer punto 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

No tengo sus datos, así que intentemos esto en el conjunto de datos de AirPassenger en R. Simplemente eliminaré un período de dos años cerca del 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)

Y ahí está tu interpolación.

salida gráfica

Por supuesto, no es perfecto. Supongo que es el resultado de que los patrones en la parte anterior de los datos son diferentes a los de la última parte (el pico de julio a agosto no es tan fuerte en los años anteriores). Pero, como puede ver en la imagen, es claramente mejor que solo el pronóstico o el casting posterior solo. Me imagino que sus datos pueden obtener resultados ligeramente menos confiables, ya que no existe una variación estacional tan fuerte.

Supongo que podría intentar esto, incluidos los intervalos de confianza también, pero no estoy seguro de la validez de hacerlo tan simple como esto.

nada101
fuente