I am comparing prediction intervals from linear regression and ARIMA for a simple AR(1) model:
p = lag(p)
The models were built on monthly data from 2003-2013 years. Predictions were made for 2014 year.
As for linear regression I am making dynamic forecasting - I mean that I don't use known values from 2014 year, instead for future predictions I use previously predicted values.
Comparing prediction intervals from this two models I noticed how different they are.
library(dplyr)
library(forecast)
data = data.frame(year = rep(seq(2003,2014), each = 12), mon = rep(1:12, 12), p = c(2.36, 2.13, 1.82, 1.50, 1.83, 1.52, 1.66, 2.17, 1.72, 1.21, 1.92, 1.82, 1.40, 2.35, 2.54, 2.52, 2.31, 2.97, 3.05, 3.05, 3.55, 3.38, 3.01, 3.06, 2.88, 2.63, 3.27, 3.48, 3.71, 3.62, 5.17, 5.51, 6.00, 6.35, 5.94, 5.89, 6.01, 5.64, 5.78, 6.03, 6.02, 5.73, 4.28, 3.78, 2.75, 2.89, 3.52, 3.68, 3.09, 2.52, 2.00, 1.78, 1.85, 1.81, 1.62, 1.14, 2.10, 2.65, 3.16, 3.16, 4.34, 5.35, 5.33, 6.00, 7.57, 8.67, 9.17, 6.52, 6.11, 3.99, 2.10, 0.46, -0.53, -0.18, -0.29, -0.99, -3.29, -3.96, -4.32, -1.06, -0.91, 0.39, 1.94, 3.46, 3.98, 3.55, 3.33, 2.96, 3.51, 3.44, 3.61, 3.32, 3.00, 2.91, 2.75, 2.98, 2.90, 2.78, 3.04, 4.11, 4.32, 4.21, 4.20, 4.33, 4.02, 4.17, 4.14, 3.42, 3.27, 3.28, 3.41, 2.53, 2.60, 2.66, 2.81, 2.76, 3.38, 3.27, 2.66, 3.55, 3.31, 3.20, 2.68, 2.49, 2.29, 2.30, 2.03, 1.63, 1.40, 1.41, 1.88, 1.66, 1.93, 1.95, 2.12, 2.48, 2.62, 2.35, 2.16, 2.12, 1.73, 1.45, 1.24, 0.62))
data$lag_p = lag(data$p)
data_mod = data[data$year %in% 2003:2013, ]
data_fore = data[data$year %in% 2014, ]
model1 = lm(p ~ lag_p, data = data_mod)
model2 = Arima(ts(data_mod$p, frequency = 12, start = c(2003, 1)), order = c(1,0,0))
# Arima forecast
arima_fore = forecast(model2, h = 12)
# Linear regression forecast
data_fore_copy = data_fore
dyn_fore = vector()
lwr = vector(); upr = vector()
for ( i in 1:nrow(data_fore_copy) ) {
data_dyn = data_fore_copy[i, ]
y_hat = predict(model1, newdata = data_dyn, interval = "prediction", level = 0.95)
dyn_fore = c(dyn_fore, y_hat[1, 1])
lwr = c(lwr, y_hat[1, 2])
upr = c(upr, y_hat[1, 3])
if ( i < nrow(data_fore_copy) ) {
data_fore_copy[i + 1, "lag_p"] = y_hat[1, 1]
}
}
my_dates = seq(as.Date("2013-12-01"), as.Date("2014-12-01"), length.out = 13)
actual_values = data_fore$p
plot_type = "Dynamic Forecast"
pre_actual = data_mod[nrow(data_mod), "p"]
plot(my_dates, c(NA, dyn_fore), type = "l", col = "red", cex.main = 0.8, ylab = "Values", xlab = "Dates",
ylim = c(-2, 8),
main = "Prediction intervals. Arima vs Linear regression")
lines(my_dates, c(pre_actual, data_fore$p), col = "blue")
lines(my_dates, c(NA, lwr), col = "magenta")
lines(my_dates, c(NA, upr), col = "magenta")
lines(my_dates, c(NA, arima_fore$mean), col = "green")
lines(my_dates, c(NA, arima_fore$lower[,2]), col = "orange")
lines(my_dates, c(NA, arima_fore$upper[,2]), col = "orange")
legend("topleft", lty = c(1, 1, 1, 1, 1), col = c("red", "green", "blue", "magenta", "orange"), legend = c("Forecast. Linear regression", "Forecast. Arima", "Actual", "Prediction interval: Linear regression", "Prediction interval: Arima"))
My question is: I guess that prediction interval from linear regression is so different from ARIMA (whereas predictions match almost exactly) because linear regression constructs prediction intervals as if it makes STATIC predictions but I actually make DYNAMIC predictions. So, what can I do in order to get prediction interval from linear regression which will be like prediction interval from ARIMA in case of DYNAMIC forecasting? How can I transform it (are there any generalizable formulas)?
I mean that if I, for example, standing in today make predictions 12 periods ahead DYNAMICLY, then uncertainty in my predictions should INCREASE with time like on the picture below:
Two videos on youtube about constructing prediction intervals: