关于入境旅游人数的时间序列分析.doc

关于入境旅游人数的时间序列分析

专业：统计学 姓名：佟虹生 指导老师：汪小英

摘要
大众旅游时代的到来，使旅游日益成为现代人类社会主要的生活方式和社会经济活动，旅游业以其强劲的势头成为全球经济产业中最具活力的“朝阳产业”。随着社会生产力不断发展，劳动生产率不断提高，以及人们生活水平的迅速提高和带薪假期的增加，旅游业将持续高速度发展，成为世界最重要的经济部门之一。中国同样如此，有数据统计2014年全国旅游及相关产业增加值占国内生产总值的比重为4.33%。
众所周知旅游业是一个存在显著季节效应的行业，如果能对旅游业的客流量作出准确的预测，将会有利于商家更好的把握商机。本文选取入境旅游的客流量作为时间序列，将简单地分析该序列的季节效应，并对序列拟合ARIMA模型，并作出简单的预测。
关键词：入境旅游；季节效应；时间序列；ARMIA；预测

一、 引言

旅游是在人的基本生活需求得到适度满足后的一种新的消费行为，一种带有浓厚文化内涵的群体活动。人们离开常住地到异国他乡访问的旅行和暂时停留所引起的各种现象和关系的总和。
我国拥有丰富的旅游资源，疆域辽阔，既有风景秀丽的江南水乡，也有粗犷豪迈的西北风情；我国拥有悠久的历史文化，目前已经公布了99个国家级历史文化古城，长城、故宫、颐和园等已经被列入世界文化遗产名录；我国还是一个拥有多个民族的国家，各个民族的习俗和风情很容易使人产生很强烈的向往之情。所有这些，都为我国旅游业的发展奠定了一个良好的基础，使得我国吸引了大量的入境游客。
入境旅游是指他国居民前来我国的旅游活动，或者是指他国居民进入本国国境以内的旅游活动，入境旅游属于国际旅游。目前入境旅游已成为构成我国旅游业的重要组成部分。
时间序列分析(Time series analysis)是一种动态数据处理的统计方法。该方法基于随机过程理论和数理统计学方法，研究随机数据序列所遵从的统计规律，以用于解决实际问题。因此本文将通过时间序列分析的方法来研究分析中国的入境旅游人数。

二、研究方法

差分运算具有强大的确定性信息提取能力，许多非平稳序列差分后会显示出平稳序列的性质，这是我们成这个非平稳序列为差分平稳序列。对差分平稳序列可以使用ARIMA模型进行拟合。
具有如下结构的模型称为求和自回归移动平均（autoregressive integrated moving average），简记为ARIMA(p，d，q)模型：
ΦB∇dxt=ΘBεt Εεt=0，Varεt=σε2，Εεtεs=0，s≠tΕxsεt=0，∀s<t
式中，∇d=(1-B)d；ΦB=1-ϕ1B-⋯-ϕpBp，为平稳可逆ARMA（p，q）模型的自回归系数多项式；ΘB=1-θ1B-⋯-θqBq，为可逆ARMA（p，q）模型的移动平滑系数多项式。
ARIMA（p，d，q）模型是指d阶差分后自相关最高阶数为p，移动平均最高阶数为q的模型，通常它包含了p+q个独立的未知系数：ϕ1，⋯，ϕp，θ1，⋯，θq。
如果该模型中有部分自相关系数ϕj(1≤j<p)或部分移动平滑系数θk(1≤k<q)为零，即原ARIMA（p，d，q）模型中有部分系数省缺了，那么该模型称为疏系数模型。
疏系数模型一般形式为
ARIMA(p1,⋯,pm,d,(q1,⋯,qn))
式中，p1，⋯，pm为非零自相关系数的阶数，q1，⋯，qn为非零移动平滑系数的阶数。
在实际操作中，疏系数模型时有应用。
ARIMA模型可以对具有季节效应的序列建模。根据季节效应提取的难易程度，可以分为简单季节模型和乘积季节模型。
简单季节模型是指序列中的季节效应和其他效应之间是加法关系，即
xt=St+Tt+It
式中，Tt代表序列的长期趋势波动；St代表序列的季节性（周期性）变化；It代表随机波动。
这时，各种效应信息的提取都非常容易。通常简单的周期步长差分即可将序列中的季节信息提取充分，简单的低阶差分即可将趋势信息提取充分，提取完季节信息和趋势信息之后的序列就是一个平稳序列，可以用ARMA模型拟合。
所以简单季节模型实际上就是通过趋势差分、季节差分将序列转换为平稳序列，再对其进行拟合。它的模型结构通常如下：
∇D∇dxt=Θ(B)Φ(B)εt
式中，
（1） D为周期步长，d为提取趋势信息所用的差分阶数。
（2） εt为白噪声序列，且Εεt=0，Varεt=σε2。
（3） ΘB=1-θ1B-⋯-θqBq，为q阶移动平均系数多项式。
（4） ΦB=1-ϕ1B-⋯-ϕpBp，为p阶自回归系数多项式。
但更为常见的情况是，序列的季节效应、长期趋势和随机波动之间有着复杂德尔交互影响关系，简单的ARIMA模型并不足以提取其中的相关关系，这时通常需要采用乘积季节模型。
乘积模型的构造原理如下：
当序列具有短期相关性时，通常可以使用低阶ARMA（p，q）模型提取。
当序列具有季节效应，季节效应本身还具有相关性时，季节相关性可以使用以周期为步长单位的ARMA（P，Q）模型提取。
由于短期相关性和季节效应之间具有乘积关系，所以拟合模型实质为ARMA（p，q）和ARMA（P，Q）的乘积。综合前面的d阶趋势差分和D阶以周期S为步长的季节差分运算，对原观察值序列拟合的乘积模型完整的结构如下：
∇d∇SDxt=Θ(B)ΘS(B)Φ(B)ΦS(B)εt
式中，ΦB=1-ϕ1B-⋯-ϕpBp
ΘB=1-θ1B-⋯-θqBq
ΦSB=1-ϕ1BS-⋯-ϕPBPS
ΘSB=1-θ1BS-⋯-θQBQS
该乘积模型简记为ARIMA(p,d,q)×(P,D,Q)S。

