「入門はじめての統計解析」をPythonで写経 Vol.21 ～７章「はじめての時系列分析」②３項移動平均、指数平滑化、自己回帰モデル

ネイピア DS

2025年7月12日 11:46

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

書籍の著者　石村貞夫先生

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

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

この記事は 移動平均、指数平滑化、自己回帰モデル を実践します。
ChatGPTの活用も継続してまいります！

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

はじめに

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

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

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

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

7.2 ３項移動平均
7.3 指数平滑化
7.4 自己回帰 AR(1) モデル

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

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

### インポート

# 数値計算
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が「時間の旅」の未知先案内人に。

🌿未来にそっと触れる分析へ

時系列データには、今日までの出来事が刻まれています。
でもそれだけではありません――
そこには「これからどうなっていくのか」という、未来へのヒントも静かに埋もれています。

今回ご紹介するのは、そんなデータの“流れ”をなめらかに整えたり、
少しだけ先の姿をやさしく予測したりするための方法です。

移動平均で波を整え、
指数平滑化で未来をそっとのぞき、
自己回帰モデルでデータが語る“自分自身とのつながり”をとらえます。

穏やかな流れをとらえることで、未来の風景が少しずつ見えてくる――
そんな感覚を、あなたと一緒に味わっていきましょう。

➡️ Section 7.2 ３項移動平均

テキストによると「移動平均とは、時系列データの変動を滑らかに変換する手法」であり、「移動平均をすることにより、時系列データのトレンドを浮かび上がらせる」ことができるとのこと。

移動平均はデータの凸凹を緩やかに（滑らかに）するものです。
３項移動平均は「ある時点とその前後の時点の３時点の平均」です。
株価の世界の移動平均は「５日移動平均線」「25日移動平均線」などがありますね！

ということで、テキストの表 7.2.1「受療者数」のデータをお借りして、図 7.2.1「受療者数」のチャートを描画しましょう。

### 受療者 p.274 表7.2.1

# データの登録
data1 = pd.DataFrame(
    {'受療者数': [576, 626, 754, 727, 823, 855, 766, 943, 926, 1005, 1092, 1105]},
    index=list(range(1994, 2006)))
data1.index.name = '年'

# 結果の表示
data1

【実行結果】
1994 年から 12 年間：年単位のデータです。

### データの可視化 p.274 図7.2.1
data1.plot(marker='D')
plt.xlim(data1.index.min(), data1.index.max())
plt.ylim(380, 1220);

【実行結果】
年々増加する傾向が見られます。

このデータの３項移動平均を計算しましょう。
pandas の rolling() メソッドを利用します。
平均する期間は window で、平均値を真ん中の時点にする center=True を設定します。
「.sum()」で３期間の合計、「.mean()」で３期間の平均を算出できます。

### 3項移動平均 p.275 表7.1.2

# rollingメソッドで窓の幅=3、中心に結果を格納する方法で合計を算出
data1['3項の合計'] = data1['受療者数'].rolling(window=3, center=True).sum()
# rollingメソッドで窓の幅=3、中心に結果を格納する方法で平均値を算出
data1['3項の平均値'] = data1['受療者数'].rolling(window=3, center=True).mean()
# 結果の表示
data1.round(2)

【実行結果】
ひとまず計算完了です。

３項移動平均の「滑らかな」時系列を可視化しましょう。

### 3項移動平均のグラフ p.275 図7.2.2

# 3項移動平均のみのグラフの描画
data1[['3項の平均値']].plot(marker='D', color='tab:orange')
plt.xlim(data1.index.min(), data1.index.max())
plt.ylim(380, 1220);

# 原系列と3項移動平均を比べるグラフの描画
data1[['受療者数', '3項の平均値']].plot(marker='D')
plt.xlim(data1.index.min(), data1.index.max())
plt.ylim(380, 1220);

【実行結果】
上がテキストの図 7.2.2「３項移動平均」に相当します。
下が元の時系列データと３項移動平均を並べた図です。

下の図を見ると元のデータの凸凹が３項移動平均では滑らかになっています。
「一貫して増加する傾向」がはっきりと分かります。

なお平均を取る過程で「両端」時点のデータが無くなることに留意しましょう。

➡️ Section 7.3 指数平滑化

テキストによると「指数平滑化は予測のための時系列分析と考えられる」とし、２つの公式を紹介しています。

