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


 A: 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)


A: 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.
A: 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.
