AR(I)MA时间序列建模过程——步骤和python代码

发表: 2017-05-27 浏览: 4320

建模分析数据分析 Python

也算是我自己的笔记，都2017年了，ARMA模型感觉已经过时了吧，用深度学习来预测的表现是不是更好，不过反正都研究了一下了，就做个记录吧。

1.异常值和缺失值的处理

这绝对是数据分析时让所有人都头疼的问题。异常和缺失值会破坏数据的分布，并且干扰分析的结果，怎么处理它们是一门大学问，而我根本还没入门。

（1）异常值

3 ways to remove outliers from your data提供了关于如何对时间序列数据进行异常值检测的方法，作者认为移动中位数的方法最好，代码如下：

from pandas import rolling_median

threshold = 3 #指的是判定一个点为异常的阈值

df['pandas'] = rolling_median(df['u'], window=3, center=True).fillna(method='bfill').fillna(method='ffill') 

#df['u']是原始数据，df['pandas'] 是求移动中位数后的结果，window指的是移动平均的窗口宽度

difference = np.abs(df['u'] - df['pandas'])

outlier_idx = difference > threshold

rolling_median函数详细说明参见pandas.rolling_median

(2)缺失值

缺失值在DataFrame中显示为nan，它会导致ARMA无法拟合，因此一定要进行处理。
a.用序列的均值代替，这样的好处是在计算方差时候不会受影响。但是连续几个nan即使这样替代也会在差分时候重新变成nan，从而影响拟合回归模型。
b.直接删除。我在很多案例上看到这样的做法，但是当一个序列中间的nan太多时，我无法确定这样的做法是否还合理。

2.平稳性检验

序列平稳性是进行时间序列分析的前提条件，主要是运用ADF检验。

from statsmodels.tsa.stattools import adfuller

def test_stationarity(timeseries):

    dftest = adfuller(timeseries, autolag='AIC')

    return dftest[1]

    #此函数返回的是p值

adfuller函数详细说明参见statsmodels.tsa.stattools.adfuller

3.不平稳的处理

（1）对数处理。对数处理可以减小数据的波动，因此无论第1步检验出序列是否平稳，都最好取一次对数。关于为什么统计、计量学家都喜欢对数的原因，知乎上也有讨论：在统计学中为什么要对变量取对数？
（2）差分。一般来说，非纯随机的时间序列经一阶差分或者二阶差分之后就会变得平稳。那差分几阶合理呢？我的观点是：在保证ADF检验的p<0.01的情况下，阶数越小越好，否则会带来样本减少、还原序列麻烦、预测困难的问题。——这是我的直觉，还没有查阅资料求证。基于这样的想法，构造了选择差分阶数的函数：

def best_diff(df, maxdiff = 8):

    p_set = {}

    for i in range(0, maxdiff):

        temp = df.copy() #每次循环前，重置

        if i == 0:

            temp['diff'] = temp[temp.columns[1]]

        else:

            temp['diff'] = temp[temp.columns[1]].diff(i)

            temp = temp.drop(temp.iloc[:i].index) #差分后，前几行的数据会变成nan，所以删掉

        pvalue = test_stationarity(temp['diff'])

        p_set[i] = pvalue

        p_df = pd.DataFrame.from_dict(p_set, orient="index")

        p_df.columns = ['p_value']

    i = 0

    while i < len(p_df):

        if p_df['p_value'][i]<0.01:

            bestdiff = i

            break

        i += 1

    return bestdiff

（3）平滑法。利用移动平均的方法来处理数据，可能可以用来处理周期性因素，我还没实践过。
（4）分解法。将时间序列分解成长期趋势、季节趋势和随机成分，同样没实践过。
对于（3）（4），参见《python时间序列分析》或者Complete guide to create a Time Series Forecast (with Codes in Python)【翻译版《时间序列预测全攻略（附带Python代码）》】

4.随机性检验

只有时间序列不是一个白噪声（纯随机序列）的时候，该序列才可做分析。（参考自：《时间序列ARIMA模型详解：python实现店铺一周销售量预测》）

from statsmodels.stats.diagnostic import acorr_ljungbox

def test_stochastic(ts):

    p_value = acorr_ljungbox(ts, lags=1)[1] #lags可自定义

    return p_value

acorr_ljungbox函数详细说明参见statsmodels.stats.diagnostic.acorr_ljungbox

5.确定ARMA的阶数

ARMA(p,q)是AR(p)和MA(q)模型的组合，关于p和q的选择，一种方法是观察自相关图ACF和偏相关图PACF, 另一种方法是通过借助AIC、BIC统计量自动确定。由于我有几千个时间序列需要分别预测，所以选取自动的方式，而BIC可以有效应对模型的过拟合，因而选定BIC作为判断标准。

from statsmodels.tsa.arima_model import ARMA

def proper_model(data_ts, maxLag): 

    init_bic = float("inf")

    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

