¿Cómo puedo producir una serie temporal desestacionalizada en R?

Question

Tengo una serie temporal de niveles de actividad por hora para un periodo de unos 2 meses (1704 observaciones). Evidentemente, hay un fuerte componente "estacional" (freq=24) en esta serie temporal, y la actividad muestra fluctuaciones diarias entre la noche y el día.

Me interesa estudiar la relación entre la actividad horaria y las variables ambientales, pero creo que primero tengo que eliminar la estacionalidad, porque de lo contrario hay una fuerte relación positiva entre la actividad y la temperatura del aire -pero eso sería sobre todo porque hace más calor durante el día y estamos más activos durante el día-, pero lo que me gustaría averiguar es si somos más activos en los días cálidos que en los fríos, y cuánto desfase hay entre el aumento de la temperatura y el aumento de la actividad.

He ejecutado algunas funciones de correlación cruzada para intentar responder a estas preguntas, pero creo que la fuerte ciclicidad de 24 horas está afectando a esos resultados. He descompuesto la serie de tiempo usando "descomponer" en R, que es limpio, pero no sé cómo utilizar esa información para dar una serie de tiempo real, desestacionalizado para trabajar.

He aquí una muestra de los datos:

[1] 24 16 40 48 50 38 24  4  4  5  3  6  4  4  4  3 12 63 55 42 56 20 10 26 45 47 66 64 59
[30] 54 24  5  6  2  4  3  6 10  6  2 13 39 26 17 24 13 19 26 17 32 54 68 58 39 20  0  3  2
[59]  8  2  4  1  5 11  5 60 57 54 40 40 53 74 40 42 57 46 46 26  9  8  4  6 14  8  5  3  2
[88]  7 19 47 53 43 53 51 55 64 48 64 57 56 52 34 22  8  5  6  4  6  3  4  7  6 27 40 48 41
[117] 43 51 50 44 56 64 68 46 49 35 16  2 14  3  7  3 13  3  3  2 14 49 62 42 41 57 52 63 32
[146] 54 59 60 68 24 12  2  2  2  2  7  6  5  9 10 26 53 50 59 28 45 47 44 48 55 59 77 86 33
[175] 18 16 10  6  9  9 14  7  9  7  9 46 57 41 33 32 34 29 39 39 27 26  4 10  9  6  6  2  4
[204]  1  2  2  4  4 17 50 47 24 27 34 26 38 20  6 20 15 25  8  2  2  3  6  4  3  3  4  4  2
[233] 18 41 63 52 37 32 32 28 48 20  6 10  9  7  5 10  4  3  4  7  4  3  4 10  8 56 47 50 27
[262] 30 22 38 38 28 33 24 18 12 14  2 10  4 21  4  5  6  4  4 20 41 46 16  8 20 24 21 16 27
[291] 10  6 14  5  6  6 12  2 10  7

Answer 1

Para desestacionalizar los datos diarios y mensuales, utilice el paquete desestacionalizar.

Para desestacionalizar los datos trimestrales, utilice causfinder::desestacionalizarQ función:

#######deseasonalizeQ: deseasonalize quarterly data #######
# Inspired by excellent work of Jason Delaney on Quarterly Deseasonalize: https://www.youtube.com/watch?v=Jr_2nj6M7L8
# Mutatis mutandis replica of Jason's logic in R 

sales <- ts(c(6,15,10,4,10,18,15,7,14,26,23,12,19,28,25,18,22,34,28,21,24,36,30,20,28,40,35,27))

deseasonalizeQ <- function (x){
x <- ts(x)
#Step1: Centered moving averages: create cma time series having the same length with the original time series x
# cma has 2 NAs on both ends.
cma <- filter(x, filter = c(1/8, 1/4, 1/4, 1/4, 1/8), sides=2)

#Step2: Ratios = Original time series / centered moving averages
ratio <- x/cma

#Step3: Unadjusted 4 seasonal indexes
unadj4si <- ts(1:4)
# floor((length(x)-4)/4)  #"-4" is 4 NA at both ends; below "-1" is due to starting "0:" in multiplication

unadj4si[1] <- mean(ratio[3+4*(0:(floor((length(x)-4)/4) - 1))])
unadj4si[2] <- mean(ratio[4+4*(0:(floor((length(x)-4)/4) - 1))])
unadj4si[3] <- mean(ratio[5+4*(0:(floor((length(x)-4)/4) - 1))])
unadj4si[4] <- mean(ratio[6+4*(0:(floor((length(x)-4)/4) - 1))])

#Step4: Adjusted 4 seasonal indexes
adj4si <- ts(1:4)
adj4si[1] <- unadj4si[1]/mean(c(unadj4si[1],unadj4si[2],unadj4si[3],unadj4si[4]))
adj4si[2] <- unadj4si[2]/mean(c(unadj4si[1],unadj4si[2],unadj4si[3],unadj4si[4]))
adj4si[3] <- unadj4si[3]/mean(c(unadj4si[1],unadj4si[2],unadj4si[3],unadj4si[4]))
adj4si[4] <- unadj4si[4]/mean(c(unadj4si[1],unadj4si[2],unadj4si[3],unadj4si[4]))

#Step5: Propogated adjusted seasonal indexes
propadjsi <- ts(1:length(x))

propadjsi[3+4*(0:(floor((length(x)-4)/4) - 1))] <- adj4si[1]
propadjsi[4+4*(0:(floor((length(x)-4)/4) - 1))] <- adj4si[2]
propadjsi[5+4*(0:(floor((length(x)-4)/4) - 1))] <- adj4si[3]
propadjsi[6+4*(0:(floor((length(x)-4)/4) - 1))] <- adj4si[4]

propadjsi[1] <- adj4si[3]
propadjsi[2] <- adj4si[4]
propadjsi[length(x)-1] <- adj4si[1]
propadjsi[length(x)] <- adj4si[2]

#Step6: Deseasonalized values
out <- x/propadjsi  # deseasonalized = x/propadjsi
out
}

deseasonalizeQ(sales)

#Time Series:Start = 1, End = 28, Frequency = 1  
[1]  6.673117 11.015814  8.941810  6.442787 11.121862 13.218976 13.412714 
[8] 11.274878 15.570607 19.094077 20.566162 19.328362 21.131538 20.562852
[15] 22.354524 28.992543 24.468097 24.969177 25.037067 33.824633 26.692469
[22] 26.437953 26.825429 32.213936 31.141214 29.375503 31.296333 43.488814

###### Plots ########
salesSA <- deseasonalizeQ(sales)
salesSAsalesORJ <- cbind(salesSA, sales)

plot(salesSAsalesORJ, plot.type="single", main="Compare", ylab="values", col=c("blue", "red"), lty=1:2)
legend(10, 40, legend=c("salesSA","sales"), col=c("blue", "red"), lty=1:2)

Answer 2

Puede utilizar el paquete STL para descomponer la estación de una serie temporal

https://stat.ethz.ch/R-manual/R-devel/library/stats/html/stl.html

Yo lo quitaría por días y horas y comprobaría el error en cada uno de ellos. También podrías comprobar si es necesaria la descomposición de tendencias.