三、对中国入境旅游人数的实证研究

本文实证研究所选取的数据是中国入境旅游人数，数据频率为月度，客流量单位为万人，时间跨度从2001年1月至2015年12月，共计180个数据，来源于中经网统计数据库。

3.1原序列以及差分序列的相关检验
3.1.1原序列的平稳检验

图3.1 原始序列时序图
首先画出原始序列的的时序图，进行观察，并结合单位根检验判断它的平稳性，图3.1为原始序列的时序图，表3.1为原始序列的3阶单位根检验结果。
表3.1 原始序列3阶单位根检验结果

增广 Dickey-Fuller 单位根检验
类型
滞后
Rho
Pr < Rho
Tau
Pr < Tau
F
Pr > F
零均值
0
-0.1237
0.6538
-0.12
0.6421

1
0.1717
0.7223
0.22
0.7505

2
0.3017
0.7548
0.63
0.8520

3
0.3173
0.7587
0.78
0.8810

单均值
0
-28.4429
0.0013
-4.02
0.0017
8.17
0.0010

1
-19.9134
0.0102
-3.44
0.0109
6.17
0.0088

2
-8.6838
0.1794
-2.31
0.1684
3.13
0.2727

3
-6.7822
0.2846
-2.08
0.2537
2.73
0.3758
趋势
0
-63.3460
0.0005
-6.19
<.0001
19.16
0.0010

1
-49.0645
0.0005
-5.06
0.0003
12.96
0.0010

2
-21.0121
0.0498
-3.20
0.0874
5.26
0.1262

3
-16.1923
0.1367
-2.77
0.2117
3.98
0.3823

从时序图中可以观察到这是一个典型的有着上升趋势的非平稳序列，同时单位根检验也印证了我们的判断，从2阶开始p值大于0.05说明原始序列非平稳。同时从时序图中，我们还可以观察到原始序列在上升趋势的同时，还具备明显的周期性，周期长度为12，与常识的认知相符合。所以接下来要考虑对原始序列进行1阶12步差分处理，1阶差分去除原始序列的趋势，1阶差分后的序列进行12步差分提取差分后序列的季节信息即周期性。

3.1.2差分序列的相关检验

图3.2 差分后序列的时序图
我们对1阶12步差分以后的序列进行平稳性和白噪声检验，图3.2为差分后序列的时序图，表3.2为差分后序列的3阶单位根检验结果，表3.3为差分后序列的白噪声检验结果。
表3.2 差分后序列的3阶单位根检验结果

