SARIMAXモデルを用いた電力価格へ影響するパラメータの考察


はじめに

私はインフラ系企業に勤めていて、独学でPythonやTableauを学んでいました。
社内データの活用の際、検討段階でデータ分析ができたり、
何ができるかなどがわかることでより良いアウトプットを出すことできるのでは?
と考えて、ここ数か月勉強をしています。
今回は、卸電力価格が何に影響を受けているか考えてみたいと思います。

概要

今回行うこと

九州エリアの前日市場価格に影響を及ぼす因子に関する分析を行います。
当該エリアは、全国の中でも一日における価格のギャップが大きいエリアです。
理由は、再エネが豊富だからです。
今後、ほかのエリアも当該エリアのようになると予想されていて、
興味を持っているので、考えてみたいと思いました。

使うデータ

九州電力エリアの前日前市場価格
出典:https://www.jepx.jp/electricpower/market-data/spot/

JEPXのウェブサイト

需要実績(需要全体、水力、原子力、隣接エリアへ流入する電力量などのデータ)
出典:https://www.kyuden.co.jp/td_area_jukyu/jukyu.html

九州電力が公開している需給実績

手順

  • データの前処理

  • SARIMAXを用いて分析を行い、
    卸電力価格に影響するパラメータは何かを見る。

モデル実装

実行環境について

  • Google Chrome

  • Google Colaboratory

  • Python3

  • Windows11

前処理

今回は2023年4月1日~2024年3月31日までのデータを集めてました。
それぞれ、①電力価格のデータは30分毎、②需要データは1時間毎のデータです。
今回は、日最高値について考えるので、それぞれの日の最大の数字のみをINPUTデータにします。

まず、使用するライブラリのインポートを行います。

import pandas as pd
import os
from datetime import datetime
import matplotlib.pyplot as plt
from matplotlib.pylab import rcParams
%matplotlib inline
from sklearn.preprocessing import MinMaxScaler
import statsmodels.api as sm
from statsmodels.tsa.seasonal import seasonal_decompose
from statsmodels.tsa.api import SARIMAX, AutoReg
from statsmodels.tsa.arima.model import ARIMA

次に、ドライブ上に保存したデータを作成した下記のようなディレクトリに保存していきます。

#Google Driveをマウントする。
from google.colab import drive
drive.mount('/content/drive')

#データを格納するフォルダを作成する
kakaku_dir = '/content/drive/My Drive/input/kakaku'
jyuyou_dir = '/content/drive/My Drive/input/jyuyou'
process_dir = '/content/drive/My Drive/process'
os.makedirs(kakaku_dir)
os.makedirs(jyuyou_dir)

#各データを保存する。
kakaku_path = os.path.join('/content/drive/My Drive', 'spot_summary_2023.csv')
jyuyou_path1 = os.path.join('/content/drive/My Drive','area_jyukyu_jisseki_2023_1Q.csv')
jyuyou_path2 = os.path.join('/content/drive/My Drive','area_jyukyu_jisseki_2023_2Q.csv')
jyuyou_path3 = os.path.join('/content/drive/My Drive','area_jyukyu_jisseki_2023_3Q.csv')
jyuyou_path4 = os.path.join('/content/drive/My Drive','area_jyukyu_jisseki_2023_4Q.csv')

kakaku_data = pd.read_csv(kakaku_path, encoding='shift-jis')
jyuyou_data1 = pd.read_csv(jyuyou_path1, encoding='shift-jis')
jyuyou_data2 = pd.read_csv(jyuyou_path2, encoding='shift-jis')
jyuyou_data3 = pd.read_csv(jyuyou_path3, encoding='shift-jis')
jyuyou_data4 = pd.read_csv(jyuyou_path4, encoding='shift-jis')

