################################################################## #simulated critical values for the limiting null distribution# ############################################################## set.seed(100) dd<-c(-.499,-(12:1)*0.04,0,(1:12)*0.04,.499) numd<-length(dd) result<-NULL n<-5000 reptim<-30000 stat<-rep(0,reptim) L0<-0.15 U0<-0.85 Ln0<-as.integer(n*L0) Un0<-as.integer(n*U0) sublength<-Un0-Ln0+1 for (h in 1:numd) { d1<-dd[h] for (j in 1:reptim) { data<-as.vector(simARMA0(n,H=d1+1/2)) substat<-rep(0,sublength) for (k in Ln0:Un0) { mean1<-mean(data[1:k]) mean2<-mean(data[(k+1):n]) inter1<-cumsum(data[1:k])-(1:k)*mean1 inter2<-cumsum(data[(k+1):n])-(1:(n-k))*mean2 substat[k-Ln0+1]<-sqrt(n)*abs(mean1-mean2)/(sqrt(sum(inter1^2)+sum(inter2^2))/n) } stat[j]<-max(substat) print(j) } result<-rbind(result,as.vector(quantile(stat,c(0.9,0.95,0.99)))) print(h) } ################################################# #long memory, reptim=10000 #n=5000, simARMA0 from longmemo package in R #d=-0.499, -(12:1)*0.04,0,(1:12)*0.04,.499 ################################################# > result h=1 [,1] [,2] [,3] [1,] 20.09091 21.99803 25.45926 [,1] [,2] [,3] [1,] 20.05159 21.92415 25.17966 [,1] [,2] [,3] [1,] 20.23206 22.04757 25.57672 [2,] 20.13137 22.14904 25.89856 [3,] 20.90437 23.09339 26.81124 [4,] 21.27919 23.71602 28.34925 [5,] 21.60397 23.93517 29.22708 [6,] 22.43455 25.07212 30.48307 [7,] 23.04724 25.84554 31.76960 [8,] 23.79224 26.91413 33.22689 [9,] 24.31893 27.81315 34.73245 [10,] 25.31452 29.06364 36.63289 [11,] 26.07800 29.93718 38.13267 [12,] 27.21747 31.36806 40.58119 [13,] 28.12480 32.22254 41.98237 [14,] 28.70614 33.42419 43.32156 [15,] 29.93047 34.52267 45.57731 [16,] 31.21991 36.46689 48.11292 [17,] 32.51348 38.05722 50.15767 [2,] 33.71634 39.24952 51.55109 [3,] 34.82833 40.43071 53.81783 [4,] 36.28219 42.84610 56.25739 [5,] 37.18398 44.04682 59.46327 [6,] 38.10746 45.33738 60.13217 [7,] 39.84096 47.17937 62.52489 [8,] 40.58492 48.16548 65.38060 [9,] 41.86039 49.65368 68.31884 [10,] 43.85135 52.92012 71.68123 [11,] 43.76838 51.58495 69.85238 h=26 [2,] 43.50805 51.91513 68.68221 h=27 [,1] [,2] [,3] [1,] 44.7115 53.50515 71.58486 [,1] [,2] [,3] [1,] 43.91084 52.57802 72.04208 ############################# #n=5000, reptim=10000 ############################# 20.05159 21.92415 25.17966 20.13137 22.14904 25.89856 20.90437 23.09339 26.81124 21.27919 23.71602 28.34925 21.60397 23.93517 29.22708 22.43455 25.07212 30.48307 23.04724 25.84554 31.76960 23.79224 26.91413 33.22689 24.31893 27.81315 34.73245 25.31452 29.06364 36.63289 26.07800 29.93718 38.13267 27.21747 31.36806 40.58119 28.12480 32.22254 41.98237 28.70614 33.42419 43.32156 29.93047 34.52267 45.57731 31.21991 36.46689 48.11292 32.51348 38.05722 50.15767 33.71634 39.24952 51.55109 34.82833 40.43071 53.81783 36.28219 42.84610 56.25739 37.18398 44.04682 59.46327 38.10746 45.33738 60.13217 39.84096 47.17937 62.52489 40.58492 48.16548 65.38060 41.86039 49.65368 68.31884 43.50805 51.91513 69.68221 44.71150 53.50515 71.58486 ############################ ############################ #local whittle estimation ############################ locwhittle<-function(d) { n<-length(ts) xv<-2*pi*(1:m)/n yv<-fft(ts[n:1]) yv<-yv[2:(m+1)] yv<-(Mod(yv))^2/2/pi/n obj1<-mean(xv^(2*d)*yv) obj2<-mean(log(xv))*2*d return(log(obj1)-obj2) } ############################################33 #not used here ############################################## Mmake.long<-function(d,ARts) #ARts is H number of time series of length n #[n,H] size, each column is a time series { L <- 10000 coef <- rep(0,L) coef[1] <- 1 for (k in 1:(L-1)) { coef[k+1] <- prod((d-1+(1:k))/(1:k)) } n<-dim(ARts)[1] H<-dim(ARts)[2] result<-matrix(0,n,H) for (k in 1:n) { result[k,1:H]<-matrix(coef[1:k],1,k)%*%ARts[k:1,1:H] } return(result[(n/2+1):n,1:H]) } ####################################### #AR(1) model with Gaussian innovations ####################################### MAR<-function(n,H,rr) { inter<-matrix(0,n+30,H) epsilon<-matrix(rnorm((n+30)*H,0,1),(n+30),H) for (j in 1:(n+29)) { inter[j+1,1:H]<-epsilon[j+1,1:H]+rr*inter[j,1:H] } return(inter[31:(n+30),1:H]) } ####################################### #MA(1) model with Gaussian innovations ####################################### MMA<-function(n,H,rr) { inter<-matrix(0,n+30,H) epsilon<-matrix(rnorm((n+31)*H,0,1),(n+31),H) for (j in 1:(n+30)) { inter[j,1:H]<-epsilon[j+1,1:H]+rr*epsilon[j,1:H] } return(inter[31:(n+30),1:H]) } ####################################3 #####################################3 #read the critical values (90%, 95%, 99%) cv<-read.table("cv3.txt") CV<-matrix(0,27,3) for (j in 1:3) { CV[,j]<-cv[,j] } ################################ #simulation setup (SIZE) ##################### n<-100 reptim<-5000 ################### #long memory ################### set.seed(100) d<-0.4 data<-matrix(0,n,reptim) for (j in 1:reptim) { data[,j]<-simARMA0(n,d+1/2) } ################# #AR(1) set.seed(100) data<-MAR(n,reptim,-0.5) ################# #################### #MA(1) set.seed(100) data<-MMA(n,reptim,0.5) #########################3 ################# #main program ################# dn<-NULL m<-as.integer(n^(2/3)) for (j in 1:reptim) { ts<-data[,j] opt3 <- nlminb(start=c(0),obj=locwhittle,lower=-0.499,upper=0.499) dn<-c(dn,opt3\$par) } critval<-matrix(0,2,reptim) seq1<-c(-0.499, -(12:1)*0.04,0,(1:12)*0.04,.499) for (j in 1:reptim) { index1<-which.max(as.numeric(dn[j]1) { #90% critval[1,j]<-((-dn[j]+seq1[index1])*CV[index1-1,1]+(dn[j]-seq1[index1-1])*CV[index1,1])/(seq1[index1]-seq1[index1-1]) #95% critval[2,j]<-((-dn[j]+seq1[index1])*CV[index1-1,2]+(dn[j]-seq1[index1-1])*CV[index1,2])/(seq1[index1]-seq1[index1-1]) } } ############### stat<-rep(0,reptim) L0<-0.15 U0<-0.85 Ln0<-as.integer(n*L0) Un0<-as.integer(n*U0) sublength<-Un0-Ln0+1 decision1<-rep(0,reptim) decision2<-rep(0,reptim) for (j in 1:reptim) { data1<-data[,j] substat<-rep(0,sublength) for (k in Ln0:Un0) { mean1<-mean(data1[1:k]) mean2<-mean(data1[(k+1):n]) inter1<-cumsum(data1[1:k])-(1:k)*mean1 inter2<-cumsum(data1[(k+1):n])-(1:(n-k))*mean2 substat[k-Ln0+1]<-sqrt(n)*abs(mean1-mean2)/(sqrt(sum(inter1^2)+sum(inter2^2))/n) } stat[j]<-max(substat) #90% decision1[j]<-as.numeric(stat[j]>critval[1,j]) #95% decision2[j]<-as.numeric(stat[j]>critval[2,j]) print(j) } #return(c(sum(decision1)/reptim,sum(decision2)/reptim,dn[1:100])) #} round(sum(decision1)/reptim*100,2) round(sum(decision2)/reptim*100,2) ##################### #Power ##################### n<-100 reptim<-5000 ############################# #series with a change point ############################## unit1<-c(rep(0,.25*n),rep(1,.75*n)) unit2<-c(rep(0,.5*n),rep(1,.5*n)) unit3<-c(rep(0,.8*n),rep(1,.2*n)) extraMat1<-matrix(rep(unit1,reptim),n,reptim) extraMat2<-matrix(rep(unit2,reptim),n,reptim) extraMat3<-matrix(rep(unit3,reptim),n,reptim) ################### #long memory ################### set.seed(100) d<-0 data<-matrix(0,n,reptim) for (j in 1:reptim) { data[,j]<-simARMA0(n,d+1/2) } dataALT<-data+extraMat2*1 ################# #AR(1) set.seed(100) data<-MAR(n,reptim,-0.5) dataALT<-data+extraMat3*0.5 ################# #################### #MA(1) set.seed(100) data<-MMA(n,reptim,0.5) dataALT<-data+extraMat3*2 #########################3 ############################## ################# #main program ################# dn<-NULL m<-as.integer(n^(2/3)) for (j in 1:reptim) { ts<-dataALT[,j] opt3 <- nlminb(start=c(0),obj=locwhittle,lower=-0.499,upper=0.499) dn<-c(dn,opt3\$par) } critval<-matrix(0,2,reptim) seq1<-c(-0.499, -(12:1)*0.04,0,(1:12)*0.04,.499) for (j in 1:reptim) { index1<-which.max(as.numeric(dn[j]1) { #90% critval[1,j]<-((-dn[j]+seq1[index1])*CV[index1-1,1]+(dn[j]-seq1[index1-1])*CV[index1,1])/(seq1[index1]-seq1[index1-1]) #95% critval[2,j]<-((-dn[j]+seq1[index1])*CV[index1-1,2]+(dn[j]-seq1[index1-1])*CV[index1,2])/(seq1[index1]-seq1[index1-1]) } } ############### stat<-rep(0,reptim) L0<-0.15 U0<-0.85 Ln0<-as.integer(n*L0) Un0<-as.integer(n*U0) sublength<-Un0-Ln0+1 decision1<-rep(0,reptim) decision2<-rep(0,reptim) for (j in 1:reptim) { data1<-dataALT[,j] substat<-rep(0,sublength) for (k in Ln0:Un0) { mean1<-mean(data1[1:k]) mean2<-mean(data1[(k+1):n]) inter1<-cumsum(data1[1:k])-(1:k)*mean1 inter2<-cumsum(data1[(k+1):n])-(1:(n-k))*mean2 substat[k-Ln0+1]<-sqrt(n)*abs(mean1-mean2)/(sqrt(sum(inter1^2)+sum(inter2^2))/n) } stat[j]<-max(substat) #90% decision1[j]<-as.numeric(stat[j]>critval[1,j]) #95% decision2[j]<-as.numeric(stat[j]>critval[2,j]) #print(j) } #return(c(sum(decision1)/reptim,sum(decision2)/reptim,dn[1:100])) #} #sum(decision1)/reptim #sum(decision2)/reptim round(sum(decision1)/reptim*100,2) round(sum(decision2)/reptim*100,2)