• Home

비선형 모형과 추세 분석(아파트 매매 데이터 기반)

지난번 포스팅을 통해서 얻은 데이터와 그 이전 포스팅의 아파트 매매가에 미치는 층수, 크기, 년도 효과에 대한 분석의 후속 분석으로 같은 데이터를 기반을 하는 분석이지만 2015년 Q1의 데이터가 어느정도 모였으니 이의 가격동향과 더불어 다중 회귀모형의 비선형적인 효과를 좀더 다른 방식의 알고리즘으로 모델링 해보고 시각화 해보는 과정을 거쳐보도록 하겠다.

참고로, 분석을 위해 2015년 2월까지의 아파트 매매 데이터를 다시 크롤링했다.

노파심에서 이야기 하지만, 이 분석은 부동산에 전혀 지식이 없는 한 개인이 내가 살고 있는 아파트 단지와 근처 몇몇 단지를 분석한 결과이므로 확대 해석은 매우 위험한 생각이라 말씀드리고 싶다. 그저 생활 데이터 분석의 하나일 뿐이다.

suppressPackageStartupMessages(
{
library(data.table)
library(mgcv)
library(party)
library(randomForest)
library(ggplot2)
library(dplyr)
library(extrafont)
library(randtests)
})

#http://freesearch.pe.kr/archives/4298 에서 만든 데이터를 로딩한다. 
load('result_sales_dt.RData')

setnames(result_sales_dt, 1:10,c('si_gun_gu', 'm_bun', 's_bun', 'dangi', 'area','cont_date', 'price', 'floor', 'year_of_construct', 'road_nm'))

result_sales_dt[,price:=as.numeric(price)]
result_sales_dt[,floor:=as.numeric(floor)]
result_sales_dt[,mm:=substr(yyyymm, 5,6)]
result_sales_dt[between(as.numeric(mm), 1, 3), qrt:='Q1']
result_sales_dt[between(as.numeric(mm), 4, 6), qrt:='Q2']
result_sales_dt[between(as.numeric(mm), 7, 9), qrt:='Q3']
result_sales_dt[between(as.numeric(mm), 10, 12), qrt:='Q4']
result_sales_dt[,yyyyqrt:=paste0(substr(yyyymm, 1,4), qrt)]
result_sales_dt[,yyyy:=factor(substr(yyyymm, 1,4))]
#이전 포스팅과 다르게 쿼터로 날짜는 자른다. 
result_sales_dt[,yyyyqrt:=factor(yyyyqrt)]


my_apt <- result_sales_dt[si_gun_gu %chin% " 서울특별시 관악구 신림동" & m_bun == 1730]

ggplot(my_apt, aes(floor, price)) + geom_point(aes(colour=yyyy)) 

plot of chunk unnamed-chunk-1

일단 층별 매매가를 보면 위와 같은 비선형적인 패턴이 보인다. 익히 알고 있는 로열층이 가격대가 높은것을 실제 데이터를 통해 확인이 가능하다.

이런 사전지식 하에서 이전 포스팅에서 했던 층변수에 2차 함수를 할당한 회귀모형 시각화를 확인해 보자.

#아웃라이어 레코드 제거 (잔차 분석 참고 ...)
my_apt <- my_apt[-c(54,67,202)]

lm_mdl <- lm(price ~ yyyyqrt + I(floor^2) + floor  + area, data=my_apt)

#지난번 하지 않았던 오차상관성 검정... 상관성은 없다고 봐도 무방하다. 
runs.test(resid(lm_mdl))
## 
##  Runs Test
## 
## data:  resid(lm_mdl)
## statistic = 1.0647, runs = 226, n1 = 214, n2 = 214, n = 428,
## p-value = 0.287
## alternative hypothesis: nonrandomness
levelss <- levels(factor(my_apt$yyyyqrt))

aggrs_fun <- function(floor, area, mdls){
  predict(mdls, data.table(yyyyqrt=factor(rep("2015Q1"),levels=levelss), floor=floor, area=area))
}


floor_seq <- seq(1,24, length=30)
area_seq <- seq(59,140, length=30)

persp(floor_seq,area_seq, z=outer(floor_seq,area_seq, aggrs_fun, lm_mdl), col="green", ticktype="detailed",shade=TRUE, theta=-25,phi = 5, xlab="floor", ylab="area", zlab = "price")

plot of chunk _lm

나쁘지 않아 보인다. 사실 다중 회귀 모형에서 2차 함수 형태를 직접 데이터 시각화를 통해서 발견해 채용한 그 과정이 흥미로웠다고 생각한다.

