03 SOIの移動平均と自己相関 | データサイエンス 2026

1. この回でやること

オリジナルの SOI には、短周期・長周期の変動が混在している。ここでは、簡単なフィルタ処理として移動平均を用い、比較的見やすい周期成分を取り出す。

121か月移動平均を計算する
元の SOI から 121か月移動平均を引き、10年以上の変動を除く
さらに 13か月移動平均をかけ、1年以内の変動を弱める
その時系列に対して自己相関係数を求める

ここで使う移動平均は、考え方を学ぶための簡単な方法である。実際の研究では、フーリエ変換やデジタルフィルタを使うことも多い。

2. Pythonでよく出る機能の意味

この回では、計算式そのものよりも、配列を作る・欠損を除く・条件で切り出す・関数を定義する といった操作が多く出てくる。ここが読めると、自己相関のコードもかなり理解しやすくなる。

書き方	意味
pd.read_csv(...)	表形式データを読み込む
dropna(subset=[...])	指定列に欠損がある行を落とす
reset_index(drop=True)	行番号を振り直す
df[(条件1) & (条件2)]	条件に合う行だけを取り出す
to_numpy()	pandas の列や表を NumPy 配列に変える
def 関数名(...):	自分で関数を定義する
pd.Series(x).rolling(...).mean()	移動平均を計算する
np.isnan(x)	欠損値 NaN かどうかを調べる
~np.isnan(x)	NaN ではないデータだけを選ぶためのマスクを作る
np.arange(a, b)	等間隔の整数列を作る
np.full(n, np.nan)	長さ n の NaN 配列を作る
np.asarray(x)	x を NumPy 配列として扱う
enumerate(...)	番号と中身を同時に取り出す

大事な見方
長いコードは、「何を計算したいか」と「そのためにデータをどんな形にしているか」に分けて読むと理解しやすい。

STEP 1

3. 準備とデータ読み込み

前回作成した SOI_timeseries_updated.txt を読み込む。自己相関は教科書に合わせて 1951–1997 年に制限してよい。

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

# データ読み込み
df = pd.read_csv("SOI_timeseries_updated.txt", sep=r"\s+")
df = df.dropna(subset=["SOI"]).reset_index(drop=True)

# 教科書に合わせるなら期間を制限
df = df[(df["YEAR"] >= 1951) & (df["YEAR"] <= 1997)].reset_index(drop=True)

time = df["TIME"].to_numpy()
soi = df["SOI"].to_numpy()

print(df.head())
print(len(soi))

このブロックで新しく出るもの

pd.read_csv：テキストファイルを表として読む
sep=r"\s+"：空白が1個以上あるところで区切る
dropna(subset=["SOI"])：SOI が欠損の行を落とす
reset_index(drop=True)：削除後に行番号を 0,1,2,... と振り直す
df[(条件)]：条件に合う行だけ残す
to_numpy()：列を NumPy 配列に変える

コードの意味

df は DataFrame で、表形式のデータそのもの
dropna をしているのは、欠損を含んだままだと後の計算や描画がややこしくなるから
df[(df["YEAR"] >= 1951) & (df["YEAR"] <= 1997)] は、「YEAR が 1951 以上かつ 1997 以下」の行だけを残すという意味
time と soi を NumPy 配列にしておくと、この後の数値計算がやりやすい

確認

TIME は何を表しているか
1951–1997 の期間にすると、データ数はいくつになるか
len(soi) は何を数えているか

STEP 2

4. 121か月移動平均

まず、121か月（約10年）移動平均を計算する。中央配置の移動平均を使う。

def moving_average(x, window):
    return pd.Series(x).rolling(window=window, center=True).mean().to_numpy()

soi_ma121 = moving_average(soi, ____)

このブロックで新しく出るもの

def moving_average(...)：自分で関数を作る
pd.Series(x)：1次元データを pandas の Series にする
rolling(window=...)：移動窓を作る
center=True：窓の中心に値を配置する
mean()：その窓の平均をとる
to_numpy()：最後に NumPy 配列へ戻す

