ある日の株価が急落したとします。翌日の株価は「前日の下落ショックの影響」を引きずるかもしれませんが、1週間後にはそのショックの影響はほとんど消えているでしょう。このように、「過去に起きたショック（ノイズ）の影響が、一定期間だけ残り、やがて消える」という現象をモデル化したいとき、MA（Moving Average: 移動平均）モデルが力を発揮します。
MAモデルは、時系列分析における最も基本的なモデルの一つであり、ARモデルと並ぶ重要な構成要素です。MAモデルを理解することで、以下のような場面で活用できます。
- 金融市場のショック分析 — ニュースや政策変更が株価に与える一時的影響の持続期間をモデル化する
- 製造工程の品質管理 — 機械の突発的な異常が後続の製品品質にどの程度影響するかを定量化する
- ARMAモデル・ARIMAモデルへの発展 — ARモデルと組み合わせることで、より柔軟な時系列モデリングが可能になる
本記事の内容
- ARモデルの復習とMAモデルの直感的な理解
- MA(q)モデルの数学的定式化
- 自己共分散と自己相関の導出（q次で打ち切りになる性質）
- ARモデルとの双対性（AR表現↔MA表現）
- 反転可能性条件の意味と導出
- Pythonでのシミュレーションと自己相関の検証
- statsmodelsによるMAモデルのフィッティング
前提知識
この記事を読む前に、以下の記事を読んでおくと理解が深まります。
ARモデルの復習
MAモデルを理解するための出発点として、まずARモデルを簡単に振り返りましょう。AR(p)モデルは、現在の値が「過去の自分自身の値」の線形結合で決まるモデルでした。
$$ x_t = \phi_1 x_{t-1} + \phi_2 x_{t-2} + \cdots + \phi_p x_{t-p} + \varepsilon_t $$
ここで $\varepsilon_t \sim \text{WN}(0, \sigma^2)$ はホワイトノイズです。AR(1)モデル $x_t = \phi_1 x_{t-1} + \varepsilon_t$ を再帰的に展開すると、
$$ x_t = \varepsilon_t + \phi_1 \varepsilon_{t-1} + \phi_1^2 \varepsilon_{t-2} + \phi_1^3 \varepsilon_{t-3} + \cdots $$
のように、過去の全てのノイズ $\varepsilon_{t-k}$ の重み付き和として表現できます。ここで重要なのは、AR(1)モデルでは過去のショックの影響が $\phi_1^k$ の割合で指数的に減衰しながらも、無限に続くということです。
しかし、現実には「ショックの影響が有限の期間だけ残り、その後は完全に消える」というケースも多くあります。例えば、ある工場で機械が一時的に故障した場合、修理後すぐに正常に戻り、数日後にはそのショックの影響は完全にゼロになるかもしれません。このような「ショックの影響が有限期間で打ち切られる」現象を直接モデル化するのがMAモデルです。
ARモデルがショックの影響を間接的に（指数減衰で）表現するのに対し、MAモデルはショックの影響を直接的に（有限個のラグで）表現します。この対比を頭に入れておくと、MAモデルの定式化が自然に理解できるでしょう。
MAモデルの直感
MAモデルのアイデアを一言でまとめると、「現在の値は、現在のショックと過去の有限個のショックの重み付き和で決まる」ということです。
日常的なアナロジーで考えてみましょう。あなたが毎日のコーヒーの満足度を記録しているとします。コーヒーの満足度は、その日のコーヒー豆の品質（今日のショック）に加え、昨日の豆の品質や一昨日の豆の品質の「残り香」にも少し影響されるかもしれません。しかし、1週間前の豆の品質は全く影響しないでしょう。この「過去のショックの影響が有限期間だけ残る」というメカニズムが、MAモデルの核心です。
もう少し技術的に言えば、MAモデルはフィルタとして理解できます。ホワイトノイズ（完全にランダムな信号）を入力し、それに有限インパルス応答（FIR: Finite Impulse Response）フィルタをかけた出力が、MA過程に従う時系列データになります。信号処理の世界では、FIRフィルタは非常に馴染みのある概念であり、MAモデルはまさにその統計学版と言えます。
この直感を踏まえて、次にMAモデルの数学的な定式化を見ていきましょう。
MA(q)モデルの定式化
定義
MA(q)モデルは、現在の値が現在のホワイトノイズと過去 $q$ 個のホワイトノイズの線形結合で表されるモデルです。
$$ x_t = \varepsilon_t + \theta_1 \varepsilon_{t-1} + \theta_2 \varepsilon_{t-2} + \cdots + \theta_q \varepsilon_{t-q} $$
これをまとめて書くと、
$$ \begin{equation} x_t = \varepsilon_t + \sum_{j=1}^{q} \theta_j \varepsilon_{t-j} \end{equation} $$
となります。ここで、$\theta_1, \theta_2, \dots, \theta_q$ はMAパラメータ（MA係数）、$\varepsilon_t \sim \text{WN}(0, \sigma^2)$ はホワイトノイズです。
この式が何を意味しているか、直感的に読み解いてみましょう。$x_t$ は、時刻 $t$ におけるショック $\varepsilon_t$ に加え、1期前のショック $\varepsilon_{t-1}$ の影響が $\theta_1$ 倍、2期前のショック $\varepsilon_{t-2}$ の影響が $\theta_2$ 倍、…、$q$ 期前のショック $\varepsilon_{t-q}$ の影響が $\theta_q$ 倍だけ残っている、と解釈できます。そして、$q+1$ 期以上前のショックの影響は完全にゼロです。
ラグ演算子による表現
ARモデルと同様に、ラグ演算子 $B$（$Bx_t = x_{t-1}$）を使うと、MA(q)モデルは非常にコンパクトに書けます。
$$ x_t = (1 + \theta_1 B + \theta_2 B^2 + \cdots + \theta_q B^q) \varepsilon_t $$
MA多項式を $\Theta(B) = 1 + \theta_1 B + \theta_2 B^2 + \cdots + \theta_q B^q$ と定義すれば、
$$ \begin{equation} x_t = \Theta(B) \varepsilon_t \end{equation} $$
と書けます。これは「ホワイトノイズに多項式フィルタ $\Theta(B)$ を適用した出力が $x_t$ である」という、先ほどのFIRフィルタの解釈をそのまま数式で表しています。
MA(1)モデルの例
最も単純なMA(1)モデルを見てみましょう。
$$ x_t = \varepsilon_t + \theta_1 \varepsilon_{t-1} $$
この式は、「現在の値は、今のショックと1期前のショックの $\theta_1$ 倍の和」を意味します。$\theta_1 > 0$ ならば、正のショックの影響が次の期にも正の方向に残ります。$\theta_1 < 0$ ならば、正のショックの影響が次の期には負の方向に作用します（振動的な振る舞い）。
MAモデルの最も重要な性質は、自己相関が有限のラグで打ち切られることです。次のセクションでこれを数学的に証明しましょう。
自己共分散と自己相関の導出
MAモデルの最も際立った特徴は、自己共分散（自己相関）がラグ $q$ で打ち切りになることです。これを数学的に示していきましょう。
MA(q)モデルの期待値
まず、MA(q)モデルの期待値を計算します。$\varepsilon_t$ はホワイトノイズで $\mathbb{E}[\varepsilon_t] = 0$ ですから、
$$ \mathbb{E}[x_t] = \mathbb{E}\left[\varepsilon_t + \sum_{j=1}^{q} \theta_j \varepsilon_{t-j}\right] = 0 + \sum_{j=1}^{q} \theta_j \cdot 0 = 0 $$
期待値は0です。定数項 $\mu$ を加えたモデル $x_t = \mu + \varepsilon_t + \sum_{j=1}^{q} \theta_j \varepsilon_{t-j}$ を考えれば $\mathbb{E}[x_t] = \mu$ となりますが、以下では簡単のため $\mu = 0$ とします。
自己共分散の導出
自己共分散 $\gamma(h) = \text{Cov}(x_t, x_{t+h})$ を導出します。便宜上 $\theta_0 = 1$ と定義しておくと、MA(q)モデルは
$$ x_t = \sum_{j=0}^{q} \theta_j \varepsilon_{t-j} $$
と書けます。ここから自己共分散を計算していきます。
$$ \gamma(h) = \text{Cov}(x_t, x_{t+h}) = \mathbb{E}[x_t \cdot x_{t+h}] $$
$x_t$ と $x_{t+h}$ をそれぞれ展開して代入すると、
$$ \gamma(h) = \mathbb{E}\left[\left(\sum_{i=0}^{q} \theta_i \varepsilon_{t-i}\right) \left(\sum_{j=0}^{q} \theta_j \varepsilon_{t+h-j}\right)\right] $$
ホワイトノイズの性質 $\mathbb{E}[\varepsilon_s \varepsilon_r] = \sigma^2 \delta_{sr}$（$s = r$ のとき $\sigma^2$、$s \neq r$ のとき $0$）を使います。$\varepsilon_{t-i}$ と $\varepsilon_{t+h-j}$ が同じ時刻のノイズであるためには $t – i = t + h – j$、すなわち $j = i + h$ が必要です。
$j = i + h$ を代入すると、$j$ が $0$ 以上 $q$ 以下である範囲、すなわち $0 \leq i + h \leq q$ つまり $0 \leq i \leq q – h$ の範囲で和をとります。
$$ \begin{equation} \gamma(h) = \begin{cases} \sigma^2 \displaystyle\sum_{i=0}^{q-h} \theta_i \theta_{i+h} & (0 \leq h \leq q) \\ 0 & (h > q) \end{cases} \end{equation} $$
この結果は非常に重要です。ラグ $h$ が $q$ を超えると自己共分散が完全にゼロになるのです。これがMAモデルの最大の特徴であり、ARモデルの自己共分散が指数的に減衰しながらも決してゼロにならないのとは対照的です。
具体例：MA(1)モデルの自己共分散
MA(1)モデル $x_t = \varepsilon_t + \theta_1 \varepsilon_{t-1}$ の場合、$q = 1$ ですので、
$h = 0$（分散）の場合:
$$ \gamma(0) = \sigma^2 \sum_{i=0}^{1} \theta_i^2 = \sigma^2(1 + \theta_1^2) $$
$h = 1$ の場合、和の範囲は $0 \leq i \leq 0$ なので $i = 0$ のみです:
$$ \gamma(1) = \sigma^2 \theta_0 \theta_1 = \sigma^2 \theta_1 $$
$h \geq 2$ の場合:
$$ \gamma(h) = 0 $$
自己相関関数（ACF）
自己相関関数 $\rho(h) = \gamma(h) / \gamma(0)$ を計算します。MA(1)の場合、
$$ \begin{equation} \rho(h) = \begin{cases} 1 & (h = 0) \\ \displaystyle\frac{\theta_1}{1 + \theta_1^2} & (h = 1) \\ 0 & (h \geq 2) \end{cases} \end{equation} $$
ラグ2以降の自己相関が完全にゼロになることに注目してください。これがMAモデルのACFの特徴的なパターンです。
MA(1)モデルのACFの最大値
$\rho(1) = \theta_1 / (1 + \theta_1^2)$ の絶対値の最大値を調べてみましょう。$f(\theta_1) = \theta_1 / (1 + \theta_1^2)$ を $\theta_1$ で微分して0と置くと、
$$ f'(\theta_1) = \frac{(1 + \theta_1^2) – \theta_1 \cdot 2\theta_1}{(1 + \theta_1^2)^2} = \frac{1 – \theta_1^2}{(1 + \theta_1^2)^2} = 0 $$
これより $\theta_1 = \pm 1$ で極値をとり、$|\rho(1)| = 1/2$ が最大です。つまり、MA(1)モデルのラグ1自己相関は $\pm 0.5$ を超えないという性質があります。
一般のMA(q)モデルのACF
一般のMA(q)モデルのACFは、
$$ \begin{equation} \rho(h) = \begin{cases} \displaystyle\frac{\sum_{i=0}^{q-h} \theta_i \theta_{i+h}}{\sum_{i=0}^{q} \theta_i^2} & (1 \leq h \leq q) \\ 0 & (h > q) \end{cases} \end{equation} $$
となります。ACFがラグ $q$ で打ち切りになるこの性質は、実データからMA次数 $q$ を同定する際の最も重要な手がかりとなります。具体的には、標本ACFを計算して「ある特定のラグ以降でACFが有意にゼロと異ならない」場合、そのラグをMA次数の候補とすることができます。
ここまでで、MAモデルの自己相関構造がARモデルと本質的に異なることがわかりました。次に、この違いをより深く理解するために、ARモデルとMAモデルの双対的な関係を見ていきましょう。
ARモデルとMAモデルの双対性
ARモデルとMAモデルは、一見全く異なるモデルに見えますが、実は深いところでつながっています。この双対性を理解することは、ARMAモデルやARIMAモデルを扱う上で非常に重要です。
AR(1)モデルのMA表現
まず、AR(1)モデル $x_t = \phi x_{t-1} + \varepsilon_t$（$|\phi| < 1$）を再帰的に展開してみましょう。
$x_{t-1} = \phi x_{t-2} + \varepsilon_{t-1}$ を元の式に代入すると、
$$ x_t = \phi(\phi x_{t-2} + \varepsilon_{t-1}) + \varepsilon_t = \phi^2 x_{t-2} + \phi \varepsilon_{t-1} + \varepsilon_t $$
さらに $x_{t-2}$ を展開して代入すると、
$$ x_t = \phi^3 x_{t-3} + \phi^2 \varepsilon_{t-2} + \phi \varepsilon_{t-1} + \varepsilon_t $$
この操作を無限に繰り返すと（$|\phi| < 1$ のとき $\phi^k x_{t-k} \to 0$）、
$$ \begin{equation} x_t = \sum_{k=0}^{\infty} \phi^k \varepsilon_{t-k} = \varepsilon_t + \phi \varepsilon_{t-1} + \phi^2 \varepsilon_{t-2} + \cdots \end{equation} $$
これは $\theta_k = \phi^k$ としたMA($\infty$)表現に他なりません。つまり、定常なAR(1)過程は、無限次のMAモデルとして表現できるのです。
MA(1)モデルのAR表現
逆の方向も見てみましょう。MA(1)モデル $x_t = \varepsilon_t + \theta \varepsilon_{t-1}$ から $\varepsilon_t$ について解くと、
$$ \varepsilon_t = x_t – \theta \varepsilon_{t-1} $$
ここで $\varepsilon_{t-1} = x_{t-1} – \theta \varepsilon_{t-2}$ を代入すると、
$$ \varepsilon_t = x_t – \theta(x_{t-1} – \theta \varepsilon_{t-2}) = x_t – \theta x_{t-1} + \theta^2 \varepsilon_{t-2} $$
さらに $\varepsilon_{t-2}$ について同様に代入を繰り返すと、
$$ \varepsilon_t = x_t – \theta x_{t-1} + \theta^2 x_{t-2} – \theta^3 x_{t-3} + \cdots $$
$|\theta| < 1$ のとき、この級数は収束し、
$$ \begin{equation} x_t = \sum_{k=1}^{\infty} (-\theta)^k x_{t-k} + \varepsilon_t \end{equation} $$
つまり、$|\theta| < 1$ のMA(1)過程は、$\phi_k = (-\theta)^k$ としたAR($\infty$)表現を持ちます。
双対性のまとめ
この双対性を表にまとめます。
| 性質 | AR(p) | MA(q) |
|---|---|---|
| 定義 | 過去の自分自身の線形結合 | 過去のノイズの線形結合 |
| ACF | 指数的減衰（打ち切りなし） | ラグ $q$ で打ち切り |
| PACF | ラグ $p$ で打ち切り | 指数的減衰（打ち切りなし） |
| 等価表現 | MA($\infty$)で表現可能 | AR($\infty$)で表現可能 |
| 条件 | 定常性条件 | 反転可能性条件 |
ここで新しく登場したPACF（偏自己相関関数）について補足します。PACFは、ラグ $k$ の偏自己相関を「$x_t$ と $x_{t-k}$ の間の相関から、$x_{t-1}, \dots, x_{t-k+1}$ の影響を除去した相関」として定義されます。AR(p)モデルのPACFはラグ $p$ で打ち切りになります。
この双対性から、ACFとPACFのパターンを見ることで、データがARモデルとMAモデルのどちらに適しているかを判断できます。
- ACFが指数的減衰、PACFがラグ $p$ で打ち切り → AR(p)モデルが適切
- ACFがラグ $q$ で打ち切り、PACFが指数的減衰 → MA(q)モデルが適切
- ACF、PACFともに指数的減衰 → ARMA(p,q)モデルが適切
この双対性を理解した上で、次にMA表現からAR表現への変換が可能であるための条件、すなわち反転可能性条件について詳しく見ていきましょう。
反転可能性条件
反転可能性とは何か
前のセクションで、MA(1)モデルをAR($\infty$)表現に変換する際、$|\theta| < 1$ という条件が必要でした。この条件を反転可能性条件（invertibility condition）と呼びます。
なぜ反転可能性が重要なのでしょうか？実は、同じ自己相関構造を持つMAモデルは一意ではないのです。MA(1)モデルの自己相関を思い出してください。
$$ \rho(1) = \frac{\theta_1}{1 + \theta_1^2} $$
この式で $\theta_1 = \theta$ と $\theta_1 = 1/\theta$ を代入してみましょう。
$\theta_1 = \theta$ のとき: $\rho(1) = \theta / (1 + \theta^2)$
$\theta_1 = 1/\theta$ のとき: $\rho(1) = (1/\theta) / (1 + 1/\theta^2) = (1/\theta) / ((θ^2 + 1)/\theta^2) = \theta / (1 + \theta^2)$
驚くべきことに、$\theta$ と $1/\theta$ はまったく同じ自己相関構造を与えます。つまり、データのACFだけからは $\theta$ と $1/\theta$ のどちらが「正しい」パラメータなのか区別できません。
反転可能性条件の役割
この曖昧さを解消するために、反転可能性条件を課します。$|\theta| < 1$ の方を選ぶことで、
- MA過程をAR($\infty$)で表現できる（パラメータ推定の実用上重要）
- パラメータの一意性が保証される（モデルの解釈が明確になる）
一般のMA(q)モデルの反転可能性条件
一般のMA(q)モデルの反転可能性条件は、MA多項式 $\Theta(z) = 1 + \theta_1 z + \theta_2 z^2 + \cdots + \theta_q z^q$ の全ての根が単位円の外側にあること、と表現されます。
$$ \begin{equation} \Theta(z) = 0 \quad \Rightarrow \quad |z| > 1 \end{equation} $$
これはARモデルの定常性条件（AR多項式 $\Phi(z) = 1 – \phi_1 z – \cdots – \phi_p z^p = 0$ の根が単位円の外側）と対称的な構造を持っています。
MA(1)の場合、$\Theta(z) = 1 + \theta_1 z = 0$ の根は $z = -1/\theta_1$ です。$|z| > 1$ であるためには $1/|\theta_1| > 1$、すなわち $|\theta_1| < 1$ が必要です。これは先ほど直接的に導いた条件と一致します。
MA(2)の場合、$\Theta(z) = 1 + \theta_1 z + \theta_2 z^2 = 0$ の2つの根 $z_1, z_2$ がともに $|z_i| > 1$ を満たす必要があります。ヴィエトの公式から $z_1 z_2 = 1/\theta_2$、$z_1 + z_2 = -\theta_1/\theta_2$ であることを使うと、具体的な条件を導くことも可能です。
反転可能性条件を理解したところで、ここまでの理論をPythonで実際に確認してみましょう。シミュレーションを通じて、MAモデルの自己相関の打ち切り性質を目で見て実感します。
PythonによるMAモデルのシミュレーション
MA(1)モデルのシミュレーション
まず、MA(1)モデルからデータを生成し、自己相関関数を確認します。理論通りにラグ1だけに有意な自己相関が現れ、ラグ2以降はゼロになるかを確認しましょう。
import numpy as np
import matplotlib.pyplot as plt
from statsmodels.tsa.stattools import acf, pacf
np.random.seed(42)
# パラメータ設定
n = 1000
theta1 = 0.7
sigma = 1.0
# ホワイトノイズの生成
eps = np.random.normal(0, sigma, n + 1)
# MA(1)モデルのシミュレーション: x_t = eps_t + theta1 * eps_{t-1}
x = np.zeros(n)
for t in range(n):
x[t] = eps[t + 1] + theta1 * eps[t]
# 自己相関・偏自己相関の計算
acf_vals = acf(x, nlags=20)
pacf_vals = pacf(x, nlags=20)
# 理論的なACF
acf_theory = np.zeros(21)
acf_theory[0] = 1.0
acf_theory[1] = theta1 / (1 + theta1**2)
fig, axes = plt.subplots(2, 2, figsize=(12, 8))
# 時系列プロット
axes[0, 0].plot(x[:200], linewidth=0.8)
axes[0, 0].set_title(f'MA(1) Process (θ = {theta1})')
axes[0, 0].set_xlabel('Time')
axes[0, 0].set_ylabel('x(t)')
# ACF
lags = np.arange(21)
axes[0, 1].bar(lags, acf_vals, width=0.3, color='steelblue', label='Sample ACF')
axes[0, 1].scatter(lags, acf_theory, color='red', zorder=5, s=30, label='Theoretical ACF')
axes[0, 1].axhline(y=1.96/np.sqrt(n), color='gray', linestyle='--', alpha=0.5)
axes[0, 1].axhline(y=-1.96/np.sqrt(n), color='gray', linestyle='--', alpha=0.5)
axes[0, 1].set_title('ACF of MA(1)')
axes[0, 1].set_xlabel('Lag')
axes[0, 1].legend()
# PACF
axes[1, 0].bar(lags, pacf_vals, width=0.3, color='coral')
axes[1, 0].axhline(y=1.96/np.sqrt(n), color='gray', linestyle='--', alpha=0.5)
axes[1, 0].axhline(y=-1.96/np.sqrt(n), color='gray', linestyle='--', alpha=0.5)
axes[1, 0].set_title('PACF of MA(1)')
axes[1, 0].set_xlabel('Lag')
# ヒストグラム
axes[1, 1].hist(x, bins=40, density=True, color='steelblue', alpha=0.7)
axes[1, 1].set_title('Distribution of MA(1)')
axes[1, 1].set_xlabel('x')
axes[1, 1].set_ylabel('Density')
plt.tight_layout()
plt.show()
上のグラフから、MA(1)モデルの重要な特徴を確認できます。
- ACF（右上）: 標本ACF（青いバー）はラグ1で有意な正の値（約0.47）を示し、ラグ2以降は95%信頼区間（灰色の破線）の内側に収まっています。赤い点で示した理論値 $\rho(1) = 0.7 / (1 + 0.49) \approx 0.47$ とよく一致しています。これがMA(1)モデルの「ラグ1で打ち切り」という特徴です。
- PACF（左下）: PACFは指数的に減衰しています。これはARモデルのACFの減衰パターンと双対的な関係にあります。
- 分布（右下）: MA(1)過程のデータは正規分布に近い分布を示しています。これはホワイトノイズが正規分布に従う場合の理論的予測と一致します。
MA(2)モデルのシミュレーション
次に、MA(2)モデルでラグ2まで自己相関が残ることを確認しましょう。
import numpy as np
import matplotlib.pyplot as plt
from statsmodels.tsa.stattools import acf
np.random.seed(42)
# パラメータ設定
n = 1000
theta1 = 0.7
theta2 = -0.4
sigma = 1.0
# ホワイトノイズの生成
eps = np.random.normal(0, sigma, n + 2)
# MA(2)モデルのシミュレーション
x = np.zeros(n)
for t in range(n):
x[t] = eps[t + 2] + theta1 * eps[t + 1] + theta2 * eps[t]
# 標本ACF
acf_vals = acf(x, nlags=20)
# 理論的ACF
gamma0 = sigma**2 * (1 + theta1**2 + theta2**2)
gamma1 = sigma**2 * (theta1 + theta1 * theta2)
gamma2 = sigma**2 * theta2
acf_theory = np.zeros(21)
acf_theory[0] = 1.0
acf_theory[1] = gamma1 / gamma0
acf_theory[2] = gamma2 / gamma0
fig, axes = plt.subplots(1, 2, figsize=(12, 4))
# 時系列プロット
axes[0].plot(x[:300], linewidth=0.8)
axes[0].set_title(f'MA(2) Process (θ₁={theta1}, θ₂={theta2})')
axes[0].set_xlabel('Time')
axes[0].set_ylabel('x(t)')
# ACF
lags = np.arange(21)
axes[1].bar(lags, acf_vals, width=0.3, color='steelblue', label='Sample ACF')
axes[1].scatter(lags, acf_theory, color='red', zorder=5, s=30, label='Theoretical ACF')
axes[1].axhline(y=1.96/np.sqrt(n), color='gray', linestyle='--', alpha=0.5)
axes[1].axhline(y=-1.96/np.sqrt(n), color='gray', linestyle='--', alpha=0.5)
axes[1].set_title('ACF of MA(2)')
axes[1].set_xlabel('Lag')
axes[1].legend()
plt.tight_layout()
plt.show()
print(f"理論的ACF: ρ(1) = {acf_theory[1]:.4f}, ρ(2) = {acf_theory[2]:.4f}")
print(f"標本ACF: ρ(1) = {acf_vals[1]:.4f}, ρ(2) = {acf_vals[2]:.4f}")
MA(2)モデルの結果を見ると、理論通りの振る舞いが確認できます。
- ラグ1とラグ2に有意な自己相関が存在 — $\rho(1)$ は正、$\rho(2)$ は負の値を取っています。$\theta_2 = -0.4 < 0$ なので、2期前のショックが現在の値を逆方向に引っ張る効果がラグ2の負の自己相関として現れています。
- ラグ3以降のACFはゼロ — 95%信頼区間の内側に収まっており、MA(2)の理論的予測と一致します。
- 理論値と標本値の一致 — サンプルサイズ $n = 1000$ では、標本ACFが理論ACFと良い精度で一致しています。
パラメータ $\theta$ の影響の可視化
MAパラメータの値を変えたときに、自己相関がどう変わるかを可視化してみましょう。
import numpy as np
import matplotlib.pyplot as plt
from statsmodels.tsa.stattools import acf
np.random.seed(42)
theta_values = [0.3, 0.7, -0.5, -0.9]
n = 2000
nlags = 15
fig, axes = plt.subplots(2, 4, figsize=(14, 6))
for idx, theta1 in enumerate(theta_values):
# シミュレーション
eps = np.random.normal(0, 1, n + 1)
x = np.array([eps[t + 1] + theta1 * eps[t] for t in range(n)])
# 理論ACF
rho1_theory = theta1 / (1 + theta1**2)
# 時系列
axes[0, idx].plot(x[:200], linewidth=0.6)
axes[0, idx].set_title(f'θ = {theta1}')
axes[0, idx].set_ylim(-5, 5)
if idx == 0:
axes[0, idx].set_ylabel('x(t)')
# ACF
acf_vals = acf(x, nlags=nlags)
axes[1, idx].bar(range(nlags + 1), acf_vals, width=0.3, color='steelblue')
axes[1, idx].axhline(y=rho1_theory, color='red', linestyle='--',
label=f'ρ(1)={rho1_theory:.3f}')
axes[1, idx].axhline(y=1.96/np.sqrt(n), color='gray', linestyle='--', alpha=0.3)
axes[1, idx].axhline(y=-1.96/np.sqrt(n), color='gray', linestyle='--', alpha=0.3)
axes[1, idx].set_ylim(-0.6, 0.6)
axes[1, idx].legend(fontsize=8)
if idx == 0:
axes[1, idx].set_ylabel('ACF')
plt.suptitle('MA(1) with Different θ Values', fontsize=13)
plt.tight_layout()
plt.show()
この比較図から、$\theta$ の値がMA(1)過程に与える影響が明確に読み取れます。
- $\theta > 0$（左2つ）: ラグ1の自己相関が正で、隣り合うデータ点が同じ方向に動きやすい（平滑化された時系列）
- $\theta < 0$（右2つ）: ラグ1の自己相関が負で、隣り合うデータ点が逆方向に動きやすい（振動的な時系列）
- $|\theta|$ が大きいほど時系列の見た目が滑らか（または振動的） — ただし $|\rho(1)| \leq 0.5$ の範囲内に制約されています
いずれの場合もラグ2以降のACFはゼロ近傍に留まっており、MA(1)の理論的性質が確認できます。次に、実際のデータに対してMAモデルをフィッティングする方法を見ていきましょう。
statsmodelsによるMAモデルのフィッティング
シミュレーションデータでの検証
まず、真のパラメータが既知のシミュレーションデータを使い、statsmodelsのARIMAクラスでMAモデルのパラメータが正しく推定できることを確認します。ARIMAモデルにおいてAR次数 $p = 0$、差分次数 $d = 0$ とすることで、純粋なMA(q)モデルが得られます。
import numpy as np
import matplotlib.pyplot as plt
from statsmodels.tsa.arima.model import ARIMA
from statsmodels.graphics.tsaplots import plot_acf, plot_pacf
np.random.seed(123)
# 真のパラメータ
n = 500
theta_true = [0.6, -0.3] # MA(2)モデル
sigma_true = 1.0
# MA(2)データの生成
eps = np.random.normal(0, sigma_true, n + 2)
x = np.array([eps[t+2] + theta_true[0]*eps[t+1] + theta_true[1]*eps[t]
for t in range(n)])
# ARIMA(0,0,2) = MA(2)でフィッティング
model = ARIMA(x, order=(0, 0, 2))
result = model.fit()
print("=" * 50)
print("MA(2) Model Estimation Results")
print("=" * 50)
print(f"True θ₁ = {theta_true[0]:.3f}, Estimated θ₁ = {result.maparams[0]:.3f}")
print(f"True θ₂ = {theta_true[1]:.3f}, Estimated θ₂ = {result.maparams[1]:.3f}")
print(f"True σ² = {sigma_true**2:.3f}, Estimated σ² = {result.params[-1]:.3f}")
print(f"AIC = {result.aic:.2f}")
print(f"BIC = {result.bic:.2f}")
print()
print(result.summary())
推定結果を確認すると、$\theta_1$ と $\theta_2$ の推定値が真の値に近いことがわかります。サンプルサイズ $n = 500$ でも、最尤推定法によるパラメータ推定は十分な精度を持っています。
残差診断
モデルの適切さを評価するために、残差診断を行います。適切なモデルであれば、残差はホワイトノイズに近い性質を示すはずです。
import numpy as np
import matplotlib.pyplot as plt
from statsmodels.tsa.arima.model import ARIMA
from statsmodels.graphics.tsaplots import plot_acf
from statsmodels.stats.diagnostic import acorr_ljungbox
from scipy import stats
np.random.seed(123)
# データ生成（前のセクションと同じ）
n = 500
theta_true = [0.6, -0.3]
sigma_true = 1.0
eps = np.random.normal(0, sigma_true, n + 2)
x = np.array([eps[t+2] + theta_true[0]*eps[t+1] + theta_true[1]*eps[t]
for t in range(n)])
# フィッティング
model = ARIMA(x, order=(0, 0, 2))
result = model.fit()
# 残差の取得
residuals = result.resid
fig, axes = plt.subplots(2, 2, figsize=(12, 8))
# 残差の時系列
axes[0, 0].plot(residuals, linewidth=0.6)
axes[0, 0].axhline(y=0, color='red', linestyle='--', alpha=0.5)
axes[0, 0].set_title('Residuals')
axes[0, 0].set_xlabel('Time')
# 残差のACF
plot_acf(residuals, lags=20, ax=axes[0, 1])
axes[0, 1].set_title('ACF of Residuals')
# Q-Qプロット
stats.probplot(residuals, dist="norm", plot=axes[1, 0])
axes[1, 0].set_title('Q-Q Plot of Residuals')
# 残差のヒストグラム
axes[1, 1].hist(residuals, bins=30, density=True, alpha=0.7, color='steelblue')
x_range = np.linspace(residuals.min(), residuals.max(), 100)
axes[1, 1].plot(x_range, stats.norm.pdf(x_range, residuals.mean(), residuals.std()),
'r-', linewidth=2, label='Normal fit')
axes[1, 1].set_title('Distribution of Residuals')
axes[1, 1].legend()
plt.tight_layout()
plt.show()
# Ljung-Box検定
lb_test = acorr_ljungbox(residuals, lags=[5, 10, 15, 20])
print("\nLjung-Box Test Results:")
print(lb_test)
残差診断の結果から、モデルの適合度を評価できます。
- 残差の時系列（左上）: 残差に明確なパターンや趨勢がなく、ランダムな変動に見えます。これはモデルがデータの構造を適切に捉えていることを示しています。
- 残差のACF（右上）: 残差のACFが全てのラグで95%信頼区間内に収まっており、残差に有意な自己相関が残っていないことが確認できます。
- Q-Qプロット（左下）: データ点がほぼ直線上に並んでおり、残差の正規性が保たれています。
- Ljung-Box検定: p値が有意水準0.05を上回れば、残差にホワイトノイズに有意な自己相関がないことを意味します。
モデル次数の選択
実際の分析では、MA次数 $q$ は未知であり、データから推定する必要があります。AICやBICを使った次数選択の例を示します。
import numpy as np
import matplotlib.pyplot as plt
from statsmodels.tsa.arima.model import ARIMA
import warnings
warnings.filterwarnings('ignore')
np.random.seed(42)
# 真のモデル: MA(2) with θ₁=0.6, θ₂=-0.3
n = 500
theta_true = [0.6, -0.3]
eps = np.random.normal(0, 1.0, n + 2)
x = np.array([eps[t+2] + theta_true[0]*eps[t+1] + theta_true[1]*eps[t]
for t in range(n)])
# MA(0)からMA(6)までのAIC/BICを計算
max_q = 6
aic_values = []
bic_values = []
q_range = range(0, max_q + 1)
for q in q_range:
try:
model = ARIMA(x, order=(0, 0, q))
result = model.fit()
aic_values.append(result.aic)
bic_values.append(result.bic)
except Exception:
aic_values.append(np.nan)
bic_values.append(np.nan)
fig, ax = plt.subplots(figsize=(8, 5))
ax.plot(list(q_range), aic_values, 'o-', color='steelblue', label='AIC', linewidth=2)
ax.plot(list(q_range), bic_values, 's-', color='coral', label='BIC', linewidth=2)
ax.axvline(x=2, color='gray', linestyle='--', alpha=0.5, label='True order (q=2)')
ax.set_xlabel('MA Order (q)')
ax.set_ylabel('Information Criterion')
ax.set_title('Model Selection: AIC and BIC')
ax.legend()
ax.set_xticks(list(q_range))
plt.tight_layout()
plt.show()
best_aic = np.nanargmin(aic_values)
best_bic = np.nanargmin(bic_values)
print(f"AICが最小のq: {best_aic}")
print(f"BICが最小のq: {best_bic}")
print(f"真のq: 2")
AIC/BICによるモデル選択の結果、$q = 2$ が最適な次数として選ばれることが確認できます。AICとBICはともに次数2で最小値を取り、真のモデル次数を正しく同定しています。$q > 2$ では余分なパラメータに対するペナルティが効いて情報量基準が増加しており、過剰適合が抑制されていることがわかります。
ここまでの分析を踏まえて、ARモデルとMAモデルの特徴を総括し、実務での使い分けについてまとめましょう。
ARモデルとMAモデルの実践的な使い分け
ここまでの理論を踏まえて、ARモデルとMAモデルの使い分けの指針を示します。
ACF/PACFによる同定
実データに対してまずACFとPACFを計算し、以下のパターンに照らし合わせます。
import numpy as np
import matplotlib.pyplot as plt
from statsmodels.tsa.arima_process import ArmaProcess
from statsmodels.graphics.tsaplots import plot_acf, plot_pacf
np.random.seed(42)
n = 1000
# AR(2)過程
ar_params = np.array([1, -0.7, 0.2]) # 1 - 0.7*B + 0.2*B^2
ma_params = np.array([1])
ar_process = ArmaProcess(ar_params, ma_params)
ar_data = ar_process.generate_sample(n)
# MA(2)過程
ar_params2 = np.array([1])
ma_params2 = np.array([1, 0.7, -0.3])
ma_process = ArmaProcess(ar_params2, ma_params2)
ma_data = ma_process.generate_sample(n)
fig, axes = plt.subplots(2, 2, figsize=(12, 8))
plot_acf(ar_data, lags=20, ax=axes[0, 0])
axes[0, 0].set_title('ACF of AR(2)')
plot_pacf(ar_data, lags=20, ax=axes[0, 1])
axes[0, 1].set_title('PACF of AR(2)')
plot_acf(ma_data, lags=20, ax=axes[1, 0])
axes[1, 0].set_title('ACF of MA(2)')
plot_pacf(ma_data, lags=20, ax=axes[1, 1])
axes[1, 1].set_title('PACF of MA(2)')
plt.tight_layout()
plt.show()
この比較図は、AR(2)とMA(2)のACF/PACFの違いを明確に示しています。
- AR(2): ACF（左上）は指数的に減衰しながら振動し、PACF（右上）はラグ2で明確に打ち切りになっています。
- MA(2): ACF（左下）はラグ2で打ち切りになり、PACF（右下）は指数的に減衰しながら振動しています。
この対称的なパターンは、ARモデルとMAモデルの双対性の直接的な反映です。
実務での考え方
実務では、ACF/PACFのパターンが教科書通りにきれいに現れるとは限りません。そのような場合は以下のアプローチが有効です。
- ARMAモデルを検討する — ACFもPACFも明確に打ち切りにならない場合、ARMA(p,q)モデルが適切かもしれません。
- 情報量基準で比較する — 複数のモデル候補をAIC/BICで比較し、最適な次数を選択します。
- 残差診断を行う — どのモデルを選んでも、最後に残差がホワイトノイズに近いことを確認することが重要です。
まとめ
本記事では、MA（移動平均）モデルの理論と実装について詳しく解説しました。
- MAモデルの核心: 現在の値が、現在と過去 $q$ 個のホワイトノイズの線形結合で決まるモデルであり、ショックの影響が $q$ 期で打ち切られる
- 自己相関の打ち切り: MA(q)モデルのACFはラグ $q$ で打ち切りになり、これがMAモデルの最大の特徴である
- ARモデルとの双対性: 定常なAR過程はMA($\infty$)で表現でき、反転可能なMA過程はAR($\infty$)で表現できる
- 反転可能性条件: MA多項式の根が全て単位円の外側にあることが必要であり、パラメータの一意性を保証する
- モデル同定: ACFの打ち切りパターンからMA次数を同定し、AIC/BICで確認する
MAモデルを理解することで、ARモデルと組み合わせたARMAモデル、さらに非定常時系列に対応したARIMAモデルへの道が開けます。
次のステップとして、以下の記事も参考にしてください。