本次比赛目标：通过时间序列模型，预测接下来一个月,俄罗斯某商超集团每件商品在各个商店的总销售额。
预测了销售额又如何呢？当然作用大大，例如，良性的库存备货提高周转率，发现异常销售以改进，制定成本预算项目，审时度势因地制宜制定有效营销策略，等等。
# 数据集大小：最大数据文件含46万条数据。
# 比赛链接：https://www.kaggle.com/c/competitive-data-science-predict-future-sales
# 本次尝试：
# 先参考kaggle大神分享 https://www.kaggle.com/jagangupta/time-series-basics-exploring-traditional-ts
# 然后遇到Facebook开源工具prophet下载不下来，加上对时间序列模型的使用尚不熟悉，就转为参考书籍《Python数据科学·技术详解与商业实践》中的第十八章“时间序列建模”。
# 所以，全文内容结构是这样的：
# 一、Kaggle商品销售数据初探索。
# 二、《Python数据科学》时间序列模型温习，含定义说明和实例演示。
# 三、回到kaggle商品销售数据预测，将第二部分所学运用于最终的预测。
# Basic packages
import numpy as np # linear algebra
import pandas as pd # data processing, CSV file I/O (e.g. pd.read_csv)
import random as rd # generating random numbers
import datetime # manipulating date formats
# Viz
import matplotlib.pyplot as plt # basic plotting
import seaborn as sns # for prettier plots
# TIME SERIES
from statsmodels.tsa.arima_model import ARIMA
from statsmodels.tsa.statespace.sarimax import SARIMAX
from pandas.plotting import autocorrelation_plot
from statsmodels.tsa.stattools import adfuller, acf, pacf,arma_order_select_ic
import statsmodels.formula.api as smf
import statsmodels.tsa.api as smt
import statsmodels.api as sm
import scipy.stats as scs
# settings
import warnings
warnings.filterwarnings("ignore")
# Import all of them
sales=pd.read_csv("D:/2018_BigData/Python/Kaggle_learning/Predict Future Sales - data science/data-sources/sales_train.csv")
item_cat=pd.read_csv("D:/2018_BigData/Python/Kaggle_learning/Predict Future Sales - data science/data-sources/item_categories.csv")
item=pd.read_csv("D:/2018_BigData/Python/Kaggle_learning/Predict Future Sales - data science/data-sources/items.csv")
shops=pd.read_csv("D:/2018_BigData/Python/Kaggle_learning/Predict Future Sales - data science/data-sources/shops.csv")
test=pd.read_csv("D:/2018_BigData/Python/Kaggle_learning/Predict Future Sales - data science/data-sources/test.csv")
# 源数据文件描述
# sales_train.csv - the training set. Daily historical data from January 2013 to October 2015.
# test.csv - the test set. You need to forecast the sales for these shops and products for November 2015.
# sample_submission.csv - a sample submission file in the correct format.
# items.csv - supplemental information about the items/products.
# item_categories.csv - supplemental information about the items categories.
# shops.csv- supplemental information about the shops.
# 数据标签说明
# ID - an Id that represents a (Shop, Item) tuple within the test set
# shop_id - unique identifier of a shop
# item_id - unique identifier of a product
# item_category_id - unique identifier of item category
# item_cnt_day - number of products sold. You are predicting a monthly amount of this measure
# item_price - current price of an item
# date - date in format dd/mm/yyyy
# date_block_num - a consecutive month number, used for convenience. January 2013 is 0, February 2013 is 1,..., October 2015 is 33
# item_name - name of item
# shop_name - name of shop
# item_category_name - name of item category
sales.head()
| date | date_block_num | shop_id | item_id | item_price | item_cnt_day | |
|---|---|---|---|---|---|---|
| 0 | 02.01.2013 | 0 | 59 | 22154 | 999.00 | 1.0 |
| 1 | 03.01.2013 | 0 | 25 | 2552 | 899.00 | 1.0 |
| 2 | 05.01.2013 | 0 | 25 | 2552 | 899.00 | -1.0 |
| 3 | 06.01.2013 | 0 | 25 | 2554 | 1709.05 | 1.0 |
| 4 | 15.01.2013 | 0 | 25 | 2555 | 1099.00 | 1.0 |
item_cat.head()
| item_category_name | item_category_id | |
|---|---|---|
| 0 | PC - Гарнитуры/Наушники | 0 |
| 1 | Аксессуары - PS2 | 1 |
| 2 | Аксессуары - PS3 | 2 |
| 3 | Аксессуары - PS4 | 3 |
| 4 | Аксессуары - PSP | 4 |
item.head()
| item_name | item_id | item_category_id | |
|---|---|---|---|
| 0 | ! ВО ВЛАСТИ НАВАЖДЕНИЯ (ПЛАСТ.) D | 0 | 40 |
| 1 | !ABBYY FineReader 12 Professional Edition Full... | 1 | 76 |
| 2 | ***В ЛУЧАХ СЛАВЫ (UNV) D | 2 | 40 |
| 3 | ***ГОЛУБАЯ ВОЛНА (Univ) D | 3 | 40 |
| 4 | ***КОРОБКА (СТЕКЛО) D | 4 | 40 |
shops.head(3)
| shop_name | shop_id | |
|---|---|---|
| 0 | !Якутск Орджоникидзе, 56 фран | 0 |
| 1 | !Якутск ТЦ "Центральный" фран | 1 |
| 2 | Адыгея ТЦ "Мега" | 2 |
test.head(3)
| ID | shop_id | item_id | |
|---|---|---|---|
| 0 | 0 | 5 | 5037 |
| 1 | 1 | 5 | 5320 |
| 2 | 2 | 5 | 5233 |
#formatting the date column correctly
sales.date=sales.date.apply(lambda x:datetime.datetime.strptime(x, '%d.%m.%Y'))
# check
print(sales.info())
print(sales.head(2))
<class 'pandas.core.frame.DataFrame'>
RangeIndex: 2935849 entries, 0 to 2935848
Data columns (total 6 columns):
date datetime64[ns]
date_block_num int64
shop_id int64
item_id int64
item_price float64
item_cnt_day float64
dtypes: datetime64[ns](1), float64(2), int64(3)
memory usage: 134.4 MB
None
date date_block_num shop_id item_id item_price item_cnt_day
0 2013-01-02 0 59 22154 999.0 1.0
1 2013-01-03 0 25 2552 899.0 1.0
# Aggregate to monthly level the required metrics
monthly_sales=sales.groupby(["date_block_num","shop_id","item_id"])[
"date","item_price","item_cnt_day"].agg({"date":["min",'max'],"item_price":"mean","item_cnt_day":"sum"})
## Lets break down the line of code here:
# aggregate by date-block(month),shop_id and item_id
# select the columns date,item_price and item_cnt(sales)
# Provide a dictionary which says what aggregation to perform on which column
# min and max on the date
# average of the item_price
# sum of the sales
monthly_sales.head(8)
| date | item_price | item_cnt_day | ||||
|---|---|---|---|---|---|---|
| min | max | mean | sum | |||
| date_block_num | shop_id | item_id | ||||
| 0 | 0 | 32 | 2013-01-03 | 2013-01-31 | 221.0 | 6.0 |
| 33 | 2013-01-03 | 2013-01-28 | 347.0 | 3.0 | ||
| 35 | 2013-01-31 | 2013-01-31 | 247.0 | 1.0 | ||
| 43 | 2013-01-31 | 2013-01-31 | 221.0 | 1.0 | ||
| 51 | 2013-01-13 | 2013-01-31 | 128.5 | 2.0 | ||
| 61 | 2013-01-10 | 2013-01-10 | 195.0 | 1.0 | ||
| 75 | 2013-01-17 | 2013-01-17 | 76.0 | 1.0 | ||
| 88 | 2013-01-16 | 2013-01-16 | 76.0 | 1.0 | ||
# number of items per cat
x=item.groupby(['item_category_id']).count()
x=x.sort_values(by='item_id',ascending=False)
x=x.iloc[0:10].reset_index()
x
| item_category_id | item_name | item_id | |
|---|---|---|---|
| 0 | 40 | 5035 | 5035 |
| 1 | 55 | 2365 | 2365 |
| 2 | 37 | 1780 | 1780 |
| 3 | 31 | 1125 | 1125 |
| 4 | 58 | 790 | 790 |
| 5 | 30 | 756 | 756 |
| 6 | 72 | 666 | 666 |
| 7 | 19 | 628 | 628 |
| 8 | 61 | 598 | 598 |
| 9 | 23 | 501 | 501 |
# plot
plt.figure(figsize=(8,4))
ax= sns.barplot(x.item_category_id, x.item_id, alpha=0.8)
plt.title("Items per Category")
plt.ylabel('# of items', fontsize=12)
plt.xlabel('Category', fontsize=12)
plt.show()
# 为什么前面 x 要赋值前十？因为10个数据值比较适合条形图可视化，而原84个值则数量太多，不宜画图。
# seaborn的条形图barplot,明显比matplotlib的条形图,好看很多。
# 还有个发现：seaborn的这个barplot，绘制的时候按照item_category_id升序排列了，估计是默认的。
# 回到分析正题：销量最高的是40号商品，达到了接近5000的量，其次是55号商品约2500销量和37号1900销量。销量第一名40号是遥遥领先。
# 大神原贴说，我们的目的是预测每个商品在每个商店的下个月销售额，属于时间序列预测类型；
# 那首先我们先找个简单的时间序列预测类型练练手，例如就预测下个月所有商品在所有商店的销售总额。
ts=sales.groupby(["date_block_num"])["item_cnt_day"].sum() # 按月周期，求总和销售额
ts.astype('float') # 转换字符格式成浮点数，方便计算。
plt.figure(figsize=(16,8)) # 画个图
plt.title('Total Sales of the company')
plt.xlabel('Time')
plt.ylabel('Sales')
plt.plot(ts); #如果此语句替换成 plt.show()，则图形结果是一张有坐标轴和标题的白板，没有任何折线或其他内容。突然想起plot是折线图。。。
# 另外，刚刚发现，上边最后一个分号“;”，原来作用类似等同于plt.show().——可减一行代码，不愧是简洁优美。
# 回到分析上来，随着时间推移，销售总额的两个峰值，答曰在第12个月和第24个月，刚好是年周期。莫非这是其中一个规律？应该就是的。
# 以12个月为周期，看看滚动变化。
# python时间序列分析之_用pandas中的rolling函数计算时间窗口数据——更详细用法可通过搜索引擎深入研究。
plt.figure(figsize=(16,6))
plt.plot(ts.rolling(window=12,center=False).mean(),label='Rolling Mean');
plt.plot(ts.rolling(window=12,center=False).std(),label='Rolling sd');
plt.legend(); # legend()函数表示显示标签，如果没有，则'Rolling Mean'和'Rolling sd'标签不显示。
# 明显看到，滚动月均销售总额在逐月下降，是经济形势逐月下降，还是公司战略和市场认可度发生了负面的变化呢？
# 作为数据分析师，发现这么要命的趋势变化，得进一步深挖原因，最好还能找到应对措施。
# 至于滚动月均标准差，隐隐约约展现了12个月一个周期，而且周期内缓慢上升；而周期外整体是下降趋势。
# 说明年初销售额下降幅度不大，到了年中甚至年尾，一年下来各月销售额下降越来越快，直到次年年初才又稳住下降趋势。
# 此处斗胆推测：
# （1）如果企业经营和选品没问题，那就是整体市场需求发生了变化，而且逐年下降。另外，此类商品在年初需求又比较大，在年中年尾需求较小。
# （2）如果整体市场需求稳定没有下降趋势，那就证明是企业经营或者是选品或者是营销宣传出了问题，需要深入排查解决。
# 下面深入分解：长期趋势Trend、季节性seasonality和随机残差residuals。
# 强行补充小知识：平稳性处理之“分解”
# 所谓分解就是将时序数据分离成不同的成分。statsmodels使用的X-11分解过程，它主要将时序数据分离成长期趋势、季节趋势和随机成分。
# 与其它统计软件一样，statsmodels也支持两类分解模型，加法模型和乘法模型，model的参数设置为"additive"（加法模型）和"multiplicative"（乘法模型）。
import statsmodels.api as sm # 导入统计建模模块
# multiplicative
res = sm.tsa.seasonal_decompose(ts.values,freq=12,model="multiplicative")
# 这里用到的.tsa.seasonal_decompose()函数，经尝试：参数ts.values时，横坐标是Time；参数ts时，横坐标是date_block_num。其他不变。
# freg这个参数容后研究，这里暂且猜测是周期12个月。
# plt.figure(figsize=(16,12))
fig = res.plot()
# fig.show() # 此句，可加可不加。
# 得到不同的分解成分，接下来可以使用时间序列模型对各个成分进行拟合。
res = sm.tsa.seasonal_decompose(ts.values,freq=12,model="additive")
#plt.figure(figsize=(16,12))
fig = res.plot()
# 这里加法模型的结果，和上面乘法模型的结果，几乎一模一样。
# Stationarity 平稳性:
# Stationarity refers to time-invariance of a series. (ie) Two points in a time series are related to each other by only how far apart they are, and not by the direction(forward/backward)
# When a time series is stationary, it can be easier to model. Statistical modeling methods assume or require the time series to be stationary.
# 当时间序列稳定的时候，它更容易被建模。统计学建模方式是假设或要求时间序列平稳的。
# There are multiple tests that can be used to check stationarity.
# 有各种各样的测试方法，可以用于平稳性检验。例如：
# ADF( Augmented Dicky Fuller Test)
# KPSS
# PP (Phillips-Perron test)
# 以上三种均属于时间序列平稳性检验方法中的“单位根检验法”。
# Let's just perform the ADF which is the most commonly used one.
# 下面我们只采用最常用的ADF检验方法。
# Stationarity tests 平稳性检验
# 写一个名为test_stationarity的统计性检验模块，以便将某些统计检验结果更加直观的展现出来。
def test_stationarity(timeseries):
#Perform Dickey-Fuller test:
print('Results of Dickey-Fuller Test:')
dftest = adfuller(timeseries, autolag='AIC')
dfoutput = pd.Series(dftest[0:4], index=['Test Statistic','p-value','#Lags Used','Number of Observations Used'])
# pd.series的index用于给特定行附上标签
for key,value in dftest[4].items():
dfoutput['Critical Value (%s)'%key] = value # 继续备注，方便查看
print (dfoutput)
test_stationarity(ts)
# adfuller函数说明：返回的结果是一个数组。
# 第一个是adt检验的结果，也就是书上说的t统计量的值。
# 第二个是t统计量的P值。
# 第三个是计算过程中用到的延迟阶数。
# 第四个是用于ADF回归和计算的观测值的个数。
# 第五个是配合第一个一起看的，是在99%，95%，90%置信区间下的临界的ADF检验的值。如果第一个值比第五个值小证明平稳，反正证明不平稳。
# 从以下结果看，第一个值 Test Statistic -2.395704 比 第五个值三个置信区间的临界ADF检验值都大，所以证明不平稳。
# 接下来就是要消除不平稳因素，获得新的满足平稳性的时间序列数据集。
Results of Dickey-Fuller Test:
Test Statistic -2.395704
p-value 0.142953
#Lags Used 0.000000
Number of Observations Used 33.000000
Critical Value (1%) -3.646135
Critical Value (5%) -2.954127
Critical Value (10%) -2.615968
dtype: float64
# to remove trend
from pandas import Series as Series
# create a differenced series 差分序列
def difference(dataset, interval=1):
diff = list()
for i in range(interval, len(dataset)):
value = dataset[i] - dataset[i - interval]
diff.append(value)
return Series(diff)
# invert differenced forecast
def inverse_difference(last_ob, value):
return value + last_ob
# 普及帖: 为什么差分时间序列数据?
# 差分是一种变换时间序列数据集的方法。
# 它可以用于消除序列对时间性的依赖性，即所谓的时间性依赖。这包含趋势和周期性的结构。
# 不同的方法可以帮助稳定时间序列的均值，消除时间序列的变化，从而消除(或减少)趋势和周期性。
# 通过从当前观察中减去先前观察值来实现差分。
# difference(t) = observation(t) - observation(t-1)
# 这样可以计算出序列差分。
# 延迟差分
# 将连续观察值之间的差值称为延迟-1差分。
# 可以调整延迟差分来适应特定的时间结构。
# 对于有周期性成分的时间序列，延迟可能是周期性的周期(宽度)。
# 差分序列
# 执行差分操作后，如非线性趋势的情况下，时间结构可能仍然存在。
# 因此，差分过程可以一直重复，直到所有时间依赖性被消除。
# 执行差分的次数称为差分序列。
# 上面那一步,是原贴中自定义差分函数def difference。
# diff（）函数
# 经搜索得知,pandas库有内置差分函数,可直接引用，自动计算差分数据集。这个diff（）函数是由Series和DataFrame对象提供。
# 就像前一节中手动定义的差分函数一样，它需要一个参数来指定间隔或延迟，在本例中称为周期（periods）。
# 优势是，代码行更少更简洁，且可以保留时间序列中时间和日期的信息。
# diff（）函数应用例子：
# from pandas import read_csv
# from pandas import datetime
# from matplotlib import pyplot
# def parser(x):
# return datetime.strptime('190'+x, '%Y-%m')
# series = read_csv('shampoo-sales.csv', header=0, parse_dates=[0], index_col=0, squeeze=True, date_parser=parser)
# diff = series.diff()
# pyplot.plot(diff)
# pyplot.show()
# 信息参考来源：http://www.atyun.com/4248.html
ts=sales.groupby(["date_block_num"])["item_cnt_day"].sum()
ts.astype('float')
plt.figure(figsize=(16,16))
plt.subplot(311)
plt.title('Original')
plt.xlabel('Time')
plt.ylabel('Sales')
plt.plot(ts) # 这图前面单个画过，这里用于与后面差分对比。
plt.subplot(312)
plt.title('After De-trend')
plt.xlabel('Time')
plt.ylabel('Sales')
new_ts=difference(ts)
plt.plot(new_ts)
plt.plot()
plt.subplot(313)
plt.title('After De-seasonalization')
plt.xlabel('Time')
plt.ylabel('Sales')
new_ts=difference(ts,12) # assuming the seasonality is 12 months long
plt.plot(new_ts)
plt.plot()
# now testing the stationarity again after de-seasonality
test_stationarity(new_ts)
Results of Dickey-Fuller Test:
Test Statistic -3.270101
p-value 0.016269
#Lags Used 0.000000
Number of Observations Used 21.000000
Critical Value (1%) -3.788386
Critical Value (5%) -3.013098
Critical Value (10%) -2.646397
dtype: float64
# 平稳性检验通过，达到了95%置信区间之内。下面开始深入预测工作。
# adding the dates to the Time-series as index
# 原数据集没有时间格式的月份，现用date_range函数以index的形式补上。
ts=sales.groupby(["date_block_num"])["item_cnt_day"].sum()
ts.index=pd.date_range(start = '2013-01-01',end='2015-10-01', freq = 'MS')
ts=ts.reset_index()
ts.head()
| index | item_cnt_day | |
|---|---|---|
| 0 | 2013-01-01 | 131479.0 |
| 1 | 2013-02-01 | 128090.0 |
| 2 | 2013-03-01 | 147142.0 |
| 3 | 2013-04-01 | 107190.0 |
| 4 | 2013-05-01 | 106970.0 |
# kaggle大神帖子用的Facebook开源工具【Prophet】，
# 由于本机一直提示网络断开安装不上（很有可能是翻墙不稳定），
# 所以转为参照书籍《Python数据科学》用另一种方法【平稳时间序列分析ARMA模型】继续。
# 在开始下一步之前，我们先通过《Python数据科学》第18章“时间序列建模”及其相关数据，对ARMA模型先做初步了解和熟悉。
# 一、平稳时间序列
# 实际中讨论的平稳时间序列是指对任意时间下，序列的均值、方差存在并为常数，且协方差函数与自相关系数只与时间间隔k有关。
# 只有平稳时间序列才可以进行统计分析，因为平稳性保证了时间序列数据都是出自同一分布，这样才可以计算均值、方差、延迟k期的协方差、延迟k期的相关系数。
# 一个独立同标准正态分布的随机序列就是平稳时间序列。
# 当然，若一个平稳时间序列的序列值之间没有相关性，那么就意味着这种数据前后没有规律，也就无法挖掘出有效的信息。这种序列称为纯随机序列。
# 在纯随机序列中，有一种序列称为白噪声序列，这种序列随机且各期的方差一致。
# 所以从这种意义上说，平稳时间序列分析在于充分挖掘时间序列之间的关系，当时间序列中的关系呗提取出来后，剩下的序列就应该是一个白噪声序列。
# 平稳时间序列模型主要有以下3种。
# （1）自回归模型（Auto Regression Model），简称 AR 模型。
# （2）移动平均模型（Moving Average Model），简称 MA 模型。
# （3）自回归移动平均模型（Auto Regression Moving Average Model），简称ARMA模型。
# 用于判断ARMA类型的自相关和偏自相关函数如下：
# （1）自相关函数（Autocorrelation Function，ACF）：描述时间序列任意两个时间间隔 k 的相关系数。
# （2）偏自相关函数（Partial Autocorrelation Function，PACF）：描述时间序列中在任意两个时间间隔 k 的时刻，去除1至k-1这个时间段中的其他数据的相关系数，这在统计学中称为偏相关系数。
# 二、ARMA模型
# 1、AR模型
# AR 模型，又称自回归模型，其认为时间序列当期观测值与前n期有线性关系，而与前n+1期无线性关系。
# 在AR(n)系统中，Xt具有 n 阶动态性。
# 以AR（1）模型为例，其中（1）表示滞后1期。
# AR（p）模型有以下重要性质。
# （1）某期观测值Xt的期望与系数序列α有关，方差有界。
# （2）自相关系数（ACF）拖尾，且值呈指数衰减（时间越近的往期观测对当期观测的影响越大）。
# （3）偏自相关函数（PACF） p 阶截尾。
# 其中ACF与PACF的性质可用于识别该平稳时间序列是适合滞后多少期的AR模型。
# 2、MA模型
# MA模型认为，如果一个系统在t时刻的相应Xt，与其以前时刻t-1,t-2,…的相应Xt-1,Xt-2，…无关，而与其以前时刻t-1,t-2,…，t-m 进入系统的扰动项εt-1，εt-2，…，εt-m存在一定的相关关系，那么这一类系统为MA（m）模型。
# 其中，εt是白噪声过程。
# MA（q）模型有以下重要性质。
# （1）t期系统扰动项εt的期望为常数，方差也为常数。
# （2）自相关系数（ACF） q 阶截尾。
# （3）偏自相关函数（PACF）拖尾。
# 其中ACF与PACF的性质可用于识别该平稳时间序列是适合滞后多少期的MA模型。
# 3、ARMA模型
# ARMA模型结合了AR模型与MA模型的共同特点，其认为序列是受到前期观测数据与系统扰动的共同影响。
# 具体来说，一个系统，如果它在时刻t的响应Xt，不仅与其以前时刻的自身值有关，还与其以前时刻进入系统的扰动项存在一定的依存关系，那么这个系统就是自回归移动平均模型。
# 对于平稳时间序列模型来说，AR模型、MA模型都属于ARMA（n,m）模型的特例。
# ARMA（p,q）模型的性质如下：
# （1）Xt的期望与系数序列α有关，方差有界。
# （2）自相关系数（ACF）拖尾。
# （3）偏自相关函数（PACF）拖尾。
# 4、ARMA模型的定阶与识别
# 之前介绍过AR模型、MA模型和ARMA模型的一些重要性质，其中自相关系数（ACF）与偏自相关系数（PACF）可以用于判断平稳时间序列数据适合哪一种模型和阶数。
# 平稳时间序列的识别方法：
# -AR模型：自相关系数拖尾，偏自相关系数截尾；
# -MA模型：自相关系数截尾，偏自相关函数拖尾；
# -ARMA模型：自相关函数和偏自相关函数均拖尾。
# 需要注意的是，一般样本长度大于50才能有一定的精确度。
# 5、使用AIC和BIC定阶与识别
# 使用ACF与PACF对ARMA模型进行定阶时，只能精确到MA模型与AR模型的阶数，而对于ARMA模型无法确定阶数，
# 而且由于估计误差的存在，实际中有时甚至很难判断AR模型与MA模型的截尾期数。
# 在实际操作中，可以使用以下方法识别ARMA模型：
# 通过AIC或者BIC准则识别，两个统计量都是越小越好。
# AIC或者BIC准则特别适用于ARMA模型，当然也适用于AR模型或者MA模型。
# 三、在Python中进行AR建模
# 在数据集ts_simu200.csv中有3份模拟的时间序列数据，分别为AR（1）模拟数据、MA（1）模拟数据和ARMA（1,1）模拟数据。
# 本节使用AR1_a序列演示在Python中建立 AR 模型。
# 1、探索平稳性
# 载入数据后，选择AR1_a列为原始数据，并将原始数据转换为时间序列数据。
ts_simu200 = pd.read_csv(r"D:\2018_BigData\Python\Python_book\18Timeseries\ts_simu200.csv")
ts_simu200.head(2)
| t | AR1_a | AR1_b | AR1_c | AR2 | ARIMA_110 | ARMA_11_a | ARMA_11_b | ARMA_22 | MA1_a | MA1_b | MA2 | EX1_a | EX1_b | EX1_c | EX1_d | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| 0 | 1 | -1.792035 | -0.230986 | -0.393871 | -0.870515 | -2.071533 | 1.027111 | 1.356889 | -5.104612 | 1.951711 | -0.519669 | -0.095012 | 0.640405 | 0.699464 | 0.389539 | 0.188441 |
| 1 | 2 | -0.743791 | -0.708602 | 0.286157 | -0.670041 | -2.189523 | -0.803688 | 2.027196 | -5.997497 | 0.407913 | 0.035999 | -0.130953 | -0.757885 | -1.696354 | -1.056063 | 0.405051 |
dates = pd.date_range(start="2017/01/01",periods=200)
ts_simu200.set_index(dates,inplace=True)
dta=ts_simu200["AR1_a"]
dta.head(3)
# 表头不见了，暂时不管。
2017-01-01 -1.792035
2017-01-02 -0.743791
2017-01-03 0.644999
Freq: D, Name: AR1_a, dtype: float64
plt.figure(figsize=(16,6))
plt.plot(dta);
# 如上图所示，初步看来时序图中的数据是平稳的，
# 因此可以进一步画出自相关系数（ACF）与偏自相关系数（PACF）的定阶图，以进一步验证平稳性并确定使用什么模型，阶数是多少。
# 2、定阶
# 在Python中可以使用acf与pacf函数画出ACF图与PACF图，如下所示：
# 自相关和偏自相关
fig = plt.figure(figsize=(16,16))
fig = sm.graphics.tsa.plot_acf(dta,lags=20) #lags表示滞后的阶数
fig = sm.graphics.tsa.plot_pacf(dta,lags=20)
plt.show()
# 代码运行效果如下图。从图中可以看出，ACF显示拖尾，而PACF显示1阶截尾，可以判断应使用AR(1)模型。
<Figure size 1152x1152 with 0 Axes>
# 3、AR 建模
# Python中的函数sm.tsa.ARMA，可以建立AR模型、MA模型、ARMA模型和带差分的ARIMA模型。
# 这里使用sm.tsa.ARMA函数，设定AR阶数为1，差分和MA的阶数都为0，即为AR（1）模型，对应的参数设置为（1,0）：
ar10 = sm.tsa.ARMA(dta,(1,0)).fit()
# 4、残差白噪声检验
# AR模型是否提取了原数据的足够信息的重要参考是AR模型的残差是否是白噪声序列，在Python中首先通过残差的自相关图和偏自相关图，代码如下：
# 检验下残差序列
resid = ar10.resid
fig = plt.figure(figsize=(12,8))
fig = sm.graphics.tsa.plot_acf(resid.values.squeeze(),lags=20)
fig = sm.graphics.tsa.plot_pacf(resid,lags=20)
plt.show()
# 残差的自相关图和偏自相关图如下图所示。可以看出，残差已经无信息可提取。
# 在ACF图中，残差滞后，各期均无显著的自相关性（ACF第0期代表与自身的相关性，其值恒为1）；
# 在PACF图中，各期也无显著的偏自相关性。可以判定，残差序列为白噪声序列。
<Figure size 864x576 with 0 Axes>
# 5、AR模型预测
# 利用下面代码预测未来20期的情况，并绘制曲线图，展示已有真实值、预测值及预测的置信区间。
predict_dta = ar10.forecast(steps=5)
import datetime
fig = ar10.plot_predict(
pd.to_datetime("2017-01-01")+datetime.timedelta(days=190),
pd.to_datetime("2017-01-01")+datetime.timedelta(days=220),
dynamic=False,plot_insample=True)
plt.show()
# 四、非平稳时间序列分析ARIMA模型
# 本节主要介绍针对非平稳时间序列使用的差分处理手段，将非平稳时间序列转换为平稳数据后再用ARMIA建模，还介绍了其在Python中的实现。
# 1、差分与ARIMA模型
# 实际上很多时间序列数据都是非平稳的，直接以平稳时间序列分析方法进行分析是不合适的。
# 可以通过差分等手段，将非平稳时间序列数据编程平稳时间序列，再采用ARIMA模型建模。
# （1）差分运算
# 差分运算是一种非常简便、有效的确定性信息提取的方法，而Cramer分解定理在理论上保证了适当阶数的差分一定可以充分提取确定性信息。
# 差分运算的实质是使用自回归的方式提取确定性信息，
# 1阶差分即用当期观测数据减去前1期的观测数据构成差分项；
# 2阶差分在1阶差分的基础上，对1阶差分的结果再进行差分。
# 以此类推，d阶差分是在d-1阶差分的基础上，对d-1阶差分的结果再进行差分。
# 适度的差分能够有效将非平稳时间序列转换为平稳时间序列。
# （2）ARIMA模型的建模步骤
# ARIMA模型适用于非平稳时间序列数据，其中的I表示差分的次数，适当的差分可使原序列称为平稳序列后，再进行ARIMA模型的建模。
# 其建模步骤与ARMA模型类似，分为以下5步。
# ①平稳：通过差分的手段，对非平稳时间序列数据进行平稳操作。
# ②定阶：确定宁ARIMA模型的阶数p、q。
# ③估计：估计未知参数。
# ④检验：检验残差是否是白噪声过程。
# ⑤预测：利用模型进行预测。
# 2、在Python中进行ARIMA建模
# （1）模拟数据ARIMA建模
# ARIMA建模与前面介绍的ARMA建模过程类似，分为5步。
# ①探索平稳性
# 载入数据后，选择ARIMA_110列为原始数据，将其转换为时间序列数据。
ts_simu200 = pd.read_csv(r"D:\2018_BigData\Python\Python_book\18Timeseries\ts_simu200.csv",index_col="t")
ts_simu200.head(2)
| AR1_a | AR1_b | AR1_c | AR2 | ARIMA_110 | ARMA_11_a | ARMA_11_b | ARMA_22 | MA1_a | MA1_b | MA2 | EX1_a | EX1_b | EX1_c | EX1_d | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| t | |||||||||||||||
| 1 | -1.792035 | -0.230986 | -0.393871 | -0.870515 | -2.071533 | 1.027111 | 1.356889 | -5.104612 | 1.951711 | -0.519669 | -0.095012 | 0.640405 | 0.699464 | 0.389539 | 0.188441 |
| 2 | -0.743791 | -0.708602 | 0.286157 | -0.670041 | -2.189523 | -0.803688 | 2.027196 | -5.997497 | 0.407913 | 0.035999 | -0.130953 | -0.757885 | -1.696354 | -1.056063 | 0.405051 |
dates = pd.date_range(start="2017/01",periods=200)
ts_simu200.set_index(dates,inplace=True)
dta=ts_simu200["ARIMA_110"]
dta.head(3)
2017-01-01 -2.071533
2017-01-02 -2.189523
2017-01-03 -2.305250
Freq: D, Name: ARIMA_110, dtype: float64
plt.figure(figsize=(16,6))
plt.plot(dta);
# 将数据序列可视化，效果如图。
# 从图中可以看出数据不同个时间段的均值差别较大，并不是一个平稳时间序列。
# 为进一步确定，可以调用Python中的包tseries中的函数adfuller进行平稳性检验。
# 该检验原假设为该时间序列是非平稳的，备择假设为时间序列是平稳的。
# 平稳性检验
result = adfuller(dta)
print("ADF Statistic:%f"% result[0])
print("p-value:%f"% result[1])
print()
print(result)
ADF Statistic:-1.897349
p-value:0.333309
(-1.8973490643758149, 0.33330922038075694, 1, 198, {'1%': -3.4638151713286316, '5%': -2.876250632135043, '10%': -2.574611347821651}, 504.6915217240421)
# 从检验结果上来看，p值为0.333，无法拒绝平稳假设，因此有足够理由认为该数据是非平稳的，所以需要对原数据进行差分。由于趋势是线性的，预先判断1阶差分就可以平稳。
# 使用函数diff进行差分，设定为1表示进行1阶差分，再画出差分后的时序图，如下所示：
# 差分序列的时序图
diff1 = dta.diff(1)
diff1.plot(figsize=(16,8));
# 1阶差分效果如下图。阶差分后，时序图中显示序列转化为了平稳时间序列。
# ②定阶
# 使用acf与pacf函数对差分数据画出ACF图与PACF图，如下。
# 差分序列的自相关ACF和偏自相关图PACF
ddta = diff1
ddta.dropna(inplace=True)
fig = plt.figure(figsize=(16,8))
fig = sm.graphics.tsa.plot_acf(ddta,lags=20)
fig = sm.graphics.tsa.plot_pacf(ddta,lags=20)
plt.show()
# 从图中可以看出，此数据的ACF拖尾，PACF为1阶截尾，可以判定应使用ARIMA（1,1,0）模型。
<Figure size 1152x576 with 0 Axes>
# ③ARIMA模型建模
# 使用Python中的函数sm.tsa.ARMA，设定AR阶数为1，差分为1,；MA模型的阶数都为0，即为ARIMA(1,1,0)模型，对应的参数设置为（1,1,0）。
arima110 = sm.tsa.ARIMA(dta,(1,1,0)).fit()
# ④残差白噪声检验
# 通过自相关和偏自相关图直观展示残差序列。
# 检验下残差序列
resid = arima110.resid
fig = plt.figure(figsize=(12,8))
fig = sm.graphics.tsa.plot_acf(resid.values.squeeze(),lags=20)
fig = sm.graphics.tsa.plot_pacf(resid,lags=20)
plt.show()
# 从图中可以看出，残差已经无信息可提取，在ACF图中，残差滞后，各期均无显著的自相关性；在PACF图中，残差没有滞后，各期也无显著的偏自相关性。
<Figure size 864x576 with 0 Axes>
# ⑤进行预测
# 利用下面代码预测未来20期的情况，并绘制曲线图，展示已有真实值、预测值及预测的置信区间。
import datetime
fig = arima110.plot_predict(pd.to_datetime('2017-01-01')+datetime.timedelta(days=1),
pd.to_datetime('2017-01-01')+datetime.timedelta(days=220),
dynamic=False, plot_insample=True)
plt.show()
# 如图：原始序列与向前20期的预测。
# （2）带季节性的ARIMA建模
# 这里使用一家拖拉机销售厂商的数据，其记录了2003年到2014年每月的拖拉机销售数据，数据带有明显的趋势性与季节性。如下：
sales_data = pd.read_csv(r"D:\2018_BigData\Python\Python_book\18Timeseries\tractor_sales.csv")
sales_data.head(15)
| Month-Year | Number of Tractor Sold | |
|---|---|---|
| 0 | 3-Jan | 141 |
| 1 | 3-Feb | 157 |
| 2 | 3-Mar | 185 |
| 3 | 3-Apr | 199 |
| 4 | 3-May | 203 |
| 5 | 3-Jun | 189 |
| 6 | 3-Jul | 207 |
| 7 | 3-Aug | 207 |
| 8 | 3-Sep | 171 |
| 9 | 3-Oct | 150 |
| 10 | 3-Nov | 138 |
| 11 | 3-Dec | 165 |
| 12 | 4-Jan | 145 |
| 13 | 4-Feb | 168 |
| 14 | 4-Mar | 197 |
sales_data.rename(columns={'Number of Tractor Sold':'Tractor-Sales'}, inplace=True)
sales_ts = sales_data['Tractor-Sales']
plt.figure(figsize=(12,5))
plt.plot(sales_ts)
plt.xlabel("Years")
plt.ylabel("Tractor Sales")
Text(0,0.5,'Tractor Sales')
# ①探索平稳性
# 需要注意的是，该数据的季节效应的波幅有明显增加，在这种情况下考虑对原始数据做对数。
plt.figure(figsize=(12,5))
plt.plot(np.log(sales_ts))
plt.xlabel("Years")
plt.ylabel("Log(Tractor Sales)")
Text(0,0.5,'Log(Tractor Sales)')
# 和原序列相比，对数变换后的时间序列的波幅基本一致了，且仍旧存在趋势性与季节性。因此对取对数的数据进行1阶差分，并查看其自相关和偏自相关函数。
sales_ts_log = np.log(sales_ts)
sales_ts_log.dropna(inplace=True)
sales_ts_log_diff = sales_ts_log.diff(periods=1)
sales_ts_log_diff.dropna(inplace=True)
fig,axes = plt.subplots(1,2,sharey=False,sharex=False)
fig.set_figwidth(12)
fig.set_figheight(4)
smt.graphics.plot_acf(sales_ts_log_diff,lags=30,ax=axes[0],alpha=0.5)
smt.graphics.plot_pacf(sales_ts_log_diff,lags=30,ax=axes[1],alpha=0.5)
plt.tight_layout()
# 通过以上分析，我们队数据的特征有了一定的认识，知道需要进行对数转换（因为趋势因素与季节因素都是累乘关系），并且需要进行原始数据前后1阶和季节前后1阶的差分。
# 在实际工作中，我们不再严格遵循使用自相关、偏自相关函数确定参数的方法，而是通过遍历所有可能的参数，以AIC统计量最小的模型作为最优模型。
# 以下是设置参数的语句：
# 设置自相关（AR）、差分（I）、移动平均（MA）的三个参数的取值范围
import itertools
p = d = q = range(0,2)
pdq=list(itertools.product(p,d,q))
# 设置季节效应的自相关（AR）、差分（I）、移动平均（MA）的三个参数的取值范围
seasonal_pdq = [(x[0],x[1],x[2],12) for x in list(itertools.product(p,d,q))]
# 下面语句中的核心函数是 sm.tsa.statespace.SARIMAX()，只有当数据有季节性时，才会使用该函数。
# 其中S代表该函数可以处理季节效应；X表示该函数可以加入外生变量，即对当前时间序列有具有预测作用的变量，类似线性回归中的X。
import sys
warnings.filterwarnings("ignore") # 忽略ARIMA模型无法估计出结果时的报警信息
best_aic = np.inf
best_pdq = None
best_seasonal_pdq = None
temp_model = None
for param in pdq:
for param_seasonal in seasonal_pdq:
try:
temp_model = sm.tsa.statespace.SARIMAX(sales_ts_log,
order = param,
seasonal_order = param_seasonal,
enforce_stationarity=True,
enforce_invertibility=True)
results = temp_model.fit()
if results.aic < best_aic:
best_aic = results.aic
best_pdq = param
best_seasonal_pdq = param_seasonal
except:
continue
print("Best SARIMAX{}x{}12 model - AIC:{}".format(best_pdq, best_seasonal_pdq, best_aic))
Best SARIMAX(0, 1, 1)x(1, 0, 1, 12)12 model - AIC:-495.2412861246029
# ②季节性ARIMA模型建模
# 根据上一步分析获得的参数，设置好SARIMAX()函数的参数，在此运行，并对残差进行检验。
best_model = sm.tsa.statespace.SARIMAX(sales_ts_log,
order=(0, 1, 1),
seasonal_order=(1, 0, 1, 12),
enforce_stationarity=True,
enforce_invertibility=True)
best_results = best_model.fit()
print(best_results.summary().tables[0])
print(best_results.summary().tables[1])
Statespace Model Results
==========================================================================================
Dep. Variable: Tractor-Sales No. Observations: 144
Model: SARIMAX(0, 1, 1)x(1, 0, 1, 12) Log Likelihood 251.621
Date: Thu, 09 May 2019 AIC -495.241
Time: 13:25:01 BIC -483.390
Sample: 0 HQIC -490.425
- 144
Covariance Type: opg
==========================================================================================
==============================================================================
coef std err z P>|z| [0.025 0.975]
------------------------------------------------------------------------------
ma.L1 -0.3574 0.069 -5.195 0.000 -0.492 -0.223
ar.S.L12 0.9933 0.006 175.826 0.000 0.982 1.004
ma.S.L12 -0.5532 0.097 -5.728 0.000 -0.742 -0.364
sigma2 0.0013 0.000 9.216 0.000 0.001 0.002
==============================================================================
best_results.plot_diagnostics(lags=30, figsize=(16,12))
plt.show()
# 下图是对残差进行的检验。可以确认服从正太分布，且不存在滞后效应。因此不能再提取任何信息。
# ③进行预测
# 利用下面代码预测未来36期的情况，并绘制曲线图，展示已有真实值、预测值及预测的置信区间。
import math
n_steps = 36
pred_uc_95 = best_results.get_forecast(steps=n_steps, alpha=0.05)
pred_pr_95 = pred_uc_95.predicted_mean
pred_ci_95 = pred_uc_95.conf_int()
idx = pd.date_range(sales_ts.index[-1], periods=n_steps, freq='MS')
fc_95 = pd.DataFrame(np.column_stack([np.power(math.e, pred_pr_95),
np.power(math.e, pred_ci_95)]), index=idx,
columns=['forecast', 'lower_ci_95', 'upper_ci_95'])
fc_95.head(2)
# 说明：
# 获取预测值的语句：best_results.get_forecast()，其中alpha表示置信区间，0.05表示置信区间为95%
# 将预测的均值提取出来：pred_uc_95.predicted_mean()
# 将预测值的置信区间提取出来：pred_uc_95.conf_inf()
# 由于之前对数据取了自然对数，现在这些预测数值需要取自然指数。
| forecast | lower_ci_95 | upper_ci_95 | |
|---|---|---|---|
| 1970-01-01 | 567.464517 | 528.213603 | 609.632119 |
| 1970-02-01 | 566.201497 | 519.958562 | 616.557084 |
sales_ts.head()
0 141
1 157
2 185
3 199
4 203
Name: Tractor-Sales, dtype: int64
# dates = pd.date_range(start='2003-01-01', freq='MS', periods=len(sales_data))
# sales_ts["date"]=dates
# # sales_ts.set_index(dates, inplace=True)
# sales_ts.head(2)
# 将预测的结果通过图形展现出来。
# axis = sales_ts.plot(label='Observed', figsize=(15, 6))
# fc_95['forecast'].plot(ax=axis, label='Forecast', alpha=0.7)
# axis.fill_between(fc_95.index, fc_95['lower_ci_95'], fc_95['upper_ci_95'],
# color='k', alpha=.25)
# axis.set_xlabel('Years')
# axis.set_ylabel('Tractor Sales')
# plt.legend(loc='best')
# plt.show()
# 报错，暂时忽略。先回到kaggle比赛。alueError: view limit minimum -36011.4 is less than 1 and is an invalid Matplotlib date value. This often happens if you pass a non-datetime value to an axis that has datetime units
# 【正题】回到kaggle比赛：商品销售预测
ts
| index | item_cnt_day | |
|---|---|---|
| 0 | 2013-01-01 | 131479.0 |
| 1 | 2013-02-01 | 128090.0 |
| 2 | 2013-03-01 | 147142.0 |
| 3 | 2013-04-01 | 107190.0 |
| 4 | 2013-05-01 | 106970.0 |
| 5 | 2013-06-01 | 125381.0 |
| 6 | 2013-07-01 | 116966.0 |
| 7 | 2013-08-01 | 125291.0 |
| 8 | 2013-09-01 | 133332.0 |
| 9 | 2013-10-01 | 127541.0 |
| 10 | 2013-11-01 | 130009.0 |
| 11 | 2013-12-01 | 183342.0 |
| 12 | 2014-01-01 | 116899.0 |
| 13 | 2014-02-01 | 109687.0 |
| 14 | 2014-03-01 | 115297.0 |
| 15 | 2014-04-01 | 96556.0 |
| 16 | 2014-05-01 | 97790.0 |
| 17 | 2014-06-01 | 97429.0 |
| 18 | 2014-07-01 | 91280.0 |
| 19 | 2014-08-01 | 102721.0 |
| 20 | 2014-09-01 | 99208.0 |
| 21 | 2014-10-01 | 107422.0 |
| 22 | 2014-11-01 | 117845.0 |
| 23 | 2014-12-01 | 168755.0 |
| 24 | 2015-01-01 | 110971.0 |
| 25 | 2015-02-01 | 84198.0 |
| 26 | 2015-03-01 | 82014.0 |
| 27 | 2015-04-01 | 77827.0 |
| 28 | 2015-05-01 | 72295.0 |
| 29 | 2015-06-01 | 64114.0 |
| 30 | 2015-07-01 | 63187.0 |
| 31 | 2015-08-01 | 66079.0 |
| 32 | 2015-09-01 | 72843.0 |
| 33 | 2015-10-01 | 71056.0 |
plt.figure(figsize=(12,5))
plt.xlabel("Time Series",size=16)
plt.ylabel("Total Monthly Sales",size=16)
# ax.set_xticklabels(ts.index) 无法显示横轴标为月份标签，暂时忽略。
plt.plot(ts.item_cnt_day);
# 如前文所述，有较为明显的12期一波峰，整体下降趋势。
# （1）探索平稳性
# 该数据的季节效应的波幅有轻微变化，无法判断其影响；在这种情况下考虑对原数据取对数观察。
plt.figure(figsize=(12,5))
plt.xlabel("Time Series",size=16)
plt.ylabel("Log(Total Monthly Sales)",size=16)
plt.plot(np.log(ts.item_cnt_day));
# 和原序列相比，对数变换后的时间序列波幅，与原序列几乎一致，当然也仍旧存在趋势性与季节性。所以原序列无需取对数即可进行下一步：1阶差分。
# 对原数据进行 1 阶差分，并查看其自相关和偏自相关函数。
ts_diff = ts.item_cnt_day.diff(periods=1)
ts_diff.dropna(inplace=True)
fig,axes = plt.subplots(1,2,sharey=False,sharex=False)
fig.set_figwidth(12)
fig.set_figheight(4)
smt.graphics.plot_acf(ts_diff,lags=30,ax=axes[0],alpha=0.5)
smt.graphics.plot_pacf(ts_diff,lags=30,ax=axes[1],alpha=0.5)
plt.tight_layout()
# 从图中可以看出，滞后几期内的ACF中只有第一期较显著，PACF则19期和30期较显著。
# 从发现经营异常的角度看，PACF中，19期和30期都是值得关注的时间段。
# 从自相关函数ACF看，以12期为一周期Q，1期（1个月）为P，可以看到：
# 前两周期Q满足P1月正数（增长），P2-3月负数（下降），4月增长，5-6月下降，7月增长，8-9月下降，10月增长，11-12月下降。
# 而第三周期，则1-4月增长，5-6月下降，7月增长。将原销售下降的2-3月，拉回成了销售增长。
# 因此，可以判断季节性的间隔周期为12期（即12个月），而且季节效应为AR（1）。
# 试试 2 阶差分效果
ts_diff = ts.item_cnt_day.diff(periods=2)
ts_diff.dropna(inplace=True)
fig,axes = plt.subplots(1,2,sharey=False,sharex=False)
fig.set_figwidth(12)
fig.set_figheight(4)
smt.graphics.plot_acf(ts_diff,lags=30,ax=axes[0],alpha=0.5)
smt.graphics.plot_pacf(ts_diff,lags=30,ax=axes[1],alpha=0.5)
plt.tight_layout()
# 从图中可以看出，滞后几期内的ACF中只有第一期较显著，PACF则25期较显著。
# 从发现经营异常的角度看，PACF中，18期和25期都是值得关注的时间段。
# 下面我们通过遍历所有可能的参数，以AIC统计量最小的模型作为最优模型。
# 设置自相关（AR）、差分（I）、移动平均（MA）的三个参数的取值范围
import itertools
p = d = q = range(0,2)
pdq=list(itertools.product(p,d,q))
# 设置季节效应的自相关（AR）、差分（I）、移动平均（MA）的三个参数的取值范围
seasonal_pdq = [(x[0],x[1],x[2],12) for x in list(itertools.product(p,d,q))]
# 下面语句中的核心函数是 sm.tsa.statespace.SARIMAX()，只有当数据有季节性时，才会使用该函数。
# 其中S代表该函数可以处理季节效应；X表示该函数可以加入外生变量，即对当前时间序列有具有预测作用的变量，类似线性回归中的X。
import sys
warnings.filterwarnings("ignore") # 忽略ARIMA模型无法估计出结果时的报警信息
best_aic = np.inf
best_pdq = None
best_seasonal_pdq = None
temp_model = None
for param in pdq:
for param_seasonal in seasonal_pdq:
try:
temp_model = sm.tsa.statespace.SARIMAX(ts.item_cnt_day,
order = param,
seasonal_order = param_seasonal,
enforce_stationarity=True,
enforce_invertibility=True)
results = temp_model.fit()
if results.aic < best_aic:
best_aic = results.aic
best_pdq = param
best_seasonal_pdq = param_seasonal
except:
continue
print("Best SARIMAX{}x{}12 model - AIC:{}".format(best_pdq, best_seasonal_pdq, best_aic))
# 通过执行该脚本，获得的最优模型为：SARIMAX(0, 1, 0)x(0, 1, 0, 12)，各参数的含义是：
# 第一个(0, 1, 0)代表该模型是ARIMA(0,1,0),即1阶差分I（1）模型；
# 第二个(0, 1, 0, 12)代表季节效应为ARIMA(0,1,0)，也是1阶差分I（1）模型；
# 最后的12代表季节效应为12期，这是我们之前自己设置的。
Best SARIMAX(0, 1, 0)x(0, 1, 0, 12)12 model - AIC:452.598370416518
# （2）季节性ARIMA模型建模
# 根据上一步分析获得的参数，设置好SARIMAX()函数的参数，在此运行，并对残差进行检验。
best_model = sm.tsa.statespace.SARIMAX(ts.item_cnt_day,
order=(0, 1, 0),
seasonal_order=(0, 1, 0, 12),
enforce_stationarity=True,
enforce_invertibility=True)
best_results = best_model.fit()
best_results.plot_diagnostics(lags=20, figsize=(16,12))
plt.show()
# 不服从正态分布，还有信息可以提取。怎么提取呢？先留个悬念。
# 随便看看模型参数和效果。
print(best_results.summary().tables[0])
print(best_results.summary().tables[1])
Statespace Model Results
==========================================================================================
Dep. Variable: item_cnt_day No. Observations: 34
Model: SARIMAX(0, 1, 0)x(0, 1, 0, 12) Log Likelihood -225.299
Date: Thu, 09 May 2019 AIC 452.598
Time: 13:25:05 BIC 453.643
Sample: 0 HQIC 452.825
- 34
Covariance Type: opg
==========================================================================================
==============================================================================
coef std err z P>|z| [0.025 0.975]
------------------------------------------------------------------------------
sigma2 1.147e+08 4.49e+07 2.555 0.011 2.67e+07 2.03e+08
==============================================================================
# （3）进行预测
# 利用下面代码预测未来12期的情况，并绘制曲线图，展示已有真实值、预测值及预测的置信区间。
n_steps = 12
pred_uc_95 = best_results.get_forecast(steps=n_steps, alpha=0.05)
pred_pr_95 = pred_uc_95.predicted_mean
pred_ci_95 = pred_uc_95.conf_int()
idx = pd.date_range(ts.index[-1], periods=n_steps, freq='MS')
fc_95 = pd.DataFrame(np.column_stack([pred_pr_95, pred_ci_95]), index=idx,
columns=['forecast', 'lower_ci_95', 'upper_ci_95'])
fc_95
# 说明：
# 获取预测值的语句：best_results.get_forecast()，其中alpha表示置信区间，0.05表示置信区间为95%
# 将预测的均值提取出来：pred_uc_95.predicted_mean()
# 将预测值的置信区间提取出来：pred_uc_95.conf_inf()
# 回到比赛原题：预测接下来一个月,俄罗斯某商超集团每件商品在各个商店的总销售额。
# 我们这里，得到了接下来不分商品不分商店的销售总额，81479。
| forecast | lower_ci_95 | upper_ci_95 | |
|---|---|---|---|
| 1970-01-01 | 81479.0 | 60491.682959 | 102466.317041 |
| 1970-02-01 | 132389.0 | 102708.451602 | 162069.548398 |
| 1970-03-01 | 74605.0 | 38253.900570 | 110956.099430 |
| 1970-04-01 | 47832.0 | 5857.365917 | 89806.634083 |
| 1970-05-01 | 45648.0 | -1281.067570 | 92577.067570 |
| 1970-06-01 | 41461.0 | -9947.217821 | 92869.217821 |
| 1970-07-01 | 35929.0 | -19598.221578 | 91456.221578 |
| 1970-08-01 | 27748.0 | -31613.096795 | 87109.096795 |
| 1970-09-01 | 26821.0 | -36140.951124 | 89782.951124 |
| 1970-10-01 | 29713.0 | -36654.723826 | 96080.723826 |
| 1970-11-01 | 36477.0 | -33130.055982 | 106084.055982 |
| 1970-12-01 | 34690.0 | -38012.198860 | 107392.198860 |
ts.index[-1]
33
# 将预测的结果通过图形展现出来。
axis = ts.item_cnt_day.plot(label='Observed', figsize=(15, 6))
# fc_95['forecast'].plot(ax=axis, label='Forecast', alpha=0.7)
# axis.fill_between(fc_95.index, fc_95['lower_ci_95'], fc_95['upper_ci_95'],
# color='k', alpha=.25)
axis.set_xlabel('Years')
axis.set_ylabel('Total Monthly Sales')
plt.title("Observed",size=20)
plt.figure(figsize=(5, 4))
fc=fc_95['forecast'].plot()
fc.set_xlabel('Years')
fc.set_ylabel('Total Monthly Sales')
plt.title("Forecast",size=20)
plt.legend(loc='best')
plt.show()
# 回到比赛原题：预测接下来一个月,俄罗斯某商超集团每件商品在各个商店的总销售额。
# 我们前面，得到了接下来不分商品不分商店的销售总额，81479。
# 下面，按每件商品每个商店重新分组，预测各商品各商店接下来一个月的总销售额。
sales.head()
| date | date_block_num | shop_id | item_id | item_price | item_cnt_day | amount | |
|---|---|---|---|---|---|---|---|
| 0 | 2013-01-02 | 0 | 59 | 22154 | 999.00 | 1.0 | 0 |
| 1 | 2013-01-03 | 0 | 25 | 2552 | 899.00 | 1.0 | 0 |
| 2 | 2013-01-05 | 0 | 25 | 2552 | 899.00 | -1.0 | 0 |
| 3 | 2013-01-06 | 0 | 25 | 2554 | 1709.05 | 1.0 | 0 |
| 4 | 2013-01-15 | 0 | 25 | 2555 | 1099.00 | 1.0 | 0 |
# 如上，查看sales数据内容，发现各条记录含单价和数量，但是没有销售额，所以接下来我们先补充一列“销售额”。
# 参考 data_4temp3["cv_balance"] = data_4temp3[["avg_balance","stdev_balance"]].apply(lambda x: x[1]/x[0],axis=1)
sales["amount"]=sales[["item_price","item_cnt_day"]].apply(lambda x:sales["item_price"]*sales["item_cnt_day"])
print(sales.info())
print(sales.head(3))
<class 'pandas.core.frame.DataFrame'>
RangeIndex: 2935849 entries, 0 to 2935848
Data columns (total 7 columns):
date datetime64[ns]
date_block_num int64
shop_id int64
item_id int64
item_price float64
item_cnt_day float64
amount float64
dtypes: datetime64[ns](1), float64(3), int64(3)
memory usage: 156.8 MB
None
date date_block_num shop_id item_id item_price item_cnt_day \
0 2013-01-02 0 59 22154 999.0 1.0
1 2013-01-03 0 25 2552 899.0 1.0
2 2013-01-05 0 25 2552 899.0 -1.0
amount
0 999.0
1 899.0
2 -899.0
# 前面已经按商品和商店分过组了，但是没有汇总销售金额，这次我们重新分组汇总。
month_sales=sales.groupby(["date_block_num","shop_id","item_id"])["item_cnt_day","amount"].agg({"item_cnt_day":"sum","amount":"sum"})
month_sales.reset_index()
month_sales.tail()
| item_cnt_day | amount | |||
|---|---|---|---|---|
| date_block_num | shop_id | item_id | ||
| 33 | 59 | 22087 | 6.0 | 714.0 |
| 22088 | 2.0 | 238.0 | ||
| 22091 | 1.0 | 179.0 | ||
| 22100 | 1.0 | 629.0 | ||
| 22102 | 1.0 | 1250.0 |
month_sales.columns
Index(['item_cnt_day', 'amount'], dtype='object')
month_sales.shape
# 总共16万条记录，我们的目标就是预测这16万条（对应单个商品单个商店34个月）中，约5000个单位商品单位商店接下来一个月的销售总额
(1609124, 2)
# (未完待续)
1万+
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
