社会科学市场博弈和价格预测之时间序列挖掘（Datawhale AI 夏令营）

深入理解赛题——探索性数据分析

        首先，我们先介绍一下什么是EDA： 探索性数据分析（Exploratory Data Analysis, EDA）是一组数据分析技术，旨在总结其主要特征，通常通过可视化手段来实现。EDA 的目标是通过数据的统计摘要和图形展示来发现数据的结构、异常值、模式、趋势、关系以及变量之间的相互作用。

为什么进行EDA？

        在现在的数据挖掘类比赛中，模型和方法选择空间往往很小，同时存在不少自动机器学习框架（如AutoGulon、AutoSKLearn）会基于一定规则，自动构造特征，采用尽可能多的模型组合来获得好分数。因此最后的关键涨分点落在了对数据的理解上，并由此构造的强特征（对结果有关键影响的变量）。

        因此在拿到一个数据集后，需要对数据做尽可能多的探索，了解数据所在的领域先验知识，数据本身的特性等，并总结为一系列有用的信息。下面，就让我们对本次比赛的数据做进一步的探索吧(ง •̀_•́)ง。

数据基本信息

import pandas as pd
import seaborn as sns
import matplotlib.pylab as plt
from pathlib import Path
import warnings
warnings.filterwarnings('ignore')
plt.style.use('ggplot')
plt.rcParams['font.sans-serif'] = ['WenQuanYi Micro Hei', "SimHei"]
plt.rcParams['axes.unicode_minus'] = False

base_path = Path("data")

# 读取数据
electricity_price = pd.read_csv(base_path / "electricity_price_parsed.csv", parse_dates=["timestamp"], index_col=0)
electricity_price.columns = ["demand", "price"]
electricity_price.head()

        通过数据基本信息的读取，我们可以确认数据的时间范围和格式。

# 创建测试集掩码，标记出所有价格为 NaN 的数据行
test_mask = electricity_price["price"].isna()

# 创建训练集掩码，标记出所有价格不为 NaN 的数据行,其中~代表布尔取反，即True和False互换。 `[True,False]`和 `~[False,True]`一致
train_mask = ~test_mask

# 打印训练集的范围和总长度
print(f"训练集范围：{electricity_price[train_mask].index.min()} --> {electricity_price[train_mask].index.max()}\t总长度{len(electricity_price[train_mask])}")

# 打印测试集的范围和总长度
print(f"测试集范围：{electricity_price[test_mask].index.min()} --> {electricity_price[test_mask].index.max()}\t总长度{len(electricity_price[test_mask])}")

        确认是否有缺失值，可以发现除了测试集部分，其他均无缺失。

electricity_price.info()

数据观察

        进一步观察数据，发现数据中的出清价格在一小时内都是相同的，但负荷会在一小时中变动。

electricity_price.head(12)

        定义一个函数来确认这一点：

def check_repeated(data, repeat_count=4):
    for i in range(0, len(data), repeat_count):
        subsequence = data[i:i + repeat_count]
        if len(set(subsequence)) != 1:
            print(f"序列数据不是元素不断重复 {repeat_count} 次")
            return
    print(f"序列数据是元素不断重复 {repeat_count} 次")

check_repeated(electricity_price[train_mask]["price"])

        从图中也能确认，在一小时中出清价格一致。

electricity_price.loc["2023-01-03"].plot(y="price", figsize=(18, 5), marker='o')
plt.title("2023年1月3日出清价格走势")

统计指标分析

        首先我们从最基本的单变量数值关系入手，探究其统计指标的表现。

electricity_price.describe()

        总负荷和数值分布：

ax = sns.displot(
    electricity_price,
    x="demand",
    aspect=1.5,
    height=5,
    kde=True
)
ax.set(title="总负荷数值分布")

        出清价格数值分布：

ax = sns.displot(electricity_price, x="price", aspect=1.5, height=5, kde=True)
ax.set(title="出清价格数值分布")

分时统计特征

        不同小时的总电力需求和电价：

electricity_price["hour"] = electricity_price.index.hour
electricity_price["month"] = electricity_price.index.month
electricity_price["day"] = electricity_price.index.day
electricity_price["weekday"] = electricity_price.index.weekday
electricity_price["year"] = electricity_price.index.year

fig, ax = plt.subplots(1, 2, figsize=(20, 9), constrained_layout=True)

