見出し画像

「入門はじめての統計解析」をPythonで写経 Vol.20 ~ 7章「はじめての時系列分析」①トレンド、周期変動、不規則変動

7章「はじめての時系列分析」

書籍の著者 石村貞夫 先生


この記事は、書籍「入門はじめての統計解析」7章「はじめての時系列分析」の Python写経活動 を取り扱います。

書籍の図・表・計算を淡々とPython化する写経シリーズです。
7章は時系列データの分析に取り組みます。

この記事は テキストが3つの基本時系列と呼ぶ トレンド周期変動不規則変動 を実践します。
ChatGPTの活用も継続してまいります!

では書籍を開いて統計解析の旅に出発です🚀

タイムマシンのイラスト:「いらすとや」さんより

はじめに


このブログシリーズは、書籍「入門はじめての統計解析」(東京図書、「テキスト」と呼びます)の Python 写経を通じて得た「統計解析の楽しさ」をご紹介します。

書籍の紹介と引用表記はリンク先の記事に掲載しています。

7章 はじめての時系列分析


この記事は7章の以下のSectionを取り扱います。

7.1 3つの基本時系列

記事に用いるデータは、テキストに掲載されたデータそのものを引用しています。
データ件数の少ないものはコード上でデータを登録し、データ件数の多いものはCSVファイル化してデータを読み込みしています。

7章で用いるライブラリをインポートします。

### インポート

# 数値計算
import statistics                           # Python標準ライブラリ
import numpy as np
import pandas as pd

# 統計
import scipy.stats as stats
import pymannkendall as mk                  # mann-kendall検定(トレンドの検定)

# 時系列分析
import statsmodels.api as sm
from scipy.signal import periodogram        # ピリオドグラム
import statsmodels.tsa.api as tsa           # 指数平滑法

# 描画
import matplotlib.pyplot as plt
plt.rcParams['font.family'] = 'Meiryo'

イントロダクション

ChatGPTが「時間の旅」へと、誘います。


✅ 未来は「今」の延長線じゃない。だから、時系列分析。

ビジネスは、時間で動いています。
昨日と今日、去年と今年――そこには見えない“流れ”があります。
時系列分析は、その流れを読み解く技術。
変化の兆しを見つけ、トレンドをつかみ、未来の形をそっと描き出します。

でも、気をつけて。
未来は、過去のコピーではありません。
そこには、予測できない揺らぎも、まだ見ぬチャンスもある。
だから、数字を見るだけじゃ足りない。
私たちは、不確実性の中にも描くのです。

時間を味方にする――その第一歩が、ここにあります。

ChatGPTで生成した時系列分析のイメージ

➡️ Section 7.1 3つの基本時系列

「時系列」にまつわる用語をテキストから引用いたします。

時系列は「時間と共に変化するデータの列」
時系列データは「時間を $${t}$$ とすると $${\{ \cdots,\ x(t-3),\ x(t-2),\ x(t-1),\ x(t)\}}$$ 」と表現できる
時系列分析は「時系列データから未来を予測する統計手法」

テキストより引用

ということで、仮想の時系列データ $${y}$$ を生成して、テキストの図 7.1.2「人工的株価変動」のチャートを描いてみましょう。

### 人工的株価変動 p.267 図7.1.2

## 設定と準備
# グラフタイトル
titles = ['$y_1=0.8x$', '$y_2=10\ \sin(x/4)$', '$y_3=10(random-0.5)$',
          '$y=y_1+y_2+y_3$']
# 乱数生成器
rng = np.random.default_rng(seed=0)

## データの作成
x = np.linspace(-50, 50, 201)
y1 = 0.8 * x
y2 = 10 * np.sin(x / 4)
y3 = 10 * (rng.random(size=len(x)) - 0.5)
y = y1 + y2 + y3
y_list = [y1, y2, y3, y]

## 描画
fig, axes = plt.subplots(2, 2, figsize=(10, 5), tight_layout=True)
for y_plot, title, ax in zip(y_list, titles, axes.flat):
    ax.plot(range(len(x)), y_plot)
    ax.set(title=title, xlabel='時間 $t$', yticks=[])