这个函数的原理是，根据设定的maxLag，通过循环输入p和q值，选出拟合后BIC最小的p、q值。
然而在statsmodels包里还有更直接的函数：

import statsmodels.tsa.stattools as st

order = st.arma_order_select_ic(timeseries,max_ar=5,max_ma=5,ic=['aic', 'bic', 'hqic'])

order.bic_min_order

timeseries是待输入的时间序列，是pandas.Series类型，max_ar、max_ma是p、q值的最大备选值。
order.bic_min_order返回以BIC准则确定的阶数，是一个tuple类型

6.拟合ARAM

from statsmodels.tsa.arima_model import ARMA

model = ARMA(timeseries, order=order.bic_min_order)

result_arma = model.fit(disp=-1, method='css')

ARMA函数详细说明参见statsmodels.tsa.arima_model.ARMA，.fit参见statsmodels.tsa.arima_model.ARMA.fit
对于差分后的时间序列，运用于ARMA时该模型就被称为ARMIA，在代码层面改写为model = ARIMA(timeseries, order=(p,d,q))，但是实际上，用差分过的序列直接进行ARMA建模更方便，之后添加一步还原的操作即可。

7.预测的y值还原

从前可知，放入模型进行拟合的数据是经过对数或（和）差分处理的数据，因而拟合得到的预测y值要经过差分和对数还原才可与原观测值比较。
暂时写了对数处理过的还原：

def predict_recover(ts):

    ts = np.exp(ts)

    return ts

8.判定拟合优度

在我学习计量经济学的时候，判断一个模型拟合效果是用一个调整R方的指标，但是似乎在机器学习领域，回归时常用RMSE(Root Mean Squared Error，均方根误差)，可能是因为调整R方衡量的预测值与均值之间的差距，而RMSE衡量的是每个预测值与实际值的差距。《均方根值（RMS）+ 均方根误差（RMSE）+标准差（Standard Deviation）》、《（转）SSE,MSE,RMSE,R-square指标讲解》提供了详细公式。

train_predict = result_arma.predict()

train_predict = predict_recover(train_predict) #还原

RMSE = np.sqrt(((train_predict-timeseries)**2).sum()/timeseries.size)

9.预测未来的值

用statsmodel这个包来进行预测，很奇怪的是我从来没成功过，只能进行下一步（之后一天）的预测，多天的就无法做到了。这里提供他们的代码，供大家自行试验，成功的话一定要留言教我啊！跪谢。

predict_ts = result_arma.predict()

predict方法详细说明参见statsmodels.tsa.arima_model.ARMAResults.predict，反正我不太懂这个方法怎么使用……

还有根据How to Create an ARIMA Model for Time Series Forecasting with Python，用来预测的代码是：

for t in range(len(test)):

    model = ARIMA(history, order=(5,1,0))

    model_fit = model.fit(disp=0)

    output = model_fit.forecast()

    yhat = output[0]

    predictions.append(yhat)

    obs = test[t]

    history.append(obs)

    print('predicted=%f, expected=%f' % (yhat, obs))

然而我真的不懂，按他写forecast方法的方式，每次循环预测的都是history样本的下一个值，因而如何用这个循环来预测history样本的之后，比如十个值？如果不用循环，直接令forecast中的参数steps=为要预测的时长，我也没成功……
forecast方法详细说明参见statsmodels.tsa.arima_model.ARIMAResults.forecast，
此外，Stackoverflow上的一个解答：ARMA out-of-sample prediction with statsmodels，又给了一个预测的写法。

10. 更方便的时间序列包：`pyflux`

好在《AR、MA及ARMA模型》提到了python的另一个包pyflux，它的文档在PyFlux 0.4.0 documentation。这个包在macOS上安装之前需要安装XCode命令行工具：

xcode-select --install

同时它的画图需要安装一个seaborn的包（如果没有Anaconda则用pip的方式。《数据可视化（三）- Seaborn简易入门》简要介绍了seaborn，它是“在matplotlib的基础上进行了更高级的API封装”。

conda install seaborn

我用这个包写了一个简略而完整的ARIMA建模：

# -*- coding: utf-8 -*-

"""

Created on Tue Jan 31 14:13:58 2017

@author: 竹间为简

@published in: 简书

"""



import pandas as pd

from statsmodels.tsa.stattools import adfuller

import statsmodels.tsa.stattools as st

import numpy as np

import pyflux as pf





daily_payment = pd.read_csv('xxx.csv',parse_dates=[0], index_col=0)



def test_stationarity(timeseries):

    dftest = adfuller(timeseries, autolag='AIC')

    return dftest[1]





def best_diff(df, maxdiff = 8):

    p_set = {}

    for i in range(0, maxdiff):

        temp = df.copy() #每次循环前，重置

        if i == 0:

            temp['diff'] = temp[temp.columns[1]]

        else:

            temp['diff'] = temp[temp.columns[1]].diff(i)

            temp = temp.drop(temp.iloc[:i].index) #差分后，前几行的数据会变成nan，所以删掉

        pvalue = test_stationarity(temp['diff'])

        p_set[i] = pvalue

        p_df = pd.DataFrame.from_dict(p_set, orient="index")

        p_df.columns = ['p_value']

    i = 0

    while i < len(p_df):

        if p_df['p_value'][i]<0.01:

            bestdiff = i

            break

        i += 1

    return bestdiff