regression tree를 동일한 방식으로 살펴보자!

ctree_mdl <- ctree(price ~ yyyyqrt + floor + area, data=my_apt)

persp(floor_seq,area_seq, z=outer(floor_seq,area_seq, aggrs_fun, ctree_mdl), col="green", ticktype="detailed",shade=TRUE, theta=-25,phi = 5, xlab="floor", ylab="area", zlab = "price")

plot of chunk tree

물론 예상했지만 거의 비선형적인 패턴을 추출하지 못한것으로 보인다. 왜냐면 6층부터 24층까지 동일한 층수 효과를 적용했기 때문이다.

rf_mdl <- randomForest(price ~ yyyyqrt + floor  + area, data=my_apt)


persp(floor_seq,area_seq, z=outer(floor_seq,area_seq, aggrs_fun, rf_mdl), col="green", ticktype="detailed",shade=TRUE, theta=-25,phi = 5, xlab="floor", ylab="area", zlab = "price")

plot of chunk rf

랜덤포레스트는 일반 regression tree보다는 비선형적인 특징이 잘 설명되는것을 볼 수 있으나 사실 눈으로 확인한 패턴에 가깝지는 않다.

gam_mdl <- gam(price ~ yyyyqrt + s(floor) + area, data=my_apt)


persp(floor_seq,area_seq, z=outer(floor_seq,area_seq, aggrs_fun, gam_mdl), col="green", ticktype="detailed",shade=TRUE, theta=-25,phi = 5, xlab="floor", ylab="area", zlab = "price")

plot of chunk gam

개인적으로는 GAM(generalized additive model)에서 Spline형태로 표현된 비선형 패턴이 가장 실제 데이터에서 보였던 패턴과 가장 가깝다고 생각한다.

그럼 몇층이 가장 가치가 높다고 볼 수 있을까?

plot(gam_mdl)

plot of chunk unnamed-chunk-2

8층에서 20층 사이는 오차 범위가 겹치니 이중에 본인이 좋아하는 층을 선택하는게 좋을것 같다.


2015년 1Q가 다 지나가고 데이터도 좀 쌓였으니 시간에 따른 매매 가격 변동 추이를 자세하게 볼 필요가 있을 것 같다.

theme_set(theme_bw(base_family = "UnBatang"))

ptbl <- data.frame(summary(gam_mdl)$p.table)
ptbl$params <- rownames(ptbl)
ptbl <- data.table(ptbl)  

ptbl_yyyy <- ptbl[params %like% "^yyyy"]

ggplot(ptbl_yyyy, aes(x=params,y=Estimate)) + geom_errorbar(aes(ymin=Estimate - 2*`Std..Error`, ymax=Estimate + 2*`Std..Error`), linetype=2) + geom_point() + geom_line(aes(group=1)) + theme(axis.text.x=element_text(angle=90)) + geom_smooth(aes(group=1), method=loess) + ggtitle("서울특별시 관악구 신림동 1730 신림 푸르지오 1차 ")

plot of chunk unnamed-chunk-3

위 그래프는 해당 단지에서 면적과 층수에 대한 효과를 제외한 시간에 따른 가격 변동을 2006년을 기준(분양 이후)으로 표현한 도식화된 그래프이다. 과연 금년 3월에 실시된 역대 최저 주택담보대출 금리의 효과가 앞으로 어떻게 반영이 될지 매우 궁금하다.

비교를 위해 근처 두 아파트 단지의 추이를 같은 방법으로 비교했다. 2007년, 2008년이 매매가가 고점이였다는 걸 감안하고 보면 된다.

my_apt2 <- result_sales_dt[si_gun_gu %chin% " 서울특별시 관악구 신림동" & m_bun == 1736]


gam_mdl2 <- gam(price ~ yyyyqrt + s(floor) + area, data=my_apt2)


ptbl <- data.frame(summary(gam_mdl2)$p.table)
ptbl$params <- rownames(ptbl)
ptbl <- data.table(ptbl)  

ptbl_yyyy <- ptbl[params %like% "^yyyy"]

ggplot(ptbl_yyyy, aes(x=params,y=Estimate)) + geom_errorbar(aes(ymin=Estimate - 2*`Std..Error`, ymax=Estimate + 2*`Std..Error`), linetype=2) + geom_point() + geom_line(aes(group=1)) + theme(axis.text.x=element_text(angle=90)) + geom_smooth(aes(group=1), method=loess) + ggtitle("서울특별시 관악구 신림동 1736 신림 푸르지오 2차")