plt.show()

【実行結果】
$${y_1}$$「直線」、$${y_2}$$「周期のある曲線」、$${y_3}$$「ランダムな線」を足し合わせたデータが右下の $${y}$$ です。
テキストは $${y}$$ を「平均株価のようなグラフ」と名付けています。

続けてテキストは上のチャートを参照しつつ、時系列データが次の「3つの基本時系列」から成ることを説明します。

  • トレンド($${y_1}$$)

  • 周期変動($${y_2}$$)

  • 不規則変動($${y_3}$$)

テキストの図 7.1.3「時系列の3つの基本パターン」で確認しましょう。
上のチャートの $${y_1, y_2, y_3}$$ を再度描画するものです。

### 時系列の3つの基本パターン p.267 図7.1.3

## 設定と準備
titles = ['(a) トレンド', '(b) 周期変動', '(c) 不規則変動']
y_list = [y1, y2, y3]

## 描画
fig, axes = plt.subplots(1, 3, figsize=(10, 2.5), tight_layout=True)
for y_plot, title, ax in zip(y_list, titles, axes.flat):
    ax.plot(range(len(x)), y_plot)
    ax.set(title=title, xlabel='時間 $t$', yticks=[])
plt.show()

【実行結果】

テキストに沿って「トレンド」「周期変動」「不規則変動」を深堀りしていきます!
時系列データ $${y}$$ はこの先も活用しますので、頭の片隅に置いておいてください。

トレンド

テキストの「トレンドで大切なポイントは?」を引用いたします。

1.時系列データにトレンドがあるといってよいか?
▶️ ケンドールによるトレンドの検定

2.時系列にトレンドがある場合、将来の値を予測しよう
▶️ 曲線の当てはめ
▶  指数平滑化

▶️を付与した2つのテーマを実践します。
指数平滑化は次回の記事で実践します。

▶️ ケンドールによるトレンドの検定
次の仮説に関するノンパラメトリック検定です。
・帰無仮説 $${H_0}$$「トレンドは存在しない」
・対立仮説 $${H_1}$$「トレンドが存在する」

Python のライブラリ pyMannKendall で「マン・ケンドール検定」(ケンドールによるトレンドの検定)を実行できます。
公式サイトはこちら!

仮想の時系列データ $${y}$$ にトレンドが存在するかどうか、マン・ケンドール検定で確かめましょう。
$${y}$$ の成分にトレンド $${y_1}$$ を含んでいるので「トレンドが存在する」はずです!

「季節性(周期性)・自己相関を考慮しないオリジナル検定」と「季節性を考慮する季節性検定」を実行します。

◆ オリジナル検定
Mann-Kendall の original_test 関数に時系列データの変数名を与えるだけです!

mk.original_test(y)

### ケンドールによるトレンドの検定 Mann-Kendall検定 p.269
# https://github.com/mmhs013/pyMannKendall

## yのトレンドの検定(オリジナル検定)
trend, h, p, z, tau, s, var_s, slope, intercept  = mk.original_test(y)
print(f'Test Results: trend={trend}, is trend?={h}, p_value={p:.4f}')

## トレンドの描画
# トレンド推定値の算出
trend_line = np.arange(len(y)) * slope + intercept
# 描画領域の設定
plt.figure(figsize=(10, 4))
# 観測値y(青)の描画
plt.plot(range(len(y)), y, label='観測値')
# トレンド真値y1(グレー)の描画
plt.plot(range(len(y)), y1, color='grey', label='トレンド真値')
# トレンド推定値(赤点線)の描画
plt.plot(range(len(y)), trend_line, color='tab:red', ls='--',
         label='トレンド推定値')
# 修飾
plt.xlabel('時間 t')
plt.yticks([])
plt.legend();

【実行結果】
デフォルトの有意水準 $${5\%}$$ で有意であり「上昇トレンドが存在する」と言えます。