コードの意味

moving_average は「配列 x と窓幅 window を受け取って、移動平均を返す関数」
関数にしておくと、121か月移動平均にも13か月移動平均にも同じ処理を再利用できる
center=True を付けないと、平均値が右寄せや左寄せになって時刻とずれてしまう
端のほうは十分な個数のデータが取れないため、移動平均の結果は NaN になる

考えてみよう

121 という数字は何を意味するか
center=True を付ける理由は何か
なぜ端の数か月分では移動平均が計算できないのか

plt.figure(figsize=(12,4))
plt.plot(time, soi, label="Original SOI", alpha=0.5)
plt.plot(time, soi_ma121, label="121-month moving average", linewidth=2)
plt.axhline(0, color="k", linewidth=0.8)
plt.xlabel("Year")
plt.ylabel("SOI")
plt.title("SOI and 121-month moving average")
plt.legend()
plt.grid(True, alpha=0.3)
plt.tight_layout()
plt.show()

描画でよく出るもの

plt.figure(figsize=(12,4))：図の大きさを決める
plt.plot(x, y)：折れ線グラフを描く
label=...：凡例に表示する名前
alpha=0.5：少し透明にする
linewidth=2：線を太くする
plt.axhline(0)：y=0 の横線を引く
plt.legend()：凡例を表示する
plt.grid(True)：グリッドを表示する
plt.tight_layout()：文字や図がはみ出さないように整える

STEP 3

5. 長周期成分の除去

次に、元の SOI から 121か月移動平均を引く。これにより、10年以上の長周期変動を取り除く。

soi_hp = soi - soi_ma121
mask_hp = ~np.isnan(soi_hp)

time_hp = time[mask_hp]
soi_hp_valid = soi_hp[mask_hp]

このブロックで新しく出るもの

np.isnan(x)：x が NaN かどうかを判定する
~：True/False を反転する
x[mask]：mask が True の場所だけ取り出す

コードの意味

soi_hp は high-pass 的な成分で、ゆっくり変わる長周期成分を差し引いた残差
ただし soi_ma121 の端は NaN なので、引き算した結果の端も NaN になる
そこで mask_hp = ~np.isnan(soi_hp) により、「NaNではない場所だけ True」の配列を作る
time_hp = time[mask_hp] は、時刻も同じ場所だけ残す処理
soi_hp_valid = soi_hp[mask_hp] は、実際に計算できた部分だけを取り出している

考えてみよう

~np.isnan(...) は何をしているか
なぜ移動平均を引いた直後のデータをそのまま使えないのか
time にも同じ mask をかける理由は何か

plt.figure(figsize=(12,4))
plt.plot(time_hp, soi_hp_valid, linewidth=1.2)
plt.axhline(0, color="k", linewidth=0.8)
plt.xlabel("Year")
plt.ylabel("SOI")
plt.title("SOI after removing >10-year component")
plt.grid(True, alpha=0.3)
plt.tight_layout()
plt.show()

STEP 4

6. 13か月移動平均

さらに 13か月移動平均を行い、1年以内の短周期変動を弱める。こうして、年々変動を中心に見やすい時系列を作る。

soi_filt = moving_average(soi_hp_valid, ____)
mask_filt = ~np.isnan(soi_filt)

time_filt = time_hp[mask_filt]
soi_filt_valid = soi_filt[mask_filt]

コードの意味

ここでも同じ moving_average 関数を再利用している
121か月移動平均では長周期を取り除き、13か月移動平均では月ごとの細かいギザギザを弱める
この結果、2〜5年程度の変動が相対的に目立ちやすくなる
13か月移動平均でも端は NaN になるため、もう一度 mask を作って有効データだけ残す

考えてみよう

13 という数字は何を落としたい処理か
121か月移動平均のあとに、なぜさらに13か月移動平均を行うのか
soi_hp_valid と soi_filt_valid はどう違うか

