IBM SPSS ARIMA 算法使用案例介绍

余温半世容颜

IBM SPSS ARIMA 算法使用案例介绍
ARIMA 算法概述
时间序列是以规律的时间间隔采集的测量值的有序集合。时间序列分析的主要目的是根据现有的历史数据来预测未来的数据。
IBM SPSS ARIMA（差分自回归移动平均值）模型是一种典型的时间序列模型，可用于预测未来值，提供比指数平滑模型更复杂的趋势和季节性成分建模方法，在模型中包含预测变量。包括自回归模型（AR），移动平均模型（MA）和自回归移动平均模型（ARMA）。
Python Notebook 环境
可以通过下面两种 notebook 环境来进行预测分析：
本地环境
本地环境需要有 IBM SPSS 算法 jar 包及其依赖的包（如果没有的话可以用后面介绍的IBM Watson Studio 环境，除了加载输入数据不同其它代码是一样的）。
首先安装 Python 到路径 C:\Python35，这里 Python 版本是 3.5.2，选择自定义安装，指定安装路径，并加入系统 PATH 变量。如果下面三个库还没有安装的话从命令行执行 pip 安装：
pip install pyspark
pip install jupyter
pip install numpy
然后把下面 IBM SPSS 两个 jar 包放到指定路径 C:\Python35\Lib\site-packages\pyspark\jars，这里 com.ibm.spss.com.ibm.json.JSON4J-1.0.jar 是 IBM SPSS 算法包依赖的库，IBM SPSS 算法包这里使用的版本号是 5.1。
spss-assembly-spark2.0-scala2.11-5.1*.jar
com.ibm.spss.com.ibm.json.JSON4J-1.0.jar
解压文件 sparkspsspy-5.1*.zip, 然后把目录 spss 复制到下面的路径即可。 C:\Python\Python35\Lib\site-packages
任选工作目录，从命令行运行下面的命令就会在当前缺省浏览器中打开 Jupyter Notebook.
C:\temp\arima\jupyter notebook
如图 1：
图 1. Jupyter Home

右键点击 New 按钮，下拉菜单中选择 Python 3，如图 2：

接下来就可以在 notebook 的 cell 里面输入 Python 代码了。
IBM Watson Studio
IBM Watson Studio 已经集成了 IBM SPSS 算法包及依赖包。登录进去后，在 Projects 里面新建一个 Project，然后在 Assets->Notebooks->New notebook，选择 run time 的环境去创建 notebook，图 3 显示了内核是 Python 3.5 with Spark 2.1 的 notebook.
图 3. IBM Watson Studio Notebook

使用案例说明
某产品分类公司想要根据过去十年的销售数据来预测来年的男装季度销售，还想知道哪些预测变量（predictor）对目标被预测变量（target）影响比较大，从而做出相应的决策。预测变量包括：邮递的产品目录数（mail）和产品目录的页数（page）、开通的订购热线数目（phone）、印刷广告投入额（print）以及客户服务代表人数（service）。
这里的使用案例的输入数据 catalog_seasfac.csv 来自文件 catalog_seasfac.sav ， 可在任何 IBM® SPSS® Modeler 安装程序的 Demos 目录中找到。
输入原始数据的转换
首先加载数据 catalog_seasfac.csv 为 spark dataframe ，清单 1 显示了在本地 notebook 环境下的代码，csv 文件在当前 notebook 路径"C:\temp\arima\"下面。
清单 1. 本地环境加载输入数据为 Spark Data Frame

1. from pyspark.sql import SparkSession
   
2. spark = SparkSession.builder.getOrCreate()
   
3. from pyspark.sql.types import *
   
4. from pyspark import SparkContext
   
5. from pyspark import SQLContext
   
6. sc = SparkContext.getOrCreate();
   
7. sqlContext=SQLContext.getOrCreate(sc);
   
8. 
   
9. df = sqlContext.read.format('com.databricks.spark.csv') \
   
10.     .options(header='true', inferschema='true') \
    
11. .load('./catalog_seasfac.csv') # this is your csv file
    
12. print(df.schema)
    
13. df.show(5)

复制代码

