d <- 150 ## The largest term in the expansion
N <- 100000 ## The number of simulations
k <- 3:d ## index
xx <- seq(-1,1,length=201) ## grid of x values from -1 to 1 at which to evaluate the process
stdev <- 1/(sqrt(gamma(k+1))) ## vector of standard deviations of Y/root(k!)
polymtrx <- outer(xx,k,"^") ## a matrix of the x's raised to the k power
scls <- exp(xx^2) - 1 - xx^2-xx^4/2 ## scaling factor
pck<- abs(xx)<0.03 ## for a grid mesh of .001, small values of xx can cause computation problems
scls[pck]<-xx[pck]^6/6+xx[pck]^8/24 ## for these values, calculate the scaling factor using the next two terms directly
zetamm <- matrix(rnorm(N*(d-2),mean=0,sd=rep(stdev,N)), nrow=d-2,ncol=N) ## a matrix of simulated normals
trac <- polymtrx %*% zetamm/sqrt(scls) ## matrix multiplication to do the proper sum
trac[xx==0,] <- 0 ## set the value of the sum at index of delta=0 equal to 0
mins <- apply(trac,2,min) ## finds the minimum of each column
rmin <- sort((abs(-zetamm[2,]*2*sqrt(6) - mins) - (-zetamm[2,]*2*sqrt(6) + mins))/2)^2
mht <-hist(mins^2,breaks=400,prob=T) ## draws a histogram of the minimums
quantile(rmin,p=0.95,type=6) ## the empirical critical value, adjusts for bias in quantile estimate
## written by Drew Carter and Doug Steigerwald, 2011
## user selections:
## d is the number of terms in the Taylor-series approximation and equals J-1 in the paper; J-1=150 is the recommended value
## N is the number of simulations; we ran 100,000
## xx is the grid of values for eta and corresponds to H in the paper, the length embeds the grid mesh, thus
## for a grid of values from (-1,1) and a grid mesh of .01, the length is 201 (as in the program above)
## for a grid of values from (-1,1) and a grid mesh of .001, the length is 2001
## for a grid mesh of .001, the values of xx near zero can lead to numeric accuracy problems in calculating scls
## because scls equals xx^6/6 + xx^8/24 + ...., we replace the scls calucation with the first two terms from this expansion for small values of xx