グレーの「トレンドの真値 $${=y_1}$$」と赤い点線のオリジナル検定によるトレンド推定値はかなり近く、若干のずれが見られる程度です。

◆ 季節性検定
Mann-Kendall の seasonal_test 関数に時系列データの変数名と周期を与えます。

mk.seasonal_test(y, period=period)

## yのトレンドの検定(季節性検定)
period = 51
trend, h, p, _, _, _, _, slope, intercept = mk.seasonal_test(y, period=period)
print(f'Test Results: trend={trend}, is trend?={h}, p_value={p:.4f}')

## トレンドの描画
# トレンド推定値の算出
trend_line = np.arange(len(y)) / period * slope + intercept
# 描画領域の設定
plt.figure(figsize=(10, 4))
# 観測値y(青)の描画
plt.plot(range(len(y)), y, label='観測値')
# トレンド真値y1(グレー)の描画
plt.plot(range(len(y)), y1, color='grey', label='トレンド真値')
# トレンド推定値(赤点線)の描画
plt.plot(range(len(y)), trend_line, color='tab:red', ls='--',
         label='トレンド推定値(季節性考慮)')
# 修飾
plt.xlabel('時間 t')
plt.yticks([])
plt.legend();

【実行結果】
こちらもデフォルトの有意水準 $${5\%}$$ で有意であり「上昇トレンドが存在する」と言えます。

トレンド推定値は真値とほぼ一致しています。
周期性を捉えているように見えます。

▶️ 曲線の当てはめ
テキストの図 7.1.6「明日の株価は⋯?」に相当する「明日の予測図」を描きましょう。
変数 $${y}$$ が時間 $${t}$$ に対して次式(曲線)で表せるとします。

$$
y(t) = \cfrac{t+t^3}{20} + 誤差
$$

### 明日の株価は…? p.269 図7.1.6

## 設定と準備
# 曲線関数の定義
f = lambda t: (t + t**3) / 20
# 描画データの作成
time = np.arange(1, 7)
y_obs = f(time) + np.array([0.2, -0.2, 0.4, -0.4, 0, 0])
# 描画用設定
blue, red, white, black = 'tab:blue', 'tab:red', 'white', 'black'

## 描画
# 描画領域の設定
fig, ax = plt.subplots(figsize=(7, 3))
# 今日までの株価に基づく推定値(青実線)の描画
ax.plot(t_val:=np.linspace(1, 5), f(t_val), color=blue)
# 明日の株価予測線(赤点線)の描画
ax.plot(t_val:=np.linspace(5, 6), f(t_val), ls='--', color=red)
# 昨日・今日・明日の垂直点線の描画
ax.vlines([4, 5, 6], -0.3, y_obs[3:], color=black, ls='--', lw=0.5)
# 今日までの株価のデータ点(青丸)の描画
ax.plot(time[:-1], y_obs[:-1], 'o', ms=8, mec=blue, mfc=white)
# 明日の株価予測点(赤丸)の描画
ax.plot(time[-1], y_obs[-1], 'o', ms=8, color=red)
# 修飾
ax.set_xticks(ticks=time, labels=['', '', '', '昨日', '今日', '明日'], fontsize=12)
ax.set_ylabel('株価', fontsize=12)
ax.set(ylim=(-0.3, y_obs[-1]*1.05), yticks=[])
plt.show()

【実行結果】
今日までの実測値(白抜きの点)を当てはめた「曲線」から、明日の株価(赤い点)を予測しています。

周期変動

テキストは次の2つの用語を使い分けしています。
・周期変動:周期的に繰り返されている時系列
・季節変動:四季や12ヶ月のように季節的に繰り返す周期変動

テキストの「周期変動で大切なポイントは?」を引用いたします。

1.どのような周期で繰り返しているのか?
▶️ スペクトル分析

2.時系列に周期変動がある場合、将来の値を予測しよう
▶  指数平滑化
▶  自己回帰モデル

3.季節変動のときの分析手法(例)
▶  12 ヶ月移動平均
▶️ 季節性の分解

▶️を付与した2つのテーマを実践します。
指数平滑化、移動平均、自己回帰モデルは次回の記事で実践します。

