はじめに

東京工業大学/株式会社Nospare の栗栖です．
この記事は「関数時系列データの主成分分析(1)」の続編です．今回は前回紹介した関数時系列データに対する functional principal component analysis (FPCA) の拡張について解説します．前回紹介したFPCAは各関数データが独立または定常な時系列データである場合に利用可能な方法でした．特に以下では非定常な関数時系列データに対するFPCAとその応用例として非定常関数時系列データの2標本検定と予測に焦点を当てて解説します．この記事の内容は主に最近の研究成果 Kurisu(2021) に基づいています．またこの記事の続編の「関数時系列データの主成分分析(3)」では (定常な) 関数時系列データの予測について具体的なデータ分析例を紹介する予定です．

非定常関数時系列データ

$X_{T} = \{X_{t,T}\}_{t \in \mathbb{Z}}$ は内積 $\langle \cdot, \cdot \rangle$ をもつヒルベルト空間 $H$ に値をとる関数時系列とします．
「関数時系列データの主成分分析(1)」とは異なり，ここでは $X_{T} = \{X_{t,T}\}_{t \in \mathbb{Z}}$ について特定の時刻のまわりでの局所的な性質について調べることを考えます．時間に関して局所的な性質に注目して統計分析を行うため，時系列データは各時刻でその性質が変化する非定常性をもっていても良いというのがこの枠組みの良い点です．例えば以下のようなモデルを考えてみましょう．

X_{t,T} = m_{t/T} + v_{t},

ここで $m_{t/T}$ は $\|m_{t/T}\|= \sqrt{\langle m_{t/T}, m_{t/T} \rangle} <\infty$ を満たす関数，$\{v_{t}\}_{t \in \mathbb{Z}}$ は $E[v_{t}] = 0$ かつ $E\left[\|v_{t}\|^{4} \right]<\infty$ を満たす定常な関数時系列とします．関数 $m_{t/T}$ は $t$ ごとに異なって良く，この場合関数時系列データ $\{X_{t,T}\}$ 自身は各時刻で平均関数 $m_{t/T}$ が変化するため非定常性をもちます．一方で，$t/T \approx u \in (0,1)$ と近似できる場合は局所的に

X_{t,T} \approx X_{t}^{(u)} = m_{u} + v_{t}

という関係が成り立ち，$u$ を固定すれば $\{X_{t}^{(u)}\}_{t \in \mathbb{Z}}$ は各時刻で共通の平均関数 $m_{u}$ をもつ定常な関数時系列とみなすことができます．このように局所的に定常な関数時系列で近似できる場合，定常な関数時系列データに対する統計分析のアイデアを拡張することで非定常関数時系列データの統計分析が可能になります.

非定常な関数時系列データの例として，特定の都市における年間気温データ (Aue and van Delft(2020)) や一日単位の金融資産価格データなどが知られています．

関数主成分分析(FPCA)

ここでは非定常関数時系列 $X_{T} = \{X_{t,T}\}_{t \in \mathbb{Z}}$ において $t/T \approx u \in (0,1)$ と近似できる場合に注目し，共分散オペレーター $C_{t/T} \approx C_{u}=E[X_{t}^{(u)} \otimes X_{t}^{(u)}]=E[\langle X_{t}^{(u)}, \cdot \rangle X_{t}^{(u)}]$ を推定することを考えてみます．さらに，前回の記事と同様に共分散オペレーター $C_{u}$ は以下のように特異値分解できることがわかります：

C_{u}(x) = \sum_{j=1}^{\infty}\lambda_{u,j} \langle v_{u,j}, x \rangle v_{u,j},\ x \in H.

ここで $\{\lambda_{u,j}\}$, $\{v_{u,j}\}$ はそれぞれ $C_{u}$ の固有値(eigenvalue)，固有関数(eigenfunction)と呼ばれ，特に $\lambda_{u,1} \geq \lambda_{u,2} \geq \dots$，

C_{u}(v_{u,j})=\lambda_j v_{u,j},\quad \|v_{u,j}\|=1,\quad \langle v_{u,j}, v_{u,k} \rangle=0\ (j \neq k)