如果是 IBM Watson Studio notebook, 需要把 catalog_seasfac.csv 上传使用。登录 IBM Watson Studio 后，新建一个 project，然后在 Assets 里面新建一个 notebook,在 notebook 右上方，点击"Find and add data"图标上传数据，然后右键选择已上传的数据 catalog_seasfac.csv ，选择 Insert to Code ->Insert SparkSession DataFrame，notebook 于是自动插入代码来创建 Spark Dataframe,这里也可以修改自动插入代码加入预先定义数据的类型 df_schema，见清单 2。
清单 2. 在 IBM Watson Studio 中加载输入数据为 Spark Data Frame

1. from pyspark.sql import SparkSession
   
2. from pyspark.sql.types import *
   
3. spark = SparkSession.builder.getOrCreate()
   
4. df_schema = StructType([
   
5.                 StructField("date",TimestampType(),True),
   
6.                 StructField("men",DoubleType(),True),
   
7.                 StructField("women",DoubleType(),True),
   
8.                 StructField("jewel",DoubleType(),True),
   
9.                 StructField("mail",DoubleType(),True),
   
10.                 StructField("page",DoubleType(),True),
    
11.                 StructField("phone",DoubleType(),True),
    
12.                 StructField("print",DoubleType(),True),
    
13.                 StructField("service",DoubleType(),True),
    
14.                 StructField("YEAR_",DoubleType(),True),
    
15.                 StructField("MONTH_",DoubleType(),True),
    
16.                 StructField("DATE_",StringType(),True),
    
17.                 StructField("Seasonal_Err_Men",DoubleType(),True),
    
18.                 StructField("Seasonal_AdjSer_Men",DoubleType(),True),
    
19.                 StructField("Seasonal_Factors_Men",DoubleType(),True),
    
20.                 StructField("Seasonal_TrendCycle_Men",DoubleType(),True)])
    
21. df = spark.read\
    
22.   .format('org.apache.spark.sql.execution.datasources.csv.CSVFileFormat')\
    
23.   .option('header', 'true')\
    
24.   .schema(df_schema)\
    
25.   .load(cos.url('catalog_seasfac.csv', '***********************************'))
    
26. print(df.schema)
    
27. df.show(5)

复制代码

图 4 是打印出的数据模式和前 5 行数据，从图中可以看到输入数据的所有字段及其类型，注意 date 这里是 TimestampType 类型。因为该产品分类公司希望预测来年的男装季度销售，而当前原始数据的时间间隔是月，所以需要对原始数据作转换把时间间隔转换为季度。
图 4. Data schema 和数据

清单 3 代码引入 TSDP, RTSDP 和 ARIMA 模型相关的包。
清单 3. TSDP, RTSDP 和 ARIMA 相关依赖包

1. from spss.ml.forecasting.timeseriesdatapreparation import TimeSeriesDataPreparation
   
2. from spss.ml.forecasting.traditional.timeseriesforecasting import TimeSeriesForecastingArima
   
3. from spss.ml.forecasting.reversetimeseriesdatapreparation import ReverseTimeSeriesDataPreparation
   
4. from spss.ml.forecasting.params.predictor import Predictor,ScorePredictor
   
5. from spss.ml.forecasting.params.targetorderlist import TargetOrderList
   
6. from spss.ml.forecasting.params.tforderlist import TFOrderList
   
7. from spss.ml.forecasting.params.temporal import ForecastEs
   
8. from spss.ml.common.wrapper import LocalContainerManager

复制代码

基于 IBM SPSS 当前对算法的设计,输入的原始数据必须由 IBM SPSS TSDP (Time Series Data Preparation) 先处理转换为一个合适的格式供 ARIMA 模型或其它 IBM SPSS 时间序列模型来消费,然后 IBM SPSS RTSDP (Reverse Time Series Data Preparation) 再将 ARIMA 或其它 IBM SPSS 时间序列模型的输出转换为同原始数据一样的格式。
在 IBM SPSS Modeler 和 IBM SPSS Statistics 中，TSDP 和 RTSDP 被嵌入在产品中对用户是透明的，所以不需要考虑数据转换的问题。但是在 notebook 中这种通过 API 调用来使用 IBM SPSS 时间序列模型的方式必须考虑数据前后的转换。
清单 4 代码即调用 TSDP，把输入数据的 dataframe 转换为供 ARIMA 消费的格式。这里可以通过 TSDP 对输入数据做不同的转换。
清单 4. 调用 TSDP

1. lcm = LocalContainerManager()
   
2. tsdp = TimeSeriesDataPreparation(lcm). \
   
3.     setMetricFieldList(["men","mail","page","phone","print","service"]). \
   
4.     setDateTimeField("date"). \
   
