I am new to statistics, so pardon any mistakes in my question.

I have two time series $X_i$ and $Y_i$. Assuming that they're stationary AR(1) processes with possibly different means, how do I test for difference of means?

I found this link (but tests only one sample) which discusses using gls: https://stat.ethz.ch/pipermail/r-help/2006-May/105495.html

So I am wondering how to extend this to multiple samples (use a dummy variable?) or if there are any other methods (preferably in R).


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)



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]))


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])


    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]))

    names(res) <-NULL

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)


This function generates N pairs of n sized sample of two AR(1) processes.

Then I ran the following simulation.

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% {

nn <- c(10,30,100,500,1000)

rrmn <- foreach(n=nn,.combine=rbind) %do% {

rrmns <- foreach(s=c(0.1,0.5,1),.combine=rbind) %do% {

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])


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. enter image description here

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$: enter image description here

Now the same graphs for t.test ($\rho=0.1$): enter image description here

And $\rho=0.9$: enter image description here

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.

  • $\begingroup$ What would your approach be when $X_i$ and $Y_i$ are sampled at different times? $\endgroup$ – user27886 Aug 31 '16 at 15:01

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy