Result of a diagnostic test of a predictive model looking too good

Question

I have created a predictive model that outputs a predictive density. I used 1000 rolling windows to estimate the model and predict one step ahead in each window. I collected the 1000 predictions and compared them to the actual realizations. I used several diagnostic tests, among them Kolmogorov-Smirnov. I saved the $p$-value of the test.

I did the same for multiple time series. Then I looked at all of the $p$-values from the different series. I found that they are 0.440, 0.579, 0.848, 0.476, 0.753, 0.955, 0.919, 0.498, 0.997. At first I was quite happy that they are much larger than 0.010, 0.050 or 0.100 (to use the standard cut-off values). But then a colleague of mine pointed out that the $p$-values should be distributed as $\text{Uniform}[0,1]$ under the null of correct predictive distribution, and so I should perhaps not be so happy.

On the one hand, the colleague must be right; the $p$-values should ideally be uniformly distributed. On the other hand, I have found that my model predicts "better" than the true model normally would; the discrepancy between the predicted density and the realized density is less than one would normally expect between the true density and the realized density. This could be an indication of overfitting if I were evaluating my model in-sample, but the model has been evaluated out of sample. What does this tell me? Should I be concerned with a diagnostic test's $p$-values being too high?

You could say this is just a small set of $p$-values (just 8 of them) so anything could happen, and you might be right. However, suppose I have a larger set of $p$-values that are closer to 1 than uniformly distributed; is that a problem? What does that tell me?

---

Update: I have adapted the code by @jbowman and played with different values of $k$ and $n$ to arrive at the following histograms of $p$-values from the Kolmogorov-Smirnov test. The relationship between the shape of the histogram and $k$ and $n$ seems nontrivial.

• For small $n$ (columns on the left), increasing $k$ (going down the rows) appears to push the $p$-values towards $0$.
• For large $n$ (columns on the right), increasing $k$ (going down the rows) appears to push the $p$-values towards $1$.
• For cases in between such as $n=30$, I do not see a clear trend.

(The number of simulation runs is $1000$ for each histogram.)

You can experiment for yourself by running the code below:

m=1e3 # number of simulation runs
n=1e3 # size of training sample (length of rolling window)
k=1e3 # size of test sample / number of rolling windows within one simulation run

ks_value <- rep(NA, m)
p_value  <- rep(NA, k)

# Runs 45  seconds for n=1e3, m=1e3
# Runs 5.5 minutes for n=1e4, m=1e3
print(paste0(Sys.time()," start"))
for (j in 1:m) {
   if(j%%100==0) print(paste0(Sys.time()," simulation run no. ",j," of ",m))
   set.seed(j); x <- rnorm(n+k)
   for (i in 1:k) {
      xbar       <- mean(x[i:(n+i-1)])   
      s_d        <- sd  (x[i:(n+i-1)])
      p_value[i] <- pnorm(x[n+i], xbar, s_d)
   }
   ks_value[j]   <- ks.test(p_value, punif)$p.value
}; print(paste0(Sys.time()," end"))

hist(ks_value)

It is still a puzzle for me why this is happening. I can see that the $p$-values of observations are autocorrelated across windows (within a simulation run) because of the overlap of the training windows that causes the estimated mean and estimated standard deviation to be autocorrelated. However, I do not see why this should skew their unconditional distribution towards $1$. Actually, deliberately setting the estimated mean and estimated standard deviation to be fixed across windows (within a simulation run) yields KS test's $p$-values that look fairly uniform when collected across the simulation runs. This appears to imply that the problem is caused by data leakage, i.e. test observations becoming training observations as the window rolls.
But why/how exactly does that cause nonuniformity?

Answer 1

That "rolling windows" in your second sentence raises the question of whether the 1000 predictive densities are actually independent of each other... with obvious implications for the validity of the diagnostic tests. Let's perform an experiment, which of course will not match your problem exactly.

We will generate independent 1100 standard normal variates, then compute a rolling sample mean and standard deviation based on windows of size 100. This will give us 1000 estimates of the distributions of $x_{i+1}$ for $i = 101, \dots, 1100$. We then calculate the predictive cumulative density function $p_{i+1}$ of each of the $x_{i+1}$, and compare the resulting collection of 1000 of these with a $U(0,1)$ distribution using the Kolmogorov-Smirnov test. (Given our sample size of 100, the use of a Normal distribution for predictive purposes seems reasonable.) We repeat the entire exercise 100 times and plot a histogram of the resulting K-S test p-values.

ks_value <- rep(0, 100)
p_value <- rep(0, 1000)

for (j in 1:100) {
   x <- rnorm(1100)
   for (i in 1:1000) {
      xbar <- mean(x[i:(i+99)])   
      s_d <- sd(x[i:(i+99)])
      p_value[i] <- pnorm(x[i+100], xbar, s_d)
   }
   ks_value[j] <- ks.test(p_value, punif)$p.value
}

with resultant histogram:

which certainly seems to match the pattern of your observed p-values better than a Uniform distribution!

Changing our code allows us to verify that in the independent sample case we get a reasonable collection of p-values from the K-S test:

for (j in 1:100) {
   for (i in 1:1000) {
      x <- rnorm(101)
      xbar <- mean(x[1:100])   
      s_d <- sd(x[1:100])
      p_value[i] <- pnorm(x[101], xbar, s_d)
   }
   ks_value[j] <- ks.test(p_value, punif)$p.value
}

which, given only 100 p-values, does not appear unreasonably far from a Uniform distribution.

One mechanism that would cause this result is as follows. Consider a given $x_{i+1}$ and the associated predictive distribution. If $x_{i+1}$ is "large", its p-value will also be large. As the rolling window moves forward, though, that "large" value will cause the sample mean and standard deviation to, all other things being equal, be slightly larger than they would have been had $x_{i+1}$ been "small". This will reduce the probability of obtaining a large p-value from the comparison between $x_{i+2}$ and its associated predictive distribution, working in a way similar to that of antithetic variables in simulation. Similarly, if $x_{i+1}$ is very near the sample mean (relative to the size of the standard deviation), the sample standard deviation will become smaller for the comparison between $x_{i+2}$ and its associated distribution, increasing the chances for a more extreme p-value. Consequently, the distribution of p-values from the predictive distribution tests would seem to be more uniform than they ought to be given their sample size, although with a large enough number of replications no doubt their non-uniformity would manifest itself.

We can see these effects in the following plots. We generate 10 standard normal variates, with the 11th being fixed at either 2.0 or 0.2. We calculate the p-value of the 11th observation and the 12th observation (given a window width of 10) and repeat 100,000 times. We then plot the histograms of the four sets of p-values.

Our code:

p_values <- matrix(0, nrow=100000, ncol=2)

for (i in 1:100000) {
   x <- c(rnorm(10), 2, rnorm(1))

   p_values[i,1] <- pt((x[11] - mean(x[1:10])) / (sqrt(1 + 1/10) * sd(x[1:10])), df=9)
   p_values[i,2] <- pt((x[12] - mean(x[2:11])) / (sqrt(1 + 1/10) * sd(x[2:11])), df=9)
}

hist(p_values[,1], main="Histogram of p-values for x[11] = 2", xlab = "p-values")
hist(p_values[,2], main="Histogram of p-values for x[12] | x[11] = 2", xlab = "p-values")

First, the case where $x_{11} = 2.0$:

And again, for the case where $x_{11} = 0.2$:

It's also important to note that the effects in the second of each of the two pairs of plots persist for the entire length of the window.