def produce_diffed_timeseries(df, diffn):

    if diffn != 0:

        df['diff'] = df[df.columns[1]].apply(lambda x:float(x)).diff(diffn)

    else:

        df['diff'] = df[df.columns[1]].apply(lambda x:float(x))

    df.dropna(inplace=True) #差分之后的nan去掉

    return df





def choose_order(ts, maxar, maxma):

    order = st.arma_order_select_ic(ts, maxar, maxma, ic=['aic', 'bic', 'hqic'])

    return order.bic_min_order





def predict_recover(ts, df, diffn):

    if diffn != 0:

        ts.iloc[0] = ts.iloc[0]+df['log'][-diffn]

        ts = ts.cumsum()

    ts = np.exp(ts)

#    ts.dropna(inplace=True)

    print('还原完成')

    return ts





def run_aram(df, maxar, maxma, test_size = 14):

    data = df.dropna()

    data['log'] = np.log(data[data.columns[0]])

    #    test_size = int(len(data) * 0.33)

    train_size = len(data)-int(test_size)

    train, test = data[:train_size], data[train_size:]

    if test_stationarity(train[train.columns[1]]) < 0.01:

        print('平稳，不需要差分')

    else:

        diffn = best_diff(train, maxdiff = 8)

        train = produce_diffed_timeseries(train, diffn)

        print('差分阶数为'+str(diffn)+'，已完成差分')

    print('开始进行ARMA拟合')

    order = choose_order(train[train.columns[2]], maxar, maxma)

    print('模型的阶数为：'+str(order))

    _ar = order[0]

    _ma = order[1]

    model = pf.ARIMA(data=train, ar=_ar, ma=_ma, target='diff', family=pf.Normal())

    model.fit("MLE")

    test = test['payment_times']

    test_predict = model.predict(int(test_size))

    test_predict = predict_recover(test_predict, train, diffn)

    RMSE = np.sqrt(((np.array(test_predict)-np.array(test))**2).sum()/test.size)

    print("测试集的RMSE为："+str(RMSE))

pyflux的predict函数就十分易用，model.predict(h = )就可。详细参见ARIMA的文档，画图起来也是十分方便。
Time Series Forecasting using ARIMA in Python也提供了利用pyflux进行建模的例子。

11. 调参

对于ARIMA来说，可调参的地方也挺多：
差分阶数。
p、q的阶数。
模型拟合的方法：MLE、OLS等，参见Bayesian Inference和Classical Inference。
预测的周期、滚动预测的周期。

12. 结合卡尔曼滤波

在时间序列中存在的噪声，会干扰序列中每个点的波动，给预测造成难度，所以人们就想出一个办法来过滤这个噪声，一个有名的办法叫”卡尔曼滤波“
这个东西我还没研究……给出参考资料：

其他参考读物：

在针对ARIMA以及时间序列分析话题搜寻资料时，还接触了以下资料：

《时间序列模型的三个重要概念》：自相关、白噪声、平稳性。
《时间序列分析之一：数据预处理》，讲了对数据进行平稳、正态、独立、周期、趋势项检验，以及提了一句对异常值的处理。
有什么好的模型可以做高精度的时间序列预测呢？提供了一些关于预测模型选择的见解。
时间序列建模问题，如何准确的建立时间序列模型？讨论了建立时间序列模型的流程。
Machine Learning Strategies for Time Series Prediction
Data transformations and forecasting models: what to use and when，顾名思义，列举了预测模型的适用范围。
Autoregressive moving average model Wikipedia
《需求预测我喜欢用ARIMA模型》，通俗、形象。
《金融时间序列入门（二）----MA & ARMA & ARIMA》
《python时间序列分析之ARIMA》的贡献在于讲到了对模型的检验，然而没说清楚做检验的用处，真要了解还得去看计量经济学书籍。
Seasonal ARIMA with Python，贡献在于讲到了seasonal difference的方法，但也没讲如何确定一个序列应该应用这种差分方法。
How to Use and Remove Trend Information from Time Series Data in Python，如果一个时间序列带有趋势的特征，如何应对的方法。
How to Make Predictions for Time Series Forecasting with Python，贡献在于讲了如何保存拟合好的模型，以及有新的数据加入的时候，更新预测模型。
Jason Brownlee这个人最近写了N篇关于TimeSeries的文章……

题外话：

<Evaluating Machine Learning Models>讲解了对各种模型使用的评价指标、调参的方法以及AB测试的陷阱，它的未完成翻译版见《机器学习模型评价(Evaluating Machine Learning Models)-主要概念与陷阱》

0 个评论

要回复文章请先登录或注册