銀行員 RとPythonに出会う

Rネタを中心に、いろいろと更新していきます

~続・不動産とファイナンス・賃貸物件入居者編(3)~「機械学習を使って東京23区のお買い得賃貸物件を探してみた」を千葉県で再度やってみる(予測編)

今回は、不動産賃貸物件の賃料予測シリーズの最終回です。

(分析のソースコードは最後に一括掲載します)

 

データクレンジングした千葉県の賃貸物件情報データセットを使って、予測に入っていきます。

 

その前に、データクレンジングについて補足があります。

スクレイピングの段階で物件の住所が抜け落ちていることに気づきましたので、後から付け足しました。

市ごとにデータを抽出していたため追加は問題なかったのですが、コードを参照される方はご注意ください。

 

毎度のことながら、こちらが参考記事です。

www.analyze-world.com

 

基礎分析

実は今回の記事は下の記事の続編でして、全体の外観は既に見ていました。

前回と今回の違いは、物件の構造や駐車場有無など新たに変数を追加した点です。

スクレイピングにすごく時間がかかりました)

なお、今回も市別物件数等の基礎統計は出していますが、大きな違いはないため重複する内容は前回に譲ります。

https://blog.hatena.ne.jp/d_s/d-s.hatenablog.com/edit?entry=10257846132614731711

 

では、最初に今回使うデータセットを載せておきます。

df(38906,19)

> str(df)

'data.frame':	38906 obs. of  19 variables:

 $ name        : Factor w/ 26511 levels "#NAME?","(仮)D-room千葉寺",..: 17354 2060 24672 6816 24832 12633 12634 23640 3630 3631 ...

 $ city        : Factor w/ 10 levels "浦安市","鎌ヶ谷市",..: 6 6 6 6 6 6 6 6 6 6 ...

 $ shikikin    : int  38000 0 40000 0 0 38000 39000 0 0 0 ...

 $ reikin      : int  38000 0 0 0 0 38000 39000 45000 0 0 ...

 $ hoshoukin   : int  0 0 0 0 0 0 0 0 0 0 ...

 $ layout      : Factor w/ 8 levels "1DK","1K","1LDK",..: 2 1 1 3 8 2 2 8 3 3 ...

 $ area        : num  21.1 34.3 33.1 42 34.3 ...

 $ direction   : Factor w/ 9 levels "-","西","東",..: 3 3 8 5 3 6 6 4 5 2 ...

 $ type        : Factor w/ 5 levels "アパート","その他",..: 4 4 5 4 4 1 1 1 4 4 ...

 $ age         : int  89 52 48 46 52 22 22 44 48 48 ...

 $ rout1       : Factor w/ 21 levels "JR外房線","JR京葉線",..: 4 4 4 15 4 4 4 4 7 1 ...

 $ station     : Factor w/ 159 levels "おゆみ野","くぬぎ山",..: 8 8 89 50 8 8 8 8 100 100 ...

 $ distance    : int  72 25 100 25 28 80 52 28 68 68 ...

 $ construction: Factor w/ 10 levels "その他","プレコン",..: 7 6 10 6 6 10 10 7 7 7 ...

 $ floor       : Factor w/ 25 levels "1","1-2","1-3",..: 16 12 1 20 12 1 12 12 12 16 ...

 $ height      : int  4 4 1 5 4 2 2 2 4 4 ...

 $ car.dum     : Factor w/ 4 levels "近隣","付無料",..: 3 1 2 4 4 3 4 4 4 4 ...

 $ free_rent   : int  0 1 0 0 1 0 0 0 0 0 ...

 $ monthly_cost: int  38000 53000 40000 27000 50000 40000 41000 45000 49000 54000 ...

 

 

 

 

はじめに、家賃と相関が高い変数を見ておきます。

家賃との相関係数が絶対値で0.5超のものを抽出してあります。 

f:id:d_s:20180906224319p:plain

 

monthly_cost:家賃、area:専有面積、age:築年ですが、家賃と専有面積の相関が1番高そうです。

感覚的にもそうだろうなと思いますよね。予測の際に重要になりそうな変数です。

 

次は、建物構造をスクレイピングしてきているので家賃との対応を見ておきます。

f:id:d_s:20180906225837p:plain

やはり木造が・・・という結果になりました。

自分もアパート暮らしを経験したことがありますが、物件の構造は気にしていました。

鉄筋だの鉄骨だの細かくこだわっていませんでしたが、木造だけは「音が気になる」という先入観を持っていたので、当時選択肢にすら入れていなかったのを覚えています。

 

以下、変数ごとにヒストグラムを描いて分布を確認して、データが寄っていたら対数変換して正規分布に近くなるか確認して・・・と続いていくわけなのですが、長々としてしまいますのでこんな例を載せておきます。

f:id:d_s:20180906231428p:plain

f:id:d_s:20180906231713p:plain

最寄り駅までの距離ですが、これを対数変換すると

f:id:d_s:20180906231547p:plain

f:id:d_s:20180906231740p:plain

幾分か正規分布、線形に近づいたと思います。

今回はこのまま重回帰に進みますが、もっとぐりぐりと探索したいところではあります。その他処理については、コードに解説を。

 (後で気づきましたが、元のデータセットで対数変換しようと思ったものを変換し忘れていました。orz)

 

次に予測に入っていきます。

今回もデータセットを学習:テストで8:2に分けて検証します。

学習データ:tarin(31111, 14)

テストデータ:test(7776, 14)

 予測1 重回帰

まずは重回帰分析です。

コードを参照していただければと思いますが、カテゴリー変数はダミー化しています。

先に述べたように、今回Step関数(AIC最小化)で変数選択していますので、下の結果は変数選択後の結果です。

> summary(lm.step)



Call:

lm(formula = monthly_cost ~ city.浦安市 + city.鎌ヶ谷市 + city.市川市 + 

    city.松戸市 + city.千葉市 + city.船橋市 + city.八千代市 + 

    layout.1DK + layout.1K + layout.1LDK + layout.1SDK + layout.1SK + 

    layout.1SLDK + area + type.アパート + type.その他 + type.テラスハウス + 

    type.マンション + age + rout1.JR外房線 + rout1.JR京葉線 + 

    rout1.JR常磐線 + rout1.JR総武線 + rout1.JR総武線快速 + 

    rout1.JR総武本線 + rout1.JR内房線 + rout1.JR武蔵野線 + 

    rout1.つくばエクスプレス + rout1.京成千原線 + rout1.京成千葉線 + 

    rout1.京成本線 + rout1.新京成線 + rout1.千葉都市モノレール + 

    rout1.都営新宿線 + rout1.東京メトロ東西線 + rout1.東武野田線 + 

    rout1.北総線 + distance + construction.軽量鉄骨 + construction.鉄筋コン + 

    construction.鉄骨 + construction.鉄骨鉄筋 + floor.1 + floor.12 + 

    floor.13 + floor.14 + floor.15 + floor.17 + floor.18 + floor.19 + 

    floor.2 + `floor.2-3` + floor.3 + floor.36 + floor.38 + floor.6 + 

    height + car.dum.近隣 + car.dum.付無料 + car.dum.敷地内, 

    data = train_dummy)



Residuals:

   Min     1Q Median     3Q    Max 

-47975  -4407   -385   3754  92307 



Coefficients:

                           Estimate Std. Error  t value             Pr(>|t|)    

(Intercept)               51597.992   1333.095   38.705 < 0.0000000000000002 ***

city.浦安市                9991.120    309.562   32.275 < 0.0000000000000002 ***

city.鎌ヶ谷市              -860.327    420.582   -2.046             0.040808 *  

city.市川市                7080.343    218.837   32.354 < 0.0000000000000002 ***

city.松戸市                 785.224    171.400    4.581   0.0000046401196611 ***

city.千葉市               -4705.804    206.609  -22.776 < 0.0000000000000002 ***

city.船橋市                2900.313    179.287   16.177 < 0.0000000000000002 ***

city.八千代市             -4384.625    295.029  -14.862 < 0.0000000000000002 ***

layout.1DK                 2178.404    202.462   10.760 < 0.0000000000000002 ***

layout.1K                  1092.225    124.968    8.740 < 0.0000000000000002 ***

layout.1LDK                3864.571    212.770   18.163 < 0.0000000000000002 ***

layout.1SDK                3438.787   1274.857    2.697             0.006992 ** 

layout.1SK                 3032.182    976.678    3.105             0.001907 ** 

layout.1SLDK               7881.395    932.678    8.450 < 0.0000000000000002 ***

area                        933.801      8.357  111.738 < 0.0000000000000002 ***

type.アパート            -16857.882   1252.718  -13.457 < 0.0000000000000002 ***

type.その他              -15205.806   2156.760   -7.050   0.0000000000018226 ***

type.テラスハウス        -14009.415   1547.312   -9.054 < 0.0000000000000002 ***

type.マンション          -15053.623   1266.946  -11.882 < 0.0000000000000002 ***

age                        -615.966      3.725 -165.344 < 0.0000000000000002 ***

rout1.JR外房線           3185.541    416.257    7.653   0.0000000000000202 ***

rout1.JR京葉線          10369.839    373.456   27.767 < 0.0000000000000002 ***

rout1.JR常磐線           3613.268    327.765   11.024 < 0.0000000000000002 ***

rout1.JR総武線          10559.045    303.565   34.783 < 0.0000000000000002 ***

rout1.JR総武線快速      12349.595    389.127   31.737 < 0.0000000000000002 ***

rout1.JR総武本線         2588.423    414.103    6.251   0.0000000004139921 ***

rout1.JR内房線           5378.499    461.783   11.647 < 0.0000000000000002 ***

rout1.JR武蔵野線         1811.428    408.226    4.437   0.0000091398650456 ***