指数平滑化の公式

① 指数平滑化の定義
時点 $${t}$$ における $${1}$$ 期先の予測値 $${\hat{x}(t, 1)}$$ は重み（平滑化係数）を $${\alpha}$$ とし、時点 $${t}$$ 以前の観測値 $${x(t), x(t-1), x(t-2), \cdots}$$ を用いて、以下のように定義されます。

$$
\hat{x}(t, 1) = \alpha \cdot x(t) + \alpha(1-\alpha) \cdot x(t-1) + \alpha(1-\alpha)^2 \cdot x(t-2) + \cdots
$$

テキストの数式を引用

重み $${\alpha=0.7}$$ のとき、各時点の全体の重みは

時点 $${t}$$ の重みは全体で $${\alpha=0.7}$$
時点 $${t-1}$$ の重みは全体で $${\alpha(1 - \alpha)=0.7 \times 0.3=0.21}$$
時点 $${t-2}$$ の重みは全体で $${\alpha(1 - \alpha)^2=0.7 \times 0.3^2=0.063}$$

のように最近の時点は重みが大きいので予測値に対する影響が大きく、過去になるほど重みが「指数関数的に」小さくなって予測値に対する影響が小さくなります。

◆

② １期前の予測値で再帰的に表現

$$
\hat{x}(t, 1)=\alpha \cdot x(t) + (1 - \alpha) \cdot \hat{x}(t-1, 1)
$$

テキストの数式を引用

こちらの式はシンプルですね！
ただ、時点 $${t-1}$$ の予測には時点 $${t-2}$$ の予測値が用いられ、時点 $${t-2}$$ の予測には・・・、と過去に「再帰的に」遡る必要があります。

② の１期前の予測値を用いる方法で「指数平滑化を用いた1期先予測関数」を書いてみましょう。

### 指数平滑化を用いた1期先予測関数
def exponential_smoothing_pred(list_obs, alpha):

    # 1期先予測関数
    def prediction(alpha, prev_obs, prev_pred):
        return alpha * prev_obs + (1 - alpha) * prev_pred

    # numpy配列、pandasシリーズをリスト化
    list_obs = list(list_obs)
    
    # 予測値リストの初期化: 最初の2期間に1期目の観測値を設定
    list_pred = [list_obs[0]] * 2

    # 観測値の3期~最終-1期の予測を実行して予測値リストに格納
    for obs in list_obs[1:-1]:
        pred = list_pred[-1]
        list_pred.append(prediction(alpha, obs, pred))
    
    # 誤差平方和の算出
    SSE = sum([(obs - pred)**2 for obs, pred in zip(list_obs, list_pred)])
    
    # 観測値の最終期の1期先の予測を実行
    one_period_forecast = prediction(alpha, list_obs[-1], list_pred[-1])

    # 戻り値: 1期先予測値、誤差平方和、観測値に対応する予測リスト
    return one_period_forecast, SSE, list_pred

指数平滑化例題 p.278

テキストの表 7.3.1「保釈率」のデータをお借りして、指数平滑化による１期先予測を実践しましょう。

■ 分析の準備
データを登録します。

### 指数平滑化 保釈率 p.278 表7.3.1

# データの登録
data2 = pd.DataFrame(
    {'保釈率': [47, 55, 41, 40, 36, 44, 32, 28, 27, 26, 27, 33]},
    index=range(1, 13))
data2.index.name = '時点t'
# 結果の表示
data2

【実行結果】
12 期のデータです。

時系列の推移を可視化しましょう。
pandas データフレームに対して「.plot」メソッドを使うと

### データの可視化
data2.plot(marker='D', ylim=(20, 60));

【実行結果】
全体的に減少傾向のように見えます。
時点 $${2, 6, 11, 12}$$ には増加も見られます。

■ 指数平滑化で１期先予測：重みを指定
指数平滑化で１期先予測をしましょう。
最初にテキストどおりの重み $${\alpha=0.8}$$ で計算します。

### 指数平滑化を用いた1期先予測 

# 関数利用

# 設定
alpha = 0.8
# 指数平滑化の実行
forecast, SSE, preds = exponential_smoothing_pred(data2['保釈率'], alpha)
# 結果の表示
print(f'第{data2.index.max()+1}期の予測値: {forecast:.6f}')
print(f'SSE: {SSE:.5f}')
pd.DataFrame(dict(予測値=preds[1:] + [forecast]), index=range(1, len(preds)+1)
            ).round(2)

