一、导入数据
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)
结论:acf拖尾,pacf拖尾,12阶以内的ACF和PACF均不截尾,
故用ARMA(1,1)模型提取差分后序列的短期自相关信息
自相关系数1阶截尾,故用 ARMA(0,1)12 拟合季节自相关信息
偏自相关系数1阶截尾,故用 ARMA(1,0)12 拟合季节自相关信息
故模型识别为:ARIMA(1,1,1)X(0,1,1)12 和 ARIMA(1,1,1)X(1,1,0)12 或者ARIMA(0,1,6)X(0,1,1)12
5.5模型拟合与检验
(1)拟合ARIMA(1,1,1)X(0,1,1)12
fit_1<-arima(analy1,order=c(1,1,1),seasonal = list(order=c(0,1,1),period=12))
fit_1
检验模型ARIMA(1,1,1)X(0,1,1)12
ts.diag(fit_1)
(2)拟合模型ARIMA(1,1,1)X(1,1,0)12
fit_2=arima(analy1,order = c(1,1,1),seasonal=list(order=c(1,1,0),period=12))
fit_2
检验模型ARIMA(1,1,1)X(1,1,0)12
ts.diag(fit_2)
结论:两个模型都不通过检验
(3)拟合模型ARIMA(0,1,6)X(0,1,1)12
fit_3=arima(analy1,order = c(0,1,6),seasonal=list(order=c(0,1,1),period=12))
fit_3
∇
∇
12
∇∇_{12}
∇∇12
X
t
X_t
Xt =
(1-0.924B -0.0793
B
2
B^{2}
B2 -0.3075
B
3
B^{3}
B3 + 0.0590
B
4
B^{4}
B4 - 0.0684
B
5
B^{5}
B5 + 0.3202
B
6
B^{6}
B6) * (1 - 0.4184
B
12
B^{12}
B12)
e
t
e_t
et
检验模型ARIMA(0,1,6)X(0,1,1)12
ts.diag(fit_3)
(4)拟合模型ARIMA(0,1,6)X(1,1,0)
fit_4=arima(analy1,order = c(0,1,6),seasonal=list(order=c(1,1,0),period=12))
fit_4
检验模型ARIMA(0,1,6)X(1,1,0)
ts.diag(fit_4)
(5)拟合模型ARIMA(3,1,2)X(1,1,0)
x4_fit<-arima(analy1,order=c(3,1,2),seasonal = list(order=c(1,1,0),period=12))
x4_fit
∇ ∇ 12 ∇∇_{12} ∇∇12 X t X_t Xt = 1 − 1.7905 B + 0.7905 B 2 1 − 0.8229 B − 0.0030 B 2 + 0.3107 B 3 \frac{1-1.7905B+0.7905B^2}{1-0.8229B-0.0030B^2+0.3107B^3} 1−0.8229B−0.0030B2+0.3107B31−1.7905B+0.7905B2*( 1 − 0.4073 B 1-{0.4073B} 1−0.4073B) e t e_t et
六、评判模型的相对优劣
data.frame(AIC(fit_1), AIC(fit_4), BIC(fit_1), BIC(fit_4))