6
$\begingroup$

I need to identify seasonality/ periodicity of a dataset so as to develop an ARMAX model. This is what the original time-series looks like enter image description here

I have plotted the periodogram of the dataset.

Ps: I used the first difference of the original time series so as to remove the trend from the original time series . This what the periodogram looks likeenter image description here

I also removed the linear trend from the data and estimated the periodogram. Below is what it looks like

enter image description here

I also estimated periodogram of the raw data set without any differencing. Clearly, it shows there is some trend in the data. I am using following R-code to generate it

A <- read.xlsx("Interpolation_Data.xlsx")
y <- A[,2]
x <- as.numeric(A[,4])
diff.y <- diff(y)
diff.x <- diff(x)
Fs<- 1/60
n=length(diff.y)
xdft <-fft(diff.y)
xdft <- xdft[1:(n/2+1)]
psdx <- (1/(Fs*n))* abs(xdft)^2
psdx[c(-1,-n)] <- 2*psdx[c(-1,-n)]
f <- seq(0,Fs/2,by = Fs/n)
plot(f,10*log10(psdx),type='l')

This is the link to data set

$\endgroup$
4
  • $\begingroup$ Isn't this already asked and answered at stats.stackexchange.com/q/16117? $\endgroup$
    – whuber
    Jul 22, 2015 at 16:13
  • 1
    $\begingroup$ @whuber Yes.However,this periodogram shows multiple periods. my question is , what is the next step after you estimate the periodogram? All the literatures that I have reviewed do not clearly specify the next step after periodogram $\endgroup$
    – Spandyie
    Jul 22, 2015 at 16:20
  • $\begingroup$ is this daily data ? $\endgroup$
    – IrishStat
    Jul 22, 2015 at 22:52
  • $\begingroup$ @IrishStat the data is collected every 60 seconds $\endgroup$
    – Spandyie
    Jul 23, 2015 at 0:37

1 Answer 1

7
$\begingroup$

Well, the periodogram after taking first differences doesn't indicate any clear periodicity. However, be aware that taking first differences amplifies high-frequency components, which should appear in the power spectrum as a quadratic trend, and in the log-power spectrum as a log-shaped trend (which is roughly compatible with what we see here).

My recommendations:

– First compute the periodogram without any preprocessing.

– Then, remove the linear trend, but not by using first differences but by fitting a line by regression and subtracting the result from the data (there should also be a "detrend" function in R).

– Use a proper spectral estimation function. I'm not familiar with R, but I'm sure it contains an implementation of Welch's modified periodogram method.


Update for the updated question:

The time series exhibits a dominant period of roughly 360 samples, which for a sampling rate of 1 per minute means 360 minutes. The dominant frequency should therefore be about 0.0028 min$^{-1}$. This seems to be consistent with the periodogram after subtracted trend. Try zooming into the low-frequency range to more precisely determine the peak location. Superimposed to that is a secondary oscillation of 180 samples period or 0.0056 min$^{-1}$ frequency.

Analyzing your data myself I get a main peak at 0.0029 min$^{-1}$ and a secondary peak at 0.0059 min$^{-1}$ (factor 2, the first harmonic). The secondary peak is smaller by a factor of about 36, or by 16 dB. I can share my code, but I am using Matlab, not R.

Btw., the frequency scale on your plots does not appear to be consistent; sometimes it goes up to 0.008, sometimes to 0.5? I think the former gives the frequency in Hz and the latter in min$^{-1}$.

$\endgroup$
8
  • 1
    $\begingroup$ @Spandy, thanks, but now I feel that this spectrum might actually be disturbed by a trend (the strong increase towards 0 is an indicator). Would you also follow my second recommendation: subtract a linear trend? $\endgroup$
    – A. Donda
    Jul 22, 2015 at 19:09
  • 1
    $\begingroup$ @Spandy, see my updated answer. $\endgroup$
    – A. Donda
    Jul 23, 2015 at 2:10
  • 1
    $\begingroup$ @Spandy, in this case normalized frequency is identical to frequency in min$^{-1}$, because the sampling rate is 1 (per minute). I think in this case that's the most useful frequency scale. – The difference in estimated peak frequency is probably due to the fact that the algorithm I use 0-padded the time series to 1024 samples before spectral analysis; with such a slow rhythm it is hard to determine the frequency exactly. $\endgroup$
    – A. Donda
    Jul 23, 2015 at 2:58
  • 1
    $\begingroup$ My code: x = detrend(x); [P, fs] = pwelch(x, numel(x), [], [], 1); semilogy(fs, P, '.-'), using pwelch from Matlab's Signal Processing Toolbox. – If my answer helped you, consider accepting it. $\endgroup$
    – A. Donda
    Jul 23, 2015 at 2:59
  • 1
    $\begingroup$ really appreciate your help $\endgroup$
    – Spandyie
    Jul 23, 2015 at 3:04

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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