How to test difference between times series - does "time series ANOVA" exist? Description: Temperature if six objects (obj_A, obj_B,...obj_F) was measured in 1 hour intervals (1,2...10)
Objects were under the influence of two treatments (A and B). Treatment A = obj_A, obj_B, obj_C; treatment B = obj_D, obj_E, obj_F.
Problem is that measured values of each object are serially correlated, therefore I can not use classical
ANOVA. How to take into account such fact?
# example data
my.data <- data.frame(object = rep(c("obj_A","obj_B","obj_C",
                                     "obj_D","obj_E","obj_F"),
                                      each = 10),
                      time = rep(c(1:10), times = 6),
                      treatment = rep(c("A","B"), each = 30),
                      value = c(4,4,7,8,8,10,8,12,14,12,
                                8,8,12,12,10,12,10,11,12,16,
                                12,12,11,13,12,16,16,14,16,20,
                                11,20,23,27,31,29,31,32,28,30,
                                12,16,17,23,22,24,33,31,31,32,
                                14,13,19,20,24,26,24,28,25,23))

# converting values to time series object
obj_A <- ts(my.data$value[my.data$object=="obj_A"],
            start = 1, end = 10, frequency = 1)
obj_B <- ts(my.data$value[my.data$object=="obj_B"],
            start = 1, end = 10, frequency = 1)
obj_C <- ts(my.data$value[my.data$object=="obj_C"],
            start = 1, end = 10, frequency = 1)
obj_D <- ts(my.data$value[my.data$object=="obj_D"],
            start = 1, end = 10, frequency = 1)
obj_E <- ts(my.data$value[my.data$object=="obj_E"],
            start = 1, end = 10, frequency = 1)
obj_F <- ts(my.data$value[my.data$object=="obj_F"],
            start = 1, end = 10, frequency = 1)

# plot -> blue = treatment A; red = treatment B
ts.plot(obj_A, obj_B, obj_C, obj_D, obj_E, obj_F,
        col=c("deepskyblue","deepskyblue1","deepskyblue2",
              "darkred","indianred","indianred1"),
        lwd = 2.5, lty = 2, xlab = "time", ylab = "temperature")


How to rigorously test whether temperature of objects differ according to used treatment, but without ignoring serial
 correlation?
 A: You notice correcly that usual ANOVA cannot handle this type of very heteroskedastic and highly dependent time point data. But also, other repeated measures or multivariate procedures would fail because you have only 6 replications but 10 time points, so called high dimensional data. In the end, it is proven that there cannot be an exact test for data like that.
However, there is a good approximative test, it's the Huynh-Feldt-procedure, see this technical report for references and for a generalization to possibly unequal covariance matrices that you may use. It is defined even for your very small data set. There are no assumptions on the variances and covariances. It is intended for normally distributed data, but I think you're fine with it. 
You can test if the curves differ (set T = diag(10)) and if there is an interaction between treatment and time point (set T = diag(10) - 1/10). It's too simple for an extra R-package. Check this: 
Tmat <- diag(10) # or Tmat <- diag(10)-1/10
Y1 <- t(matrix(my.data[my.data$treatment=="B","value"],10,3)) %*% Tmat
    Y2 <- t(matrix(my.data[my.data$treatment=="A","value"],10,3)) %*% Tmat

Mean1 <- apply(Y1,2,mean)
Mean2 <- apply(Y2,2,mean)
S1 <- cov(Y1)
S2 <- cov(Y2)

F <- sum((Mean1-Mean2)**2) / (sum(diag(S1/3 + S2/3))) # see eq. 3.18
trS1sq <- 3/2 * ( (sum(diag(S1)))**2 - 2/3 * sum(S1**2)) # see eq. 3.26
trS2sq <- 3/2 * ( (sum(diag(S2)))**2 - 2/3 * sum(S2**2)) 
trS1S1 <- 3/2 * ( sum(S1**2) - (sum(diag(S1)))**2/2)  # 3.27
trS2S2 <- 3/2 * ( sum(S2**2) - (sum(diag(S2)))**2/2)  # 3.27

f <- (trS1sq + trS2sq + 2*sum(diag(S1))*sum(diag(S2))) / 
    (trS1S1+trS2S2+2*sum(S1*S2))  # see p. 16
f0 <- (trS1sq + trS2sq + 2*sum(diag(S1))*sum(diag(S2))) / 
    (trS1S1/2 + trS2S2/2)

p.value <- 1-pf(F,f,f0)

A: I believe your problem is a textbook case of longitudinal data analysis for which there is an extensive statistical methodology developed. In particular, analysis of response profiles is an approach that I would employ in this case. 
To motivate this, I will provide some background from Applied Longitudinal Analysis by G.M. Fitzmaurice:

1.2 Longitudinal and clustered data
The defining feature of longitudinal studies is that measurements of
  the same individuals are taken repeatedly through time, thereby
  allowing the direct study of change over the time. The primary goal of
  a longitudinal study is to characterize the change in response over
  the time and the factors that influence change. (...)
A distinctive feature of longitudinal data is that they are
  clustered. In longitudinal studies the clusters are composed of the repeated measurements obtained from a single individual at different
  occasions. Observations within a cluster will typically exhibit
  positive correlation, and this correlation must be accounted for in
  the analysis. Longitudinal data also have a temporal order; the first
  measurement within a cluster necessarily comes before the second
  measurement, and so on. The ordering of the repeated measures has
  important implications for analysis.

and regarding analyzing response profiles choice: 

Methods for analyzing response profiles are appealing when there is a
  single categorical covariate (perhaps denoting different treatment or
  exposure groups) and when no specific a priori pattern for the
  differences in the response profiles between groups can be specified.
  When repeated measures are obtained at the same sequence of occasions,
  the data can be summarized by the estimated mean response at each
  occasion, stratified by levels of a group factor. At any given level
  of the group factor, the sequence of means over time is referred to as
  the mean response profile.

The nlme R library contains a glm function to evaluate such model, giving you a longitudinal "equivalent" of Anova:
library(nlme)

# Model assuming the same variance for each time point
gls.fit <- 
  gls(value ~ factor(time) + treatment, 
      data = my.data,
      corr = corSymm(form = ~ 1 | object),
      control = glsControl(tolerance = 0.01, msTol = 0.01, 
                           maxIter = 1000000, msMaxIter = 1000000))
summary(gls.fit)

# Model allowing for different variance structure for each time point 
gls.fit.diff.var <- 
  gls(value ~ factor(time) + treatment, 
      data = my.data,
      corr = corSymm(form = ~ 1 | object),
      weights = varIdent(form = ~ 1 | factor(time)),
      control = glsControl(tolerance = 0.01, msTol = 0.01, 
                           maxIter = 1000000, msMaxIter = 1000000))
summary(gls.fit.diff.var)

Unfortunately, in this case model estimation does not coverage (even though I changed control parameters for more convenient). I am afraid there is just too few cases in your data set to estimate parameters. 
