How is a prediction interval generated for many input observations?

Question

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")

Answer 1

Ok, so as long as I can can see what you are trying to do it's not good. In the first place you are trying to do a forecast on observed data, this is wrong. If you have observed data you can do a prediction interval for its mean, and its variance. Doing a prediction interval for observed data makes no sense. As you can see in the first graph your plot is just the plot of the observed data and the "Prediction interval" is just a bound defined by its standard deviation times two.

In the second place you have not done a time series model. In order to do so, you have to define your dependent variable in a recursive way, remember that a Time Series model is:

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

Talking about the prediction intervals, you don't need new variables because your dependent variable is the only one. Brief, the prediction intervals depend upon the coefficients of your model, and its variance.

I conseil you to read more about time series before trying to do a model. You will have a bad time if you don't do it.