恐怖袭击等级预测量化与ARMIA时间序列建模的例子

一.恐怖袭击的全球分布量化图:(量化分类由k-means算法得)

# coding:utf-8

import pandas as pd
import mpl_toolkits.basemap             #地图只在Spyder中加载是成功的!!!
import matplotlib.pyplot as plt
import seaborn as sns

plt.style.use('classic')  # 此句设置出雨林绿和中天蓝
#sns.set()

data_file1 = "../data/附件1.xlsx"
data_file2 = "../temp/q1_聚类分级结果.xlsx"
df1 = pd.read_excel(data_file1)    #内含经纬度数据
df2 = pd.read_excel(data_file2)    #量化分级的数据F值和分级
dfdf = pd.merge(df1, df2)
columns1 = dfdf.columns.tolist()
dfdf = dfdf[["eventid", "latitude", "longitude", "F值", "分级"]]

plt.subplots(figsize=(20, 9))
basemap = mpl_toolkits.basemap.Basemap()
basemap.drawcoastlines()
basemap.drawcountries(linewidth=1.5)

# cm = plt.cm.get_cmap('gist_rainbow')
# 直接将聚类的簇号(按照簇的大小来量化分级,簇越大袭击的等级越低)赋给颜色作为能量渐变值
colors = {1: 'red', 2: 'orange', 3: 'yellow', 4: 'cyan', 5: 'green'}
dengji = {1: '一级恐怖袭击', 2: '二级恐怖袭击', 3: '三级恐怖袭击', 4: '四级恐怖袭击', 5: '五级恐怖袭击'}
daxiao = {1: 50, 2: 40, 3: 30, 4: 10, 5: 3}
for i in range(len(colors)):
    px = dfdf['longitude'][dfdf["分级"] == i + 1]
    py = dfdf['latitude'][dfdf["分级"] == i + 1]
    plt.scatter(px, py, c=colors[i + 1], vmin=0, vmax=20000, s=daxiao[i + 1], label=dengji[i + 1])
    # plt.scatter(dfdf['longitude'], dfdf['latitude'], c=dfdf["分级"].apply(lambda x:colors[x]), vmin=0, vmax=200, s=9)
    plt.rcParams['font.sans-serif'] = ['SimHei']  # 标题不能显示汉字,这么处理
    plt.title('1998-2017世界恐怖袭击案发地分布图')
    plt.legend()  # 这里怎么写???
    plt.savefig('../img/world'+str(i)+'.png')
    #plt.show()

效果图:

二.ARMIA时间序列例程

例1：CO2回归预测

"""例1:时间序列建模:ARIMA(差分自回归移动平均模型)"""          
import warnings
import itertools
import pandas as pd
import numpy as np
import statsmodels.api as sm
import matplotlib.pyplot as plt
plt.style.use('fivethirtyeight')
data = sm.datasets.co2.load_pandas()
y = data.data 
y.plot(figsize=(15, 6))
plt.show()      

#1.数据预处理。 每周数据可能很棘手，因为它是一个很短的时间，所以让我们使用每月平均值。 我们将使用resample函数进行转换。 为了简单起见，我们还可以使用fillna()函数来确保我们的时间序列中没有缺少值 
# The 'MS' string groups the data in buckets by start of the month
y = y['co2'].resample('MS').mean()    #将数据按月分组取平均(这是一种累加或者累减的思想)
# The term bfill means that we use the value before filling in missing values
y = y.fillna(y.bfill()) #使用缺失值的前一个来做填充
print(y)  
# 注：y.shape从(2284,1) -> (526,)

#2.探索这个时间序列e作为数据可视化：
y.plot(figsize=(15, 6))
plt.show()

#3.ARIMA时间序列模型调参
# Define the p, d and q parameters to take any value between 0 and 2
p = d = q = range(0, 2)
 
# Generate all different combinations of p, q and q triplets
pdq = list(itertools.product(p, d, q))
"""
	[(0, 0, 0),
	 (0, 0, 1),
	 (0, 1, 0),
	 (0, 1, 1),
	 (1, 0, 0),
	 (1, 0, 1),
	 (1, 1, 0),
	 (1, 1, 1)]
""" 
#3.ARIMA时间序列模型调参
# Define the p, d and q parameters to take any value between 0 and 2
p = d = q = range(0, 2)
 
# Generate all different combinations of p, q and q triplets
pdq = list(itertools.product(p, d, q))
 
# Generate all different combinations of seasonal p, q and q triplets
#产生000--111的8个二进制编码
seasonal_pdq = [(x[0], x[1], x[2], 12) for x in list(itertools.product(p, d, q))]
 
print('Examples of parameter combinations for Seasonal ARIMA...')
print('SARIMAX: {} x {}'.format(pdq[1], seasonal_pdq[1]))
print('SARIMAX: {} x {}'.format(pdq[1], seasonal_pdq[2]))
print('SARIMAX: {} x {}'.format(pdq[2], seasonal_pdq[3]))
print('SARIMAX: {} x {}'.format(pdq[2], seasonal_pdq[4]))