▶️ スペクトル分析
ChatGPTがスペクトル分析をざっくりまとめてくれました。


🔍 スペクトル分析とは?

時系列データを 周波数成分(frequency components)に分解する手法です。

  • 周期性の強いデータは、ある周波数成分が 大きなパワー(エネルギー) を持つ

  • その周波数の逆数が、データの周期を表します

📘 ざっくり流れ(Python × scipyのピリオドグラムで作成したグラフ)

🎯 グラフの読み方のポイント

グラフの ピークが大きいほど、その周波数成分がデータに強く現れています。
周期は:

$$
周期 = \cfrac{1}{周波数}
$$

たとえば、ピークが 0.142857 ≒ 1/7 にあれば、「7日周期」の成分が強いことを意味します。

🌟 どんなときに役立つ?

  • 売上の曜日パターンや月次サイクル

  • 季節的なイベントや定期的な行動

  • 周期性があるかどうかの探索的分析


仮想の時系列データ $${y}$$ の周期をスペクトル分析の1手法「ピリオドグラム」で調べましょう。
分析に入る前に… ChatGPTが教えてくれた注意点を確認します。


トレンド(傾きやゆるやかな変化)を含んだ時系列データにピリオドグラムをそのまま適用するときには、注意が必要です。以下に詳しく解説します。

🎯 なぜ注意が必要なのか?

✅ 理由:トレンド成分は「低周波成分」として現れるため

  • ピリオドグラムは周波数ごとのエネルギー(パワー)を可視化する方法

  • 緩やかな上昇/下降トレンドは、非常に低い周波数(≒周期が長い)として現れます

  • その結果、「0に近い周波数」で大きなピークが出てしまうことがあります

このような「トレンド由来の低周波ピーク」は、周期性ではないにもかかわらず強く出てしまうため、本来検出したい周期のピークを見えにくくしてしまうのです。

✅ 解決策は…
トレンド除去 or 前処理

🎼 イメージで例えるなら

  • トレンドは重低音のドローン

  • あなたが聴きたいのはその上で踊るリズムやメロディ(=周期成分)

  • ドローン音を一度下げてから、リズムを聴き取るような作業が「トレンド除去」

✅ 結論


■ トレンド成分を除去した時系列データでスペクトル分析

さきほどマン・ケンドール検定でトレンド成分を推定しました。
仮想の時系列データ $${y}$$ から推定したトレンド成分を除去して、ピリオドグラムを実施しましょう。

まずはトレンド除去です。

### スペクトル分析の準備:トレンド除去

# 仮想の時系列データからトレンド成分を除去 ※trend_line はマン・ケンドール検定で推定
y_detrend = y - trend_line

# 結果の可視化
plt.figure(figsize=(10, 4))
plt.plot(range(len(y)), y_detrend)
plt.title('仮想の時系列データ:トレンド除去後')
plt.xlabel('時間 $t$', fontsize=12);

【実行結果】
上昇トレンドがなくなり、一定の周期で繰り返す感じになりました。

いよいよピリオドグラムで周期を推定します!
scipy の periodogram を利用します。
引数にトレンド除去後の仮想データを与えます。

### スペクトル分析

## ピリオドグラムを用いてパワースペクトル密度を推定
# ピリオドグラム法の実行
freqs, power = periodogram(y_detrend)
# 最大のパワースペクトルとその周波数の取得
peak_power = np.max(power)
peak_freq = freqs[np.argmax(power)]
# 推定周期の算出
period_est = 1 / peak_freq

## プロット
# 描画領域の設定
plt.figure(figsize=(10, 4))
# スペクトル分析の結果の描画
plt.plot(freqs, power)
# 最大のパワースペクトルの点の描画
plt.plot(peak_freq, peak_power, 'o', color='tab:red')
# 修飾
plt.xlabel('周波数 [1/t]')
plt.ylabel('パワー')
plt.title(f'スペクトル分析:パワー最大の周波数 {peak_freq:.5f}, 周期 {period_est}')
plt.grid(True)
plt.show()

