Since $X_i$ and $Y_i$ are AR(1) processes, we can treat $(X_i,Y_i)$ as a VAR(1) process. Then testing the difference of means means testing the restrictions on VAR coefficients. For this simple case there is a formula derived in Hamilton's book "Time series analysis". The formula in the link is the one before the start of new section. Unfortunately it is not implemented in R, as far as I know. I checked only vars package, but I can say with some confidence that if it is not implemented there, it is not implemented in other packages.
Here is some code I've cobbled to calculate this statistic for simulated AR(1) processes:
n <- 1000
e1 <- rnorm(n,2,sd=1)
e2 <- rnorm(n,2.5,sd=1)
xm <- arima.sim(n,model=list(ar=0.5),innov=e1)
ym <- arima.sim(n,model=list(ar=0.6),innov=e2)
eq1<-dynlm(xm~L(xm)+L(ym))
eq2<-dynlm(ym~L(xm)+L(ym))
res<-rbind(residuals(eq1),residuals(eq2))
omega <- 1/n*tcrossprod(res,res)
xx <- rbind(1,xm[-1000],ym[-1000])
xxinv <- solve(tcrossprod(xx,xx))
c1 <- coef(eq1)[1]
c2 <- coef(eq2)[1]
chi <- (c1-c2)^2/(xxinv[1,1]*(omega[1,1]-2*omega[1,2]+omega[2,2]))
1-pchisq(chi,1)
After playing with different values of means, it seems that statistic is working more or less. Of course the results should be tested for different values of variance of errors and different AR specifications and different sample sizes. On the other hand from theoretical point of view statistic is quite nice, since it allows unequal variances.
Update Out of curiosity I did some MC simulations. For this I wrote a bit more efficient version of test:
ham.test.lm2 <- function(xm,ym) {
n <- length(xm)
df <- cbind(x=xm[-1],y=ym[-1],lx=xm[-n],ly=ym[-n])
eq1<-lsfit(df[,3:4],df[,1])
eq2<-lsfit(df[,3:4],df[,2])
res<-rbind(residuals(eq1),residuals(eq2))
omega <- 1/n*tcrossprod(res,res)
xx <- rbind(1,xm[-n],ym[-n])
ch <- chol(tcrossprod(xx,xx))
xxinv <- chol2inv(ch)
c1 <- coef(eq1)[1]
c2 <- coef(eq2)[1]
chi <- (c1-c2)^2/(xxinv[1,1]*(omega[1,1]-2*omega[1,2]+omega[2,2]))
res<-1-pchisq(chi,1)
names(res) <-NULL
res
}
This runs about 60 times faster than the original. Then I wrote special function for generating two AR(1) processes given their parameters.
gen.arima <-function(n,N,mn1,mn2,rho1,rho2,sd1=1,sd2=2,burn.in=50){
foreach(i=1:N,.combine=c) %dopar% {
e1 <- rnorm(n+burn.in,mn1,sd=sd1)
e2 <- rnorm(n+burn.in,mn2,sd=sd2)
xm <- arima.sim(n+burn.in,model=list(ar=rho1),innov=e1)
ym <- arima.sim(n+burn.in,model=list(ar=rho2),innov=e2)
list(window(cbind(xm,ym),start=burn.in+1))
}
}
This function generates N pairs of n sized sample of two AR(1) processes.
Then I ran the following simulation.
library("iterators")
library("foreach")
library("multicore")
library("doMC")
registerDoMC(16) ##the simulations ran on server with 2 quad-core Xeon processors with HT
rho <- seq(0.1,0.9,by=0.1)
nr <- length(rho)
rr <- numeric()
for(i in 1:nr) {
for(j in (i):nr)
rr <- rbind(rr,cbind(rho[i],rho[j]))
}
mu <- seq(0,2,by=0.1)
rrm <- foreach(m=mu,.combine=rbind) %do% {
cbind(rr,m)
}
nn <- c(10,30,100,500,1000)
rrmn <- foreach(n=nn,.combine=rbind) %do% {
cbind(rrm,n)
}
rrmns <- foreach(s=c(0.1,0.5,1),.combine=rbind) %do% {
cbind(rrmn,s)
}
colnames(rrmns) <- c("rho1","rho2","mudiff","slen","sigma")
cc <- system.time(res <- foreach(i=1:nrow(rrmns),.combine=rbind) %dopar% {
a <- rrmns[i,]
aa <- gen.arima(a[4],500,2,2+a[3],a[1],a[2],a[5],a[5])
hm<-sapply(aa,function(x)ham.test.lm2(x[,1],x[,2]))
tt<-sapply(aa,function(x)t.test(x[,1],x[,2])$p.value)
h1<-sum(hm<0.05)/length(hm)
t1<-sum(tt<0.05)/length(tt)
c(a,h1,t1)
})
write.csv(file="hammc.csv",res)
This simulation calculates the number of times the statistic rejects the null hypothesis that the means are zero for various sample sizes and various combinations of parameters of AR(1) processes. I used the maximum burn in for $\rho=0.9$ when generating AR(1) processes. I also compare the performance of the statistic to simple t-test. The resulting data.frame can be downloaded here.
Here is a few graphs illustrating the perfomance.
This is for AR(1) processes with $\rho=0.1$, 0.1 and 1 indicate the variance of the errors. mudiff is the true difference of the means.
Here is the similar graph for $\rho=0.9$:
Now the same graphs for t.test ($\rho=0.1$):
And $\rho=0.9$:
Here is the code which produced the graphs:
mc <- read.csv("hammc.csv")
colnames(mc)[7] <- "Hamilton"
colnames(mc)[8] <- "t.test"
colnames(mc)[5] <- "sample.size"
qplot(x=sample.size,y=Hamilton,data=mc[mc$sigma!=0.5 & mc$rho1==0.1 & mc$rho2==0.1,],group=mudiff,geom="line",color=mudiff)+facet_wrap(~sigma)
qplot(x=sample.size,y=Hamilton,data=mc[mc$sigma!=0.5 & mc$rho1==0.9 & mc$rho2==0.9,],group=mudiff,geom="line",color=mudiff)+facet_wrap(~sigma)
qplot(x=sample.size,y=t.test,data=mc[mc$sigma!=0.5 & mc$rho1==0.1 & mc$rho2==0.1,],group=mudiff,geom="line",color=mudiff)+facet_wrap(~sigma)
qplot(x=sample.size,y=t.test,data=mc[mc$sigma!=0.5 & mc$rho1==0.9 & mc$rho2==0.9,],group=mudiff,geom="line",color=mudiff)+facet_wrap(~sigma)
Conclusions Note that the graphs measure the power of the statistic for mudiff greater than zero and size of the test for zero mudiff, so ideally the the blue line for mudiff=0 should be straight line at 0.05, and all the other lines should be straight lines at 0.95. Since the statistic is assymptotic this does not happen. What it is clear that t.test is not suitable for higher values of $\rho$, it performs quite nice for $\rho=0.1$, but very bad for $\rho=0.9$. Hamilton's test on the other hand has low power even for larger sample sizes when the true difference in means is not very large. Also the results are the same for error variance 0.1 and 1, so I do not rule out the error in my code.