【実行結果】
１期からコツコツ予測値を積み上げています。
そうじゃないと13期の予測ができないので…
SSE は観測値と予測値の誤差の二乗和（二乗誤差）です。

観測値と予測値を比べてみましょう。

### データの可視化
data2.plot(marker='D', ylim=(20, 60))
plt.plot(range(2, 14), preds[1:] + [forecast], marker='D', label='予測値')
plt.legend();

【実行結果】
オレンジの予測値は青い観測値を１期ずらした感じに見えますね。
もしかして「前期と同じ値を予測する」でＯＫかも！？

Python のライブラリで指数平滑化を実践しましょう。
statsmodels を利用します。

🖲️statsmodels の SimpleExpSmoothing()

# statsmodels利用

# 設定
alpha = 0.8

# 指数平滑化の実行
fit = tsa.SimpleExpSmoothing(data2).fit(smoothing_level=alpha, optimized=False)

# 結果の表示
forecast = fit.forecast(1)
print(f'第{forecast.index[0]}期の予測値: {forecast.values[0]:.6f}')
print(f"α: {fit.params['smoothing_level']:.5f}")
print(f'SSE: {fit.sse:.5f}')
fit.level.rename('予測値').to_frame().round(2)

【実行結果】
自作関数の結果と一致しています！

◆

■ 指数平滑化で１期先予測：重みを推定
statsmodels の指数平滑化は重み $${\alpha}$$ の最適値も計算してくれます。
早速、予測してみましょう。

# statsmodels利用 αの最適値を推定

# 指数平滑化の実行
fit = tsa.SimpleExpSmoothing(data2).fit(optimized=True)

# 結果の表示
forecast = fit.forecast(1)
print(f'第{forecast.index[0]}期の予測値: {forecast.values[0]:.6f}')
print(f"α: {fit.params['smoothing_level']:.5f}")
print(f'SSE: {fit.sse:.5f}')
fit.level.rename('予測値').to_frame().round(2)

【実行結果】
$${\alpha \approx 0.69}$$ となりました。
13 期の予測値は $${\alpha=0.8}$$ のときよりも小さくなっています。

可視化しましょう。

### データの可視化 statsmodelsのα最適値
data2.plot(marker='D', ylim=(20, 60))
plt.plot(range(2, 14), fit.level, marker='D', label='予測値')
plt.legend();

【実行結果】
$${\alpha=0.8}$$ のときと大差ない感じです。

自作関数を使って重み $${\alpha}$$ の最適値を探索することもできます。
$${\alpha = [0.0, 0.1, 0.2, \ldots, 0.9, 1.0]}$$ を探索して二乗誤差 SSE が最小になる $${\alpha}$$ を見つけます。

### 誤差平方和の最小となるαの探索 p.277

# 設定
alphas = np.arange(0, 1.1, 0.1)  # 探索するαの値
forecasts, SSEs = [], []         # リストの初期化

# 指数平滑化を用いた予測の実行
for alpha in alphas:
    forecast, SSE, pred_list = exponential_smoothing_pred(data2['保釈率'], alpha)
    forecasts.append(forecast)
    SSEs.append(SSE)

# 結果の表示
(pd.DataFrame({'α': alphas, '誤差平方和': SSEs, '予測値': forecasts})
 .sort_values(['誤差平方和'])
 .reset_index(drop=True)
).round(2)

【実行結果】
$${\alpha=0.7}$$ のとき二乗誤差（誤差平方和）最小となりました。

理解チェック指数平滑化 p.279

木材の出荷量の時系列データを利用して、重み $${\alpha=0.7}$$ の指数平滑化で１期先予測を行います。

コードを淡々と書きます。
テキストのデータをお借りします。

### 指数平滑化 ブナの丸太の出荷量 p.279 表7.3.3

# データの登録
data3 = pd.DataFrame(
    {'採伐量': [2206, 2406, 2259, 2407, 2718, 2267, 2089, 1868, 1778, 1577,
               1486, 1999, 1059, 1122, 1034, 960, 938, 854, 767, 805]},
     index=range(1, 21))
data3.index.name = '時点t'

