日本の平均気温偏差の変動データのモデリングと分析
気象庁で 1898年から現在までの全国の平均気温の気温偏差を公開していたのでStanで時系列モデルを作り、その性能の評価を試みました。
ここでの結果には実際の気候変動の原因となっていると考えられている諸現象の効果は含まれておらず、大まかな傾向のみを取り出しているものであることにご注意ください。
使用データ
若干整形したもの(kion_hensa.csv)
気候偏差とは"観測された月平均気温からから、1971～2000年の30年平均値を差し引いたもの"でこの処理によって観測地点間の差異を平均しているようです。
コード
２つのモデルを試行しました。
model0 (月別成分なし)
model{
real m[N];
for(i in 3:N)
trend[i]~normal(2*trend[i-1]-trend[i-2],s_trend);
for(i in 3:N)
ar[i] ~ normal(c_ar[1]*ar[i-1] + c_ar[2]*ar[i-2], s_ar);
for(i in 1:N)
m[i]<-trend[i]+ar[i];
for(i in 1:N)
X[i]~normal(m[i],s_tot);
}
https://gist.github.com/xiangze/9451645#file-model0-stan
model1 (月別成分あり)
model{
real m[N];
for(i in 3:N)
trend[i]~normal(2*trend[i-1]-trend[i-2],s_trend);
for(i in 3:N)
ar[i] ~ normal(c_ar[1]*ar[i-1] + c_ar[2]*ar[i-2], s_ar);
for(i in 12:N)
mm[i]~normal(-mm[i-1]-mm[i-2]-mm[i-3]-mm[i-4]-mm[i-5]-mm[i-6]-mm[i-7]-mm[i-8]-mm[i-9]-mm[i-10]-mm[i-11],s_mm);
for(i in 1:N)
m[i]<-trend[i]+mm[i]+ar[i];
for(i in 1:N)
X[i]~normal(m[i],s_tot);
}
https://gist.github.com/xiangze/9451645#file-model1-stan
(gist)
周期成分の取り入れ方が
mm[i]~normal(-mm[i-1]-mm[i-2]-mm[i-3]-mm[i-4]-mm[i-5]-mm[i-6]-mm[i-7]-mm[i-8]-mm[i-9]-mm[i-10]-mm[i-11],s_mm);
と少し変わって見えますが、これについては『予測にいかす統計モデリングの基本』のpage 11で詳しく説明されています。RStanで『予測にいかす統計モデリングの基本』の売上データの分析をトレースしてみたでも説明されています。
"Random variable is nan:0, but must not be nan!"というエラーに対しては変数の値の範囲の制約をつけることで解決させました。制約はparameter blockに定義された変数のみに付けることが出来ますが、そうするとその変数は出力されることになってしまいます*1。
結果
計算時間の制約からを1973年〜2013年のデータに関する結果を示します。1898年以降の全データに対しても計算を行いましたが8000 stepsを経てもRhatが収束に至りませんでした。
グラフ作成に関してはStanでポアソン分布の状態空間モデルを参考にさせていただきました。matplotlibを用いたコードはこちらです。
model0 (月別成分なし)
model1 (月別成分あり)
長期的な傾向を示すtrend項からはその性質通り、一貫した変化が見られます。グラフを見た限りでは2つのモデルのtrend項に大きな差はありません。月別の効果mmには(これも当然ですが)１年周期のパターンが現れています。
2つのモデルのどちらが妥当かという問題に関しては 情報量基準によって判別されるべきと考えられます。(これから勉強します 続く...)
考察
一般に言われているように気温が上昇していく傾向をtrendの項の動きから読み取れるようになっています。しかしこの結果がなんらかの物理現象を反映していると主張することは出来ません。
世界的には米国海洋大気庁気候データセンター（NCDC）などで気温データを経緯度にして1度間隔で取得しており、
http://www.ncdc.noaa.gov/ghcnm/v3.php
で公開されているので空間情報をモデルにいれることができます(世界の平均気温の偏差の算出方法)。平均化したデータを用いるよりは既存の情報を入れることでよいモデルになると考えられます。一方で物理モデルなしに階層的なモデリングを行うのは難しいと考えています。
http://www.data.kishou.go.jp/climate/cpdinfo/temp/clc_wld.html
一般には気候変動の研究、予測には物理現象のモデル化と観測データをデータ同化という手法で統合して用いているそうです。
Reference
地球温暖化が止まった？ －近未来の気候変動予測－
http://ipccwg1.restec.or.jp/2013sympo/report/pdf/02_kimoto.pdf
気温偏差場の長期変動
http://www.dpac.dpri.kyoto-u.ac.jp/mukou/meeting/Report/06/sakai.pdf
世界と日本の気候変動(2012年)
http://www.jma-net.go.jp/fukuoka/chosa/report2012/report2012_download/pdf/2012_1.pdf
データ同化入門
以前データ同化についてTokyoRで発表しました。
http://xiangze.hatenablog.com/entry/20120922/1348334113
*1:この制約が結果にもたらす影響に関しては個人的にまだ良く分かっていません
