Rで周波数解析や成分分解

今回は、まずこんなデータを考えてみる。
期間は2006/1/4〜2015/11/20。

Rのspec.pgramを使用してピリオドグラムを求める。
spec.pgramはデフォルトでは対数グラフを作成するが、今回は対数化は行わないので、log="no"を指定する。

> spec.pgram(y,log="no")

このピリオドグラムから周波数0.1のところにピークがあることがわかる。つまり、原系列は周波数0.1（周期10）のデータであることがわかる。
このように、時系列データに隠されている周期性を解析できたわけだが、このような手法がスペクトル分析だ（理系学部を卒業した人には、常識中の常識ですね・・・）。
今回はスペクトル分析を、Rのパッケージstatsのspec.pgram関数を用いて行ったということ。

ちなみに、既にわかっていると思うが、原系列のデータは次のようにして作成した周期0.1の三角関数である。

> x<-1:100
> T<-10
> y<-ts(cos(2*pi*x/T))

さて、今度は次のようなデータに対してスペクトル分析を行う。
期間は2006/1/4〜2015/11/20。

> plot(tpxc)

> spec.pgram(tpxc,log="no")

うーん、周波数小さすぎる。
見づらいので対数化してみる。

周波数が小さいのでstlを使用しても「 引数 'x' の評価中にエラーが起きました (関数 'plot' に対するメソッドの選択時): stl(tpxc.ts, s.window = "per") でエラー: 系列は周期的でないか、2 未満の周期しか持ちません」とでそうだが、stlで成分分解を行ってみる。

まずは前処理的なやつ。

> tpxc.ts<-ts(as.numeric(tpxc))

周期2でstl。

> tpxc.ts<-ts(as.numeric(tpxc.ts),frequency=2)
> plot(stl(tpxc.ts,s.window="per"))

周期200でstl。

> tpxc.ts<-ts(as.numeric(tpxc.ts),frequency=200)
> plot(stl(tpxc.ts,s.window="per"))

周期100でstl。

> tpxc.ts<-ts(as.numeric(tpxc.ts),frequency=100)
> plot(stl(tpxc.ts,s.window="per"))

周期240でstl。

> tpxc.ts<-ts(as.numeric(tpxc.ts),frequency=240)
> plot(stl(tpxc.ts,s.window="per"))

周期240でのstlの結果を見ながら考えると、大局観として、例えば次のような直近3ヶ月の見通しを立てることもできるだろう。

・季節要因は、残り1/4周期（凡そ60営業日。つまり凡そ3ヶ月。）かけて上昇するだろう。つまり、12月〜2月までの季節要因は上昇。
・トレンドは、2007年3月と6月頃のピークに近づきつつある。そのため、そろそろ反転の可能性があるが、ピークをつけるまでは上昇する可能性もある。
・残差は、マイナスからプラスに転じようとしている傾向がある。少なくとも増加傾向にあるように見える。
以上の3点から、向こう3ヶ月で、上昇するだろう。

向こう3ヶ月で上昇するという想定を立てられたが、来月は12月。年末さ。もう少し細かな想定を立てたい。
そういうわけで、周期80でstlした結果。

これを見ると、今度は次のようなシナリオを作れるだろう。

・季節要因は、残り半周期（凡そ40営業日。つまり凡そ2ヶ月。）かけて下落するだろう。つまり、12月〜1月の季節要因は下落。
・トレンドは、2007年3月と6月頃のピークに到達することなく既に下落を開始している。トレンドが簡単に変わることは考えづらいので、12月〜1月も下落しそう。
・残差は、上限近辺にありそうで、そろそろ下落しそう。
以上の3点から、向こう2ヶ月で、下落するだろう。

あれれ、周期240でやったときは向こう3ヶ月で上昇しそうというシナリオだったのに、周期80でやったときは向こう2ヶ月で下落しそうというシナリオになった。ということは、12月〜1月で下落して、2月に上昇するのか？。こんな考え方もできるだろう。

もちろん、大雑把なチャートを見てチャーチスト的な適当な？見通しなので、これだけでは駄目だが、形はどうあれ見通し（何らかの想定）を持つことは大切かと。