Statistical Programming

Scalaで自己共分散と自己相関の実装

時系列解析入門 - Stats Hack!!

自己共分散も自己相関も時系列データ元々の分布(関数)の時間軸をずらした時にどれほど元のグラフと相関があるかについて計る指標。自己相関は自己共分散の値を使って求め、自己共分散と違い値の大きさが単位に依存しない(分散と標準偏差みたいな関係。本質的には同じだけど標準化されているかそうでないかの違い)。

自己共分散=((元の関数の各値ー平均)×(ずらした関数の各値ー平均))の平均。

f:id:kyu9q:20140223112931p:plain

def autocovariance(raw:Seq[Double],lag:Int):Double={

  lazy val mean=meanf(raw)

  lazy val lag_raw=raw.drop(lag)

  if(raw.length<=lag) 0  

  else lag_raw.zip(raw).map{x=>(x._1-mean)*(x._2-mean)}.reduce{(a,b)=>a+b}/raw.length

  }

 実装ではまず初めに引数のデータの平均と、引数のlagの分だけ先頭から要素を削除したSeqを各変数に代入している。その変数がlazy(遅延評価)になっているのはその後の行を見たらわかるが、もしlagがデータ数以上だった場合自己共分散は常に0になり計算がもったいないから。んでもしlagが必要以上に大きくなかったら、時間軸をずらしたラグデータと元のデータをzipする。次にmapメソッドで対となる要素ー平均を掛け合わせたものの合計/元のデータ数をする。なぜここで同じ平均を使っているのかというと確かに自己共分散を求めるにはそれぞれのデータの平均、上の式で言えばµt,µsを使う必要があるが本来ラグデータは時間軸をずらしただけのデータであり要素数も平均も元のデータと同じ。実装の都合上やりやすいからあたかも要素を削除しているように見えるけどそれは単純に見かけの話。

次に自己相関。自己相関は自己共分散の値を利用して求めるだけ。

ρ(h)= r(h) / r(0)

自己相関(lag h)=自己共分散(h)/ 自己共分散(0)

 普通のRとかでは自己相関関数(acf)を実行すると20くらいまでそれぞれのlag hに応じた自己相関を作ってくれるから、それも真似てみる。

def autocorrelation(raw:Seq[Double])

:IndexedSeq[Double]={

  val r0=autocovariance(raw,0)

  if (raw.length<20) (0 to raw.length).map(x=>autocovariance(raw,x)/r0)

  else (0 to 20).map(x=>autocovariance(raw,x)/r0)

}