季节效应序列分析

原创

已于 2022-06-07 09:54:55 修改 · 3.1k 阅读

36 ·

CC 4.0 BY-SA版权

文章标签：

#r语言

于 2022-04-20 23:28:26 首次发布

本文通过分析‘季节效应.csv’数据，展示了时序图并确定了序列的长期趋势、季节性和随机波动。选用乘法模型进行分解，提取了趋势和季节效应。进一步运用Holt-Winters三参数指数平滑法进行预测，比较了加法和乘法模型的效果，发现乘法模型更优。接着利用ARIMA模型进行差分和模型识别，尝试了多种ARIMA配置并进行了模型检验，最终得出ARIMA(0,1,6)X(0,1,1)12模型可能适合该序列。

一、导入数据

analy <- read.csv("D:\\季节效应.csv",header = T,row.names = 1,encoding = "UTF-8")
analy

二、时序图

analy1 <- ts(analy,start = c(1963,1),end = c(1976,12),frequency = 12)
plot(analy1)

在这里插入图片描述
根据时序图，该序列存在长期向上的趋势、季节性变化和随机波动的影响。随着趋势的递增，序列呈现出和季节变化相关的稳定周期性波动，应采用乘法模型。

三、趋势效应与季节效应的提取

Result1=decompose(analy1,type = "multiplicative")  #乘法模型
plot(Result1)

在这里插入图片描述

# 提取趋势效应
plot(Result1$trend)
abline(lm(Result1$trend~time(analy1)),col="2")   #选择2号颜色为趋势线标线

在这里插入图片描述
结论：总体递增

# 提取季节效应
sea=as.vector(Result1$seasonal)
sea=sea[1:12]
plot(sea,type = 'l')

在这里插入图片描述
结论：如图所示，2月份为最低点，6、7、8月份达到高峰值，存在明显季节性。
因此指数平滑模型选择Holt-Winters三参数指数平滑。

四、Holt-Winters三参数指数平滑

(加法模型)

H=HoltWinters(analy1,seasonal = "additive")  #seasonal = "additive"可省略，默认加法模型
H

在这里插入图片描述
结论：Xt+k=875.583001+1.971655k+Sj，k>=1

#手动预测：
xx1977.1=875.583001+1.971655*1-53.694536
xx1977.2=875.583001+1.971655*2-116.436891
xx1977.3=875.583001+1.971655*3-105.404840
xx1977.1
xx1977.2
xx1977.3

#系统预测
fore_H=forecast(H1, h=36)
fore_H

在这里插入图片描述

(乘法模型)

H1=HoltWinters(analy1,seasonal =  "multiplicative")  #乘法模型
H1

在这里插入图片描述
结论：xt+k=(875.5123349+1.9568538k)Sj

#手动预测
x1977.1=(875.5123349+1.9568538*1)*0.9301067
x1977.2=(875.5123349+1.9568538*2)*0.8613974
x1977.3=(875.5123349+1.9568538*3)*0.8741871
x1977.1
x1977.2
x1977.3

#系统预测
fore_H1=forecast(H1, h=36)
fore_H1

在这里插入图片描述
结论：乘法模型的预测效果优于加法模型。

plot(fore_H1)
lines(fore_H1$fitted,col=3)

在这里插入图片描述

五、ARIMA模型

5.1差分平稳化:

因原序列存在季节性和趋势性，对原序列做差分平稳化，使其转为平稳性模型。

Difference1=diff(analy1)   #一阶差分
plot(Difference1)

在这里插入图片描述
对原序列做1阶差分后，消除了递增的趋势。

D=diff(diff(analy1),12)    #一阶十二步差分
plot(D)

在这里插入图片描述
然后对1阶差分后的序列进行12步差分，消除了季节趋势。

5.2对差分后序列作平稳性检验

adf.test(D,nlag=3)

在这里插入图片描述
结论：p=0.01<a，拒绝原假设，为平稳序列

5.3对差分后序列作白噪声检验

lag = c()  
Q.value = c()  
p.value = c()  
for(i in 1:12){  
  test = Box.test(D,lag=i,type="Ljung-Box")  
  lag[i] = i  
  Q.value[i] = test$statistic  
  p.value[i] = test$p.value  
}  
data.frame(lag=lag,Q.value=Q.value,p.value=p.value)