【時系列解析】1次自己回帰過程の記憶（自己相関）はいつ消える？

私の研究室の学生は，私があまりにしつこく1次自己回帰過程の話をするのでうんざりしているかもしれません．さらに，しつこくやっていきます．

1次自己回帰過程

　 $y_{n}=a y_{n-1} + w_n$ 　( $-1 < a < 1$ ， $w_n$ は平均0，分散 $\sigma^2$ の白色ノイズ)

　 $\displaystyle S(f)=\frac{\sigma^2}{1+a^2-2 a \cos (2 \pi f)}=\frac{\sigma^2}{(1-a)^2-2 a( \cos (2 \pi f)-1)}$

です．今回は，このパワースペクトルのグラフの構造と時系列の記憶（自己相関）が消える時間の長さ（ラグ）の関係について説明します．

aの値を変えてグラフを描いたもの（両対数目盛だぞ）が以下です．ここでは，aは1に近いとします．

パワースペクトルの特徴は，低い周波数で平らな丘があり，高い周波数側で下り坂があります．aを1に近づけていく（0.9 -> 0.99 -> 0.999）と，高い周波数側の下り坂の領域が広がり， $f^{-2}$ に比例する構造が見えてきます．

1次自己回帰過程は，短期記憶過程の典型例です．時系列解析の分野では，「記憶」とは「自己相関（自己共分散）」のことです．短期記憶過程では，特徴的な時間の長さがあり，その時間以上では記憶が消える，つまり，自己相関が0とみなせます．では，この自己相関が0とみなせる特徴的な時間の長さを，どのようにすれば見積もることができるのでしょうか．

指数関数的に減衰する関数では，時定数とか，緩和時間とか呼ばれる特徴量があります．
　 $x(t) \propto e^{-t/\tau_0}$
の形のとき，時定数は $\tau_0$ です．肩に乗っている部分が-1になるとして，tを求めれば，それが時定数です．

1次自己回帰過程の自己共分散関数 $C(k)$ は，

　 $\displaystyle C(k)=\frac{\sigma^2}{1-a^2} a^{k} \propto e^{(\log a) k}$

ですので，時定数は

　 $\displaystyle -\frac{1}{\log a}$

です（ $\displaystyle \log a$ はマイナスの値です）．時定数は値がほぼ0とみなせる時間だと，理解している人がいると思います．では，

　 $k > \displaystyle -\frac{1}{\log a}$

で，記憶（自己相関）が消えるのでしょうか．違います．オーダーとしては正しい（10倍も違わない）ですが，無相関とみなせる時間としてもっと正確な範囲は，

　 $\displaystyle k > \frac{2 \pi \sqrt{a}}{1-a}$ ，あるいは， $\displaystyle -\frac{2 \pi}{\log a}$

です．

なぜかを説明するために，パワースペクトルを考えます．パワースペクトルでは，自己相関が0とみなせる低周波数の領域で，平らな丘が現れます（両対数目盛でのグラフだぞ）．平らな丘は，白色ノイズと同じ構造です．これは，1次自己回帰過程のパワースペクトルに見られる低周波数側の平地のことです．ですので，パワースペクトルの低周波数側で，平らな丘に達する周波数を見積もってやれば，記憶が消えるスケールを見積もることができます．

解析的にこの周波数を計算するために，低周波数側の丘の部分と，高周波数側の坂道の近似関数を求めます．近似のために，aは1に非常に近いことを仮定します．さらに，

　 $\displaystyle \cos x = 1 - \frac{x^2}{2!}+ \cdots$

を使います．fは，1/2までの値をとるので，fの2乗まで考えます．

このとき，パワースペクトルは

　 $\displaystyle S(f) \approx \frac{\sigma^2}{(1-a)^2+ 4 \pi^2 a f^2}$

となり，分母の $(1-a)^2$ が， $4 \pi^2 a f^2$ よりも大きいときは，

　 $\displaystyle S(f) \approx \frac{\sigma^2}{(1-a)^2}$

分母の $4 \pi^2 a f^2$ が， $(1-a)^2$ よりも大きいときは，

　 $\displaystyle S(f) \approx \frac{\sigma^2}{ 4 \pi^2 a f^2}$

と近似できそうです．数学的な厳密さとか言われても私は数学者ではないので，数値的に検証します．検証した結果が下の図の右上です．水色の破線（近似）と赤の線（近似なし）が一致しているので，正しそうです．

坂から丘に変化する点を，上の2つの近似式の交点を使って見積もると，

　 $\displaystyle f = \frac{1-a}{2 \pi \sqrt{a}}$

になります．これを，特徴的な周波数 $\displaystyle f_c$ とします．この逆数が特徴的な時間になるので，これを $k_c$ と表せば，

　 $\displaystyle k_c = \frac{2 \pi \sqrt{a}}{1-a}$

となります．今日は疲れたので説明しませんが，

　 $\displaystyle k_c = -\frac{2 \pi}{\log a}$

もほぼ同じ値になります（正直，どうやって導いたか今は思い出せません）．どちらも， $1-a$ でテーラー展開すると，