plt.figure(figsize=(12,4))
plt.plot(time_filt, soi_filt_valid, linewidth=2)
plt.axhline(0, color="k", linewidth=0.8)
plt.xlabel("Year")
plt.ylabel("SOI")
plt.title("Filtered SOI")
plt.grid(True, alpha=0.3)
plt.tight_layout()
plt.show()

121か月移動平均で長周期を落とし、13か月移動平均で短周期を弱めるので、結果として数年スケールの変動が目立ちやすくなる。

STEP 5

7. 自己相関係数

自己相関とは、同じ時系列を前後にずらして比較したときの相関である。ここでは教科書に合わせて、内積型の定義で自己相関係数を求める。

def autocorrelation_fortran_style(x, max_lag):
    x = np.asarray(x)
    lags = np.arange(-max_lag, max_lag + 1)
    acf = np.full(len(lags), np.nan)
    n = len(x)

    for i, lag in enumerate(lags):
        if lag < 0:
            x1 = x[-lag:]
            x2 = x[:n + lag]
        elif lag > 0:
            x1 = x[:n - lag]
            x2 = x[lag:]
        else:
            x1 = x
            x2 = x

        num = np.sum(__________)
        den = np.sqrt(np.sum(x1**2) * np.sum(x2**2))
        if den > 0:
            acf[i] = num / den

    return lags, acf

lags, acf = autocorrelation_fortran_style(soi_filt_valid, ____)

このブロックで新しく出るもの

np.asarray(x)：x を NumPy 配列として扱う
np.arange(-max_lag, max_lag + 1)：lag の一覧を作る
np.full(len(lags), np.nan)：結果を入れる NaN 配列を作る
for i, lag in enumerate(lags)：番号 i と lag の値を同時に回す
x[a:b]：配列の一部を切り出す
np.sum(...)：合計を計算する
np.sqrt(...)：平方根をとる

コードの意味

lags は何か月ずらして比較するかの一覧で、たとえば -120 から +120 まで作る
acf は各 lag に対応する自己相関係数を保存する配列
lag < 0 のときは、一方を前へ、もう一方を後ろへずらした形で比較する
lag > 0 のときは、その逆向きにずらして比較する
lag = 0 なら、自分自身との比較になる
x1 と x2 は、比較する2本の同じ長さの配列
num は内積、den は大きさの正規化で、これにより値が -1 から 1 の範囲に入りやすくなる

特に大事なところ

x[-lag:] は、末尾側を切り出す書き方
x[:n + lag] は、先頭から n+lag 個分を切り出す書き方
この2つを組み合わせると、「ずらしたあとで重なっている部分だけ」を比較できる
つまり、端のそろわない部分は自動的に除外している

考えてみよう

lag = 0 のとき、なぜ自己相関は 1 になるのか
なぜ左右ほぼ対称な形になるのか
num = np.sum(x1 * x2) は何を意味しているか

plt.figure(figsize=(10,5))
plt.plot(lags, acf, color="k", linewidth=1.5)
plt.axhline(0, color="k", linestyle="--", linewidth=0.8)
plt.axvline(0, color="k", linewidth=0.8)
plt.xlim(-120, 120)
plt.ylim(-1, 1)
plt.xlabel("Month")
plt.ylabel("Auto-Correlation")
plt.title("Autocorrelation of filtered SOI")
plt.tight_layout()
plt.show()

この描画で新しく出るもの

plt.axvline(0)：x=0 の縦線を引く
plt.xlim(-120, 120)：横軸の範囲を決める
plt.ylim(-1, 1)：縦軸の範囲を決める
linestyle="--"：破線にする

8. 結果の見方

自己相関図では、数年おきに山や谷が現れるかを見る。SOI は完全な正弦波ではないので、3年・4年・5年で必ず明確な正のピークが出るとは限らない。

図	見るポイント
121か月移動平均	10年以上のゆっくりした変動を表しているか
差し引き後の SOI	長周期成分が弱まっているか
13か月移動平均後	数年スケールの変動が見やすくなっているか
自己相関	約 48 か月前後に山があるか、左右対称か