【実行結果】
周波数 $${0.0199}$$ がパワー最大となり、この周波数の逆数=周期は $${50.25}$$ です。
仮想データはだいたい $${50}$$ の周期で繰り返しているようですね!

▶️ 季節性の分解
季節成分のコレログラムを見てみましょう。
コレログラムは自己相関・偏自己相関と呼ばれる「時系列データとその時系列データの◯時点前との相関」を可視化するものです。
時点のずれは「ラグ」と呼ばれ、1時点前は「ラグ1」、2時点前は「ラグ2」です。

仮想の時系列データ $${y}$$ の季節成分 $${y_2}$$ をコレログラムで可視化します。
statsmodels の plot_acf() と plot_pacf() を利用します。

### 季節性の分解 p.271

# 周期変動の自己相関・偏自己相関の可視化
fig, (ax1, ax2) = plt.subplots(2, 1, figsize=(8, 6), tight_layout=True)
sm.graphics.tsa.plot_acf(y2, lags=100, auto_ylims=True, ax=ax1)
sm.graphics.tsa.plot_pacf(y2, lags=100, auto_ylims=True, ax=ax2);

【実行結果】
上が自己相関、下が偏自己相関です。
横軸は「ラグ」、青く塗られた領域は 95% 信頼区間です。
信頼区間を超えている点(白い領域)が統計的に有意となります。

自己相関の方に周期性が見られますね!
きれいな波が描かれています。
ラグ 1~9 とラグ 18~29 あたりまでが有意です。

偏自己相関は比較対象のラグだけの影響を抽出しています。
他のラグの影響を取り除いています。
ラグ 1 ~ 7 あたりが有意のようです。

■ 季節成分の分離

時系列データの趨勢・傾向など分析する際に季節成分を取り除くことがあります。
仮想の時系列データ $${y}$$ から季節成分を分離してみましょう。
同時にトレンド成分、不規則変動成分も分離します。
statsmodels の seasonal_decompose() を利用します。

# 季節性の分解
result = tsa.seasonal_decompose(y, period=int(period_est), two_sided=True)
fig = result.plot()
fig.set_size_inches(8, 8);

【実行結果】
上から「仮想データそのもの(原系列)」「トレンド成分」「季節成分」「不規則変動成分」です。
トレンド:長期的な傾向が直線的な上昇傾向であることが分かります。

記事の最初の $${y_1, y_2, y_3}$$ のプロット(真のトレンド・季節成分・不規則変動)と比べると…

  • 直線的なトレンド成分をうまく分離できているように見えます。

  • 季節成分はギザギザしています。

不規則変動

テキストの「不規則変動で大切なポイントは?」を引用いたします。

1.時系列データは本当に不規則なのか?
▶️ 連の総数によるランダムの検定

2.時系列データがホワイトノイズか?
▶️ ボックス・リュングの検定

▶️を付与した2つのテーマを実践します。

▶️ 連の検定
連の検定は「データがランダムなのか」に関する統計的検定です。
仮説は:
帰無仮説 $${H_0}$$「2つに分割したデータの並び方はランダムである」
対立仮説 $${H_1}$$「2つに分割したデータの並び方はランダムではない」
ChatGPTが連の検定をざっくりまとめてくれました。


🔍 連の検定ってなに?
「データに並び方の偏りがないか」を調べる検定です。

🎲 たとえば…
コインを投げて「表(H)・裏(T)」の並びを観察したとき、
「HHHHTTTT」みたいに偏った並びになっていたら、

「ほんとにランダムなの?」
「何か規則性があるんじゃない?」

…と疑いたくなりますよね。
そんなときに使えるのが「連の検定」です。

🧩「連」とは?
連(run)=同じ値が連続して続くかたまり
(例)H H H T T H T T T
→ 「HHH」「TT」「H」「TTT」 → 連の数は 4つ

📊 どんなことがわかる?

  • 連が少なすぎる → 偏り(グルーピング)があるかも?

  • 連が多すぎる → 過剰な切り替え(バラつきすぎ)かも?

