2
$\begingroup$

I have a time series of hourly activity levels for a period of about 2 months (1704 observations). There is obviously a strong "seasonal" component (freq=24) to this time series, with activity showing daily fluctuations between night and day.

I'm interested in looking at the relationship between hourly activity and environmental variables, but I'm thinking I need to remove the seasonality first, because otherwise there is a strong positive relationship between activity and air temperature - but that would mostly be because it's warmer during the day and we're more active during the day - but what I would like to find out is if we are more active on warm days than on cold days, and how much of a lag there is between increasing temp and increasing activity.

I ran some cross-correlation functions to try and address these questions, but I think the strong 24 hour cyclicity is affecting those results. I've decomposed the time series using "decompose" in R, which is neat, but I don't know how to use that information to give an actual, deseasonalized time series to work with.

Here is a sample of the data:

[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
Improve this question
$\endgroup$

2 Answers 2

Reset to default
1
$\begingroup$

You can use the STL package to decompose the season of a timeseries

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

I would remove it by days and hours and check the error in each of them. You might also want to check if trend decomposition is necessary.

Improve this answer
$\endgroup$
2
  • $\begingroup$ Thanks for the reply! I definitely like that stl seems more flexible than decompose. For anyone else reading, this article was really helpful: wessa.net/download/stl.pdf I still have a couple of questions: 1) From the above article, I figured out some starting values to use for my arguments in stl. I just can't figure out how to translate between the terminology used in the paper - n(p), n(i), n(s) - and the arguments in stl. 2) Will stl give me a time series that can be used for cross-correlation? I see the 3 components in time.series but not sure what to do with them! $\endgroup$
    – sitka
    Apr 5, 2015 at 17:25
  • $\begingroup$ To anyone else reading, I don't know if this is correct or rigorous, but it works. I basically followed the code from this blog post: align-alytics.com/blog/… that uses stl and explains how to extract the decomposed series: decomposed <- stl(time.series, s.window="periodic") seasonal <- decomposed$time.series[,1] trend <- decomposed$time.series[,2] remainder <- decomposed$time.series[,3] $\endgroup$
    – sitka
    Apr 6, 2015 at 18:15
2
$\begingroup$

For deseasonalization of Daily and Monthly data, use deseasonalize package.

For deseasonalization of Quarterly data, use causfinder::deseasonalizeQ function:

#######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)
Improve this answer
$\endgroup$
3
  • $\begingroup$ Hi, thank you for your reply! I've looked into both those options, but do you know if they'll work with hourly data on a 24 hour cycle? It looks like for deseasonalize you have to choose either daily or monthly as an argument - I tried setting s=12 in getds for funsies but it still looks very seasonal. Can you adjust for frequency in deseasonalizeQ? $\endgroup$
    – sitka
    Apr 5, 2015 at 17:41
  • $\begingroup$ @sitka, causfinder::deseasonalizeQ makes use of "centered moving averages" during deseasonalization. So, you can change the filter "c(1/8, 1/4, 1/4, 1/4, 1/8)" acc. to your needs. I suggest you to repeat Jason's work in Excel completely. Then, copy .xlsx this time to work for your hourly/daily business. Now, repeat Jason's logic mutatis mutandis for your case. Exctract what the filter coefficients can be. Without doing that Excel job, I am not good enough to say you the relevant filter coefficients in your case of hourly/daily data. There is also another option I will say in next post. $\endgroup$
    – Erdogan CEVHER
    Apr 5, 2015 at 18:48
  • $\begingroup$ @sitka, You may try to resort to non-R options for deseasonalization since that phase in fact merely is of PRE-PROCESSING data. Hence, there are better ways of deseasonalization outside of R: "TRAMO-SEATS" > "X-13" > "centered moving averages". Work this: eviews.com/EViews8/ev8ecx13_n.html Note that TRAMOSEATS is better than X-11 or X-13 since it is model based. Hence, perform complete deseasonalization outside of R. Then do whatever you need to do in R. Eviews' .wf1 to .csv to R's dataframe! To my knowledge, R has no TRAMO-SEATS function to perform deseasonalization. $\endgroup$
    – Erdogan CEVHER
    Apr 5, 2015 at 18:54

Your Answer

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

Not the answer you're looking for? Browse other questions tagged or ask your own question.