# Combining two time-series by averaging the data points

I would like to combine the forecasted and backcasted (viz. the predicted past values) of a time-series data set into one time-series by minimizing the Mean Squared Prediction Error.

Say I have time series from 2001-2010 with a gap for the year 2007. I have been able to forecast 2007 using the 2001-2007 data (red line - called it $Y_f$) and to backcast using the 2008-2009 data (light blue line - call it $Y_b$).

I would like to combine the data points of $Y_f$ and $Y_b$ into a imputed data point Y_i for each month. Ideally I would like to obtain the weight $w$ such that it minimizes the Mean Squared Prediction Error (MSPE) of $Y_i$. If this is not possible, how would I just find the average between the two time-series' data points?

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

As a quick example:

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


I would like to get (just showing the averaging... Ideally minimizing the 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 • What is the forecasting model (arima, ets, some other)? (+1) for the approach suggestion, I once thought about such way, but stayed within Expectation-Maximization after interpolation. In principle the learning period could matter, to give the higher weight for the model based on larger info (red forecast in picture). Some accuracy criteria could also be potentially useful to make weights, to be not so deterministically linked to the time series lengths. Aug 8, 2011 at 12:06
• Sorry about leaving out the forecasting model. The one above is simply using the predict function of the forecast package. However, I think that I'm going to use the HoltWinters forecasting model to predict and backcast. I have time series with little <50 counts, and tried Poisson regression forecasting - but for some reason to very weak predictions. Aug 8, 2011 at 14:36
• The data for counts seems to have a break exactly at the place you show, the forecasts and back-casts also illustrate the same thing. In Poisson you made regression of $\log(counts)$ on time trend $t$? Aug 8, 2011 at 14:46
• Do you have just counts or some additional related time series without NA values? It seems that making learning period MSPE could be misleading since the sub-periods' are well described by linear tendencies, but in the missed period a drop down somewhere occurs, and it actually could be any point. Note also that since the forecasts are collinear in trend, their average will introduce two structural breaks instead of seemingly one. Aug 9, 2011 at 8:59
• Sorry for only getting back now @Dmitij. What is this 'break' you are talking about? I did do the log(counts) for the GLM regression. And there is a subset of the count data that have counts less than <6 which will force me to use it. I only have the counts. If you look at this question you will get an idea of the data that I have. The above counts are only for the '15up' AgeGroup. If this makes sense? Aug 9, 2011 at 13:28

Assuming you have the Squared Prediction Errors for the forecast and backcast individually I would recommend this: Let w be a vector of length 12, let m be the month that you are interested in.

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


Now w is the weight for the forecast and 1-w is the weight for the backcast.

• This appears to just weight the lowest value more highly (to the point that negative numbers can end up having weights > 1). What's the point? Also, line two s/w/m/ Apr 5, 2012 at 1:04
• How would you get negative squared prediction errors? Apr 9, 2012 at 8:15

Your purpose is to perform a fixed interval (FI) smoothing of the time series. The smoothed value of the observation at time $t$ is defined as a conditional expectation $$\widehat{Y}_{t} := \mathbb{E}(Y_t|\mathbf{Y}_{1:r},\,\mathbf{Y}_{s:n})$$ where the notation $\mathbf{Y}_{u:v} := [Y_u,\,Y_{u+1}, \, \dots,\,Y_v]$ is for the vector of the observations from time $u$ to time $v$. Above, the gap is assumed to be the interval ranging from time $r+1$ to $s-1$, and $n$ is the whole series' length. The time $t$ is in the gap and the expectation could be written $\widehat{Y}_{t|1:r, s:n}$ to recall its conditional nature.

The smoothed value does not have the simple form you guess. For a gaussian stationnary time series with known covariance structure, the estimated $\widehat{Y}_{t}$ for $t$ in the gap can be found by solving a linear system.

When the time series model can be put in State Space (SS) form, the FI smoothing is a standard operation based on Kalman filtering and it can be done e.g. using available R functions. You simply need to specify that the values in the gap are missing. The smoothing algorithm estimates the hidden state $\boldsymbol{\alpha}_t$ which contains all the relevant information about $Y_t$ for $t$ in the gap. ARIMA models can be put in SS form.

Interestingly, the FI smoothing can be written as a combination of two filters: one forward and one backward, leading to a formula of the kind you expected, but for the hidden state estimation $\boldsymbol{\alpha}_t$ (forecast and backcast), but not for the observation $Y_t$. This is known as Rauch-Tung-Striebel filtering.

At least in the multiplicative versions, 'ad hoc' forecasting procedures like Holt-Winters rely on stochastic models with no simple FI algorithms since they can not be put in SS form. The smoothing formula can probably be approximated by using SS model, but it is much simpler then to use Structural Time Series models with log transformations. The 'KalmanSmooth', 'tsSmooth' and 'StructTS' functions of the R stats package can do the job. You should have a look at the books by Harvey or by Durbin and Koopman cited in the R help pages. The smoothing algorithm can provide a conditional variance for the estimated $Y_t$ and can be used to build smoothing intervals, which usually tend to be larger in the middle of the gap. Note however that the estimation of Stuctural Models can be difficult.

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) I find your suggested approach, of taking the means of the fore- and back-casts, interesting.

One thing that might be worth pointing out is that in any system exhibiting chaotic structure the forecasts are likely to be more accurate over shorter periods. That isn't the case for all systems, for example a damped pendulum could be modelled by a function with the wrong period, in which case all the medium term forecasts are likely to be wrong, while the long term ones are all going to be very accurate, as the system converges to zero. But it looks to me, from the graph in the question, that this might be a reasonable assumption to make here.

That implies that we might be better off relying more on the forecast data for the earlier part of the missing period, and more on the back-cast data for the latter part. The simplest way to do this would be to have using a linearly decreasing weight for the forecast, and the opposite for the back-cast:

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


This gives a little weight of the backcast on the first element. You could also use n-1, without the subscripts at the end, if you wanted to use only the forecast value on the first interpolated point.

> w
 0.92307692 0.84615385 0.76923077 0.69230769 0.61538462 0.53846154
 0.46153846 0.38461538 0.30769231 0.23076923 0.15384615 0.07692308


I don't have your data, so let's try this on the AirPassenger dataset in R. I'll just remove a two-year period near the centre:

> 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)$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)$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) 