增广 Dickey-Fuller 单位根检验
类型
滞后
Rho
Pr < Rho
Tau
Pr < Tau
F
Pr > F
零均值
0
-212.748
0.0001
-17.23
<.0001

1
-209.293
0.0001
-10.18
<.0001

2
-337.999
0.0001
-9.06
<.0001

3
-697.217
0.0001
-8.03
<.0001

单均值
0
-212.743
0.0001
-17.18
<.0001
147.63
0.0010

1
-209.295
0.0001
-10.15
<.0001
51.51
0.0010

2
-338.006
0.0001
-9.03
<.0001
40.78
0.0010

3
-697.566
0.0001
-8.01
<.0001
32.07
0.0010
趋势
0
-212.752
0.0001
-17.13
<.0001
146.76
0.0010

1
-209.397
0.0001
-10.12
<.0001
51.22
0.0010

2
-338.152
0.0001
-9.00
<.0001
40.53
0.0010

3
-700.571
0.0001
-7.99
<.0001
31.89
0.0010
表3.3 差分后序列的白噪声检验结果

白噪声的自相关检查
至滞后
卡方
自由度
Pr > 卡方
自相关
6
20.15
6
0.0026
-0.281
0.085
-0.155
0.001
-0.073
0.053
12
67.60
12
<.0001
-0.038
0.022
0.053
-0.017
0.155
-0.481
18
92.44
18
<.0001
0.193
-0.191
0.177
-0.066
0.117
-0.104
24
98.50
24
<.0001
0.083
-0.041
0.047
-0.115
0.079
-0.032

从1阶12步差分后的序列时序图中不难看出，已经充分地提取了原始序列的趋势，同时差分后的序列也没有明显的周期性，当然差分后序列的3阶单位根检验也印证了我们的判断。再观察差分后序列的白噪声检验结果，发现无论是6阶、12阶、18阶或是24阶，p值都小于0.05，说明差分后的序列为一个非白噪声序列，所以综上，我们可以判断原始序列经过1阶12步差分以后得到一个非白噪声的平稳序列，有分析的价值，考虑到季节效应，可以尝试对差分以后的序列拟合ARIMA模型。

3.2拟合ARIMA模型
3.2.1考虑简单季节模型

观察1阶12步差分后序列的自相关系数和偏自相关系数图，图3.3为差分后序列的自相关系数，图3.4为差分后序列的偏自相关系数。

图3.3 差分后序列的自相关系数

图3.4 差分后序列的偏自相关系数

自相关图显示延迟12阶自相关系数显著大于2倍标准差范围，这说明差分后序列中仍然蕴含着非常显著的季节效应。延迟1阶、2阶的自相关系数也大于2倍标准差，这说明差分后的序列还具有短期相关性。观察偏自相关图，得到结论基本一致。根据自相关系数和偏自相关系数的性质，可以尝试拟合疏系数模型AR（1,12）、MA（1,2,12）、ARMA（（1,12），（1,2,12）），表3.4、表3.5、表3.6分别为3种模型的残差白噪声检验结果。
表3.4 AR（1,12）的残差白噪声检验结果

残差的自相关检查
至滞后
卡方
自由度
Pr > 卡方
自相关
6
5.34
4
0.0017
-0.016
-0.078
-0.120
-0.072
-0.064
0.031
12
12.94
10
0.0753
-0.055
0.008
0.040
-0.018
0.111
-0.157
表3.5 MA（1,2,12）的残差白噪声检验结果

残差的自相关检查
至滞后
卡方
自由度
Pr > 卡方
自相关
6
6.24
3
0.0013
-0.067
0.020
-0.136
-0.058
-0.076
0.060
12
11.38
9
0.0802
-0.044
-0.011
-0.009
-0.096
0.041
-0.124

表3.6 ARMA（（1,12），（1,2,12））的残差白噪声检验结果

残差的自相关检查
至滞后
卡方
自由度
Pr > 卡方
自相关
6
10.38
1
0.0013
-0.136
0.081
-0.156
-0.036
-0.081
0.057
12
28.54
7
0.0002
-0.039
0.005
0.025
-0.085
0.072
-0.292

