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

转载自:https://www.jianshu.com/p/cced6617b423

侵删

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_armax_mapq值的最大备选值。
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这个包来进行预测,很奇怪的是我从来没成功过,只能进行下一步(之后一天)的预测,多天的就无法做到了。可以啊,根据How to Create an ARIMA Model for Time Series Forecasting with Python,输入前100天的数据,则可以预测任意时间长度的值,只需要

predict_ts = result_arma.predict(101,110,...)
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))

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

11. 调参

  • 对于ARIMA来说,可调参的地方也挺多:
  • 差分阶数。
  • pq的阶数。
  • 模型拟合的方法:MLE、OLS等,参见Bayesian InferenceClassical Inference
  • 预测的周期、滚动预测的周期。

12. 结合卡尔曼滤波

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

其他参考读物:

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

题外话:

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



作者:竹间为简
链接:https://www.jianshu.com/p/cced6617b423
來源:简书
著作权归作者所有。商业转载请联系作者获得授权,非商业转载请注明出处。

阅读更多
想对作者说点什么?

博主推荐

换一批

没有更多推荐了,返回首页