plot of chunk unnamed-chunk-4

my_apt3 <- result_sales_dt[si_gun_gu %chin% " 서울특별시 관악구 신림동" & m_bun == 1735]


gam_mdl3 <- gam(price ~ yyyyqrt + s(floor) + area, data=my_apt3)


ptbl <- data.frame(summary(gam_mdl3)$p.table)
ptbl$params <- rownames(ptbl)
ptbl <- data.table(ptbl)  

ptbl_yyyy <- ptbl[params %like% "^yyyy"]

ggplot(ptbl_yyyy, aes(x=params,y=Estimate)) + geom_errorbar(aes(ymin=Estimate - 2*`Std..Error`, ymax=Estimate + 2*`Std..Error`), linetype=2) + geom_point() + geom_line(aes(group=1)) + theme(axis.text.x=element_text(angle=90)) + geom_smooth(aes(group=1), method=loess) + ggtitle("서울특별시 관악구 신림동 1735 관악산 휴먼시아2차")

plot of chunk unnamed-chunk-4

세단지 모두 최근의 추세가 반등을 하느냐, 못하느냐 가장 중요한 시점임에는 분명한것 같다. 정부의 정책과 더불어 앞으로의 매매 추이가 매우 궁금해진다.

국토교통부 실거래가 데이터 크롤링 코드

국토교통부 실거래가 데이터 스크래핑 코드를 공유한다. 블로그에 썼던 데이터 분석에 사용한 데이터는 친분이 있는 분으로 부터 받은 데이터인데, 새로운 매매 데이터가 올라가면서 매번 요청하기 힘들것 같아서 코드를 작성해 봤고, 아래와 같이 간단하게 스크래핑 코드를 만들 수 있었다.

데이터 스크래핑 코드는 항상 그렇듯이 임시방편적이고, 코드가 깨끗하지 않다. 그래서 좀 색다르게 magrittr 방식으로 코딩해 봤으나 그다지 나아 보이지 않는다. ;;

수백만건의 데이터를 아래와 같은 간단한 코드로 수집할 수 있다는게 데이터 크롤링의 재미가 아닐까 생각하며, 다시한번 rvest, stringi, XLConnect, data.table과 같은 유용한 패키지의 존재에 감사한다.

library(rvest)
library(httr)
library(stringi)
library(XLConnect)
library(data.table)

## 엑셀 데이터 다운로딩 부분
# 데이터는 working dir아래 data 디렉토리에 쌓인다. 따리서 data 디렉토리를 만들어 주자!

su <- file("succ.txt", "w")

agent_nm <- 
  "Mozilla/5.0 (Macintosh; Intel Mac OS X 10.10; rv:35.0) Gecko/20100101 Firefox/35.0"

maxidx <- html('http://rt.molit.go.kr/rtFile.do?cmd=list') %>% 
  html_nodes('.notiWhite1 .notiBorad01') %>%
  html_text %>% as.numeric %>% max

for(i in maxidx:1){
  urls_view <- sprintf("http://rt.molit.go.kr/rtFile.do?cmd=view&seqNo=%d", i)
  r <- GET(urls_view,
           user_agent(agent_nm))
  htxt <- html(r, "text")

  if((html_nodes(htxt, "td.notiBorad14")[[2]]  %>% 
        html_text()  %>% 
        stri_trim_both) == '첨부파일이 존재하지 않습니다.') next

  download_tags <- html_nodes(htxt, "td.notiBorad14")[[2]] %>%
    html_nodes('a[href^="javascript:jsDown"]')  

  for(dtag in download_tags){

    dtag %>% html_attr("href") %>% 
      stri_match_all_regex(pattern="javascript:jsDown\\('([0-9]+)','([0-9]+)'\\);") %>%
      .[[1]]  %>%  
          { 
            f_idx <<- .[2] %>% as.numeric
            s_idx <<- .[3] %>% as.numeric
          }

    f_nm <- dtag %>% html_text

    urls <- sprintf("http://rt.molit.go.kr/rtFile.do?cmd=fileDownload&seq_no=%d&file_seq_no=%d", f_idx,s_idx)
    r <- GET(urls,
             user_agent(agent_nm))
    bin <- content(r, "raw")

    #1kb 미만의 데이터는 버림(에러 페이지?) 
    if(length(bin) < 1000) next
    writeBin(bin, sprintf("data/%s",f_nm))
    cat(sprintf("%d, %d\n", f_idx,s_idx), file = su)
    print(sprintf("%d, %d", f_idx,s_idx))
  }
}

