# How is a prediction interval generated for many input observations?

I am aware that the "prediction interval", as defined in most textbooks on linear models, is focused on the uncertainty in the model being fit and is used to estimate an output prediction range for exact inputs. However, as is normally the case, what about situations where the input isn't exact? How do you calculate a "prediction interval" that accounts for both the uncertainty in the model and the uncertainty in the input data?

For example,

Assume the model is: y = a0 + (a1 * x1) + (a2 * x2)

Where y, x1, and x2 are time series vectors.

I only have one observation of the y time series with its associated x1 and x2 time series. That data is used to fit the model. However, I also have 1000 additional observations for x1 and x2. I can easily calculate individual "prediction intervals", using the model and each x1 x2 pair, however I want to estimate a "prediction interval" that allows for all x1 and x2 observations. To be more specific, I generated the example below. The questionable code is below the #====== comment line.

Basically, I used lm(..) to fit the y vector to its associated x1 and x2 vectors. Next, I used predict.lm(...interval="prediction", level=alpha) to generate its typical fit, lwr, and upr vectors for each additional x1 x2 pair (see Graph 1, fit is the solid line, upr and lwr are the dashed lines for two observations). I then collected these vectors for all x1 x2 pairs and used alpha to extract the upper and lower tails of all of the lwr and upr vectors (see Graph 2).

Is this scheme valid? Does the alpha that was used in predict.lm(... level=alpha) apply directly to counting-up/sorting the results to generate a "prediction interval" that allows for all x1 x2 pairs? Can a 5% range from a competing model (for example an ARIMA model) be compared to this 5% "prediction interval"?

I'm fairly sure that the following scheme isn't right, but so far, I haven't figured out what I need to fix.

set.seed(1)

numpoi <- 10 #Number of data points in a time series vector

#First independent variable "x1", first observation
x1mea <- 0.03
x1sta <- 0.05
x1 <- cumsum(rnorm(numpoi, mean=x1mea, sd=x1sta))

#Second independent variable "x2", first observation
x2mea <- -0.01
x2sta <- 0.1
x2 <- cumsum(rnorm(numpoi, mean=x2mea, sd=x2sta))

#Dependent variable "y", first observation
a0 <- 3
a1 <- 2
a2 <- 1
noimea <- 0.0
noista <- 0.1
y <- a0 + (a1 * x1) + (a2 * x2) + rnorm(numpoi, mean=noimea, sd=noista)

#Build a data frame of the "first observation" data
datfra <- data.frame(y=y, x1=x1, x2=x2)

#Fit the model for the "first observation"
mod <- lm(y ~ x1 + x2, data=datfra)
summary(mod)

#Set up desired alpha value
alpha <- 0.95
onetai <- (1 - alpha)/2 #Convert the two tail "alpha" to a one tail value for use later

#Generate some new data "a" for a second observation of "x1" and "x2".
x1a <- cumsum(rnorm(numpoi, mean=x1mea, sd=x1sta))
x2a <- cumsum(rnorm(numpoi, mean=x2mea, sd=x2sta))
datfraa <- data.frame(y=rep(NA, numpoi), x1=x1a, x2=x2a)
modprea <- predict(mod, newdata=datfraa, interval="prediction", level=alpha)

#Generate some new data "b" for a third observation of "x1" and "x2".
x1b <- cumsum(rnorm(numpoi, mean=x1mea, sd=x1sta))
x2b <- cumsum(rnorm(numpoi, mean=x2mea, sd=x2sta))
datfrab <- data.frame(y=rep(NA, numpoi), x1=x1b, x2=x2b)
modpreb <- predict(mod, newdata=datfrab, interval="prediction", level=alpha)

#Plot the results for both new data "a" and "b"
plot(modprea[, 1], type="l", ylim=c(min(modprea, modpreb), max(modprea, modpreb)), main="Graph 1 - Second and Third Observations for x1 and x2", lwd=2, col="red")
lines(modprea[, 2], lwd=2, lty=2, col="red")
lines(modprea[, 3], lwd=2, lty=2, col="red")
lines(modpreb[, 1], lwd=2,  col="green")
lines(modpreb[, 2], lwd=2, lty=2, col="green")
lines(modpreb[, 3], lwd=2, lty=2, col="green")

#===========================================================================
#The code below is where my question lies.    Is this the appropriate method
#to account for all observations?

#Run the above calculation scheme on "all" observations.
numtri <- 1000 #All observations
modprecfit <- matrix(0, nrow=numpoi, ncol = numtri) #Matrix to hold all "fit" vectors
modpreclwr <- matrix(0, nrow=numpoi, ncol = numtri) #Matrix to hold all "lwr" vectors
modprecupr <- matrix(0, nrow=numpoi, ncol = numtri) #Matrix to hold all "upr" vectors

#Fill up the modprecfit, modpreclwr, and modprecupr matricies.
for (i in 1:numtri) {

#Generate the new data and put it in a dataframe
x1c <- cumsum(rnorm(numpoi, mean=x1mea, sd=x1sta))
x2c <- cumsum(rnorm(numpoi, mean=x2mea, sd=x2sta))
datfrac <- data.frame(y=rep(NA, numpoi), x1=x1c, x2=x2c)

#Predict the new "y" for new input data
modprec <- predict(mod, newdata=datfrac, interval="prediction", level=alpha)

#Store the "ft", "lwr", and "upr" vectors so they can be processed later
modprecfit[, i] <- modprec[, 1]
modpreclwr[, i] <- modprec[, 2]
modprecupr[, i] <- modprec[, 3]

}

#Extract the average "fit", the lower quantile "lwr", and the upper quantile "upr"
modprecfitfin <- apply(modprecfit, 1, quantile, 0.5)
modpreclwrfin <- apply(modpreclwr, 1, quantile, onetai)
modprecuprfin <- apply(modprecupr, 1, quantile, (1 - onetai))

plot(modprecfitfin, type="l", ylim=c(min(modpreclwrfin), max(modprecuprfin)), main="Graph 2 - All Observations for x1 and x2", lwd=2, col="red")
lines(modpreclwrfin, lwd=2, lty=2, col="red")
lines(modprecuprfin, lwd=2, lty=2, col="red")  $$Y_t = \phi_0 + Y_{t-1}+Y_{t-2}+\epsilon_t$$
In consequence, you can't have three variables, because when you are doing a time series model, all you need is one. And, you can't compute the model through lm, you have the arima function, which is the one for this kind of situations.