# Confidence Intervals for Mean Difference between Actuals and Forecast

I am trying to assess the influence of a product change compared to a hypothetical forecast. I use a time series model to predict a realistic future and compare it to the observed actuals. I am interested in the mean absolute and percentual difference between the prediction and the actuals.

How may I correctly estimate confidence intervals in this scenario and account for the specific time series structure of the data?

library(forecast)
library(ggplot2)

train   <- subset(AirPassengers, end = 84)
actuals <- subset(AirPassengers, start = 85)

m1 <- ets(train, "AAN")
fc <- forecast(m1, h = 12 * 5)

fc %>%
autoplot() +
autolayer(actuals) # Current Approach: T-Test for Confidence Intervals
mean(actuals - fc$mean) #>  60.60125 t.test(actuals, fc$mean)  # possibly use paired = TRUE
#>
#>  Welch Two Sample t-test
#>
#> data:  actuals and fc$mean #> t = 5.3794, df = 78.309, p-value = 7.533e-07 #> alternative hypothesis: true difference in means is not equal to 0 #> 95 percent confidence interval: #> 38.17509 83.02742 #> sample estimates: #> mean of x mean of y #> 396.4333 335.8321 uplift_pct <- (actuals - fc$mean) / fc$mean mean(uplift_pct) #>  0.1767661 t.test(uplift_pct) #> #> One Sample t-test #> #> data: uplift_pct #> t = 7.4578, df = 59, p-value = 4.535e-10 #> alternative hypothesis: true mean is not equal to 0 #> 95 percent confidence interval: #> 0.1293381 0.2241941 #> sample estimates: #> mean of x #> 0.1767661  • You put an "N" in for the seasonal component, so naturally your prediction doesn't have one. Since virtually all the time series structure is due to the predictable annual variation, you're missing out on it. – jbowman Jan 14 '18 at 22:45 • I actually produced the "bad" forecast deliberately to show a clear distinction between the forecast and the actuals, because I am mostly interested in learning about the confidence intervals for their mean differences. (I probably should have thought of a more realistic example, but this was the first, that came to my mind) – stats-hb Jan 15 '18 at 19:23 • Oh... so the problem is that you don't know how to get the pointwise confidence intervals from the forecast object? – jbowman Jan 15 '18 at 19:25 • not exactly - I do not know, if my t-test-approach for the confidence intervals is valdid in a time series context (with potential autocorrelations in the data) and when one of the time series is not "manifes" but an estimation itself (with its own uncertainty) – stats-hb Jan 15 '18 at 19:32 ## 1 Answer This is easily done using the general structure you've provided. Unfortunately, you left out the annual seasonal term, and consequently the algorithm was unable to estimate the seasonal component that contains most of the time series structure. Note also that for this data the series seems to increase in a nonlinear fashion, as does the variability due to seasonal swings. This indicates that a Box-Cox transform may help. The following code does what you want, with a little extra added for the Box-Cox transform. The confidence intervals, as plotted, do not contain the additional variability due to the estimation of the Box-Cox parameter$\lambda$, which, as it happens, equals$0$- i.e., a log transform performed best according to a minimum AIC criterion. library(forecast) library(ggplot2) train <- subset(AirPassengers, end = 84) actuals <- subset(AirPassengers, start = 85) aic <- ets(train, "ZZZ", lambda=0, biasadj=TRUE)$aic
min_lambda <- 0
for (lambda in seq(0.05, 1, by=0.05)) {
m2 <- ets(train, "ZZZ", lambda=lambda, biasadj=TRUE)
if (m2$aic < aic) { aic <- m2$aic
min_lambda <- lambda
}
}
m2 <- ets(train, "ZZZ", lambda=min_lambda, biasadj=TRUE)
fc2 <- forecast(m2, h=60)
fc2 %>% autoplot() + autolayer(actuals)


