Чтобы продать что-нибудь ненужное, нужно сначала купить что-нибудь ненужное, а у нас денег нет.
— Трое из Простоквашино

Введение

Так получилось, что я живу в своей квартире (или кондо по-местному) в Монреале. И однажды, примерно год назад меня посетила мысль что неплохо бы перебраться в собственный дом. Некоторый опыт покупки и продажи жилья у меня уже был и, в принципе, можно было бы подойти к этому вопросу просто, как поступает большинство местных обывателей: нанять риэлтора и предоставить ему разобраться со всеми вопросами, но это было бы скучно и неинтересно.

Поэтому я решил подойти к этому делу научно. Есть задача: надо разобраться сколько примерно стоит то что у меня есть, и где находится то что я могу себе позволить. Ну и попутный вопрос — понять куда дует ветер. И изучить гео-пространственные вычисления в R.

В принципе, сразу было понятно, что просто отдельный дом ( single family house по-местному) я не потяну, если хочу остаться в цивилизованной местности и ударять ежедневным велопробегом по глобальному потеплению. Другой распространённый местный вариант — купить дуплекс или триплекс, т.е. дома, где две или три квартиры: в одной живешь, в остальных разводишь кроликов остальные сдаются жильцам. Тут появляется ещё одна неизвестная величина — доход от аренды.

Поэтому захотелось сделать карту города с ценами жилья на продажу, ценами на аренду, ну и ещё иметь возможность отслеживать, как это всё изменяется со временем.

Сервисы типа zillow, показывающие красивые картинки с графиками расчетной цены жилья, для нашей местности не работают, и изучать можно только регулярные отчёты местной ассоциации риэлторов, которые регулярно рапортуют, что цены идут вверх, аж удержу нет: https://apciq.ca/en/real-estate-market/. Цены на аренду тоже никто нигде массово не анализирует, можно просто помониторить доски объявлений, но это занятие довольно скучное.

В общем, первым делом я прошерстил интернет на тему готовых открытых данных о ценах на местную недвижимость, почти ничего не нашёл, зато обнаружил, что местный сервис листинга жилья на продажу имеет слабо документированный интерфейс, по которому можно таскать информацию о том, что сейчас доступно на рынке: https://github.com/Froren/realtorca

После этого дело оставалось за малым — был написан скрипт на питоне, который каждый день проходится по острову Монреаль и окрестностям и собирает информацию о рынке, дополнительно был сделан веб-скраппер на основе requests и beatifulsoap, он тоже ежедневно собирает объявления о жилье на съем.

Оба сервиса дают ограниченную информацию — риэлторский листинг показывает геолокацию, цену которую просят, тип жилья, количество спален, иногда размер жилья и земельного участка, этажность; доска объявлений о жилье на аренду показывает только цену, геолокацию и количество спален остальное опционально.

К счастью, нам на помощь приходит портал открытых данных города и openstreet map, через которые можно заполнить некоторые пробелы.

Обработка данных

Вебскраппер — довольно примитивный скрипт, каждый день он делает дамп всего, что может найти в sqlite базу данных, таким образом можно отслеживать как недвижимость появляется на рынке, как меняется её цена и когда она исчезает. Чего нельзя сделать, так это узнать, за сколько именно её продали, сколько было посетителей, в каком состоянии дом и т.д.

После скраппера данные обрабатываются в R, с использованием пакетов из набора tidy-verse, Simple Features for R, довольно ценный ресурс о том как это всё использовать — он-лайн книжка Geocomputation with R, графики нарисованы с помощью ggplot2 ( входит в tidyverse), и tmap.

Сперва я делаю всякий банальные вещи, типа преобразования систем координат в местную топологическую систему, и пространственное пересечение (join?) с базой данных недвижимости от города.

Отслеживание локальной цены квартиры

Для начала, я сделал простенькую статистику по ценам квартир в моём квартале, с помощью dplyr отбираем только кондоминимумы, а потом делаем пересечение с географическим описанием моего квартала:

Для начала надо сказать R что мы имеем дело с географическими данными, и каждая квартира это сферический конь в вакууме:

library(tidyverse)
library(sf)

property<-read_csv("....") %>% 
 st_as_sf(coords=c("lng","lat"), crs=4326) %>% 
 st_transform(crs=32188)

Выбираем район города в котором я живу:

neighbourhood<-geojson_sf("quartierreferencehabitation.geojson") %>%
 st_transform(32188) %>% 
 filter(nom_qr %in% c("Saint-Louis", "Milton-Parc")) %>% 
 summarize() %>% 
 st_buffer(dist=0)

И делаем пересечение со списком квартир:

neighbors <- st_join(property, neighbourhood, left=F)

Теперь можно взять карту с openstreetmap и показать с чем мы имеем дело:

osm_neighbourhood<-read_osm(st_bbox(neighbourhood%>%st_transform(4326)), ext=1.5, type="esri")

Изображаем карту и квартиры с помощью пакета tmap :

library(tmap)
library(tmaptools)

tm_shape(osm_neighbourhood) + tm_rgb(alpha=0.7)+
  tm_shape(neighbourhood) + tm_borders(col='red',alpha=0.8)  + 
  tm_shape(neighbors) + tm_symbols(shape=3,size=0.2,alpha=0.8) +
  tm_shape(ref_home) + tm_symbols(col='red',shape=4,size=0.5,alpha=0.8)+
  tm_compass(position=c("right", "bottom"))+
  tm_scale_bar(position=c("right", "bottom"))

Теперь посмотрим на цены, в зависимости от площади и наличия парковочного места:

Если применить простую линейную регрессию ( решается методом наименьших квадратов):

lm(price ~ parking:area_interior)

Получается следующий результат:

## Coefficients:
##                            Estimate Std. Error t value Pr(>|t|)    
## (Intercept)                33776.10   22175.97   1.523    0.129    
## parkingFALSE:area_interior   444.28      23.54  18.876   <2e-16 ***
## parkingTRUE:area_interior    523.01      19.65  26.614   <2e-16 ***

Т.е. в районе моего дома каждый квадратный фут в доме без парковки добавляет 444$ к базовой цене в 33к, а с парковкой +523$.

Если подставить мои параметры, получаем ожидаемую стоимость 443k$, с доверительным интервалом [433k$ — 453k$]

Однако, если посмотреть на разницу между ценой, которую предсказывает модель и тем что есть на рынке получается такая картина:

Т.е. имеет ошибка зависит от предсказанной величины, т.е. мы имеем дело с нарушением одного из условий использования простой линейном модели. В принципе, на это дело можно забить и использовать и так. Но, если изучить литературу, можно обнаружить, что для предсказания цен с иногда используют generalized linear model с распределением ошибки inverse Gaussian distribution и логарифмической функцией связи, по-русски звучит коряво, но стоит попробовать:

И оценка стоимости: 435k$, с 95% диапазоном [419k$ — 450k$] — видно, что с практической точки зрения оценка осталась такой же.

Заметим, что в обоих моделях я игнорирую количество комнат, этаж и расположение, зато результаты модели очень просто воспринимать — всего три параметра.
Если хочется учитывать больше параметров в модели, возникает проблема, что многие из них на самом деле коррелируют друг с другом — т.е. в квартире с большей площадью обычно больше комнат и больше шансов, что будет парковка, да и если просто комбинировать большое количество дискретных факторов (комнаты X парковка X этаж ) возникает слишком много параметров.

Поэтому, можно модель упростить и пытаться вывести самый главный параметр (площадь) за скобки, и оценивать цену за квадратный фут (или метр, но в попугаях я длиннее).

Если опять использовать generalized linear model с распределением ошибки inverse Gaussian distribution и логарифмической функцией связи:

glm(price_sqft ~ parking + bedrooms,family=inverse.gaussian(link="log")

То получается следующая решение:

## (Intercept)    parkingTRUE   bedrooms2   bedrooms3   bedrooms4 
## 503.1981961   1.1215828   0.9720589   0.9662187   0.8325715

Т.е. квадратный фут стоит, в среднем 503$, парковка добавляет 12% процентов к цене, две спальни — уменьшают цену на 2.8%, 3 спальни — на 3.3%, 4 спальни на 17%, по сравнению с квартирой на одну спальню той же площади.
Итоговая оценка моей квартиры получается 430k$ [ 413k$ — 448k$]

Временная динамика

Все предыдущие оценки не принимали во внимание как цена меняется во времени. Было бы интересно отслеживать как как цена меняется в зависмости от фазы луны со временем.
Если просто изобразить цены квартир по времени и пытаться что-то сгладить, получается какой-то тренд, но не понятно — это из-за того что цены растут, или просто в этом году стало больше квартир с парковками продаваться?
Для начала, сглаживаем цены по всем квартирам функцией loess.

Если просто разбить квартиры по группам, получается каша — потому-что в каждой группе слишком данные слишком разреженные.

Понятно, что квартир не так много чтобы простым сглаживанием можно было вычленнить временную составляющую. Нужно, чтобы можно было одновременно учитывать постоянную разницу между квартирами с разными параметрами ( а желательно ещё и положением) и общим временным трендом.

Тут нам на помощь приходит моделирование с помощью “обобщенных аддитивных моделей” Generalized additive model

Если не вдаваться особенно в детали, то это позволяет взять обобщенную линейную модель и добавить к ней некую гладкую функцию известных переменных для улучшения предсказательных свойств. В R это реализовано в пакете mgcv через функцию gam:

gam(price_sqft ~ parking + bedrooms + s(start_date, k=24), family=inverse.gaussian(link="log"))

Таким образом мы строим модель, опять используя логарифмическое преобразование, с распределением ошибки inverse Gaussian distribution, но с дополнительной гладкой функцией, которая зависит от времени, в этом случае будут сплайн с 24 узлами. Интересная особенность метода gam — автоматическая оценка гладкости данных, в результате параметр k не очень сильно влияет на конечный результат.

В результате получается такая аппроксимацию данных ( для 2х комнатной квартиры с парковкой):

В результате, получается такая оценка цены: 429k [413k-447k], т.е на самом деле тоже самое что и раньше. Но теперь можно надеяться что цена растёт, а не падает со временем.

Сколько времени уйдёт на продажу

Тут всё довольно просто, если использовать метод для анализа выживаемости В принципе по этому методу можно выяснить что разные виды квартир продаются с разной скоростью, но разница не существенная.

Видно, что половина квартир исчезает с рынка за 60 дней. На самом деле это означает что они продаются ещё быстрее, просто их не сразу снимают с сайта.

Оцениваем дома

Аналогично, когда я собирался смотреть чего там мне хотят продать риэлторы, я для начала смотрел сколько стоят дома в округе. Легко сделать выборку по списку всей недвижимости, чтобы найти чего продавалось в радиусе 1км, например для триплекса которым я как-то интересовался:

# индекс недвижимости на сайте листингов
selected_mls=17758383 
# максимальная дистанция поиска 2км
max_distance=2000  
# не смотрим на односемейные дома и квартиры
plex_pe<-prop_geo_p %>% filter(type!='Apartment', type!='House')
ref<-plex_pe%>%filter(mls==selected_mls) 

# делаем круглую маску вокруг дома
search_roi <- st_buffer(ref, max_distance) 
# и отфильтровываем всякие левые записи , где риелтор внёс какую-то чушь 
result <- st_intersection(plex_pe %>% filter(mls!=selected_mls), search_roi) %>% 
filter(area_interior<10000, area_interior>100,area_land>0,price<1e7,price>100 ) 

Получаем результат:

И смотрим на цены:

Видно, что за дом просят слегка больше чем можно было бы рассчитывать, но видно что вариабельность цен гораздо больше чем у квартир, что в общем-то понятно — значение имеет не только площадь и количество спален, но и гораздо больше состояние ( все эти плексы были построены в основном в начале XX века), какие там жильцы и сколько они платят и т.д.
Если использовать модель через оценку стоимости квадратного фута, то получается оценка 523k$, с диапазоном [ 570k$ — 620k$]

Рисуем карты со средними ценами

Когда я только начинал заниматься всеми этими упражнениями, я хотел нарисовать карту распределения цен в Монреале. Для начала, можно просто усреднить цены на аренду по каждому кварталу. С помощь функций из пакета sf это делается элементарно:

Выбираем квартиры с двумя спальнями, берём только одну интересующую нас переменную (цена), делаем пространственное объединение с картой кварталов, а потом считаем медиану по каждому кварталу:

aggregate(filter(kijiji_geo_p,bedrooms==2)%>%dplyr::select(price), mtl_p, median, join = st_contains)

Получилась карта с прямыми линиями по границам кварталов, что немного противоестественно. Вряд ли с одной цены улицы цена аренды будет так сильно отличатся от другой. Нужно применить метод который позволяет моделировать цену плавно.
Применяем обобщенную линейную модель, только теперь плавная функция будет зависеть от пространственных координат:

gam(price_sqft ~ type + bedrooms + parking + s(x,y,k=100), family=inverse.gaussian(link="log"))

Интерполируем результаты на растр, с шагом 100м:

pred_rent_whole <- raster(extent(mtl_land),res=100)
crs(pred_rent_whole)<-crs(mtl_land)
my_predict<-function(...) predict(...,type="response")
pred_rent_whole<- raster::interpolate(pred_rent_whole, model_rent_geo_whole, fun=my_predict, xyOnly=T,const=data.frame(bedrooms=2))

# убираем информацию экстраполированную за пределы острова
pred_rent_whole <- mask(pred_rent_whole, mtl_land)

И рисуем это всё с помощью tmap:

tm_shape(osm_mtl)+tm_rgb(alpha=0.6)+
  tm_shape(mtl_arr) + tm_borders(alpha=0.8, col='black')+
  tm_shape(pred_rent_whole)+tm_raster(style="cont",alpha=0.7, title='$')+  tm_shape(subway_stop_p%>%dplyr::select(stop_name))+tm_symbols(col='blue',alpha=0.2,size=0.03)+
  tm_shape(subway_p)+tm_lines(col='blue',alpha=0.2)+
  tm_compass(position=c("right", "bottom"))+
  tm_scale_bar(position=c("left", "bottom"))+
  tm_layout(scale=1.5)

Синенькие линии — местное метро.

Аренда в центральной части, которая мне более интересна.

Аналогично, можно построить распределение цен на дома.

Ну и как всякий хороший капиталист, я могу теперь оценить прибыльность покупки триплекса ( цена аренды за год/(цена за квадратный фут * площадь).
Строим распределение средняя площади триплексов с тремя спальнами и парковкой.

И рассчитываем прибыльность ( аренда за год/общая цена ).

Ну и заодно, чтобы два раза не вставать, можно отслеживать как изменяется цена со временем (я выбрал только районы которые мне были интересны).

Заключение

Вот так, с помощью пары скриптов на питоне и R можно нарисовать пару красивых графиков оценить что, когда и где покупать или продавать. Но жизнь штука более сложная, в реальном применении не хватает знания реальной цены продажи (в наших краях это доступно только зарегистрированным риэлторам). Так что не стоит ожидать, что полученные прогнозы совпадут с действительностью на 100%. В общем кто не спрятался — я не виноват.

Исходный код

Все данные и исходный код выложен в репозитории. Покупайте наших слонов!

Бонус для тех, кто дочитал до конца

Интерактивная карта с результатами: http://www.ilmarin.info/re_mtl/