# 使用网格搜索，我们已经确定了为我们的时间序列数据生成最佳拟合模型的参数集。 我们可以更深入地分析这个特定的模型。 
warnings.filterwarnings("ignore") # specify to ignore warning messages
for param in pdq:
    for param_seasonal in seasonal_pdq:
        try:
            mod = sm.tsa.statespace.SARIMAX(y,
                                            order=param,
                                            seasonal_order=param_seasonal,
                                            enforce_stationarity=False,
                                            enforce_invertibility=False)
 
            results = mod.fit()
 
            print('ARIMA{}x{}12 - AIC:{}'.format(param, param_seasonal, results.aic))
        except:
            continue
# 我们的代码的输出表明， SARIMAX(1, 1, 1)x(1, 1, 1, 12)产生最低的AIC值为277.78。 因此，我们认为这是我们考虑过的所有模型中的最佳选择。 

#我们首先将最佳参数值插入到新的SARIMAX模型中
mod = sm.tsa.statespace.SARIMAX(y,
                                order=(1, 1, 1),
                                seasonal_order=(1, 1, 1, 12),
                                enforce_stationarity=False,
                                enforce_invertibility=False)
 
results = mod.fit()
print(results.summary().tables[1])
# coef列显示每个特征的重量（即重要性）以及每个特征如何影响时间序列
# P>|z| 列通知我们每个特征重量的意义。 这里，每个重量的p值都低于或接近0.05 ，所以在我们的模型中保留所有权重是合理的

#在适合季节性ARIMA模型（以及任何其他模型）的情况下，运行模型诊断是非常重要的，以确保没有违反模型的假设。 plot_diagnostics对象允许我们快速生成模型诊断并调查任何异常行为
results.plot_diagnostics(figsize=(15, 12))
plt.show()

# 我们的主要关切是确保我们的模型的残差是不相关的，并且平均分布为零。 如果季节性ARIMA模型不能满足这些特性，这是一个很好的迹象，可以进一步改善。 
#在这种情况下，我们的模型诊断表明，模型残差正常分布如下：
#在右上图中，我们看到红色KDE线与N(0,1)行（其中N(0,1) ）是正态分布的标准符号，平均值0 ，标准偏差为1 ） 。 这是残留物正常分布的良好指示。
#左下角的qq图显示，残差（蓝点）的有序分布遵循采用N(0, 1)的标准正态分布采样的线性趋势。 同样，这是残留物正常分布的强烈指示。
#随着时间的推移（左上图）的残差不会显示任何明显的季节性，似乎是白噪声。 这通过右下角的自相关（即相关图）来证实，这表明时间序列残差与其本身的滞后版本具有低相关性。 
#这些观察结果使我们得出结论，我们的模型产生了令人满意的合适性，可以帮助我们了解我们的时间序列数据和预测未来价值。
#虽然我们有一个令人满意的结果，我们的季节性ARIMA模型的一些参数可以改变，以改善我们的模型拟合。 例如，我们的网格搜索只考虑了一组受限制的参数组合，所以如果我们拓宽网格搜索，我们可能会找到更好的模型。 

#我们已经获得了我们时间序列的模型，现在可以用来产生预测
#我们首先将预测值与时间序列的实际值进行比较，这将有助于我们了解我们的预测的准确性。 get_prediction()和conf_int()属性允许我们获得时间序列预测的值和相关的置信区间。
pred = results.get_prediction(start=pd.to_datetime('1998-01-01'), dynamic=False)
pred_ci = pred.conf_int()
# 上述规定需要从1998年1月开始进行预测。
#dynamic=False参数确保我们产生一步前进的预测，这意味着每个点的预测都将使用到此为止的完整历史生成。
#我们可以绘制二氧化碳时间序列的实际值和预测值，以评估我们做得如何。 注意我们如何在时间序列的末尾放大日期索引。 
ax = y['1990':].plot(label='observed')
pred.predicted_mean.plot(ax=ax, label='One-step ahead Forecast', alpha=.7)
 
ax.fill_between(pred_ci.index,
                pred_ci.iloc[:, 0],
                pred_ci.iloc[:, 1], color='k', alpha=.2)
 
ax.set_xlabel('Date')
ax.set_ylabel('CO2 Levels')
plt.legend()
plt.show()
#总体而言，我们的预测与真实价值保持一致，呈现总体增长趋势。

# 量化我们的预测的准确性也是有用的,使用MSE（均方误差）.
y_forecasted = pred.predicted_mean
y_truth = y['1998-01-01':]
 
# Compute the mean square error
mse = ((y_forecasted - y_truth) ** 2).mean()
print('The Mean Squared Error of our forecasts is {}'.format(round(mse, 2)))

# 然而，使用【动态预测】可以获得更好地表达我们的真实预测能力。 在这种情况下，我们只使用时间序列中的信息到某一点，之后，使用先前预测时间点的值生成预测。 
pred_dynamic = results.get_prediction(start=pd.to_datetime('1998-01-01'), dynamic=True, full_results=True)
pred_dynamic_ci = pred_dynamic.conf_int(