sns.lineplot(
    electricity_price.groupby(["hour", "year"])["price"].mean().reset_index(),
    x="hour",
    y="price",
    ax=ax[0],
    marker="o",
    hue="year",
    palette="Set1"
)
ax[0].set_title("不同年份的分时电价")
ax[0].set_xticks(range(24))

sns.lineplot(
    electricity_price.groupby(["hour", "year"])["demand"].mean().reset_index(),
    x="hour",
    y="demand",
    ax=ax[1],
    marker="o",
    hue="year",
    palette="Set1"
)
ax[1].set_title("不同年份的分时总负荷")
ax[1].set_xticks(range(24))

透视表（Pivot Table）

pivot = pd.pivot_table(
    electricity_price,
    values="price",
    index="month",
    columns="hour"
)
pivot = pivot.astype(int)
pivot

plt.figure(figsize=(17, 10))
sns.heatmap(
    pivot,
    cmap="coolwarm",
    linewidths=0.5,
    annot=True,
    fmt=".0f",
    annot_kws={"size": 12, "weight": "bold", "color": "black"}
)
plt.title("不同月份和时间下的电价")

负电价

minus_mask = electricity_price["price"] < 0

plt.figure(figsize=(12, 7))
ax = sns.countplot(electricity_price[minus_mask], x="hour")
ax.set(title="负电量频数出现的小时分布")

plt.figure(figsize=(12, 7))
ax = sns.countplot(electricity_price[minus_mask], x="month")
ax.set(title="负电量频数出现的月分布")

(
    electricity_price[minus_mask]
    .groupby(["month", "day"])["price"]
    .size()
    .reset_index()
    .sort_values("price", ascending=False)
    .head(15)
)

典型的负电价趋势：

plt.figure(figsize=(20, 8))
ax = sns.lineplot(electricity_price.loc["2022-05-01":"2022-05-09"]["price"], color="black")

for i in range(1, 10):
    ax.axvspan(f"2022-05-0{i} 10:00:00", f"2022-05-0{i} 15:00:00", color='yellow', alpha=0.2)

plt.axhline(y=0, color="red", linestyle="--")
plt.setp(ax.get_xticklabels(), rotation=45, ha='right')
plt.title("一段比较典型的负电价趋势：2022年5月1日 - 2022年5月9日")

假期对负电价的影响：

plt.figure(figsize=(20, 8))
ax = sns.lineplot(electricity_price.loc["2023-01-21":"2023-02-02"]["price"], color="black")
plt.axhline(y=0, color="red", linestyle="--")
plt.setp(ax.get_xticklabels(), rotation=45, ha='right')
plt.title("假期对负电价的影响 2022年1月28日 - 2022年2月6日")

高电价

upper_threshold = (
    electricity_price["price"].mean() +
    3 * electricity_price["price"].std()
)
high_abnormal_mask = (
    electricity_price["price"] > upper_threshold
)

plt.figure(figsize=(12, 7))
ax = sns.countplot(electricity_price[high_abnormal_mask], x="hour")
ax.set(title="高电量频数出现的小时分布")

plt.figure(figsize=(12, 7))
ax = sns.countplot(electricity_price[high_abnormal_mask], x="month")
ax.set(title="高电量频数出现的月分布")

(
    electricity_price[high_abnormal_mask]
    .groupby(["month", "day"])["price"]
    .size()
    .reset_index()
    .sort_values("price", ascending=False)
    .head(15)
)

典型高电价时间段：

plt.figure(figsize=(20, 8))
ax = sns.lineplot(electricity_price.loc["2022-08-01":"2022-08-08"]["demand"], color="black")
plt.setp(ax.get_xticklabels(), rotation=45, ha='right')
plt.axvspan("2022-08-03", "2022-08-06", color="yellow", alpha=0.2)
plt.title("2022年8月1日 - 2022年8月8日 总需求")

plt.figure(figsize=(20, 8))
ax = sns.lineplot(electricity_price.loc["2022-08-01":"2022-08-08"]["price"], color="black")
plt.axhline(y=upper_threshold, color="red", linestyle="--")
plt.setp(ax.get_xticklabels(), rotation=45, ha='right')
plt.title("典型高电价时间段 2022年8月1日 - 2022年8月8日 中的3-6日")

双变量分析

electricity_price[["demand", "price"]].corr()

