【時系列解析】非整数積分 (和分)があれば，パワースペクトルはもっと自由：ARFIMA(0, d, 0)の半歩手前から大きく外れてみる

　今回はパワースペクトル $S(f)$ が， $\beta$ を定数として，

　 $S(f) \propto f^{-\beta}$

の形になる，差分方程式で記述される確率過程を考えます．前の記事
chaos-kiyono.hatenablog.com
で，m階積分 (和分)過程のパワースペクトル（もどき）を計算しました．つまり，

分散 $\sigma^2$ の白色ノイズ $\{x [n]\}_{n=1}^N$ の $m$ 階積分 (和分)過程で，長さ $N$ の時系列 $\{y [n]\}_{n=1}^N$ のパワースペクトル $S_y(f_k)$ を時系列のフーリエ変換を使って計算すると，

　 $\displaystyle S_y(f_k) = \frac{\sigma^2}{\left(2- 2 \cos 2 \pi f_k \right)^{m}} = \frac{\sigma^2}{\left( 2 \sin \pi f \right)^{2m}}$

になるはず．

ということを書きました．「パワースペクトル（もどき）」と表現したのは， $m \ge 0.5$ のとき非定常なので，自己共分散 (自己相関)関数のフーリエ変換としてのパワースペクトルが存在しないからです．でも，有限の長さの時系列サンプルがあって，その時系列のフーリエ変換を使ってパワースペクトルを計算すると，上の関数になるはずです．私は数学者ではないので，定常の範囲とか，そんなこと全く気にしません．ですので，今回は $m \ge 0.5$ の世界も考えます．

　ということで，前の記事では， $m$ は整数でした． $m$ が正のとき， $m$ 階積分 (和分)過程，負の値のとき， $m$ 階差分過程になります．

　これらの過程では，

　 $\displaystyle S(f) = \frac{\sigma^2}{\left( 2 \sin \pi f \right)^{2m}} \approx 4^{-m} \sigma^2 (\pi f)^{-2m} \propto f^{-2m}$ 　

になるので， $S(f) \propto f^{-\beta}$ の形になってます．でも残念なことに， $\beta$ は，２の倍数の整数しかとれません．

　何とか，いろんな実数をとれるようにできないでしょうか

　「そんなの簡単だろ， $m$ を整数から実数に拡張すれば， $\beta$ はどんな値でもとれるだろ」と考えたあなたは，枠に収まらない人間です．

　問題は，対応する差分方程式がどうなるかです．この問題は既に考えられていて，答えがあります．

　 $\begin{eqnarray} (1-L)^m y[n] &=& x[n] \\ y[n] &=& \frac{1}{(1-L)^m}x[n] \\ &=& \sum_{j=0}^{\infty} \frac{ \Gamma (j+m)}{\Gamma (j+1) \Gamma (m)} x[n-j] \end{eqnarray}$ $\displaystyle$

という風に， $\frac{1}{(1-L)^m}$ の部分を展開してやることで， $\{x[n] \}$ の重み付き和として， $y[n]$ の値を計算することができます．

　ARFIMA $(0, d, 0)$ では，ここでの $m$ が， $d$ の値で，

　 $\displaystyle - \frac{1}{2} < d < \frac{1}{2}$

の範囲を考えます．これは，定常がどうとかいうことで範囲が決まっています．

　しかし私には，こんな範囲の制限を守る気はないので， $m$ の値をもっと広い範囲で変えて数値実験してみます．つまり， $\{x [n]\}_{n=1}^N$ を平均０，分散 $\sigma^2$ の白色ノイズとして，次の式を使って $\{y [n]\}_{n=1}^N$ を計算します．

　 $\begin{eqnarray} y[n] &=& \sum_{j=0}^{\infty} \frac{ \Gamma (j+m)}{\Gamma (j+1) \Gamma (m)} x[n-j] \end{eqnarray}$ $\displaystyle$

以下は，ARFIMA $(0, d, 0)$ では，やってはダメな領域の結果ですので，決して真似しないでください．責任はとれません．

出発点の１階の積分過程を確認

　 $m=1$ が，１階の積分過程です．やってみた結果が下の図です．図は，左が $\frac{ \Gamma (j+m)}{\Gamma (j+1) \Gamma (m)}$ の値，中がサンプル時系列，右が200個のサンプル時系列から推定したパワースペクトルの平均です．

問題なく，

　 $S(f) \propto f^{-2}$

