检测季节性R中
问题描述:
问题:检测使用周期图和FFT在日常数据的周期性图案R.检测季节性R中
的问题是如何在R代码里面的周期图检测月度,季度,半年度,annual..etc数据中的循环模式。换句话说,我需要检测为低频率(即:11个月=> 2 * PI/365,6个月=> 4 * PI/365等)的周期性图案的存在
重复的例子:
library(weatherData)
w2009=getWeatherForYear("sfo",2009)
w2010=getWeatherForYear("sfo",2010)
w2011=getWeatherForYear("sfo",2011)
w2012=getWeatherForYear("sfo",2012)
w2013=getWeatherForYear("sfo",2013)
w2014=getWeatherForYear("sfo",2014)
w=rbind(w2009,w2010); w=rbind(w,w2011); w=rbind(w,w2012)
w=rbind(w,w2013); w=rbind(w,w2014)
# Next we analyze the periodograms
# This is IMAGE 1
TSA::periodogram(w$Max_TemperatureF)
# Next: I dont really know to use this information
GeneCycle::periodogram(w$Max_TemperatureF)
# Next THIS IS IMAGE 2
stats::spectrum(w$Max_TemperatureF)
# I also tried . This is IMAGE 3
f.data <- GeneCycle::periodogram(tmax)
harmonics <- 1:365
plot(f.data$freq[harmonics]*length(tmax),]
f.data$spec[harmonics]/sum(f.data$spec),
xlab="Harmonics (Hz)", ylab="Amplitute Density", type="h")
阅读的答案后,我所做的:
per <- TSA::periodogram(w$Max_TemperatureF,lwd = 1)
x <- which(per$freq < 0.01)
plot(x = per$freq[x], y = per$spec[x], type="s")
我的问题是到底有什么意思呢?我们有季节性循环吗?
答
如果你正在寻找一个长期(365天),你会发现它非常低频
> 1/365
[1] 0.002739726
下,实际上,你可以在这个值上看到你的第一个图像左侧的峰值。过滤器,以较低的频率,如果你想放大:
per <- TSA::periodogram(w$Max_TemperatureF,lwd = 1)
x <- which(per$freq < 0.01)
plot(x = per$freq[x], y = per$spec[x], type="s")
另一种方式来寻找周期为自相关估计(acf
):
acf(w$Max_TemperatureF, lag.max = 365*3)
参见季节性分解:
ts1 <- ts(data = w$Max_TemperatureF, frequency = 365)
plot(stl(ts1, s.window = "periodic"))
谢谢,我在低频上添加了周期图。但我认为要找到一个年度,频率是2 * pi/365,为什么1/365足够了?所以频率是每个$ freq有一个尖峰的权利?那么你如何计算振幅?谢谢 !!! – Oniropolo 2015-03-03 02:59:29
对。如果你有频率“可疑”,你可以去季节性分解:价值分解为周期部分(有幅度),趋势和剩余(见'stl')。 – bergant 2015-03-03 03:06:05
你的数据是每日时间序列 - 这就是为什么年度频率是1/365。 – bergant 2015-03-03 03:11:53