close(su)




## 엑셀 데이터에서 테이블을 추출해 하나의 아파트 매매 데이터로 통합하는 코드 
## 연립 다세대, 단독 다가 데이터도 간단한 코드 변환으로 통합할 수 있다. 


f_list <- list.files('data') %>%  
  stri_trans_nfc  %>% 
  .[stri_detect_fixed(.,'매매아파트')]

total_list <- list()

#파일의 각 시트(지역)를 data.table로 변환하고 필요 필드 추가함 
for(xlsf in f_list){
  wb <- loadWorkbook(paste0('data/',xlsf))
  sells <- list()
  fname <- stri_replace_last_fixed(xlsf, '.xls',"")
  yyyymm <- substring(fname, 1, 6)
  typenm <- substring(fname, 9)
  for(nm in getSheets(wb)) {
    df <- data.table(readWorksheet(wb, sheet = nm, header = TRUE))
    df[,`:=`(region=nm, yyyymm=yyyymm, typenm=typenm)]
    df[,`거래금액.만원.`:= stri_replace_all_fixed(`거래금액.만원.`, ',', '')]
    sells[[nm]] <- df
  }
  total_list[[paste0(yyyymm,typenm)]] <- rbindlist(sells)
}

result_sales_dt <- rbindlist(total_list)

#결과 데이터 저장 
save(result_sales_dt, file='result_sales_dt.RData')

실거래가 페이지가 바뀌지 않는 이상 동작을 할 것이고, result_sales_dt객체에 2015년 1월 기준 4,964,004건의 매매데이터가 쌓여 있는 것을 확인해 볼 수 있을 것이다.

사족이나 그동안 스크래핑은 Python과 같은 다른 언어를 사용하거나 webkit 라이브러리를 이용해 브라우저 객체를 통해 스크래핑 하곤 했는데, 이젠 여러 패키지 덕분에 R에서도 매우 간단해졌다. 10년 전만해도 HTTP 프로토콜 책을 옆에 끼고 크롤러 코딩을 했었는데, 격세지감을 느낄 따름이다.

윈도우에서 동작여부가 매우 궁금하나, 이런 부분에 대해서는 case by case라 시간 소모하고 싶지 않다.

그러나 참고를 위해서 코드가 동작하는 환경 정보는 추가한다.

## R version 3.1.2 (2014-10-31)
## Platform: x86_64-apple-darwin13.4.0 (64-bit)
## 
## locale:
## [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
## 
## attached base packages:
## [1] stats     graphics  grDevices utils     datasets  methods   base     
## 
## other attached packages:
## [1] RWordPress_0.2-3    knitr_1.9           XLConnect_0.2-10   
## [4] XLConnectJars_0.2-9 stringi_0.4-1       rvest_0.2.0        
## [7] data.table_1.9.4   
## 
## loaded via a namespace (and not attached):
##  [1] bitops_1.0-6    chron_2.3-45    digest_0.6.8    evaluate_0.5.5 
##  [5] formatR_1.0     htmltools_0.2.6 httr_0.6.1      magrittr_1.5   
##  [9] markdown_0.7.4  plyr_1.8.1      Rcpp_0.11.4     RCurl_1.95-4.5 
## [13] reshape2_1.4.1  rJava_0.9-6     rmarkdown_0.5.1 stringr_0.6.2  
## [17] tools_3.1.2     XML_3.98-1.1    XMLRPC_0.3-0    yaml_2.1.13

내가 사는 아파트 매매 가격 분석

아직 무주택자인 관계로 국내 부동산 매매 시장에 매우 관심이 많다. 그러다 몇번의 국토교통부 실거래가를 수집하려 틈틈히 노력을 기울였으나 대부분 전수 데이터를 얻는덴 실패하고 그러다 잊고 지내기를 몇번을 거쳐왔는데 최근 실거래가 조회 페이지에 자료제공항목이 생기면서 이런 저런 과정을 통해 원천 전수 데이터를 얻을 수 있게 되었다.

데이터를 얻어서 할 수 있는 일은 역시 분석을 통한 정보의 습득이다. 그리고 그 분석이 내가 지금 살고 있는 아파트와 관련이 있으면 또한 재미가 있을 것인데 이유는 지금 전세를 살고 있는 아파트의 환경이 꽤 맘에 들기 때문이다.