rout1.つくばエクスプレス   3203.856    365.079    8.776 < 0.0000000000000002 ***

rout1.京成千原線           1333.445    440.695    3.026             0.002482 ** 

rout1.京成千葉線           8878.187    445.089   19.947 < 0.0000000000000002 ***

rout1.京成本線             2570.059    295.439    8.699 < 0.0000000000000002 ***

rout1.新京成線            -1117.875    336.274   -3.324             0.000887 ***

rout1.千葉都市モノレール   1237.903    420.995    2.940             0.003280 ** 

rout1.都営新宿線          11057.461    793.120   13.942 < 0.0000000000000002 ***

rout1.東京メトロ東西線     7552.354    352.766   21.409 < 0.0000000000000002 ***

rout1.東武野田線          -1299.222    355.462   -3.655             0.000258 ***

rout1.北総線               1446.186    497.302    2.908             0.003639 ** 

distance                   -279.290      4.316  -64.716 < 0.0000000000000002 ***

construction.軽量鉄骨      2811.936    119.929   23.447 < 0.0000000000000002 ***

construction.鉄筋コン      2424.085    229.301   10.572 < 0.0000000000000002 ***

construction.鉄骨          1085.687    197.038    5.510   0.0000000361614633 ***

construction.鉄骨鉄筋      -779.408    371.952   -2.095             0.036139 *  

floor.1                   -3087.492    175.923  -17.550 < 0.0000000000000002 ***

floor.12                   2906.244   1011.471    2.873             0.004065 ** 

floor.13                   2868.653   1216.872    2.357             0.018410 *  

floor.14                   7940.393   1781.355    4.458   0.0000083210633342 ***

floor.15                  30143.089   4101.146    7.350   0.0000000000002032 ***

floor.17                  11181.280   5029.571    2.223             0.026216 *  

floor.18                  33601.725   7074.130    4.750   0.0000020438072318 ***

floor.19                  34882.878   7073.643    4.931   0.0000008206846043 ***

floor.2                   -1399.758    172.474   -8.116   0.0000000000000005 ***

`floor.2-3`                8776.983   4071.860    2.156             0.031129 *  

floor.3                    -701.876    177.975   -3.944   0.0000804200069929 ***

floor.36                  45586.183   4128.933   11.041 < 0.0000000000000002 ***

floor.38                  12905.202   7082.366    1.822             0.068440 .  

floor.6                    -744.190    365.900   -2.034             0.041974 *  

height                     1311.134     24.931   52.591 < 0.0000000000000002 ***

car.dum.近隣               -868.668    116.089   -7.483   0.0000000000000747 ***

car.dum.付無料             4812.216    625.234    7.697   0.0000000000000144 ***

car.dum.敷地内             -179.877    101.314   -1.775             0.075834 .  

---

Signif. codes:  0***0.001**0.01*0.05 ‘.’ 0.1 ‘ ’ 1



Residual standard error: 7036 on 31050 degrees of freedom

Multiple R-squared:  0.843,	Adjusted R-squared:  0.8427 

F-statistic:  2778 on 60 and 31050 DF,  p-value: < 0.00000000000000022

Intercept(切片)のEstimate(偏回帰係数)が51597ですので、解釈としては51597円をベースとして、加算減算し家賃が求まる式ができました。

例えば市で見ると、浦安の偏回帰係数が9991、千葉が-4705ですので、その物件が浦安であれば51597円+9991円、千葉であれば51597円-4705円といった具合で計算されていきます。

市を見ると、東京から遠ざかるにつれ係数が下がっていくのでこれは納得ですね。

area(専有面積)については、偏回帰係数が933と出ています。

1平米あたり933円という解釈になりますが、今回のデータセットの平米家賃単価(平均家賃÷平均専有面積)は2193円ですので、その差額が他の変数で説明されている部分と捉えられそうです。

今回の重回帰分析の決定係数は0.84なのでなかなか良い結果と言えるのではないでしょうか。

試しに、テストデータで予測を確認しておきます。

f:id:d_s:20180906234229p:plain

重回帰でそれなりに予測できていそうです。

こうなってくると、次のランダムフォレストでの結果も楽しみです。

 

予測2 ランダムフォレスト

比較しやすいように、先に結果を載せてしまいます。

f:id:d_s:20180906234927p:plain

さらに精度がよくなったようです。

実は今回、ランダムフォレストの学習時間短縮のために、学習データをさらに6割に削減したので不利なはずでした。

それなのにこの結果とは。。。改めてランダムフォレストの強力さを実感しました。

(テストデータは重回帰とまったく一緒です)

 

なお、今回のパラメータは例のごとくtuneRF関数でmtry(特徴量の数)を求めました。

mtryは8が最適のようです。

f:id:d_s:20180906235127p:plain

変数重要度は以下の通りアウトプットされました。

ランダムフォレストでは築年数、専有面積が重要な要因のようです。

f:id:d_s:20180906235042p:plain

以上、重回帰とランダムフォレストをざっとみてきましたが、ここで精度を比較しておきます。

RMSE(平均)という「平均化された誤差」という指標で比較したいと思うのですが、

> ##### RMSE #####

> sqrt( sum((test$monthly_cost - lm_pred)^2) / length(test$monthly_cost) )

[1] 7095.295

> sqrt( sum((test$monthly_cost - rf_pred)^2) / length(test$monthly_cost) )

[1] 5133.426

 上記の通り、重回帰:7095、ランダムフォレスト5133と、ランダムフォレストのほうがRMSEが小さい(誤差が相対的に低い)ので、ランダムフォレストのほうが良い結果を出すことができました。

 

予測3 お得物件を見つける

やっとここまでたどり着きました。

以下はランダムフォレストで予測した家賃をもとに、実際の家賃との差額からお得係数を算出してソートしたものです。

(お得係数:僕が勝手に考えた係数で、差分÷家賃共益費で算出)

 

物件 間取 面積 築年 最寄駅距離 構造 家賃+共益費 予測家賃 差分 お得係数
景中寮 流山市 ワンルーム 18.00 31 東武野田線 9 木造 12000 31058.24 19058.24 1.5881867
カルフォルニアハウス江戸川台202号室 流山市 1K 7.30 14 東武野田線 6 鉄筋コン 20000 47337.43 27337.43 1.3668713
カルフォルニアハウス江戸川台201号室 流山市 1K 6.60 14 東武野田線 6 鉄筋コン 20000 46878.56 26878.56 1.3439278
カルフォルニアハウス江戸川台101号室 流山市 1K 7.70 14 東武野田線 6 鉄筋コン 20000 46210.84 26210.84 1.3105422
JR総武線津田沼駅2階建築26年 船橋市 1K 20.00 26 JR総武線 15 木造 22000 47256.80 25256.80 1.1480366
ヴィレッジダイドー 千葉市 1LDK 42.00 46 千葉都市モノレール 25 鉄筋コン 27000 57002.69 30002.69 1.1112108
ビレッジハウス二和1号棟305号室 船橋市 1K 45.36 57 新京成線 4 鉄筋コン 34000 66637.96 32637.96 0.9599399
アーバンハイム金丸 習志野市 ワンルーム 18.00 19 京成本線 6 木造 27000 48302.86 21302.86 0.7889947
大和田ハイツI-1 八千代市 1K 20.00 50 京成本線 10 木造 19000 33909.16 14909.16 0.7846928
ハイホーム田中305号室 市川市 1K 19.80 48 JR総武線快速 7 鉄筋コン 31000 55022.72 24022.72 0.7749265

 

この中から、気になる物件を見てみます。

まず、お得度第一位の物件   

f:id:d_s:20180906221406p:plain

f:id:d_s:20180906221506p:plain

 おお、男子寮がランクインしてきました。

学生向けの物件でしょうか。ちょっと前提から外れる物件な気もしなくはないですが、破格の値段です。

寮なので、シェアハウスのごとく部屋以外は共同のようですが、写真で見る限り中は綺麗でした。

 

つぎは、個人的に目に留まった物件で、総武線市川駅徒歩7分です。

お得ランキングとしては10位です。築古なので年季が入っているようですが・・・  

f:id:d_s:20180906221327p:plain

f:id:d_s:20180906221444p:plain

家賃共益費込みで31,000円!!!!!!!

 

築古ワンルームですが、東京の近さと、最寄り駅の近さを勘案するとすごくお得な物件ではないでしょうか?おまけに角部屋、鉄筋コンクリートです。

 

ちなみに、中の写真を見てみると・・・

 

f:id:d_s:20180906221616p:plain

f:id:d_s:20180906221635p:plain

 

 築古だけとリフォームしてある感じで綺麗!!!

 

そしてお風呂が見当たらない!!!!!!!!!

 

なるほど、そういうことだったんですね笑

 

もしかしたら、お得度上位物件は、変数で考慮されていない特徴ある訳あり物件が並んでいるのかもしれません。笑

 

逆に言えば、一般的な物件は予測精度が良いのかもしれません。

風呂有り無し等、スクレイピングレベルで再度試すのは面倒ですが、今回の施行は面白いものとなりました。

 

次に

賃貸物件の賃料予測は無事行うことができましたが、次回は収益物件の値段の予測問題に入っていきます。収益物件の価格の妥当性を検証したいわけです。

実は、収益物件売買情報掲載サイトの「楽待」をスクレイピングして、既にデータセットを作成してあります。

楽待でも、今回の賃料予測に対応して千葉県の収益物件を取ってきていますので、収支予測を今回の結果をベースに行い、DCFベースの物件価格評価を目指します。

DCFをやるからには、割引率は収益物件の想定利回り(投資家が求める期待リターンとも言えるのか)の構造を重回帰等で割り出して算出してみたら面白いんじゃないかと考えています。

 

賃貸物件の賃料予測は、これにて終わりです。

 