plt.figure(figsize=(8, 10))
sns.regplot(
    data=electricity_price.loc["2022"],
    x="demand",
    y="price",
    scatter_kws={"s": 0.5, "alpha": 0.6, "color": "black"},
    color="red",
    lowess=True
)

时间序列模型构建

        在完成了以上的探索性数据分析之后，我们对数据有了较为深入的了解。接下来，我们将构建一个时间序列模型来进行电价预测。

构建时间序列模型

        时间序列模型用于处理和分析按时间顺序排列的数据点。在本次任务中，我们将使用一些常见的时间序列建模方法，如ARIMA（自回归积分滑动平均模型）和Prophet模型。

        首先，我们需要将数据整理为时间序列格式，并划分训练集和测试集。

      

from statsmodels.tsa.arima.model import ARIMA
from sklearn.metrics import mean_squared_error
import numpy as np

# 准备数据
train_data = electricity_price[train_mask]["price"]
test_data = electricity_price[test_mask]["price"]

# 构建并训练ARIMA模型
model = ARIMA(train_data, order=(5,1,0)) # (p,d,q)分别表示自回归阶数、差分阶数、移动平均阶数
model_fit = model.fit(disp=0)

# 预测
predictions = model_fit.forecast(steps=len(test_data))[0]

# 计算均方误差
mse = mean_squared_error(test_data, predictions)
print(f'Mean Squared Error: {mse}')

# 绘制预测结果
plt.figure(figsize=(12, 6))
plt.plot(test_data.index, test_data, label='Actual Price')
plt.plot(test_data.index, predictions, color='red', label='Predicted Price')
plt.title('Electricity Price Prediction')
plt.xlabel('Time')
plt.ylabel('Price')
plt.legend()
plt.show()

        在上述代码中，我们使用了ARIMA模型对电价进行了预测，并计算了预测结果的均方误差。同时，我们将实际价格和预测价格进行了可视化展示。

        除了ARIMA模型，我们还可以尝试使用Prophet模型进行预测。

          通过以上步骤，我们完成了对数据的探索性分析，发现了很多有趣的现象和规律，为后续的模型构建和特征提取提供了重要的参考。

from fbprophet import Prophet

# 准备数据
train_data = electricity_price[train_mask].reset_index()[["timestamp", "price"]]
train_data.columns = ["ds", "y"]

# 构建并训练Prophet模型
model = Prophet()
model.fit(train_data)

# 生成未来的日期数据
future = model.make_future_dataframe(periods=len(test_data), freq='15T')

# 进行预测
forecast = model.predict(future)

# 绘制预测结果
model.plot(forecast)
plt.title('Electricity Price Prediction with Prophet')
plt.xlabel('Time')
plt.ylabel('Price')
plt.show()

        通过上述代码，我们可以使用Prophet模型对电价进行预测，并将预测结果进行可视化展示。

玻尔兹曼财富模型

        接下来，我们将探讨另一种方法，即使用Agent-Based Modeling（ABM）来模拟市场中的博弈行为。这里我们介绍一种经典的模型——玻尔兹曼财富模型。

玻尔兹曼财富模型

import numpy as np
import pandas as pd
import seaborn as sns
from tqdm import tqdm
import matplotlib.pylab as plt
plt.style.use('ggplot')
plt.rcParams['font.sans-serif'] = ['SimHei']
plt.rcParams['axes.unicode_minus'] = False

from mesa.model import Model
from mesa.time import RandomActivation
from mesa.agent import Agent
from mesa.datacollection import DataCollector

def compute_gini(model):
    agent_wealths = [agent.wealth for agent in model.schedule.agents]
    x = sorted(agent_wealths)
    N = len(model.agents)
    B = sum(xi * (N - i) for i, xi in enumerate(x)) / (N * sum(x))
    return 1 + (1 / N) - 2 * B

class WealthModel(Model):
    def __init__(self):
        super().__init__()
        self.schedule = RandomActivation(self)
        self.datacollector = DataCollector(model_reporters={"gini": compute_gini}, agent_reporters={"Wealth": "wealth"})

    def step(self):
        self.datacollector.collect(self)
        self.schedule.step()

    def add_agents(self, agents: list[Agent]):
        for agent in agents:
            self.schedule.add(agent)