になってくれました．積分過程では，上左図のように，和をとる重み $\frac{ \Gamma (j+m)}{\Gamma (j+1) \Gamma (m)}$ は１で一定になります．昔の値を積分 (和分)するので，当然です．

1/fノイズ生成に挑戦

　1/fノイズは， $m=0.5$ ( $d = 0.5$ に相当)とすれば，生成できそうです．やってみた結果が下の図です．

期待通り，

　 $S(f) \propto f^{-1}$

になってくれました．上左図のように，和をとる重み $\frac{ \Gamma (j+m)}{\Gamma (j+1) \Gamma (m)}$ は，単調減少しますが，0には，なかなか，なりません．

失敗例

　私の数値実験では， $m \le 1$ の領域において，期待通り $S(f) \propto f^{-2m}$ になりました．しかし．それよりも大きな値を設定すると，うまくいきませんでした．下の図は， $m=1.2$ の結果です．明らかにずれています．

$m > 2$ の領域では，和をとる重み $\frac{ \Gamma (j+m)}{\Gamma (j+1) \Gamma (m)}$ が，過去に戻るほど，大きな値になっています．これは，さすがに不自然ですし，値が異常に大きくなると思います．

まとめ

　今回は，ARFIMA $(0, d, 0)$ の説明をしようと思いましたが，やめました．ARFIMA $(0, d, 0)$ では想定していない外れの領域の数値実験をしました． $0.5 \le m \le 2$ の領域も，非定常ですが，期待通りのパワースペクトルを示してくれました．

おまけ

　使用したRスクリプトです．

###########################
# 平均をとるサンプル時系列数
N.samps <- 100
##########################
# 時系列の長さ (時系列は奇数長になります)
N <- 1001
# 奇数に変換
N <- round(N/2)*2+1
# 半分
M <- (N-1)/2
#########################
# パラメタの設定
m <- 0.5
sig2 <- 1
sig <- sqrt(sig2)
#########################
# 数値実験
j.max <- N
j <- 0:j.max
coef <- exp(lgamma(j+m)-lgamma(j+1)-lgamma(m))
for(i in 1:N.samps){
　#########################
　# 白色ノイズの生成
　w <- rnorm(N+j.max,0,sig)
　#########################
 if(m == 0){
  x.sim <- w[1:N]
 }else{
  x.sim <- c()
  for(k in 1:N){
   x.sim[k] <- w[(k+j.max):k] %*% coef
  }
 }
　#########################
　# PSDの推定
　tmp <- fft(x.sim)
　PSD <- abs(tmp)^2/N
 if(i == 1){
  PSD.samps <- data.frame(PSD)
 }else{
  PSD.samps <- cbind(PSD.samps,data.frame(PSD))
 }
}
# サンプルの平均
PSD <- rowMeans(PSD.samps)
#########################
# 結果の描画
# モデルのパワースペクトルを与える
PSD.model <- function(f,m,sig2){
 return((f!=0)*sig2/(2*abs(sin(pi*f)))^(2*m))
}
#########################
# f = 0を避けて計算
f <- c((1:M)/N,-(M:1)/(N))
# f = 0のとき，パワースペクトルを0に設定
PSD.model <- c(0,PSD.model(f,m,sig2))
f <- c(0,(1:M)/N,-(M:1)/(N))

par(mfrow=c(1,3),cex.lab=1.4)
plot(0:j.max,coef,"l",col=4,lwd=2,xlab="j",ylab="Weight",xaxs="i",las=1)
abline(h=0,lty=2,col=8)
plot(0:(N-1),x.sim,"l",xaxs="i",col=4,xlab="i",ylab="x[i]",main=paste("Synthetic Time Series: m =",m),lwd=2,las=1)
#
shift.PSD <- mean(PSD[f > 0]/PSD.model[f > 0])
plot(f[f > 0],PSD[f > 0],"l",log="xy",col=4,pch=16,xlim=c(f[2],0.5),xaxs="i",xlab="f",ylab="",main="Power Spectrum",lwd=4,las=1)
lines(f[f > 0],PSD.model[f > 0]*shift.PSD,col=2,lwd=2,lty=2)
legend("topright", legend=c(paste("Target PSD: beta =",2*m)), col=c(2), lty=c(2),lwd=c(2))

axis(1,at=10^(-15:15)%x%(1:9),label=FALSE,tck=-0.015)
axis(2,at=10^(-15:15)%x%(1:9),label=FALSE,tck=-0.015)