# 結果の表示
data3

【実行結果】

時系列推移を可視化します。

# 可視化
data3.plot(marker='D', ylim=(0, 2900));

【実行結果】
減少傾向が見られます。
時点 $${2, 5, 12, 14}$$ のときに増加が見られます。

指数平滑化による１期先予測を行います。
まずは自作関数利用です。

### 指数平滑化

# 関数利用
# 設定
alpha = 0.7
# 指数平滑化の実行
forecast, SSE, preds = exponential_smoothing_pred(data3['採伐量'], alpha)
# 結果の表示
print(f'第{data3.index.max()+1}期の予測値: {forecast:.6f}')
print(f'SSE: {SSE:.5f}')

# 回答の表を作成(t=20時点まで)
data3_pred = data3.copy()
data3_pred['予測値'] = [np.nan] + preds[1:]
data3_pred.round(2)

【実行結果】
予測値は $${804}$$。減少傾向を反映している感じです。

観測値と予測値を重ねて可視化しましょう。

### データの可視化
data3.plot(marker='D', ylim=(0, 2900))
plt.plot(range(2, 22), preds[1:] + [forecast], marker='D', label='予測値')
plt.legend();

【実行結果】
観測値の増加・減少の影響をかなり受けている感じです。
１期前の観測値が予測値にスライドしているような印象です。

statsmodels を利用して重み $${\alpha}$$ の推定も行います。

🖲️statsmodels の SimpleExpSmoothing()

# statsmodels利用 αの最適値を推定

# 指数平滑化の実行
fit = tsa.SimpleExpSmoothing(data3).fit()
# 結果の表示
forecast = fit.forecast(1)
print(f'第{forecast.index[0]}期の予測値: {forecast.values[0]:.6f}')
print(f"α: {fit.params['smoothing_level']:.5f}")
print(f'SSE: {fit.sse:.5f}')
fit.level.rename('予測値').to_frame().round(2)

【実行結果】
$${\alpha}$$ の推定値は $${0.732}$$、21 期の予測値は $${803}$$ です。

可視化しましょう。

### データの可視化 statsmodelsのα最適値
data3.plot(marker='D', ylim=(0, 2900))
plt.plot(range(2, 22), fit.level, marker='D', label='予測値')
plt.legend();

【実行結果】

➡️ Section 7.4 自己回帰 AR(1) モデル

自己回帰モデルはある時点の観測値が、重み付けした過去時点の観測値にホワイトノイズを加えたもの、という考え方をしています。
過去の「自分」を変数にして回帰分析する感じなので「自己回帰」。

AR(1) モデルは「１期前の観測値にホワイトノイズを加える」自己回帰（Auto Regression）です。
カッコの中の数値は「次数」と呼ばれ、AR(1) の場合は「次数１」です。
ちなみに AR(2) や AR(3) などの次数の自己回帰モデルもあります！

テキストの「時系列 AR(1) モデルの定義」を引用します。

時系列データを $${\{\cdots, x(t-3), x(t-2), x(t-1), x(t)\}}$$ としたとき、

$$
x(t) = \alpha_1 \cdot x(t-1) + u(t)
$$

を、自己回帰 AR(1)モデルという。
ただし、$${u(t)}$$ はホワイトノイズ。
このとき、1期先の最適な予測値 $${\hat{x}(t, 1)}$$は、

$$
\hat{x}(t, 1) = \alpha_1 \cdot x(t)
$$

となる。

１期前の観測値に重み＝自己回帰係数 $${\alpha_1}$$ を掛けて、ホワイトノイズを足しています。
シンプルな公式です。

仮想データによる例題

テキストに例題がないので、この記事限定の仮想データを作成して、AR(1) モデルを構築しましょう。

■ 分析の準備
AR(1) モデルの数式 $${x(t)=0.72 \cdot x(t_1) + ホワイトノイズ}$$ で仮想データを作成します。
ホワイトノイズは標準正規分布乱数です。

### 時系列データの作成

## 設定
T = 100                               # 時間の数
rng = np.random.default_rng(seed=42)  # 乱数生成器

## データの生成
# ノイズの生成：標準正規分布乱数
noise = rng.standard_normal(size=T)
# 観測値obsの作成：自己回帰 x(t) = 0.72 x(t-1) + noise(t)
obs = [noise[0]]
for t in range(1, T):
    obs.append(0.72 * obs[t-1] + noise[t])