ここでは「きれいな単一周期」を探すというより、SOI が 2〜5 年程度の時間スケールで変動しているかを確認することが大事である。

8. 解答の表示

まずは上の穴埋めコードを自分で考えてみること。確認したい場合だけ、パスワードを入力して解答と完全スクリプトを表示する。

穴埋め解答

soi_ma121 = moving_average(soi, 121)
soi_filt = moving_average(soi_hp_valid, 13)
num = np.sum(x1 * x2)
lags, acf = autocorrelation_fortran_style(soi_filt_valid, 120)

完全スクリプト

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

# 1. Read data
df = pd.read_csv("SOI_timeseries_updated.txt", sep=r"\s+")
df = df.dropna(subset=["SOI"]).reset_index(drop=True)
df = df[(df["YEAR"] >= 1951) & (df["YEAR"] <= 1997)].reset_index(drop=True)

time = df["TIME"].to_numpy()
soi = df["SOI"].to_numpy()

# 2. Moving average
def moving_average(x, window):
    return pd.Series(x).rolling(window=window, center=True).mean().to_numpy()

soi_ma121 = moving_average(soi, 121)
soi_hp = soi - soi_ma121
mask_hp = ~np.isnan(soi_hp)

time_hp = time[mask_hp]
soi_hp_valid = soi_hp[mask_hp]

soi_filt = moving_average(soi_hp_valid, 13)
mask_filt = ~np.isnan(soi_filt)

time_filt = time_hp[mask_filt]
soi_filt_valid = soi_filt[mask_filt]

plt.figure(figsize=(12,4))
plt.plot(time, soi, label="Original SOI", alpha=0.5)
plt.plot(time, soi_ma121, label="121-month moving average", linewidth=2)
plt.axhline(0, color="k", linewidth=0.8)
plt.xlabel("Year")
plt.ylabel("SOI")
plt.title("SOI and 121-month moving average")
plt.legend()
plt.grid(True, alpha=0.3)
plt.tight_layout()
plt.show()

plt.figure(figsize=(12,4))
plt.plot(time_hp, soi_hp_valid, linewidth=1.2)
plt.axhline(0, color="k", linewidth=0.8)
plt.xlabel("Year")
plt.ylabel("SOI")
plt.title("SOI after removing >10-year component")
plt.grid(True, alpha=0.3)
plt.tight_layout()
plt.show()

plt.figure(figsize=(12,4))
plt.plot(time_filt, soi_filt_valid, linewidth=2)
plt.axhline(0, color="k", linewidth=0.8)
plt.xlabel("Year")
plt.ylabel("SOI")
plt.title("Filtered SOI")
plt.grid(True, alpha=0.3)
plt.tight_layout()
plt.show()

def autocorrelation_fortran_style(x, max_lag):
    x = np.asarray(x)
    lags = np.arange(-max_lag, max_lag + 1)
    acf = np.full(len(lags), np.nan)
    n = len(x)

    for i, lag in enumerate(lags):
        if lag < 0:
            x1 = x[-lag:]
            x2 = x[:n + lag]
        elif lag > 0:
            x1 = x[:n - lag]
            x2 = x[lag:]
        else:
            x1 = x
            x2 = x

        num = np.sum(x1 * x2)
        den = np.sqrt(np.sum(x1**2) * np.sum(x2**2))
        if den > 0:
            acf[i] = num / den

    return lags, acf

lags, acf = autocorrelation_fortran_style(soi_filt_valid, 120)

plt.figure(figsize=(10,5))
plt.plot(lags, acf, color="k", linewidth=1.5)
plt.axhline(0, color="k", linestyle="--", linewidth=0.8)
plt.axvline(0, color="k", linewidth=0.8)
plt.xlim(-120, 120)
plt.ylim(-1, 1)
plt.xlabel("Month")
plt.ylabel("Auto-Correlation")
plt.title("Autocorrelation of filtered SOI")
plt.tight_layout()
plt.show()