One-tailed F-test with one restriction

beta <- c(1,4,2,3)       # The coefficients
n <- 20                  # The number of observations
combo <- c(0,1,1,-1,-1)  # The contrast (starting with an intercept coefficient)
sigma <- 2               # The error SD
n.sim <- 1e3             # Number of iterations of the simulation
#
# Prepare to simulate.
#
set.seed(17)
k <- length(beta)
combo.name <- paste0("Contrast(", paste(combo, collapse=","), ")")
#
# Simulate with many different sets of regressors and responses, but all
# using the same `beta`.
#
p.values <- replicate(n.sim, {
  #
  # Create regressors.
  #
  x <- matrix(rnorm(n*k), n)
  colnames(x) <- paste0("X", 1:k)
  y <- x %*% beta + rnorm(n, 0, sigma)
  #
  # Test the combination.
  #
  fit <- lm(y ~ ., as.data.frame(x))
  beta.hat <- crossprod(coef(fit), combo)
  beta.hat.se <- sqrt(combo %*% vcov(fit) %*% combo)
  t.stat <- beta.hat / beta.hat.se
  p <- pt(t.stat, n-k-1) # A one-sided p-value
  #
  # When running this "for real," uncomment the following lines to
  # see the results.
  #
#   stats <- rbind(coef(summary(fit)), c(beta.hat, beta.hat.se, t.stat, p))
#   rownames(stats)[k+2] <- combo.name
#   print(stats, digits=3)
  #
  # Return the p-value for the (two-sided) t-test
  #
  p
})
#
# Display the simulation results.
#
n.bins <- floor(n.sim / 100)
(p.dist <- chisq.test(table(floor(n.bins*p.values)))$p.value) #$
hist(p.values, freq=FALSE, breaks=n.bins,
     sub=paste("p (uniform) =", format(p.dist, digits=3)))
abline(h = 1, col="Gray", lwd=2, lty=3)