5.     setInputTimeInterval("MONTH"). \
   
6.     setOutTimeInterval("QUARTER")
   
7. tsdpOutput=tsdp.transform(df)

复制代码

LocalContainerManager 用来创建容器（container）的对象, 通常包含数据转换或建模过程中产生的 Json, PMML 或 StatXML 等与模型相关的文件。
setMetricFieldList 指定变量字段，这里的变量应该是连续变量，包括被预测变量字段和预测变量字段。
setDateTimeField 指定时间日期字段的名字，这里是 date。
setInputTimeInterval 指定输入数据的时间间隔是天（DAY），星期（WEEK），月（MONTH），季度（QUARTER）或年（YEAR），通过 setOutTimeInterval 指定转换后的时间间隔。本文指定的输入时间间隔是月（MONTH），输出时间间隔是季度（QUARTER）。这里的时间间隔根据原始数据的不同，还可以是时间间隔秒（SECOND），分钟（MINUTE）或小时（HOUR），甚至非标准时间间隔循环周期（CYCLIC_PERIOD）。
在选择时间序列的模型建模之前，通常需要分析一下时间序列曲线。清单 5 的代码用来把 TSDP 转换后的格式再通过 RTSDP 转换为和输入原始数据一样的格式，这样就可以把历史每个季度的男装销售曲线图画出来，便于进一步的分析。
清单 5. 调用 RTSDP

1. rtsdpRawdata = ReverseTimeSeriesDataPreparation(lcm). \
   
2.     setInputContainerKeys([tsdp.uid])
   
3. rtsdpRawdataOutput=rtsdpRawdata.transform(tsdpOutput)
   
4. print(rtsdpRawdataOutput.schema)
   
5. rtsdpRawdataOutput.show(5)

复制代码

从图 5 中可以看到，经过 TSDP 处理之后，输入数据的时间间隔已经从月变成了季度，这里对预测变量作了聚集，现在每个季度的值是之前三个月的值的和。
图 5. TSDP, RTSDP 转换后的数据

清单 6 代码画出过去十年（1989 - 1999）的男装季度销售额曲线图：
清单 6. 画出时间序列图

1. import matplotlib.pyplot as plt
   
2. import matplotlib.pyplot as plt
   
3. %matplotlib inline
   
4. fig = plt.figure(figsize=(16,8))
   
5. x=rtsdpRawdataOutput.select('date').rdd.map(lambda row : row[0]).collect()
   
6. y=rtsdpRawdataOutput.select('men').rdd.map(lambda row : row[0]).collect()
   
7. plt.plot(x, y,color="blue", linewidth=1.0, linestyle="-", label="men")
   
8. plt.xlabel('date')
   
9. plt.legend(loc=0)

复制代码

从图 6 中可以看到，该时间序列整体呈现上升趋势，即序列值趋向于随时间变化而增加，并维持线性趋势。此序列还有一个明显的季节模式，即年度高点在第四季度。
图 6. 时间序列图

因为该产品分类公司除了想预测来年每个季度的销售额，还想知道哪些变量对该销售额影响最大，所以考虑使用 IBM SPSS ARIMA 模型来进行预测分析。
ARIMA 算法建模
清单 7 的代码用来创建 ARIMA 模型。
清单 7. 创建 ARIMA 模型

1. arima = TimeSeriesForecastingArima(lcm). \
   
2.     setInputContainerKeys([tsdp.uid]). \
   
3.     setTargetPredictorList([Predictor(targetList = [["men"]],
   
4.                                      predictorIncludeList=[["mail"],["page"],["phone"],["print"],["service"]])]).\
   
5.     setTargetOrderList([TargetOrderList(
   
6.      targetList=[["men"]],
   
7.      nonSeasonal=[1,0,1],
   
8.      seasonal=[1,0,1],
   
9.      transType="none"
   
10.      )]). \
    
11.     setTFOrderList([TFOrderList(
    
12.      targetList=[["men"]],
    
13.      predictorList=[["mail"],["page"],["phone"],["print"],["service"]],   
    
14.      nonSeasonal=[1,0,0],
    
15.      seasonal=[0,0,0],
    
16.      transType="none"
    
17.      )]). \
    
18.     setOutInputData(True). \
    
19.     setCalcPI(True)
    
20. arima_model=arima.fit(tsdpOutput)

复制代码

