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.
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. $\endgroup$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. $\endgroup$