MATLABで『時系列解析入門』９章

状態空間モデル

モデル

$\ell$ 変量の時系列 $y_n$ を表現するための状態空間モデルをいかに示す.

変数の下に示す数字は次元である.

例１：AR(3)

AR(3)を状態空間モデルの形で表してみる. AR(3)は，

ゆえ，$x_n = (y_n, y_{n-1}, y_{n-1})$ とおいて，

が $x_n = F_n x_{n-1} + G_n v_n$ に対応する. また，

が $y_n = H x_n + w_n$ に対応する. もちろん，$Q = \sigma^2$, $R = 0$ となる.

フィルタ，予測，平滑化

フィルタ，予測，平滑化のアルゴリズムは教科書参照のこと. 以下のコードでは，可能な限り関数名や変数を教科書に合わせている.

ARモデルを用いた予測をするための準備

ARモデルの状態空間表現

上に述べたように，AR(m)モデルは以下のように各行列を設定したモデルとして表現できる ( ただし一意ではない ) .

Kalman filter

予測にしろ平滑化にしろまずは，カルマンフィルタを用いて，

一期先予測
• $x_{n|n-1}$：$y_{n-1}$ までの観測に基づく $x_n$ の状態の予測の期待値
• $V_{n|n-1}$：$x_{n|n-1}$ の分散共分散行列

フィルタ
• $x_{n|n}$：$y_{n}$ が観測されたことによって更新された $x_n$ の期待値
• $V_{n|n}$：$x_{n|n}$ の分散共分散行列

を求める必要がある. 以下で求まる

１期先予測

\begin{align} \large x_{n|n-1} &\large= F_n x_{n-1|n-1} \label{ahead} \\[5pt] \large V_{n|n-1} &\large= F_n V_{n-1|n-1} F_n^{\mathrm{T}} + G_n Q_n G_n^{\mathrm{T}} \notag \end{align}

フィルタ

\begin{align} \large K_n &\large= V_{n|n-1} H_n^{\mathrm{T}} (H_n V_{n|n-1} H_n^{\mathrm{T}} + R_n)^{-1} \label{filter} \\[5pt] \large x_{n|n} &\large= x_{n|n-1} + K_n (y_n - H_n x_{n|n-1}) \notag \\[5pt] \large V_{n|n} &\large= (I - K_n H_n) V_{n|n-1} \notag \end{align}

長期予測

カルマンフィルタのプロセスの終点から，1期先予測を続けていけば長期予測が行える. アルゴリズムは以下となる. 式\eqref{ahead} と見比べれば同じことを繰り返しているだけだとわかる.

\begin{align} \large x_{n+i|n} &\large= F_{n+i} x_{n+i-1|n} \label{superahead} \\[5pt] \large V_{n+i|n} &\large= F_{n+i} V_{n+i-1|n} F_{n+i}^{\mathrm{T}} + G_{n+i} Q_{n+i} G_{n+i}^{\mathrm{T}} \notag \end{align}

初期状態

上のアルゴリズムを回すためには，$x_{0|0}$, $V_{0|0}$ を用意する必要がある. $x_{0|0}$ は，観測値がなにもないときの $x_0 = (\vecseq{y}{0}{-m+1})^{\mathrm{T}}$ の期待値なので $x_{0|0} = (0, \dots, 0)^{\mathrm{T}}$ である ( データから平均は事前に引くと仮定 ) . $V_{0|0}$ は，同様に観測値がなにもないときの $x_0$ の分散共分散行列なので，$\mathrm{Cov}(x_0, x_0)$ であるが，これはすなわち，

\begin{align}\large{ V_{0|0} = \begin{pmatrix} C_0 & C_1 & \cdots & C_{1-m} \\ C_1 & C_0 & \cdots & C_{2-m} \\ \vdots & \vdots & \ddots & \vdots \\ C_{1-m}& C_{2-m}& \cdots & C_{0} \\ \end{pmatrix} \notag }\end{align}

ということである. 以下のプログラムでは，第六章 で作ったautocovARMA函数を使ってYule-Walker方程式を解くことで各 $C_i$ を求める.

ARモデルを用いた食料産業従事者数データの長期予測

使用データ