⇒ 統計的に「ランダムに見えるかどうか」を判断できます!

✅ よくある使いどころ

  • コイン投げやランダム性の検証

  • 品質検査で「合格・不合格」の並びの偏り確認

  • 時系列データの“切り替わりパターン”をチェックするとき


仮想の時系列データ $${y}$$ の $${y_3}$$(不規則変動成分)がランダムなのか(本当に不規則なのか)について、連の検定で確かめましょう。

自作関数の作成に当たり、次のWebサイトの情報を利用いたしました。
ありがとうございます!

では実装に入ります。
まずは連の検定関数の定義から。

### 連の検定関数 p.273

# z検定関数
def z_test(z_val, alpha=0.05, alternative='two-sided'):

    # 標準正規分布の設定 ※scipy.stats利用
    std_normal_dist = stats.norm(loc=0, scale=1)
    
    # 棄却限界値c_valとp値p_valの算出
    p_val = std_normal_dist.sf(x=abs(z_val))
    if alternative == 'two-sided':
        c_val = std_normal_dist.ppf(q=1 - alpha/2)
        c_val = -c_val, c_val
        p_val = p_val * 2
    elif alternative=='less':
        c_val = std_normal_dist.ppf(q=alpha)
    elif alternative=='greater':
        c_val = std_normal_dist.ppf(q=1 - alpha)
    
    # 戻り値:
    return {'alternative': alternative, 'c_value': c_val, 'z_value': z_val,
            'alpha': alpha, 'p_value': p_val}


# 連の検定関数
def runs_test(list_x, cutoff='mean', alpha=0.05, alternative='two-sided'):

    ### 基準値を用いてデータを二値分類 ※numpy利用
    # 標本サイズの算出
    N = len(list_x)
    # 基準値の設定 ※デフォルト:mean(平均), 選択肢:median(中央値)
    cutoff_value = (
        statistics.median(list_x) if cutoff=='median' else sum(list_x) / N)
    # データxを0,1に分類 ※x>cutoff_value: 0 as A, x<=cutoff_value: 1 as B
    bins_x = [int(x <= cutoff_value) for x in list_x]
    
    ### 連の処理
    # 基準値を超える0の個数na, 基準値以下1の個数nbの算出
    na, nb = N - sum(bins_x), sum(bins_x)
    # 連の数Kの算出
    K = sum([x != y for x, y in zip(list(bins_x)[1:], list(bins_x)[:-1])]) + 1
    
    ### 検定処理
    # 連の数Kが従う正規分布の平均の算出
    mu_K = 2 * na * nb / N + 1
    # 連の数Kが従う正規分布の分散の算出
    sigma2_K = (2 * na * nb * (2 * na * nb - N)) / (N**2 * (N - 1))
    # 検定統計量zの算出
    z = (K - mu_K) / sigma2_K**(1/2)
    # z検定処理
    result = z_test(z, alpha, alternative)
    
    ## 戻り値
    return result | {'N': N, 'N_A': na, 'N_B': nb,
            'cutoff method': cutoff, 'cutoff_value': cutoff_value, 'runs_K': K}

では仮想の時系列データ $${y}$$ の $${y_3}$$(不規則成分)で連の検定を実行します。
データの中央値で2分割して、2つのデータ間の並び方がランダムかどうかを調べます。
$${y_3}$$ は一様分布乱数なので「ランダム・不規則」であって欲しいです!

### 連の検定 p.273
# 帰無仮説:2つに分割したデータの並び方はランダムである

# 関数利用
result = runs_test(y3, cutoff='median')  # cutoff: データを2分割する際の境界値
result['z_value'], result['p_value']

【実行結果】
検定統計量と $${p}$$ 値です。

$${p}$$ 値 $${0.1325 > 0.05}$$ であり、有意水準 $${5\%}$$ で帰無仮説を棄却できません。
「2つに分割したデータの並び方はランダムではない」とは言えない、ので、$${y_3}$$(不規則成分)はランダムだろう、という感じです。

