#code for simulating critical values ################ #no trimming ################ simulate3mod<-function(n,p,reptim) { #reptim<-3000 stat<-rep(0,reptim) sublength<-n-1 scale<-(1:(n-1))/n^2*((n-1):1) for (j in 1:reptim) { data<-matrix(rnorm(n*p,0,1),p,n) substat<-rep(0,sublength) for (k in 1:(n-1)) { if (p==1) { mean1<-mean(data[1,1:k]) mean2<-mean(data[1,(k+1):n]) } if (p>1) { mean1<-apply(matrix(data[,1:k],ncol=k),1,mean) mean2<-apply(matrix(data[,(k+1):n],ncol=n-k),1,mean) } inter1<-NULL inter2<-NULL for (h in 1:p) { inter1<-rbind(inter1,cumsum(data[h,1:k])-(1:k)*mean1[h]) inter2<-rbind(inter2,cumsum(data[h,n:(k+1)])-(1:(n-k))*mean2[h]) } M1<-inter1%*%t(inter1)/n^2 M2<-inter2%*%t(inter2)/n^2 #Sigmahat<-matrix(1,p,p) substat[k]<-n*matrix(scale[k]*(mean1-mean2),1,p)%*%solve(M1+M2)%*%matrix(scale[k]*(mean1-mean2),p,1) } stat[j]<-max(abs(substat)) print(j) } return(quantile(stat,c(0.9,0.95,0.99,0.995,0.999))) } #Generate AR(1) model 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]) } ########################## #change in mean, comparison #between fixed bandwidth, data-dependent #bandwidth, and self-normalizated test statistic ################################################## n<-100 reptim<-3000 rho<-0.4 set.seed(100) ts<-MAR(n,reptim,rho) statFB<-rep(0,reptim) statDDB<-rep(0,reptim) statSN<-rep(0,reptim) inter1<-NULL inter2<-NULL scale<-(1:(n-1))/n^2*((n-1):1) ln<-ceiling(n^(1/3)) for (j in 1:reptim) { data<-ts[,j] substat<-rep(0,n-1) for (k in 1:(n-1)) { mean1<-mean(data[1:k]) mean2<-mean(data[(k+1):n]) inter1<-cumsum(data[1:k])-(1:k)*mean1 inter2<-cumsum(data[n:(k+1)])-(1:(n-k))*mean2 M1<-sum(inter1^2)/n^2 M2<-sum(inter2^2)/n^2 substat[k]<-n*(scale[k]*(mean1-mean2))^2/(M1+M2) } statSN[j]<-max(abs(substat)) autocov<-acf(data,type="covariance",plot=F,lag.max=n-1)[[1]][,1,1] #Bartlett kernel, fixed bandwidth sigmasqFB<-sum(c(autocov[ln:1],autocov[2:ln])*(1-abs(seq(-ln+1,ln-1)/ln))) statFB[j]<-max(abs(cumsum(data)-(1:n)*mean(data)))/sqrt(n)/sqrt(sigmasqFB) #data dependent bandwidth demeandata<-data-mean(data) rhohat<-sum(demeandata[2:n]*demeandata[1:(n-1)])/sum((demeandata[1:(n-1)])^2) alphahat<-4*rhohat^2/(1-rhohat^2)^2 lnDDB<-ceiling(1.1447*(alphahat*n)^(1/3)) sigmasqDDB<-sum(c(autocov[lnDDB:1],autocov[2:lnDDB])*(1-abs(seq(-lnDDB+1,lnDDB-1)/lnDDB))) statDDB[j]<-max(abs(cumsum(data)-(1:n)*mean(data)))/sqrt(n)/sqrt(sigmasqDDB) print(j) } ############ #rejection rates ######################## #1.36 corresponds to 5% significance level ############################################# rejFB<-sum(statFB>1.36)/reptim rejDDB<-sum(statDDB>1.36)/reptim rejSN<-sum(statSN>40.13805)/reptim round(c(rejFB,rejDDB,rejSN),3)*100 ################### #function ################### rejmeanshiftSIZE<-function(ts) { n<-dim(ts)[1] reptim<-dim(ts)[2] statFB<-rep(0,reptim) statDDB<-rep(0,reptim) statDDB2<-rep(0,reptim) statSN<-rep(0,reptim) inter1<-NULL inter2<-NULL scale<-(1:(n-1))/n^2*((n-1):1) ln<-ceiling(n^(1/3)) for (j in 1:reptim) { data<-ts[,j] substat<-rep(0,n-1) diff<-rep(0,n-1) for (k in 1:(n-1)) { mean1<-mean(data[1:k]) mean2<-mean(data[(k+1):n]) diff[k]<-mean1-mean2 inter1<-cumsum(data[1:k])-(1:k)*mean1 inter2<-cumsum(data[n:(k+1)])-(1:(n-k))*mean2 M1<-sum(inter1^2)/n^2 M2<-sum(inter2^2)/n^2 substat[k]<-n*(scale[k]*(mean1-mean2))^2/(M1+M2) } khat<-which.max(abs(diff)*sqrt(scale)) residual<-c(data[1:khat]-mean(data[1:khat]),data[(khat+1):n]-mean(data[(khat+1):n])) autocov2<-acf(residual,type="covariance",plot=F,lag.max=n-1)[[1]][,1,1] statSN[j]<-max(abs(substat)) autocov<-acf(data,type="covariance",plot=F,lag.max=n-1)[[1]][,1,1] #Bartlett kernel, fixed bandwidth sigmasqFB<-sum(c(autocov[ln:1],autocov[2:ln])*(1-abs(seq(-ln+1,ln-1)/ln))) statFB[j]<-max(abs(cumsum(data)-(1:n)*mean(data)))/sqrt(n)/sqrt(sigmasqFB) #data dependent bandwidth demeandata<-data-mean(data) rhohat<-sum(demeandata[2:n]*demeandata[1:(n-1)])/sum((demeandata[1:(n-1)])^2) alphahat<-4*rhohat^2/(1-rhohat^2)^2 lnDDB<-ceiling(1.1447*(alphahat*n)^(1/3)) sigmasqDDB<-sum(c(autocov[lnDDB:1],autocov[2:lnDDB])*(1-abs(seq(-lnDDB+1,lnDDB-1)/lnDDB))) statDDB[j]<-max(abs(cumsum(data)-(1:n)*mean(data)))/sqrt(n)/sqrt(sigmasqDDB) rhohat2<-sum(residual[2:n]*residual[1:(n-1)])/sum((residual[1:(n-1)])^2) alphahat2<-4*rhohat2^2/(1-rhohat2^2)^2 lnDDB2<-ceiling(1.1447*(alphahat2*n)^(1/3)) sigmasqDDB2<-sum(c(autocov2[lnDDB2:1],autocov2[2:lnDDB2])*(1-abs(seq(-lnDDB2+1,lnDDB2-1)/lnDDB2))) statDDB2[j]<-max(abs(cumsum(data)-(1:n)*mean(data)))/sqrt(n)/sqrt(sigmasqDDB2) print(j) } rejFB<-sum(statFB>1.36)/reptim rejDDB<-sum(statDDB>1.36)/reptim rejDDB2<-sum(statDDB2>1.36)/reptim rejSN<-sum(statSN>40.13805)/reptim qFB<-quantile(statFB,0.95) qDDB<-quantile(statDDB,0.95) qDDB2<-quantile(statDDB2,0.95) qSN<-quantile(statSN,0.95) result<-rbind(round(c(rejFB,rejDDB,rejDDB2,rejSN),3)*100,round(c(qFB,qDDB,qDDB2,qSN),3)) return(result) } set.seed(100) ts<-MAR(100,5000,0.4) rejmeanshiftSIZE(ts)