ARIMA模型预测产品销量

任务：做出一个时间序列预测模型，从而实现对一系列产品销量的月预测

你好！ 这是你第一次使用 Markdown编辑器 所展示的欢迎页。如果你想学习如何使用Markdown编辑器, 可以仔细阅读这篇文章，了解一下Markdown的基本语法知识。

方法：针对数据特征，采用了ARIMA模型。并成功实现了任务

在这里预设了8种数据预处理模型，从而使数据平稳，保证尽可能预测更多产品

`

# -*- coding: utf-8 -*-
import os
import math
import pandas as pd
import numpy as np
from statsmodels.graphics.tsaplots import plot_acf, plot_pacf
from statsmodels.tsa.stattools import adfuller
from statsmodels.tsa import arima_model
from statsmodels.stats.diagnostic import acorr_ljungbox
import warnings
warnings.filterwarnings("ignore")


def ARIMA_M0(series,n=3):
    '''
    只讨论ARIMA模型，预测，数字索引从1开始
    series:时间序列,值以月为单位输出，原始数据应该超过160天，不然可能报错
    n:需要往后预测的月个数
    '''
    #数据预处理与清洗
    series.index=pd.to_datetime(series.index)    
    series=series.resample('D').asfreq().fillna(0)
    series[(series/series.mean())>10]=int(4*series.mean())
    series=series.resample('M').sum()
    series[series<=0]=1
    #对清洗后的数据处理
    series = np.array(series)
    series = pd.Series(series.reshape(-1))
    currentDir = os.getcwd()#获得当前工作路径
    #差分数据处理
    ser_sqrt=np.sqrt(series)
    ser_sqrt1=ser_sqrt.diff(1)[1:]
    ser_sqrt2=ser_sqrt1.diff(1)[1:]
    ser_1=series.diff(1)[1:]
    ser_2=ser_1.diff(1)[1:]
    ser_pow=series**2
    ser_pow1=ser_pow.diff(1)[1:]
    ser_pow2=ser_pow1.diff(1)[1:]
    ser_log=np.log(series)
    ser_log1=ser_log.diff(1).dropna()
    ser_log2=ser_log1.diff(1).dropna()
    #创建数据列表
    SERIES={'原数1差':ser_1,'原数2差':ser_2,'原对1差':ser_log1,'原对2差':ser_log2,'原开1差':ser_sqrt1,'原开2差':ser_sqrt2,'原乘1差':ser_pow1,'原乘2差':ser_pow2}
    cont0=[]
    cont=[]
    for i in SERIES:
        #plot_acf(SERIES[i]).savefig(currentDir+'/一阶差分自相关图.png')
        #plot_pacf(SERIES[i]).savefig(currentDir+'/一阶差分偏自相关图.png')
        #单位根检验+白噪声检验
        unitP = adfuller(SERIES[i])
        noiseP = max(acorr_ljungbox(SERIES[i], lags=3)[-1])
        if unitP[1]<0.05 and unitP[0]<min(unitP[4].values()) and noiseP<=0.05:
            unitAssess = '单位根检验中p值为%f，小于0.05，白噪声检验中p值为%f，小于0.05，认为该序列判断为非白噪声平稳序列'%(unitP[1],noiseP)
            cont0.append(i)
        else:
            unitAssess = '单位根检验中p值为%f，白噪声检验中p值为%f，认为该序列不适合做ARIMA时间预测序列'%(unitP[1],noiseP)
            cont.append((i,unitAssess))
            
            #BIC准则确定p、q值
    for i in cont0:
        bics = list()
        for p in range(8):
            tmp = list()
            for q in range(5):
                if i=='原数1差':
                    try:
                        tmp.append(arima_model.ARIMA(series, (p, 1, q)).fit().bic)
                    except Exception as e:
                        #print(str(e))
                        tmp.append(1e+10)#加入一个很大的数
                elif i=='原开1差':
                    try:
                        tmp.append(arima_model.ARIMA(ser_sqrt, (p, 1, q)).fit().bic)
                    except Exception as e:
                        #print(str(e))
                        tmp.append(1e+10)#加入一个很大的数
                elif i=='原开2差':
                    try:
                        tmp.append(arima_model.ARIMA(ser_sqrt, (p, 2, q)).fit().bic)
                    except Exception as e:
                        #print(str(e))
                        tmp.append(1e+10)#加入一个很大的数
                elif i=='原乘1差':
                    try:
                        tmp.append(arima_model.ARIMA(ser_pow, (p, 1, q)).fit().bic)
                    except Exception as e:
                        #print(str(e))
                        tmp.append(1e+10)#加入一个很大的数
                elif i=='原乘2差':
                    try:
                        tmp.append(arima_model.ARIMA(ser_pow, (p, 2, q)).fit().bic)
                    except Exception as e:
                        #print(str(e))
                        tmp.append(1e+10)#加入一个很大的数
                elif i=='原对1差':
                    try:
                        tmp.append(arima_model.ARIMA(ser_log, (p, 1, q)).fit().bic)
                    except Exception as e:
                        #print(str(e))
                        tmp.append(1e+10)#加入一个很大的数
                elif i=='原对2差':
                    try:
                        tmp.append(arima_model.ARIMA(ser_log, (p, 2, q)).fit().bic)
                    except Exception as e:
                        #print(str(e))
                        tmp.append(1e+10)#加入一个很大的数
                else:
                    try:
                        tmp.append(arima_model.ARIMA(series, (p, 2, q)).fit().bic)
                    except Exception as e:
                        #print(str(e))
                        tmp.append(1e+10)#加入一个很大的数        
            bics.append(tmp)
        bics = pd.DataFrame(bics)
        p, q = bics.stack().idxmin()
        #print('BIC准则下确定p,q为%s,%s'%(p,q))
        #建模
        if i=='原数1差':
            model = arima_model.ARIMA(series,order=(p, 1, q)).fit()
        elif i=='原开1差':
            model = arima_model.ARIMA(ser_sqrt,order=(p, 1, q)).fit()
        elif i=='原开2差':   
            model = arima_model.ARIMA(ser_sqrt,order=(p, 2, q)).fit()
        elif i=='原乘1差':
            model = arima_model.ARIMA(ser_pow,order=(p, 1, q)).fit()
        elif i=='原乘2差':   
            model = arima_model.ARIMA(ser_pow,order=(p, 2, q)).fit()
        elif i=='原对1差':
            model = arima_model.ARIMA(ser_log,order=(p, 1, q)).fit()
        elif i=='原对2差':   
            model = arima_model.ARIMA(ser_log,order=(p, 2, q)).fit()
        else:
            model = arima_model.ARIMA(series,order=(p, 2, q)).fit()
        #白噪声残差检测
        noiseR = min(acorr_ljungbox(model.resid, lags=3)[-1])
        if noiseR>0.05:
            noiseAsse = '残差白噪声检验%f，大于0.05，认为该模型残差序列为白噪声.模型可用'%noiseR
            if i=='原开1差' or i=='原开2差':
                predict = (model.forecast(n)[0])**2
            elif i=='原对1差' or i=='原对2差':
                predict = np.exp(model.forecast(n)[0])
            elif i=='原乘1差' or i=='原乘2差':
                pdt_pow=model.forecast(n)[0]
                pdt_pow[pdt_pow<=0]=1
                predict =np.sqrt(pdt_pow)
            else:
                predict = model.forecast(n)[0]
            cont.append((i,'AR模型阶数p',p,'MA模型阶数q',q,'往后预测%d个月的序列'%(n),predict))
        else:
            noiseAsse = '残差白噪声检验中p值为%f，小于0.05，认为该模型残差序列为非白噪声，模型有待改进'%noiseR
            cont.append((i,'AR模型阶数p',p,'MA模型阶数q',q,noiseAsse))   
    return cont