Combining two time-series by averaging the data points

Question

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

Answer 1

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.

Answer 2

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)

Answer 3

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
 [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

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

And there's your interpolation.

Of course, it's not perfect. I guess that's a result of the patterns in the earlier part of the data being different to those in the latter part (the Jul-Aug peak is not so strong in earlier years). But as you can see from the image, it's clearly better than just the forecasting or the back casting alone. I would imagine that your data may get slightly less reliable results, as there is not such a strong seasonal variation.

My guess would be that you could try this including the confidence intervals too, but I'm not sure of the validity of doing it as simply as this.