前两篇主要介绍了一些线性模型AR/MA/ARMA/ARIMA。这篇主要介绍时间序列的季节模型。
参考文献1.《金融时间序列分析》 第2版 Ruey S.Tsay著 王辉、潘家柱 译2.TimeSeriesAnalysisTsa http://nipy.bic.berkeley.edu/nightly/statsmodels/doc/html/tsa.html3.时间序列预测全攻略（附带Python代码）
六、季节模型6.1 季节性时间序列有些时间序列呈现出一定的循环或周期性，这样的时间序列叫季节性时间序列。 下面我们以月入境旅游人数数据为例，看看序列是不是具有季节性
# 载入相关的库
from scipy import stats
import statsmodels.api as sm # 统计相关的库
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
temp = DataAPI.EcoDataProGet('2180000195','20100101','20151231')
temp = temp.set_index('periodDate')
temp = temp.sort()
temp['dataValue'].plot(figsize=(15,6),grid=True)
可以看出，序列看上去有着季节性，我们在看看每年的情况
compare={}
for index in temp['dataValue'].index:
year = index[:4]
month = index[5:7]
if year not in compare.keys():
compare.update({year:{month:temp['dataValue'][index]}})
else:
compare[year].update({month:temp['dataValue'][index]})
compare = pd.DataFrame(compare)
compare.plot(figsize=(20,6))
plt.legend(loc=0)
除了部分月份有一定偏差，序列整体的走势非常相似，那么可以说明有季节性么？我们先了解下季节指数
6.1.1 季节指数所谓季节指数就是用简单平均法计算的周期内各时期季节性影响的相对数。
首先计算周期内各期平均数
然后计算总平均数
再计算季节指数
季节指数反映了该季度与总平均值之间的一种比较稳定的关系：
如果比值大于1，说明该季度的值常常会高于总平均值；
如果比值小于1，说明该季度的值常常低于总平均值；
如果序列的季节指数都近似为1，就说明该序列没有明显的季节性。
我们来计算一下上面数据的季节指数：
ave_k = compare.mean(axis=1)
ave = ave_k.mean()
compare['S'] = ave_k/ave
compare
可以看出，虽然看上去具有季节性，但根据季节指数来看，季节指数除了2月份以外几乎全都接近1，序列的季节性很弱，不利于后面的建模。
好吧，看来只好我自己造一组数据了。
T = np.array([4,7,8,4,1,1,3])
data = T
for i in range(1,20):
x = T i
data = np.append(data,x)
plt.figure(figsize=(10,6))
plt.plot(data)
再来看看它的季节指数：
ave = data.mean()
term = len(T)
numOfTerm = len(data)/len(T)
S = []
for i in range(term):
s = 0
for j in range(numOfTerm):
s = data[i j*term]
S.append(s/numOfTerm)
S = S/ave
S
array([ 0.96296296, 1.18518519, 1.25925926, 0.96296296, 0.74074074, 0.74074074, 0.88888889])
可以看出，这次的季节指数比刚刚明显多了。
6.2 季节性差分在有些应用中，季节性是次要的，我们需要把它从数据中消除，这个过程叫季节调整，其中季节性差分化是一种常见的方法。
在之前我们介绍过差分（正规差分化），其形式为：
我们将它推广，如果一个序列具有周期性，且周期为s，则季节性差分化为：
我们以第一个上序列为例（旅游的），周期为12，我们进行季节性差分试试：
dataDiff = temp['dataValue'].diff(12)[12:]
dataDiff.plot(figsize=(15,6))
可以看出，季节性基本上被消除，但是序列不够平稳，我们不妨再做1阶差分看看
diffdiff = dataDiff.diff()[1:]
diffdiff.plot(figsize=(15,6)
t = sm.tsa.stattools.adfuller(diffdiff) # ADF检验
print "p-value: ",t[1]
p-value: 2.01727492597e-12
经过ADF检验，我们发现p-value小于显著性水平0.05，序列是平稳的。
6.3 多重季节性模型上一节中，经过了季节性差分和正规差分后，序列成为了平稳时间序列，则我们可以用MA模型对多重差分后的模型建模。则我们有模型：
其中，B为向后推移算子，代表前一时刻的值。B的s次方则为前s时刻的值。{at}为白噪声序列，s为序列的周期。上述模型为航空模型。
另外，序列经过了季节性差分和正规差分消除序列相关性（序列平稳），则可以认为是间隔s和间隔1的周期共同影响，因此该模型称为多重季节性MA模型
6.3.1 建模由于statsmodels中没有找到直接对序列建立季节模性的函数。我们使用参数估计的方式来建模。我们以一开始的旅游数据为例
如何进行参数估计？首先记：
根据白噪声性质，我们计算序列的方差与协方差如下：
接着，我们计算{ωt}的这三个参数:
data = diffdiff[:-5] # 最后5个用来预测
test = diffdiff[-5:]
var = data.var()
length = len(data)
acf,q,p = sm.tsa.acf(data,nlags =length,qstat=True) ## 计算自相关系数 及p-value
out = np.c_[range(1,length), acf[1:], q, p]
output=pd.DataFrame(out, columns=['lag', "AC", "Q", "P-value"])
output = output.set_index('lag')
print 'ACF1',output.ix[1]['AC']
print 'ACF12',output.ix[12]['AC']
ACF1 -0.550397467822
ACF12 -0.331684718252
我们知道，我们的s为12，根据公式：
两个自相关系数为：
cov_1 = output.ix[1]['AC']*var
cov_12 = output.ix[12]['AC']*var
然后我们根据之前推出的3个等式，代入Var,Cov1,Cov12，解3个方程的方程组。最终可归结为求两个二元一次方程：
# 求二元一次方程，返回绝对值大于1的解
def calc(a,b,c):
if b**2-4*a*c<0:
return np.nan
elif b**2-4*a*c==0:
return -b/(2*a)
else:
if abs((-b (b**2-4*a*c)**(1/2))/(2*a))<1 :
return (-b (b**2-4*a*c)**(1/2))/(2*a)
elif abs((-b-(b**2-4*a*c)**(1/2))/(2*a))<1:
return (-b-(b**2-4*a*c)**(1/2))/(2*a)
else:
return np.nan
print var/cov_1,var/cov_12
theta = calc(1,var/cov_1,1)
Theta = calc(1,var/cov_12,1)
print theta,Theta
-1.81686882383 -3.01491128463
nan nan
未能求出满足要求的θ和Θ的值。。。也就是序列难以用以上模型表达，即使用其他方法求出最优估计，应该满足的3条性质也无法满足
另一方面，在求季节指数时我们发现该序列季节性不强，这也是原因之一。
下面我们介绍其他模型。
6.4 综合分析模型（加法模型、乘法模型、混合模型）通常时间序列包括3个因素：
趋势因素T
季节性因素S
不规则因素I
有时候也会加入周期要素C,我们这里将它归为季节性因素
常用的综合分析模型有：
加分模型
乘法模型
混合模型
我们之前分析过，旅游数据的季节指数显示没有明显的季节性。因此我们这里以自己构造的数据为例。
T = np.array([4,7,8,4,1,1,3])
data = T
for i in range(1,20):
temp = T i
data = np.append(data,temp)
plt.figure(figsize=(10,6))
plt.plot(data)
我们选用加法模型来试试：
先分离季节性,我们可以看出周期为7，用移动平均来滤波
seasonRemove = pd.rolling_mean(data,7)
st = data - seasonRemove
plt.figure(figsize=(10,5))
plt.plot(seasonRemove,label = 'seasonRemove')
plt.plot(st,label = 'st')
plt.legend(loc=0)
对于去除季节性后的序列，我们用最小二乘来拟合：
from sklearn import linear_model
train = np.mat(range(7,len(seasonRemove))).T
clf = linear_model.LinearRegression()
clf.fit(train,seasonRemove[7:])
k = clf.coef_
print 'k=',k
x = range(7,len(seasonRemove))
y = k*x
b = (seasonRemove[7:] - y).mean()
print 'b=',b
y = y b
plt.figure(figsize=(10,6))
plt.plot(x,y,label='Tt')
plt.plot(seasonRemove,label='seasonal removed')
plt.legend(loc=0)
k = [ 0.14285714]
b = 3.14285714286
有
拟合效果很好， 接着根据残差计算
It=[]
for i in range(len(seasonRemove[7:])):
It.append(y[i]-seasonRemove[7 i])
plt.figure(figsize=(10,6))
plt.plot(It,label='It')
plt.legend()
看残差我们发现，虽然很小，量级为e-15，但不是平稳的，随着序列继续推进，误差渐渐变大。不过也可能是计算性能上的误差。
可以看出，加法模型很适合这个序列，可以有效分为趋势和季节两种特征。
6.4.2 分解利用sm.tsa.seasonal_decompose，可以直接将时间序列分解为趋势、季节和残差，并有加法模型和乘法模型两种模式可选。以上面的数据为例
result = sm.tsa.seasonal_decompose(np.array(data),model='additive',freq=1) #加法模型
plt.figure(figsize=(15,8))
plt.subplot(411)
plt.plot(result.observed, label='Original')
plt.legend(loc=0)
plt.subplot(412)
plt.plot(result.trend, label='Trend')
plt.legend(loc=0)
plt.subplot(413)
plt.plot(result.seasonal,label='Seasonality')
plt.legend(loc=0)
plt.subplot(414)
plt.plot(result.resid, label='Residuals')
plt.legend(loc=0)
本篇简单介绍了季节模型，由于书上太偏理论而且讲的内容太少，所以这篇写起来相当费劲儿，python可用的现成模型也少。都是网上查阅资料写的。
个人感觉季节模型最重要的是如何分离出季节特征，其中加法模型最常用，本篇也只举了简单的例子。想要深入还得继续研究。
,
