#はじめに
東京工業大学/株式会社Nospare の栗栖です．この記事では
「関数時系列データの主成分分析(1)」
「関数時系列データの主成分分析(2)」
に引き続き，関数時系列データに対する主成分分析 (functional principal component analysis, FPCA) の応用例について紹介します．今回は特に関数時系列データの予測に焦点を当て，Rでの分析例を紹介します．

#関数時系列データの例

以下の図は R のパッケージ ftsaに含まれるオーストラリアで1921年から2015年までに得られた各年齢における出生率のデータです．縦軸が出生数，横軸が年齢(15~49歳)に対応しています．図の中のそれぞれの曲線が各年の出生率に対応しています．出生率のデータ以外にも，各都市における PM10 の観測値(Aue et al.(2015)) や一日の中での株価の変化(Li et al.(2020))など，様々な分野で関数時系列データの応用例が知られています．

出生率データの分析

実際に関数時系列データの予測を行う前にデータにFPCAを適用してその結果を眺めてみましょう．
まずは関数データの平均関数の推定をしてみます．下の図の黒い太線がデータの平均です．

library(ftsa)
library(fda)
#Australian Fertility Data
yearbasis35 <- create.fourier.basis(c(14,49),35)
# harmonic acceleration differential operator
harmLfd = vec2Lfd(c(0,(2*pi/35)^2,0), c(14, 49))
AFRfdPar = fdPar(yearbasis35,harmLfd,
                 1e0) #control smoothness
AFRfd = smooth.basis(15:49,Australiafertility$y,AFRfdPar)
AFRpca = pca.fd(AFRfd$fd,nharm=4)
plot(Australiasmoothfertility,
     col = rainbow(length(Australiasmoothfertility$y[1,])), xlab = "Age",
     ylab = "Count of live birth (per 1,000 females)",
     main = "Australia fertility rates (1921-2015)")
lines(AFRpca$meanfd,lwd=5,col=1)

次に関数時系列データの固有値を推定してみましょう．下の図が推定された固有値を大きい順にプロットした結果です．
この結果を見ると第3主成分まで利用すればデータの情報のロスをかなり抑えつつ分析できそうであるということがわかります．

par(mfrow=c(1,1),mar = c(8, 8, 4, 2))
plot(AFRpca$values[1:8],xlab='component',ylab='eigenvalue',col="black",
     cex.lab=1,cex.axis=1,cex=1)

次に固有関数の推定を行ってみます．一例として第1-第4主成分までの推定値をプロットした結果が下の図です．今回の出生率のデータでは第1，第2主成分が昔と最近の大まかな傾向を捉えています．

# functional principal components
harmfd = AFRpca$harmonics
harmvals = eval.fd(15:49,harmfd)
# plot all 4 FPCs
par(mfrow=c(1,1),mar = c(8, 8, 4, 2))
matplot(15:49,harmvals,xlab='year',ylab='PCs',
        lwd=2,lty=1,cex.lab=1,cex.axis=1,type='l')
legend(40,-0.2,c('FPC1','FPC2','FPC3','FPC4'),col=1:4,lty=1,lwd=2)
title('Fertility Rate Principle Component Functions')

##予測

Hyndman and Shang(2009) では FPCA を用いた関数時系列データの予測方法が提案されています．ここでは彼らのアイデアを簡単に紹介します．

$X=\{X_{t}\}_{t \in \mathbb{Z}}$ を関数時系列データとし，$X_{1},\cdots,X_{T}$ を観測値とします．また $\{\lambda_j\}$, $\{v_j\}$ を $X$ の共分散オペレーター $C$ の固有値，固有関数とします．さらに $\hat{C}$, $\{\hat{\lambda}_j\}_{j=1}^{q}$, $\{\hat{v}_j\}_{j=1}^{q}$ をそれぞれ $C$, $\{\lambda_j\}_{j=1}^{q}$, $\{v_j\}_{j=1}^{q}$ の推定値とします．これらは FPCA を用いて計算できます．観測値 $X_{1},\cdots,X_{T}$ をもとに $X_{T+h}$ を予測する ($h$-step ahead prediction) 方法は以下の通りです：

(Step1) $X_{t}$ の FPC スコア

\begin{align*}
\hat{\boldsymbol{X}}_t &= (\langle X_{t}, \hat{v}_1 \rangle,\cdots,\langle X_{t}, \hat{v}_q \rangle)'= :(X_{t}^{(1)},\cdots,X_{t}^{(q)})',\ t=1,\cdots,T
\end{align*}

を計算．