kakaku_path_rev = os.path.join(kakaku_dir, 'spot_summary_2023.csv')
jyuyou_data1_rev = os.path.join(jyuyou_dir,'area_jyukyu_jisseki_2023_1Q.csv')
jyuyou_data2_rev = os.path.join(jyuyou_dir,'area_jyukyu_jisseki_2023_2Q.csv')
jyuyou_data3_rev = os.path.join(jyuyou_dir,'area_jyukyu_jisseki_2023_3Q.csv')
jyuyou_data4_rev = os.path.join(jyuyou_dir,'area_jyukyu_jisseki_2023_4Q.csv')

kakaku_data.to_csv(kakaku_path_rev, index=False)
jyuyou_data1.to_csv(jyuyou_data1_rev, index=False)
jyuyou_data2.to_csv(jyuyou_data2_rev, index=False)
jyuyou_data3.to_csv(jyuyou_data3_rev, index=False)
jyuyou_data4.to_csv(jyuyou_data4_rev, index=False)

需要データを下記のようにきれいに整形します。

'jyuyou'という名前のディレクトリに入れていて、
九州電力のウェブサイトからとってきたファイルは4半期ごとに分かれています。
→1つに結合する

また、1,2行目にカラム名が分かれています。
→1行にします

’DATE_TIME’というカラムに入っているのは、'yyyy/mm/dd hh:mm'です。
→'yyyy/mm/dd’に変えます。

さいごに日付ごとの最大値のみを抽出して、csvファイルに保存します。

#'input'フォルダ内のすべてのCSVファイルを取得
csv_files = [f for f in os.listdir(jyuyou_dir) if f.endswith('.csv')]

# CSVファイルを読み込み、結合するためのリストを作成
jyuyou = []
for csv_file in csv_files:
    file_path = os.path.join(jyuyou_dir, csv_file)
    df = pd.read_csv(file_path, encoding='shift-jis', header=None)
    # 1行目と2行目を読み込む
    first_row = df.iloc[0]
    second_row = df.iloc[1]

    # 1行目が空欄だったら2行目のセルのみを返し、それ以外は1行目と2行目を結合
    new_columns = [second if pd.isna(first) or first == '' else f'{first}_{second}' for first, second in zip(first_row, second_row)]
    
    # 新しいカラム名を設定
    df.columns = new_columns
    # 1行目と2行目を削除
    df = df.drop([0, 1]).reset_index(drop=True)
    jyuyou.append(df)

daily_jyuyou = pd.concat(jyuyou, ignore_index=True)
daily_jyuyou['DATE_TIME'] = pd.to_datetime(daily_jyuyou['DATE_TIME']).dt.strftime('%Y/%m/%d')
daily_jyuyou = daily_jyuyou.groupby('DATE_TIME').max()
daily_jyuyou_file = os.path.join(process_dir, 'daily_jyuyou.csv')
daily_jyuyou.to_csv(daily_jyuyou_file)

次に価格データについてです。
こちらは九州エリア以外の価格などいろいろ入っているので、
必要なもののみを残します。

#価格データ読み込み
#必要なものを残す
kakaku = pd.read_csv(kakaku_dir, encoding='shift-jis')
kakaku = kakaku[['受渡日', '売り入札量(kWh)', '買い入札量(kWh)', 'エリアプライス九州(円/kWh)',]]
kakaku = kakaku.rename(columns={'受渡日': 'DATE_TIME'})

#MAX,MIN,MEANのカラムを作る。
daily_max = kakaku.groupby('DATE_TIME').max()

output = os.path.join(process_dir, 'combined.csv')
daily_max.to_csv(output, index=False, encoding='shift-jis')

最後に需要と価格データを結合して、
INPUTデータが完成しました。

dataset = pd.merge(daily_max,daily_jyuyou, on='DATE_TIME')

データの確認

いったん価格データの1年間の推移をみてみます。
夏や冬において日最高値が跳ねる傾向が見られます。

rcParams['figure.figsize'] = 15, 6
y = dataset["エリアプライス九州(円/kWh)"].astype("float")
y.plot()
九州エリアの日最高値

自己相関係数と偏自己相関係数についてもみてみます。