を満たします．

関数時系列 $X_{T} = \{X_{t,T}\}_{t \in \mathbb{Z}}$ について，$X_{1,T},\dots,X_{t,T}$ が観測値として得られているとします．このとき $C_{u}$ は

\hat{C}_{u} = {1 \over Th}\sum_{t=1}^{T}K\left({u-t/T \over h}\right) (X_{t,T} - \bar{X}^{(u)}_{T}) \otimes (X_{t,T} - \bar{X}^{(u)}_{T})

で推定できます．ここで，

\bar{X}^{(u)}_{T} = {1 \over Th}\sum_{t=1}^{T}K\left({u-t/T \over h}\right)X_{t,T}

は $t/T \approx u \in (0,1)$ と近似できるような局所的な時刻における標本平均，
$h$ はバンド幅と呼ばれ，$T \to \infty$ で $h \to 0$ となる定数列，
$K$ はカーネル関数と呼ばれ, 特に $\int_{\mathbb{R}}K(x)dx=1$ を満たします．
実際，$u \in (0,1)$ に対して適当な条件のもとで一致性をもつことが示せます：

\|\hat{C}_{u}-C_{u}\|_{\mathcal{L}} = O_{p}\left(h + {1 \over Th} + {1 \over Th^{2}}\right).

($\|\cdot\|_{\mathcal{L}}$ は $H$ 上の有界線形作用素の作用素ノルム，$\hat{C}_{u}$ と $C_{u}$ の近さを測る指標)
共分散オペレーター $C_{u}$ が推定できると $\lambda_{u,j}$, $v_{u,j}$ も推定することができます．$\{\hat{\lambda}_{u,j}\}$, $\{\hat{v}_{u,j}\}$ をそれぞれ $\hat{C}_{u}$ の固有値，固有関数とすると，適当な条件のもとで

\begin{align}
\max_{1 \leq j \leq q}|\hat{\lambda}_{u,j} - \lambda_{u,j}| &= O_{p}\left(h + \sqrt{{1 \over Th}} + {1 \over Th^2}\right),\\ 
\max_{1 \leq j \leq q}\|\hat{c}_{u,j}\hat{v}_{u,j} - v_{u,j}\| &= O_{p}\left(h + \sqrt{{1 \over Th}} + {1 \over Th^{2}}\right).
\end{align}

($q \geq 1$ は固定)
であることが示せます．ここで $\hat{c}_{u,j} = \text{sign}(\langle \hat{v}_{u,j},v_{u,j}\rangle)$ です.

前回の記事同様，上記の結果より，$\hat{v}_{u,j}$ で $v_{u,j}$ を推定する場合，向き(sign)については正しく推定できないため，この結果を統計分析に用いる場合，$v_{u,j}$ の向きに依存しない方法を使う必要があることに注意しておきます．

2標本検定

独立な2つの非定常関数時系列データ $\{X_{t,T_{1}}\}_{t \in T_{1}}$, $\{X_{t,T_{2}}\}_{t \in T_{2}}$ (${T_1 \over T_1 + T_2} \to \theta \in (0,1)$, $\min\{T_1, T_2\} \to \infty$) がそれぞれロケーションモデル

\begin{align}
X_{t,T_{1}} &= m_{t/T_{1}}^{(1)} + v_{1,t},\\
X_{t,T_{2}} &= m_{t/T_{2}}^{(2)} + v_{2,t}
\end{align}

に従う場合に

\mathbb{H}_0: \|m_{u}^{(1)} - m_{u}^{(2)}\|=0\ \text{vs.}\ \mathbb{H}_1: \|m_{u}^{(1)} - m_{u}^{(2)}\|>0