Python のライブラリも試してみましょう。
statsmodels の連の検定です。
なお、執筆現在、sandbox 扱いのため、将来に実装方法が変わる可能性があります。

🖲️statsmodels の runstest_1samp()

# statsmodels利用 (z値, p値)
sm.stats.runstest_1samp(y3, cutoff='median')

【実行結果】

▶️ ホワイトノイズ
テキストのホワイトノイズの定義をお借りします。

確率変数の列$${\{\cdots, X(t-2), X(t-1), X(t), X(t+1)\}}$$が次の性質を満たすとき、この列をホワイトノイズという。

$$
\begin{align*}
平均\ &E(X(t)) = 0 \\
分散\ &\text{Var}(X(t)) = \sigma^2 \\
共分散\ &\text{Cov}(X((t)), X(t-1)) = 0 \quad (s=\cdots, -2, -1, 1, 2, \cdots) \\
\end{align*}
$$

テキストの数式を引用

ホワイトノイズはざっくり、どの時点でも「平均0」「分散一定」であり、異なる時点間の「共分散=0」となります。

ホワイトノイズデータを作ってみましょう。
標準正規分布 $${\text{Normal}\ (0, 1^2)}$$ に従う乱数です。
平均0、分散1です。

### ホワイトノイズの例 (標準正規分布乱数)

# 設定と準備
T = 100                               # 時間の数
rng = np.random.default_rng(seed=14)  # 乱数生成器の設定:14, 19, 27, 28, 30

# ホワイトノイズの生成:正規分布 N(0, σ²=1)に従う乱数なので正規ホワイトノイズ
white_noise = rng.standard_normal(size=T)

# 可視化
plt.figure(figsize=(10, 4))
plt.plot(white_noise)
plt.title(f'正規ホワイトノイズ $\sim Normal\ (0, 1^2)$')
plt.xlabel('時間 $t$');

【実行結果】

■ ホワイトノイズと自己相関

ホワイトノイズは「自己相関がない」性質を持ちます。
コレログラムで確認しましょう。

# ホワイトノイズの自己相関・偏自己相関の可視化
fig, (ax1, ax2) = plt.subplots(2, 1, figsize=(8, 6), tight_layout=True)
sm.graphics.tsa.plot_acf(white_noise, lags=20, auto_ylims=True, ax=ax1)
sm.graphics.tsa.plot_pacf(white_noise, lags=20, auto_ylims=True, ax=ax2);

【実行結果】
ラグ1以降、有意な自己相関、偏自己相関はありません。
ホワイトノイズが自己相関を持たないことが分かりました。

▶️ リュング・ボックス検定(ボックス・リュングの検定)
リュング・ボックス検定は「時系列データがホワイトノイズなのか」に関する統計的検定です。
仮説は:
帰無仮説 $${H_0}$$「データはホワイトノイズである」
対立仮説 $${H_1}$$「データはホワイトノイズではない」

さきほど作ったホワイトノイズデータでリュング・ボックス検定を実行しましょう。
statsmodels の acorr_ljungbox() を利用します。
ラグ10 まで確認します。

# ホワイトノイズのリュング・ボックス検定
# statsmodels利用 ※index:ラグ,lb_stat列:検定統計量, lb_pvalue:p値
# 帰無仮説:データはホワイトノイズである

result = sm.stats.acorr_ljungbox(white_noise)
result.index.name = 'ラグ'
result.columns = ['検定統計量', 'p値']
result.round(3)

【実行結果】
ラグ10までのすべてで $${p}$$ 値 $${> 0.05}$$ となっており、有意水準 $${5\%}$$ で有意とは言えず、帰無仮説を棄却できません。
帰無仮説「データはホワイトノイズである」を受容します。

続いて仮想の時系列データ $${y}$$ の不規則成分 $${y_3}$$ がホワイトノイズかどうかをリュング・ボックス検定で確かめましょう。

# 仮想データの不規則成分 y3 のリュング・ボックス検定
# statsmodels利用 ※index:ラグ,lb_stat列:検定統計量, lb_pvalue:p値
# 帰無仮説:データはホワイトノイズである

