ノンパラメトリック回帰モデルの推定
横浜国立大学/株式会社Nospare リサーチャーの栗栖です．
この記事は前回の記事「ノンパラメトリック回帰モデルの適応的推定（２）」 から引き続き，深層ニューラルネットワークを利用したノンパラメトリック回帰モデルの推定について解説します．今回からは複数回にわたって最近の研究成果「ノンパラメトリック時系列回帰モデルに対する適応的深層学習」Kurisu, Fukami and Koike (2022) について解説を行います．前回まではi.i.d.データに対するノンパラメトリック回帰モデルの推定に関して紹介してきましたが，今回はそれらの結果の時系列データへの拡張について紹介します．
ノンパラメトリック時系列回帰モデル
この研究では以下のようなノンパラメトリック時系列回帰モデルを考えます：
Y_t = m(X_t) + \eta(X_t)v_t,\ t=1,\dots, T.
ここで，$m$ はノンパラメトリック回帰モデルにおける回帰関数，$\eta^2$ は 分散関数，$\eta(X_t)v_t$ は観測誤差に対応します．
以上のモデルに従うデータ $(Y_t,X_t) \in \mathbb{R} \times \mathbb{R}^d$，$t=1,\dots,T$ が観測されるとします．特にこの研究ではデータを用いて回帰関数 $m$ の $[0,1]^d$ 上での値の推定 (即ち，$f_0 := m\mathbf{1}_{[0,1]^d}$) の推定を考えます．
上記のモデルは以下のような，時系列データ分析でよく利用される多くのモデルを含みます．
(例1) AR(p)-ARCH(q) モデル
Y_t = m(Y_{t-1},\dots,Y_{t-p}) + \eta(Y_{t-1},\dots, Y_{t-q})v_t
(例2) 一般化加法的 AR(d) モデル
Y_t = \phi(m_1(Y_{t-1}) + \dots + m_d(Y_{t-d})) + \eta(Y_{t-1},\dots, Y_{t-d})v_t
(例3) Threshold AR(1) モデル
\begin{align*}
Y_t =
\begin{cases}
a_1Y_{t-1} + v_t & \text{if $Y_{t-1} \leq r$}, \\
a_2Y_{t-1} + v_t & \text{if $Y_{t-1} > r$}.
\end{cases}
\end{align*}
(例4) Time-varying AR(p)-ARCH(q) モデル
Y_{t,T} = m\left({t \over T},Y_{t-1,T},\dots,Y_{t-p,T}\right) + \eta\left({t \over T},Y_{t-1,T},\dots, Y_{t-q,T}\right)v_t,\ t=1,\dots,T
深層ニューラルネットワーク
この研究で考える深層ニューラルネットワーク (deep neural network, DNN) について説明しておきます．
$L \in \mathbb{N}$ を隠れ層の数 (or DNNの深さ)，$\mathbf{p}=(d,p_1,\dots,p_{L},1)\in \mathbb{N}^{L+2}$ をDNNの幅パラメータとします．またDNNの活性化関数を$\sigma:\mathbb{R} \to \mathbb{R}$ とします．活性化関数が $\sigma$ でネットワーク構造 $(L,\mathbf{p})$ をもつDNN $f:\mathbb{R}^{d} \to \mathbb{R}$ は以下で定義されます．
f(x) = A_{L+1} \circ \sigma_L \circ A_L \circ \sigma_{L-1} \circ \cdots \sigma_1 \circ A_1(x) .
ここで，$A_{\ell}:\mathbb{R}^{p_{\ell-1}} \to \mathbb{R}^{p_\ell}$ は $A_{\ell}(x) = W_{\ell}x + \mathbf{b}_\ell$で 表される線形変換です．$A_{\ell}$ の表現に出てくる $W_\ell$ は $p_{\ell-1}\times p_\ell$の重み行列で，$\mathbf{b}_{\ell} \in \mathbb{R}^{p_\ell}$ はシフトベクトル，$\sigma_\ell:\mathbb{R}^{p_\ell} \to \mathbb{R}^{p_\ell}$ は成分ごとに活性化関数を並べたもの $\sigma_\ell=(\sigma(z_1),\dots,\sigma(z_{p_\ell}))'$ です．活性化関数のクラスとしては $C$-Lipschitz なものを考えます．これはある正の定数 $C>0$ が存在して任意の $x_1, x_2 \in \mathbb{R}$ に対して $|\sigma(x_1) - \sigma(x_2)| \leq C|x_1 - x_2|$ が成り立つようなクラスで，rectifier linear unit (ReLU) 活性化関数 $\sigma(x) = \max\{x,0\}$ や sigmoid 活性化関数 $\sigma(x) = 1/(1+e^{-x})$ などが含まれます．
さらにDNN $f$ に対して，そのパラメータをまとめたベクトルを
\theta(f) := (vec(W_1)',\mathbf{b}'_1,\dots,vec(W_{L+1})',\mathbf{b}'_{L+1})'
とします．$vec(W)$ は行列 $W$ の成分を縦にならべて列ベクトルにする変換を表します．またDNN $f$ に対して depth($f$) を$f$の隠れ層の数，width($f$) を最大幅 $\max_{0 \leq \ell \leq L}p_\ell$ とします．更に$d$-次元の入力と$1$ 次元の出力をもち，活性化関数 $\sigma$ をもつ DNN の集合を $\mathcal{F}_{\sigma, d,1}$ と書くことにします．以上の記号の下で，本研究で考えるDNNの集合を以下のように定義します：
\begin{align*}
\mathcal{F}_\sigma(L,N,B,F)&:= \{f\mathbf{1}_{[0,1]^d}: f \in \mathcal{F}_{\sigma,d,1}, \text{depth}(f)\leq L, \\
&\quad \quad \quad \text{width}(f)\leq N, \|\theta(f)\|_{\infty} \leq B, \|f\|_{\infty} \leq F\}\\
\mathcal{F}_\sigma(L,N,B,F,S) &:= \{f \in \mathcal{F}_{\sigma}(L,N,B,F): \|\theta(f)\|_0 \leq S\}.
\end{align*}
ここで $x \in \mathbb{R}^p$ に対して $\|x\|_{\infty}=\max_{1 \leq j \leq p}|x_j|$，$\|x\|_0=\sum_{j=1}^p \mathbf{1}_{\{x_j \neq 0\}}$，$\|f\|_{\infty} = \sup_{x \in [0,1]^{d}}|f(x)|$ とします．
$\mathcal{F}_{\sigma}(L,N,B,F,S)$ において $S$ はスパースレベルと呼ばれ，$\mathcal{F}_{\sigma}(L,N,B,F)$ において非ゼロパラメータ数に関して制約を課したDNNのクラスなのでスパース制約DNNと呼ぶことにします．
DNN推定量
本研究ではノンパラメトリック時系列回帰モデルの回帰関数 $f_0 = m\mathbf{1}_{[0,1]^d}$ の DNN 推定量として，以下の2種類の推定量を考えます：
\begin{align*}
\hat{f}_{T,np} &\in \text{argmin}_{f \in \mathcal{F}_{\sigma}(L,N,B,F,S)}\left({1 \over T}\sum_{t=1}^{T}(Y_t - f(X_t))^2\right).\\
\hat{f}_{T,sp} &\in \text{argmin}_{f \in \mathcal{F}_{\sigma}(L,N,B,F)}\left({1 \over T}\sum_{t=1}^{T}(Y_t - f(X_t))^2 + J_T(f)\right).
\end{align*}
ここで，$J_T(f) = \lambda_T\|\theta(f)\|_{\text{clip},\tau_T}$，$\lambda_T>0$，また $\theta \in \mathbb{R}^p$ に対して
\|\theta\|_{\text{clip},\tau} := \sum_{j=1}^p\left({|\theta_j| \over \tau} \wedge 1\right), \tau>0
です．上記の2種類の推定量において，$\hat{f}_{T,np}$ をDNN推定量，$\hat{f}_{T,sp}$ をスパース制約DNN推定量と呼ぶことにします．各推定量のパフォーマンスを測る指標としては以下の汎化誤差 (or 予測誤差)を考えます：$\hat{f} \in \{\hat{f}_{T,np}, \hat{f}_{T,sp}\}$ として
R(\hat{f},f_0) := E\left[{1 \over T}\sum_{t=1}^{T}(\hat{f}(X_t^{\ast}) - f_0(X_t^{\ast}))^2\right].
ここで，$\{X_1^{\ast},\dots,X_T^{\ast}\}$ は $\{X_1,\dots, X_T\}$と独立かつ同じ分布をもつデータです．
データ$X = \{X_t\}_{t=1}^{T}$と観測誤差 $\{v_t\}_{t=1}^{T}$ に対して適当な時系列構造 (exponential $\beta$-mixing + 適当なフィルトレーションに関する可測性) を仮定することで $\hat{f}_{T,np}$，$\hat{f}_{T,sp}$ の汎化誤差のレートを導出することができます．
- 時系列データが exponential $\beta$-mixing (+ 定常) であるための十分条件は Chen and Chen(2000) などいくつかの研究があり，既に挙げた AR(d)-ARCH(d) モデルや一般化加法的 AR(d) モデル, threshold AR モデルなどはその条件を満たす特別な場合であることが知られています．
- $\{X_t\}_{t=1}^{T}$ が exponential $\beta$-mixing の十分条件は $v_t$ が独立同分布で，
(1) 回帰関数 $m$ と $\eta$ が任意の有界集合上での有界性：任意の $K>0$ に対して
\sup_{\|x\|\leq K}|m(x)|<\infty,\ 0<\inf_{\|x\|\leq K}\eta(x) \leq \sup_{\|x\|\leq K}\eta(x)<\infty,
(2) $X_t=(Y_{t-1},\dots,Y_{t-d})'$ の各成分が大きな値をとる場合における制約：ある $M>0$，$c_{m,i}>0$，$c_{\eta,i}>0$，$i=0,\dots, d$ が存在して
\begin{align*}
|m(x)| &\leq c_{m,0} + \sum_{i=1}^{d}c_{m,i}|x_i|,\ \text{for $|x| \geq M$,}\\
\eta(x) &\leq c_{\eta,0} + \sum_{i=1}^{d}c_{\eta,i}|x_i|,\ \text{for $|x| \geq M$,}\\
\sum_{i=1}^{d}&(c_{m,i} + c_{\eta,i}E[|v_1|])<1
\end{align*}
であれば満たされるので，$[0,1]^d$ 上での $m$ の関数形についてはほとんど制約がなく，DNN で近似できる様々な回帰関数 $m$ に対して本研究の手法が適用可能であることがわかります．本研究における主要な結果は以下の2つです．特に以下の結果において時系列データの定常性は必要ない (非定常であっても良い) 点に注意しておきます．
DNN推定量の汎化誤差バウンド
$\mathcal{F} = \mathcal{F}_{\sigma}(L,N,B,F,S)$，$F \geq 1$ とします．このとき，任意の $\rho \in (0,1)$ に対してある正の定数 $C_\rho$ が存在して
R(\hat{f}_{T,np},f_0) \leq \rho \inf_{f \in \mathcal{F}}R(f,f_0) + C_{\rho}F^2{S(L+1)\log ((L+1)(N+1)BT)(\log T) \over T}.
が成り立ちます．
スパース制約DNN推定量の汎化誤差バウンド
$\mathcal{F} = \mathcal{F}_{\sigma}(L,N,B,F)$，$F \geq 1$ とします．さらに正の定数 ($C_L, C_B, C_N, C_\lambda,C_\tau,\nu_0, \nu_1, \nu_2$) に対して
L \leq C_L\log^{\nu_0} T,\ N \leq C_N T^{\nu_1},\ 1 \leq B \leq C_B T^{\nu_2},
また $J_T(f)$ のチューニングパラメータ $\lambda_T$，$\tau_T$ は
\lambda_T = C_{\lambda}F^2{\log ^{3+\nu_0}T \over T},\ \tau_T(L+1)((N+1)B)^{L+1} \leq C_\tau T^{-1}
を満たすとします．このとき，$(\eta, C_L,C_N,C_B,C_{\lambda}, C_\tau,\nu_0,\nu_1,\nu_2)$ に依存するある正の定数 $C$ が存在して
R(\hat{f}_{T,sp},f_0) \leq 6 \inf_{f \in \mathcal{F}}\left(R(f,f_0) + J_T(f)\right) + CF^2\left({1 + \log T \over T}\right).
が成り立ちます．
結果に関する考察
-
この研究では独立同分布のデータの下でのアプローチと異なる証明のアプローチをとることで一般の (非定常な) 時系列データや確率過程のモデルにもDNNの理論が拡張可能であることを示しました．
-
前回の記事で既にスパース制約を課さない推定量 $\hat{f}_{T,np}$ は $f_0$ の属する関数空間に関して適応的な推定にならない ($f_0$ の滑らかさなどのパラメータに依存してDNNの集合のパラメータ$L,N,B,S$を設定しなければ理論上最適な汎化誤差のレートが達成できない) ことを紹介しました．$f_0$ がどのような関数であるかは事前にはわからないため，DNNのチューニングパラメータの設定においてこのことは実用上の問題となります．DNN推定量の汎化誤差に関する上記の結果から，時系列データに対しても同様の問題が起こることがわかります．
-
一方で，スパース制約DNN推定量の汎化誤差に関する上記の結果から，スパース制約を課した推定量はDNNのチューニングパラメータを $f_0$ の情報を用いずに設定でき，さらにスパース制約DNN推定量は ($\log T$ のべき乗の違いを無視すれば) Schmidt-Hieber(2020) などで議論されているような様々な関数に対して理論上最適な推定レートを達成可能であることを示すことができます．この点については次回の記事で解説する予定です．
まとめ
この記事では，Kurisu, Fukami and Koike (2022) における深層ニューラルネットワークを用いたノンパラメトリック時系列回帰モデルの適応的推定について紹介しました．株式会社Nospareでは時系列解析に限らず，統計学の様々な分野を専門とする研究者が所属しております．統計アドバイザリーやビジネスデータの分析につきましては株式会社Nospare までお問い合わせください．
参考文献
[1] Chen, M and Chen, G. (2000) Geometric ergodicity of nonlinear autoregressive models with changing conditional variances. Canadian Journal of Statistics 3, 605-613.
[2] Kurisu, D., Fukami, R. and Koike, Y. (2022) Adaptive deep learning for nonparametric time series regression. arXiv:2207.02546.
[3] Schmidt-Hieber, J. (2020) Nonparametric regression using deep neural networks with ReLU activation function. Annals of Statistics 48, 1875-1897.