という検定を考えてみましょう．ただし，$\{v_{1,t}\}$, $\{v_{2,t}\}$ は平均ゼロの独立な(定常)関数時系列，$m_{t/T}^{(1)}$，$m_{t/T}^{(2)} \in H$ は $X_{1,t}$, $X_{2,t}$ の平均関数です．
まずいくつか記号を導入しておきます.
$C_{u,1}$, $C_{u,2}$ をそれぞれ $\{X_{t,T_{1}}\}$, $\{X_{t,T_{2}}\}$ の共分散オペレーターとし，$D_{u} = (1-\theta)C_{u,1} + \theta C_{u,2}$, さらに $\{\eta_{u,j}\}$, $\{\nu_{u,j}\}$ を $D_{u}$ の固有値，固有関数とします．
このとき $C_{u,1}$, $C_{u,2}$ は

\hat{C}_{u,j}={1 \over T_{j}h}\sum_{t=1}^{T_{j}}K\left({u-t/T_{j} \over h}\right) (X_{t,T_{j}}-\bar{X}_{j,T}^{(u)}) \otimes (X_{t,T_{j}}-\bar{X}_{j,T}^{(u)}),\quad  j=1,2

で推定できます. $\hat{\theta} = {T_1 \over T_1 + T_2}$ として $\{\hat{\eta}_{u,j}\}$, $\{\hat{\nu}_{u,j}\}$ を $\hat{D}_{u} = (1-\hat{\theta})\hat{C}_{u,1} + \hat{\theta}\hat{C}_{u,2}$ の固有値，固有関数とします．また

\bar{X}_{j,T}^{(u)}={1 \over T_{j}h}\sum_{t=1}^{T_{j}}K\left({u-t/T_{j} \over h}\right)X_{t,T_{j}},\quad  j=1,2

を $t/T \approx u \in (0,1)$ と近似できるような局所的な時刻における標本平均として，検定統計量を

U_{T_1,T_2}^{(u)}:={T_1 T_2 h \over T_{1} + T_{2}}\|\bar{X}_{1,T}^{(u)} - \bar{X}_{2,T}^{(u)}\|^2 = {T_1 T_2 h \over T_{1} + T_{2}}\sum_{j=1}^{\infty}\langle \bar{X}_{1,T}^{(u)} - \bar{X}_{2,T}^{(u)}, \nu_{u,j} \rangle^{2}

で定義します．このとき適当な条件の下で以下が成り立ちます:

\begin{align}
U_{T_1,T_2}^{(u)}&
\begin{cases}
\stackrel{d}{\to} \sum_{j=1}^{\infty}\eta_{u,j}N_{j}^{2}=: U_{\infty,\infty}^{(u)} & \text{under $\mathbb{H}_0$},\\
\stackrel{p}{\to} \infty &  \text{under $\mathbb{H}_1$}.
\end{cases}
\end{align}

ここで $\{N_{j}\}$ は独立な標準正規分布に従う確率変数列です．検定を実行するにあたり，$U_{\infty,\infty}^{(u)}$ の分位点を推定する必要がありますが，漸近分布 $U_{\infty,\infty}^{(u)}$ は無限和で表されるため厳密にその分位点を計算することができません．実用上は上記の無限和の表現をFPCAを用いて有限次元に還元した検定統計量とその漸近分布

\bar{U}_{T_1,T_2}^{(u)}:={T_1 T_2 h \over T_{1} + T_{2}}\sum_{j=1}^{q}\hat{\eta}_{u,j}^{-1}\langle \bar{X}_{1,T}^{(u)} - \bar{X}_{2,T}^{(u)}, \hat{\nu}_{u,j}\rangle^2 
\begin{cases}
\stackrel{d}{\to} \sum_{j=1}^{q}N_{j}^{2} \sim \chi^{2}(q) & \text{under $\bar{\mathbb{H}}_0$},\\
\stackrel{p}{\to} \infty & \text{under $\bar{\mathbb{H}}_1$}.
\end{cases}

を利用します．ここで $\chi^{2}(q)$ は自由度 $q$ のカイ2乗分布で,

\bar{\mathbb{H}}_0: \max_{1 \leq j \leq q}\langle m_{u}^{(1)} - m_{u}^{(2)}, \nu_{u,j} \rangle^{2} = 0\ \text{vs.}\ \bar{\mathbb{H}}_1: \max_{1 \leq j \leq q}\langle m_{u}^{(1)} - m_{u}^{(2)}, \nu_{u,j} \rangle^{2}>0

