R言語所蔵のサンプルデータセットをABC順（大文字）優先で拝見させていただいとります。今回はsunspots、太陽黒点数データです。毎度サンプルデータセットのデータは古いよなどとブー垂れてましたが今回は古いことにも意義があります。3つのデータセットをあわせると1700年以来2014年までのデータが含まれている、と。
※「データのお砂場」投稿順Indexはこちら
遥か古代、私メが高校へ行っていたころ、昼休み毎、天文ドームがあった（しかし据え付けの望遠鏡などはとうに失われていた）部室で先輩が太陽黒点の観測をしてました。立派な仕事だとリスペクトしつつも、「地味な観測だな、黒い点の数を数えるだけじゃん、俺はやりたくね～」などと思ってました。すみません。継続は力なりであります。高校生でもできる地味な観測、その上世界のどこで観測しても同じ結果が得られる筈。1700年、日本でいったら元禄年間、この数年後に、赤穂浪士が本所松坂町に討ち入りの時代のデータから含まれております。
今回のサンプルデータセット
サンプルデータセットの解説ページが以下に。
Monthly Sunspot Numbers, 1749–1983
今回太陽黒点のデータについては、似た名前のサンプルデータセットが３つ存在します。ほれこのように。
上記の解説ページは、一番下の sunspots というデータセットの説明です。毎月のデータ、1749年（将軍様は徳川家重公）から1983年までです。
ひとつのデータセットをロードしたつもりが、ずるずると似た名前の他のデータセットもロードされてくるパターンが多い（前回など）のですが、今回はそんなことはありません。上記３つは全て別々にロードしないとなりません。対象期間が微妙にずれてます。似たようなデータなのですが解説ページも独立に存在します。
なおsunspotsのデータでいうと、1960年まではスイス、チューリッヒの天文台のデータ、それ以降は東京天文台のデータみたいっす。
今回は全体を見通すということで３つ全部をざっと眺めてみたいと思います。
まずは生データ
３つ全てをロードするとこんな感じです。すべて時系列データですが、長さがすべて異なります。また、観測データのソースも異なるみたいです。
最初は、一番短い、年単位の sunspot.year をプロットしてみます。犬公方様（徳川綱吉公）の元禄時代から1988年の「昭和元禄」時代までです。皆さんよくご存じの「太陽黒点の１１年周期」見えてるような気がします。
つづいて、一番長い sunspot.month データをプロットしてみます。上記年単位のプロットが「ローパス・フィルタ」かかってスムースであるのに対して、下記の月単位は毎月の変動分、動きが激しいです。
さて、sunspotsの解説ページには処理例が記載されとりますが、結局時系列プロットでしかありません。ぶっちゃけ、上記の素のプロットにプラスしてタイトルをつけただけ。
周期を推定するために自己相関をとってみる
とりあえず周期を推定するためにacf(自己相関）をとってみます。当然、念頭にあるのは「１１年周期」であります。まずは sunspot.year データから。処理の指定はこんな感じ。
acf(sunspot.year)
結果は以下のようです。汚い字で１１yearsと書いてありますが、10年周期か11年周期か微妙。
さて、もっとも長い sunspot.month データをacfに掛けてみます。ただし「月次データ」なのでデフォルト値でかけると短い周期のLagで打ち切られてしまうので、lag.max=15*12 などとしてしまいました。１年12か月で15年分ね。
acf(sunspot.month, lag.max=15*12)
上の結果が以下です。見事に滑らかなグラフになりました。多数回の蓄積があればこそ。継続は力なりのお陰なんであります。
上記のように、自己相関のピークは10.6~10.7年くらいであるように見えます。次回もsunspot.monthか？