python时间序列分析（ARIMA模型）

原文地址：https://blog.csdn.net/u011596455/article/details/78650458

转载请注明出处。

什么是时间序列

时间序列简单的说就是各时间点上形成的数值序列，时间序列分析就是通过观察历史数据预测未来的值。在这里需要强调一点的是，时间序列分析并不是关于时间的回归，它主要是研究自身的变化规律的（这里不考虑含外生变量的时间序列）。

为什么用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
  9  
 10 class arima_model:
 11  
 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
 21  
 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)
 27  
 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()
 47  
 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'
 61  
 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
 73  
 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)
 87  
 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
 91  
 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
 99  
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)
105  
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']
110  
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)
119  
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')
127  
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)
138  
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
   6  
   7 from __future__ import absolute_import
   8 from statsmodels.compat.python import string_types, range
   9 # for 2to3 with extensions
  10  
  11 from datetime import datetime
  12  
  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
  19  
  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
  35  
  36 _armax_notes = """
  37  
  38         Notes
  39         -----
  40         If exogenous variables are given, then the model that is fit is
  41  
  42         .. math::
  43  
  44            \\phi(L)(y_t - X_t\\beta) = \\theta(L)\epsilon_t
  45  
  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 """
  55  
  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."""
  65  
  66 _arma_model = "Autoregressive Moving Average ARMA(p,q) Model"
  67  
  68 _arima_model = "Autoregressive Integrated Moving Average ARIMA(p,d,q) Model"
  69  
  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."""
  79  
  80 _predict_notes = """
  81         Notes
  82         -----
  83         Use the results predict method instead.
  84 """
  85  
  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::
  94  
  95          1970Q1, 1
  96          1970Q2, 1.5
  97          1970Q3, 1.25
  98          1970Q4, 2.25
  99          1971Q1, 1.2
 100          1971Q2, 4.1
 101  
 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 """
 108  
 109 _predict = """
 110         %(Model)s model in-sample and out-of-sample prediction
 111  
 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
 138  
 139         Returns
 140         -------
 141         %(returns)s
 142         %(extra_section)s
 143 """
 144  
 145 _predict_returns = """predict : array
 146             The predicted values.
 147  
 148 """
 149  
 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}
 157  
 158 _arma_results_predict = _predict % {"Model" : "ARMA", "params" : "",
 159                                     "extra_params" : "",
 160                                     "returns" : _predict_returns,
 161                                     "extra_section" : _results_notes}
 162  
 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'}
 167  
 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}
 173  
 174 _arima_results_predict = _predict % {"Model" : "ARIMA",
 175                                      "params" : "",
 176                                      "extra_params" :
 177                                      """typ : str {'linear', 'levels'}
 178  
 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}
 185  
 186 _arima_plot_predict_example = """        Examples
 187         --------
 188         >>> import statsmodels.api as sm
 189         >>> import matplotlib.pyplot as plt
 190         >>> import pandas as pd
 191         >>>
 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()
 200  
 201         .. plot:: plots/arma_predict_plot.py
 202 """
 203  
 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                       }
 219  
 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                       }
 240  
 241  
 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
 249  
 250  
 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)
 258  
 259  
 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
 277  
 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)
 307  
 308     endog = np.zeros(p + steps - 1)
 309  
 310     if p and start:
 311         endog[:p] = y[start-p:start]
 312     elif p:
 313         endog[:p] = y[-p:]
 314  
 315     return endog, resid, mu
 316  
 317  
 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]
 338  
 339     if q:
 340         i = 0  # if q == 1
 341     else:
 342         i = -1
 343  
 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
 349  
 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
 354  
 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
 358  
 359  
 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
 368  
 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]
 374  
 375  
 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)
 383  
 384     return start
 385  
 386  
 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
 397  
 398  
 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
 403  
 404  
 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
 415  
 416  
 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
 431  
 432  
 433 def _check_estimable(nobs, n_params):
 434     if nobs <= n_params:
 435         raise ValueError("Insufficient degrees of freedom to estimate")
 436  
 437  
 438 class ARMA(tsbase.TimeSeriesModel):
 439  
 440     __doc__ = tsbase._tsa_doc % {"model" : _arma_model,
 441                                  "params" : _arma_params, "extra_params" : "",
 442                                  "extra_sections" : _armax_notes %
 443                                  {"Model" : "ARMA"}}
 444  
 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
 461  
 462     def _fit_start_params_hr(self, order):
 463         """
 464         Get starting parameters for fit.
 465  
 466         Parameters
 467         ----------
 468         order : iterable
 469             (p,q,k) - AR lags, MA lags, and number of exogenous variables
 470             including the constant.
 471  
 472         Returns
 473         -------
 474         start_params : array
 475             A first guess at the starting parameters.
 476  
 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.
 483  
 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
 534  
 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.")
 550  
 551         # check MA coefficients
 552         # print start_params
 553         return start_params
 554  
 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
 571  
 572     def score(self, params):
 573         """
 574         Compute the score function at params.
 575  
 576         Notes
 577         -----
 578         This is a numerical approximation.
 579         """
 580         return approx_fprime_cs(params, self.loglike, args=(False,))
 581  
 582     def hessian(self, params):
 583         """
 584         Compute the Hessian at params,
 585  
 586         Notes
 587         -----
 588         This is a numerical approximation.
 589         """
 590         return approx_hess_cs(params, self.loglike, args=(False,))
 591  
 592     def _transparams(self, params):
 593         """
 594         Transforms params to induce stationarity/invertability.
 595  
 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)
 603  
 604         # just copy exogenous parameters
 605         if k != 0:
 606             newparams[:k] = params[:k]
 607  
 608         # AR Coeffs
 609         if k_ar != 0:
 610             newparams[k:k+k_ar] = _ar_transparams(params[k:k+k_ar].copy())
 611  
 612         # MA Coeffs
 613         if k_ma != 0:
 614             newparams[k+k_ar:] = _ma_transparams(params[k+k_ar:].copy())
 615         return newparams
 616  
 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)
 629  
 630         # MA coeffs
 631         if k_ma != 0:
 632             newparams[k+k_ar:k+k_ar+k_ma] = _ma_invtransparams(macoefs)
 633         return newparams
 634  
 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
 655  
 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)
 659  
 660     def geterrors(self, params):
 661         """
 662         Get the errors of the ARMA process.
 663  
 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         """
 672  
 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
 678  
 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)
 684  
 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])
 695  
 696             k_ar = self.k_ar
 697             k_ma = self.k_ma
 698  
 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()
 710  
 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)
 714  
 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")
 720  
 721         endog = self.endog
 722         resid = self.geterrors(params)
 723         k_ar = self.k_ar
 724  
 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))
 732  
 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
 744  
 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
 758  
 759     def loglike(self, params, set_sigma2=True):
 760         """
 761         Compute the log-likelihood for ARMA(p,q) model
 762  
 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)
 774  
 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)
 780  
 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:]
 803  
 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
 810  
 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.
 816  
 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.
 863  
 864         Returns
 865         -------
 866         statsmodels.tsa.arima_model.ARMAResults class
 867  
 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
 873  
 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')
 880  
 881         """
 882         k_ar = self.k_ar
 883         k_ma = self.k_ma
 884  
 885         # enforce invertibility
 886         self.transparams = transparams
 887  
 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'
 891  
 892         # (re)set trend and handle exogenous variables
 893         # always pass original exog
 894         k_trend, exog = _make_arma_exog(endog, self.exog, trend)
 895  
 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.")
 901  
 902         # check again now that we know the trend
 903         _check_estimable(len(endog), k_ar + k_ma + k_exog + k_trend)
 904  
 905         self.k_trend = k_trend
 906         self.exog = exog    # overwrites original exog from __init__
 907  
 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
 912  
 913         # choose objective function
 914         if k_ma == 0 and k_ar == 0:
 915             method = "css"  # Always CSS when no AR or MA terms
 916  
 917         self.method = method = method.lower()
 918  
 919         # adjust nobs for css
 920         if method == 'css':
 921             self.nobs = len(self.endog) - k_ar
 922  
 923         if start_params is not None:
 924             start_params = np.asarray(start_params)
 925  
 926         else:  # estimate starting parameters
 927             start_params = self._fit_start_params((k_ar, k_ma, k), method)
 928  
 929         if transparams:  # transform initial parameters to ensure invertibility
 930             start_params = self._invtransparams(start_params)
 931  
 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
 942  
 943         if transparams:  # transform parameters back
 944             params = self._transparams(params)
 945  
 946         self.transparams = False  # so methods don't expect transf.
 947  
 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)
 954  
 955  
 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):
 960  
 961     __doc__ = tsbase._tsa_doc % {"model" : _arima_model,
 962                                  "params" : _arima_params, "extra_params" : "",
 963                                  "extra_sections" : _armax_notes %
 964                                  {"Model" : "ARIMA"}}
 965  
 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
 975  
 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?
 997  
 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
