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