setInputContainerKeys 用来设置 TSDP 转换后的容器相关信息。
setTargetPredictorList 用来设置被预测变量（target）和预测变量（predictor）。这里被预测变量就是要预测的男装销售额 men，预测变量包括 mail, page, phone, print 和 service。
setTargetOrderList 用来设置 ARIMA 的阶数，针对被预测变量(target)。对于 ARIMA 的阶数，所有的值都必须是非负整数。对于自回归和移动平均值来说，该值表示最大阶数。所有低于该值的正阶数都将包括在模型中。例如，如果指定 2，那么模型包括的阶数是 2 和 1。阶数依次包括下面三个值,本文中阶数的设置为实验值。
自回归（Autoregressive, p）：模型中的自回归阶数。自回归阶数指定序列中哪些以前的值用于预测当前值。例如，自回归阶数 1 指定序列中过去一个时限的值用于预测当前值。
差分（Difference, d）：指定在估计模型之前应用于序列的差分的阶。当趋势出现时（具有趋势的序列通常是不稳定的，而 ARIMA 建模时假定是稳定的），差分是必需的并可用于去除这些趋势的影响。差分阶数与序列趋势度一致，1 即一阶差分表示线性趋势，2 即二阶差分表示二次趋势，0 表示没有差分。
移动平均值（Moving Average, q）：模型中移动平均值阶数的值。移动平均值阶数指定如何使用与序列以前值均值之间的偏差来预测当前值。例如，移动平均值阶数 1 指定在预测序列的当前值时，可考虑与序列（来自过去一个时限中的每一个）均值之间的偏差。
setTFOrderList 用来设置转换函数的阶数,针对预测变量（predictor）。对于转换函数的阶数，所有的值都必须是非负整数。对于分子和分母来说，该值表示最大阶数。所有低于该值的正阶数都将包括在模型中。阶数依次包括下面三个值,本文中阶数的设置为实验值。
分子（Numerator）: 转换函数的分子阶数指定选定的独立（预测变量）序列中哪些先前的值用于预测相依序列的当前值。例如，分子阶数 1 指定独立序列过去一个时限的值以及独立序列的当前值可用于预测每个相关序列的当前值。
差分（Difference）:指定在估计模型之前应用于所选独立（预测变量）序列的差分的阶数。当趋势出现时，差分是必需的并可用于去除这些趋势的影响。
分母（Denominator）: 转换函数的分母阶数指定如何使用与选定独立（预测变量）序列的先前值均值之间的偏差来预测相依序列的当前值。例如，分母阶数 1 指定当预测每个相关序列的当前值时，可考虑与独立序列过去一个时限的均值之间的偏差。
nonSeasonal 指非季节性的阶数，取值可以参考上面阶数的定义。
seasonal 指季节性的阶数，取值可以参考上面阶数的定义。
注意：季节性自回归成分、移动平均值成分和差分成分与其非季节性对应成分起着相同的作用（季节性分子、分母和差分成分与其非季节性对应成分起着相同的作用）。但是对于季节阶数，当前的序列值会受到由一个或多个季节周期分隔的以前序列值的影响。例如，对于以月为时间单位的数据（季节周期为 12），季节阶数 1 表示当前序列值会受到当前序列之前的 12 个周期内的序列值的影响。因此，对于以月为时间单位数据，将季节阶数指定为 1 相当于将非季节阶数指定为 12。
setOutInputData 用来将输入数据包含在预测的输出中，从而可以比较输入值和输出预测值。
setCalcPI 用来计算预测变量的重要性,只有在存在预测变量的时候才有效。
ARIMA 模型解析及预测变量重要性
建模后还需要对模型相关文件进行解析，清单 8 的代码用来输出容器中所有模型相关的文件的列表：
清单 8. 输出模型相关文件列表

1. arima_container = lcm.importContainerSeq(arima.uid)
   
2. arima_container.entriesNames()

复制代码

图 7 为输出文件列表的例子，这里 TSDPOutput.json 是 TSDP 的输出，datamodel.xml 是数据模型文件，containers[0]中的 *.es.xml 和 0.xml 是 PMML 文件，containers[1]中的 StatXML.xml 和 0.xml 都是 StatXML。StatXML.xml 中可以看到构建模型的各种设置及其值，0.xml 可以查看模型的各个统计量的信息进而判断模型的好坏。
图 7. 模型相关文件列表

清单 9 中的代码输出 PMML 文件，数据模型文件和 StatXML 文件，其它文件可仿照例子输出。

