时间序列分析实战（三）：时序因素分解法

🍉CSDN小墨&晓末:https://blog.csdn.net/jd1813346972

   个人介绍: 研一｜统计学｜干货分享
         擅长Python、Matlab、R等主流编程软件
         累计十余项国家级比赛奖项，参与研究经费10w、40w级横向

1 目的

  有某欧洲小镇1963年1月至1976年12月每月旅馆入住的房间数构成时间序列$x_t$，选择适当模型对序列进行因素分解，同时考察序列的规律。

2 读取数据

  运行程序：

data=scan("F:/时间序列分析/数据.txt")
data1=ts(data,start = c(1963,1),frequency = 12)
data1

   运行结果：

##       Jan  Feb  Mar  Apr  May  Jun  Jul  Aug  Sep  Oct  Nov  Dec
## 1963  501  488  504  578  545  632  728  725  585  542  480  530
## 1964  518  489  528  599  572  659  739  758  602  587  497  558
## 1965  555  523  532  623  598  683  774  780  609  604  531  592
## 1966  578  543  565  648  615  697  785  830  645  643  551  606
## 1967  585  553  576  665  656  720  826  838  652  661  584  644
## 1968  623  553  599  657  680  759  878  881  705  684  577  656
## 1969  645  593  617  686  679  773  906  934  713  710  600  676
## 1970  645  602  601  709  706  817  930  983  745  735  620  698
## 1971  665  626  649  740  729  824  937  994  781  759  643  728
## 1972  691  649  656  735  748  837  995 1040  809  793  692  763
## 1973  723  655  658  761  768  885 1067 1038  812  790  692  782
## 1974  758  709  715  788  794  893 1046 1075  812  822  714  802
## 1975  748  731  748  827  788  937 1076 1125  840  864  717  813
## 1976  811  732  745  844  833  935 1110 1124  868  860  762  877

3 原始时序特征

  运行程序：

plot(data1,sub='图1 1963-1976年旅馆入住房间数时序图',ylab="入住房间数")

  运行结果：

  从图1可以看出，该时间序列具有明显的线性递增趋势以及以年为周期的季节效应，同时，不具有大的经济周期循环和特殊交易日的信息，所以可以确定这个序列受到三个因素的影响：长期趋势、季节效应和随机波动。由于随着趋势的递增，每个季节的振幅也在增大，这说明季节效应受到趋势的影响，因此选择乘法模型进行因素分解。即：

$$x_t=T_t\times S_t\times I_t$$

4 趋势效应提取

4.1 复合移动平均法

  运行程序：

#2×12复合移动平均
m12=filter(data1/12,rep(1,12))
m2_12=filter(m12/2,rep(1,2),sides=1)  
m=data.frame(data1,m12,m2_12)
m[1:12,]

  运行结果：

##    data1      m12    m2_12
## 1    501       NA       NA
## 2    488       NA       NA
## 3    504       NA       NA
## 4    578       NA       NA
## 5    545       NA       NA
## 6    632 569.8333       NA
## 7    728 571.2500 570.5417
## 8    725 571.3333 571.2917
## 9    585 573.3333 572.3333
## 10   542 575.0833 574.2083
## 11   480 577.3333 576.2083
## 12   530 579.5833 578.4583

4.2 序列趋势效应

  运行程序：

#绘制移动平均效果图
plot(data1,lty=2,sub='图2 入住房间数序列趋势效应效果图')
lines(m2_12)

  运行结果：

4.3 去除序列趋势效应

  运行程序：

#绘制残差序列图
y_t=data1/m2_12
plot(y_t,sub='图3 入住房间数序列消除趋势效应效果图')

  运行结果：

  对原序列去除趋势效应后，效果图见图3，接下来提取季节效应及随机效应并查看其特征。

5 季节效应提取

5.1 序列矩阵化

  运行程序：

y_t=matrix(y_t,ncol=12,byrow=T)
y_t

  运行结果：