です．ここでのポイントはFPCAを用いることにより検定統計量の漸近分布が未知パラメータに依存せず，さらにオリジナルのデータの情報のロスを抑えつつ検定が実行可能になっている点です．以上の内容は Horvath et al.(2013) の非定常な場合への拡張になっています．

予測

関数時系列データ $\{X_{t,T}\}_{t=1}^{T_{1}}$ ($T_{1}<T$) を利用して $k$ 期先予測を考えてみましょう．ここでは $X_{T_{1}+k}$ の代わりに，有限次元で近似した

\tilde{X}_{T_{1}+k,T} = \sum_{j=1}^{q}\langle X_{T_{1},T}, v_{u,j} \rangle v_{u,j}

を予測することを考えます．ここで $\{v_{u,j}\}$ は共分散オペレータ $C_{u} = E[X_{t}^{(u)} \otimes X_{t}^{(u)}]$ の固有関数です．$C_{u}$ は時刻 $T_{1}$ までのデータを利用して以下で推定できます：

\bar{C}_{u} = {2 \over Th}\sum_{t=1}^{T_{1}}K\left({u-t/T \over h}\right) (X_{t,T}-\tilde{X}_{T}^{(u)}) \otimes (X_{t,T}-\tilde{X}_{T}^{(u)}),

ここで，

\tilde{X}_{T}^{(u)} = {2 \over Th}\sum_{t=1}^{T_{1}}K\left({u-t/T \over h}\right) X_{t,T}

です．$\bar{C}_{u}$ の固有関数を $\{\bar{v}_{u,j}\}$ とすると，$X_{T_{1}+k}$ は以下で予測することができます：

\hat{X}_{T_{1}+k,T} = \sum_{j=1}^{q}\langle X_{T_{1},T}, \bar{v}_{u,j} \rangle \bar{v}_{u,j}.

適当な条件の下で以下が成り立つことがわかります：

 \|\hat{X}_{T_{1}+k,T} - \tilde{X}_{T_{1}+k,T}\| = O_{p}\left({k \over T} + h + \sqrt{{1 \over T_{1}h}} + {1 \over T_{1}h^{2}} + {1 \over h}\left|u-{T_{1} \over T}\right|\right).

この結果を見ると，近似誤差の第1項は予測期間が短いほど近似が良くなることを意味しており，直感的にも整合的です．第2-4項は $\{v_{u,j}\}_{j=1}^{q}$ の推定誤差，第5項は $T_{1}/T \approx u$ の近似精度に対応しています．またこの結果は Aue et al.(2015) の方法の拡張になっています．

まとめ

この記事では Kurisu(2021) に基づいて非定常な関数時系列データに対する主成分分析 (FPCA) とそれを利用した2標本検定と予測について解説しました．関数データ分析では他にも分類 (classification) や回帰 (parametric/nonparametric regression) も重要なトピックですが，これらについてはまた別の機会に記事にしようと思います．この記事の最初にも書きましたが，続編の「関数時系列データの主成分分析(3)」では関数時系列データの予測について具体的なデータ分析例を紹介する予定です．

株式会社Nospareには統計学の様々な分野を専門とする研究者が所属しております．統計アドバイザリーやビジネスデータの分析につきましては株式会社Nospare までお問い合わせください．

参考文献
[1] Aue, A., Dubart Norinho, D. and Hormann, S. (2015). On the prediction of stationary functional time series. Journal of the American Statistical Association. 110, 378-392.
[2] Aue, A. and van Delft, A. (2020). Testing for stationarity of functional time series in the frequency domain. Annals of Statistics. 48, 2505-2547.
[3] Horvath, L., Kokoszka, P. and Reeder, R. (2013). Estimation of the mean of functional time series and a two-sample problem. Journal of the Royal Statistical Society Series B. 75, 103-122.
[4] Kurisu, D. (2021). On the estimation of locally stationary functional time series. arXiv:2105.11873.