which leads to the following graph: Clearly this is not a very good forecast at the longer horizons, but in real life we'd probably forecast over a shorter horizon to start with and use out-of-sample predictions to select $\lambda$. In any case, the objective is to show you how to capture the annual structure and the confidence intervals of the forecasts in the model.

Now, on to something more like the question you really asked! The forecast object fc2 contains, as a default, 80% prediction intervals and 95% prediction intervals in the first and second columns respectively of fc2$lower and fc2$upper. We can extract these intervals from the forecast object as follows:

> ci <- data.frame(cbind(lower_95=fc2$lower[,2], upper_95=fc2$upper[,2]))
lower_95 upper_95
Jan 1956 263.7051 301.8168
Feb 1956 264.3147 307.9777
Mar 1956 302.9796 358.7300
Apr 1956 290.3735 348.8713
May 1956 285.0566 347.1562
Jun 1956 319.1490 393.6367


These are pointwise prediction intervals, in the sense that for each period taken individually there is an $\alpha$% probability of falling outside the interval (given all the assumptions.) However, forecast errors are autocorrelated, so the probability of falling outside the interval in, say, Feb. 1956 is not independent of whether or not we fell outside the interval in Jan. 1956; consequently, a t-test based on the errors over the whole forecast horizon would be biased due to the failure of independence.

Since forecast errors tend to be positively correlated for forecast models that take into account short-term dependencies in the variate of interest, we can apply Hochberg's step-up procedure (see here) for a joint test.

ci <- data.frame(cbind(obsno=1:60, lower_95=fc2$lower[,2], upper_95=fc2$upper[,2]))
ci$p_value <- 1-pnorm(abs((actuals-fc2$mean))/((ci$upper-ci$lower)/(2*1.96)))
ci <- ci[order(ci$p_value),] ci$hochberg_cutoff <- 0.05/(61-1:60)


with result:

> head(ci)
obsno lower_95 upper_95    p_value  hochberg_cutoff
51    51 421.3037 713.9750 0.03590013     0.0008333333
38    38 332.3203 525.6013 0.05502230     0.0008474576
28    28 336.5221 500.7919 0.06149721     0.0008620690
27    27 349.7245 517.0156 0.06288973     0.0008771930
36    36 323.5269 505.8235 0.06568278     0.0008928571
26    26 303.6777 445.9398 0.07688071     0.0009090909


which in this instance would not reject any time period as having fallen outside the Hochberg-corrected 95% test interval, as the minimum observed p-value is $>$ the Hochberg cutoff level.

There are plenty of other multiple comparison corrections, see for example the Wikipedia page, that may be more or less useful to you.

• Nice example showing how to choose an optimal lambda for the transformation. I don"t accept your answer right now though, because I"m looking for the confidence intervals for the difference between actuals and forecast. – stats-hb Jan 15 '18 at 19:27
• No worries, I misunderstood your sticking point. I'll update the answer with something hopefully more on target! – jbowman Jan 15 '18 at 19:29
• Thank you very much for your update. The way I understand your approach is, that you look you compute pointwise p_values based on the prediction interval on the forecast (really elegant!) and then correct for multiple testing. This is basically testing for any significant deviation between actuals and forecast. In my use case, I need to quantify the mean and/or the cumulative effect of an intervention. The logic is very much comparable to the CausalImpact-package from Google. This is why I am so interested in the confidence intervals for the mean effect for a given period (say 4 years or so). – stats-hb Jan 18 '18 at 23:20
• It sounds like what you'd ideally do in your real application is have an external regressor mixed in with the ETS model to capture the intervention effect and be able to test its significance in a more straightforward manner? – jbowman Jan 18 '18 at 23:24
• Possibly, yes. I look for a way to summarise the differences (with confidence intervals). I actually use the prophet-package to generate the forecast and compare this to the actuals. My main point is, that a t-test on the differences (act - fc) for a certain time period does not account for autocorrelation and the uncertainty of the forecast and might potentially be biased. – stats-hb Jan 19 '18 at 8:23