부동산 데이터는 멀티레벨 분석에 최적화된 데이터 특징을 가지고 있다. 이는 지역별 고정효과가 있고 지역의 아파트 브랜드에 따른 랜덤효과가 분명 있을 것이기 때문이다. 더 나아가서는 층별효과 등이 존재할 것이다.

궁극적으로 하고 싶은 분석은 멀티레벨 분석이나 멀티레벨의 시작점은 역시 회귀모형이니 작은 단위로 모델링이 가능한지 확인하고 간단한 추론 과정도 거쳐보도록 하겠다.

2006년부터 2014년 12월까지 국토교통부에 등록된 전국 아파트 매매 데이터는 약 500만건이다. 이 중에서 현재 살고 있는 신림푸르지오1차 아파트의 전체 매매 데이터 423건을 가져와 분석을 시작해 보겠다.

suppressPackageStartupMessages({
  library(sjPlot)
  library(dplyr)
})

load('myapt.RData')

#대략 어떤 포맷인지?
glimpse(myaptdf)
## Observations: 423
## Variables:
## $ id          (dbl) 1693, 43839, 43840, 43841, 105791, 105792, 105793,...
## $ yyyymm      (chr) "200601", "200602", "200602", "200602", "200603", ...
## $ region      (chr) "서울", "서울", "서울", "서울", "서울", "서울", "서울", "서울", "서...
## $ district    (chr) "서울특별시 관악구 신림동", "서울특별시 관악구 신림동", "서울특별시 관악구 신림동",...
## $ street_main (chr) "1730", "1730", "1730", "1730", "1730", "1730", "1...
## $ street_sub  (chr) "0", "0", "0", "0", "0", "0", "0", "0", "0", "0", ...
## $ name        (chr) "신림푸르지오", "신림푸르지오", "신림푸르지오", "신림푸르지오", "신림푸르지오", ...
## $ dealtype    (chr) "BUY", "BUY", "BUY", "BUY", "BUY", "BUY", "BUY", "...
## $ renttype    (chr) "NA", "NA", "NA", "NA", "NA", "NA", "NA", "NA", "N...
## $ area_excl   (dbl) 138.74, 84.79, 114.84, 114.84, 84.79, 84.79, 84.79...
## $ date        (chr) "11~20", "21~28", "1~10", "21~28", "1~10", "11~20"...
## $ property    (dbl) 63500, 41200, 49000, 49700, 44800, 39000, 40500, 4...
## $ deposit     (dbl) 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,...
## $ fee_month   (dbl) 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,...
## $ floor       (dbl) 22, 16, 6, 14, 11, 24, 16, 5, 10, 9, 19, 1, 5, 14,...
## $ build_year  (chr) "2005", "2005", "2005", "2005", "2005", "2005", "2...
## $ roadname    (chr) "남부순환로", "남부순환로", "남부순환로", "남부순환로", "남부순환로", "남부순환...
## $ yyyy        (chr) "2006", "2006", "2006", "2006", "2006", "2006", "2...

myaptdf 데이터의 yyyy 필드는 년도별 효과를 보기위해 필자가 생성한 변수이며, 원본의 yyyymm에서 추출한 필드이다.

모델 설명변수로는 층수, 매매년도, 면적을 넣는다.

fit <- lm(property ~ floor + yyyy + area_excl, data=myaptdf)

summary(fit)$adj.r.squared
## [1] 0.7752066
plot(fit, 1)

plot of chunk _mdl

plot(fit, 2)

plot of chunk _mdl

shapiro.test(resid(fit))
## 
##  Shapiro-Wilk normality test
## 
## data:  resid(fit)
## W = 0.974, p-value = 7.509e-07

오차의 등분산성과 오차의 정규성 그래프를 확인해보면 세 레코드의 아웃라이어 데이터가 존재하는 것을 볼 수 있다. 원인이 뭔지 모르겠지만 예측값에 비해 매우 높은 가격을 주고 매매가 이뤄진 경우라 볼 수 있는데, 이 데이터는 제거하고 다시 모델링을 해본다(가장 고점이였던 2007년 2008년을 고려해 보더라도 이 세 매매건은 말도 안되는 가격에 거래가 된 것이라 볼 수 있다. 예측 가격 대비 적게는 1억 많게는 2.5억 더 주고 거래가 된 케이스니..). 게다가 모델은 오차 정규성 검정도 실패를 했다는 것을 알 수 있다.

fit2 <- lm(property ~ floor + yyyy + area_excl, data=myaptdf[-c(54,67,69),])

summary(fit2)$adj.r.squared
## [1] 0.7923766
AIC(fit2)
## [1] 8191.741
plot(fit2, 1)