Benilton, novamente obrigado pela resposta e tempo. Antes de executar o seu código, porém, eu tentei rodar o seguinte:

##
P <- 2000:2010
Rep <- letters[1:5]
TF <- 1:2
CC <- paste("H-",101:299,sep="")
a <- expand.grid(P=P,R=Rep,TF=TF,CC=CC)
a$y <- rnorm(nrow(a))
a$z <- runif(nrow(a))

# Alternativa 3
library(snow)
ncpu<-4
cl<-makeSOCKcluster(rep("localhost",ncpu))
clusterEvalQ(cl,library(zoo))
clusterExport(cl,'a')

ind<-split(seq(nrow(a)),list(a$R,a$CC,a$TF))
a.2<-parLapply(cl,ind,function(i,k){
cbind(ma.y=c(rep(NA,k-1),rollmean(a$y[i],k)),
ma.z=c(rep(NA,k-1),rollmean(a$z[i],k)))
},k=3)
a.2.1<-cbind(a,do.call(rbind,a.2))
a.2.2<-na.omit(a.2.1)


Eu fiz um split dos índices para depois juntar as médias móveis 'a.2' na base 'a'. Tudo vai muito bem até que ao juntar as duas partes, quando para algumas linhas os valores ficam incorretos.

Uma outra tentativa foi:

##
# Alternativa 4
ma.fun<-function(x,k){
  c(rep(NA,k-1),rollmean(x,k))
}

library(data.table)
a.3<-data.table(a)
a.3.1<-a.3[,list(P,y=ma.fun(y,k=3),z=ma.fun(z,k=3)),by=list(R,CC,TF)]
a.3.2<-na.omit(a.3.1)


Mas na base original tem algumas combinações de níveis inexistentes e dá erro, avisando que estes 'tratamentos' tem comprimento zero. Como eu poderia lidar com isso na função 'ma.fun'?

Grato,
Rubem