一、导入数据
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)
结论:差分后的序列为平稳非白噪声序列。
5.4模型识别
acf(D,lag.max=60)
pacf(D,lag.max=60)