文档章节

AR(I)MA时间序列建模过程--with python

易野
 易野
发布于 2017/07/12 15:53
字数 3155
阅读 48
收藏 0
点赞 0
评论 0

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这个包来进行预测,很奇怪的是我从来没成功过,只能进行下一步(之后一天)的预测,多天的就无法做到了。这里提供他们的代码,供大家自行试验,成功的话一定要留言教我啊!跪谢。

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)-主要概念与陷阱

结语:

信息过载是当今时代的毛病,搜一个主题可以得到浩如烟海的资料,以及浩如烟海的相关主题资料……他们之间还是互相抄的,所以真的累死了。

本文转载自:www.jianshu.com/p

共有 人打赏支持
易野
粉丝 3
博文 142
码字总数 114058
作品 0
深圳
用ARIMA模型做需求预测

本文结构: 时间序列分析? 什么是ARIMA? ARIMA数学模型? input,output 是什么? 怎么用?-代码实例 常见问题? 时间序列分析? 时间序列,就是按时间顺序排列的,随时间变化的数据序列。...

aliceyangxi1987 ⋅ 2017/05/02 ⋅ 0

【python数据挖掘课程】二十三.时间序列金融数据预测及Pandas库详解

这是《Python数据挖掘课程》系列文章,也是我上课内容及书籍中的一个案例。本文主要讲述时间序列算法原理,Pandas扩展包基本用法以及Python调用statsmodels库的时间序列算法。由于作者数学比...

eastmount ⋅ 05/09 ⋅ 0

手把手教你用Prophet快速进行时间序列预测(附Prophet和R代码)

简介 对于任何业务而言,基于时间进行分析都是至关重要的。库存量应该保持在多少?你希望商店的客流量是多少?多少人会乘坐飞机旅游?类似这样待解决的问题都是重要的时间序列问题。 这就是时...

技术小能手 ⋅ 06/12 ⋅ 0

R语言 时间序列ARIMA模型方法

原理什么的百度一搜一堆,看不明白,先学会用这个工具吧!   ARIMA:全称为自回归积分滑动平均模型(Autoregressive Integrated Moving Average Model,简记ARIMA),是由博克思(Box)和詹金斯...

gdyflxw ⋅ 2017/02/17 ⋅ 0

时序预测(网络流量预测)方法调研总结