从以上模型残差的检验结果可以看出，三种模型的残差在滞后6阶的白噪声检验中p值都小于0.05，都不能通过，说明残差在短期还有信息没有提取完全，所以模型拟合效果均不理想。考虑到该序列既具有短期相关性又具有季节效应，而且短期相关性和季节效应使用加法模型无法充分、有效提取，可以认为该序列的季节效应和短期相关性之间具有复杂的关联性。这时，假定短期相关性和季节效应之间具有乘积关系，尝试使用乘积模型来拟合序列的发展。

3.2.2考虑乘积季节模型

还是观察1阶12步差分之后序列的自相关系数和偏自相关系数图，两者显示12阶以内的自相关系数和偏自相关系数均不截尾，所以尝试使用ARMA（1,1）模型提取差分后序列的短期自相关信息。
在考虑季节自相关特征，这时考察延迟12阶、24阶等以周期长度为单位的自相关系数和偏自相关系数的特征，自相关系数图（图3.3）显示延迟12阶自相关系数显著非零，但是延迟24阶自相关系数落入2倍标准差范围。而偏自相关系数图（图3.4）显示延迟12阶和延迟24阶的偏自相关系数都显著非零。所以可以认为季节自相关特征是自相关系数截尾，偏自相关系数拖尾，这时以12步为周期的ARMA（0,1）12模型提取差分后序列的季节自相关信息。
综合前面的差分信息，要拟合的乘积模型为ARIMA（1,1,1）×（0,1,1）12
∇∇12=1-θ1B1-ϕ1B1-θ12B12εt
表3.7为该模型条件最小二乘下的参数估计。
表3.8 参数估计

条件最小二乘估计
参数
估计
标准
误差
t 值
近似
Pr > |t|
滞后
MU
-0.19927
0.30542
-0.65
0.5150
0
MA1,1
0.84881
0.07839
10.83
<.0001
1
MA2,1
0.80888
0.05088
15.90
<.0001
12
AR1,1
0.57604
0.12083
4.77
<.0001
1

可以看到，该模型除去常数项不显著外，其余系数都显著非零，所以重新拟合不含常数项的该模型，并观察残差序列的白噪声检验。表3.8为不含常数项的参数估计，表3.9为不含常数项下的残差序列白噪声检验结果。
表3.8 不含常数项的参数估计

条件最小二乘估计
参数
估计
标准
误差
t 值
近似
Pr > |t|
滞后
MA1,1
0.83316
0.08461
9.85
<.0001
1
MA2,1
0.80609
0.05095
15.82
<.0001
12
AR1,1
0.56010
0.12669
4.42
<.0001
1

表3.9 不含常数项的残差白噪声检验结果

残差的自相关检查
至滞后
卡方
自由度
Pr > 卡方
自相关
6
2.22
3
0.5278
-0.014
0.034
-0.014
-0.014
-0.028
0.101
12
4.35
9
0.8870
-0.042
-0.007
0.026
-0.073
0.063
0.009
18
13.76
15
0.5439
0.027
-0.180
0.081
-0.090
0.053
0.009

从表3.8可以看出不含常数项的模型中所有系数均显著，表3.9则说明不含常数项的模型残差序列均不拒绝为白噪声序列，所以模型拟合成功，得到最终的模型的口径为：
∇∇12=1-0.83316B1-0.56010B1-0.80609B12εt

图3.5 拟合值和观察值联合图
将序列拟合值和序列观察值联合作图。图3.5为拟合值和观察值联合图。
图中分别作出模型拟合值以及模型拟合在95%置信水平下的置信上下限，均以折现连接，图中星星标记未连接的即为实际观测值，可以直观地看出该季节乘积模型对原序列拟合效果良好。

3.3预测

利用该模型对2016年的入境旅游人数进行预测，并绘制预测序列的时序图。表3.10为2016年每个月度入境旅游人数的预测值（单位：万人）。

表3.10 2016年每个月度入境旅游人数预测（单位：万人）

以下变量的预测:x
date
预测
95% 置信限
date
预测
95% 置信限
2016年1月
1112.679
1030.0292
1195.3287
2016年7月
1141.4522
1007.6716
1275.2327
2016年2月
986.6735
884.4935
1088.8535
2016年8月
1163.4328
1025.8239
1301.0417
2016年3月
1159.3093
1046.6533
1271.9652
2016年9月
1131.0711
989.826
1272.3163
2016年4月
1191.4926
1071.8249
1311.1603
2016年10月
1188.1258
1043.3837
1332.868
2016年5月
1124.4214
999.3415
1249.5012
2016年11月
1140.8166
992.6861
1288.9472
2016年6月
1099.2116
969.5473
1228.8759
2016年12月
1181.1752
1029.7462
1332.6041

