1. 程式人生 > >【時間序列】ARIMA模型在鞋服行業銷售預測中的運用

【時間序列】ARIMA模型在鞋服行業銷售預測中的運用

大綱:

-資料處理
-模型構建
-擬合效果


1.資料處理

真實業務資料。來源於特步四川分公司。

資料按照地區可以劃分為:成都/樂山/南充/綿陽等;按品類可以劃分為羽絨服/板鞋/短袖POLO等等。

資料時間跨度:2014年1月~2017年10月

樣本:成都地區跑鞋銷量預測

#讀取資料
library(readxl)
Shoes2_Chengdu <- read_excel("Desktop/跑鞋區域資料.xlsx", sheet = "成都", col_types = c("date", "numeric"))
plot(Shoes2_Chengdu,type="line")
attach(Shoes2_Chengdu)
Quantity=as.numeric(Quantity)
detach(Shoes2_Chengdu)

繪製趨勢圖,對整體趨勢有一個直觀的感受:


節假日因素調整:因為特步在每年的春節,五一,十一等節假日都會舉行促銷活動。除了春節外,其他的節假日都在固定時間。我們所需要調整的是每年春節的不固定日期。

#調整春節因素,2014年春節在一月
a=Shoes2_Chengdu$Quantity[1]
Shoes2_Chengdu$Quantity[1]=Shoes2_Chengdu$Quantity[2]
Shoes2_Chengdu$Quantity[2]=a

如果有過小的值,代表是反季商品促銷。(比如羽絨服在夏季的銷量)直接做0值處理,不做預測(因為對銷售決策制定沒有指導意義)。

#剔除過小的值
for (i in 1:47)
{if(Shoes2_Chengdu$Quantity[i]<20)
  Shoes2_Chengdu$Quantity[i]<-0

}

載入時間序列包,將資料整理成時間序列格式。

#載入時間序列程式包
library(tseries)
library(forecast)
library(dplyr)
library(stats)
#轉換成時間序列
Shoes2_Chengduts<-ts(Shoes2_Chengdu$Quantity[1:47],fre=12,start=c(2014,01))

通過箱形圖檢測異常值:箱形圖可以用來觀察資料整體的分佈情況,利用中位數,25/%分位數75/%分位數,上邊界,下邊界等統計量來來描述資料的整體分佈情況。通過計算這些統計量,生成一個箱體圖,箱體包含了大部分的正常資料,而在箱體上邊界和下邊界之外的,就是異常資料。

(因為樣本量跨越年份較少,所以基本檢測不出異常值)。如果有異常值,通過插值法調整。

#箱型圖檢測異常值
for(i in 1:12){
  bp[,i]<-(Shoes2_Chengduts[c(i,i+12,i+24)])}
bptest=boxplot(bp)
bptest$out



2.模型構建

選擇ARIMA模型。

首先先檢測模型的平穩性。

#檢測一下平穩性
adf.test(Shoes2_Chengduts)
	Augmented Dickey-Fuller Test


data:  Shoes2_Chengduts
Dickey-Fuller = -7.0728, Lag order = 3, p-value = 0.01
alternative hypothesis: stationary

Warning message:
In adf.test(Shoes2_Chengduts) : p-value smaller than printed p-value

可以看出P值很小,我們可以認為模型整體沒有向上或向下的趨勢,基本平穩。

接下來通過ACF圖和PACF圖確定模型階數。


tsdisplay(Shoes2_Chengduts)

可以看出ACF圖並沒有明顯的截尾的表現,明視訊記憶體在季節性。實際上我們用常識也可以判斷出來,銷售每年都是有周期性的。所以我們加入季節性因子,考慮12階次差分。