本日のコード

library(ggplot2)

library(tidyverse)

library(caret)

library(randomForest)

library(ggrepel)

library(psych)

library(kernlab)

library(knitr)

library(corrplot)



options(scipen=100)



# データセットの読み込み

df <- read.csv(データセットを置いてあるファイルパスを入力,

               fileEncoding = "cp932", sep = ",")



# 後の工程を考えて全角中黒は変換

df$type <- df$type %>% 

  str_replace('テラス・タウンハウス','テラスハウス') %>%

  as.factor()



# データチェック

df %>%

  str()



# rent1とrent2を足してmonthly_costを作成、rent1・rent2は削除(rent1:家賃、rent2:共益費)

df <- df %>%

  mutate(monthly_cost=rent1+rent2) %>%

  select(-c(3,4))



# 相関チェック

numericVars <- which(sapply(df, is.numeric)) #numeric型の変数を抽出する

numericVarNames <- names(numericVars)

all_numVar <- df[, numericVars]

cor_numVar <- cor(all_numVar, use="pairwise.complete.obs") 

# monthly_cost(家賃)との相関が高い順位ソート

cor_sorted <- as.matrix(sort(cor_numVar[,'monthly_cost'], decreasing = TRUE))

# 相関係数が絶対値で0.5超のものを抽出

CorHigh <- names(which(apply(cor_sorted, 1, function(x) abs(x)>0.5)))

cor_numVar <- cor_numVar[CorHigh, CorHigh]

# 相関図プロット

corrplot.mixed(cor_numVar, tl.col="black", tl.pos = "lt")







# 物件数が多い順位に並び替え

df1 <- df %>%

  select(city) %>%

  table %>%

  as.data.frame() %>%

  arrange(Freq)

df1$city <- with(df1,reorder(., Freq))



# 千葉県市別物件数のプロット

ggplot(df1,aes(city,n,fill = city)) +

  geom_bar(stat="identity", alpha=0.8)+

  ggtitle("千葉県市別物件数")+xlab("")+ylab("物件数")+ 

  theme(axis.text.x = element_text(angle = 180, hjust =1))+

  theme_bw(base_family = "HiraKakuProN-W3") +

  coord_flip() +

  ylim(0,25000) + scale_y_continuous(labels = scales::comma)



# 家賃の分布を見る

hist(df$monthly_cost,xlim = c(0, 200000),breaks = 30, col="#993435")



# 市別家賃BoxPlot

cost_median <- with(df, reorder(city, monthly_cost, median))

par(las=1, cex.axis=0.7, family = "HiraKakuProN-W3")

boxplot(monthly_cost ~ cost_median, data = df,

        xlab = "賃料+共益費", ylab = "",

        main = "市別家賃", varwidth = TRUE, horizontal=TRUE)



par(las=1, cex.axis=0.7, family = "HiraKakuProN-W3")

boxplot(monthly_cost ~ cost_median, data = df,

        xlab = "賃料+管理費", ylab = "",

        main = "市別家賃(〜20万円)", varwidth = TRUE, horizontal=TRUE, ylim = c(0,200000), outline=FALSE)



#市別築年数

age_med <- with(df, reorder(city, age, median))

par(las=1, cex.axis=0.7, family = "HiraKakuProN-W3")

boxplot(age ~ age_med, data = df,

        xlab = "築年数", ylab = "",

        main = "市別築年数", varwidth = TRUE, horizontal=TRUE, ylim = c(0,70), outline=FALSE)



# 建物構造別家賃

cost_median <- with(df, reorder(construction, monthly_cost, median))

par(las=1, cex.axis=0.7, family = "HiraKakuProN-W3")

boxplot(monthly_cost ~ cost_median, data = df,

        xlab = "賃料+共益費", ylab = "",

        main = "建物構造別家賃", varwidth = TRUE, horizontal=TRUE,ylim = c(0,200000))



# 家賃 VS 専有面積

plot(df$area,df$monthly_cost,

     xlim = c(0, 100),

     ylim = c(0, 200000),pch=".")



# 家賃 VS 専有面積

plot(df$age,df$monthly_cost,

     xlim = c(0, 60),

     ylim = c(0, 200000),pch=".")



# 家賃 VS 最寄り駅までの距離

plot(df$distance,df$monthly_cost,

     xlim = c(0, 100),

     ylim = c(0, 200000),pch=".")



# 家賃 VS 最寄り駅までの距離

plot(df$floor,df$monthly_cost,

     xlim = c(0, 10),

     ylim = c(0, 200000),pch=".")



# rout(最寄り駅)

df_rout <- df %>%

  select(rout1) %>%

  table %>%

  as.data.frame() %>%

  arrange(Freq)



# 構造

df_const <- df %>%

  select(construction) %>%

  table %>%

  as.data.frame() %>%

  arrange(Freq)

df_const



hist(df2$area, main = '専有面積')

hist(df2$age, main = '築年数')

hist(log(df2$distance), main = '最寄り駅距離(分)')

hist(df2$floor, main = '階')

hist(df2$height, main = '物件高さ(階建)')

hist(df2$free_rent, main = 'フリーレント実施有無')

hist(df2$monthly_cost, main = '賃料/月')



qqnorm(df2$area, main = '専有面積')

qqnorm(df2$age, main = '築年数')

qqnorm(log(df2$distance), main = '最寄り駅距離(分)')

qqnorm(df2$floor, main = '階')

qqnorm(df2$height, main = '物件高さ(階建)')

qqnorm(df2$free_rent, main = 'フリーレント実施有無')

qqnorm(df2$monthly_cost, main = '賃料/月')





#monthly_cost(家賃+共益費)を20万円以下に絞る

df2 <- df %>%

  filter(monthly_cost <= 200000) %>%

  filter(rout1 != '成田スカイアクセス') %>%

  filter(construction != 'ブロック') %>%

  filter(construction != '鉄骨プレ') %>%

  select(-c(shikikin, reikin, hoshoukin, station, direction))



df3 <- dummyVars(~.,data=df2[-1]) %>%

  predict(df2) %>%

  as.data.frame()



# train, test分割

index <- createDataPartition(df2$monthly_cost, p=.8, list=F)

train <- df2[index,]

test <- df2[-index,]

# dammy

train_dummy <- df3[index,]

test_dummy <- df3[-index,]

# randomforest学習用(上記trainだと時間がかかりすぎるため)

index2 <- createDataPartition(train$monthly_cost, p=.6, list=F)

train_rf <- train[index2,]



##### Linear Regression #####

lm.1 <- lm(monthly_cost~., data = train_dummy)

lm.step <- step(lm.1, trace = 1)

summary(lm.step)

lm_pred <- predict(lm.step, newdata = test_dummy)

lm_pred_csv <- data.frame(lm_pred)

write.csv(lm_pred_csv, "お好みのファイルパス")



test_sort <- test_dummy %>%

  arrange((monthly_cost))

lm_pred_sort <- predict(lm.step, test_sort)



ggplot(test_sort) + 

  geom_line(aes(x=1:nrow(test_sort),y=lm_pred_sort,color="red"))+ 

  geom_line(aes(x=1:nrow(test_sort),y=monthly_cost,color="blue"))+ 

  #geom_line(aes(x=1:nrow(test),y=df_lm_pred,color="green"))+

  ggtitle("テストデータ vs 予測データ(重回帰)")+xlab("index")+ylab("家賃")+ 

  scale_color_hue(name = "", labels = c("テストデータ","重回帰")) +

  theme(axis.text.x = element_text(angle = 180, hjust =1))+theme_bw(base_family = "HiraKakuProN-W3")+

  ylim(0,200000)



##### Random Forest #####



# tuneRFでグリッドサーチ。個々の決定木を作成する際に使用する特徴量の数mtryを求める

# system.time(

#   tuneRF(train_rf %>% select(-c(name,monthly_cost)), train_rf$monthly_cost,

#          doBest=TRUE, trace=TRUE, plot=TRUE )

#   )# 結果 → 8

rf <- randomForest(monthly_cost ~ ., data=train_rf[-1], mtry=8)

saveRDS(rf, file="rf_0906")

summary(rf)

plot(rf)

print(rf$importance)

varImpPlot(rf, main="Contribution",n.var = 10)



# 予測を行う

rf_pred <- predict(rf, test[-1]) #rent_admを削除

rf_pred_csv <- data.frame(rf_pred)

write.csv(rf_pred_csv, "お好みのファイルパス")



# テストデータと予測データをソートしておく

test_sort <- test[-1] %>%

  arrange((monthly_cost))

rf_pred_sort <- predict(rf, test_sort)



ggplot(test_sort) + 

  geom_line(aes(x=1:nrow(test_sort),y=rf_pred_sort,color="red"))+ 

  geom_line(aes(x=1:nrow(test_sort),y=monthly_cost,color="blue"))+ 

  #geom_line(aes(x=1:nrow(test),y=df_lm_pred,color="green"))+

  ggtitle("テストデータ vs 予測データ(ランダムフォレスト)")+xlab("index")+ylab("家賃")+ 

  scale_color_hue(name = "", labels = c("テストデータ","ランダムフォレスト") ) +

  theme(axis.text.x = element_text(angle = 180, hjust =1))+theme_bw(base_family = "HiraKakuProN-W3")+

  ylim(0,200000)



##### Random Forest end #####





##### RMSE #####

sqrt( sum((test$monthly_cost - lm_pred)^2) / length(test$monthly_cost) )

sqrt( sum((test$monthly_cost - rf_pred)^2) / length(test$monthly_cost) )





##### 最終予測 #####



df2$pred_rf <- predict(rf, newdata = df2)

df2$diff <- (df2$pred_rf - df2$monthly_cost)

df2$diff.monthly_cost <- df2$diff / df2$monthly_cost