(Step2) FPCによって得られた $q$-次元のベクトル $\hat{\boldsymbol{X}}_t$, $t=1,\cdots, T$ の各成分を互いに独立な1次元の時系列データと仮定して時系列データの予測手法を適用する．ここでは指数平滑化 (exponential smoothing) による予測を考えます．指数平滑化では一期先の予測は以下のように計算されます．

\begin{align*}
X_{T+1|T}^{(k)} &= \alpha X_{T}^{(k)} + \alpha(1-\alpha)X_{T-1}^{(k)}+\cdots + \alpha(1-\alpha)^{T-1}X_{1}^{(k)},\ \text{($1$ 期先予測)}\\
& \vdots\\
X_{T+h|T}^{(k)} &= \alpha X_{T+h-1|t}^{(k)} + \cdots + \alpha(1-\alpha)^{T+h-1}X_{T+1|T}^{(k)}\\
&\quad + \alpha(1-\alpha)^{T+h}X_{T} + \cdots + \alpha(1-\alpha)^{T+h-1}X_{1}^{(k)}\ \text{($h$ 期先予測)}
\end{align*}

$0<\alpha<1,\ t=p,\cdots,T$. 他にもAue et al.(2015)では $q$ 個の時系列データを独立とせず，$q$ 変量のベクトル値自己回帰 (vetor autoregressive, VAR) モデルを用いて予測をする方法を提案しています．ARモデルを用いた予測については FitAR などの R のパッケージで Durbin-Levinson アルゴリズムを用いて $h$ 期先の予測値 ${X}_{T+h|T}^{(k)}$ を計算することができます．

(Step3) (Step2) で求めた ${X}_{T+h|T}^{(k)}$ と $\{\hat{v}_j\}_{j=1}^{q}$ から，Karhunen-Loeve 展開を用いて最終的な関数データの予測値 $X_{T+h|T}$ を計算：

X_{T+h|T} = \sum_{k=1}^{q}\langle X_{T+h|T}^{(k)}, \hat{v}_j \rangle \hat{v}_j.

[](
Aue et al.(2015)では他にも VAR モデルのラグ $p$ と予測に利用するFPCスコアの数 $q$ の選び方についても議論されています．)

上記の方法を利用してオーストラリアの出生率データを分析してみましょう．1921年から2015年までのデータを利用して2016年から2035年まで最長20期先の予測を行ってみます．予測の方法は以下の2通りの方法を考えます．

1つ目は(Step2)で $k$ 期先予測を行う際に $k$ 期先予測を用いて予測を行う. (multiple prediction)
2つ目は(Step2)で $k$ 期先予測を行う際に1期先予測を繰り返して予測を行う．(iterative prediction)

今回は予測に利用するFPCの数は $3$ として分析を行いました．
下の図は multiple prediction の結果です．灰色の線が予測の計算に利用した1921年-2015年のデータ，赤い曲線が2016年の予測に対応し，それ以降はグラデーションで色を変えて2035年(青曲線)までの予測をプロットしています．

# Plot the historical data in gray
nFPC <-3 #set number of FPC 
plot(Australiasmoothfertility,
     col = gray(0.8), xlab = "Age",
     ylab = "Count of live birth (per 1,000 females)",
     main = "Forecasted fertility rates (2007-2026)")
# Plot the forecasts in rainbow color 
plot(forecast(ftsm(Australiasmoothfertility,
                   order = nFPC), #order: number of FPC used
              h = 20), add = TRUE) 
legend("topright", c("2007", "2026"),
       col = c("red", "blue"), lty = 1)

下の図は iterative prediction の結果です．multiple prediction と同様に赤い曲線が2016年の予測に対応し，2035年(青曲線)までの予測をプロットしています．$k$ 期先の予測を直接行う multiple predction と逐次的に行う iterative prediction で結果に違いが出ることがわかります．


plot(Australiasmoothfertility,
     col = gray(0.8), xlab = "Age",
     ylab = "Count of live birth (per 1,000 females)",
     main = "Forecasted fertility rates (2007-2026)")
# Plot the forecasts in rainbow color 
plot(ftsmiterativeforecasts(
  Australiasmoothfertility, components = nFPC, #components: number of FPC used
  iteration = 20), add = TRUE)
legend("topright", c("2007", "2026"),
       col = c("red", "blue"), lty = 1)

#まとめ
この記事では関数時系列データに対する主成分分析 (FPCA) を利用した予測について紹介しました．今後も関数時系列データ解析の関連テーマについて紹介記事を公開予定です．

株式会社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] Hyndman, R.J. and Shang, H.L. (2009). Forecasting functional time series. Journal of the Korean Statistical Society. 38, 199-201.
[3] Li, D., Robinson, P.M., and Shang, H.L. (2020). Long-range dependent curve time series. Journal of the American Statistical Association. 115, 957-971.