1026  
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
1035  
1036         return end - self.k_diff, out_of_sample
1037  
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.
1043  
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.
1090  
1091         Returns
1092         -------
1093         `statsmodels.tsa.arima.ARIMAResults` class
1094  
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
1100  
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')
1107  
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)
1118  
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)
1144  
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
1196  
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
1209  
1210         else:  # pragma : no cover
1211             raise ValueError("typ %s not understood" % typ)
1212  
1213     predict.__doc__ = _arima_predict
1214  
1215  
1216 class ARMAResults(tsbase.TimeSeriesModelResults):
1217     """
1218     Class to hold results from fitting an ARMA model.
1219  
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.
1230  
1231     Returns
1232     --------
1233     **Attributes**
1234  
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 = {}
1319  
1320     #TODO: use this for docstring when we fix nobs issue
1321  
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
1343  
1344     @cache_readonly
1345     def arroots(self):
1346         return np.roots(np.r_[1, -self.arparams])**-1
1347  
1348     @cache_readonly
1349     def maroots(self):
1350         return np.roots(np.r_[1, self.maparams])**-1
1351  
1352     @cache_readonly
1353     def arfreq(self):
1354         r"""
1355         Returns the frequency of the AR roots.
1356  
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)
1364  
1365     @cache_readonly
1366     def mafreq(self):
1367         r"""
1368         Returns the frequency of the MA roots.
1369  
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)
1377  
1378     @cache_readonly
1379     def arparams(self):
1380         k = self.k_exog + self.k_trend
1381         return self.params[k:k+self.k_ar]
1382  
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:]
1388  
1389     @cache_readonly
1390     def llf(self):
1391         return self.model.loglike(self.params)
1392  
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)))
1400  
1401     def cov_params(self):  # add scale argument?
1402         params = self.params
1403         hess = self.model.hessian(params)
1404         return -inv(hess)
1405  
1406     @cache_readonly
1407     def aic(self):
1408         return -2 * self.llf + 2 * self._ic_df_model
1409  
1410     @cache_readonly
1411     def bic(self):
1412         nobs = self.nobs
1413         return -2 * self.llf + np.log(nobs) * self._ic_df_model
1414  
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
1419  
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
1438  
1439     @cache_readonly
1440     def resid(self):
1441         return self.model.geterrors(self.params)
1442  
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
1448  
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
1452  
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)
1457  
1458         fcasterr = np.sqrt(sigma2 * np.cumsum(ma_rep**2))
1459         return fcasterr
1460  
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]
1465  
1466         return conf_int
1467  
1468     def forecast(self, steps=1, exog=None, alpha=.05):
1469         """
1470         Out-of-sample forecasts
1471  
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) %
1483  
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))
1508  
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
1515  
1516         # compute the standard errors
1517         fcasterr = self._forecast_error(steps)
1518         conf_int = self._forecast_conf_int(forecast, fcasterr, alpha)
1519  
1520         return forecast, fcasterr, conf_int
1521  
1522     def summary(self, alpha=.05):
1523         """Summarize the Model
1524  
1525         Parameters
1526         ----------
1527         alpha : float, optional
1528             Significance level for the confidence intervals.
1529  
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.
1535  
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))
1556  
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                     ]
1570  
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])]
1578  
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)
1583  
1584         # Make the roots table
1585         from statsmodels.iolib.table import SimpleTable
1586  
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"])
1617  
1618             smry.tables.append(roots_table)
1619         return smry
1620  
1621     def summary2(self, title=None, alpha=.05, float_format="%.4f"):
1622         """Experimental summary function for ARIMA Results
1623  
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
1633  
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.
1639  
1640         See Also
1641         --------
1642         statsmodels.iolib.summary2.Summary : class to hold summary
1643             results
1644  
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))
1659  
1660         k_ar, k_ma = self.k_ar, self.k_ma
1661  
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 = []
1681  
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
1688  
1689         # Summary
1690         from statsmodels.iolib import summary2
1691         smry = summary2.Summary()
1692  
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))
1701  
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)
1709  
1710         return smry
1711  
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)
1717  
1718  
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)
1724  
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)
1730  
1731  
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)
1738  
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)
1744  
1745         if plot_insample:
1746             ax.plot(x[:end + 1 - start], self.model.endog[start:end+1],
1747                     label=self.model.endog_names)
1748  
1749         ax.legend(loc='best')
1750  
1751         return fig
1752     plot_predict.__doc__ = _plot_predict
1753  
1754  
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)
1763  
1764  
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
1770  
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)
1775  
1776         fcerr = np.sqrt(np.cumsum(cumsum_n(ma_rep, self.k_diff)**2)*sigma2)
1777         return fcerr
1778  
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
1783  
1784     def forecast(self, steps=1, exog=None, alpha=.05):
1785         """
1786         Out-of-sample forecasts
1787  
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) %
1799  
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
1808  
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)
1827  
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:]
1832  
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
1837  
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)
1843  
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)
1849  
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)
1855  
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)
1862  
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)
1868  
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)
1877  
1878         ax.legend(loc='best')
1879  
1880         return fig
1881  
1882     plot_predict.__doc__ = _arima_plot_predict
1883  
1884  
1885 class ARIMAResultsWrapper(ARMAResultsWrapper):
1886     pass
1887 wrap.populate_wrapper(ARIMAResultsWrapper, ARIMAResults)
1888  
1889  
1890 if __name__ == "__main__":
1891     import statsmodels.api as sm
1892  
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))
1898  
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))
1904  
1905     # test CSS
1906     arma22_css = ARMA(y_arma22)
1907     res22css = arma22_css.fit(trend='nc', order=(2, 2), method='css')
1908  
1909     data = sm.datasets.sunspots.load()
1910     ar = ARMA(data.endog)
1911     resar = ar.fit(trend='nc', order=(9, 0))
1912  
1913     y_arma31 = arma_generate_sample([1, -.75, -.35, .25], [.1],
1914                                     nsample=1000)
1915  
1916     arma31css = ARMA(y_arma31)
1917     res31css = arma31css.fit(order=(3, 1), method="css", trend="nc",
1918                              transparams=True)
1919  
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')
1924  
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')
1930  
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')
1935  
1936     # ARIMA Model
1937     from statsmodels.datasets import webuse
1938     dta = webuse('wpi1')
1939     wpi = dta['wpi']
1940  
1941     mod = ARIMA(wpi, (1, 1, 1)).fit()

posted @ 2018-12-16 16:43 hou永胜阅读(56559) 评论(5) 编辑收藏举报

会员力量，点亮园子希望

刷新页面返回顶部

hou永胜

python时间序列分析（ARIMA模型）

公告