机器学习大作业：多元时间序列预测Sarima
提前说明，由于图片太多了，我不可能一个一个上传，所以我整合并上传至了优快云。
初步处理数据和数据的可视化
原始数据分别保存成两个CSV文件
原本的xlsx文件包括两个工作文件，我们在转化为csv文件时出现问题，所以只能将原本的origindata分别做成两个:
annual_data.csv与monthly_data.csv
注意
我们应该如何分开保存？
原本的表格文件内含两个工作簿，简单地“另存为”会导致非常难以理解的错误，出现无法修正的BUG
我们不应该简单地在原本的xlax进行“另存为”
而是应该另外开一个工作簿，然后将相应的工作簿复制粘贴到新的csv中，进行保存
保存的方式如下（记得设置成UTF-8方式）：
初步判断csv文件中数据形式和形状（以monthly_data为例）
import numpy as np
import pandas as pd
data = pd.read_csv('monthlydata.csv')
print(data.shape)
data.head()
(594, 10)
| Month | Coal Electric Power Sector CO2 Emissions | Natural Gas Electric Power Sector CO2 Emissions | Distillate Fuel, Including Kerosene-Type Jet Fuel, Oil Electric Power Sector CO2 Emissions | Petroleum Coke Electric Power Sector CO2 Emissions | Residual Fuel Oil Electric Power Sector CO2 Emissions | Petroleum Electric Power Sector CO2 Emissions | Geothermal Energy Electric Power Sector CO2 Emissions | Non-Biomass Waste Electric Power Sector CO2 Emissions | Total Energy Electric Power Sector CO2 Emissions | |
|---|---|---|---|---|---|---|---|---|---|---|
| 0 | 1973 January | 73.112 | 12.163 | 2.397 | 0.128 | 23.698 | 26.223 | Not Available | Not Available | 111.499 |
| 1 | 1973 February | 65.369 | 11.696 | 2.080 | 0.106 | 19.886 | 22.071 | Not Available | Not Available | 99.137 |
| 2 | 1973 March | 65.006 | 13.979 | 1.181 | 0.083 | 18.851 | 20.115 | Not Available | Not Available | 99.100 |
| 3 | 1973 April | 61.717 | 14.612 | 1.031 | 0.130 | 15.783 | 16.944 | Not Available | Not Available | 93.273 |
| 4 | 1973 May | 62.686 | 17.326 | 0.957 | 0.167 | 16.920 | 18.044 | Not Available | Not Available | 98.057 |
把文件中每个头进行翻译：
煤电力部门二氧化碳排放
天然气电力部门二氧化碳排放
蒸馏燃料，包括煤油型喷气燃料、石油电力部门二氧化碳排放
石油焦电力部门二氧化碳排放
残余燃料油电力部门二氧化碳排放
石油电力部门二氧化碳排放
地热能电力部门二氧化碳排放
非生物质废物电力部门二氧化碳排放
总能源电力部门二氧化碳排放
试验：对monthlydata的数据第一列"Coal Electric Power Sector CO2 Emissions"进行趋势分解 出现了错误
import statsmodels
from statsmodels.tsa.seasonal import seasonal_decompose
from dateutil.parser import parse
import pandas as pd
import matplotlib.pyplot as plt
# 导入数据
df = pd.read_csv('monthlydata.csv')
dates = pd.DatetimeIndex([parse(d).strftime('%Y %B') for d in df['Month']])
df.set_index(dates, inplace=True)
# 分解模型--乘法
result = seasonal_decompose(df['Coal Electric Power Sector CO2 Emissions'],
model='multiplicative')
# 画图
plt.rc('font', family='SimHei', size=13)
plt.rcParams.update({
'figure.figsize': (10, 10)})
result.plot().suptitle('')
plt.show()
---------------------------------------------------------------------------
ValueError Traceback (most recent call last)
Cell In[2], line 13
10 df.set_index(dates, inplace=True)
12 # 分解模型--乘法
---> 13 result = seasonal_decompose(df['Coal Electric Power Sector CO2 Emissions'], model='multiplicative')
15 # 画图
16 plt.rc('font', family='SimHei', size=13)
File ~\anaconda3\envs\pytorch\lib\site-packages\statsmodels\tsa\seasonal.py:153, in seasonal_decompose(x, model, filt, period, two_sided, extrapolate_trend)
150 nobs = len(x)
152 if not np.all(np.isfinite(x)):
--> 153 raise ValueError("This function does not handle missing values")
154 if model.startswith("m"):
155 if np.any(x <= 0):
ValueError: This function does not handle missing values
出现了错误，提示我们不能读取空白的数据
这里需要注意在csv文件的最下方是丢失了数据，需要我们预测的，所以只能读取到2021 December
正确地趋势分解 Coal Electric Power Sector CO2 Emissions 列 和 Natural Gas Electric Power Sector CO2 Emissions 列
import statsmodels
from statsmodels.tsa.seasonal import seasonal_decompose
from dateutil.parser import parse
import pandas as pd
import matplotlib.pyplot as plt
# 导入数据
df = pd.read_csv('monthlydata.csv')
# 将"Month"列转换为日期时间对象
df['Month'] = pd.to_datetime(df['Month'], format='%Y %B')
# 筛选出截止到"2021 December"的数据
df = df[df['Month'] <= '2021-12-01']
# 将"Month"列设置为索引
df.set_index('Month', inplace=True)
# 分解模型--乘法
result = seasonal_decompose(df['Coal Electric Power Sector CO2 Emissions'],
model='multiplicative')
# 画图
plt.rc('font', family='SimHei', size=13)
plt.rcParams.update({
'figure.figsize': (20, 40)})
#plt.rcParams.update({'figure.figsize': (10,10)})
result.plot().suptitle('')
plt.show()
import statsmodels
from statsmodels.tsa.seasonal import seasonal_decompose
from dateutil.parser import parse
import pandas as pd
import matplotlib.pyplot as plt
# 导入数据
df = pd.read_csv('monthlydata.csv')
# 将"Month"列转换为日期时间对象
df['Month'] = pd.to_datetime(df['Month'], format='%Y %B')
# 筛选出截止到"2021 December"的数据
df = df[df['Month'] <= '2021-12-01']
# 将"Month"列设置为索引
df.set_index('Month', inplace=True)
# 分解模型--乘法
result = seasonal_decompose(
df['Natural Gas Electric Power Sector CO2 Emissions'],
model='multiplicative')
# 画图
plt.rc('font', family='SimHei', size=13)
plt.rcParams.update({
'figure.figsize': (20, 40)})
#plt.rcParams.update({'figure.figsize': (10,10)})
result.plot().suptitle('')
plt.show()
初步地分析了两列数据，可以看到Trend清晰可见，而季节性特征非常明显
初步解决问题：哪些指标之间又相互关联的关系？
引入库（start）
import warnings
import itertools
import pandas as pd
import numpy as np
from sklearn.preprocessing import StandardScaler
import statsmodels
from statsmodels.tsa.seasonal import seasonal_decompose
from dateutil.parser import parse
import statsmodels.api as sm
import seaborn as sns
import matplotlib.pyplot as plt
plt.style.use('fivethirtyeight')
#在这里，使用了 "fivethirtyeight" 样式表。"fivethirtyeight" 是一种图形设计风格，得名于美国著名的政治和经济新闻网站 FiveThirtyEight。
#这种风格通常以清晰、易读、注重数据可视化的特点而闻名。
处理monthlydata.csv并返回一个data.csv
def handle_monlydata_csv(file_path):
data = pd.read_csv(file_path)
# 将"Month"列转换为日期时间对象
data['Month'] = pd.to_datetime(data['Month'], format='%Y %B')
# 筛选出截止到"2021 December"的数据
data = data[data['Month'] <= '2021-12-01']
text = '''Coal Electric Power Sector CO2 Emissions 煤电能源部门（二氧化碳排放） 替换为 Emission_1
Natural Gas Electric Power Sector CO2 Emissions 天然气电力部门二氧化碳排放 替换为 Emission_2
Distillate Fuel, Including Kerosene-Type Jet Fuel, Oil Electric Power Sector CO2 Emissions 蒸馏燃料，包括喷气式燃料、石油电力部门二氧化碳排放 替换为Emission_3
Petroleum Coke Electric Power Sector CO2 Emissions 石油焦电力部门二氧化碳排放 替换为Emission_4
Residual Fuel Oil Electric Power Sector CO2 Emissions 残余燃油电力部门二氧化碳排放 替换为Emission_5
Petroleum Electric Power Sector CO2 Emissions 石油电力部门二氧化碳排放 替换为 Emission_6
Geothermal Energy Electric Power Sector CO2 Emissions 地热能电力部门二氧化碳排放 替换为Emission_7
Non-Biomass Waste Electric Power Sector CO2 Emissions 非生物质废物电力部门二氧化碳排放 替换为Emission_8
Total Energy Electric Power Sector CO2 Emissions 总能源电力部门二氧化碳排放 替换为Emission_9'''
print(text)
data.rename(
columns={
'Coal Electric Power Sector CO2 Emissions': 'Emission_1'},
inplace=True)
data.rename(columns={
'Natural Gas Electric Power Sector CO2 Emissions': 'Emission_2'
},
inplace=True)
data.rename(columns={
'Distillate Fuel, Including Kerosene-Type Jet Fuel, Oil Electric Power Sector CO2 Emissions':
'Emission_3'
},
inplace=True)
data.rename(columns={
'Petroleum Coke Electric Power Sector CO2 Emissions':
'Emission_4'
},
inplace=True)
data.rename(columns={
'Residual Fuel Oil Electric Power Sector CO2 Emissions':
'Emission_5'
},
inplace=True)
data.rename(columns={
'Petroleum Electric Power Sector CO2 Emissions': 'Emission_6'
},
inplace=True)
data.rename(columns={
'Geothermal Energy Electric Power Sector CO2 Emissions':
'Emission_7'
},
inplace=True)
data.rename(columns={
'Non-Biomass Waste Electric Power Sector CO2 Emissions':
'Emission_8'
},
inplace=True)
data.rename(columns={
'Total Energy Electric Power Sector CO2 Emissions':
'Emission_9'
},
inplace=True)
#替换不可用数据
data.replace('Not Available', np.nan, inplace=True)
# 将"Month"列设置为数据框的索引
data.set_index('Month', inplace=True)
data = data.astype(float)
#前向填充，使用中位数
data.ffill(inplace=True)
data['Emission_7'] = data['Emission_7'].fillna(data['Emission_7'].median())
data['Emission_8'] = data['Emission_8'].fillna(data['Emission_8'].median())
print(data.index)
#打印索引，我们发现索引正确，从1973到2021
print(data)
#检验成功 我们截取了1973年1月到2021年12月的数据
return data
file_path = 'monthlydata.csv'
handle_monlydata_csv = handle_monlydata_csv(file_path)
handle_monlydata_csv.to_csv('handle.csv')
Coal Electric Power Sector CO2 Emissions 煤电能源部门（二氧化碳排放） 替换为 Emission_1
Natural Gas Electric Power Sector CO2 Emissions 天然气电力部门二氧化碳排放 替换为 Emission_2
Distillate Fuel, Including Kerosene-Type Jet Fuel, Oil Electric Power Sector CO2 Emissions 蒸馏燃料，包括喷气式燃料、石油电力部门二氧化碳排放 替换为Emission_3
Petroleum Coke Electric Power Sector CO2 Emissions 石油焦电力部门二氧化碳排放 替换为Emission_4
Residual Fuel Oil Electric Power Sector CO2 Emissions 残余燃油电力部门二氧化碳排放 替换为Emission_5
Petroleum Electric Power Sector CO2 Emissions 石油电力部门二氧化碳排放 替换为 Emission_6
Geothermal Energy Electric Power Sector CO2 Emissions 地热能电力部门二氧化碳排放 替换为Emission_7
Non-Biomass Waste Electric Power Sector CO2 Emissions 非生物质废物电力部门二氧化碳排放 替换为Emission_8
Total Energy Electric Power Sector CO2 Emissions 总能源电力部门二氧化碳排放 替换为Emission_9
DatetimeIndex(['1973-01-01', '1973-02-01', '1973-03-01', '1973-04-01',
'1973-05-01', '1973-06-01', '1973-07-01', '1973-08-01',
'1973-09-01', '1973-10-01',
...
'2021-03-01', '2021-04-01', '2021-05-01', '2021-06-01',
'2021-07-01', '2021-08-01', '2021-09-01', '2021-10-01',
'2021-11-01', '2021-12-01'],
dtype='datetime64[ns]', name='Month', length=588, freq=None)
Emission_1 Emission_2 Emission_3 Emission_4 Emission_5 \
Month
1973-01-01 73.112 12.163 2.397 0.128 23.698
1973-02-01 65.369 11.696 2.080 0.106 19.886
1973-03-01 65.006 13.979 1.181 0.083 18.851
1973-04-01 61.717 14.612 1.031 0.130 15.783
1973-05-01 62.686 17.326 0.957 0.167 16.920
... ... ... ... ... ...
2021-08-01 101.985 68.814 0.346 0.875 0.510
2021-09-01 80.443 54.323 0.251 0.732 0.416
2021-10-01 64.520 51.508 0.286 0.722 0.346
2021-11-01 59.326 48.144 0.301 0.919 0.324
2021-12-01 62.391 48.368 0.339 0.696 0.342
Emission_6 Emission_7 Emission_8 Emission_9
Month
1973-01-01 26.223 0.032 0.896 111.499
1973-02-01 22.071 0.032 0.896 99.137
1973-03-01 20.115 0.032 0.896 99.100
1973-04-01 16.944 0.032 0.896 93.273
1973-05-01 18.044 0.032 0.896 98.057
... ... ... ... ...
2021-08-01 1.731 0.040 0.892 173.463
2021-09-01 1.400 0.039 0.863 137.067
2021-10-01 1.353 0.040 0.892 118.313
2021-11-01 1.544 0.039 0.863 109.915
2021-12-01 1.377 0.040 0.892 113.068
[588 rows x 9 columns]
检验结果
我们从打印的结果可以看出，我们成功设置了handle.csv
handle.csv除了对原本的monthlydata.csv的NOT AVAILABLE的数据进行了修改，经过前向替换成中位数的数据
下面我们读取并处理handle.csv的数据
data = pd.read_csv('handle.csv')
data['Month'] = pd.to_datetime(data['Month'])
data.set_index('Month', inplace=True)
数据标准化并可视化
data_extra = data
# 使用StandardScaler对选定的列进行标准化
scaler = StandardScaler()
data_extra_standardized = pd.DataFrame(scaler.fit_transform(data_extra),
columns=data_extra.columns)
import itertools
# 获取所有可能的两两组合
column_combinations = list(
itertools.combinations(data_extra_standardized.columns, 2))
# 遍历组合并绘制图形
for column1, column2 in column_combinations:
# 创建一个新的图形
plt.figure(figsize=(30, 15))
# 绘制两个数据列
plt.plot(data_extra_standardized[column1], label=column1)
plt.plot(data_extra_standardized[column2], label=column2)
# 添加图例
plt.legend()
# 添加标题和轴标签
plt.title('Visualization of {} and {}'.format(column1, column2))
plt.xlabel('Index')
plt.ylabel('Values')
# 显示图形
plt.show()
print(data.index)
#打印索引，我们发现索引正确，从1973到2021
print(data)
#检验成功 我们截取了1973年1月到2021年12月的数据
DatetimeIndex(['1973-01-01', '1973-02-01', '1973-03-01', '1973-04-01',
'1973-05-01', '1973-06-01', '1973-07-01', '1973-08-01',
'1973-09-01', '1973-10-01',
...
'2021-03-01', '2021-04-01', '2021-05-01', '2021-06-01',
'2021-07-01', '2021-08-01', '2021-09-01', '2021-10-01',
'2021-11-01', '2021-12-01'],
dtype='datetime64[ns]', name='Month', length=588, freq=None)
Emission_1 Emission_2 Emission_3 Emission_4 Emission_5 \
Month
1973-01-01 73.112 12.163 2.397 0.128 23.698
1973-02-01 65.369 11.696 2.080 0.106 19.886
1973-03-01 65.006 13.979 1.181 0.083 18.851
1973-04-01 61.717 14.612 1.031 0.130 15.783
1973-05-01 62.686 17.326 0.957 0.167 16.920
... ... ... ... ... ...
2021-08-01 101.985 68.814 0.346 0.875 0.510
2021-09-01 80.443 54.323 0.251 0.732 0.416
2021-10-01 64.520 51.508 0.286 0.722 0.346
2021-11-01 59.326 48.144 0.301 0.919 0.324
2021-12-01 62.391 48.368 0.339 0.696 0.342
Emission_6 Emission_7 Emission_8 Emission_9
Month
1973-01-01 26.223 0.032 0.896 111.499
1973-02-01 22.071 0.032 0.896 99.137
1973-03-01 20.115 0.032 0.896 99.100
1973-04-01 16.944 0.032 0.896 93.273
1973-05-01 18.044 0.032 0.896 98.057
... ... ... ... ...
2021-08-01 1.731 0.040 0.892 173.463
2021-09-01 1.400 0.039 0.863 137.067
2021-10-01 1.353 0.040 0.892 118.313
2021-11-01 1.544 0.039 0.863 109.915
2021-12-01 1.377 0.040 0.892 113.068
[588 rows x 9 columns]
分析对比，初步得出结论
分析标准化的数据，我们可以初步得出结论：
Emission_1和Emission_9可能存在关联
Emission_3和Emission_5可能存在关联
Emission_3和Emission_6可能存在关联
Emission_5和Emission_6可能存在关联
其中1和9,5和6应该是强关联
# 分解模型--乘法 并单独分解出某一列数据的趋势
result1 = seasonal_decompose(data['Emission_1'], model='multiplicative')
result2 = seasonal_decompose(data['Emission_9'], model='multiplicative')
result3 = seasonal_decompose(data['Emission_3'], model='multiplicative')
result4 = seasonal_decompose(data['Emission_5'], model='multiplicative')
result5 = seasonal_decompose(data['Emission_6'], model='multiplicative')
trend1 = result1.trend
trend2 = result2.trend
trend3 = result3.trend
trend4 = result4.trend
trend5 = result5.trend
# 画图
plt.rc('font', family='SimHei', size=13)
plt.rcParams.update({
'figure.figsize': (30, 15)})
# 绘制第一个数据列的趋势
plt.plot(trend1, label='Trend - Emission_1')
# 绘制第二个数据列的趋势
plt.plot(trend2, label='Trend - Emission_9')
# 添加图例
plt.legend()
# 添加标题和轴标签
plt.title('Trend_VS')
plt.xlabel('Date')
plt.ylabel('Trend Values')
# 显示图形
plt.show()
# 绘制第一个数据列的趋势
plt.plot(trend3, label='Trend - Emission_3')
# 绘制第二个数据列的趋势
plt.plot(trend4, label='Trend - Emission_5')
# 添加图例
plt.legend()
# 添加标题和轴标签
plt.title('Trend_VS')
plt.xlabel('Date')
plt.ylabel('Trend Values')
# 显示图形
plt.show()
# 绘制第一个数据列的趋势
plt.plot(trend3, label='Trend - Emission_3')
# 绘制第二个数据列的趋势
plt.plot(trend5, label='Trend - Emission_6')
# 添加图例
plt.legend()
# 添加标题和轴标签
plt.title('Trend_VS')
plt.xlabel('Date')
plt.ylabel('Trend Values')
# 显示图形
plt.show()
# 绘制第一个数据列的趋势
plt.plot(trend4, label='Trend - Emission_5')
# 绘制第二个数据列的趋势
plt.plot(trend5, label='Trend - Emission_6')
# 添加图例
plt.legend()
# 添加标题和轴标签
plt.title('Trend_VS')
plt.xlabel('Date')
plt.ylabel('Trend Values')
# 显示图形
plt.show()
再次进行比对
Emission_3&Emission_4存在相关情况
Emission_5&Emission_6存在相关情况
数据的预测
引入库
import warnings
import itertools
import pandas as pd
import numpy as np
from sklearn.preprocessing import StandardScaler
import statsmodels
from statsmodels.tsa.seasonal import seasonal_decompose
from dateutil.parser import parse
import statsmodels.api as sm
import seaborn as sns
import matplotlib.pyplot as plt
from statsmodels.tsa.stattools import adfuller
plt.rc('font', family='SimHei', size=13)
plt.style.use('fivethirtyeight')
#在这里，使用了 "fivethirtyeight" 样式表。"fivethirtyeight" 是一种图形设计风格，得名于美国著名的政治和经济新闻网站 FiveThirtyEight。
#这种风格通常以清晰、易读、注重数据可视化的特点而闻名。
warnings.filterwarnings("ignore")
读取handle.csv
data = pd.read_csv('handle.csv')
data['Month'] = pd.to_datetime(data['Month'])
data.set_index('Month', inplace=True)
print(data.index)
#打印索引，我们发现索引正确，从1973到2021
print(data)
#检验成功 我们截取了1973年1月到2021年12月的数据
DatetimeIndex(['1973-01-01', '1973-02-01', '1973-03-01', '1973-04-01',
'1973-05-01', '1973-06-01', '1973-07-01', '1973-08-01',
'1973-09-01', '1973-10-01',
...
'2021-03-01', '2021-04-01', '2021-05-01', '2021-06-01',
'2021-07-01', '2021-08-01', '2021-09-01', '2021-10-01',
'2021-11-01', '2021-12-01'],
dtype='datetime64[ns]', name='Month', length=588, freq=None)
Emission_1 Emission_2 Emission_3 Emission_4 Emission_5 \
Month
1973-01-01 73.112 12.163 2.397 0.128 23.698
1973-02-01 65.369 11.696 2.080 0.106 19.886
1973-03-01 65.006 13.979 1.181 0.083 18.851
1973-04-01 61.717 14.612 1.031 0.130 15.783
1973-05-01 62.686 17.326 0.957 0.167 16.920
... ... ... ... ... ...
2021-08-01 101.985 68.814 0.346 0.875 0.510
2021-09-01 80.443 54.323 0.251 0.732 0.416
2021-10-01 64.520 51.508 0.286 0.722 0.346
2021-11-01 59.326 48.144 0.301 0.919 0.324
2021-12-01 62.391 48.368 0.339 0.696 0.342
Emission_6 Emission_7 Emission_8 Emission_9
Month
1973-01-01 26.223 0.032 0.896 111.499
1973-02-01 22.071 0.032 0.896 99.137
1973-03-01 20.115 0.032 0.896 99.100
1973-04-01 16.944 0.032 0.896 93.273
1973-05-01 18.044 0.032 0.896 98.057
... ... ... ... ...
2021-08-01 1.731 0.040 0.892 173.463
2021-09-01 1.400 0.039 0.863 137.067
2021-10-01 1.353 0.040 0.892 118.313
2021-11-01 1.544 0.039 0.863 109.915
2021-12-01 1.377 0.040 0.892 113.068
[588 rows x 9 columns]
数据平稳性检验
这里我们借2.0.4趋势分解得到的图进行分析：
结果我们可以看到时间序列的趋势、季节和随机效应分别从数据中分离出来。煤电力部门二氧化碳排放在一定时间前上升，然后下降。且具有周期长度为12个月的季节变动。观察序列的随机效应，可以发现提取出趋势和季节效应后的残差基本稳定。常用的对时间序列进行平稳性检验的方法有：
绘制移动平均值和标准差
检验时间序列平稳性的一种正式的方法是绘制移动平均值和标准差，通过观察时间序列的平均值和标准差是否随时间变化来判断序列是否平稳。在Python中定义函数TestStationaryPlot来检验序列的平稳性。
def TestStationaryPlot(ts, figsize=(30, 18), linewidth=2):
rol_mean = ts.rolling(window=12, center=False).mean()
rol_std = ts.rolling(window=12, center=False).std()
plt.figure(figsize=figsize) # 设置图的大小
plt.plot(ts, color='blue', label=u'原始数据')
plt.plot(rol_mean, color='red', linestyle='-.', label=u'移动平均')
plt.plot(rol_std, color='black', linestyle='--', label=u'标准差')
plt.xticks(fontsize=25)
plt.yticks(fontsize=25)
plt.xlabel(u'时间（年）', fontsize=25)
plt.ylabel(u'总排放量', fontsize=25)
plt.legend(loc='best', fontsize=18)
plt.title(u'移动平均和标准差', fontsize=27)
plt.show(block=True)
在这里，window=12 是指滚动窗口的大小，用于计算移动平均和移动标准差。滚动窗口是一个在时间序列数据上滑动的窗口，其大小由window参数指定。在这个情况下，窗口大小为 12，即每次计算移动平均和移动标准差时，都使用过去的 12 个数据点。
具体来说，对于移动平均，窗口内的数据点的平均值被计算，然后窗口向前滑动一个时间步长，再计算下一组数据点的平均值，依此类推。
同样，对于移动标准差，窗口内的数据点的标准差被计算，然后窗口向前滑动一个时间步长，再计算下一组数据点的标准差。
这种技术常用于平滑时间序列数据，以便更好地观察数据的趋势和周期性。在这里，窗口大小为 12 意味着每个计算点都基于前面 12 个数据点的信息。
简要来说，这里的滚动窗口设置12个点，意思是设置了通过过去的12个月的数据，预测未来
ADF检验
检验序列平稳性的第二种方法是使用统计学中的ADF检验。检验的零假设是时间序列是非平稳的。检验结果比较了不同置信水平下的检验统计量和临界值。如果“检验统计量”小于“临界值”，我们可以拒绝零假设，认为在给定的显著性水平下这个序列是平稳的。在Python中定义TestStationaryAdfuller函数，对序列进行ADF检验。
def TestStationaryAdfuller(ts, cutoff=0.01):
ts_test = adfuller(ts, autolag='AIC')
ts_test_output = pd.Series(ts_test[0:4],
index=[
'Test Statistic', 'p-value', '#Lags Used',
'Number of Observations Used'
])
for key, value in ts_test[4].items():
ts_test_output['Critical Value (%s)' % key] = value
print(ts_test_output)
if ts_test[1] <= cutoff:
print(u"拒绝原假设，即数据没有单位根，序列是平稳的。")
else:
print(u"不能拒绝原假设，即数据存在单位根，数据是非平稳序列。")
TestStationaryPlot(data['Emission_1'])
TestStationaryAdfuller(data['Emission_1'])
Test Statistic -0.937859
p-value 0.775249
#Lags Used 18.000000
Number of Observations Used 569.000000
Critical Value (1%) -3.441895
Critical Value (5%) -2.866633
Critical Value (10%) -2.569482
dtype: float64
不能拒绝原假设，即数据存在单位根，数据是非平稳序列。
TestStationaryPlot(data['Emission_2'])
TestStationaryAdfuller(data['Emission_2'])
Test Statistic 1.214420
p-value 0.996084
#Lags Used 18.000000
Number of Observations Used 569.000000
Critical Value (1%) -3.441895
Critical Value (5%) -2.866633
Critical Value (10%) -2.569482
dtype: float64
不能拒绝原假设，即数据存在单位根，数据是非平稳序列。
TestStationaryPlot(data['Emission_3'])
TestStationaryAdfuller(data['Emission_3'])
Test Statistic -2.742302
p-value 0.067000
#Lags Used 16.000000
Number of Observations Used 571.000000
Critical Value (1%) -3.441854
Critical Value (5%) -2.866615
Critical Value (10%) -2.569473
dtype: float64
不能拒绝原假设，即数据存在单位根，数据是非平稳序列。
TestStationaryPlot(data['Emission_4'])
TestStationaryAdfuller(data['Emission_4'])
Test Statistic -1.350412
p-value 0.605805
#Lags Used 19.000000
Number of Observations Used 568.000000
Critical Value (1%) -3.441915
Critical Value (5%) -2.866642
Critical Value (10%) -2.569487
dtype: float64
不能拒绝原假设，即数据存在单位根，数据是非平稳序列。
TestStationaryPlot(data['Emission_5'])
TestStationaryAdfuller(data['Emission_5'])
Test Statistic -1.560655
p-value 0.503322
#Lags Used 15.000000
Number of Observations Used 572.000000
Critical Value (1%) -3.441834
Critical Value (5%) -2.866606
Critical Value (10%) -2.569468
dtype: float64
不能拒绝原假设，即数据存在单位根，数据是非平稳序列。
TestStationaryPlot(data['Emission_6'])
TestStationaryAdfuller(data['Emission_6'])
Test Statistic -1.702096
p-value 0.430033
#Lags Used 14.000000
Number of Observations Used 573.000000
Critical Value (1%) -3.441814
Critical Value (5%) -2.866597
Critical Value (10%) -2.569463
dtype: float64
不能拒绝原假设，即数据存在单位根，数据是非平稳序列。
TestStationaryPlot(data['Emission_7'])
TestStationaryAdfuller(data['Emission_7'])
Test Statistic -0.488247
p-value 0.894290
#Lags Used 19.000000
Number of Observations Used 568.000000
Critical Value (1%) -3.441915
Critical Value (5%) -2.866642
Critical Value (10%) -2.569487
dtype: float64
不能拒绝原假设，即数据存在单位根，数据是非平稳序列。
TestStationaryPlot(data['Emission_8'])
TestStationaryAdfuller(data['Emission_8'])
Test Statistic -2.584044
p-value 0.096340
#Lags Used 13.000000
Number of Observations Used 574.000000
Critical Value (1%) -3.441794
Critical Value (5%) -2.866588
Critical Value (10%) -2.569459
dtype: float64
不能拒绝原假设，即数据存在单位根，数据是非平稳序列。
TestStationaryPlot(data['Emission_9'])
TestStationaryAdfuller(data['Emission_9'])
Test Statistic -1.535025
p-value 0.516119
#Lags Used 14.000000
Number of Observations Used 573.000000
Critical Value (1%) -3.441814
Critical Value (5%) -2.866597
Critical Value (10%) -2.569463
dtype: float64
不能拒绝原假设，即数据存在单位根，数据是非平稳序列。
由此可见，所有数据列是非平稳序列，需要进行序列平稳化
序列平稳化
消除时间序列的趋势效应和季节性效应最常用方法之一是差分法。对原始的序列作1阶12步差分来提取原序列的趋势效应和季节效应，差分后的序列图及检验结果如下所示。
消除趋势和季节性，进行一阶差分和季节性一阶差分
#消除消除趋势和季节性
# 一阶差分
data_first_difference = data - data.shift(1)
data_first_difference.dropna(inplace=True) # 删除包含NaN的行
# 季节性一阶差分（假设月份为季节性周期）
data_seasonal_first_difference = data.diff(12)
data_seasonal_first_difference.dropna(inplace=True) # 删除包含NaN的行
TestStationaryPlot(data_seasonal_first_difference['Emission_1'])
TestStationaryAdfuller(data_seasonal_first_difference['Emission_1'])
Test Statistic -5.472154
p-value 0.000002
#Lags Used 12.000000
Number of Observations Used 563.000000
Critical Value (1%) -3.442019
Critical Value (5%) -2.866687
Critical Value (10%) -2.569511
dtype: float64
拒绝原假设，即数据没有单位根，序列是平稳的。
TestStationaryPlot(data_seasonal_first_difference['Emission_2'])
TestStationaryAdfuller(data_seasonal_first_difference['Emission_2'])
Test Statistic -6.651603e+00
p-value 5.104387e-09
#Lags Used 1.200000e+01
Number of Observations Used 5.630000e+02
Critical Value (1%) -3.442019e+00
Critica
最低0.47元/天 解锁文章
扫一扫
591
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
