ARIMA模型---时间序列分析---温度预测

（图片来自百度）

数据

分析数据第一步还是套路------画图

数据看上去比较平整，但是由于数据太对看不出具体情况，于是将只取前300个数据再此画图

这数据看上去很不错，感觉有隐藏周期的意思

代码

#coding:utf-8
import csv
import matplotlib.pyplot as plt

def read_csv_data(aim_list_1, aim_list_2, file_name):
    i = 0
    csv_file = csv.reader(open(file_name,'r'))
    for data in csv_file:
        if (i == 0):
            i += 1
            continue
        aim_list_1.append(float(data[1]))
        aim_list_2.append(data[3])
    return
def plot_picture(x, y):
    plt.xlabel('x')
    plt.ylabel('y')
    plt.plot(x, y)
    plt.show()
    return
if __name__ == '__main__':
    temp = []
    tim = []
    file_name = 'C:/Users/lichaoxing/Desktop/testdata.csv'
    read_csv_data(temp, tim, file_name)
    plot_picture(tim[:300], temp[:300])

使用ARIMA模型（ARMA）

第一步观察数据是否是平稳序列，通过上图可以看出是平稳的

如果不平稳，则需要进行预处理，方法有   对数变换    差分

对于平稳的时间序列可以直接使用ARMA(p, q)模型进行拟合

ARMA (p, q) :  AR(p) + MA(q)

此时参数p和q的确定可以通过观察ACF和PACF图来确定

通过观察PACF图可以看出，阶数为9也就是p=9，这里ACF图看出自相关呈现震荡下降收敛，但是怎么决定出q，我没太明白，这里姑且拍脑袋才一个吧就q=3

但是这里我遇到了一个问题，没有搞懂，就是平稳的序列，如果我进行一阶差分后应该仍然是平稳的序列，但是这个时候我又画了一个ACF与PACF图，竟然是下图这样，lag的范围是-0.04到0.04（不懂）

---

lag的范围是-0.04到0.04的问题原因（修改于再次使用此模型）

原因：

当时，我使用的是一阶差分，也就是让数据的后一个值减去前一个值得到新的值，这样就会导致第一个值变为缺失值（下面的数据是再此使用此模型时的数据，与原博客数据无关）

就是因为此处的值为缺失值，导致绘制ACF与PACF时数据有问题而无法成功显示

解决办法，在绘制上述图形前，将第一个数据去除：

dta= dta.diff(1)
dta = dta.truncate(before= ym[1])#删除第一个缺失值

---

其实还有就是使用ADF检验，得到的结果如图，这个p值很小===》平稳

画图代码

def acf_pacf(temp, tim):
    x = tim
    y = temp
    dta = pd.Series(y, index = pd.to_datetime(x))
    fig = plt.figure(figsize=(9,6))
    ax1 = fig.add_subplot(211)
    fig = sm.graphics.tsa.plot_acf(dta,lags=50,ax=ax1)
    ax2 = fig.add_subplot(212)
    fig = sm.graphics.tsa.plot_pacf(dta,lags=50,ax=ax2)
    show()

ADF检验代码

def test_stationarity(timeseries):
    dftest = adfuller(timeseries, autolag='AIC')
    return dftest[1]

这里先使用ARMA(9,3)来实验测试一下效果，取前300个数据中的前250个作为train，后面的作为test

效果

可以说这个模型是真的强大，预测的还是十分准确的

代码

def test_300(temp, tim):

    x = tim[0:300]
    y = temp[0:300]
    dta = pd.Series(y[0:249], index = pd.to_datetime(x[0:249]))

    fig = plt.figure(figsize=(9,6))
    ax1 = fig.add_subplot(211)
    fig = sm.graphics.tsa.plot_acf(dta,lags=30,ax=ax1)
    ax2 = fig.add_subplot(212)
    fig = sm.graphics.tsa.plot_pacf(dta,lags=30,ax=ax2)

    arma_mod = sm.tsa.ARMA(dta, (9, 3)).fit(disp = 0)

    predict_sunspots = arma_mod.predict(x[200], x[299], dynamic=True)

    fig, ax = plt.subplots(figsize=(9, 6))
    ax = dta.ix[x[0]:].plot(ax=ax)
    predict_sunspots.plot(ax=ax)

    show()