df_sort <- df2 %>%

  arrange(desc(diff.monthly_cost))

kable(head(df_sort[c(1,2,3,4,6,7,8,9,14,15,16,17)],20))

 

~続・不動産とファイナンス・賃貸物件入居者編(2)~「機械学習を使って東京23区のお買い得賃貸物件を探してみた」を千葉県で再度やってみる(データクレンジング編)

今回は前回記事の続きで、SUUMOから取得してきたデータのクレンジングを行います。

d-s.hatenablog.com

 

 なお、このテーマはshokosakaさんの記事が参考となっています。

www.analyze-world.com

 

参考記事ではPythonでデータクレンジング(前処理)を行っていますが、私はR で実施しています。

Pythonでも似たようなコードになると思うので、Pythonユーザーの方でも問題ないと思いますが、ご希望の方はコメントくださればPythonのコードも書いて追記します。

 

では、さっそく内容に入っていきます。

前回の記事でスクレイピングしてきたデータをクレンジングした結果を先に見ておきます。

以下のようになりました。

> t(head(suumo_data,3))

             1                                 2                             3           

name         "ホーリーグランド稲毛山王305号室" "JR総武線稲毛駅4階建築52年" "大宮町貸家"

rent1        "30000"                           "50000"                       "40000"     

rent2        "0"                               "3000"                        "0"         

shikikin     "30000"                           "0"                           "40000"     

reikin       "30000"                           "0"                           "0"         

hoshoukin    "0"                               "0"                           "0"         

layout       "1K"                              "1DK"                         "1DK"       

area         "21.06"                           "34.27"                       "33.12"     

direction    "東"                              "東"                          "北西"      

type         "マンション"                      "マンション"                  "一戸建て"  

age          "89"                              "52"                          "48"        

rout1        "JR総武線"                      "JR総武線"                 "JR総武線"

station      "稲毛"                            "稲毛"                        "千葉"      

distance     " 72"                             " 25"                         "100"       

construction "鉄骨"                            "鉄筋コン"                    "木造"      

floor        "3"                               "2"                           "1"         

height       "4"                               "4"                           "1"         

car.dum      "敷地内"                          "近隣"                        "付無料"    

free_rent    "0"                               "1"                           "0"         

変数名は日本語をローマ字で充てたものと、簡単な英単語なので何かわかると思います。

補足をしておきますと以下の通りです。

area: 専有面積、age: 築年数、rout1: 最寄駅、distance: 最寄駅までの徒歩時間、

floor: その部屋が何階か、height: その物件は何階建か、car.dum: 駐車場有無、

free_rent: フリーレントの実施有無(期間は問わず)

 

以下が今回作成したRのコードです。

コードを解説しながら進めていこうと思ったのですが、コメントアウトを見ていただければそれなりにわかると思いますので、説明を省略してコードを載せてしまいます。

注意点としては、変数「distance(最寄り駅までの分数)」の取り扱いがあります。抽出した時点のデータセットでは徒歩とバスでの所要分数が入り混じっていました。

今回は、「最寄り駅まで徒歩で何分か」という基準で統一することにしたので、バス~分という場合は、バス~分に×4をして徒歩~分に置き換えています。

具体的には、バス10分であれば徒歩40分に置き換える。そんな感じです。

置き換えの根拠も載せてありますが、あくまで私なりの置き換え方法なので、何かアドバイスがある方はぜひ教えてください。

 

次回は、予測と物件リコメンドに入ります。

library(tidyverse)



# ファイルのPATHを指定して読み込む

strDirPath = 'スクレイピングしてきたデータが格納されているファイルPATHを各自記入'



# ファイルを一気に読み込んで一つのデータフレームにする

lf <- list.files(path = strDirPath, full.names = T)

data <- lapply(lf, read.csv, fileEncoding='cp932')

suumo_data <- do.call(rbind, data)



# 念のため、以下の列で重複がないか調べておく。

# 今回はデータ抽出後にLibreOfficeCalcで重複削除していたため重複はなかった。

suumo_data %>% distinct('name', 'rent1', 'rent2', 'layout', 'area','floor_height', .keep_all = TRUE)



# shiki_rei, access1, access2, access3, floor_heightについて分割処理を行う

suumo_data <- suumo_data %>%

  separate(shiki_rei, c('shikikin', 'reikin'), '/', remove = TRUE) %>%

  separate(access1, c('rout1', 'distance1'), '/', remove = TRUE) %>%

  separate(floor_heigth, c('floor', 'height'), '/', remove = TRUE)



# distance1はさらにstationtとdistanceに分割

suumo_data <- suumo_data %>%

  separate(distance1, c('station', 'distance'), '駅', remove = TRUE)



# 最寄駅にバス会社を選択している物件はdistanceがNAになる

# 最寄がバス会社で登録されている物件は登録数が非常に少ないかつ前提が大きく異なるような物件のためドロップ

suumo_data <- suumo_data[!is.na(suumo_data$distance),]



# rent2に含まれる文字列を0に置き換える

suumo_data$rent2 <- suumo_data$rent2 %>%

  str_replace('※定期借家', '') %>% 

  str_replace('-', '0')

# 数値で扱うデータは'-'を0にしておく

for (row in c(4, 5, 6, 7)){

  suumo_data[, row] <- suumo_data[, row] %>%

    str_replace('-', '0')

}



# rent1,shikikin,reikin,hoshoukin,shikibiki_shoukyakuに×10,000円して円単位にする

suumo_data$rent1 <- as.integer(suumo_data$rent1) * 10000

suumo_data$shikikin <- as.integer(suumo_data$shikikin) * 10000

suumo_data$reikin <- as.integer(suumo_data$reikin) * 10000

suumo_data$hoshoukin <- as.integer(suumo_data$hoshoukin) * 10000

suumo_data$shikibiki_shoukyaku <- as.integer(suumo_data$shikibiki_shoukyaku) * 10000



# rent2に含まれる文字列を0に置き換える

suumo_data$rent2 <- suumo_data$rent2 %>%

  str_replace('※定期借家', '') %>% 

  str_replace('-', '0')



# 数値で扱うデータは'-'を0にしておく

for (row in c(4, 5, 6, 7)){

  suumo_data[, row] <- suumo_data[, row] %>%

    str_replace('-', '0')

}



# 新築のageは「新」と表示されるため、0に置き換える

suumo_data$age[suumo_data$age == '新'] <- 0



# 平屋の物件はheightがNAになるため、1階建にする

suumo_data$height[is.na(suumo_data$height)] <- '1階建'



# floorの'平屋'は1階に変換→その後'階'を全部削除

# また、データ内の'0'と'建'と今回の'B'は誤りだと思われるため削除

suumo_data$floor <- suumo_data$floor %>%

  str_replace('平屋', '1階') %>% 

  str_replace('階', '') %>% 

  str_replace('0', '') %>% 

  str_replace('建', '') %>% 

  str_replace('B', '')



# 建物高さを抽出

# 地下は削除

suumo_data$height <- suumo_data$height %>%

  str_replace('地下1地上', '') %>% 

  str_replace('地下2地上', '') %>% 

  str_replace('地下3地上', '') %>% 

  str_replace('地下4地上', '') %>% 

  str_replace('地下5地上', '') %>% 

  str_replace('地下6地上', '') %>% 

  str_replace('地下7地上', '') %>% 

  str_replace('地下8地上', '') %>% 

  str_replace('地下9地上', '') %>% 

  str_replace('階建', '')



# kosuu,car,joukenの

						と

を除去(スクレイピング処理で誤って拾ってきた文字列の除去)

