100年間の気温変化

こんな話を見かけた。
「昔の夏はもっと涼しかっただろ!」と思い、50年前の日本の気温を調べてみたところ、まさかの結果が…! | netgeek
これについて言及していた人がいて、データソースがおかしい、ということのようだ。
気象庁 | 日本の年平均気温
変化がない、と言ってるほうも気象庁を参照している、といっている。どこからデータをとってきたのだか…
これをみて、実際に気象庁からデータを取ってきてプロットしたという人がいた。

 
ここでよくあるのが、気温というデータは典型的な自己回帰がある時系列データなので、単純に線形回帰するのはどうなのか、という話。
 
自分は統計素人なので時系列データの解析は詳しくない。状態空間モデルでrstanを使って適当にやってみる。
あと、上のツイートのデータ、どこから取ってきたのかが分からない。最高気温が39度の年があるが、気象庁のどのリンクをたどって取得できたのかがわからない。
気象庁の ホーム > 各種データ・資料 > 過去の気象データ検索 から、東京・東京を選択して、8月の100年間の平均・最高・最低気温データを使うことにした。
データリンク
 
モデルとしては、平均・最高・最低気温の状態は2階差分で正規分布で決まり、実際の観測気温は適当にコーシー分布からサンプリングした。
3つの気温状態のサンプリングに使う分散パラメータは、平均・最高・最低気温で共通としていたが\hat{R} の収束が悪く、3つ別々にしたら改善したので結局別々にしている。
 
結果としては、1910年までは最高・平均気温の変化はほとんどないが、最低気温は常に上昇している傾向である。
1920年から1940年頃は、最高・平均気温が上昇している。最高気温については、1940年代移行、横ばいであるが、2000年移行また上昇している。
平均気温は緩やかに上昇傾向。
最低気温は最高・平均気温に比べてかなり上昇幅が大きいが、2010年くらいから頭打ち、なのかもしれない。

 

# モデル
data{
  int N_year;
  vector<lower=0, upper=40>[3] Temp[N_year];
}
parameters{
  vector<lower=0, upper=40>[3] S[N_year]; // state 
  real<lower=0> s_s[3]; // state variance
  real<lower=0> s_o[3]; // observation variance
}
model{
  for(i in 3:N_year){
    for(j in 1:3){
      S[i, j] ~ normal(2*S[i-1, j] - S[i-2, j], s_s[j]);
    }
  }
  for(i in 1:3){
    Temp[, i] ~ cauchy(S[, i], s_o[i]);
  }
}
library(rstan)
rstan_options(auto_write = TRUE)
options(mc.cores = parallel::detectCores())

dat <- read.table("temperature.txt", header=F, row=1)
m02 <- stan_model("tokyo02.stan")
standata <- list(N_year=nrow(dat), Temp=dat)
fit <- sampling(m02, standata, iter=1500)
ex <- extract(fit)