adf.test(diff(Shoes2_Chengduts,12)

結論同上,平穩。

	Augmented Dickey-Fuller Test

data:  diff(Shoes2_Chengduts, 12)
Dickey-Fuller = -4.005, Lag order = 3, p-value = 0.02093
alternative hypothesis: stationary

做ACF和PACF圖

tsdisplay(diff(Shoes2_Chengduts,12))

ACF在1階處衰減,PACF在1階處截尾。初步確定Arima階數為(0,0,0)(1,0,1)[12]

arima1<-arima(Shoes2_Chengduts,order=c(0,0,0),seasonal=list(order=c(2,0,1),period=12))

我們再來看一下R中自帶的auto命令給出的最優引數。

> auto1<-auto.arima(Shoes2_Chengduts,trace=T)


 ARIMA(2,1,2)(1,0,1)[12] with drift         : Inf
 ARIMA(0,1,0)            with drift         : 830.1212
 ARIMA(1,1,0)(1,0,0)[12] with drift         : 826.1885
 ARIMA(0,1,1)(0,0,1)[12] with drift         : Inf
 ARIMA(0,1,0)                               : 827.9735
 ARIMA(1,1,0)            with drift         : 825.481
 ARIMA(1,1,0)(0,0,1)[12] with drift         : 826.216
 ARIMA(1,1,0)(1,0,1)[12] with drift         : Inf
 ARIMA(2,1,0)            with drift         : 824.9653
 ARIMA(2,1,1)            with drift         : Inf
 ARIMA(3,1,1)            with drift         : Inf
 ARIMA(2,1,0)                               : 822.7944
 ARIMA(2,1,0)(1,0,0)[12]                    : 819.5303
 ARIMA(2,1,0)(1,0,1)[12]                    : 820.9409
 ARIMA(1,1,0)(1,0,0)[12]                    : 824.0301
 ARIMA(3,1,0)(1,0,0)[12]                    : 821.3261
 ARIMA(2,1,1)(1,0,0)[12]                    : 820.4318
 ARIMA(3,1,1)(1,0,0)[12]                    : 820.7753
 ARIMA(2,1,0)(1,0,0)[12] with drift         : 821.4325


 Best model: ARIMA(2,1,0)(1,0,0)[12] 
究竟兩個模型哪個更好?

比較模型的擬合效果,通常使用赤池資訊法則(即AIC)來衡量模型的優劣。為了避免過度擬合帶來的偏差,AIC中增加了對多餘變數的懲罰項。我們可以計算AIC的值,越小的AIC的值說明模型的擬合效果最好。

現在我們需要比較兩個模型的AIC值。

arima1$aic
[1] 811.2276
> auto1$aic
[1] 819.5303

R函式選擇的模型比我們認為比較出來的模型AIC值更高。當然這並不意味著哪個模型“最優”,auto1 的AIC值較高的原因可能是因為加入更多引數,但由於引數值比較大,所以對“過擬合風險”的懲罰項較大。有時候我們需要在衡量過擬合和準確性之間作出抉擇。可以最終兩者都嘗試一下,作出判斷。

3.擬合效果

(1)ARIMA1

fit1<-forecast(arima1,h=1,level = c(99.5))
plot(fit1)
lines(fit1$fitted,col="green")
lines(Shoes2_Chengduts,col='black')

可以看出在16年之前的擬合曲線有滯後現象。明顯感覺有某些因素沒考慮到哦。

再來比較2017年的預測精度。

fitvsture<-data.frame(floor(fit1$fitted[37:47]),Shoes2_Chengdu$Quantity[37:47])
#調整春節因素,2017年春節在1月
a<-fitvsture[1,1]
fitvsture[1,1]=fitvsture[2,1]
fitvsture[2,1]=a
#精確度校驗
d<-Shoes2_Chengdu$Quantity[1:10]
for(i in 1:10){
  d[i]<-1-(abs(fitvsture[i,1]-fitvsture[i,2])/max(fitvsture[i,1],fitvsture[i,2]))
}
print(fitvsture)
print(d)
error=sum(d[1:10])/10
print(error)
> print(fitvsture)
   floor.fit1.fitted.37.47.. Shoes2_Chengduts.37.47.
1                       3889                    6860
2                       3393                    3609
3                       3935                    4162
4                       5646                    4529
5                       4699                    4080
6                       4724                    3453
7                       3749                    2958
8                       3216                    3368
9                       3281                    3078
10                      3115                    4198
11                      3933                     883
> print(d[1:10])
 [1] 0.5669096 0.9401496 0.9454589 0.8021608 0.8682698 0.7309483 0.7890104
 [8] 0.9548694 0.9381286 0.7420200
> error=sum(d[1:10])/10
> print(error)
[1] 0.8277926
#11月份資料不足,不做預測。不計入精度排名。預測總精度基本上可以達80%。

(2)AUTO1

fit1<-forecast(auto1,h=1,level = c(99.5))
plot(fit1)
lines(fit1$fitted,col="green")
lines(Shoes2_Chengduts,col='black')
fitvsture<-data.frame(floor(fit1$fitted[37:47]),Shoes2_Chengduts[37:47])

#調整春節因素,2017年春節在1月
a<-fitvsture[1,1]
fitvsture[1,1]=fitvsture[2,1]
fitvsture[2,1]=a
fitvsture
#精確度校驗
d<-Shoes2_Chengdu$Quantity[1:10]
for(i in 1:10){
  d[i]<-1-(abs(fitvsture[i,1]-fitvsture[i,2])/max(fitvsture[i,1],fitvsture[i,2]))
}
print(fitvsture)
print(d[3:10])
error=sum(d[3:10])/7
print(error)

> print(fitvsture)
   floor.fit1.fitted.37.47.. Shoes2_Chengduts.37.47.
1                       3393                    6860
2                       3889                    3609
3                       3935                    4162
4                       5646                    4529
5                       4699                    4080
6                       4724                    3453
7                       3749                    2958
8                       3216                    3368
9                       3281                    3078
10                      3115                    4198
11                      3933                     883
> print(d[1:10])
 [1] 0.4946064 0.9280021 0.9454589 0.8021608 0.8682698 0.7309483 0.7890104
 [8] 0.9548694 0.9381286 0.7420200
> error=sum(d[1:10])/10
> print(error)
[1] 0.8193475


結論:我們可看出,R軟體自動生成的模型AUTO1 的整體擬合效果較ARIMA1較差一些,但兩者相差很小,精確到小數點後兩位後幾乎可以忽略不計。為了降低過擬合風險和提高預測精度,我們最終還是選擇ARIMA1。