## データのまとめ
# データフレーム化
data4 = pd.DataFrame(
    {'測定値': obs},
    index=pd.date_range(start='2016-01-01', periods=100, freq='MS'))

## データの可視化
data4.plot(figsize=(7, 4), xlabel='年月', grid=True);

【実行結果】
2016 年 1 月から 100 か月の時系列データです。

■ モデルの構築
AR(1) モデルを構築します。
Python の statsmodels ライブラリの ARIMA を利用します。
次のように書きます。

model = tsa.ARIMA(時系列データ, order=(ARの次数, 0, 0), trend=なし).fit()

### AR(1)モデルの構築

fitted_model = tsa.ARIMA(data4, order=(1, 0, 0), trend='n').fit()
fitted_model.summary()

【実行結果】
こちらはモデルのサマリーです。
中段の「ar.L1」行の「coef」が推定した自己回帰係数 $${0.7319}$$ です。
データ生成に用いた $${\alpha=0.72}$$ に近似しています。
「sigma2」はホワイトノイズの分散 $${0.5988}$$ です。
データ生成に用いた $${\sigma^2=1^2}$$ と少々離れている感じです。

■ 将来予測
AR(1) モデルで予測を行ってみましょう。
12 か月先予測をやってみます。次のように書きます。

forecast = model.get_forecast(steps=予測する時点の数)

### AR(1)モデルで予測

## 予測
# 予測の実行
forecasts = fitted_model.get_forecast(steps=12)
# 予測値（平均）と信頼区間の取得
forecasts_mean = forecasts.predicted_mean
conf_ints = forecasts.conf_int()

## 描画
# 描画領域の設定
plt.figure(figsize=(7, 4))
# 観測値の描画（青色）
plt.plot(data4, label='観測値')
# 予測値（平均）の描画（赤色）
plt.plot(forecasts_mean, color='tab:red', label='予測値')
# 予測値の95%信頼区間の塗りつぶし
plt.fill_between(conf_ints.index, conf_ints.iloc[:, 0], conf_ints.iloc[:, 1], 
                 color='tomato', alpha=0.2, label='95%信頼区間')
# 観測値の最終時点に垂直点線を描画
plt.axvline(data4.index[-1], color='gray', lw=3, ls='--')
# 修飾
plt.xlabel('年月')
plt.grid()
plt.legend();

【実行結果】
赤い線は予測値の平均値、薄赤色のエリアは予測値の 95% 信頼区間です。
観測値のギザギザな動きが削ぎ落とされて「滑らかに」なっています。

時系列データの「その先」を学べる書籍の紹介

📗 実践Data Scienceシリーズ Pythonではじめる時系列分析入門
時系列分析を Python で基礎から応用まで学べる良書です。
時系列分析の学びを進めたい方におすすめしたいです。

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

📘 ChatGPTのひとこと：

ゆるやかな変動をすくい取り、
昨日と今日の延長線の先に、そっと明日を描く。
移ろうデータのなかにも、微かな規則が眠っていることを――
そんな気づきを、今日もまたそっとノートに書きとめて。

今回の写経は以上です。

この記事がシリーズ最終回です。
長い間お付き合いいただきまして、ありがとうございます！

ChatGPTと一緒に締めくくりたいと思います。

「わからない」を乗り越えるたびに、見える世界は少しずつ広がってきました。
きっと今のあなたは、はじめて統計と出会ったあの日よりも、ずっと遠くまで来ています。

データを読む力は、未来をつくる力。
小さな数の動きの裏に、大きなヒントが隠れていることを、あなたはもう知っています。

次に扉を開くのは、どんな世界でしょう？
この先もきっと、あなたらしいリズムで学びを続けていけます。

統計を味方に、未来を描こう。
あなたのストーリーは、データとともにこれからも、静かに紡がれていくでしょう。

おわり

シリーズの記事

前の記事

ブログの紹介

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

１．のんびり統計

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

２．実験！たのしいベイズモデリング１＆２をPyMC Ver.5で

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

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

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

４．楽しい写経ベイズ・Python等

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

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

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

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

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

７．Python機械学習プログラミング実践記

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

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

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

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

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

買ってよかったもの

12万件

私のイチオシ

7.3万件