I have written a short script to estimate the size of a test of non-constant
mean in an AR(1) time series. When I run the simulation on my PC (R version
2.7.1), I get the expected result: an empirical size much larger than the
nominal size. On the Red Hat machine (R version 2.7.2) in my department,
however, I get p-values > .45 for every simulated test. Below is my
simulation code (uses functions from the library waveslim - both computers
have the most current version). I'm trying to reconcile the difference
between these two sets of results, as the Linux machine is significantly
faster and I would much prefer to run my simulations there. Does anyone have
an idea of the source of the discrepancy? Thank you in advance for any help
you can provide, and please excuse the inelegant coding!
# These functions define the test.
fitphi <- function(data,x,sub){
y<-1:length(x)
y<-(y-.5)/length(y)
(1/length(data))*sum(data*cos(pi*sub*y))
}
mos <- function(data, x){
rdata <- rank(data)/(length(x)+1)
Rn <- rep(0, length(x))
coeff <- Rn
for (i in 1:length(x)){
coeff[i] <- fitphi(data,x,i)
}
for (j in 1:length(x)){
Rn[j]<-(24/j)*length(data)*sum((coeff[1:j])^2)
}
max(Rn)
}
#The function used to find the wavelet bootstrap replicate of the original
series
wavestrap <- function (y, wf = "d8", J = log(length(y), 2) - 1, infl =
sqrt(1.1)){
y.dwt <- dwt(y, wf, n.levels = J)
resample.dwt <- y.dwt
for (i in 1:(length(y.dwt)-1)) {
resample.dwt[[i]] <- sample(y.dwt[[i]]*infl, replace = TRUE)
}
idwt(resample.dwt)
}
#The simulation:
arcoeff<-.9 #AR(1) Coefficient
T<-128 #Realization Length
bsreps<-300 #Number of Bootstrap Replicates used for test
nosims<-20 #Number of Simulations
wspval<-rep(0,nosims) #initialize vector of p-values
#outer loop for recording calculated p-values
for (j in 1:nosims){
series<-arima.sim(T, model=list(ar=arcoeff)) #series would be the
observed
AR(1) series
test<-mos(series,1:length(series)) #observed value of test statistic
wsdist<-rep(0,bsreps) #initialize vector for wavelet bootstrap null
distribution
#innner loop for finding null distribution via wavelet bootstrap
for (i in 1:bsreps){
wsrep<-wavestrap(series, wf = "d4", infl=1)
wsdist[i]<-mos(wsrep,1:length(series))
}
wspval[j]<-sum(wsdist>=test)/length(wsdist) #find p-value (one-sided
test)
}
--
View this message in context:
http://www.nabble.com/Wavelet-Bootstrap-Size-Simulation-tp22191678p22191678.html
Sent from the R help mailing list archive at Nabble.com.
______________________________________________
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
