気象庁の時系列データをRで遊んでみた

はじめに

最近時系列データの解析に興味があって、いろいろと調べてたけれども、実データでないといまいちピンとこないので、いろいろとRの使いかたを調べつつ遊んでみた。emacsでもessをセッティングしているからいくらでもできるのだけど、せっかくなのでRStudioを使ってみた。

現場ですぐ使える時系列データ分析 ~データサイエンティストのための基礎知識~

作者: 横内大介,青木義充
出版社/メーカー: 技術評論社
発売日: 2014/02/18
メディア: 単行本（ソフトカバー）
この商品を含むブログ (6件) を見る

元データ

気象庁が過去の天気データを公開してくれているので、こちらを利用した。
いろいろとデータの種別やら期間やらあるけど、今回は1877年から2006年までの東京の月平均気温を使っている。

気象庁|過去の気象データ・ダウンロード

とりあえず読み込む

> tokyo <- read.csv("data.csv")
> View(tokyo)
> tokyo
        年月 平均気温
1     1877/1      3.2
2     1877/2      3.6
3     1877/3      6.2
4     1877/4     13.6
5     1877/5     16.5
6     1877/6     22.0
7     1877/7     26.5
8     1877/8     25.9
9     1877/9     21.3
10   1877/10     15.9
11   1877/11      9.6
12   1877/12      5.8
13    1878/1      2.3
(後略)

本来は、csvにコメント部分や他の列もあるけど、そこは先に加工したものを取り込んだ。read.csvでそのへんはカバーできるので実は不要であるが。

とりあえずデータテーブルだと扱いづらいので、平均気温を年ごとに扱えるように分割しておく。
はじめ2次元配列に分割しようとして、なんかうまくいかないと悩んでたけど、3次元配列にしてみたらうまくいったのでよかった。

> splited_tokyo = array(tokyo$平均気温, dim=c(1,12,140))
> splited_tokyo[1, 1:12, 1:3]
      [,1] [,2] [,3]
 [1,]  3.2  2.3  3.2
 [2,]  3.6  2.5  5.4
 [3,]  6.2  7.2  8.0
 [4,] 13.6 11.5 12.6
 [5,] 16.5 18.3 18.0
 [6,] 22.0 20.0 21.4
 [7,] 26.5 26.0 26.1
 [8,] 25.9 24.6 26.6
 [9,] 21.3 22.8 21.3
[10,] 15.9 15.8 15.0
[11,]  9.6  9.7  9.7
[12,]  5.8  5.1  8.0

とりあえずグラフ化してみる

複数の時系列データをグラフ化するなら、ts.plotを使うのが簡単

> ts.plot(splited_tokyo[,,], type ="l")

f:id:grugrut:20170226202340p:plain

140年分のデータをまとめてみると、やっぱり毎年同じように気温は変化する(1月2月が一番寒くて、8月が一番暑い)ということがわかるね。

20年ごとの相関も見てみた。

> pairs(splited_tokyo[1,1:12,c(1,21,41,61,81,101,121,140)])

f:id:grugrut:20170226203601p:plain

細かいところをみると違うところもあるけど、ほぼ一直線上にプロットされており、かなり相関係数は1に近いといえそう(計算して出せよ、というのは正しいがやってない)

過去と比べてみる

温暖化温暖化っていうけど、本当に昔に比べて気温はあがってるのだろうか。

とりあえず時系列データということで、tsで時系列データに変換する

> ts_tokyo <- ts(data=tokyo$平均気温, start=c(1877, 1), frequency=12)
> plot(ts_tokyo)

f:id:grugrut:20170226205023p:plain

なんとなく上昇傾向にある気はするけど、よくわからない

> plot(stl(ts_tokyo, s.window="periodic"))

f:id:grugrut:20170226205313p:plain

トレンドは確かに上昇傾向にあるけど、残差も多いのではなかろうか、これ。ただ、残差はプラスマイナス両方に散らばってるので、全体としては上昇傾向にある気がする。

ARIMAモデルによる推定もしてみた。

> library("forecast")
> arima_tokyo <- auto.arima(ts_tokyo, trace=T, stepwise=F, seasonal=T)
> plot(forecast(arima_tokyo, c(95), h=480))

[f:id:grugrut:20170226210851p:plain]

これは、2050年には1年中16度ぐらいの過ごしやすかちょっと肌寒いぐらいの気温になる、ということであろうか。

* 別アプローチ

1877年から1976年までの平均、標準偏差に対し、1977年から2016年までの値がどれぐらいの範囲にいるかを調べてみた。

>|r|
> m <- apply(splited_tokyo[1,1:12,1:100], 1, function(v){return (mean(v))})[1:12]
> m
 [1]  3.482  4.183  7.277 12.951 17.226 20.815 24.694 26.057
 [9] 22.367 16.403 10.982  5.873
> s <- apply(splited_tokyo[1,1:12,1:100], 1, function(v){return (sd(v))})[1:12]
> s

> ts.plot(splited_tokyo[,,101:140], type ="l")
> lines(m-3*s, type = "l", col="red", lwd="3")
> lines(m-2*s, type = "l", col="yellow", lwd="3")
> lines(m-s, type = "l", col="green", lwd="3")
> lines(m, type = "l", col="blue", lwd="3")
> lines(m+s, type = "l", col="green", lwd="3")
> lines(m+2*s, type = "l", col="yellow", lwd="3")
> lines(m+3*s, type = "l", col="red", lwd="3")