1. containers = arima_container.containers();
   
2. print("containers[0] 0.xml:\n"+containers[0].containerEntry('0.xml').stringContent()+"\n")
   
3. print("containers[0] datamodel.xml:\n"+containers[0].containerEntry('datamodel.xml').stringContent()+"\n")
   
4. print("containers[1] StatXML.xml:\n"+containers[1].containerEntry('StatXML.xml').stringContent()+"\n")
   
5. print("containers[1] 0.xml:\n"+containers[1].containerEntry('0.xml').stringContent()+"\n")

复制代码

从 containers[0]中的 PMML 文件 0.xml 可以得到各个预测变量的重要性，见图 8：
图 8. 预测变量重要性

从数据模型 datamode.xml 中可以得到对应的变量的编码，１对应的是 mail，3 对应的是 phone，见图 9。因为这两个预测变量的 importance 值最大，所以这两个预测变量对被预测变量 men 的影响最大。
图 9. 数据模型中的变量

从 containers[1]中的 StatXML 文件 0.xml 中可以看到模型相关的统计量，可以对模型的好坏做出评估，见图 10：
图 10. 模型相关统计量

ARIM 算法预测
清单 10 的代码用来作预测及设置预测参数:
清单 10. 预测及参数设置

1. transformer = arima_model. \
   
2.   setTargets(ScorePredictor(targetIncludedList=[["men"]])). \
   
3.   setForecast(ForecastEs(outForecast = True, forecastSpan = ４,outCI = True))
   
4. predictions = transformer.transform(tsdpOutput)

复制代码

setTargets 用来指定被预测变量 men。
setForecast 用来指定预测相关设置，outForecast 指定是否作预测，forecastSpan 指定预测的时间间隔，这里时间间隔是季度，所以 4 指定预测来年的 4 个季度，outCI 指定是否计算输出预测置信区间 LCI (Lower Confidence Interval) 和 UCI (Upper Confidence Interval)。
清单 11 的代码用来将预测结果通过 RTSDP 转换为和输入原始数据一样的格式。对于输出的 dataframe, 为了便于显示没有选取所有字段而通过 select 选取了部分字段。
清单 11. RTSDP 转换

1. rtsdp = ReverseTimeSeriesDataPreparation(lcm). \
   
2.     setInputContainerKeys([tsdp.uid]). \
   
3.     setDeriveFutureIndicatorField(True)
   
4. 
   
5. score = rtsdp.transform(predictions).select("date","$FutureFlag","men","$TS-men","$TSLCI-men","$TSUCI-men")
   
6. score.show(score.count())

复制代码

setInputContainerKeys 用来设置 TSDP 转换后的容器相关信息。
setDeriveFutureIndicatorField 是否输出未来指示字段。
图 11 为预测的输出，这里$FutureFlag 列为 0 显示历史时间，为 1 显示未来的时间。$TS-men 列为男装销量的预测值，$TSLCI-men 列和$TSUCI-men 列表明了预测值的置信区间。
图 11. 预测输出

清单 12 的代码用来画出男装销量的历史值和预测值的比较图：
清单 12. 画出预测图

1. from datetime import datetime,time
   
2. x=score.select('date').rdd.map(lambda row : row[0]).collect()
   
3. yObserve=score.select('men').rdd.map(lambda row : row[0]).collect()
   
4. yPredict=score.select('$TS-men').rdd.map(lambda row : row[0]).collect()
   
5. 
   
6. fig = plt.figure(figsize=(16,8))
   
7. plt.plot(x, yObserve,color="blue", linewidth=1.0, linestyle="-", label="men")
   
8. plt.plot(x, yPredict,color="red", linewidth=1.0, linestyle="-", label="$TS-men")
   
9. plt.xlabel('date')
   
10. plt.legend(loc=0)

复制代码

图 12 展示了男装季度销售额 men 和其预测值的比较图$TS-men，从图中可以看到 1999 年的整体销售情况预测较好，1998 年第四季度高峰期后 1999 年上半年再次回复正常，下半年一直保持平稳的上升趋势。
图 12. 预测图

结束语
本文通过在 notebook 中转换原始数据，设置模型参数，创建模型，并根据预测变量的重要性（Predictor Importance）找到对预测影响最大的预测变量，做出预测，完整地展示了如何通过 Python 调用 IBM SPSS ARIMA API 进行建模、预测及分析。