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

巴特西

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

最新文章

热门文章