不動産価格データを使った相場分析ー世田谷区一棟マンション
このマガジンの一部記事をまとめて、Kindle本を出版しました。Kindle Unlimitedにご加入の方は読み放題で読めます。
GAFAデータサイエンティストが実践する不動産投資術
前回は世田谷区のモデルを作ってみた。地域ごとのモデルまで作ったところで、地域ごとに築年の影響とか、取引時点の年の影響とか、データを見ると若干異なるように見える、と言うところで終わっていた。
今回は前回のモデルをさらに発展させて、地域ごとに築年と取引時点の影響を変化させるようなモデルを考える。これはちょっと大変で、単純な線形モデルではできない。何を使うかというと、階層ベイズモデルを使う。
まずは前回のモデルをベースに考える。前回のモデルは、
価格＝3315 + 50 x 延べ床面積 - 155 x 築年 + 1748 x 構造 ＋ 取引年次効果 ＋ 地域効果
と言うものだった。構造というのは鉄筋コンクリートなどRC構造の場合は１で、それ以外は0と言う変数で、要するにRC構造みたいにガッチリしたマンションは1784万円を足す、というものだ。取引年次効果は、取引された年がいつかによって変わるもので、例えば2010年だったら1051万円マイナスする、というものだ。地域効果は物件がどこになるかによって変わるもので、例えば喜多見だったら-9064万円、玉川だったら8076万円足す、というものだった。
マンション価格モデルのベイズ化
まずこれまでのモデルをベイズモデルで書き換えないといけない。ベイズモデルはRだとRStanというパッケージがあり、これを使う。詳細は参考資料に任せるとして、Markov Chain Monte Carlo (MCMC)を使って乱数を発生させ、モデルフィットをする。モデルの構造をかなり柔軟に変えられる一方、乱数計算なのでモデルフィットには時間がかかる。また、別のファイルにモデルの式を書き出さないといけないくて、それも面倒臭い。
さて、前回モデルはStanで書くとこうなる。
data {
int<lower=0> N;
int<lower=0> K; // year effect index
int<lower=0> J; // chiiki effect index
vector[N] price;
vector[N] floor_size;
vector[N] age;
vector[N] is_rc_pc;
int<lower=1, upper=K> year_effect_index[N]; // 1 = 2005, 2 = 2006, 3 = 2007, etc.
int<lower=1, upper=J> chiiki_effect_index[N]; //
}
parameters {
// real year_b0; // year effect base
real b0; // intercept
real b1; // coefficient for the floor size
real b2; // coefficient for the age
real b3; // coefficient for the structure dummy (is_rc_pc)
vector[K] year_b;
vector[J] chiiki_b;
real<lower=0> sigma; // sigma for the total rent
real<lower=0> sigma_year; // sigma for the year effect
real<lower=0> sigma_chiiki; // sigma for the chiiki effect
}
model {
// prior
sigma ~ cauchy(0, 100);
sigma_year ~ cauchy(0, 100);
sigma_chiiki ~ cauchy(0, 100);
b0 ~ normal(0, 10000);
year_b ~ normal(0, sigma_year);
chiiki_b ~ normal(0, sigma_chiiki);
price ~ normal(b0 + year_b[year_effect_index] +
chiiki_b[chiiki_effect_index] +
b1*floor_size + b2*age + b3*is_rc_pc, sigma);
}これをRでモデルフィットするには、以下を実行。
library(rstan)
# 高速化のためのおまじない
options(mc.cores = parallel::detectCores())
rstan_options(auto_write = TRUE)
options(buildtools.check = function(action) TRUE )
setagaya_cleaned <- setagaya_cleaned %>%
mutate(year_dummy = transaction_year - 2004, # the smallest year, 2005 will be 1, and 2006 is 2, etc.
chiiki_dummy = as.factor(chikumei),
chiiki_dummy_index = as.integer(chiiki_dummy) ) # add index number for each area
data_list_chiiki <- list(
N = dim(setagaya_cleaned)[1],
K = max(setagaya_cleaned$year_dummy),
J = max(setagaya_cleaned$chiiki_dummy_index),
year_effect_index = setagaya_cleaned$year_dummy,
chiiki_effect_index = setagaya_cleaned$chiiki_dummy_index,
price = setagaya_cleaned$price,
age = setagaya_cleaned$age,
floor_size = setagaya_cleaned$nobeyuka_menseki_number,
is_rc_pc = setagaya_cleaned$is_rc_pc
)
setagaya_chiiki_fit <- stan(
file = "./stan/chiiki_lm.stan",
data = data_list_chiiki,
seed = 1
)取引年次で2004を引いているのは、データが2005年以降存在して、2005年を1、2006年を2、というように番号をつけるのが目的。パソコンのスペックにもよるが、モデルフィットに数分かかる。
モデル結果を見てみる。
setagaya_chiiki_fit_mcmc_sample <- rstan::extract(setagaya_chiiki_fit, permuted = FALSE)
print(setagaya_chiiki_fit, probs=c(0.05, 0.95))
mean se_mean sd 5% 95% n_eff Rhat
b0 3737.24 17.68 471.49 2939.73 4489.37 711 1.01
b1 49.65 0.01 0.58 48.70 50.60 4145 1.00
b2 -114.72 0.09 6.25 -124.86 -104.82 5080 1.00
b3 1437.99 4.68 295.97 952.51 1918.76 4006 1.00ベイズモデルだと便利なのは、信頼区間をそのまま確率として理解できる点。例えばこのモデルだと築年の係数（b2）が-125 から -104となっているが、これは築年の係数がこの間になるのが90％の確率、というように解釈できる。この区間が0を含んでいるようだと、係数が０なんじゃないか、つまりマンション価格には影響を与えていないんではないか、と解釈できる。
ベイズモデルにすると単純な線形モデルとは結果が若干だけ異なるが、ほとんど前回モデルと同じような結果になっている。
前回モデル：価格＝3315 + 50 x 延べ床面積 - 155 x 築年 + 1748 x 構造 ＋ 取引年次効果 ＋ 地域効果
ベイズモデル：価格＝3337 + 50 x 延べ床面積 - 115 x 築年 + 1438 x 構造 ＋ 取引年次効果 ＋ 地域効果
築年の係数を地域ごとに変える
前回記事の最後で、平米単価と築年の関係をグラフ化して、どうも築年の関係が地域ごとに異なるように見える、という点について触れた。その時のグラフを再掲。
前回のモデルは、築年はどの地域にあろうが一律1年につき115万円減らす、というものだった。これを、例えば玉川なら100万、東玉川なら120万、というように変えるというモデルを考える。
これはStanで以下のように書ける。
data {
int<lower=0> N;
int<lower=0> K; // year effect index
int<lower=0> J; // chiiki effect index
vector[N] price;
vector[N] floor_size;
vector[N] age;
vector[N] is_rc_pc;
int<lower=1, upper=K> year_effect_index[N]; // 1 = 2005, 2 = 2006, 3 = 2007, etc.
int<lower=1, upper=J> chiiki_effect_index[N]; //
}
parameters {
// real year_b0; // year effect base
real b0;
real b1;
// real b2;
real b3;
vector[K] year_b;
vector[J] chiiki_b;
vector[J] chiiki_age_b; // the slope for age by chiiki
real<lower=0> sigma; // sigma for the total rent
real<lower=0> sigma_year; // sigma for the year effect
real<lower=0> sigma_chiiki; // sigma for the chiiki effect
real<lower=0> sigma_chiiki_age; // sigma for the slope of age by area
}
model {
// prior
sigma ~ cauchy(0, 100);
sigma_year ~ cauchy(0, 100);
sigma_chiiki ~ cauchy(0, 100);
sigma_chiiki_age ~ cauchy(0, 100);
b0 ~ normal(0, 10000);
b1 ~ normal(0, 1000);
b3 ~ normal(0, 1000);
year_b ~ normal(0, sigma_year);
chiiki_b ~ normal(0, sigma_chiiki);
chiiki_age_b ~ normal(0, sigma_chiiki_age);
price ~ cauchy(b0 + year_b[year_effect_index] +
chiiki_b[chiiki_effect_index] +
b1*floor_size + chiiki_age_b[chiiki_effect_index] .* age + b3*is_rc_pc, sigma);
}変わったのは、chiiki_age_bという地域ごとの築年の係数を定義したところ、その係数を生み出すための標準偏差sigma_chiiki_ageを定義したところ、さらにモデルにchiiki_age_bを追加した点。
// model formula
price ~ cauchy(b0 + year_b[year_effect_index] +
chiiki_b[chiiki_effect_index] +
b1*floor_size + chiiki_age_b[chiiki_effect_index] .* age + b3*is_rc_pc, sigma);
あと、前回のデータでかなり高い物件が散見されたので、正規分布ではなくてコーシー分布でモデル化するように変更した。この式の書き方は「StanとRでベイズ統計モデリング」を参考にした。このモデルをRでフィットする。
data_list_chiiki <- list(
N = dim(setagaya_cleaned)[1],
K = max(setagaya_cleaned$year_dummy),
J = max(setagaya_cleaned$chiiki_dummy_index),
year_effect_index = setagaya_cleaned$year_dummy,
chiiki_effect_index = setagaya_cleaned$chiiki_dummy_index,
price = setagaya_cleaned$price,
age = setagaya_cleaned$age,
floor_size = setagaya_cleaned$nobeyuka_menseki_number,
is_rc_pc = setagaya_cleaned$is_rc_pc
)
setagaya_chiiki_age_fit <- stan(
file = "./stan/chiiki_random_age.stan",
data = data_list_chiiki,
seed = 1,
iter = 2000,
control = list(max_treedepth = 15)
)これは５、6分かかった。結果を見ると、
setagaya_chiiki_age_fit_mcmc_sample <- rstan::extract(setagaya_chiiki_age_fit, permuted = FALSE)
print(setagaya_chiiki_age_fit, probs=c(0.05, 0.95))
mean se_mean sd 5% 95% n_eff Rhat
b0 2584.51 4.21 152.95 2338.86 2839.98 1322 1.00
b1 47.37 0.01 0.65 46.28 48.33 3476 1.00
b3 -294.50 1.06 86.92 -435.80 -148.34 6740 1.00
(一部省略)
chiiki_b[1] 451.63 3.57 261.32 16.31 879.16 5353 1.00
chiiki_b[2] 711.47 3.65 281.30 241.48 1153.54 5948 1.00
chiiki_b[3] -736.95 3.13 195.15 -1057.90 -422.18 3892 1.00
chiiki_b[4] 212.54 4.36 334.38 -336.80 771.47 5872 1.00
chiiki_b[5] -1049.01 3.19 207.32 -1391.58 -709.06 4222 1.00
(一部省略)
chiiki_age_b[1] -66.20 0.14 12.90 -86.64 -44.87 8112 1.00
chiiki_age_b[2] -81.95 0.15 13.21 -102.12 -58.61 7752 1.00
chiiki_age_b[3] -60.86 0.09 8.32 -74.82 -47.18 9517 1.00
chiiki_age_b[4] -57.85 0.13 12.15 -78.65 -38.81 8794 1.00
chiiki_age_b[5] -70.24 0.13 12.50 -89.18 -48.00 9029 1.00
chiiki_age_b[6] -37.37 0.29 16.10 -68.98 -17.40 3102 1.00
chiiki_age_b[7] -62.53 0.08 7.50 -74.83 -50.37 8756 1.00このモデルは、築年が古くなるほど安くなる、と言う関係が世田谷区で一律-108万とするのではなくて、例えばこの地域は高級住宅地なので築年が古くなろうが値下がりが少ない、と言うようなモデルの作りにしている。
まず、ベースは2585万（係数b0）。次に、延べ床面積が増えるごとに、47万円高くなる（係数b1）。200平米なら9400万増えることという意味。RCは295万安くなる、というよくわからない結果（係数b3）。まあそういうものなのかもしれない。築年の係数はchiiki_age_bで表されていて、地域を１から順番に番号をつけている（59まである）。築年の影響は大体どこもマイナスだが、その影響が小さいところ、大きいところがある。地域の影響は、chiiki_bという係数で表されている。前回のモデルよりも全体的に小さくなっている。これは築年でも地域差を考慮に入れたことから、地域差のうち築年で説明できる部分は、そちらの係数に加わったのが理由と考えられる。とは言え、係数の大小は前回の結果に似てはいる。例えばマイナスが大きい地域は
chiiki_age_model_estimate <- as.data.frame(summary(setagaya_chiiki_age_fit)$summary[,c("mean","2.5%","97.5%")])
chiiki_age_model_estimate <- chiiki_age_model_estimate %>%
mutate(param_name = rownames(chiiki_age_model_estimate))
rownames(chiiki_age_model_estimate) <- NULL
colnames(chiiki_age_model_estimate) <- c("mean","p025","p975","param_name")
chiiki_age_model_estimate %>%
filter(grepl("chiiki_b", param_name )) %>%
mutate(chikumei = levels(setagaya_chikumei_factor)) %>%
arrange(mean) %>%
dplyr::select(chikumei,mean,p025,p975,param_name) %>%
filter(mean<0 & p975<0) %>%
head(5)
#chikumei mean p025 p975
#1 八幡山 -1190.416 -1751.713 -606.98445
#2 喜多見 -1153.425 -2227.353 -22.35367
#3 北烏山 -1138.565 -1601.493 -699.05486
#4 給田 -1115.215 -1759.956 -398.05418
#5 上祖師谷 -1049.011 -1456.259 -645.20624など、前回値段が安いと出た地域（喜多見、宇奈根、玉川田園調布、北烏山、八幡山）と被っている地域が多い。逆に値段が高い地域は、
chiiki_age_model_estimate %>%
filter(grepl("chiiki_b", param_name )) %>%
mutate(chikumei = levels(setagaya_chikumei_factor)) %>%
arrange(mean) %>%
dplyr::select(chikumei,mean,p025,p975,param_name) %>%
filter(mean>0 & p025>0) %>%
tail(5)
chikumei mean p025 p975 param_name
10 奥沢 768.8007 325.1533 1219.416 chiiki_b[22]
11 下馬 815.7732 464.0872 1161.707 chiiki_b[8]
12 北沢 859.7191 473.2526 1260.956 chiiki_b[14]
13 代沢 1160.1169 465.9408 1783.669 chiiki_b[11]
14 等々力 1193.0586 710.3462 1714.448 chiiki_b[47]高い地域は少し結果が異なる。玉川が入っていない。ひょっとしたら、玉川は再開発があったので新しい物件が多かったのかもしれない。だから築年を考慮したら、相場が高いというわけではないのかもしれない。
これらの推定値の幅をグラフ化する。
p2 <- mcmc_intervals(setagaya_chiiki_age_fit_mcmc_sample,
pars = c("chiiki_b[13]", "chiiki_b[18]", "chiiki_b[15]",
"chiiki_b[50]", "chiiki_b[5]"))
p2 <- p2 + scale_y_discrete(
labels = c("chiiki_b[13]" = levels(setagaya_chikumei_factor)[13],
"chiiki_b[18]" = levels(setagaya_chikumei_factor)[18],
"chiiki_b[15]" = levels(setagaya_chikumei_factor)[15],
"chiiki_b[50]" = levels(setagaya_chikumei_factor)[50],
"chiiki_b[5]" = levels(setagaya_chikumei_factor)[5]))
print(p2)
p3 <- mcmc_intervals(setagaya_chiiki_age_fit_mcmc_sample,
pars = c("chiiki_b[22]", "chiiki_b[8]", "chiiki_b[14]",
"chiiki_b[11]", "chiiki_b[47]"))
p3 <- p3 + scale_y_discrete(
labels = c("chiiki_b[22]" = levels(setagaya_chikumei_factor)[22],
"chiiki_b[8]" = levels(setagaya_chikumei_factor)[8],
"chiiki_b[14]" = levels(setagaya_chikumei_factor)[14],
"chiiki_b[11]" = levels(setagaya_chikumei_factor)[11],
"chiiki_b[47]" = levels(setagaya_chikumei_factor)[47]))
print(p3)
地域ごとの築年の係数を分析
では新たにモデルに追加した、地域ごとの築年の係数を見てみる。築年の係数が大きい地域と（マイナスが大きい）、小さい地域（マイナスがすくない）を見てみる。
chiiki_age_model_estimate %>%
filter(grepl("chiiki_age_b", param_name )) %>%
mutate(chikumei = levels(setagaya_chikumei_factor)) %>%
arrange(mean) %>%
# ０以下のものを取り出す
dplyr::select(chikumei,mean,p025,p975,param_name) %>%
filter(mean<0 & p975<0) %>%
head(5)
#chikumei mean p025 p975 param_name
#1 宇奈根 -108.06893 -139.7200 -74.80530 chiiki_age_b[23]
#2 大蔵 -104.88362 -158.5871 -45.62960 chiiki_age_b[20]
#3 梅丘 -96.20876 -117.1069 -75.16085 chiiki_age_b[36]
#4 祖師谷 -90.71967 -116.5335 -58.91660 chiiki_age_b[46]
#5 喜多見 -88.52433 -130.3374 -47.76930 chiiki_age_b[18]
これらの地域が築年が1年古くなると価格が下がりやすい地域といえる。例えば宇奈根は1年古くなるごとに108万安くなるので、築10年だと1080万円安くなる地域といえる。推定値の範囲はこちら。
逆に築年の係数が小さい、築年が古くなってもそれほど安くならない地域のトップ５はこちら。
chiiki_age_model_estimate %>%
filter(grepl("chiiki_age_b", param_name )) %>%
mutate(chikumei = levels(setagaya_chikumei_factor)) %>%
arrange(mean) %>%
# ０以下のものを取り出す
dplyr::select(chikumei,mean,p025,p975,param_name) %>%
# filter(mean>0) %>%
tail(5)
#chikumei mean p025 p975 param_name
#55 上野毛 -37.37407 -77.53828 -14.08838 chiiki_age_b[6]
#56 尾山台 -35.79353 -85.47200 22.01942 chiiki_age_b[25]
#57 玉川台 -14.73926 -51.76928 11.33140 chiiki_age_b[42]
#58 桜新町 10.07616 -40.98163 50.54783 chiiki_age_b[35]
#59 玉川 15.09538 -75.34670 93.84717 chiiki_age_b[41]
築年の影響が少ないエリアは上記。玉川は推定値が15万とプラスになっているが、推定値の区間を見ると-75から93と0を含んでいるので、築年の影響がほぼないといえる。ただし、先ほども書いたように、単に築年が新しい物件が多いだけなのかもしれない。上野毛は-37万と、宇奈根に比べると築年の影響がかなり小さい。10年古くなっても370万円しか安くならない。推定値の幅は下図。
p5 <- mcmc_intervals(setagaya_chiiki_age_fit_mcmc_sample,
pars = c("chiiki_age_b[6]", "chiiki_age_b[25]", "chiiki_age_b[42]",
"chiiki_age_b[35]", "chiiki_age_b[41]"))
p5 <- p5 + scale_y_discrete(
labels = c("chiiki_age_b[6]" = levels(setagaya_chikumei_factor)[6],
"chiiki_age_b[25]" = levels(setagaya_chikumei_factor)[25],
"chiiki_age_b[42]" = levels(setagaya_chikumei_factor)[42],
"chiiki_age_b[35]" = levels(setagaya_chikumei_factor)[35],
"chiiki_age_b[41]" = levels(setagaya_chikumei_factor)[41]))
print(p5)
これらの５つの地域は、他の地域よりも築年による値段の下がりが小さい地域といえる。特に玉川、桜新町、玉川台、尾山台は推定区間が0を跨いでいるので、築年の影響がほぼない古くなっても新しくても値段が変わらない、という地域と言える。
モデルから見た世田谷区の不動産事情
ここまでの結果をまとめると、こんなことが言えそう。
地域による差が結構大きい
前回のモデルだと、地域によって相場に違いがあって、全体的に価格が高い地域（二子玉川のあたり）と安い地域（狛江市に近いあたり）があるというのがわかった。今回のモデルでさらにわかったこととして、築年の効果も地域によって違いがあるという点。例えば玉川地域は再開発があったのが大きいのだろうが、そんなに相場に変化が見られなかった。
古くなっても値段がそれほど下がらない地域がある
築年の係数も地域によって大きく異なるようだった。例えば玉川地域は＋15万円となっていて、古くなっても値段が下がるどころか値上がりするという関係になっていた。もちろんこれは本当に値段が上がるわけではなくて、再開発などその地域の相場の変化があっただけなのかもしれない。実際、玉川の築年ごとの取引物件の数を調べてみると、築20年内の物件が他の地域と比べると多いことがわかる。さらにそれより新しい物件の取引データがほとんどないのも特徴。
hist_age_high_price5 <- setagaya_cleaned %>% filter(chikumei %in% c("玉川","桜新町","玉川台","尾山台","上野毛")) %>%
ggplot(aes(x=age)) + geom_histogram(binwidth = 5) +
facet_wrap(~chikumei)
plot(hist_age_high_price5)
それ以外でも上野毛や尾山台など、築年による価格の変化が36万や37万円と推測される地域もあれば、宇奈根、大倉のように108万、105万と３倍近い差があった。このように、地域による全体的な価格相場の違いがあるだけでなく、築年の影響についても地域差があるとわかった。
まとめ
国土交通省のデータを使ってモデルを作ってみた。案外これだけでも面白いことが分かる。地域ごとに相場が異なるというのはモデルでも明確にわかり、ここが安め、高め、築年の影響あまり受けない、受ける、ということがこのデータだけでも十分にわかる。単なるデータ集計ではここまではわからないので、モデルの力と言えそう。
２３区の別の地域とか、地方都市のモデルも見てみると面白そう。特に取引年次による相場の変化が、東京よりも地方都市のが遅れてくるのではないかとか、築年の影響がもっと大きい・小さい、ということがあるのかもしれない。
参考文献
StanとRでベイズ統計モデリング
モデリングの基本的なところはこれを参照。Stan記法のベクトル化、階層モデルの実装も参考になる。
実践Data Scienceシリーズ RとStanではじめる ベイズ統計モデリングによるデータ分析入門
mcmc_intervalsを使ったグラフの書き方を参考にした。
Statistical Rethinking: A Bayesian Course with Examples in R and Stanm
モデルの作り方について参考にした。今調べたら第二版が出ている。