其实，可以通过代码来自动的选择p和q的值，依据BIC准则，目标就是bic越小越好

代码

def proper_model(timeseries, maxLag):
    init_bic = 100000000
    init_properModel = None
    for p in np.arange(maxLag):
        for q in np.arange(maxLag):
            model = ARMA(timeseries, order=(p, q)) 
            try:
                results_ARMA = model.fit(disp = 0, method='css')
            except:
                continue
            bic = results_ARMA.bic
            if bic < init_bic:
                init_properModel = results_ARMA
                init_bic = bic
    return init_properModel

遇到的问题，预测时predict函数没怎么使用明白

当写于某些预测区间的时候，会报 “start”或“end”的相关错误，还有一个函数forcast，这个函数使用就是forcast(N)：预测后面N个值

返回的是预测值（array型）标准误差（array型）置信区间（array型）

还有：

对于构造时间序列，时间可以是时间格式：如 “2018-01-01”  或者就是个时间戳，在用时间戳的时候，其实在序列里它会自动识别时间戳，并加上起始时间1970-01-01 00:00:01

形式

附录（代码）

预测一序列中某一点的值

#coding:utf-8
import csv
import time
import pandas as pd
import numpy as np
from statsmodels.tsa.arima_model import ARMA
import argparse
import warnings

warnings.filterwarnings('ignore')
def timestamp_datatime(value):
    value = time.localtime(value)
    dt = time.strftime('%Y-%m-%d %H:%M',value)
    return dt

def time_timestamp(my_date):
    my_date_array = time.strptime(my_date,'%Y-%m-%d %H:%M')
    my_date_stamp = time.mktime(my_date_array)
    return my_date_stamp

def read_csv_data(aim_list_1, aim_list_2, file_name):
    i = 0
    csv_file = csv.reader(open(file_name,'r'))
    for data in csv_file:
        if (i == 0):
            i += 1
            continue
        aim_list_1.append(float(data[1]))   #1：温度  2：湿度
        dt = int(data[3])
        aim_list_2.append(dt)
    return

def proper_model(timeseries, maxLag):
    init_bic = 100000000
    init_properModel = None
    for p in np.arange(maxLag):
        for q in np.arange(maxLag):
            model = ARMA(timeseries, order=(p, q))   #bug
            try:
                results_ARMA = model.fit(disp = 0, method='css')
            except:
                continue
            bic = results_ARMA.bic
            if bic < init_bic:
                init_properModel = results_ARMA
                init_bic = bic
    return init_properModel

def test_300(temp, tim, time_in):
    
    x = []
    y = []
    end_index = len(tim)
    for i in range(0, len(tim)):
        if (time_in - (tim[i]) < 300):
            end_index = i
            break
    if (end_index < 100):
        x = tim[0: end_index]
        y = temp[0: end_index]
    else:
        x = tim[end_index - 100: end_index]
        y = temp[end_index - 100: end_index]
    
    tidx = pd.DatetimeIndex(x, freq='infer')
    dta = pd.Series(y, index = tidx)
    print(dta)
    arma_mod = proper_model(dta, 9)

    predict_sunspots = arma_mod.forecast(1)
    return predict_sunspots[0]
def predict_temperature(file_name, time_in):

    temp = []
    tim = []
    read_csv_data(temp, tim, file_name)
    
    result_temp = test_300(temp, tim, time_in)
    
    return result_temp
if __name__ == '__main__':
    
    parser = argparse.ArgumentParser()
    parser.add_argument('-f', action='store', dest='file_name')
    parser.add_argument('-t', action='store', type = int, dest='time_')

    args = parser.parse_args()
    file_name = args.file_name
    time_in = args.time_

    result_temp = predict_temperature(file_name, time_in)

    print ('the temperature is %f ' % result_temp)

在上面的代码中，预测某一点的值我采用序列中此点的前100个点作为训练集

如果给出待预测的多个点，由于每次都要计算模型的p和q以及拟合模型，时间会很慢，于是考虑将给定的待预测时间点序列切割成小段，使每一段中最大与最小的时间间隔在某一范围内

