FAO食料価格と原油価格のラグ相関をPythonで再現する（実装編）

🧭 概要

本記事では、前回の分析記事

👉 原油価格の変動は、何か月遅れて食料価格に波及するのか？

で示した相関構造・ラグ分析・回帰モデルを、Pythonで実際に再現する方法を紹介します。

使用データは以下の2系列：

期間：1990年〜2025年（共通月のみ抽出）

ラグ相関（Lag Correlation）の計算

import pandas as pd
import numpy as np
import matplotlib.pyplot as plt

# データ読み込み
fao = pd.read_csv("fao_food_price_monthly_clean.csv", parse_dates=["Date"])
energy = pd.read_csv("worldbank_energy_monthly_clean.csv", parse_dates=["Date"])

# 対象列を抽出
fao = fao[["Date", "Cereals"]]
energy = energy[["Date", "Crude oil, Brent"]]

# 結合
merged = pd.merge(fao, energy, on="Date", how="inner").dropna()
merged = merged.sort_values("Date")

# ラグ相関の算出
lags = range(-6, 7) # -6〜+6か月
corrs = []

for k in lags:
shifted = merged["Crude oil, Brent"].shift(k)
corr = merged["Cereals"].corr(shifted)
corrs.append(corr)

# 結果をDataFrame化
lag_df = pd.DataFrame({"Lag": lags, "Correlation": corrs})
best_lag = lag_df.loc[lag_df["Correlation"].idxmax()]

# 可視化
plt.figure(figsize=(8, 4))
plt.plot(lag_df["Lag"], lag_df["Correlation"], marker="o")
plt.axvline(best_lag["Lag"], color="red", linestyle="--", label=f"Peak Lag = {best_lag['Lag']} months")
plt.title("Lag Correlation: Brent Oil vs FAO Cereals Index")
plt.xlabel("Lag (months)")
plt.ylabel("Correlation coefficient")
plt.legend()
plt.grid(True)
plt.tight_layout()
plt.savefig("lag_correlation_brent_cereals.png")
plt.show()

print(f"最大相関: {best_lag['Correlation']:.3f}（ラグ {best_lag['Lag']} ヶ月）")

結果

最大相関：r = 0.879

ピークラグ：−1か月
→ 原油が約1か月先行して食料に影響。

単回帰モデルで確認（1か月遅れ効果）

import pandas as pd
import statsmodels.api as sm

# データ読み込み
fao = pd.read_csv("fao_food_price_monthly_clean.csv", parse_dates=["Date"])
energy = pd.read_csv("worldbank_energy_monthly_clean.csv", parse_dates=["Date"])

# 必要列を抽出
df = pd.merge(
fao[["Date", "Cereals"]],
energy[["Date", "Crude oil, Brent"]],
on="Date",
how="inner"
).dropna()

# 原油価格を1か月前にシフト
df["Brent_lag1"] = df["Crude oil, Brent"].shift(1)

# 欠損削除
df = df.dropna(subset=["Cereals", "Brent_lag1"])

# 回帰分析
X = sm.add_constant(df["Brent_lag1"]) # 定数項
y = df["Cereals"]

model = sm.OLS(y, X).fit()
print(model.summary())

結果

係数 β = 0.8468（p < 0.001）

R² = 0.756

→ 原油価格が1か月遅れで穀物価格を強く押し上げる傾向。

上昇／下落の非対称性を検証（NARDLの簡易版）

import pandas as pd
import statsmodels.api as sm

# データ読み込み
fao = pd.read_csv("fao_food_price_monthly_clean.csv", parse_dates=["Date"])
energy = pd.read_csv("worldbank_energy_monthly_clean.csv", parse_dates=["Date"])

# 結合
df = pd.merge(
fao[["Date", "Cereals"]],
energy[["Date", "Crude oil, Brent"]],
on="Date",
how="inner"
).dropna()

# 月次変化（ΔX, ΔY）
df["dCereals"] = df["Cereals"].diff()
df["dBrent"] = df["Crude oil, Brent"].diff()

# 1か月ラグを取る
df["dBrent_lag1"] = df["dBrent"].shift(1)

# 上昇分と下落分に分解
df["Brent_up"] = df["dBrent_lag1"].apply(lambda x: x if x > 0 else 0)
df["Brent_down"] = df["dBrent_lag1"].apply(lambda x: x if x < 0 else 0)

# 回帰
df = df.dropna()
X = sm.add_constant(df[["Brent_up", "Brent_down"]])
y = df["dCereals"]

model = sm.OLS(y, X).fit()
print(model.summary())

出力