link之家
链接快照平台
- 输入网页链接，自动生成快照
- 标签化管理网页链接
相关文章推荐
|
纯真的熊猫 · Unable to deploy ...· 1 年前 · |
|
|
坚强的伤痕 · 黎明之后（烟若雨创作的系列小说作品）_百度百科· 1 年前 · |
|
|
慷慨的蚂蚁 · 台北君悅酒店/凱悅/GRAND HYATT ...· 1 年前 · |
|
|
眼睛小的野马 · | Hitomi.la· 1 年前 · |
|
|
乖乖的生姜 · 讲座回顾 | ...· 1 年前 · |
转载请注明出处。
什么是时间序列
时间序列简单的说就是各时间点上形成的数值序列，时间序列分析就是通过观察历史数据预测未来的值。在这里需要强调一点的是，时间序列分析并不是关于时间的回归，它主要是研究自身的变化规律的（这里不考虑含外生变量的时间序列）。
为什么用python
用 两个字总结“情怀”，爱屋及乌，个人比较喜欢python，就用python撸了。能做时间序列的软件很多，SAS、R、SPSS、Eviews甚至matlab等等，实际工作中应用得比较多的应该还是SAS和R，前者推荐王燕写的《应用时间序列分析》，后者推荐“ 基于R语言的时间序列建模完整教程 ”这篇博文（ 翻译版 ）。python作为科学计算的利器，当然也有相关分析的包:statsmodels中tsa模块，当然这个包和SAS、R是比不了，但是python有另一个神器：pandas！pandas在时间序列上的应用，能简化我们很多的工作。
python推荐直接装Anaconda，它集成了许多科学计算包，有一些包自己手动去装还是挺费劲的。statsmodels需要自己去安装，这里我推荐使用0.6的稳定版，0.7及其以上的版本能在github上找到，该版本在安装时会用C编译好，所以修改底层的一些代码将不会起作用。
时间序列分析
1.基本模型
自回归移动平均模型(ARMA(p，q))是时间序列中最为重要的模型之一，它主要由两部分组成： AR代表p阶自回归过程，MA代表q阶移动平均过程，其公式如下：
依据模型的形式、特性及自相关和偏自相关函数的特征，总结如下：
在时间序列中，ARIMA模型是在ARMA模型的基础上多了差分的操作。
2.pandas时间序列操作
大熊猫真的很可爱，这里简单介绍一下它在时间序列上的可爱之处。和许多时间序列分析一样，本文同样使用航空乘客数据（AirPassengers.csv）作为样例。
数据读取：
# -*- coding:utf-8 -*-
import numpy as np
import pandas as pd
from datetime import datetime
import matplotlib.pylab as plt
# 读取数据，pd.read_csv默认生成DataFrame对象，需将其转换成Series对象
df = pd.read_csv('AirPassengers.csv', encoding='utf-8', index_col='date')
df.index = pd.to_datetime(df.index) # 将字符串索引转换成时间索引
ts = df['x'] # 生成pd.Series对象
# 查看数据格式
ts.head()
ts.head().index
查看某日的值既可以使用字符串作为索引，又可以直接使用时间对象作为索引
ts['1949-01-01']
ts[datetime(1949,1,1)]
两者的返回值都是第一个序列值：112
如果要查看某一年的数据，pandas也能非常方便的实现
ts['1949']
切片操作：
ts['1949-1' : '1949-6']
注意时间索引的切片操作起点和尾部都是包含的，这点与数值索引有所不同
pandas还有很多方便的时间序列函数，在后面的实际应用中在进行说明。
3. 平稳性检验
我们知道序列平稳性是进行时间序列分析的前提条件，很多人都会有疑问，为什么要满足平稳性的要求呢？在大数定理和中心定理中要求样本同分布（这里同分布等价于时间序列中的平稳性），而我们的建模过程中有很多都是建立在大数定理和中心极限定理的前提条件下的，如果它不满足，得到的许多结论都是不可靠的。以虚假回归为例，当响应变量和输入变量都平稳时，我们用t统计量检验标准化系数的显著性。而当响应变量和输入变量不平稳时，其标准化系数不在满足t分布，这时再用t检验来进行显著性分析，导致拒绝原假设的概率增加，即容易犯第一类错误，从而得出错误的结论。
平稳时间序列有两种定义：严平稳和宽平稳
严平稳顾名思义，是一种条件非常苛刻的平稳性，它要求序列随着时间的推移，其统计性质保持不变。对于任意的τ，其联合概率密度函数满足：
严平稳的条件只是理论上的存在，现实中用得比较多的是宽平稳的条件。
宽平稳也叫弱平稳或者二阶平稳（均值和方差平稳），它应满足：
平稳性检验：观察法和单位根检验法
基于此，我写了一个名为test_stationarity的统计性检验模块，以便将某些统计检验结果更加直观的展现出来。
# -*- coding:utf-8 -*-
from statsmodels.tsa.stattools import adfuller
import pandas as pd
import matplotlib.pyplot as plt
import numpy as np
from statsmodels.graphics.tsaplots import plot_acf, plot_pacf
# 移动平均图
def draw_trend(timeSeries, size):
f = plt.figure(facecolor='white')
# 对size个数据进行移动平均
rol_mean = timeSeries.rolling(window=size).mean()
# 对size个数据进行加权移动平均
rol_weighted_mean = pd.ewma(timeSeries, span=size)
timeSeries.plot(color='blue', label='Original')
rolmean.plot(color='red', label='Rolling Mean')
rol_weighted_mean.plot(color='black', label='Weighted Rolling Mean')
plt.legend(loc='best')
plt.title('Rolling Mean')
plt.show()
def draw_ts(timeSeries):
f = plt.figure(facecolor='white')
timeSeries.plot(color='blue')
plt.show()
Unit Root Test
The null hypothesis of the Augmented Dickey-Fuller is that there is a unit
root, with the alternative that there is no unit root. That is to say the
bigger the p-value the more reason we assert that there is a unit root
def testStationarity(ts):
dftest = adfuller(ts)
# 对上述函数求得的值进行语义描述
dfoutput = pd.Series(dftest[0:4], index=['Test Statistic','p-value','#Lags Used','Number of Observations Used'])
for key,value in dftest[4].items():
dfoutput['Critical Value (%s)'%key] = value
return dfoutput
# 自相关和偏相关图，默认阶数为31阶
def draw_acf_pacf(ts, lags=31):
f = plt.figure(facecolor='white')
ax1 = f.add_subplot(211)
plot_acf(ts, lags=31, ax=ax1)
ax2 = f.add_subplot(212)
plot_pacf(ts, lags=31, ax=ax2)
plt.show()
观察法，通俗的说就是通过观察序列的趋势图与相关图是否随着时间的变化呈现出某种规律。所谓的规律就是时间序列经常提到的周期性因素，现实中遇到得比较多的是线性周期成分，这类周期成分可以采用差分或者移动平均来解决，而对于非线性周期成分的处理相对比较复杂，需要采用某些分解的方法。下图为航空数据的线性图，可以明显的看出它具有年周期成分和长期趋势成分。平稳序列的自相关系数会快速衰减，下面的自相关图并不能体现出该特征，所以我们有理由相信该序列是不平稳的。
单位根检验：ADF是一种常用的单位根检验方法，他的原假设为序列具有单位根，即非平稳，对于一个平稳的时序数据，就需要在给定的置信水平上显著，拒绝原假设。ADF只是单位根检验的方法之一，如果想采用其他检验方法，可以安装第三方包arch，里面提供了更加全面的单位根检验方法，个人还是比较钟情ADF检验。以下为检验结果，其p值大于0.99，说明并不能拒绝原假设。
3. 平稳性处理
由前面的分析可知，该序列是不平稳的，然而平稳性是时间序列分析的前提条件，故我们需要对不平稳的序列进行处理将其转换成平稳的序列。
a. 对数变换
对数变换主要是为了减小数据的振动幅度，使其线性规律更加明显（我是这么理解的时间序列模型大部分都是线性的，为了尽量降低非线性的因素，需要对其进行预处理，也许我理解的不对）。对数变换相当于增加了一个惩罚机制，数据越大其惩罚越大，数据越小惩罚越小。这里强调一下，变换的序列需要满足大于0，小于0的数据不存在对数变换。
ts_log = np.log(ts)
test_stationarity.draw_ts(ts_log)
b. 平滑法
根据平滑技术的不同，平滑法具体分为移动平均法和指数平均法。
移动平均即利用一定时间间隔内的平均值作为某一期的估计值，而指数平均则是用变权的方法来计算均值
test_stationarity.draw_trend(ts_log, 12)
从上图可以发现窗口为12的移动平均能较好的剔除年周期性因素，而指数平均法是对周期内的数据进行了加权，能在一定程度上减小年周期因素，但并不能完全剔除，如要完全剔除可以进一步进行差分操作。
c. 差分
时间序列最常用来剔除周期性因素的方法当属差分了，它主要是对等周期间隔的数据进行线性求减。前面我们说过，ARIMA模型相对ARMA模型，仅多了差分操作，ARIMA模型几乎是所有时间序列软件都支持的，差分的实现与还原都非常方便。而statsmodel中，对差分的支持不是很好，它不支持高阶和多阶差分，为什么不支持，这里引用作者的说法：
作者大概的意思是说预测方法中并没有解决高于2阶的差分，有没有感觉很牵强，不过没关系，我们有pandas。我们可以先用pandas将序列差分好，然后在对差分好的序列进行ARIMA拟合，只不过这样后面会多了一步人工还原的工作。
diff_12 = ts_log.diff(12)
diff_12.dropna(inplace=True)
diff_12_1 = diff_12.diff(1)
diff_12_1.dropna(inplace=True)
test_stationarity.testStationarity(diff_12_1)
从上面的统计检验结果可以看出，经过12阶差分和1阶差分后，该序列满足平稳性的要求了。
d. 分解
所谓分解就是将时序数据分离成不同的成分。statsmodels使用的X-11分解过程，它主要将时序数据分离成长期趋势、季节趋势和随机成分。与其它统计软件一样，statsmodels也支持两类分解模型，加法模型和乘法模型，这里我只实现加法，乘法只需将model的参数设置为"multiplicative"即可。
from statsmodels.tsa.seasonal import seasonal_decompose
decomposition = seasonal_decompose(ts_log, model="additive")
trend = decomposition.trend
seasonal = decomposition.seasonal
residual = decomposition.resid
得到不同的分解成分后，就可以使用时间序列模型对各个成分进行拟合，当然也可以选择其他预测方法。我曾经用过小波对时序数据进行过分解，然后分别采用时间序列拟合，效果还不错。由于我对小波的理解不是很好，只能简单的调用接口，如果有谁对小波、傅里叶、卡尔曼理解得比较透，可以将时序数据进行更加准确的分解，由于分解后的时序数据避免了他们在建模时的交叉影响，所以我相信它将有助于预测准确性的提高。
4. 模型识别
在前面的分析可知，该序列具有明显的年周期与长期成分。对于年周期成分我们使用窗口为12的移动平进行处理，对于长期趋势成分我们采用1阶差分来进行处理。
rol_mean = ts_log.rolling(window=12).mean()
rol_mean.dropna(inplace=True)
ts_diff_1 = rol_mean.diff(1)
ts_diff_1.dropna(inplace=True)
test_stationarity.testStationarity(ts_diff_1)
观察其统计量发现该序列在置信水平为95%的区间下并不显著，我们对其进行再次一阶差分。再次差分后的序列其自相关具有快速衰减的特点，t统计量在99%的置信水平下是显著的，这里我不再做详细说明。
ts_diff_2 = ts_diff_1.diff(1)
ts_diff_2.dropna(inplace=True)
数据平稳后，需要对模型定阶，即确定p、q的阶数。观察上图，发现自相关和偏相系数都存在拖尾的特点，并且他们都具有明显的一阶相关性，所以我们设定p=1, q=1。下面就可以使用ARMA模型进行数据拟合了。这里我不使用ARIMA(ts_diff_1, order=(1, 1, 1))进行拟合，是因为含有差分操作时，预测结果还原老出问题，至今还没弄明白。
from statsmodels.tsa.arima_model import ARMA
model = ARMA(ts_diff_2, order=(1, 1))
result_arma = model.fit( disp=-1, method='css')
5. 样本拟合
模型拟合完后，我们就可以对其进行预测了。由于ARMA拟合的是经过相关预处理后的数据，故其预测值需要通过相关逆变换进行还原。
predict_ts = result_arma.predict()
# 一阶差分还原
diff_shift_ts = ts_diff_1.shift(1)
diff_recover_1 = predict_ts.add(diff_shift_ts)
# 再次一阶差分还原
rol_shift_ts = rol_mean.shift(1)
diff_recover = diff_recover_1.add(rol_shift_ts)
# 移动平均还原
rol_sum = ts_log.rolling(window=11).sum()
rol_recover = diff_recover*12 - rol_sum.shift(1)
# 对数还原
log_recover = np.exp(rol_recover)
log_recover.dropna(inplace=True)
我们使用均方根误差（RMSE）来评估模型样本内拟合的好坏。利用该准则进行判别时，需要剔除“非预测”数据的影响。
ts = ts[log_recover.index] # 过滤没有预测的记录
plt.figure(facecolor='white')
log_recover.plot(color='blue', label='Predict')
ts.plot(color='red', label='Original')
plt.legend(loc='best')
plt.title('RMSE: %.4f'% np.sqrt(sum((log_recover-ts)**2)/ts.size))
plt.show()
观察上图的拟合效果，均方根误差为11.8828，感觉还过得去。
6. 完善ARIMA模型
前面提到statsmodels里面的ARIMA模块不支持高阶差分，我们的做法是将差分分离出来，但是这样会多了一步人工还原的操作。基于上述问题，我将差分过程进行了封装，使序列能按照指定的差分列表依次进行差分，并相应的构造了一个还原的方法，实现差分序列的自动还原。
# 差分操作
def diff_ts(ts, d):
global shift_ts_list
# 动态预测第二日的值时所需要的差分序列
global last_data_shift_list
shift_ts_list = []
last_data_shift_list = []
tmp_ts = ts
for i in d:
last_data_shift_list.append(tmp_ts[-i])
print last_data_shift_list
shift_ts = tmp_ts.shift(i)
shift_ts_list.append(shift_ts)
tmp_ts = tmp_ts - shift_ts
tmp_ts.dropna(inplace=True)
return tmp_ts
# 还原操作
def predict_diff_recover(predict_value, d):
if isinstance(predict_value, float):
tmp_data = predict_value
for i in range(len(d)):
tmp_data = tmp_data + last_data_shift_list[-i-1]
elif isinstance(predict_value, np.ndarray):
tmp_data = predict_value[0]
for i in range(len(d)):
tmp_data = tmp_data + last_data_shift_list[-i-1]
else:
tmp_data = predict_value
for i in range(len(d)):
try:
tmp_data = tmp_data.add(shift_ts_list[-i-1])
except:
raise ValueError('What you input is not pd.Series type!')
tmp_data.dropna(inplace=True)
return tmp_data
现在我们直接使用差分的方法进行数据处理，并以同样的过程进行数据预测与还原。
diffed_ts = diff_ts(ts_log, d=[12, 1])
model = arima_model(diffed_ts)
model.certain_model(1, 1)
predict_ts = model.properModel.predict()
diff_recover_ts = predict_diff_recover(predict_ts, d=[12, 1])
log_recover = np.exp(diff_recover_ts)
是不是发现这里的预测结果和上一篇的使用12阶移动平均的预测结果一模一样。这是因为12阶移动平均加上一阶差分与直接12阶差分是等价的关系，后者是前者数值的12倍，这个应该不难推导。
对于个数不多的时序数据，我们可以通过观察自相关图和偏相关图来进行模型识别，倘若我们要分析的时序数据量较多，例如要预测每只股票的走势，我们就不可能逐个去调参了。这时我们可以依据BIC准则识别模型的p, q值，通常认为BIC值越小的模型相对更优。这里我简单介绍一下BIC准则，它综合考虑了残差大小和自变量的个数，残差越小BIC值越小，自变量个数越多BIC值越大。个人觉得BIC准则就是对模型过拟合设定了一个标准（过拟合这东西应该以辩证的眼光看待）。
def proper_model(data_ts, maxLag):
init_bic = sys.maxint
init_p = 0
init_q = 0
init_properModel = None
for p in np.arange(maxLag):
for q in np.arange(maxLag):
model = ARMA(data_ts, order=(p, q))
try:
results_ARMA = model.fit(disp=-1, method='css')
except:
continue
bic = results_ARMA.bic
if bic < init_bic:
init_p = p
init_q = q
init_properModel = results_ARMA
init_bic = bic
return init_bic, init_p, init_q, init_properModel
相对最优参数识别结果：BIC: -1090.44209358 p: 0 q: 1 ， RMSE:11.8817198331。我们发现模型自动识别的参数要比我手动选取的参数更优。
7.滚动预测
所谓滚动预测是指通过添加最新的数据预测第二天的值。对于一个稳定的预测模型，不需要每天都去拟合，我们可以给他设定一个阀值，例如每周拟合一次，该期间只需通过添加最新的数据实现滚动预测即可。基于此我编写了一个名为arima_model的类，主要包含模型自动识别方法，滚动预测的功能，详细代码可以查看附录。数据的动态添加：
from dateutil.relativedelta import relativedelta
def _add_new_data(ts, dat, type='day'):
if type == 'day':
new_index = ts.index[-1] + relativedelta(days=1)
elif type == 'month':
new_index = ts.index[-1] + relativedelta(months=1)
ts[new_index] = dat
def add_today_data(model, ts, data, d, type='day'):
_add_new_data(ts, data, type) # 为原始序列添加数据
# 为滞后序列添加新值
d_ts = diff_ts(ts, d)
model.add_today_data(d_ts[-1], type)
def forecast_next_day_data(model, type='day'):
if model == None:
raise ValueError('No model fit before')
fc = model.forecast_next_day_value(type)
return predict_diff_recover(fc, [12, 1])
现在我们就可以使用滚动预测的方法向外预测了，取1957年之前的数据作为训练数据，其后的数据作为测试，并设定模型每第七天就会重新拟合一次。这里的diffed_ts对象会随着add_today_data方法自动添加数据，这是由于它与add_today_data方法中的d_ts指向的同一对象，该对象会动态的添加数据。
ts_train = ts_log[:'1956-12']
ts_test = ts_log['1957-1':]
diffed_ts = diff_ts(ts_train, [12, 1])
forecast_list = []
for i, dta in enumerate(ts_test):
if i%7 == 0:
model = arima_model(diffed_ts)
model.certain_model(1, 1)
forecast_data = forecast_next_day_data(model, type='month')
forecast_list.append(forecast_data)
add_today_data(model, ts_train, dta, [12, 1], type='month')
predict_ts = pd.Series(data=forecast_list, index=ts['1957-1':].index)
log_recover = np.exp(predict_ts)
original_ts = ts['1957-1':]
动态预测的均方根误差为：14.6479，与前面样本内拟合的均方根误差相差不大，说明模型并没有过拟合，并且整体预测效果都较好。
8. 模型序列化
在进行动态预测时，我们不希望将整个模型一直在内存中运行，而是希望有新的数据到来时才启动该模型。这时我们就应该把整个模型从内存导出到硬盘中，而序列化正好能满足该要求。序列化最常用的就是使用json模块了，但是它是时间对象支持得不是很好，有人对json模块进行了拓展以使得支持时间对象，这里我们不采用该方法，我们使用pickle模块，它和json的接口基本相同，有兴趣的可以去查看一下。我在实际应用中是采用的面向对象的编程，它的序列化主要是将类的属性序列化即可，而在面向过程的编程中，模型序列化需要将需要序列化的对象公有化，这样会使得对前面函数的参数改动较大，我不在详细阐述该过程。
与其它统计语言相比，python在统计分析这块还显得不那么“专业”。随着numpy、pandas、scipy、sklearn、gensim、statsmodels等包的推动，我相信也祝愿python在数据分析这块越来越好。与SAS和R相比，python的时间序列模块还不是很成熟，我这里仅起到抛砖引玉的作用，希望各位能人志士能贡献自己的力量，使其更加完善。实际应用中我全是面向过程来编写的，为了阐述方便，我用面向过程重新罗列了一遍，实在感觉很不方便。原本打算分三篇来写的，还有一部分实际应用的部分，不打算再写了，还请大家原谅。实际应用主要是具体问题具体分析，这当中第一步就是要查询问题，这步花的时间往往会比较多，然后再是解决问题。以我前面项目遇到的问题为例，当时遇到了以下几个典型的问题：1.周期长度不恒定的周期成分，例如每月的1号具有周期性，但每月1号与1号之间的时间间隔是不相等的；2.含有缺失值以及含有记录为0的情况无法进行对数变换；3.节假日的影响等等。
1 # -*-coding:utf-8-*-
2 import pandas as pd
3 import numpy as np
4 from statsmodels.tsa.arima_model import ARMA
5 import sys
6 from dateutil.relativedelta import relativedelta
7 from copy import deepcopy
8 import matplotlib.pyplot as plt
10 class arima_model:
12 def __init__(self, ts, maxLag=9):
13 self.data_ts = ts
14 self.resid_ts = None
15 self.predict_ts = None
16 self.maxLag = maxLag
17 self.p = maxLag
18 self.q = maxLag
19 self.properModel = None
20 self.bic = sys.maxint
22 # 计算最优ARIMA模型，将相关结果赋给相应属性
23 def get_proper_model(self):
24 self._proper_model()
25 self.predict_ts = deepcopy(self.properModel.predict())
26 self.resid_ts = deepcopy(self.properModel.resid)
28 # 对于给定范围内的p,q计算拟合得最好的arima模型，这里是对差分好的数据进行拟合，故差分恒为0
29 def _proper_model(self):
30 for p in np.arange(self.maxLag):
31 for q in np.arange(self.maxLag):
32 # print p,q,self.bic
33 model = ARMA(self.data_ts, order=(p, q))
34 try:
35 results_ARMA = model.fit(disp=-1, method='css')
36 except:
37 continue
38 bic = results_ARMA.bic
39 # print 'bic:',bic,'self.bic:',self.bic
40 if bic < self.bic:
41 self.p = p
42 self.q = q
43 self.properModel = results_ARMA
44 self.bic = bic
45 self.resid_ts = deepcopy(self.properModel.resid)
46 self.predict_ts = self.properModel.predict()
48 # 参数确定模型
49 def certain_model(self, p, q):
50 model = ARMA(self.data_ts, order=(p, q))
51 try:
52 self.properModel = model.fit( disp=-1, method='css')
53 self.p = p
54 self.q = q
55 self.bic = self.properModel.bic
56 self.predict_ts = self.properModel.predict()
57 self.resid_ts = deepcopy(self.properModel.resid)
58 except:
59 print 'You can not fit the model with this parameter p,q, ' \
60 'please use the get_proper_model method to get the best model'
62 # 预测第二日的值
63 def forecast_next_day_value(self, type='day'):
64 # 我修改了statsmodels包中arima_model的源代码，添加了constant属性，需要先运行forecast方法，为constant赋值
65 self.properModel.forecast()
66 if self.data_ts.index[-1] != self.resid_ts.index[-1]:
67 raise ValueError('''The index is different in data_ts and resid_ts, please add new data to data_ts.
68 If you just want to forecast the next day data without add the real next day data to data_ts,
69 please run the predict method which arima_model included itself''')
70 if not self.properModel:
71 raise ValueError('The arima model have not computed, please run the proper_model method before')
72 para = self.properModel.params
74 # print self.properModel.params
75 if self.p == 0: # It will get all the value series with setting self.data_ts[-self.p:] when p is zero
76 ma_value = self.resid_ts[-self.q:]
77 values = ma_value.reindex(index=ma_value.index[::-1])
78 elif self.q == 0:
79 ar_value = self.data_ts[-self.p:]
80 values = ar_value.reindex(index=ar_value.index[::-1])
81 else:
82 ar_value = self.data_ts[-self.p:]
83 ar_value = ar_value.reindex(index=ar_value.index[::-1])
84 ma_value = self.resid_ts[-self.q:]
85 ma_value = ma_value.reindex(index=ma_value.index[::-1])
86 values = ar_value.append(ma_value)
88 predict_value = np.dot(para[1:], values) + self.properModel.constant[0]
89 self._add_new_data(self.predict_ts, predict_value, type)
90 return predict_value
92 # 动态添加数据函数，针对索引是月份和日分别进行处理
93 def _add_new_data(self, ts, dat, type='day'):
94 if type == 'day':
95 new_index = ts.index[-1] + relativedelta(days=1)
96 elif type == 'month':
97 new_index = ts.index[-1] + relativedelta(months=1)
98 ts[new_index] = dat
100 def add_today_data(self, dat, type='day'):
101 self._add_new_data(self.data_ts, dat, type)
102 if self.data_ts.index[-1] != self.predict_ts.index[-1]:
103 raise ValueError('You must use the forecast_next_day_value method forecast the value of today before')
104 self._add_new_data(self.resid_ts, self.data_ts[-1] - self.predict_ts[-1], type)
106 if __name__ == '__main__':
107 df = pd.read_csv('AirPassengers.csv', encoding='utf-8', index_col='date')
108 df.index = pd.to_datetime(df.index)
109 ts = df['x']
111 # 数据预处理
112 ts_log = np.log(ts)
113 rol_mean = ts_log.rolling(window=12).mean()
114 rol_mean.dropna(inplace=True)
115 ts_diff_1 = rol_mean.diff(1)
116 ts_diff_1.dropna(inplace=True)
117 ts_diff_2 = ts_diff_1.diff(1)
118 ts_diff_2.dropna(inplace=True)
120 # 模型拟合
121 model = arima_model(ts_diff_2)
122 # 这里使用模型参数自动识别
123 model.get_proper_model()
124 print 'bic:', model.bic, 'p:', model.p, 'q:', model.q
125 print model.properModel.forecast()[0]
126 print model.forecast_next_day_value(type='month')
128 # 预测结果还原
129 predict_ts = model.properModel.predict()
130 diff_shift_ts = ts_diff_1.shift(1)
131 diff_recover_1 = predict_ts.add(diff_shift_ts)
132 rol_shift_ts = rol_mean.shift(1)
133 diff_recover = diff_recover_1.add(rol_shift_ts)
134 rol_sum = ts_log.rolling(window=11).sum()
135 rol_recover = diff_recover*12 - rol_sum.shift(1)
136 log_recover = np.exp(rol_recover)
137 log_recover.dropna(inplace=True)
139 # 预测结果作图
140 ts = ts[log_recover.index]
141 plt.figure(facecolor='white')
142 log_recover.plot(color='blue', label='Predict')
143 ts.plot(color='red', label='Original')
144 plt.legend(loc='best')
145 plt.title('RMSE: %.4f'% np.sqrt(sum((log_recover-ts)**2)/ts.size))
146 plt.show()
修改的arima_model代码
1 #
Note: The information criteria add 1 to the number of parameters
2 # whenever the model has an AR or MA term since, in principle,
3 # the variance could be treated as a free parameter and restricted
4 # This code does not allow this, but it adds consistency with other
5 # packages such as gretl and X12-ARIMA
7 from __future__ import absolute_import
8 from statsmodels.compat.python import string_types, range
9 # for 2to3 with extensions
11 from datetime import datetime
13 import numpy as np
14 from scipy import optimize
15 from scipy.stats import t, norm
16 from scipy.signal import lfilter
17 from numpy import dot, log, zeros, pi
18 from numpy.linalg import inv
20 from statsmodels.tools.decorators import (cache_readonly,
21 resettable_cache)
22 import statsmodels.tsa.base.tsa_model as tsbase
23 import statsmodels.base.wrapper as wrap
24 from statsmodels.regression.linear_model import yule_walker, GLS
25 from statsmodels.tsa.tsatools import (lagmat, add_trend,
26 _ar_transparams, _ar_invtransparams,
27 _ma_transparams, _ma_invtransparams,
28 unintegrate, unintegrate_levels)
29 from statsmodels.tsa.vector_ar import util
30 from statsmodels.tsa.ar_model import AR
31 from statsmodels.tsa.arima_process import arma2ma
32 from statsmodels.tools.numdiff import approx_hess_cs, approx_fprime_cs
33 from statsmodels.tsa.base.datetools import _index_date
34 from statsmodels.tsa.kalmanf import KalmanFilter
36 _armax_notes = """
38 Notes
39 -----
40 If exogenous variables are given, then the model that is fit is
42 .. math::
44 \\phi(L)(y_t - X_t\\beta) = \\theta(L)\epsilon_t
46 where :math:`\\phi` and :math:`\\theta` are polynomials in the lag
47 operator, :math:`L`. This is the regression model with ARMA errors,
48 or ARMAX model. This specification is used, whether or not the model
49 is fit using conditional sum of square or maximum-likelihood, using
50 the `method` argument in
51 :meth:`statsmodels.tsa.arima_model.%(Model)s.fit`. Therefore, for
52 now, `css` and `mle` refer to estimation methods only. This may
53 change for the case of the `css` model in future versions.
54 """
56 _arma_params = """\
57 endog : array-like
58 The endogenous variable.
59 order : iterable
60 The (p,q) order of the model for the number of AR parameters,
61 differences, and MA parameters to use.
62 exog : array-like, optional
63 An optional arry of exogenous variables. This should *not* include a
64 constant or trend. You can specify this in the `fit` method."""
66 _arma_model = "Autoregressive Moving Average ARMA(p,q) Model"
68 _arima_model = "Autoregressive Integrated Moving Average ARIMA(p,d,q) Model"
70 _arima_params = """\
71 endog : array-like
72 The endogenous variable.
73 order : iterable
74 The (p,d,q) order of the model for the number of AR parameters,
75 differences, and MA parameters to use.
76 exog : array-like, optional
77 An optional arry of exogenous variables. This should *not* include a
78 constant or trend. You can specify this in the `fit` method."""
80 _predict_notes = """
81 Notes
82 -----
83 Use the results predict method instead.
84 """
86 _results_notes = """
87 Notes
88 -----
89 It is recommended to use dates with the time-series models, as the
90 below will probably make clear. However, if ARIMA is used without
91 dates and/or `start` and `end` are given as indices, then these
92 indices are in terms of the *original*, undifferenced series. Ie.,
93 given some undifferenced observations::
95 1970Q1, 1
96 1970Q2, 1.5
97 1970Q3, 1.25
98 1970Q4, 2.25
99 1971Q1, 1.2
100 1971Q2, 4.1
102 1970Q1 is observation 0 in the original series. However, if we fit an
103 ARIMA(p,1,q) model then we lose this first observation through
104 differencing. Therefore, the first observation we can forecast (if
105 using exact MLE) is index 1. In the differenced series this is index
106 0, but we refer to it as 1 from the original series.
107 """
109 _predict = """
110 %(Model)s model in-sample and out-of-sample prediction
112 Parameters
113 ----------
114 %(params)s
115 start : int, str, or datetime
116 Zero-indexed observation number at which to start forecasting, ie.,
117 the first forecast is start. Can also be a date string to
118 parse or a datetime type.
119 end : int, str, or datetime
120 Zero-indexed observation number at which to end forecasting, ie.,
121 the first forecast is start. Can also be a date string to
122 parse or a datetime type. However, if the dates index does not
123 have a fixed frequency, end must be an integer index if you
124 want out of sample prediction.
125 exog : array-like, optional
126 If the model is an ARMAX and out-of-sample forecasting is
127 requested, exog must be given. Note that you'll need to pass
128 `k_ar` additional lags for any exogenous variables. E.g., if you
129 fit an ARMAX(2, q) model and want to predict 5 steps, you need 7
130 observations to do this.
131 dynamic : bool, optional
132 The `dynamic` keyword affects in-sample prediction. If dynamic
133 is False, then the in-sample lagged values are used for
134 prediction. If `dynamic` is True, then in-sample forecasts are
135 used in place of lagged dependent variables. The first forecasted
136 value is `start`.
137 %(extra_params)s
139 Returns
140 -------
141 %(returns)s
142 %(extra_section)s
143 """
145 _predict_returns = """predict : array
146 The predicted values.
148 """
150 _arma_predict = _predict % {"Model" : "ARMA",
151 "params" : """
152 params : array-like
153 The fitted parameters of the model.""",
154 "extra_params" : "",
155 "returns" : _predict_returns,
156 "extra_section" : _predict_notes}
158 _arma_results_predict = _predict % {"Model" : "ARMA", "params" : "",
159 "extra_params" : "",
160 "returns" : _predict_returns,
161 "extra_section" : _results_notes}
163 _arima_predict = _predict % {"Model" : "ARIMA",
164 "params" : """params : array-like
165 The fitted parameters of the model.""",
166 "extra_params" : """typ : str {'linear', 'levels'}
168 - 'linear' : Linear prediction in terms of the differenced
169 endogenous variables.
170 - 'levels' : Predict the levels of the original endogenous
171 variables.\n""", "returns" : _predict_returns,
172 "extra_section" : _predict_notes}
174 _arima_results_predict = _predict % {"Model" : "ARIMA",
175 "params" : "",
176 "extra_params" :
177 """typ : str {'linear', 'levels'}
179 - 'linear' : Linear prediction in terms of the differenced
180 endogenous variables.
181 - 'levels' : Predict the levels of the original endogenous
182 variables.\n""",
183 "returns" : _predict_returns,
184 "extra_section" : _results_notes}
186 _arima_plot_predict_example = """ Examples
187 --------
188 >>> import statsmodels.api as sm
189 >>> import matplotlib.pyplot as plt
190 >>> import pandas as pd
192 >>> dta = sm.datasets.sunspots.load_pandas().data[['SUNACTIVITY']]
193 >>> dta.index = pd.DatetimeIndex(start='1700', end='2009', freq='A')
194 >>> res = sm.tsa.ARMA(dta, (3, 0)).fit()
195 >>> fig, ax = plt.subplots()
196 >>> ax = dta.ix['1950':].plot(ax=ax)
197 >>> fig = res.plot_predict('1990', '2012', dynamic=True, ax=ax,
198 ... plot_insample=False)
199 >>> plt.show()
201 .. plot:: plots/arma_predict_plot.py
202 """
204 _plot_predict = ("""
205 Plot forecasts
206 """ + '\n'.join(_predict.split('\n')[2:])) % {
207 "params" : "",
208 "extra_params" : """alpha : float, optional
209 The confidence intervals for the forecasts are (1 - alpha)%
210 plot_insample : bool, optional
211 Whether to plot the in-sample series. Default is True.
212 ax : matplotlib.Axes, optional
213 Existing axes to plot with.""",
214 "returns" : """fig : matplotlib.Figure
215 The plotted Figure instance""",
216 "extra_section" : ('\n' + _arima_plot_predict_example +
217 '\n' + _results_notes)
218 }
220 _arima_plot_predict = ("""
221 Plot forecasts
222 """ + '\n'.join(_predict.split('\n')[2:])) % {
223 "params" : "",
224 "extra_params" : """alpha : float, optional
225 The confidence intervals for the forecasts are (1 - alpha)%
226 plot_insample : bool, optional
227 Whether to plot the in-sample series. Default is True.
228 ax : matplotlib.Axes, optional
229 Existing axes to plot with.
""",
230 "returns" : """fig : matplotlib.Figure
231 The plotted Figure instance""",
232 "extra_section" : ('\n' + _arima_plot_predict_example +
233 '\n' +
234 '\n'.join(_results_notes.split('\n')[:3]) +
235 ("""
236 This is hard-coded to only allow plotting of the forecasts in levels.
237 """) +
238 '\n'.join(_results_notes.split('\n')[3:]))
239 }
242 def cumsum_n(x, n):
243 if n:
244 n -= 1
245 x = np.cumsum(x)
246 return cumsum_n(x, n)
247 else:
248 return x
251 def _check_arima_start(start, k_ar, k_diff, method, dynamic):
252 if start < 0:
253 raise ValueError("The start index %d of the original series "
254 "has been differenced away" % start)
255 elif (dynamic or 'mle' not in method) and start < k_ar:
256 raise ValueError("Start must be >= k_ar for conditional MLE "
257 "or dynamic forecast. Got %d" % start)
260 def _get_predict_out_of_sample(endog, p, q, k_trend, k_exog, start, errors,
261 trendparam, exparams, arparams, maparams, steps,
262 method, exog=None):
263 """
264 Returns endog, resid, mu of appropriate length for out of sample
265 prediction.
266 """
267 if q:
268 resid = np.zeros(q)
269 if start and 'mle' in method or (start == p and not start == 0):
270 resid[:q] = errors[start-q:start]
271 elif start:
272 resid[:q] = errors[start-q-p:start-p]
273 else:
274 resid[:q] = errors[-q:]
275 else:
276 resid = None
278 y = endog
279 if k_trend == 1:
280 # use expectation not constant
281 if k_exog > 0:
282 #TODO: technically should only hold for MLE not
283 # conditional model. See #274.
284 # ensure 2-d for conformability
285 if np.ndim(exog) == 1 and k_exog == 1:
286 # have a 1d series of observations -> 2d
287 exog = exog[:, None]
288 elif np.ndim(exog) == 1:
289 # should have a 1d row of exog -> 2d
290 if len(exog) != k_exog:
291 raise ValueError("1d exog given and len(exog) != k_exog")
292 exog = exog[None, :]
293 X = lagmat(np.dot(exog, exparams), p, original='in', trim='both')
294 mu = trendparam * (1 - arparams.sum())
295 # arparams were reversed in unpack for ease later
296 mu = mu + (np.r_[1, -arparams[::-1]] * X).sum(1)[:, None]
297 else:
298 mu = trendparam * (1 - arparams.sum())
299 mu = np.array([mu]*steps)
300 elif k_exog > 0:
301 X = np.dot(exog, exparams)
302 #NOTE: you shouldn't have to give in-sample exog!
303 X = lagmat(X, p, original='in', trim='both')
304 mu = (np.r_[1, -arparams[::-1]] * X).sum(1)[:, None]
305 else:
306 mu = np.zeros(steps)
308 endog = np.zeros(p + steps - 1)
310 if p and start:
311 endog[:p] = y[start-p:start]
312 elif p:
313 endog[:p] = y[-p:]
315 return endog, resid, mu
318 def _arma_predict_out_of_sample(params, steps, errors, p, q, k_trend, k_exog,
319 endog, exog=None, start=0, method='mle'):
320 (trendparam, exparams,
321 arparams, maparams) = _unpack_params(params, (p, q), k_trend,
322 k_exog, reverse=True)
323 # print 'params:',params
324 # print 'arparams:',arparams,'maparams:',maparams
325 endog, resid, mu = _get_predict_out_of_sample(endog, p, q, k_trend, k_exog,
326 start, errors, trendparam,
327 exparams, arparams,
328 maparams, steps, method,
329 exog)
330 # print 'mu[-1]:',mu[-1], 'mu[0]:',mu[0]
331 forecast = np.zeros(steps)
332 if steps == 1:
333 if q:
334 return mu[0] + np.dot(arparams, endog[:p]) + np.dot(maparams,
335 resid[:q]), mu[0]
336 else
:
337 return mu[0] + np.dot(arparams, endog[:p]), mu[0]
339 if q:
340 i = 0 # if q == 1
341 else:
342 i = -1
344 for i in range(min(q, steps - 1)):
345 fcast = (mu[i] + np.dot(arparams, endog[i:i + p]) +
346 np.dot(maparams[:q - i], resid[i:i + q]))
347 forecast[i] = fcast
348 endog[i+p] = fcast
350 for i in range(i + 1, steps - 1):
351 fcast = mu[i] + np.dot(arparams, endog[i:i+p])
352 forecast[i] = fcast
353 endog[i+p] = fcast
355 #need to do one more without updating endog
356 forecast[-1] = mu[-1] + np.dot(arparams, endog[steps - 1:])
357 return forecast, mu[-1] #Modified by me, the former is return forecast
360 def _arma_predict_in_sample(start, end, endog, resid, k_ar, method):
361 """
362 Pre- and in-sample fitting for ARMA.
363 """
364 if 'mle' in method:
365 fittedvalues = endog - resid # get them all then trim
366 else:
367 fittedvalues = endog[k_ar:] - resid
369 fv_start = start
370 if 'mle' not in method:
371 fv_start -= k_ar # start is in terms of endog index
372 fv_end = min(len(fittedvalues), end + 1)
373 return fittedvalues[fv_start:fv_end]
376 def _validate(start, k_ar, k_diff, dates, method):
377 if isinstance(start, (string_types, datetime)):
378 start = _index_date(start, dates)
379 start -= k_diff
380 if 'mle' not in method and start < k_ar - k_diff:
381 raise ValueError("Start must be >= k_ar for conditional "
382 "MLE or dynamic forecast. Got %s" % start)
384 return start
387 def _unpack_params(params, order, k_trend, k_exog, reverse=False):
388 p, q = order
389 k = k_trend + k_exog
390 maparams = params[k+p:]
391 arparams = params[k:k+p]
392 trend = params[:k_trend]
393 exparams = params[k_trend:k]
394 if reverse:
395 return trend, exparams, arparams[::-1], maparams[::-1]
396 return trend, exparams, arparams, maparams
399 def _unpack_order(order):
400 k_ar, k_ma, k = order
401 k_lags = max(k_ar, k_ma+1)
402 return k_ar, k_ma, order, k_lags
405 def _make_arma_names(data, k_trend, order, exog_names):
406 k_ar, k_ma = order
407 exog_names = exog_names or []
408 ar_lag_names = util.make_lag_names([data.ynames], k_ar, 0)
409 ar_lag_names = [''.join(('ar.', i)) for i in ar_lag_names]
410 ma_lag_names = util.make_lag_names([data.ynames], k_ma, 0)
411 ma_lag_names = [''.join(('ma.', i)) for i in ma_lag_names]
412 trend_name = util.make_lag_names('', 0, k_trend)
413 exog_names = trend_name + exog_names + ar_lag_names + ma_lag_names
414 return exog_names
417 def _make_arma_exog(endog, exog, trend):
418 k_trend = 1 # overwritten if no constant
419 if exog is None and trend == 'c': # constant only
420 exog = np.ones((len(endog), 1))
421 elif exog is not None and trend == 'c': # constant plus exogenous
422 exog = add_trend(exog, trend='c', prepend=True)
423 elif exog is not None and trend == 'nc':
424 # make sure it's not holding constant from last run
425 if exog.var() == 0:
426 exog = None
427 k_trend = 0
428 if trend == 'nc':
429 k_trend = 0
430 return k_trend, exog
433 def _check_estimable(nobs, n_params):
434 if nobs <= n_params:
435 raise ValueError("Insufficient degrees of freedom to estimate")
438 class ARMA(tsbase.TimeSeriesModel):
440 __doc__ = tsbase._tsa_doc % {"model" : _arma_model,
441 "params" : _arma_params, "extra_params" : "",
442 "extra_sections" : _armax_notes %
443 {"Model" : "ARMA"}}
445 def __init__(self, endog, order, exog=None, dates=None, freq=None,
446 missing='none'):
447 super(ARMA, self).__init__(endog, exog, dates, freq, missing=missing)
448 exog = self.data.exog # get it after it's gone through processing
449 _check_estimable(len(self.endog), sum(order))
450 self.k_ar = k_ar = order[0]
451
self.k_ma = k_ma = order[1]
452 self.k_lags = max(k_ar, k_ma+1)
453 self.constant = 0 #Added by me
454 if exog is not None:
455 if exog.ndim == 1:
456 exog = exog[:, None]
457 k_exog = exog.shape[1] # number of exog. variables excl. const
458 else:
459 k_exog = 0
460 self.k_exog = k_exog
462 def _fit_start_params_hr(self, order):
463 """
464 Get starting parameters for fit.
466 Parameters
467 ----------
468 order : iterable
469 (p,q,k) - AR lags, MA lags, and number of exogenous variables
470 including the constant.
472 Returns
473 -------
474 start_params : array
475 A first guess at the starting parameters.
477 Notes
478 -----
479 If necessary, fits an AR process with the laglength selected according
480 to best BIC. Obtain the residuals. Then fit an ARMA(p,q) model via
481 OLS using these residuals for a first approximation. Uses a separate
482 OLS regression to find the coefficients of exogenous variables.
484 References
485 ----------
486 Hannan, E.J. and Rissanen, J. 1982. "Recursive estimation of mixed
487 autoregressive-moving average order." `Biometrika`. 69.1.
488 """
489 p, q, k = order
490 start_params = zeros((p+q+k))
491 endog = self.endog.copy() # copy because overwritten
492 exog = self.exog
493 if k != 0:
494 ols_params = GLS(endog, exog).fit().params
495 start_params[:k] = ols_params
496 endog -= np.dot(exog, ols_params).squeeze()
497 if q != 0:
498 if p != 0:
499 # make sure we don't run into small data problems in AR fit
500 nobs = len(endog)
501 maxlag = int(round(12*(nobs/100.)**(1/4.)))
502 if maxlag >= nobs:
503 maxlag = nobs - 1
504 armod = AR(endog).fit(ic='bic', trend='nc', maxlag=maxlag)
505 arcoefs_tmp = armod.params
506 p_tmp = armod.k_ar
507 # it's possible in small samples that optimal lag-order
508 # doesn't leave enough obs. No consistent way to fix.
509 if p_tmp + q >= len(endog):
510 raise ValueError("Proper starting parameters cannot"
511 " be found for this order with this "
512 "number of observations. Use the "
513 "start_params argument.")
514 resid = endog[p_tmp:] - np.dot(lagmat(endog, p_tmp,
515 trim='both'),
516 arcoefs_tmp)
517 if p < p_tmp + q:
518 endog_start = p_tmp + q - p
519 resid_start = 0
520 else:
521 endog_start = 0
522 resid_start = p - p_tmp - q
523 lag_endog = lagmat(endog, p, 'both')[endog_start:]
524 lag_resid = lagmat(resid, q, 'both')[resid_start:]
525 # stack ar lags and resids
526 X = np.column_stack((lag_endog, lag_resid))
527 coefs = GLS(endog[max(p_tmp + q, p):], X).fit().params
528 start_params[k:k+p+q] = coefs
529 else:
530 start_params[k+p:k+p+q] = yule_walker(endog, order=q)[0]
531 if q == 0 and p != 0:
532 arcoefs = yule_walker(endog, order=p)[0]
533 start_params[k:k+p] = arcoefs
535 # check AR coefficients
536 if p and not np.all(np.abs(np.roots(np.r_[1, -start_params[k:k + p]]
537 )) < 1):
538 raise ValueError("The computed initial AR coefficients are not "
539 "stationary\nYou should induce stationarity, "
540 "choose a different model order, or you can\n"
541 "pass your own start_params.")
542 # check MA coefficients
543 elif q and not np.all(np.abs(np.roots(np.r_[1, start_params[k + p:]]
544 )) < 1):
545 return np.zeros(len(start_params)) #modified by me
546 raise ValueError("The computed initial MA coefficients are not "
547 "invertible\nYou should induce invertibility, "
548 "choose a different model order, or you can\n"
549 "pass your own start_params.")
551 # check MA coefficients
552 # print start_params
553 return start_params
555 def _fit_start_params(self, order, method):
556 if method != 'css-mle': # use Hannan-Rissanen to get start params
557 start_params = self._fit_start_params_hr(order)
558 else: # use CSS to get start params
559
func = lambda params: -self.loglike_css(params)
560 #start_params = [.1]*(k_ar+k_ma+k_exog) # different one for k?
561 start_params = self._fit_start_params_hr(order)
562 if self.transparams:
563 start_params = self._invtransparams(start_params)
564 bounds = [(None,)*2]*sum(order)
565 mlefit = optimize.fmin_l_bfgs_b(func, start_params,
566 approx_grad=True, m=12,
567 pgtol=1e-7, factr=1e3,
568 bounds=bounds, iprint=-1)
569 start_params = self._transparams(mlefit[0])
570 return start_params
572 def score(self, params):
573 """
574 Compute the score function at params.
576 Notes
577 -----
578 This is a numerical approximation.
579 """
580 return approx_fprime_cs(params, self.loglike, args=(False,))
582 def hessian(self, params):
583 """
584 Compute the Hessian at params,
586 Notes
587 -----
588 This is a numerical approximation.
589 """
590 return approx_hess_cs(params, self.loglike, args=(False,))
592 def _transparams(self, params):
593 """
594 Transforms params to induce stationarity/invertability.
596 Reference
597 ---------
598 Jones(1980)
599 """
600 k_ar, k_ma = self.k_ar, self.k_ma
601 k = self.k_exog + self.k_trend
602 newparams = np.zeros_like(params)
604 # just copy exogenous parameters
605 if k != 0:
606 newparams[:k] = params[:k]
608 # AR Coeffs
609 if k_ar != 0:
610 newparams[k:k+k_ar] = _ar_transparams(params[k:k+k_ar].copy())
612 # MA Coeffs
613 if k_ma != 0:
614 newparams[k+k_ar:] = _ma_transparams(params[k+k_ar:].copy())
615 return newparams
617 def _invtransparams(self, start_params):
618 """
619 Inverse of the Jones reparameterization
620 """
621 k_ar, k_ma = self.k_ar, self.k_ma
622 k = self.k_exog + self.k_trend
623 newparams = start_params.copy()
624 arcoefs = newparams[k:k+k_ar]
625 macoefs = newparams[k+k_ar:]
626 # AR coeffs
627 if k_ar != 0:
628 newparams[k:k+k_ar] = _ar_invtransparams(arcoefs)
630 # MA coeffs
631 if k_ma != 0:
632 newparams[k+k_ar:k+k_ar+k_ma] = _ma_invtransparams(macoefs)
633 return newparams
635 def _get_predict_start(self, start, dynamic):
636 # do some defaults
637 method = getattr(self, 'method', 'mle')
638 k_ar = getattr(self, 'k_ar', 0)
639 k_diff = getattr(self, 'k_diff', 0)
640 if start is None:
641 if 'mle' in method and not dynamic:
642 start = 0
643 else:
644 start = k_ar
645 self._set_predict_start_date(start) # else it's done in super
646 elif isinstance(start, int):
647 start = super(ARMA, self)._get_predict_start(start)
648 else: # should be on a date
649 #elif 'mle' not in method or dynamic: # should be on a date
650 start = _validate(start, k_ar, k_diff, self.data.dates,
651 method)
652 start = super(ARMA, self)._get_predict_start(start)
653 _check_arima_start(start, k_ar, k_diff, method, dynamic)
654 return start
656 def _get_predict_end(self, end, dynamic=False):
657 # pass through so predict works for ARIMA and ARMA
658 return super(ARMA, self)._get_predict_end(end)
660 def geterrors(self, params):
661 """
662 Get the errors of the ARMA process.
664 Parameters
665 ----------
666 params : array-like
667 The fitted ARMA parameters
668 order : array-like
669 3 item iterable, with the number of AR, MA, and exogenous
670 parameters, including the trend
671 """
673 #start = self._get_predict_start(start) # will be an index of a date
674 #end, out_of_sample = self._get_predict_end(end)
675 params = np.asarray(params)
676 k_ar, k_ma = self.k_ar, self.k_ma
677 k = self.k_exog + self.k_trend
679 method = getattr(self, 'method', 'mle')
680 if 'mle' in method: # use KalmanFilter to get errors
681 (y, k, nobs, k_ar, k_ma, k_lags, newparams, Z_mat, m, R_mat,
682 T_mat, paramsdtype) = KalmanFilter._init_kalman_state(params,
683 self)
685 errors = KalmanFilter.geterrors(y, k, k_ar, k_ma, k_lags, nobs,
686 Z_mat, m, R_mat, T_mat,
687 paramsdtype)
688 if
isinstance(errors, tuple):
689 errors = errors[0] # non-cython version returns a tuple
690 else: # use scipy.signal.lfilter
691 y = self.endog.copy()
692 k = self.k_exog + self.k_trend
693 if k > 0:
694 y -= dot(self.exog, params[:k])
696 k_ar = self.k_ar
697 k_ma = self.k_ma
699 (trendparams, exparams,
700 arparams, maparams) = _unpack_params(params, (k_ar, k_ma),
701 self.k_trend, self.k_exog,
702 reverse=False)
703 b, a = np.r_[1, -arparams], np.r_[1, maparams]
704 zi = zeros((max(k_ar, k_ma)))
705 for i in range(k_ar):
706 zi[i] = sum(-b[:i+1][::-1]*y[:i+1])
707 e = lfilter(b, a, y, zi=zi)
708 errors = e[0][k_ar:]
709 return errors.squeeze()
711 def predict(self, params, start=None, end=None, exog=None, dynamic=False):
712 method = getattr(self, 'method', 'mle') # don't assume fit
713 #params = np.asarray(params)
715 # will return an index of a date
716 start = self._get_predict_start(start, dynamic)
717 end, out_of_sample = self._get_predict_end(end, dynamic)
718 if out_of_sample and (exog is None and self.k_exog > 0):
719 raise ValueError("You must provide exog for ARMAX")
721 endog = self.endog
722 resid = self.geterrors(params)
723 k_ar = self.k_ar
725 if out_of_sample != 0 and self.k_exog > 0:
726 if self.k_exog == 1 and exog.ndim == 1:
727 exog = exog[:, None]
728 # we need the last k_ar exog for the lag-polynomial
729 if self.k_exog > 0 and k_ar > 0:
730 # need the last k_ar exog for the lag-polynomial
731 exog = np.vstack((self.exog[-k_ar:, self.k_trend:], exog))
733 if dynamic:
734 #TODO: now that predict does dynamic in-sample it should
735 # also return error estimates and confidence intervals
736 # but how? len(endog) is not tot_obs
737 out_of_sample += end - start + 1
738 pr, ct = _arma_predict_out_of_sample(params, out_of_sample, resid,
739 k_ar, self.k_ma, self.k_trend,
740 self.k_exog, endog, exog,
741 start, method)
742 self.constant = ct
743 return pr
745 predictedvalues = _arma_predict_in_sample(start, end, endog, resid,
746 k_ar, method)
747 if out_of_sample:
748 forecastvalues, ct = _arma_predict_out_of_sample(params, out_of_sample,
749 resid, k_ar,
750 self.k_ma,
751 self.k_trend,
752 self.k_exog, endog,
753 exog, method=method)
754 self.constant = ct
755 predictedvalues = np.r_[predictedvalues, forecastvalues]
756 return predictedvalues
757 predict.__doc__ = _arma_predict
759 def loglike(self, params, set_sigma2=True):
760 """
761 Compute the log-likelihood for ARMA(p,q) model
763 Notes
764 -----
765 Likelihood used depends on the method set in fit
766 """
767 method = self.method
768 if method in ['mle', 'css-mle']:
769 return self.loglike_kalman(params, set_sigma2)
770 elif method == 'css':
771 return self.loglike_css(params, set_sigma2)
772 else:
773 raise ValueError("Method %s not understood" % method)
775 def loglike_kalman(self, params, set_sigma2=True):
776 """
777 Compute exact loglikelihood for ARMA(p,q) model by the Kalman Filter.
778 """
779 return KalmanFilter.loglike(params, self, set_sigma2)
781 def loglike_css(self, params, set_sigma2=True):
782 """
783 Conditional Sum of Squares likelihood function.
784 """
785 k_ar = self.k_ar
786 k_ma = self.k_ma
787 k = self.k_exog + self.k_trend
788 y = self.endog.copy().astype(params.dtype)
789 nobs = self.nobs
790 # how to handle if empty?
791 if self.transparams:
792 newparams = self._transparams(params)
793 else:
794 newparams = params
795 if k > 0:
796 y -= dot(self.exog, newparams[:k])
797 # the order of p determines how many zeros errors to set for lfilter
798 b, a = np.r_[1, -newparams[k:k + k_ar]], np.r_[1, newparams[k + k_ar:]]
799 zi = np.zeros((max(k_ar, k_ma)), dtype=params.dtype)
800 for i in range(k_ar):
801 zi[i] = sum(-b[:i + 1][::-1] * y[:i + 1])
802 errors = lfilter(b, a, y, zi=zi)[0][k_ar:]
804 ssr = np.dot(errors, errors)
805 sigma2 = ssr/nobs
806 if set_sigma2:
807 self.sigma2 = sigma2
808 llf = -nobs/2.*(log(2*pi) + log(sigma2)) - ssr/(2*sigma2)
809 return llf
811 def fit(self, start_params=None, trend='c', method="css-mle",
812 transparams=True, solver='lbfgs', maxiter=50, full_output=1,
813 disp=5, callback=None, **kwargs):
814 """
815 Fits ARMA(p,q) model using exact maximum likelihood via Kalman filter.
817 Parameters
818 ----------
819 start_params : array-like, optional
820 Starting parameters for ARMA(p,q). If None, the default is given
821 by ARMA._fit_start_params. See there for more information.
822 transparams : bool, optional
823 Whehter or not to transform the parameters to ensure stationarity.
824 Uses the transformation suggested in Jones (1980). If False,
825 no checking for stationarity or invertibility is done.
826 method : str {'css-mle','mle','css'}
827 This is the loglikelihood to maximize. If "css-mle", the
828 conditional sum of squares likelihood is maximized and its values
829 are used as starting values for the computation of the exact
830 likelihood via the Kalman filter. If "mle", the exact likelihood
831 is maximized via the Kalman Filter. If "css" the conditional sum
832 of squares likelihood is maximized. All three methods use
833 `start_params` as starting parameters. See above for more
834 information.
835 trend : str {'c','nc'}
836 Whether to include a constant or not. 'c' includes constant,
837 'nc' no constant.
838 solver : str or None, optional
839 Solver to be used. The default is 'lbfgs' (limited memory
840 Broyden-Fletcher-Goldfarb-Shanno). Other choices are 'bfgs',
841 'newton' (Newton-Raphson), 'nm' (Nelder-Mead), 'cg' -
842 (conjugate gradient), 'ncg' (non-conjugate gradient), and
843 'powell'. By default, the limited memory BFGS uses m=12 to
844 approximate the Hessian, projected gradient tolerance of 1e-8 and
845 factr = 1e2. You can change these by using kwargs.
846 maxiter : int, optional
847 The maximum number of function evaluations. Default is 50.
848 tol : float
849 The convergence tolerance. Default is 1e-08.
850 full_output : bool, optional
851 If True, all output from solver will be available in
852 the Results object's mle_retvals attribute. Output is dependent
853 on the solver. See Notes for more information.
854 disp : bool, optional
855 If True, convergence information is printed. For the default
856 l_bfgs_b solver, disp controls the frequency of the output during
857 the iterations. disp < 0 means no output in this case.
858 callback : function, optional
859 Called after each iteration as callback(xk) where xk is the current
860 parameter vector.
861 kwargs
862 See Notes for keyword arguments that can be passed to fit.
864 Returns
865 -------
866 statsmodels.tsa.arima_model.ARMAResults class
868 See also
869 --------
870 statsmodels.base.model.LikelihoodModel.fit : for more information
871 on using the solvers.
872 ARMAResults : results class returned by fit
874 Notes
875 ------
876 If fit by 'mle', it is assumed for the Kalman Filter that the initial
877 unkown state is zero, and that the inital variance is
878 P = dot(inv(identity(m**2)-kron(T,T)),dot(R,R.T).ravel('F')).reshape(r,
879 r, order = 'F')
881 """
882 k_ar = self.k_ar
883 k_ma = self.k_ma
885 # enforce invertibility
886 self.transparams = transparams
888 endog, exog = self.endog, self.exog
889 k_exog = self.k_exog
890 self.nobs = len(endog) # this is overwritten if method is 'css'
892 # (re)set trend and handle exogenous variables
893 # always pass original exog
894 k_trend, exog = _make_arma_exog(endog, self.exog, trend)
896 # Check has something to estimate
897 if k_ar == 0 and k_ma == 0 and k_trend == 0 and k_exog == 0:
898 raise ValueError("Estimation requires the inclusion of least one "
899 "AR term, MA term, a constant or an exogenous "
900 "variable.")
902 # check again now that we know the trend
903 _check_estimable(len(endog), k_ar + k_ma + k_exog + k_trend)
905 self.k_trend = k_trend
906 self.exog = exog # overwrites original exog from __init__
908 # (re)set names for this model
909 self.exog_names = _make_arma_names(self.data, k_trend, (k_ar, k_ma),
910 self.exog_names)
911 k = k_trend + k_exog
913 # choose objective function
914 if k_ma == 0 and k_ar == 0:
915 method = "css" # Always CSS when no AR or MA terms
917 self.method = method = method.lower()
919 # adjust nobs for css
920 if method == 'css':
921 self.nobs = len(self.endog) - k_ar
923 if start_params is not None:
924 start_params = np.asarray(start_params)
926 else: # estimate starting parameters
927 start_params = self._fit_start_params((k_ar, k_ma, k), method)
929 if transparams: # transform initial parameters to ensure invertibility
930 start_params = self._invtransparams(start_params)
932 if solver == 'lbfgs':
933 kwargs.setdefault('pgtol', 1e-8)
934 kwargs.setdefault('factr', 1e2)
935
kwargs.setdefault('m', 12)
936 kwargs.setdefault('approx_grad', True)
937 mlefit = super(ARMA, self).fit(start_params, method=solver,
938 maxiter=maxiter,
939 full_output=full_output, disp=disp,
940 callback=callback, **kwargs)
941 params = mlefit.params
943 if transparams: # transform parameters back
944 params = self._transparams(params)
946 self.transparams = False # so methods don't expect transf.
948 normalized_cov_params = None # TODO: fix this
949 armafit = ARMAResults(self, params, normalized_cov_params)
950 armafit.mle_retvals = mlefit.mle_retvals
951 armafit.mle_settings = mlefit.mle_settings
952 armafit.mlefit = mlefit
953 return ARMAResultsWrapper(armafit)
956 #NOTE: the length of endog changes when we give a difference to fit
957 #so model methods are not the same on unfit models as fit ones
958 #starting to think that order of model should be put in instantiation...
959 class ARIMA(ARMA):
961 __doc__ = tsbase._tsa_doc % {"model" : _arima_model,
962 "params" : _arima_params, "extra_params" : "",
963 "extra_sections" : _armax_notes %
964 {"Model" : "ARIMA"}}
966 def __new__(cls, endog, order, exog=None, dates=None, freq=None,
967 missing='none'):
968 p, d, q = order
969 if d == 0: # then we just use an ARMA model
970 return ARMA(endog, (p, q), exog, dates, freq, missing)
971 else:
972 mod = super(ARIMA, cls).__new__(cls)
973 mod.__init__(endog, order, exog, dates, freq, missing)
974 return mod
976 def __init__(self, endog, order, exog=None, dates=None, freq=None,
977 missing='none'):
978 p, d, q = order
979 if d > 2:
980 #NOTE: to make more general, need to address the d == 2 stuff
981 # in the predict method
982 raise ValueError("d > 2 is not supported")
983 super(ARIMA, self).__init__(endog, (p, q), exog, dates, freq, missing)
984 self.k_diff = d
985 self._first_unintegrate = unintegrate_levels(self.endog[:d], d)
986 self.endog = np.diff(self.endog, n=d)
987 #NOTE: will check in ARMA but check again since differenced now
988 _check_estimable(len(self.endog), p+q)
989 if exog is not None:
990 self.exog = self.exog[d:]
991 if d == 1:
992 self.data.ynames = 'D.' + self.endog_names
993 else:
994 self.data.ynames = 'D{0:d}.'.format(d) + self.endog_names
995 # what about exog, should we difference it automatically before
996 # super call?
998 def _get_predict_start(self, start, dynamic):
999 """
1000 """
1001 #TODO: remove all these getattr and move order specification to
1002 # class constructor
1003 k_diff = getattr(self, 'k_diff', 0)
1004 method = getattr(self, 'method', 'mle')
1005 k_ar = getattr(self, 'k_ar', 0)
1006 if start is None:
1007 if 'mle' in method and not dynamic:
1008 start = 0
1009 else:
1010 start = k_ar
1011 elif isinstance(start, int):
1012 start -= k_diff
1013 try: # catch when given an integer outside of dates index
1014 start = super(ARIMA, self)._get_predict_start(start,
1015 dynamic)
1016 except IndexError:
1017 raise ValueError("start must be in series. "
1018 "got %d" % (start + k_diff))
1019 else: # received a date
1020 start = _validate(start, k_ar, k_diff, self.data.dates,
1021 method)
1022 start = super(ARIMA, self)._get_predict_start(start, dynamic)
1023 # reset date for k_diff adjustment
1024 self._set_predict_start_date(start + k_diff)
1025 return start
1027 def _get_predict_end(self, end, dynamic=False):
1028 """
1029 Returns last index to be forecast of the differenced array.
1030 Handling of inclusiveness should be done in the predict function.
1031 """
1032 end, out_of_sample = super(ARIMA, self)._get_predict_end(end, dynamic)
1033 if 'mle' not in self.method and not dynamic:
1034 end -= self.k_ar
1036
return end - self.k_diff, out_of_sample
1038 def fit(self, start_params=None, trend='c', method="css-mle",
1039 transparams=True, solver='lbfgs', maxiter=50, full_output=1,
1040 disp=5, callback=None, **kwargs):
1041 """
1042 Fits ARIMA(p,d,q) model by exact maximum likelihood via Kalman filter.
1044 Parameters
1045 ----------
1046 start_params : array-like, optional
1047 Starting parameters for ARMA(p,q). If None, the default is given
1048 by ARMA._fit_start_params. See there for more information.
1049 transparams : bool, optional
1050 Whehter or not to transform the parameters to ensure stationarity.
1051 Uses the transformation suggested in Jones (1980). If False,
1052 no checking for stationarity or invertibility is done.
1053 method : str {'css-mle','mle','css'}
1054 This is the loglikelihood to maximize. If "css-mle", the
1055 conditional sum of squares likelihood is maximized and its values
1056 are used as starting values for the computation of the exact
1057 likelihood via the Kalman filter. If "mle", the exact likelihood
1058 is maximized via the Kalman Filter. If "css" the conditional sum
1059 of squares likelihood is maximized. All three methods use
1060 `start_params` as starting parameters. See above for more
1061 information.
1062 trend : str {'c','nc'}
1063 Whether to include a constant or not. 'c' includes constant,
1064 'nc' no constant.
1065 solver : str or None, optional
1066 Solver to be used. The default is 'lbfgs' (limited memory
1067 Broyden-Fletcher-Goldfarb-Shanno). Other choices are 'bfgs',
1068 'newton' (Newton-Raphson), 'nm' (Nelder-Mead), 'cg' -
1069 (conjugate gradient), 'ncg' (non-conjugate gradient), and
1070 'powell'. By default, the limited memory BFGS uses m=12 to
1071 approximate the Hessian, projected gradient tolerance of 1e-8 and
1072 factr = 1e2. You can change these by using kwargs.
1073 maxiter : int, optional
1074 The maximum number of function evaluations. Default is 50.
1075 tol : float
1076 The convergence tolerance. Default is 1e-08.
1077 full_output : bool, optional
1078 If True, all output from solver will be available in
1079 the Results object's mle_retvals attribute. Output is dependent
1080 on the solver. See Notes for more information.
1081 disp : bool, optional
1082 If True, convergence information is printed. For the default
1083 l_bfgs_b solver, disp controls the frequency of the output during
1084 the iterations. disp < 0 means no output in this case.
1085 callback : function, optional
1086 Called after each iteration as callback(xk) where xk is the current
1087 parameter vector.
1088 kwargs
1089 See Notes for keyword arguments that can be passed to fit.
1091 Returns
1092 -------
1093 `statsmodels.tsa.arima.ARIMAResults` class
1095 See also
1096 --------
1097 statsmodels.base.model.LikelihoodModel.fit : for more information
1098 on using the solvers.
1099 ARIMAResults : results class returned by fit
1101 Notes
1102 ------
1103 If fit by 'mle', it is assumed for the Kalman Filter that the initial
1104 unkown state is zero, and that the inital variance is
1105 P = dot(inv(identity(m**2)-kron(T,T)),dot(R,R.T).ravel('F')).reshape(r,
1106 r, order = 'F')
1108 """
1109 arima_fit = super(ARIMA, self).fit(start_params, trend,
1110 method, transparams, solver,
1111 maxiter, full_output, disp,
1112 callback, **kwargs)
1113 normalized_cov_params = None # TODO: fix this?
1114 arima_fit = ARIMAResults(self, arima_fit._results.params,
1115 normalized_cov_params)
1116 arima_fit.k_diff = self.k_diff
1117 return ARIMAResultsWrapper(arima_fit)
1119 def predict(self, params, start=None, end=None, exog=None, typ='linear',
1120 dynamic=False):
1121 # go ahead and convert to an index for easier checking
1122 if isinstance(start, (string_types, datetime)):
1123 start = _index_date(start, self.data.dates)
1124 if typ == 'linear':
1125 if not dynamic or (start != self.k_ar + self.k_diff and
1126 start is not None):
1127 return super(ARIMA, self).predict(params, start, end, exog,
1128 dynamic)
1129 else:
1130 # need to assume pre-sample residuals are zero
1131 # do this by a hack
1132 q = self.k_ma
1133 self.k_ma = 0
1134 predictedvalues = super(ARIMA, self).predict(params, start,
1135 end, exog,
1136 dynamic)
1137 self.k_ma = q
1138 return predictedvalues
1139 elif typ == 'levels':
1140 endog = self.data.endog
1141 if not dynamic:
1142 predict = super(ARIMA, self).predict(params, start, end,
1143 dynamic)
1145 start = self._get_predict_start(start, dynamic)
1146 end, out_of_sample = self._get_predict_end(end)
1147 d = self.k_diff
1148 if 'mle' in self.method:
1149 start += d - 1 # for case where d == 2
1150 end += d - 1
1151 # add each predicted diff to lagged endog
1152 if out_of_sample:
1153 fv = predict[:-out_of_sample] + endog[start:end+1]
1154 if d == 2: #TODO: make a general solution to this
1155 fv += np.diff(endog[start - 1:end + 1])
1156 levels = unintegrate_levels(endog[-d:], d)
1157 fv = np.r_[fv,
1158 unintegrate(predict[-out_of_sample:],
1159 levels)[d:]]
1160 else:
1161 fv = predict + endog[start:end + 1]
1162 if d == 2:
1163 fv += np.diff(endog[start - 1:end + 1])
1164 else:
1165 k_ar = self.k_ar
1166 if out_of_sample:
1167 fv = (predict[:-out_of_sample] +
1168 endog[max(start, self.k_ar-1):end+k_ar+1])
1169 if d == 2:
1170 fv += np.diff(endog[start - 1:end + 1])
1171 levels = unintegrate_levels(endog[-d:], d)
1172 fv = np.r_[fv,
1173 unintegrate(predict[-out_of_sample:],
1174 levels)[d:]]
1175 else:
1176 fv = predict + endog[max(start, k_ar):end+k_ar+1]
1177 if d == 2:
1178 fv += np.diff(endog[start - 1:end + 1])
1179 else:
1180 #IFF we need to use pre-sample values assume pre-sample
1181 # residuals are zero, do this by a hack
1182 if start == self.k_ar + self.k_diff or start is None:
1183 # do the first k_diff+1 separately
1184 p = self.k_ar
1185 q = self.k_ma
1186 k_exog = self.k_exog
1187 k_trend = self.k_trend
1188 k_diff = self.k_diff
1189 (trendparam, exparams,
1190 arparams, maparams) = _unpack_params(params, (p, q),
1191 k_trend,
1192 k_exog,
1193 reverse=True)
1194 # this is the hack
1195 self.k_ma = 0
1197 predict = super(ARIMA, self).predict(params, start, end,
1198 exog, dynamic)
1199 if not start:
1200 start = self._get_predict_start(start, dynamic)
1201 start += k_diff
1202 self.k_ma = q
1203 return endog[start-1] + np.cumsum(predict)
1204 else:
1205 predict = super(ARIMA, self).predict(params, start, end,
1206 exog, dynamic)
1207 return endog[start-1] + np.cumsum(predict)
1208 return fv
1210 else: # pragma : no cover
1211 raise ValueError("typ %s not understood" % typ)
1213 predict.__doc__ = _arima_predict
1216 class ARMAResults(tsbase.TimeSeriesModelResults):
1217 """
1218 Class to hold results from fitting an ARMA model.
1220 Parameters
1221 ----------
1222 model : ARMA instance
1223 The fitted model instance
1224 params : array
1225 Fitted parameters
1226 normalized_cov_params : array, optional
1227 The normalized variance covariance matrix
1228 scale : float, optional
1229 Optional argument to scale the variance covariance matrix.
1231 Returns
1232 --------
1233 **Attributes**
1235 aic : float
1236 Akaike Information Criterion
1237 :math:`-2*llf+2* df_model`
1238 where `df_model` includes all AR parameters, MA parameters, constant
1239 terms parameters on constant terms and the variance.
1240 arparams : array
1241 The parameters associated with the AR coefficients in the model.
1242 arroots : array
1243 The roots of the AR coefficients are the solution to
1244 (1 - arparams[0]*z - arparams[1]*z**2 -...- arparams[p-1]*z**k_ar) = 0
1245 Stability requires that the roots in modulus lie outside the unit
1246 circle.
1247 bic : float
1248 Bayes Information Criterion
1249 -2*llf + log(nobs)*df_model
1250 Where if the model is fit using conditional sum of squares, the
1251 number of observations `nobs` does not include the `p` pre-sample
1252 observations.
1253 bse : array
1254 The standard errors of the parameters. These are computed using the
1255 numerical Hessian.
1256 df_model : array
1257 The model degrees of freedom = `k_exog` + `k_trend` + `k_ar` + `k_ma`
1258 df_resid : array
1259 The residual degrees of freedom = `nobs` - `df_model`
1260 fittedvalues : array
1261 The predicted values of the model.
1262 hqic : float
1263 Hannan-Quinn Information Criterion
1264 -2*llf + 2*(`df_model`)*log(log(nobs))
1265 Like `bic` if the model is fit using conditional sum of squares then
1266 the `k_ar` pre-sample observations are not counted in `nobs`.
1267 k_ar : int
1268 The number of AR coefficients in the model.
1269 k_exog : int
1270 The number of exogenous variables included in the model. Does not
1271 include the constant.
1272 k_ma : int
1273 The number of MA coefficients.
1274 k_trend : int
1275 This is 0 for no constant or 1 if a constant is included.
1276 llf : float
1277 The value of the log-likelihood function evaluated at `params`.
1278 maparams : array
1279 The value of the moving average coefficients.
1280 maroots : array
1281 The roots of the MA coefficients are the solution to
1282 (1 + maparams[0]*z + maparams[1]*z**2 + ... + maparams[q-1]*z**q) = 0
1283 Stability requires that the roots in modules lie outside the unit
1284 circle.
1285 model : ARMA instance
1286 A reference to the model that was fit.
1287 nobs : float
1288 The number of observations used to fit the model. If the model is fit
1289 using exact maximum likelihood this is equal to the total number of
1290 observations, `n_totobs`. If the model is fit using conditional
1291 maximum likelihood this is equal to `n_totobs` - `k_ar`.
1292 n_totobs : float
1293 The total number of observations for `endog`. This includes all
1294 observations, even pre-sample values if the model is fit using `css`.
1295 params : array
1296 The parameters of the model. The order of variables is the trend
1297 coefficients and the `k_exog` exognous coefficients, then the
1298 `k_ar` AR coefficients, and finally the `k_ma` MA coefficients.
1299 pvalues : array
1300 The p-values associated with the t-values of the coefficients. Note
1301 that the coefficients are assumed to have a Student's T distribution.
1302 resid : array
1303 The model residuals. If the model is fit using 'mle' then the
1304 residuals are created via the Kalman Filter. If the model is fit
1305 using 'css' then the residuals are obtained via `scipy.signal.lfilter`
1306 adjusted such that the first `k_ma` residuals are zero. These zero
1307 residuals are not returned.
1308 scale : float
1309 This is currently set to 1.0 and not used by the model or its results.
1310 sigma2 : float
1311 The variance of the residuals. If the model is fit by 'css',
1312 sigma2 = ssr/nobs, where ssr is the sum of squared residuals. If
1313 the model is fit by 'mle', then sigma2 = 1/nobs * sum(v**2 / F)
1314 where v is the one-step forecast error and F is the forecast error
1315 variance. See `nobs` for the difference in definitions depending on the
1316 fit.
1317 """
1318 _cache = {}
1320 #TODO: use this for docstring when we fix nobs issue
1322 def __init__(self, model, params, normalized_cov_params=None, scale=1.):
1323 super(ARMAResults, self).__init__(model, params, normalized_cov_params,
1324 scale)
1325 self.sigma2 = model.sigma2
1326 nobs = model.nobs
1327 self.nobs = nobs
1328 k_exog = model.k_exog
1329 self.k_exog = k_exog
1330 k_trend = model.k_trend
1331 self.k_trend = k_trend
1332 k_ar = model.k_ar
1333 self.k_ar = k_ar
1334 self.n_totobs = len(model.endog)
1335 k_ma = model.k_ma
1336 self.k_ma = k_ma
1337 df_model = k_exog + k_trend + k_ar + k_ma
1338 self._ic_df_model = df_model + 1
1339 self.df_model = df_model
1340 self.df_resid = self.nobs - df_model
1341 self._cache = resettable_cache()
1342 self.constant = 0 #Added by me
1344 @cache_readonly
1345 def arroots(self):
1346 return np.roots(np.r_[1, -self.arparams])**-1
1348 @cache_readonly
1349 def maroots(self):
1350 return np.roots(np.r_[1, self.maparams])**-1
1352 @cache_readonly
1353 def arfreq(self):
1354 r"""
1355 Returns the frequency of the AR roots.
1357 This is the solution, x, to z = abs(z)*exp(2j*np.pi*x) where z are the
1358 roots.
1359 """
1360 z = self.arroots
1361 if not z.size:
1362 return
1363 return np.arctan2(z.imag, z.real) / (2*pi)
1365 @cache_readonly
1366 def mafreq(self):
1367 r"""
1368 Returns the frequency of the MA roots.
1370 This is the solution, x, to z = abs(z)*exp(2j*np.pi*x) where z are the
1371 roots.
1372 """
1373 z = self.maroots
1374 if not z.size:
1375 return
1376 return np.arctan2(z.imag, z.real) / (2*pi)
1378 @cache_readonly
1379 def arparams(self):
1380 k = self.k_exog + self.k_trend
1381 return self.params[k:k+self.k_ar]
1383 @cache_readonly
1384 def maparams(self):
1385 k = self.k_exog + self.k_trend
1386 k_ar = self.k_ar
1387 return self.params[k+k_ar:]
1389 @cache_readonly
1390 def llf(self):
1391 return self.model.loglike(self.params)
1393 @cache_readonly
1394 def bse(self):
1395 params = self.params
1396 hess = self.model.hessian(params)
1397 if len(params) == 1: # can't take an inverse, ensure 1d
1398 return np.sqrt(-1./hess[0])
1399 return np.sqrt(np.diag(-inv(hess)))
1401 def cov_params(self): # add scale argument?
1402 params = self.params
1403 hess = self.model.hessian(params)
1404 return -inv(hess)
1406 @cache_readonly
1407 def aic(self):
1408 return -2 * self.llf + 2 * self._ic_df_model
1410 @cache_readonly
1411 def bic(self):
1412 nobs = self.nobs
1413 return -2 * self.llf + np.log(nobs) * self._ic_df_model
1415 @cache_readonly
1416 def hqic(self):
1417 nobs = self.nobs
1418 return -2 * self.llf + 2 * np.log(np.log(nobs)) * self._ic_df_model
1420 @cache_readonly
1421 def fittedvalues(self):
1422 model = self.model
1423
endog = model.endog.copy()
1424 k_ar = self.k_ar
1425 exog = model.exog # this is a copy
1426 if exog is not None:
1427 if model.method == "css" and k_ar > 0:
1428 exog = exog[k_ar:]
1429 if model.method == "css" and k_ar > 0:
1430 endog = endog[k_ar:]
1431 fv = endog - self.resid
1432 # add deterministic part back in
1433 #k = self.k_exog + self.k_trend
1434 #TODO: this needs to be commented out for MLE with constant
1435 #if k != 0:
1436 # fv += dot(exog, self.params[:k])
1437 return fv
1439 @cache_readonly
1440 def resid(self):
1441 return self.model.geterrors(self.params)
1443 @cache_readonly
1444 def pvalues(self):
1445 #TODO: same for conditional and unconditional?
1446 df_resid = self.df_resid
1447 return t.sf(np.abs(self.tvalues), df_resid) * 2
1449 def predict(self, start=None, end=None, exog=None, dynamic=False):
1450 return self.model.predict(self.params, start, end, exog, dynamic)
1451 predict.__doc__ = _arma_results_predict
1453 def _forecast_error(self, steps):
1454 sigma2 = self.sigma2
1455 ma_rep = arma2ma(np.r_[1, -self.arparams],
1456 np.r_[1, self.maparams], nobs=steps)
1458 fcasterr = np.sqrt(sigma2 * np.cumsum(ma_rep**2))
1459 return fcasterr
1461 def _forecast_conf_int(self, forecast, fcasterr, alpha):
1462 const = norm.ppf(1 - alpha / 2.)
1463 conf_int = np.c_[forecast - const * fcasterr,
1464 forecast + const * fcasterr]
1466 return conf_int
1468 def forecast(self, steps=1, exog=None, alpha=.05):
1469 """
1470 Out-of-sample forecasts
1472 Parameters
1473 ----------
1474 steps : int
1475 The number of out of sample forecasts from the end of the
1476 sample.
1477 exog : array
1478 If the model is an ARMAX, you must provide out of sample
1479 values for the exogenous variables. This should not include
1480 the constant.
1481 alpha : float
1482 The confidence intervals for the forecasts are (1 - alpha) %
1484 Returns
1485 -------
1486 forecast : array
1487 Array of out of sample forecasts
1488 stderr : array
1489 Array of the standard error of the forecasts.
1490 conf_int : array
1491 2d array of the confidence interval for the forecast
1492 """
1493 if exog is not None:
1494 #TODO: make a convenience function for this. we're using the
1495 # pattern elsewhere in the codebase
1496 exog = np.asarray(exog)
1497 if self.k_exog == 1 and exog.ndim == 1:
1498 exog = exog[:, None]
1499 elif exog.ndim == 1:
1500 if len(exog) != self.k_exog:
1501 raise ValueError("1d exog given and len(exog) != k_exog")
1502 exog = exog[None, :]
1503 if exog.shape[0] != steps:
1504 raise ValueError("new exog needed for each step")
1505 # prepend in-sample exog observations
1506 exog = np.vstack((self.model.exog[-self.k_ar:, self.k_trend:],
1507 exog))
1509 forecast, ct = _arma_predict_out_of_sample(self.params,
1510 steps, self.resid, self.k_ar,
1511 self.k_ma, self.k_trend,
1512 self.k_exog, self.model.endog,
1513 exog, method=self.model.method)
1514 self.constant = ct
1516 # compute the standard errors
1517 fcasterr = self._forecast_error(steps)
1518 conf_int = self._forecast_conf_int(forecast, fcasterr, alpha)
1520 return forecast, fcasterr, conf_int
1522 def summary(self, alpha=.05):
1523 """Summarize the Model
1525 Parameters
1526 ----------
1527 alpha : float, optional
1528 Significance level for the confidence intervals.
1530 Returns
1531 -------
1532 smry : Summary instance
1533 This holds the summary table and text, which can be printed or
1534 converted to various output formats.
1536 See Also
1537 --------
1538 statsmodels.iolib.summary.Summary
1539 """
1540 from statsmodels.iolib.summary import Summary
1541 model = self.model
1542 title = model.__class__.__name__ + ' Model Results'
1543 method = model.method
1544 # get sample TODO: make better sample machinery for estimation
1545 k_diff = getattr(self, 'k_diff', 0)
1546 if 'mle' in method:
1547 start = k_diff
1548 else:
1549 start = k_diff + self.k_ar
1550 if self.data.dates is
not None:
1551 dates = self.data.dates
1552 sample = [dates[start].strftime('%m-%d-%Y')]
1553 sample += ['- ' + dates[-1].strftime('%m-%d-%Y')]
1554 else:
1555 sample = str(start) + ' - ' + str(len(self.data.orig_endog))
1557 k_ar, k_ma = self.k_ar, self.k_ma
1558 if not k_diff:
1559 order = str((k_ar, k_ma))
1560 else:
1561 order = str((k_ar, k_diff, k_ma))
1562 top_left = [('Dep. Variable:', None),
1563 ('Model:', [model.__class__.__name__ + order]),
1564 ('Method:', [method]),
1565 ('Date:', None),
1566 ('Time:', None),
1567 ('Sample:', [sample[0]]),
1568 ('', [sample[1]])
1569 ]
1571 top_right = [
1572 ('No. Observations:', [str(len(self.model.endog))]),
1573 ('Log Likelihood', ["%#5.3f" % self.llf]),
1574 ('S.D. of innovations', ["%#5.3f" % self.sigma2**.5]),
1575 ('AIC', ["%#5.3f" % self.aic]),
1576 ('BIC', ["%#5.3f" % self.bic]),
1577 ('HQIC', ["%#5.3f" % self.hqic])]
1579 smry = Summary()
1580 smry.add_table_2cols(self, gleft=top_left, gright=top_right,
1581 title=title)
1582 smry.add_table_params(self, alpha=alpha, use_t=False)
1584 # Make the roots table
1585 from statsmodels.iolib.table import SimpleTable
1587 if k_ma and k_ar:
1588 arstubs = ["AR.%d" % i for i in range(1, k_ar + 1)]
1589 mastubs = ["MA.%d" % i for i in range(1, k_ma + 1)]
1590 stubs = arstubs + mastubs
1591 roots = np.r_[self.arroots, self.maroots]
1592 freq = np.r_[self.arfreq, self.mafreq]
1593 elif k_ma:
1594 mastubs = ["MA.%d" % i for i in range(1, k_ma + 1)]
1595 stubs = mastubs
1596 roots = self.maroots
1597 freq = self.mafreq
1598 elif k_ar:
1599 arstubs = ["AR.%d" % i for i in range(1, k_ar + 1)]
1600 stubs = arstubs
1601 roots = self.arroots
1602 freq = self.arfreq
1603 else: # 0,0 model
1604 stubs = []
1605 if len(stubs): # not 0, 0
1606 modulus = np.abs(roots)
1607 data = np.column_stack((roots.real, roots.imag, modulus, freq))
1608 roots_table = SimpleTable(data,
1609 headers=[' Real',
1610 ' Imaginary',
1611 ' Modulus',
1612 ' Frequency'],
1613 title="Roots",
1614 stubs=stubs,
1615 data_fmts=["%17.4f", "%+17.4fj",
1616 "%17.4f", "%17.4f"])
1618 smry.tables.append(roots_table)
1619 return smry
1621 def summary2(self, title=None, alpha=.05, float_format="%.4f"):
1622 """Experimental summary function for ARIMA Results
1624 Parameters
1625 -----------
1626 title : string, optional
1627 Title for the top table. If not None, then this replaces the
1628 default title
1629 alpha : float
1630 significance level for the confidence intervals
1631 float_format: string
1632 print format for floats in parameters summary
1634 Returns
1635 -------
1636 smry : Summary instance
1637 This holds the summary table and text, which can be printed or
1638 converted to various output formats.
1640 See Also
1641 --------
1642 statsmodels.iolib.summary2.Summary : class to hold summary
1643 results
1645 """
1646 from pandas import DataFrame
1647 # get sample TODO: make better sample machinery for estimation
1648
k_diff = getattr(self, 'k_diff', 0)
1649 if 'mle' in self.model.method:
1650 start = k_diff
1651 else:
1652 start = k_diff + self.k_ar
1653 if self.data.dates is not None:
1654 dates = self.data.dates
1655 sample = [dates[start].strftime('%m-%d-%Y')]
1656 sample += [dates[-1].strftime('%m-%d-%Y')]
1657 else:
1658 sample = str(start) + ' - ' + str(len(self.data.orig_endog))
1660 k_ar, k_ma = self.k_ar, self.k_ma
1662 # Roots table
1663 if k_ma and k_ar:
1664 arstubs = ["AR.%d" % i for i in range(1, k_ar + 1)]
1665 mastubs = ["MA.%d" % i for i in range(1, k_ma + 1)]
1666 stubs = arstubs + mastubs
1667 roots = np.r_[self.arroots, self.maroots]
1668 freq = np.r_[self.arfreq, self.mafreq]
1669 elif k_ma:
1670 mastubs = ["MA.%d" % i for i in range(1, k_ma + 1)]
1671 stubs = mastubs
1672 roots = self.maroots
1673 freq = self.mafreq
1674 elif k_ar:
1675 arstubs = ["AR.%d" % i for i in range(1, k_ar + 1)]
1676 stubs = arstubs
1677 roots = self.arroots
1678 freq = self.arfreq
1679 else: # 0, 0 order
1680 stubs = []
1682 if len(stubs):
1683 modulus = np.abs(roots)
1684 data = np.column_stack((roots.real, roots.imag, modulus, freq))
1685 data = DataFrame(data)
1686 data.columns = ['Real', 'Imaginary', 'Modulus', 'Frequency']
1687 data.index = stubs
1689 # Summary
1690 from statsmodels.iolib import summary2
1691 smry = summary2.Summary()
1693 # Model info
1694 model_info = summary2.summary_model(self)
1695 model_info['Method:'] = self.model.method
1696 model_info['Sample:'] = sample[0]
1697 model_info[' '] = sample[-1]
1698 model_info['S.D. of innovations:'] = "%#5.3f" % self.sigma2**.5
1699 model_info['HQIC:'] = "%#5.3f" % self.hqic
1700 model_info['No. Observations:'] = str(len(self.model.endog))
1702 # Parameters
1703 params = summary2.summary_params(self)
1704 smry.add_dict(model_info)
1705 smry.add_df(params, float_format=float_format)
1706 if len(stubs):
1707 smry.add_df(data, float_format="%17.4f")
1708 smry.add_title(results=self, title=title)
1710 return smry
1712 def plot_predict(self, start=None, end=None, exog=None, dynamic=False,
1713 alpha=.05, plot_insample=True, ax=None):
1714 from statsmodels.graphics.utils import _import_mpl, create_mpl_ax
1715 _ = _import_mpl()
1716 fig, ax = create_mpl_ax(ax)
1719 # use predict so you set dates
1720 forecast = self.predict(start, end, exog, dynamic)
1721 # doing this twice. just add a plot keyword to predict?
1722 start = self.model._get_predict_start(start, dynamic=False)
1723 end, out_of_sample = self.model._get_predict_end(end, dynamic=False)
1725 if out_of_sample:
1726 steps = out_of_sample
1727 fc_error = self._forecast_error(steps)
1728 conf_int = self._forecast_conf_int(forecast[-steps:], fc_error,
1729 alpha)
1732 if hasattr(self.data, "predict_dates"):
1733 from pandas import TimeSeries
1734 forecast = TimeSeries(forecast, index=self.data.predict_dates)
1735 ax = forecast.plot(ax=ax, label='forecast')
1736 else:
1737 ax.plot(forecast)
1739 x = ax.get_lines()[-1].get_xdata()
1740 if out_of_sample:
1741 label = "{0:.0%} confidence interval".format(1 - alpha)
1742 ax.fill_between(x[-out_of_sample:], conf_int[:, 0], conf_int[:, 1],
1743 color='gray', alpha=.5, label=label)
1745 if plot_insample:
1746 ax.plot(x[:end + 1 - start], self.model.endog[start:end+1],
1747 label=self.model.endog_names)
1749 ax.legend(loc='best')
1751 return fig
1752 plot_predict.__doc__ = _plot_predict
1755 class ARMAResultsWrapper(wrap.ResultsWrapper):
1756 _attrs = {}
1757 _wrap_attrs = wrap.union_dicts(tsbase.TimeSeriesResultsWrapper._wrap_attrs,
1758 _attrs)
1759
_methods = {}
1760 _wrap_methods = wrap.union_dicts(tsbase.TimeSeriesResultsWrapper._wrap_methods,
1761 _methods)
1762 wrap.populate_wrapper(ARMAResultsWrapper, ARMAResults)
1765 class ARIMAResults(ARMAResults):
1766 def predict(self, start=None, end=None, exog=None, typ='linear',
1767 dynamic=False):
1768 return self.model.predict(self.params, start, end, exog, typ, dynamic)
1769 predict.__doc__ = _arima_results_predict
1771 def _forecast_error(self, steps):
1772 sigma2 = self.sigma2
1773 ma_rep = arma2ma(np.r_[1, -self.arparams],
1774 np.r_[1, self.maparams], nobs=steps)
1776 fcerr = np.sqrt(np.cumsum(cumsum_n(ma_rep, self.k_diff)**2)*sigma2)
1777 return fcerr
1779 def _forecast_conf_int(self, forecast, fcerr, alpha):
1780 const = norm.ppf(1 - alpha/2.)
1781 conf_int = np.c_[forecast - const*fcerr, forecast + const*fcerr]
1782 return conf_int
1784 def forecast(self, steps=1, exog=None, alpha=.05):
1785 """
1786 Out-of-sample forecasts
1788 Parameters
1789 ----------
1790 steps : int
1791 The number of out of sample forecasts from the end of the
1792 sample.
1793 exog : array
1794 If the model is an ARIMAX, you must provide out of sample
1795 values for the exogenous variables. This should not include
1796 the constant.
1797 alpha : float
1798 The confidence intervals for the forecasts are (1 - alpha) %
1800 Returns
1801 -------
1802 forecast : array
1803 Array of out of sample forecasts
1804 stderr : array
1805 Array of the standard error of the forecasts.
1806 conf_int : array
1807 2d array of the confidence interval for the forecast
1809 Notes
1810 -----
1811 Prediction is done in the levels of the original endogenous variable.
1812 If you would like prediction of differences in levels use `predict`.
1813 """
1814 if exog is not None:
1815 if self.k_exog == 1 and exog.ndim == 1:
1816 exog = exog[:, None]
1817 if exog.shape[0] != steps:
1818 raise ValueError("new exog needed for each step")
1819 # prepend in-sample exog observations
1820 exog = np.vstack((self.model.exog[-self.k_ar:, self.k_trend:],
1821 exog))
1822 forecast, ct = _arma_predict_out_of_sample(self.params, steps, self.resid,
1823 self.k_ar, self.k_ma,
1824 self.k_trend, self.k_exog,
1825 self.model.endog,
1826 exog, method=self.model.method)
1828 #self.constant = ct
1829 d = self.k_diff
1830 endog = self.model.data.endog[-d:]
1831 forecast = unintegrate(forecast, unintegrate_levels(endog, d))[d:]
1833 # get forecast errors
1834 fcerr = self._forecast_error(steps)
1835 conf_int = self._forecast_conf_int(forecast, fcerr, alpha)
1836 return forecast, fcerr, conf_int
1838 def plot_predict(self, start=None, end=None, exog=None, dynamic=False,
1839 alpha=.05, plot_insample=True, ax=None):
1840 from statsmodels.graphics.utils import _import_mpl, create_mpl_ax
1841 _ = _import_mpl()
1842 fig, ax = create_mpl_ax(ax)
1844 # use predict so you set dates
1845 forecast = self.predict(start, end, exog, 'levels', dynamic)
1846 # doing this twice. just add a plot keyword to predict?
1847 start = self.model._get_predict_start(start, dynamic=dynamic)
1848 end, out_of_sample = self.model._get_predict_end(end, dynamic=dynamic)
1850 if out_of_sample:
1851 steps = out_of_sample
1852 fc_error = self._forecast_error(steps)
1853 conf_int = self._forecast_conf_int(forecast[-steps:], fc_error,
1854 alpha)
1856 if hasattr(self.data, "predict_dates"):
1857 from pandas import TimeSeries
1858 forecast = TimeSeries(forecast, index=self.data.predict_dates)
1859 ax = forecast.plot(ax=ax, label='forecast')
1860 else:
1861 ax.plot(forecast)
1863 x = ax.get_lines()[-1].get_xdata()
1864 if out_of_sample:
1865 label = "{0:.0%} confidence interval".format(1 - alpha)
1866 ax.fill_between(x[-out_of_sample:], conf_int[:, 0], conf_int[:, 1],
1867 color='gray', alpha=.5, label=label)
1869 if plot_insample:
1870 import re
1871 k_diff = self.k_diff
1872 label = re.sub("D\d*\.", "", self.model.endog_names)
1873 levels = unintegrate(self.model.endog,
1874 self.model._first_unintegrate)
1875 ax.plot(x[:end + 1 - start],
1876 levels[start + k_diff:end + k_diff + 1], label=label)
1878 ax.legend(loc='best')
1880 return fig
1882 plot_predict.__doc__ = _arima_plot_predict
1885 class ARIMAResultsWrapper(ARMAResultsWrapper):
1886 pass
1887 wrap.populate_wrapper(ARIMAResultsWrapper, ARIMAResults)
1890 if __name__ == "__main__":
1891 import statsmodels.api as sm
1893 # simulate arma process
1894 from statsmodels.tsa.arima_process import arma_generate_sample
1895 y = arma_generate_sample([1., -.75], [1., .25], nsample=1000)
1896 arma = ARMA(y)
1897 res = arma.fit(trend='nc', order=(1, 1))
1899 np.random.seed(12345)
1900 y_arma22 = arma_generate_sample([1., -.85, .35], [1, .25, -.9],
1901 nsample=1000)
1902 arma22 = ARMA(y_arma22)
1903 res22 = arma22.fit(trend='nc', order=(2, 2))
1905 # test CSS
1906 arma22_css = ARMA(y_arma22)
1907 res22css = arma22_css.fit(trend='nc', order=(2, 2), method='css')
1909 data = sm.datasets.sunspots.load()
1910 ar = ARMA(data.endog)
1911 resar = ar.fit(trend='nc', order=(9, 0))
1913 y_arma31 = arma_generate_sample([1, -.75, -.35, .25], [.1],
1914 nsample=1000)
1916 arma31css = ARMA(y_arma31)
1917 res31css = arma31css.fit(order=(3, 1), method="css", trend="nc",
1918 transparams=True)
1920 y_arma13 = arma_generate_sample([1., -.75], [1, .25, -.5, .8],
1921 nsample=1000)
1922 arma13css = ARMA(y_arma13)
1923 res13css = arma13css.fit(order=(1, 3), method='css', trend='nc')
1925 # check css for p < q and q < p
1926 y_arma41 = arma_generate_sample([1., -.75, .35, .25, -.3], [1, -.35],
1927 nsample=1000)
1928 arma41css = ARMA(y_arma41)
1929 res41css = arma41css.fit(order=(4, 1), trend='nc', method='css')
1931 y_arma14 = arma_generate_sample([1, -.25], [1., -.75, .35, .25, -.3],
1932 nsample=1000)
1933 arma14css = ARMA(y_arma14)
1934 res14css = arma14css.fit(order=(4, 1), trend='nc', method='css')
1936 # ARIMA Model
1937 from statsmodels.datasets import webuse
1938 dta = webuse('wpi1')
1939 wpi = dta['wpi']
1941 mod = ARIMA(wpi, (1, 1, 1)).fit()
推荐文章
|
|
坚强的伤痕 · 黎明之后（烟若雨创作的系列小说作品）_百度百科 1 年前 |
|
|
眼睛小的野马 · | Hitomi.la 1 年前 |
|
|
乖乖的生姜 · 讲座回顾 | 民航业的发展趋势与中国航空法律制度的完善 1 年前 |