主要分为线性时间序列预测模型、非线性时间序列预测模型、神经网络时间序列预测模型、Boosting预测模型、GM预测模型等。 线性时间序列模型 2 (一)自回归模型(AR(p)) 2 (二)滑动平均模型(...

qq_36558948 ⋅ 03/16 ⋅ 0

[DataAnalysis]时间序列分析

一、平稳性 1、严平稳与宽平稳的定义,一般我们都用二阶宽平稳 2、为什么要研究平稳性:若对非平稳时间序列使用现有的方法估计,则会得到虚假回归,估计模型无效。 3、ADF与DF统计量检验时间...

tomocat ⋅ 01/18 ⋅ 0

Python发展迅速,成为学术界新主流

如果说2018年以前R是数据学术界的主流,但是现在Python正在慢慢取代R在学术界的地位。 Python与R相比速度要快。Python可以直接处理上G的数据;R不行,R分析数据时需要先通过数据库把大数据转...

Python燕大侠 ⋅ 05/07 ⋅ 0

Python 用于金融数据分析第9课 Part 2-----时间序列分析模型ARIMA实现

ARIMA模型实现预测的流程 流程图 一、导入数据 我们先在终端看一下牛奶产量这个文件的前几行再决定要怎么样读取这个csv文件。 文件概览 我们可以看到有两列,第一列是时间,第二列是产量,因...

九日照林 ⋅ 03/23 ⋅ 0

7. Python3源码—Dict对象

7.1. 散列表 散列表的基本思想,是通过一定的函数将需搜索的键值映射为一个整数,将这个整数视为索引值去访问某片连续的内存区域。理论上,在最优情况下,散列表能提供O(1)复杂度的搜索效率。...

whj0709 ⋅ 06/06 ⋅ 0

团队拙作《Python机器学习实战》

之前看国内外的 Python 机器学习的书,鲜有将机器学习到底怎么做人脸识别、怎么做风险控制、怎么做 OCR 算法模型列出的,并且真正的一个 Python 应用,不止是从机器学习库中导入一下配置一下...

yijun2018 ⋅ 04/20 ⋅ 0

没有更多内容

加载失败,请刷新页面

加载更多

下一页

tcp/ip详解-链路层

简介 设计链路层的目的: 为IP模块发送和接收IP数据报 为ARP模块发送ARP请求和接收ARP应答 为RARP模块发送RARP请求和接收RARP应答 TCP/IP支持多种链路层协议,如以太网、令牌环往、FDDI、RS-...

loda0128 ⋅ 54分钟前 ⋅ 0

spring.net aop代码例子

https://www.cnblogs.com/haogj/archive/2011/10/12/2207916.html

whoisliang ⋅ 今天 ⋅ 0

发送短信如何限制1小时内最多发送11条短信

发送短信如何限制1小时内最多发送11条短信 场景: 发送短信属于付费业务,有时为了防止短信攻击,需要限制发送短信的频率,例如在1个小时之内最多发送11条短信. 如何实现呢? 思路有两个 截至到当...

黄威 ⋅ 昨天 ⋅ 0

mysql5.7系列修改root默认密码

操作系统为centos7 64 1、修改 /etc/my.cnf,在 [mysqld] 小节下添加一行:skip-grant-tables=1 这一行配置让 mysqld 启动时不对密码进行验证 2、重启 mysqld 服务:systemctl restart mysql...

sskill ⋅ 昨天 ⋅ 0

Intellij IDEA神器常用技巧六-Debug详解

在调试代码的时候,你的项目得debug模式启动,也就是点那个绿色的甲虫启动服务器,然后,就可以在代码里面断点调试啦。下面不要在意,这个快捷键具体是啥,因为,这个keymap是可以自己配置的...

Mkeeper ⋅ 昨天 ⋅ 0

zip压缩工具、tar打包、打包并压缩

zip 支持压缩目录 1.在/tmp/目录下创建目录(study_zip)及文件 root@yolks1 study_zip]# !treetree 11└── 2 └── 3 └── test_zip.txt2 directories, 1 file 2.yum...

蛋黄Yolks ⋅ 昨天 ⋅ 0

聊聊HystrixThreadPool

序 本文主要研究一下HystrixThreadPool HystrixThreadPool hystrix-core-1.5.12-sources.jar!/com/netflix/hystrix/HystrixThreadPool.java /** * ThreadPool used to executed {@link Hys......

go4it ⋅ 昨天 ⋅ 0

容器之上传镜像到Docker hub

Docker hub在国内可以访问,首先要创建一个账号,这个后面会用到,我是用126邮箱注册的。 1. docker login List-1 Username不能使用你注册的邮箱,要用使用注册时用的username;要输入密码 ...

汉斯-冯-拉特 ⋅ 昨天 ⋅ 0

SpringBoot简单使用ehcache

1,SpringBoot版本 2.0.3.RELEASE ①,pom.xml <parent><groupId>org.springframework.boot</groupId><artifactId>spring-boot-starter-parent</artifactId><version>2.0.3.RELE......

暗中观察 ⋅ 昨天 ⋅ 0

Spring源码解析(八)——实例创建(下)

前言 来到实例创建的最后一节,前面已经将一个实例通过不同方式(工厂方法、构造器注入、默认构造器)给创建出来了,下面我们要对创建出来的实例进行一些“加工”处理。 源码解读 回顾下之前...

MarvelCode ⋅ 昨天 ⋅ 0

没有更多内容

加载失败,请刷新页面

加载更多

下一页

返回顶部
顶部