fig,ax = plt.subplots(2,1,figsize=(12,8))
fig = sm.graphics.tsa.plot_acf(y, lags=100, ax=ax[0], color="darkgoldenrod")
fig = sm.graphics.tsa.plot_pacf(y, lags=100, ax=ax[1], color="darkgoldenrod")
plt.show()
自己相関係数と偏自己相関係数のコレログラム

上図(自己相関係数)を見ると、規則的な山が見られます。
下図(偏自己相関係数)を見ると、過去データの影響を取り除くことができていると考えられます。

SARIMAXモデル実装

SARIMAXモデルを実装します。
パラメータに関しては最適なものを見つけてきます。

#データの用意
X = dataset.drop("エリアプライス九州(円/kWh)",axis=1).astype("float")
y = dataset["エリアプライス九州(円/kWh)"].astype("float")

#  orderの最適化関数
def selectparameter(DATA, s):
    p = d = q = range(0, 2)
    pdq = list(itertools.product(p, d, q))
    seasonal_pdq = [(x[0], x[1], x[2], s) for x in list(itertools.product(p, d, q))]
    parameters = []
    BICs = np.array([])
    for param in pdq:
        for param_seasonal in seasonal_pdq:
            try:
                mod = sm.tsa.statespace.SARIMAX(DATA,
                                                order=param,
                                                seasonal_order=param_seasonal)
                results = mod.fit()
                parameters.append([param, param_seasonal, results.bic])
                BICs = np.append(BICs, results.bic)
            except:
                continue
    return parameters[np.argmin(BICs)]


#モデルの構築
best_params = selectparameter(y, 12)
SARIMAX_price = sm.tsa.statespace.SARIMAX(y, exog=X, order=best_params[0], seasonal_order=best_params[1], trend="c").fit()
print(SARIMAX_price.summary())

出力結果は下記です。
coefの数字についてみてみると、風力‗抑制量[MWh]が0.0114と、最も大きく寄与していそうです。

 ================================================================================
                   coef    std err          z      P>|z|      [0.025      0.975]
--------------------------------------------------------------------------------
intercept        0.0107   1.96e-05    545.131      0.000       0.011       0.011
売り入札量(kWh)   -1.927e-07   5.48e-08     -3.513      0.000      -3e-07   -8.52e-08
買い入札量(kWh)    4.253e-07   1.44e-07      2.961      0.003    1.44e-07    7.07e-07
エリア需要〔MWh〕       0.0012      0.000      4.752      0.000       0.001       0.002
原子力〔MWh〕         0.0010      0.001      1.068      0.285      -0.001       0.003
火力〔MWh〕          0.0002      0.000      0.465      0.642      -0.001       0.001
水力〔MWh〕         -0.0019      0.003     -0.752      0.452      -0.007       0.003
地熱〔MWh〕         -0.0149      0.016     -0.951      0.342      -0.046       0.016
バイオマス〔MWh〕      -0.0049      0.008     -0.620      0.536      -0.020       0.011
太陽光_実績〔MWh〕  -7.487e-05   7.14e-05     -1.049      0.294      -0.000     6.5e-05
太陽光_抑制量〔MWh〕     0.0003      0.000      1.011      0.312      -0.000       0.001
風力_実績〔MWh〕       0.0029      0.003      0.929      0.353      -0.003       0.009
風力_抑制量〔MWh〕      0.0114      0.004      3.028      0.002       0.004       0.019
揚水等〔MWh〕         0.0016      0.000      3.575      0.000       0.001       0.003
連系線〔MWh〕        -0.0006      0.000     -2.218      0.027      -0.001    -7.5e-05

最後に

今回、日最高値に関して考えましたが、
日最安値については、日最高値と異なるパラメータが寄与すると思われます。
日最高値は夕方につけることが多いのに対して、
日最安値は昼につけることが多いです。
なので、太陽光などのパラメータによる寄与が大きいと予想されます。


この記事が気に入ったらサポートをしてみませんか?