result = sm.stats.acorr_ljungbox(y3)
result.index.name = 'ラグ'
result.columns = ['検定統計量', 'p値']
result.round(3)

【実行結果】
ラグ10までのすべてで $${p}$$ 値 $${> 0.05}$$ となっており、有意水準 $${5\%}$$ で有意とは言えず、帰無仮説を棄却できません。
帰無仮説「データはホワイトノイズである」を受容します。
ですが、ホワイトノイズデータと比べると検定統計量の値が大きく、$${p}$$ 値は小さくなっています。


記事の最後をChatGPTに締めくくってもらいましょう!
「線の物語」はこれからも続いていくことでしょう。

📘 ChatGPTのひとこと:

データは静かに、でも確かに歌っています。
トレンドが旋律を描き、周期がリズムを刻み、不規則な変動さえも自由な即興演奏のように。
時系列分析とは、その楽譜を読み解くこと。
このメロディが続く先に、しばらく耳を澄ませてみましょう。
未来の音色がするかもしれません。

今回の写経は以上です。


シリーズの記事

次の記事

前の記事

目次

ブログの紹介


note で7つのシリーズ記事を書いています。
ぜひ覗いていってくださいね!

1.のんびり統計

統計検定2級の問題集を手がかりにして、確率・統計をざっくり掘り下げるブログです。
雑談感覚で大丈夫です。ぜひ覗いていってくださいね。
統計検定2級公式問題集CBT対応版に対応しています。
Python、EXCELのサンプルコードの配布もあります。

2.実験!たのしいベイズモデリング1&2をPyMC Ver.5で

書籍「たのしいベイズモデリング」・「たのしいベイズモデリング2」の心理学研究に用いられたベイズモデルを PyMC Ver.5で描いて分析します。
この書籍をはじめ、多くのベイズモデルは R言語+Stanで書かれています。
PyMCの可能性を探り出し、手軽にベイズモデリングを実践できるように努めます。
身近なテーマ、イメージしやすいテーマですので、ぜひぜひPyMCで動かして、一緒に楽しみましょう!

3.実験!岩波データサイエンス1のベイズモデリングをPyMC Ver.5で

書籍「実験!岩波データサイエンスvol.1」の4人のベイジアンによるベイズモデルを PyMC Ver.5で描いて分析します。
この書籍はベイズプログラミングのイロハをざっくりと学ぶことができる良書です。
楽しくPyMCモデルを動かして、ベイズと仲良しになれた気がします。
みなさんもぜひぜひPyMCで動かして、一緒に遊んで学びましょう!

4.楽しい写経 ベイズ・Python等

ベイズ、Python、その他の「書籍の写経活動」の成果をブログにします。
主にPythonへの翻訳に取り組んでいます。
写経に取り組むお仲間さんのサンプルコードになれば幸いです🍀

5.RとStanではじめる心理学のための時系列分析入門 を PythonとPyMC Ver.5 で

書籍「RとStanではじめる心理学のための時系列分析入門」の時系列分析をPythonとPyMC Ver.5 で実践します。
この書籍には時系列分析のテーマが盛りだくさん!
時系列分析の懐の深さを実感いたしました。
大好きなPythonで楽しく時系列分析を学びます。

6.データサイエンスっぽいことを綴る

統計、データ分析、AI、機械学習、Pythonのコラムを不定期に綴っています。
統計・データサイエンス書籍にまつわる記事が多いです。
「統計」「Python」「数学とPython」「R」のシリーズが生まれています。

7.Python機械学習プログラミング実践記

書籍「Python機械学習プログラミング PyTorch & scikit-learn編」を学んだときのさまざまな思いを記事にしました。
この書籍は、scikit-learnとPyTorchの教科書です。
よかったらぜひ、お試しくださいませ。

最後までお読みいただきまして、ありがとうございました。

いいなと思ったら応援しよう!

ネイピア DS 応援ありがとうございます。これからもがんばって記事を作成します!

この記事が参加している募集