class WealthAgent(Agent):
    def __init__(self, unique_id, model):
        super().__init__(unique_id, model)
        self.wealth = 1

    def step(self):
        if self.wealth > 0:
            other_agent = self.random.choice(self.model.schedule.agents)
            other_agent.wealth += 1
            self.wealth -= 1

steps = 200
env = WealthModel()
agents = [WealthAgent(i, env) for i in range(1000)]
env.add_agents(agents)

for i in tqdm(range(steps)):
    env.step()

df = env.datacollector.get_agent_vars_dataframe()
df

last_wealth = df.xs(steps-1, level="Step")
ax = sns.histplot(last_wealth, discrete=True)
ax.set(title="最终的财富分布", xlabel="财富量")

model_df = env.datacollector.get_model_vars_dataframe()
plt.figure(figsize=(12, 7))
ax = sns.lineplot(model_df)
ax.set(title='基尼系数变化', ylabel="基尼系数", xlabel="模拟时间步")

        通过上述代码，我们可以模拟玻尔兹曼财富模型，并观察不同个体的财富分布和基尼系数的变化。

富二代模型

        在玻尔兹曼财富模型的基础上，我们可以进一步模拟不同社会阶层之间的财富转移，例如模拟“富二代”的情况。

class WealthAgent(Agent):
    def __init__(self, unique_id, model):
        super().__init__(unique_id, model)
        self.wealth = 50

    def step(self):
        if self.wealth > 0:
            other_agent = self.random.choice(self.model.schedule.agents)
            if type(other_agent) == RichSecondAgent:
                if self.wealth >= 2:
                    self.wealth -= 2
                    other_agent.wealth += 2
                    return
            other_agent.wealth += 1
            self.wealth -= 1

class RichSecondAgent(Agent):
    def __init__(self, unique_id, model):
        super().__init__(unique_id, model)
        self.wealth = 250

    def step(self):
        if self.wealth > 0:
            other_agent = self.random.choice(self.model.schedule.agents)
            other_agent.wealth += 2
            self.wealth -= 2

steps = 5000
env = WealthModel()
agents = [WealthAgent(i, env) for i in range(280)]
agents += [RichSecondAgent(i, env) for i in range(280, 300)]
env.add_agents(agents)

for i in tqdm(range(steps)):
    env.step()

agent_df = env.datacollector.get_agent_vars_dataframe()
wealths = agent_df[agent_df.index.get_level_values("AgentID").isin([2, 3, 4, 280, 285, 297])]
plt.figure(figsize=(20, 9))
sns.lineplot(wealths, hue="AgentID", x="Step", y="Wealth", palette="Set1")

        通过上述代码，我们可以模拟“富二代”在社会中的财富转移过程，并观察不同个体的财富变化。

结语

        通过对数据的探索性分析、时间序列模型的构建以及Agent-Based Modeling的模拟，我们可以更好地理解和预测电力市场的出清价格。探索性数据分析帮助我们发现了数据中的模式和规律，揭示了市场需求、出清价格在不同时间和条件下的变化趋势。时间序列模型为我们提供了预测未来价格的有效工具，而玻尔兹曼财富模型和富二代模型则通过模拟市场博弈行为，展示了不同主体在市场中的互动和其对整体市场价格的影响。

        在实际应用中，这些方法和模型可以相互补充，提供更全面和精准的市场价格预测。例如，通过结合时间序列分析和Agent-Based Modeling，我们不仅能够预测未来的市场出清价格，还能模拟不同市场参与者的行为和策略，评估其对市场的影响。此外，EDA过程中发现的关键特征和模式，为模型的构建和优化提供了重要的参考和指导。

        展望未来，随着更多外部数据（如天气、经济指标等）的引入，以及对市场规则和参与者行为的更深入研究，我们有望构建出更为精确和可靠的价格预测模型。这不仅有助于电力市场的稳定运行，也为其他类似市场的研究和应用提供了宝贵的经验和方法。

        总之，通过系统的探索和分析，我们不仅提升了对电力市场的理解，也为进一步的研究和应用奠定了坚实的基础。希望这些方法和经验能为未来的研究者和从业者提供有价值的参考和启示。

如果你觉得这篇博文对你有帮助，请点赞、收藏、关注我，并且可以打赏支持我！

欢迎关注我的后续博文，我将分享更多关于人工智能、自然语言处理和计算机视觉的精彩内容。

谢谢大家的支持！