在使用forcast(n)函数一次预测多点，然后在预测值中找到与待预测的时间值相近的值，速度大大提升，思路如图

代码

#coding:utf-8
import csv
#import time
import pandas as pd
import numpy as np
from statsmodels.tsa.arima_model import ARMA
import warnings

warnings.filterwarnings('ignore')

def proper_model(timeseries, maxLag):
    init_bic = 1000000000
    init_p = 1
    init_q = 1
    for p in np.arange(maxLag):
        for q in np.arange(maxLag):
            model = ARMA(timeseries, order=(p, q))
            try:
                results_ARMA = model.fit(disp = 0, method='css')
            except:
                continue
            bic = results_ARMA.bic
            if bic < init_bic:
                init_p = p
                init_q = q
                init_bic = bic
    return init_p, init_q

def read_csv_data(file_name, clss = 1):
    i = 0
    aim_list_1 = [] #temperature(1) or humidity(2)
    aim_list_2 = [] #time
    csv_file = csv.reader(open(file_name,'r'))
    for data in csv_file:
        if (i == 0):
            i += 1
            continue
        aim_list_1.append(float(data[clss]))
        dt = int(data[3])
        aim_list_2.append(dt)
    
    tidx = pd.DatetimeIndex(aim_list_2, freq = None)
    dta = pd.Series(aim_list_1, index = tidx)
    init_p, init_q = proper_model(dta[:aim_list_2[100]], 9)
    return init_p, init_q, aim_list_2, dta

def for_kernel(p, q, tim, dta, tmp_time_list, result_dict):
    interval = 20
    end_index = len(tim) - 1
    for i in range(0, len(tim)):
        if (tmp_time_list[0]["time"] - tim[i] < tim[1] - tim[0]):
            end_index = i
            break

    if (end_index < 100):
        dta = dta.truncate(after = tim[end_index])
    else:
        dta = dta.truncate(before= tim[end_index - 101], after = tim[end_index])

    arma_mod = ARMA(dta, order=(p, q)).fit(disp = 0, method='css')
    #为未来interval天进行预测， 返回预测结果， 标准误差， 和置信区间
    predict_sunspots = arma_mod.forecast(interval)
    ####################################
    for tim_i in tmp_time_list:
        for tim_ in tim:
            if tim_i["time"] - tim_ >= 0 and tim_i["time"] - tim_ < tim[1] - tim[0]:
                result_dict[tim_i["time"]] = predict_sunspots[0][tim.index(tim_) - end_index]
    return

def kernel(p, q, tim, dta, time_in_list):
    interval = 20
    time_first = time_in_list[0]
    det_time = tim[1] - tim[0]
    result_dict = {}
    tmp_time_list = []
    for time_ in time_in_list:
        if time_first["time"] + det_time * interval > time_["time"]:
            tmp_time_list.append(time_)
            continue

        time_first = time_
        for_kernel(p, q, tim, dta, tmp_time_list, result_dict)
        tmp_time_list = []
        tmp_time_list.append(time_first)
    
    for_kernel(p, q, tim, dta, tmp_time_list, result_dict)
    
    return result_dict

def predict_temperature(file_name, time_in_list, clss = 1):
    p, q, tim, dta = read_csv_data(file_name, clss)

    result_temp_dict = kernel(p, q, tim, dta, time_in_list)
    
    return result_temp_dict

def predict_humidity(file_name, time_in_list, clss = 2):
    p, q, tim, dta = read_csv_data(file_name, clss)
    result_humi_dict = kernel(p, q, tim, dta, time_in_list)

    return result_humi_dict

if __name__ == '__main__':
    

    file_name = "testdata.csv"
    time_in = [{"time":1530419271,"temp":"","humi":""},{"time":1530600187,"temp":"","humi":""},{"time":1530825809,"temp":"","humi":""}]

    #time_in = [{"time":1530600187,"temp":"","humi":""},]
    result_temp = predict_temperature(file_name, time_in)
    print(result_temp)

由于后续又改动了需求，需要预测温度以及湿度，完成了项目在github

https://github.com/xinglicha0/armamodel

---