　 $\displaystyle \frac{2 \pi \sqrt{a}}{1-a} \approx (1-a) + \frac{1}{2} (1-a)^2+\cdots$
　 $\displaystyle -\frac{2 \pi}{\log a} \approx (1-a) + \frac{1}{2} (1-a)^2+\cdots$

のように2乗まで一致するので，ほとんど違いはありません．

上の図の右下の図は自己共分散のグラフです．垂直な灰色破線が時定数の位置，青色破線が $\displaystyle k_c = 1/f_c$ の位置です．灰色破線の時間の長さ（ラグ）だと値が0には見えないけど，青色破線の時間ではほぼ0に見えます．

どのくらいの時間で記憶が消えるか，わかりましたか．その式を覚えるよりも，パワースペクトルの読み方を理解することが大切です．

おまけ

図の作成に使用したRスクリプト

a <- 0.98
par(mfrow=c(2,2))
par(pty="m",cex.axis=1.3,cex.lab=1.5,mar=c(5,5,2,2))
####################
curve(sig2/(1-2*a*cos(2*pi*x)+a^2),log="xy",col=2,lwd=2,n=1000,xlim=c(0.00001,1/2),ylim=c(0.2,10000),xaxs="i",xlab="f",ylab="PSD S(f)",xaxt="n",yaxt="n")
# 対数目盛を描く
axis(1,at=10^(-5:5)%x%(1:9),label=FALSE,tck=-0.02)
tmp<-paste(paste(sep="","expression(10^",-5:5,")"),collapse=",")
v.label<-paste("axis(1,las=1,at=10^(-5:5),label=c(",tmp,"),tck=-0.03)",sep="")
eval(parse(text=v.label))
axis(2,at=10^(-5:5)%x%(1:9),label=FALSE,tck=-0.02)
tmp<-paste(paste(sep="","expression(10^",-5:5,")"),collapse=",")
v.label<-paste("axis(2,las=1,at=10^(-5:5),label=c(",tmp,"),tck=-0.03)",sep="")
eval(parse(text=v.label))
####################
curve(sig2/(1-2*a*cos(2*pi*x)+a^2),log="xy",col=2,lwd=2,n=1000,xlim=c(0.00001,1/2),ylim=c(0.2,10000),xaxs="i",xlab="f",ylab="PSD S(f)",xaxt="n",yaxt="n")
# 対数目盛を描く
axis(1,at=10^(-5:5)%x%(1:9),label=FALSE,tck=-0.02)
tmp<-paste(paste(sep="","expression(10^",-5:5,")"),collapse=",")
v.label<-paste("axis(1,las=1,at=10^(-5:5),label=c(",tmp,"),tck=-0.03)",sep="")
eval(parse(text=v.label))
axis(2,at=10^(-5:5)%x%(1:9),label=FALSE,tck=-0.02)
tmp<-paste(paste(sep="","expression(10^",-5:5,")"),collapse=",")
v.label<-paste("axis(2,las=1,at=10^(-5:5),label=c(",tmp,"),tck=-0.03)",sep="")
eval(parse(text=v.label))
########################
abline(h=sig2/(1-a)^2,col=4,lwd=2,lty=2)
curve(sig2/(pi*x)^2/(4*a),add=TRUE,col=4,lwd=2,lty=2,xlim=c(0.00001,1/2))
#curve(sig2/(pi*x)^2/(4*a)*(pi*x)^2/sin(pi*x)^2,add=TRUE,col=5,lwd=2,lty=2,xlim=c(0.00001,1/2))
abline(v=(1-a)/(2*pi*sqrt(a)),col=4,lwd=2,lty=2)
#abline(v=-log(a)/(2*pi),col=1,lwd=2,lty=2)
####################################################
# 解析的に計算した自己共分散関数
curve(sig2/(1-a^2)*a^x,col=2,lwd=2,lty=1,xlim=c(1,10000),log="x",xaxs="i",xlab="k",ylab="Autocovariance C(k)",xaxt="n")
axis(1,at=10^(-5:5)%x%(1:9),label=FALSE,tck=-0.02)
tmp<-paste(paste(sep="","expression(10^",-5:5,")"),collapse=",")
v.label<-paste("axis(1,las=1,at=10^(-5:5),label=c(",tmp,"),tck=-0.03)",sep="")
eval(parse(text=v.label))
####################################################
# 解析的に計算した自己共分散関数
curve(sig2/(1-a^2)*a^x,col=2,lwd=2,lty=1,xlim=c(1,10000),log="x",xaxs="i",xlab="k",ylab="Autocovariance C(k)",xaxt="n")
axis(1,at=10^(-5:5)%x%(1:9),label=FALSE,tck=-0.02)
tmp<-paste(paste(sep="","expression(10^",-5:5,")"),collapse=",")
v.label<-paste("axis(1,las=1,at=10^(-5:5),label=c(",tmp,"),tck=-0.03)",sep="")
eval(parse(text=v.label))
abline(h=0,col=8,lwd=2,lty=2)
abline(v=-1/log(a),col=8,lwd=2,lty=2)
abline(v=-2*pi/log(a),col=4,lwd=2,lty=2)