s_s <- apply(ex$s_s, 2, median)
s_o <- apply(ex$s_o, 2, median)
S <- mapply(function(z) apply(ex$S[,,z], 2, median), 1:3)
cols <- c("black", "red", "blue")
year <- as.numeric(rownames(dat))
yl <- range(dat)
par(mar=c(4, 4.5, 2, 2), cex.lab=1.5)
matplot(year, S, type="l", col=cols, lwd=3, las=1, lty=1, ylim=yl, xlab="year", ylab="temperature")
for(i in 1:3) points(year, dat[,i], col=cols[i], pch=16, cex=0.5)
abline(h=0:40, lty=3, col=grey(0.7))
abline(v=seq(1800, 2030, by=10), lty=3, col=grey(0.7))
# legend("bottomright", legend=sprintf("%s気温", c("平均", "最高", "最低")), col=cols, pch=15, cex=1.4)
1876 26.6 31.6 21
1877 25.9 30.5 21.1
1878 24.6 28.5 21.2
1879 26.6 31.1 22.5
1880 25.5 29.5 21.4
1881 26.7 31.1 22.7
1882 24.7 28.5 21.2
1883 25.1 29.3 21
1884 24.1 28.1 20.4
1885 25.4 29 22
1886 26.5 31.2 22.5
1887 25.3 29.3 22.2
1888 25.6 29.7 22.2
1889 25.8 30 22.3
1890 25.4 29 22.8
1891 25.5 29.6 22
1892 26.4 31.1 22.5
1893 26.2 30.6 22.7
1894 27 31.5 23.4
1895 25.5 29.7 22.8
1896 25.9 30.4 22.5
1897 25 29.2 21.7
1898 26.1 30.6 22.7
1899 26.1 30.2 22.8
1900 26.1 30.5 22.6
1901 25.1 29.9 21.6
1902 22.9 26.2 20.2
1903 25.7 31 21.6
1904 25.1 29.9 21.2
1905 22.2 25.7 19.4
1906 24.5 28.5 21.1
1907 25.8 29.6 23
1908 25.4 29.5 22.3
1909 25.2 29.7 21.7
1910 24.1 27.8 21.2
1911 25.6 29.7 22.7
1912 25.2 30.1 21.6
1913 23.8 28.8 20
1914 26.4 31.3 22.7
1915 25.7 29.7 22.9
1916 25 29.2 22.1
1917 25 29.4 21.5
1918 26.1 30.6 22.5
1919 25 29.3 21.5
1920 25.7 29.3 23
1921 25.3 29.2 22.2
1922 27.3 32.2 23.6
1923 27.2 31.7 23.6
1924 26.2 30.7 22.4
1925 25.7 29.7 22.5
1926 26.3 30.8 22.5
1927 26.6 31.4 22.8
1928 24.1 27.7 21
1929 27.1 31.4 23.2
1930 26.8 30.6 23.4
1931 26.4 30.4 22.8
1932 26.7 30.6 23.6
1933 27.5 31.6 24.1
1934 26.1 30.4 22.5
1935 24.8 28.2 21.9
1936 26.5 30.8 22.8
1937 28.2 32.8 24.3
1938 26 29.7 23.2
1939 25.8 29.8 22.4
1940 24.9 29 21.4
1941 25.4 29.3 22.5
1942 27 31.3 23.5
1943 27.4 31.9 23.7
1944 27.5 32.1 23.9
1945 26.7 31.3 23.1
1946 26.7 31.4 23.3
1947 28 33.3 23.9
1948 25.4 29.7 22.4
1949 26.6 31.4 22.9
1950 26.2 30.4 23.2
1951 26.7 31.6 23.5
1952 26.8 31.4 23.7
1953 25 28.9 22.4
1954 27 31 24.2
1955 26.3 30.2 23.3
1956 25.4 29.7 22.2
1957 27.3 31.2 24.6
1958 25.8 29.7 23
1959 26.7 30.7 23.7
1960 26.4 30.4 23.4
1961 26.8 31.2 23.7
1962 28.1 32.9 24.4
1963 26.6 30.9 23.6
1964 27.8 32.2 24.3
1965 26.7 30.8 23.1
1966 26.9 30.5 24
1967 28 32 24.8
1968 26.6 30.2 23.4
1969 27.2 31.4 23.8
1970 27.4 31.7 24
1971 26.7 30.6 23.7
1972 26.6 30.7 23.3
1973 28.5 32.7 25.5
1974 27.1 31.1 24.1
1975 27.3 31.5 24
1976 25.1 29.1 21.9
1977 25 28.4 22.3
1978 28.9 33 25.6
1979 27.4 31 24.6
1980 23.4 26.6 20.7
1981 26.2 30 23.4
1982 27.1 30.2 24.6
1983 27.5 31 24.6
1984 28.6 32.4 25.7
1985 27.9 31.6 25.1
1986 26.8 30.4 24.1
1987 27.3 30.8 24.3
1988 27 30.2 24.4
1989 27.1 30.5 24.5
1990 28.6 32.4 25.7
1991 25.5 29.1 22.6
1992 27 30.4 24.3
1993 24.8 28 22.3
1994 28.9 32.9 25.9
1995 29.4 33.7 26.1
1996 26 30 23
1997 27 30.7 23.9
1998 27.2 31 24.4
1999 28.5 32.3 25.8
2000 28.3 32.4 25.4
2001 26.4 30 23.6
2002 28 32.1 24.8
2003 26 29.5 23.1
2004 27.2 31 24.1
2005 28.1 31.8 25.1
2006 27.5 31.1 24.8
2007 29 33 25.9
2008 26.8 30.7 24.1
2009 26.6 30.1 23.8
2010 29.6 33.5 27
2011 27.5 31.2 24.6
2012 29.1 33.1 26.3
2013 29.2 33.2 26
2014 27.7 31.2 24.8
2015 26.7 30.5 23.9
2016 27.1 31.6 23.9