##            [,1]      [,2]      [,3]      [,4]      [,5]     [,6]     [,7]
##  [1,]        NA        NA        NA        NA        NA       NA 1.275980
##  [2,] 0.8930393 0.8403867 0.9041741 1.0212403 0.9709315 1.115059 1.244719
##  [3,] 0.9157786 0.8596083 0.8726676 1.0202661 0.9759282 1.109517 1.252512
##  [4,] 0.9197109 0.8605388 0.8903480 1.0161385 0.9606873 1.086375 1.221869
##  [5,] 0.8959224 0.8442748 0.8785510 1.0126904 0.9957624 1.088024 1.242261
##  [6,] 0.9227351 0.8142831 0.8768527 0.9573189 0.9898708 1.104536 1.275082
##  [7,] 0.9227468 0.8442783 0.8752808 0.9712128 0.9585319 1.088477 1.274262
##  [8,] 0.8967675 0.8334583 0.8282040 0.9738454 0.9672337 1.116629 1.268038
##  [9,] 0.8918692 0.8387205 0.8672606 0.9855716 0.9683954 1.091391 1.237236
## [10,] 0.9025797 0.8429484 0.8486417 0.9476738 0.9601540 1.069592 1.266978
## [11,] 0.9058259 0.8176428 0.8213450 0.9499142 0.9588015 1.103778 1.327046
## [12,] 0.9255660 0.8650297 0.8707124 0.9580547 0.9627160 1.080460 1.264940
## [13,] 0.8918033 0.8680851 0.8848580 0.9749484 0.9269225 1.101435 1.260261
## [14,] 0.9401082 0.8471814 0.8611058 0.9744083 0.9598156 1.071736       NA
##           [,8]      [,9]     [,10]     [,11]     [,12]
##  [1,] 1.269054 1.0221316 0.9439083 0.8330320 0.9162285
##  [2,] 1.270391 1.0062683 0.9792854 0.8262677 0.9244788
##  [3,] 1.258572 0.9791653 0.9673674 0.8480735 0.9435516
##  [4,] 1.290490 1.0014880 0.9965773 0.8508010 0.9318895
##  [5,] 1.257314 0.9768400 0.9893975 0.8732710 0.9592255
##  [6,] 1.274656 1.0164604 0.9834062 0.8281801 0.9408390
##  [7,] 1.312950 1.0026954 0.9980671 0.8409741 0.9435850
##  [8,] 1.336960 1.0091432 0.9911783 0.8335667 0.9368604
##  [9,] 1.308971 1.0267872 0.9977543 0.8446172 0.9545976
## [10,] 1.321614 1.0276278 1.0058133 0.8755799 0.9619668
## [11,] 1.285051 0.9995384 0.9682856 0.8458796 0.9542404
## [12,] 1.299224 0.9786572 0.9871404 0.8560296 0.9597128
## [13,] 1.313549 0.9808787 1.0082170 0.8341655 0.9438854
## [14,]       NA        NA        NA        NA        NA

5.2 序列总均值

  运行程序：

m=mean(y_t,na.rm = T)#计算序列总均值
m

  运行结果：

## [1] 0.9996314

5.3 序列月度均值

  运行程序：

ms=0
for(k in 1:12){
  ms[k]=mean(y_t[,k],na.rm = T)
}
ms

  运行结果：

##  [1] 0.9095733 0.8443412 0.8676924 0.9817910 0.9658270 1.0943851 1.2623986
##  [8] 1.2922152 1.0021293 0.9858768 0.8454183 0.9439278

5.4 序列季节指数

1.序列季节指数计算

  运行程序：

#求乘法模型的季节指数
S=ms/m
S

  运行结果：

##  [1] 0.9099087 0.8446526 0.8680124 0.9821531 0.9661832 1.0947887 1.2628642
##  [8] 1.2926918 1.0024989 0.9862404 0.8457301 0.9442759

2.序列季节指数可视化

  运行程序：

#绘制季节指数图
month=seq(1:12)
plot(month,S,type="o",sub='图4 入住房间数序列季节指数图',ylab="季节指数")

  运行结果：

6 随机效应提取

  运行程序：

#绘制随机效应示意图
I=data1/m2_12/S
plot(I,sub='图5 入住房间数序列随机效应示意图')

  运行结果：

  通过对原序列剔除趋势效应后，赋值给变量$y_t$，然后基于$y_t$分别求序列的总均值m和各季度均值ms。乘法模型中，各季度均值除以总均值得到季节指数S。从季节指数图（图4）可以清晰地看到，此欧洲小镇旅馆入住总序列具有6-8月为旺季的特征，随机波动序列特征如图5所示,在1附近无规律波动。