从表中不难看出关于2016年中国的入境旅游预测，2月是一个淡季，4月、10月、12月将是旺季，而且整体来看2016年下半年的形式将好于上半年。

四、 结论

本文最终选择的基于季节乘积模型的ARIMA(1,1,1)×(0,1,1)12可以较好的拟合中国入境旅游人数，所以商家可以使用该模型对未来中国的入境旅游人数进行预测，把握住商机选择在旺季进行相关的投资，以牟取利润。相关管理部门也可以依据预测值进行管理制度的规划和制定。

参考文献
[1]John C.hull.Risk Management and Financial Institutions[M].3rd,New York:Wiley,2012.257-436.
[2]王燕.应用时间序列分析[M].3rd，北京:中国人民大学出版社，2012.150-202.
[3]贾俊平.统计学[M].5th,北京:中国人民大学出版社,2014.71-121.
[4]朱世武.SAS编程技术教程[M].2nd,北京:清华大学出版社,2013.83-392.
[5]朱世武.金融计算与建模[M].北京:清华大学出版社,2007.15-64.
[6]George E.P.Box,Gwilym M.Jenkins,Gregory C.Reinsel.Time Series Analysis:Forecasting and Control[M].4th,New York:Wiley,2011,123-187.

附录
原始数据：

date
number
date
numebr
date
number
2001年1月
717.15
2003年7月
776.94
2006年1月
998.85
2001年2月
611.41
2003年8月
884.4
2006年2月
871.44
2001年3月
753.26
2003年9月
807.58
2006年3月
1002.98
2001年4月
779.7
2003年10月
854.99
2006年4月
1097.09
2001年5月
707.87
2003年11月
828.42
2006年5月
1002.39
2001年6月
711.78
2003年12月
876.68
2006年6月
1000.45
2001年7月
745.38
2004年1月
808.73
2006年7月
1090.4
2001年8月
809.14
2004年2月
753.3
2006年8月
1114.77
2001年9月
741.32
2004年3月
855.16
2006年9月
1043.77
2001年10月
757.73
2004年4月
954.71
2006年10月
1138.31
2001年11月
749.71
2004年5月
877.47
2006年11月
1042.09
2001年12月
815.09
2004年6月
893.39
2006年12月
1092.51
2002年1月
740.7
2004年7月
959.03
2007年1月
1022.98
2002年2月
711.04
2004年8月
971.56
2007年2月
933.09
2002年3月
839.54
2004年9月
917.45
2007年3月
1089.8
2002年4月
807.77
2004年10月
988.18
2007年4月
1151.68
2002年5月
787.92
2004年11月
935.21
2007年5月
1072.55
2002年6月
795.77
2004年12月
989.63
2007年6月
1066.68
2002年7月
849.37
2005年1月
938.06
2007年7月
1149.14
2002年8月
890.68
2005年2月
855.85
2007年8月
1156.82
2002年9月
865.12
2005年3月
1027.77
2007年9月
1123.13
2002年10月
869.57
2005年4月
1024.76
2007年10月
1160.14
2002年11月
841.91
2005年5月
995.14
2007年11月
1109.1
2002年12月
866.12
2005年6月
989.42
2007年12月
1152.22
2003年1月
848.43
2005年7月
1076.11
2008年1月
1080.95
2003年2月
737.6
2005年8月
1067.41
2008年2月
990.77
2003年3月
785.14
2005年9月
988.77
2008年3月
1188.44
2003年4月
564.92
2005年10月
1054.79
2008年4月
1126.51
2003年5月
543.6
2005年11月
990.78
2008年5月
1087.5
2003年6月
652.6
2005年12月
1018.31
2008年6月
1017.21