使用データはblsallfoodデータセットである. 函数にまとめず，逐次的に計算結果を示していこうと思う. このデータの長さ ( 点数 ) は156であるが，予測精度を計るために最初の120点を使う. データを確認しておく.

d = blsallfood(1:120);
plot(d)
xlabel('年'); ylabel('人数')

モデルの引数

ARモデルを当てはめる. \eqref{arSys}と\eqref{arObs}を見ればわかるように， AR次数，AR係数，分散 ( イノベーション分散 $\sigma^2$ ) を与える必要がある. これらは，今までやってきたようにAICを用いて次数を決定する. 結果，AR(15)がAIC最小となった.

M = floor(2*sqrt(length(d))); % TSSSのデフォルト値に倣った. 
method = 1;
[aic, Arcoef, s2, maic] = myAr(d, M, method);
bar(0:M, aic-min(aic), 0.5, 'LineStyle', 'none')
xlabel('AR次数'); ylabel('AIC')
ylim([0 100])

では，引数を入れていく. TSSSパッケージの函数tsmoothと同じ引数を設定する.

% observed data 教科書に合わせてyとした. 
y = d;
y = y - mean(y);

% f : 状態の遷移行列F
m = maic;
f = zeros(m, m);
arcoef = Arcoef.(['order' num2str(m)]);
f(1, :) = arcoef;
f(2:end, 1:end-1) = eye(m-1);

% g : 行列Gn
g = [1; zeros(1, m-1)'];

% h : 行列Hn
h = [1 zeros(1, m-1)];

% q : システムノイズの分散共分散行列Qn
q = s2(maic+1);

% r : 観測ノイズの分散共分散行列R
r = 0.0;

% x00 : 初期状態ベクトル x_{0|0}
x00 = zeros(m, 1);

% V00 : 初期状態の分散共分散行列 V_{0|0}
[C, k] = autocovARMA(arcoef, 0, s2(maic+1), m);
V00 =  zeros(m);
for I = 1:m-1
    tmp = eye(m-I)*C(I+1);
    V00(1:m-I, I+1:end) = V00(1:m-I, I+1:end) + tmp;
end
V00 = V00 + V00' + eye(m)*C(1);

% filterEnd : フィルタの終了時点
filterEnd = 120;

% predictEnd : 予測の終了時点
predictEnd = 156;

% minmax : 観測値の上限及び下限 ( 範囲外は異常値として扱う ) 
% 今回は設定せず. 

フィルタリング

まず，フィルタリング. 可能な限り式\eqref{ahead}, \eqref{filter} と変数の文字を揃えた.

% preallocation
xnn_1 = zeros(size(x00, 1), filterEnd);
Vnn_1 = zeros(length(V00), length(V00), filterEnd);
xnn   = zeros(size(x00, 1), filterEnd);
Vnn   = zeros(length(V00), length(V00), filterEnd);

% step 1
xnn_1(:, 1)     = f * x00;                                        % x_{1|0}
Vnn_1(:, :, 1)  = f*V00*f' + g*q*g';                              % V_{1|0}
Kn              = Vnn_1(:, :, 1) * h' / (h*Vnn_1(:, :, 1)*h'+r);  % K_1  
xnn(:, 1)       = xnn_1(:, 1) + Kn * (y(1) - h*xnn_1(:, 1));      % x_{1|1}
Vnn(:, :, 1)    = (eye(m) - Kn*h) * Vnn_1(:, :, 1);               % V_{1|1}

% step 2 - filterend
for n = 2:filterEnd
    xnn_1(:, n)     = f * xnn(:, n-1);                                    % x_{n|n-1}
    Vnn_1(:, :, n)  = f*Vnn(:, :, n-1)*f' + g*q*g';                       % V_{n|n-1}
    Kn              = Vnn_1(:, :, n) * h' / (h*Vnn_1(:, :, n)*h'+r);      % K_n  
    xnn(:, n)       = xnn_1(:, n) + Kn * (y(n) - h*xnn_1(:, n));          % x_{n|n}
    Vnn(:, :, n)    = (eye(m) - Kn*h) * Vnn_1(:, :, n);                   % V_{n|n}
end

1期先予測を確認しておく. いいかんじ.