[R] Wavelet Bootstrap Size Simulation

PudgeDavid PudgeDavid at gmail.com
Tue Feb 24 23:03:55 CET 2009


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.




More information about the R-help mailing list