date
number
date
number
date
number
2008年7月
1105.72
2011年1月
1082.77
2013年7月
1061.56
2008年8月
1076.24
2011年2月
961.05
2013年8月
1080.65
2008年9月
1056.4
2011年3月
1133.78
2013年9月
1053.59
2008年10月
1122.82
2011年4月
1227.67
2013年10月
1135.59
2008年11月
1055.09
2011年5月
1122.18
2013年11月
1079.28
2008年12月
1095.09
2011年6月
1099.69
2013年12月
1085.95
2009年1月
1033.19
2011年7月
1176.52
2014年1月
1062.86
2009年2月
939.81
2011年8月
1156.28
2014年2月
889.95
2009年3月
1054.07
2011年9月
1109.65
2014年3月
1065.4
2009年4月
1151.44
2011年10月
1176.39
2014年4月
1125.06
2009年5月
1033
2011年11月
1126.27
2014年5月
1062.83
2009年6月
994.72
2011年12月
1170.13
2014年6月
1025.28
2009年7月
1072.04
2012年1月
1060.89
2014年7月
1061.36
2009年8月
1109.19
2012年2月
1020.52
2014年8月
1095.54
2009年9月
1006.37
2012年3月
1150.18
2014年9月
1076.48
2009年10月
1113.54
2012年4月
1178.94
2014年10月
1131.87
2009年11月
1039.56
2012年5月
1097.61
2014年11月
1094.51
2009年12月
1100.67
2012年6月
1080.64
2014年12月
1158.7
2010年1月
1069.71
2012年7月
1105.63
2015年1月
1092
2010年2月
940.34
2012年8月
1128.29
2015年2月
927.13
2010年3月
1132.85
2012年9月
1095.37
2015年3月
1118.58
2010年4月
1186.15
2012年10月
1120.65
2015年4月
1199.27
2010年5月
1134.65
2012年11月
1083.11
2015年5月
1120.6
2010年6月
1088.3
2012年12月
1118.69
2015年6月
1078.62
2010年7月
1150.13
2013年1月
1079.87
2015年7月
1092.15
2010年8月
1171.41
2013年2月
943.08
2015年8月
1119.87
2010年9月
1103.69
2013年3月
1185.39
2015年9月
1133.46
2010年10月
1163.64
2013年4月
1081.34
2015年10月
1181.52
2010年11月
1089.34
2013年5月
1065.25
2015年11月
1146.67
2010年12月
1146.02
2013年6月
1056.23
2015年12月
1171.9

SAS实现过程：
proc import datafile="G:\trip.xlsx" dbms=excel out=a;/*2001年01月-2015年12月*/
sheet="sheet1";
getnames=yes;
run;
data a;
set a;
rename date=t number=x;
run;
proc gplot data=a;
plot x*t;
symbol i=join;
run;
quit;
proc arima data=a;
identify var=x stationarity=(adf=3);
run;
quit;
data b;
set a;
x_1=dif(x);
keep x_1 t;
run;
proc gplot data=b;
plot x_1*t;
symbol i=join;
run;
quit;
proc arima data=b;
identify var=x_1 stationarity=(adf=3);
run;
quit;
data c;
set b;
x_1_12=dif12(x_1);
keep x_1_12 t;
run;
proc gplot data=c;
plot x_1_12*t;
symbol i=join;
run;
quit;
proc arima data=c;
identify var=x_1_12 stationarity=(adf=3);
run;
quit;

proc arima data=c;
identify var=x_1_12;
estimate p=(1,12);
run;
quit;
proc arima data=c;
identify var=x_1_12;
estimate q=(1,2,12);
run;
quit;
proc arima data=c;
identify var=x_1_12;
estimate p=(1,12) q=(1,2,12);
run;
quit;

proc arima data=c;
identify var=x_1_12;
estimate p=1 q=(1)(12);/*p=(x,y)指疏系数模型，p=(x)(y)指第一步是AR(x)第二步是AR(y)*/
run;
quit;

proc arima data=c;
identify var=x_1_12;
estimate p=1 q=(1)(12) noint;
run;
quit;

proc arima data=a;
identify var=x(1,12);
estimate p=1 q=(1)(12) noint;
forecast lead=0 id=t out=p;
run;
quit;
proc gplot data=p;
plot x*t=1 forecast*t=2 l95*t=3 u95*t=3/overlay;
symbol1 c=black v=star i=none;
symbol2 c=red v=none i=join;
symbol3 c=green v=none i=join;
run;
quit;

proc arima data=a;
identify var=x(1,12);
estimate p=1 q=(1)(12) noint;
forecast lead=12 id=t;
run;
quit;