suumo_data$kosuu <- suumo_data$kosuu %>%

  str_replace('

						', '') %>% 

  str_replace('

', '')

suumo_data$car <- suumo_data$car %>%

  str_replace('

						', '') %>% 

  str_replace('

', '')

suumo_data$jouken <- suumo_data$jouken %>%

  str_replace('

						', '') %>% 

  str_replace('

', '') %>%

  str_replace('
						', '') %>% 

  str_replace('
', '')





# carを敷地内、近隣、付無料、無しに分類(金額情報は処理しない)

temp <- c()

for (row in suumo_data$car){

  if (str_detect(row, '敷地内')) {

    temp <- c(temp, '敷地内')

  } else if (str_detect(row, '近隣')){

    temp <- c(temp, '近隣')

  } else if (str_detect(row, '付無料')){

    temp <- c(temp, '付無料') 

  } else {

    temp <- c(temp, '無し') 

  }

}

suumo_data$car.dum <- temp



# フリーレントの実施有無

temp <- c()

for (row in suumo_data$jouken){

  if (str_detect(row, 'フリーレント')) {

    temp <- c(temp, 1)

  } else {

    temp <- c(temp, 0) 

  }

}

suumo_data$free_rent <- temp



# distanceの処理

# 最寄までバス~分または車~分の場合、徒歩~分に揃えるため×5する

# 一般に、徒歩80m/分。バスの場合、停車も加わるため時速20kmを分速換算し

# バス(車も)333m/分とする。333/80=4であるため、×4とする。家からバス停までの距離は無視。

suumo_data$distance_copy <- suumo_data$distance

suumo_data <- suumo_data %>%

  separate(distance, c('distance', 'unused'), '分', remove = TRUE)

suumo_data$distance <- suumo_data$distance %>%

  str_replace('バス', '') %>%

  str_replace('歩', '') %>%

  str_replace('車', '')

temp <- c()

for (i in rep(1:length(suumo_data$distance_copy))){

  if (str_detect(suumo_data$distance_copy[i],'バス')) {

    temp <- c(temp, as.integer(suumo_data$distance[i]) * 4)

  }  else if (str_detect(suumo_data$distance_copy[i], '車')) {

    temp <- c(temp, as.integer(suumo_data$distance[i]) * 4)

  } else {

    temp <- c(temp, as.integer(suumo_data$distance[i]))

  }

}

suumo_data$distance <- temp


# モデル作成で使いたい変数のみ残す tar <- c('name', 'rent1', 'rent2', 'shikikin', 'reikin', 'hoshoukin', 'layout', 'area', 'direction', 'type', 'age', 'rout1', 'station', 'distance', 'construction', 'floor', 'height', 'car.dum', 'free_rent' ) suumo_data <- suumo_data[tar] write.csv(suumo_data, 'お好みのファイル名', row.names = F, fileEncoding = 'cp932')

~続・不動産とファイナンス・賃貸物件入居者編(1)~「機械学習を使って東京23区のお買い得賃貸物件を探してみた」を千葉県で再度やってみる(Rでスクレイピング)

 前回の記事で、千葉県の賃貸物件の賃料予測を実行するところまで行いました。

d-s.hatenablog.com

 

しかし、データ取得のところで課題がいくつかあったため、今回データ取得(スクレイピング)のやり直しです。

コードはまとめて最下部に記載します。

前回はshokosakaさんの参考記事にのっとりPythonスクレイピングしましたが、予告通りRで完結させました。

 

データ取得

今回も、スクレイピング対象は下記のとおりです。

千葉県内で任意に選んだ市のSUUMO掲載賃貸物件データを取得

市川市 ・浦安市 ・柏市 ・鎌ヶ谷市 白井市 ・千葉市

流山市 ・習志野市 ・船橋市 ・松戸市 ・八千代市

 

前回は白井市のデータも取得していましたが、物件数が少なかったため今回は対象から外します。

また、間取りは単身者向けを前提に、ワンルーム・1K・1DK・1LDKでフィルターをかけます。

その他フィルター条件は特にかけていません。

 

追加事項

前回と大きく違うところは、「物件の構造」「向き」「駐車場有無」等物件にかかる詳細情報を多数盛り込んだことです。

おさらいになりますが、下の写真の通り「詳細を見る」をクリックした先にある・・・

f:id:d_s:20180828225321p:plain

 

この画面、詳細情報にアクセスしてデータを取得してきています。

f:id:d_s:20180828225236p:plain

f:id:d_s:20180828225502p:plain

実はまだ松戸市のデータのみスクレイピング真っ最中なのですが、いかんせん物件に1件1件アクセスしてデータを取得する関係で、スクレイピングマナーとしてスリープ時間を設けるととてつもない時間がかかります。

 

それで、取得したデータは以下のようなものになります。

サンプルで3物件見てみます。

 (変数名は適当でかっこ悪いですが。。。)

1 2 3
name アスピリアベルフルール Leaf(リーフ)103号室 Leaf01030号室
rent1 6.7 7.3 7.6
rent2 4200 3800 3500
shiki_rei -/6.7 -/7.3 -/7.6
hoshoukin - - -
shikibiki_shoukyaku - - -
layout 1LDK 1LDK 1LDK
area 44.67 44.13 44.13
direction 南西 南西 南西
type アパート アパート アパート
age 6 6 6
access1 新京成線/前原駅歩3分 JR武蔵野線/船橋法典駅歩10分 JR武蔵野線/船橋法典駅歩10分
access2 JR総武線/津田沼駅歩16分 JR総武線/西船橋駅バス17分(バス停)上山町歩6分 -
access3 京成本線/京成津田沼駅歩29分 東京メトロ東西線/西船橋駅バス17分(バス停)上山町歩6分 -
detail_layout 洋7.6LDK10.7 洋13LDK6.0 洋6LDK13.0
construction 木造 鉄骨 鉄骨
floor_heigth 1階/2階建 1階/2階建 1階/2階建
car 敷地内8640円 敷地内7560円 敷地内7560円
jouken 二人入居可 二人入居可/子供可/事務所利用不可 二人入居可
kosuu - - 8戸

 

ちゃんと、詳細情報が取れていますね。

あとは、これの前処理をしていって、再度基礎分析を行ってから予測に入っていきます。

 

本日のコード 

urlは適宜打ち換えて試してみてください。

時間がかなりかかるので、小規模で試してみると良いかもしれません。

また、スクレイピングの途中http Error500、http Error400に見舞われることがありましたが、気にせず再度実行して情報を取ってきています。

更に、データ取得中にあるページでHTMLの構造が変わる?のか、2パターンのHTMLを考慮に入れる必要があることに気づきました。

おそらくPythonで欠損になった物件はこれが原因かもしれません。

ですが、根本はよくわかっていません。

原因、回避の仕方等わかる方がいらっしゃいましたらご教授下さい。

(スクレイピングは初心者のためハードコーディングです)

library(rvest)

library(tidyverse)

library(readr)



##### 念のため文字列連結を簡単にする関数を定義しておく

"+" <- function(e1, e2) {

  if (is.character(c(e1, e2))) {

    paste(e1, e2, sep = "")

  } else {

    base::"+"(e1, e2)

  }

}



#URL(ここは各自入れ替えて。例は千葉市の賃貸住宅情報 検索結果の1ページ目・フィルター:ワンルーム。1K、1DK、1LDK)

#url <- 'https://suumo.jp/jj/chintai/ichiran/FR301FC001/?ar=030&bs=040&pc=30&smk=&po1=25&po2=99&shkr1=03&shkr2=03&shkr3=03&shkr4=03&ta=12&sa=01&cb=0.0&ct=9999999&md=01&md=02&md=03&md=04&et=9999999&mb=0&mt=9999999&cn=9999999&fw2='



fix_url <- 'https://suumo.jp'



df_all <- data_frame()





#データ取得

page <- read_html(url)



#ページ数を取得

pages <-  html_nodes(page, 'body') %>%

  html_nodes('div.pagination.pagination_set-nav') %>%

  html_text(trim = TRUE) %>%

  str_split('
|
|	')

pages <- as.integer(pages[[1]][4])



#URLを入れる配列

urls <- c()



#1ページ目を格納

urls <- c(urls, url)



#2ページ目から最後のページまでを格納

for (i in rep(1:pages)){

  pg = as.character(i+1)

  url_page = url + '&pn=' + pg

  urls <- c(urls, url_page)

}



count=0

#各ページで以下の動作をループ

for (url in urls) {

  # if (count==2) {

  #   break

  # }

  count <- count+1

  print(count)

  

  # データ取得

  tar_html <- read_html(url, encoding = 'utf-8')

  tar_url_list <- 

    tar_html %>% 

    html_nodes(xpath = "//a") %>% # aタグに格納されている

    html_attr("href") %>% # href属性のデータを取り出す

    as_data_frame() %>% 

    filter(grepl("/chintai/jnc_", .$value)) # 「詳細を見る」は"/chintai/jnc_"~で構成

  

  ## 詳細結果の数

  l <- nrow(tar_url_list)

  ## ここから各詳細結果へアクセスし、データを取得する 

  count2=0

  for (j in 1:l) {

    Sys.sleep(5) 

    count2 <- count2+1

    print(l+'-'+count2)

    

    ## 詳細結果へのURLを指定し、データを取得

    tar_url <- fix_url + tar_url_list$value[j]

    tar_html_tmp <- read_html(tar_url)

    

    name <- tar_html_tmp %>% 

      html_node('.section_h1-header-title') %>%

      html_text() %>%

      str_replace("
			","")

    

    rent1 <- tar_html_tmp %>%

      html_node(xpath = '//*[@id="js-view_gallery"]/div[1]/div[2]/div[2]/div/div[2]/div/div[1]/div/div[1]') %>%

      html_text() %>%

      str_replace("
										",'') %>%

      str_replace("万円",'')

    

    rent2 <- tar_html_tmp %>%

      html_node(xpath = '//*[@id="js-view_gallery"]/div[1]/div[2]/div[2]/div/div[2]/div/div[1]/div/div[2]/div/div[2]') %>%

      html_text() %>%

      str_replace('
												','') %>%

      str_replace("円",'')

    

    shiki_rei <- tar_html_tmp %>%

      html_node(xpath = '//*[@id="js-view_gallery"]/div[1]/div[2]/div[2]/div/div[2]/div/div[2]/ul/li[1]/div/div[2]') %>%

      html_text() %>%

      str_replace('
													
													','') %>%

      str_replace('

													','') %>%

      str_replace('

													
													','') %>%

      str_replace('
												','') %>%

      str_replace("万円",'')



    hoshoukin <- tar_html_tmp %>%

      html_node(xpath = '//*[@id="js-view_gallery"]/div[1]/div[2]/div[2]/div/div[2]/div/div[2]/ul/li[2]/div/div[2]') %>%

      html_text() %>%

      str_replace('
													
													','') %>%

      str_replace('

													','') %>%

      str_replace('

													
													','') %>%

      str_replace('
												','') %>%

      str_replace("万円",'')

    

    shikibiki_shoukyaku <-  tar_html_tmp %>%

      html_node(xpath = '//*[@id="js-view_gallery"]/div[1]/div[2]/div[2]/div/div[2]/div/div[2]/ul/li[3]/div/div[2]') %>%

      html_text() %>%

      str_replace('
													
													','') %>%

      str_replace('

													','') %>%

      str_replace('

													
													','') %>%

      str_replace('
												','') %>%

      str_replace("万円",'')



    if (is.na(rent1)){

      rent_table <- tar_html_tmp %>% 

        html_nodes('div.property_view_note-list') %>%

        html_nodes('span') %>%

        html_text()

      

      rent1 <- rent_table[1] %>%

        str_replace('万円','')

      

      rent2 <- rent_table[2] %>%

        str_replace('管理費・共益費:','') %>%

        str_replace('円','')

      

      shiki_rei <- (rent_table[3] +'/'+ rent_table[4]) %>%

        str_replace('敷金:','') %>%

        str_replace('礼金:','') %>%

        str_replace('万円','')

      

      hoshoukin <- (rent_table[5]) %>%

        str_replace('保証金:','') %>%

        str_replace('万円','')

      

      shikibiki_shoukyaku <- (rent_table[6]) %>%

        str_replace('敷引・償却:','') %>%

        str_replace('万円','')

      

    }

    

    

    detail_tbl1 <- tar_html_tmp %>%

      html_nodes('.property_view-detail') %>%

      html_nodes('td') %>%

      html_text()

    address <- detail_tbl1[1]

    layout <- detail_tbl1[3]

    area <- detail_tbl1[4] %>% 

      str_replace("m2",'')

    floor <-  detail_tbl1[6] %>% 

      str_replace("階",'')

    direction <- detail_tbl1[7]

    type <- detail_tbl1[8] 

    age <- detail_tbl1[5] %>% 

      str_replace("築",'') %>%

      str_replace("年",'') 



    if (is.na(address)){

    address <- tar_html_tmp %>% 

      html_nodes(xpath = '//*[@id="js-view_gallery"]/div[1]/div[2]/div[3]/div[2]/div[2]/div/div[2]/div') %>%

      html_text() %>%

      str_replace("
									","")

    layout <- tar_html_tmp %>% 

      html_nodes(xpath = '//*[@id="js-view_gallery"]/div[1]/div[2]/div[3]/div[1]/div/div[2]/ul/li[1]/div/div[2]') %>%

      html_text() %>%

      str_replace("
											","")

    area <- tar_html_tmp %>% 

      html_nodes(xpath = '//*[@id="js-view_gallery"]/div[1]/div[2]/div[3]/div[1]/div/div[2]/ul/li[2]/div/div[2]') %>%

      html_text() %>%

      str_replace("
											","") %>%

      str_replace("m2","")

    direction <- tar_html_tmp %>% 

      html_nodes(xpath = '//*[@id="js-view_gallery"]/div[1]/div[2]/div[3]/div[1]/div/div[2]/ul/li[3]/div/div[2]') %>%

      html_text() %>%

      str_replace("
											","")

    type <- tar_html_tmp %>% 

      html_nodes(xpath = '//*[@id="js-view_gallery"]/div[1]/div[2]/div[3]/div[1]/div/div[2]/ul/li[4]/div/div[2]') %>%

      html_text() %>%

      str_replace("
											","")

    age <- tar_html_tmp %>% 

      html_nodes(xpath = '//*[@id="js-view_gallery"]/div[1]/div[2]/div[3]/div[1]/div/div[2]/ul/li[5]/div/div[2]') %>%

      html_text() %>%

      str_replace("
											","") %>%

      str_replace("築","") %>%

      str_replace("年","")

    }



    access_tbl <- tar_html_tmp %>%

      html_nodes('.property_view-detail') %>%

      html_nodes('td') %>%

      html_nodes('div') %>%

      html_text()

    if (sum(is.na(access_tbl))==1|length(access_tbl)==0) {

      access_tbl <- tar_html_tmp %>% 

        html_node('.property_view_detail.property_view_detail--train') %>%

        html_nodes('.property_view_detail-text') %>%

        html_text()

    }

    if (length(access_tbl) == 0){

      access1 <- '-'

      access2 <- '-'

      access3 <- '-'

    } else if (length(access_tbl) == 1){

      access1 <- access_tbl[1]

      access2 <- '-'

      access3 <- '-'

    }else if (length(access_tbl) == 2){

      access1 <- access_tbl[1]

      access2 <- access_tbl[2]

      access3 <- '-'

    }else if (length(access_tbl) == 3){

      access1 <- access_tbl[1]

      access2 <- access_tbl[2]

      access3 <- access_tbl[3]

    }  

    

    detail_tbl2 <- tar_html_tmp %>% 

      html_nodes('.data_table.table_gaiyou') %>%

      html_nodes('td')

  

    detail_layout <- detail_tbl2[1] %>% 

      

      

      html_text() %>%

      str_replace("
						","")

    

    construction <- detail_tbl2[2] %>% 

      html_text() %>%

      str_replace("
						","")

    

    floor_heigth <- detail_tbl2[3] %>% 

      html_text()

    

    car <- detail_tbl2[6] %>% 

      html_text()

    

    jouken <- detail_tbl2[9] %>% 

      html_text()

    

    kosuu <- detail_tbl2[12] %>% 

      html_text()

    

    temp.DF <- data.frame(name, rent1, rent2, shiki_rei, hoshoukin,

                  shikibiki_shoukyaku, layout, area, direction,

                  type, age, access1, access2, access3, detail_layout,

                  construction, floor_heigth, car, jouken, kosuu)

    df_all <- rbind(df_all, temp.DF)

  }  

}

write.csv(df_all, "各々のPATHを入力", quote=T, row.names = F)

 

 

Rの【woeBinning】パッケージ がとても便利

分類問題をロジスティック回帰で予測しようと思ったらこのパッケージが便利そうです。

説明変数をWOE(Weight ou Evidence)ベースでビン化してIV(Information Value)

を算出してくれます。

「WOE」というよりクレジットスコアリングモデルの説明ですが、以下が参考になると思います。

https://www.worldprogramming.com/jp/blog/credit_scoring_pt5

 

説明変数のビン化(クラス化)についてはやり方がたくさんあるのでしょうが、春に受けたSASのクレジットリスクセミナーでもWOEについて時間を割いて解説されていました。

PythonでもWOEベースのビン化パッケージがあるのですが、Rもpythonも日本語で解説しているページがあまり見当たらないので紹介しておきます。

 

組み込みのGerman Creditデータセットを使用します。

Good/Badで結果および20の変数がついた1000行21列のデータフレームです。

library(woeBinning)

# German creditデータセットを読み込む
data(germancredit)
df <- germancredit
summary(df)
str(df)
kable(head(df))

# 説明変数をビン化、目的変をセット
binning <- woe.binning(df, target.var = 'creditability', pred.var = df)

# ビン化した変数をプロット
woe.binning.plot(binning, multiple.plots=F)

f:id:d_s:20180829221031p:plain

f:id:d_s:20180829230729p:plain

実行すると、IVベースでの変数ランキングと、各変数ごとの変数プロットが出てきます。

変数分出てくるので、個別変数としてはduration in monthだけ載せておきます。

このBinning関数には渡せるパラメーターまだがいくつかあって、

stop.limit:ビン化により減少するIVの減少幅にリミットを設ける

min.perc.total:各ビンに入る目的変数の数の割合に下限を設ける

など柔軟にチューニングできそうです。

その他はHelpをご参照ください。

 

duration in monthは以下の通りビン化されています。

binning[3,2]
[[1]]
                 woe cutpoints.final cutpoints.final[-1] iv.total.final bad good col.perc.a col.perc.b
(-Inf,6]  -124.59370            -Inf                   6      0.2537678   9   73  0.0300000  0.1042857
(6,15]     -36.53869               6                  15      0.2537678  80  269  0.2666667  0.3842857
(15,30]     10.83411              15                  30      0.2537678 128  268  0.4266667  0.3828571
(30, Inf]   76.63288              30                 Inf      0.2537678  83   90  0.2766667  0.1285714
Missing           NA             Inf             Missing      0.2537678   0    0  0.0000000  0.0000000
              iv.bins
(-Inf,6]  0.092555320
(6,15]    0.042976457
(15,30]   0.004746374
(30, Inf] 0.113489646
Missing            NA
hist(df$duration.in.month)

 

このdurationですが、binning関数にかける前のヒストグラムはこんな感じです。

 

f:id:d_s:20180829223624p:plain

 

近々、KaggleでみつけたHumanResourceAnalytics(会社をやめる人を予測する)データセットで、このビン化モジュールで前処理、ロジスティック回帰で予測を実施したいと思います。

 

その際に詳しく解説します。

 

 

 

~不動産とファイナンス・賃貸物件入居者編(2)~「機械学習を使って東京23区のお買い得賃貸物件を探してみた」を千葉県でやってみる

今回は前回の続きをやっていきます。

引き続き参考はこちらです。

www.analyze-world.com

 

参考のままですが、コードは追記しておきました。

では、まず前回見ていなかった築年数から見ていきます。

f:id:d_s:20180827221357p:plain

市川、浦安という東京(江戸川区)の隣接地域の築年数が高く、流山・鎌ヶ谷築年数が低いという結果が出ました。

流山、鎌ヶ谷については近年開発されている地域なので納得感がありますが、千葉は意外でした。

よくよく考えると、物件情報が掲載されている時点で、空室・または退去見込みの物件が掲載されているわけです。

よって、市川・浦安等の東京へのアクセスが良いエリアは、築浅物件ほど埋まっているのかもしれません。

「年数の経った物件が残りやすい→掲載物件ベースでは築年数が高く出やすい」という流れがあるのでしょうか?

こんな感じで物件掲載情報での集計だとバイアスがありそうなので、解釈には注意が必要だと思います。

 

予測

ここからが一番ワクワクします。 

はじめに結論ですが、まずまず予測できたと感じるものの。。。

データセット内で、物件名と物件内容が一致していないものを多数見つけました!!

前回の分析中で、欠損を含むデータをドロップして気がかりで・・・と話しましたが、やはり嫌な予感が的中しました。

どうやら、スクレイピングの最中に何らかの原因で物件名と物件の内容にズレが起きていたようです。

物件内容自体は間違い無いと思うのですが、やはりやり直した方が良さそうです。

 

ということで、スクレイピングするところからやり直すこととします。

とは言ってもせっかくなので、一応の結果は挙げておきます。参考程度ということで。

 以下ランダムフォレストの結果です。

f:id:d_s:20180828001255p:plain

パラメーターは、tuneRF関数でグリッドサーチしました。個々の決定木を作成する際に使用する特徴量の数「mtry」を求めた結果、6となりましたので、mtry=6でモデルを作成しました。

急速に収束していく様子がわかります。

次にランダムフォレストで学習した結果の変数(特徴量)の重要度を見てみます。

f:id:d_s:20180828001309p:plain

other_cost(敷金、礼金、保証金等)とarea(専有面積)、age(築年数)の重要度が高いことがわかります。以降、最寄駅、市、間取り1LDK、階、最寄り駅までの距離、間取り1K・・・と続きます。

間取りはダミー化してしまっているので不利な気もしますが、全体的に納得できるような腑に落ちない何かがあるような。。。

なお、RandomForestの特徴量の重要度算出について以下の記事が参考になります。

Random Forestで計算できる特徴量の重要度 - なにメモ

 

つづいて、予測の結果です。

f:id:d_s:20180828002941p:plain

データセット全体を、学習用:テスト用に8:2に分けて検証を行っています。

赤線のテストデータに対し、ランダムフォレストの結果が青ですので、それなりに予測できていそうですが、先に述べたようにデータセットに留意が必要です。

物件名は一応伏せてありますが、アウトプットの中ざっと見た中で気になった物件を一つ挙げてみます。

f:id:d_s:20180828141747p:plain

どうでしょうか?

こちらの物件は築古の1Kですが、津田沼駅徒歩7分の立地で、リフォームしてあるようです。

敷金礼金なしで、管理費共益費込41,000円は安いのでないでしょうか。

再チャレンジが必要になりましたが、一連のアウトプットを経て勉強になりました。

まとめ

今回の分析を通じて、どうしても追加したい思った変数がありました。 

 

 

物件の構造です。

 

 

特徴量の重要度を見て何か腑に落ちない気がしたのは、まだ追加したいと思う変数が複数あったからです。

木造か軽鉄か、または鉄筋コンクリートなのかって物件選びで大事なポイントですよね。

また、物件の方角や駐車場有無なんかもあった方がより予測精度が上がるかもしれません。

総戸数がわかれば、その物件の入居率も変数として投入できます。

 

ということで、改めてsuumoを見てみます。

f:id:d_s:20180828225321p:plain

 

 

そういえば、今回「詳細を見る」の画面を気にもしていませんでした。

スクレイピングしてきたのは、この一覧画面で得られる情報に限定されていたわけですので、「詳細を見る」の中を見てみます。

f:id:d_s:20180828225236p:plain

 

f:id:d_s:20180828225502p:plain

あるじゃないですか!!

物件の向きから構造まで、そうそうこれが欲しかったんです。

実行時間としては相当な時間を覚悟しますが、物件詳細をループで処理していけば何とかなる気がしてきました。

 

よ~し、どうせスクレイピングし直すなら全部入れてしまえ!!

ということで、全項目取得します。

 

前回、Rでのスクレイピングを途中で断念してしまったので、あえて苦しみますがRでやり直そうと思います。

Pythonに逃げるかもしれませんが、少しでもRに貢献したく。。。)

まだ途中ですが、ざっとこんな感じのデータフレームができそうです。

f:id:d_s:20180828232000p:plain

おおちゃくしてRstudioのキャプチャをのっけてしまい横が収まっていないですが、建物構造やその物件にかかる条件なんかも取得できています。

その他取得できる項目を全部変数化できればより面白くなりそうです。

まだ恥ずかしながらエラー地獄に陥っているためコードを公開できる状態ではありませんが、完成したら続編記事内で載せたいと思います。

スクレイピングはド素人ですが、Rで奮闘します。

 

 

やり出すと、何だかやり込みたくなってきました。

 

 

 次回につづきます。

 

 

#市別築年数
age_med <- with(df, reorder(city, age, median))
par(las=1, cex.axis=0.7, family = "HiraKakuProN-W3")
boxplot(age ~ age_med, data = df,
        xlab = "築年数", ylab = "",
        main = "市別築年数", varwidth = TRUE, horizontal=TRUE, ylim = c(0,70), outline=FALSE)

# モデルに投入前する変数を選択
select_var1 <- c('names',"city", 'other_cost', "age", "floor", "area",
                 "monthly_cost", "rout1",'distance1','layout')
df2 <- df %>%
  filter(monthly_cost <= 170000) %>%
  select(select_var1)

# 相関を見ておく
cor(df2[,-c(1,2,8,10)])

# train, test分割
index <- createDataPartition(df2$monthly_cost, p=.8, list=F)
train <- df2[index,]
test <- df2[-index,]

# tuneRFでグリッドサーチ。個々の決定木を作成する際に使用する特徴量の数mtryを求める
tuneRF(train[,-1], train[,-1], doBest=T) # 結果→6
rf <- randomForest(monthly_cost~., data=train[,-1], mtry=6)
saveRDS(rf, file="rf")
summary(rf)
plot(rf)
print(rf$importance)
varImpPlot(rf, main="Contribution")

# 予測を行う
rf_pred <- predict(rf, test[-1]) 
rf_pred_csv <- data.frame(rf_pred)
write.csv(rf_pred_csv, "C:\\Users\\daich\\Documents\\R_src\\realestate\\rf_pred.csv")

# テストデータと予測データをソートしておく
test_sort <- test %>%
  arrange((monthly_cost))
rf_pred_sort <- predict(rf, test_sort[-1])

ggplot(test_sort) + 
  geom_line(aes(x=1:nrow(test_sort),y=rf_pred_sort,color="red"))+ 
  geom_line(aes(x=1:nrow(test_sort),y=monthly_cost,color="blue"))+ 
  #geom_line(aes(x=1:nrow(test),y=df_lm_pred,color="green"))+
  ggtitle("テストデータ vs 予測データ(ランダムフォレスト)")+xlab("index")+ylab("家賃")+ 
  scale_color_hue(name = "", labels = c("テストデータ","ランダムフォレスト") ) +
  theme(axis.text.x = element_text(angle = 180, hjust =1))+theme_bw(base_family = "HiraKakuProN-W3")+
  ylim(0,200000)


df2$pred_rf <- predict(rf, newdata = df2)
df2$diff <- df2$pred_rf - df2$monthly_cost
df_sort <- df2 %>%
  arrange(desc(diff))
df_sort <- df_sort[,c(1,2,10,7,11,12,3,8,9,4,5,6)]
df_sort$names <- rep(1:length(df_sort$names))
colnames(df_sort) <- c('物件名','市','間取り','1.家賃','2.予想家賃','予想家賃-家賃',
                       '敷金/礼金/保証金','最寄り駅','最寄駅距離','築年数','階','専有面積')

~不動産とファイナンス・賃貸物件入居者編(1)~「機械学習を使って東京23区のお買い得賃貸物件を探してみた」を千葉県でやってみる

今回は前々回の記事で触れた通り、以下のエントリを参考に千葉県のお買い得賃貸物件を分析して探してみます。

shokosakaさんのブログですが、大変勉強になりました。

このブログをなぞる形で進めてみます。

www.analyze-world.com

 

その前に、ブログを始めてわかったのですが、高校時代にサッカー部で共に汗を流したK君がこんな面白い記事を書いていたので合わせて紹介しておきます。

K君は僕と違い物理方面の科学者(FaceBookによるとナノサイエンス?)なわけですが、経済・ファイナンス界隈の人が一度は経験するラグランジュ未定乗数法を直感的なプロットで説明してくれています。彼はいつの間にかフランスの大学で研究をやっていてフランスからレスポンスしてくれました。

また、matplotlibでCDジャケットを作ってしまうという暴挙も見せています。

koideforest.hatenadiary.com

koideforest.hatenadiary.com

 

 

さて、本題に入ります。

 

 

賃貸物件情報の取得条件

まず初めに、例にならってsuumoから賃貸物件情報を取得してきます。

もちろん千葉県を対象にするのですが、さすがに全エリアを対象にするわけにはいきませんので、以下の市に絞ってデータを取得してきています。(東京に近いところから選択)

市川市 ・浦安市 ・柏市 ・鎌ヶ谷市 ・白井市 ・千葉市

流山市 ・習志野市 ・船橋市 ・松戸市 ・八千代市

また、ある程度前提を絞りたいので、さらに単身者向け物件に限定します。

検索条件として、間取りを「ワンルーム、1K、1DK、1LDK」でフィルターしています。その他の条件は特に付与していません。

 

 賃貸物件情報の取得

僕は前処理からRで完結したいと思ったのでrvestを使ってスクレイピングしてこようと思ったのですが、何やら途中うまくいかなくなったので結局pythonで処理しました。

今回はBeautifulSoup4を使ってスクレイピングしています。

 

 結果、欠損やら何やらざっと処理した結果、最終的に28,445件の物件情報になりました。

データ取得面での前処理はざざっと終わらせましたが、正直欠損をドロップする時点で気がかりなこともありました。

いろいろ追求したいところですが、試行ですので先に進みます。

 

基礎統計

ここからが分析のスタートです。

実務では、先のデータ取得、データのクレンジング、基礎統計で外観を把握、問題あれば再度データのクレンジング・・・と途方もないループがあって一番時間がかかるところではあります。

まず、市別にSUMMO掲載物件数を見ていきます。

f:id:d_s:20180826231338p:plain

最大は千葉市で約9400件、最低は白井市で125件です。

shokosakaさんの記事のまんまですが、言わずもがな人口の多い都市順に並んでいそうです。

では、物件数と人口でプロットして確かめます。

(人口は2015年)

f:id:d_s:20180826233324p:plain

おお、綺麗な回帰直線が引けました。

やはり都市人口と物件掲載数は強い相関があることがわかりますね。

習志野市だけ信頼区間から飛び出ていますが、物件の供給過剰なのでしょうか?

昨今、相続対策と消費税増税を謳ったアパート建築が流行ったので、築年数を分析すれば見えてくるのかもしれません。

 

さて、入居者としては一番の関心事である家賃を見ていきます。

家賃の中央値が高い順に並べます。

f:id:d_s:20180827000041p:plain

これは納得ですね。どうやら都心に近い順に並んでいそうです。

流山が3番目に来たのは 意外でしたが、つくばTX沿線の開発で新しい物件が多いのかもしれません。

以前、流山おおたかの森SCのTOHOシネマズによく映画を見に行きましたが、住環境としては文句のつけようのないような場所だと感じていました。

逆に、流山市の家賃の最低価格は相当安いので旧市街・新市街で相当な差があるのでしょう。

また、市川、千葉、松戸についてはどうやら変なデータが紛れ込んでいそうです。

単身向け物件にしては相当な値段の物件が見られますが、間取りの設定ミスか値段の設定ミスでしょう。

モデルを作成する際には注意しなければなりません。

 

次に、外れ値を除外して同様に見てみます。

f:id:d_s:20180827000107p:plain

先ほどよりだいぶ見やすくなりました。

白井市については、そもそも物件数が125件程度なのでなんとも言えませんが、流山市はレンジが広いですね。

そして今更気づきましたが、不覚にも柏市を入れるのを忘れていました。。。

(知らない方のために言いますと、柏市柏レイソルのホームタウンで千葉県内でも人気のある都市です)

 

やっとここまできましたが、だいぶ遅い時間になってしまいました。

その他、基礎分析を一気に行ってしまいます。

Bar Chart(by Frequency)

f:id:d_s:20180827001308p:plain

グラフが小さくと申し訳ありませんが、要約すると

・間取り

1K→1LDK→ワンルーム?→1DK→その他

・最寄り駅

JR総武線(圧倒的)→東京メトロ東西線→JR常磐線京成本線新京成線→JR京葉線→その他

・最寄り駅徒歩(分)

10分→5分→7分→8分→3分→その他

(その他、築年数等の分析を行いたいところですが今日はここまでということで。)

念のため変数とその欠損割合を確認して今日は終わりにしたいと思います。

f:id:d_s:20180827005950p:plain

次回はいよいよモデルを作成して割安物件を探していきます。

 

※コード追記(こちらのブログのほぼまんまです)

機械学習を使って東京23区のお買い得賃貸物件を探してみた - データで見る世界

 

library(ggplot2)
library(tidyverse)
library(caret)
library(randomForest)
library(ggrepel)
library(psych)
library(makedummies)
library(DT)


options(scipen=100)

# データセットの読み込み
df <- read.csv("C:\\Users\\daichi_saito\\Documents\\realestate\\dataset_chibapre2.csv",
               fileEncoding = "cp932", sep = ",")
colnames(df) <- c('names', 'city', 'adress', 'age', 'layout','height',
                  'floor', 'area', 'monthly_cost', 'other_cost', 'rout1',
                  'station1', 'distance1', 'rout2', 'station2', 'distance2',
                  'rout3', 'station3', 'distance3')
df <- cbind(df,makedummies(df, basal_level = FALSE, col = "layout"))
df$distance1 <- as.integer(str_replace(df$distance1,'分',''))
df <- na.omit(df)
df$city <- as.character(df$city)
df$city <- factor(df$city, levels=unique(df$city))
df$layout <- as.character(df$layout)
df$layout <- factor(df$layout, levels=unique(df$layout))

# 物件数が多い順位に並び替え
df1 <- df %>%
  select(city) %>%
  group_by(city) %>%
  tally() %>%
  arrange((n))
df1$city <- with(df1,reorder(city, n))

# 千葉県市別物件数のプロット
ggplot(df1,aes(city,n,fill = city)) +
  geom_bar(stat="identity", alpha=0.8)+
  ggtitle("千葉県市別物件数")+xlab("")+ylab("物件数")+ 
  theme(axis.text.x = element_text(angle = 180, hjust =1))+
  theme_bw(base_family = "HiraKakuProN-W3") +
  coord_flip() +
  ylim(0,25000) + scale_y_continuous(labels = scales::comma)

# 人口データを追加して物件数/人口のプロット
df1$pop <- c(62,109,164,193,174,168,483,482,623,972)
ggplot(df1,aes(pop,n)) +
  geom_point()+
  ggtitle("物件数/人口のプロット")+xlab("人口(千人)")+ylab("物件数")+ geom_smooth(method = "lm") +
  geom_text_repel(label = df1$city, family = "HiraKakuProN-W3", size = 3) +
  theme_bw(base_family = "HiraKakuProN-W3") +
  ylim(0,25000) + scale_y_continuous(labels = scales::comma) +  scale_x_continuous(labels = scales::comma)
cor(df1$n, df1$pop)

# 市別家賃BoxPlot
cost_median <- with(df, reorder(city, monthly_cost, median))
par(las=1, cex.axis=0.7, family = "HiraKakuProN-W3")
boxplot(monthly_cost ~ cost_median, data = df,
        xlab = "賃料+共益費", ylab = "",
        main = "市別家賃", varwidth = TRUE, horizontal=TRUE)

par(las=1, cex.axis=0.7, family = "HiraKakuProN-W3")
boxplot(monthly_cost ~ cost_median, data = df,
        xlab = "賃料+管理費", ylab = "",
        main = "市別家賃(〜15万円)", varwidth = TRUE, horizontal=TRUE, ylim = c(0,150000), outline=FALSE)

Rと効率的フロンティアとCAPM(2)

前回の続きです。

早速効率的フロンティアを描いてみます。

dt2 <- 

  dt[,c(1,2,5)] %>% na.omit() %>%

  tidyr::spread(ticker, ret) %>%

  select(2,3,4)



# 期待リターン(平均)の算出

er_x <- mean(dt2$SOFTBANK)

er_y <- mean(dt2$TAKEDA)

er_z <- mean(dt2$TOYOTA)



# ボラティリティ(標準偏差)の算出

sd_x <- sd(dt2$SOFTBANK)

sd_y <- sd(dt2$TAKEDA)

sd_z <- sd(dt2$TOYOTA)



# 共分散の算出

cov_xy <- cov(dt2$SOFTBANK, dt2$TAKEDA)

cov_xz <- cov(dt2$SOFTBANK, dt2$TOYOTA)

cov_yz <- cov(dt2$TAKEDA, dt2$TOYOTA)



# ポートフォリオのウェイトを出しておく

x_weights <- seq(from = 0, to = 1, length.out = 1000)



# ウェイト付けしたデータフレームを生成

three_assets <- data.table(wx = rep(x_weights, each = length(x_weights)),

                           wy = rep(x_weights, length(x_weights)))



three_assets[, wz := 1 - wx - wy]



three_assets[, ':=' (er_p = wx * er_x + wy * er_y + wz * er_z,

                     sd_p = sqrt(wx^2 * sd_x^2 +

                                   wy^2 * sd_y^2 +

                                   wz^2 * sd_z^2 +

                                   2 * wx * wy * cov_xy +

                                   2 * wx * wz * cov_xz +

                                   2 * wy * wz * cov_yz))]



# 今回はショート戦略はなしで考えたいので0以上のみ

three_assets <- three_assets[wx >= 0 & wy >= 0 & wz >= 0]



# plot

ggplot() +

  geom_point(data = three_assets, aes(x = sd_p, y = er_p, color = wx - wz)) +

  geom_point(data = data.table(sd = c(sd_x, sd_y, sd_z), mean = c(er_x, er_y, er_z)),

             aes(x = sd, y = mean), color = "red", size = 3, shape = 18) +

  theme_bw() + ggtitle("効率的フロンティア") +

  xlab("Volatility") + ylab("Expected Returns") +

  scale_y_continuous(label = percent) +

  scale_x_continuous(label = percent) +

  scale_color_gradientn(colors = c("red", "blue", "green"),

                        name = expression(omega[x] - omega[z]), labels = percent)

これを実行すると以下のプロットが出ます。

f:id:d_s:20180825020057p:plain

 一応ですが出ました。

 3つの株式の組み合わせで得られるリスク、期待リターンが色付きの箇所になります。

これに、銘柄を限定せず全銘柄で効率的フロンティアを作りCAPM(資本資産価格)というものを組み合わせると無リスク資産(国債等、日本の現状では現預金と言い換えても?)を加えた戦略が得られます。

今回は例が悪かったので他から図を拝借してきますが以下のような感じです。

 

http://words.equity-investment.info/img/%E6%9C%89%E5%8A%B9%E3%83%95%E3%83%AD%E3%83%B3%E3%83%86%E3%82%A3%E3%82%A22.jpg

参照:http://words.equity-investment.info/portfolio.theory.html

効率的(有効)フロンティアのカーブが最小になる点を最小分散ポートフォリオ、資本市場線との接点を市場ポートフォリオと言います。

この市場ポートフォリオの正体こそが厳しい厳しい仮定の下得られたCAPMでして、

すべての投資家はこの市場ポートフォリオと無リスク資産を持つとしています。

市場ポートフォリオは、市場の関連情報すべてを折り込んでいるため、投資家は市場ポートフォリオ保有することで最適なポートフォリオを保持できるというものです。 

これについてはもちろんいろんな意見がありCAPMも拡張されたりするわけです。

 

まとめ

長くなりましたが、これまでやってきた流れを乱暴な言い方で締めくくると、

「現金とTOPIX上場投信で資産配分考えておけばOK!!!!!!」と言えるのでしょうか。

あくまで理論的な話ですが。

次回以降は、前回冒頭あげたような身近なネタで分析していきます。