房价影响因素的定性探索
经济领域的变量相互依赖相互影响，要找出房价的影响因素必须要有大量的数据支持。从以下几个大方向寻找影响房价的因素：人口、宏观经济、房地产投资规模、政策和利率、区域差异。本章从大量的数据出发，定性地探索上述变量与房价的关系，并分析相关程度有多大。
商品住宅销售价格波动
数据集data3，data2（数据来源：中宏产业库）中存的是分别全国及各地区商品住宅销售额、商品住宅销售面积，商品住宅投资额、施工面积、竣工面积以及土地购置费用的月度数据。它们的时间跨度为2005年2月至2018年10月。用商品住宅销售额除以商品住宅销售面积并转换单位得出相应时间的商品住宅房价（元/平方米），并将全国平均房价保存在数据集house_price中。这样计算的房价与实际房价有一定的差异，如果商品住宅销售额和商品住宅销售面积的数据可靠的话，房价数据也是可靠的。纵向地，我们看全国平均房价的时间变化。
data3=data3%>%mutate(price=housesale/area*10000)
house_price=xts(data3$price,as.Date(data3$date,format='%Y/%m/%d'))
plot(as.xts(house_price))
可见全国平均房价整体上是上升的，可以看到每年初房价都有一个下降，房价可能存在季节效应。其中，2009年至2010年保持较快的增长速度；2013年至2016年房价基本持平，大约保持在6100元/平方米。2016年后全国平均房价保持持续增长。
人口因素
衣食住行是人的基本生活要素，住房是人的必需品，人口总量和城乡结构对房价必然有很大的影响。城市人口的自然出生率将影响未来几十年房价的需求，外来人口在城市的迁入迁出则对房屋当下的需求（主要是消费型需求）有影响。人口总量是影响房屋消费的双刃剑，一方面人口增加加剧了资源稀缺性，使资源人均占有率下降，减缓区域经济增长，人均收入下降，房屋购买力下降。另一方面，人口增长伴随着城镇化进程，城镇人口增加大大激发了房屋需求，从而拉动房地产及房价的上升。此外，人口的年龄结构对房地产业有一定影响，20～30岁年龄层的青年面临结婚成家，他们购房主要是结婚所需，中年层的消费者生活基本稳定，没有婚姻的压力，他们购房或者是寻求更加舒适的居住条件或者是将购房作为一种投资方式。数据集data1（数据来源：国家统计局）保存了全国各地区2005至2017年度人口总量、城镇人口数量、人均GDP、居民消费价格指数和城镇居民人均可支配收入。
distribution=data1%>%fill(region)%>%
group_by(year)%>%
mutate(pop=sum(population))%>%
mutate(urban=sum(urbanpop))%>%
mutate(rural=pop-urban)%>%
filter(region=="beijing")%>%
gather("urban","rural",key="region",value = "popul")%>%
ggplot()+
geom_bar(aes(year,weight=popul,fill=region),position = "dodge")
distribution
由上图可以看出，2005-2017年间，全国城镇人口从5.67亿人增加到8.14亿人，平均增长率为3.42%；而农村人口逐年下降，并且在2010年全国城镇人口超过农村人口。城镇化水平也从2005年的0.43持续增长至2017年的0.59。
data1中的数据是年度数据，取data2中每年12月的数据作为该年年度数据。分别比较各地区房价与人口总量、房价与城镇化水平的关系。
data1=data1%>%fill(region)%>%
mutate(urblevel=urbanpop/population)%>%
mutate(rank=c(1:403))
data2=data2%>%fill(region,year)%>%
mutate(price=housesale/area*10000)
price=data2%>%filter(month==12)%>%
select(price)%>%
mutate(rank=c(1:403))%>%
left_join(data1,by="rank")%>%
select(region,year,population,urblevel,price)%>%
group_by(region)%>%
nest()
pop=price%>%
mutate(model=map(data,~lm(price~population,data=.)))%>%
mutate(resids=map2(data,model,add_residuals))%>%
unnest(resids)
ggplot(pop,aes(year,resid))+
geom_line(aes(group=region),alpha=1/3)+
geom_smooth(se=F)
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
urb_level=price%>%
mutate(model=map(data,~lm(price~urblevel,data=.)))%>%
mutate(resids=map2(data,model,add_residuals))%>%
unnest(resids)
ggplot(urb_level,aes(year,resid))+
geom_line(aes(group=region),alpha=1/3)+
geom_smooth(se=F)
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
从以上两图可以看出，人口、城镇化水平与房价的均有几条波动很大的残差曲线，这说明房价可能存在着区域差异。由于城镇化水平的数据是由人口数量计算而得，为避免多重共线性，选其一作为房价影响因素的解释变量。又除了波动较大的几个地区之外，城镇化水平与房价的残差图形更紧密的围绕0波动，故选择城镇化水平作为解释变量。
宏观经济因素
房地产是具有投资与消费双重属性的特殊商品。在国内宏观经济的大环境中，房地产的发展必然受到宏观经济因素的制约或促进。经济发展水平、工业化水平、居民收入水平、消费价格指数等等都是宏观经济状态的指标。理论上认为，经济发展状况良好时，投资和消费高涨，必然存在通货膨胀，物价上涨，房价自然不例外。供求关系是相互制约相互影响的，房地产投资增加必然导致房地产供给增加，给房地产业带来一个阻力，使之发展放缓。但是房价与投资、物价、收入水平的关系仍然复杂。
本文从以下几个方面探索宏观经济发展状况：地区生产总值（GDP）、居民消费价格指数（CPI）、工业化水平。
经济增长、居民收入与支出
数据集data4（数据来源：国家统计局）保存了从2005年到2017年人均国内生产总值，工业化水平和居民消费指数的年度数据。取data3中每年12月的房价数据作为当年房价数据。
2005到2017年间，人均GDP、工业化水平逐年递增。人均GDP从2005年的14368元增加到2017年的59660元，增长了4.15倍。工业化水平是第二第三产业总产值除以国内生产总值得到的数值，它是国家工业化程度的度量，它从2005年的0.88增加到2017年的0.92。
data=data3%>%
filter(month(ymd(date))==12)%>%
mutate(year=year(ymd(date)))%>%
select(price,year)%>%
right_join(data4,by="year")
ggplot(data,aes(year,price))+
geom_point(aes(y=pgdp))
ggplot(data,aes(year,price))+
geom_point(aes(y=indlevel))
ggplot(data,aes(year,price))+
geom_point(aes(y=cpi))
观察人均GDP、工业化水平、CPI对房价的散点图。他们中与房价可能有相关性的是工业化水平、人均GDP。再比较他们与房价的相关系数矩阵。
co=cor(data[complete.cases(data),-2])
knitr::kable(co,caption="相关系数" ,align="c",digits = 2)
相关系数
| price |
1.00 |
0.99 |
0.95 |
-0.26 |
| pgdp |
0.99 |
1.00 |
0.95 |
-0.19 |
| indlevel |
0.95 |
0.95 |
1.00 |
-0.16 |
| cpi |
-0.26 |
-0.19 |
-0.16 |
1.00 |
可以看出，除了CPI，人均GDP、工业化水平都与房价有明显的正相关关系。
房地产投资规模
本文主要从商品住宅投资额、施工面积、竣工面积等方面来认识房地产业的发展规模。由数据集data2，与房价走势类似，全国商品住宅投资额表现为再波动中稳步上升。并且都在年初处于低位，此后按月增加，到年底达到最大值，次年一月投资额下降到较低水平，此后进行下一轮循环。
inv=data2%>%
filter(!is.na(investment))%>%
group_by(region)
ggplot(inv,aes(region,investment))+
geom_boxplot()+
theme(axis.text.x = element_text(angle = 90))
从图中可以看出，全国商品住宅投资呈现明显的从东部沿海到西部内陆递减的趋势。从2005年2月到2017年6月全国商品住宅投资额最大的五个省份分别是江苏、广东、浙江、山东、辽宁；均值最小的五个省份分别是西藏、青海、宁夏、甘肃、新疆。
由于土地购置费是商品住宅投资额的一部分，土地购置费的时间变化也与商品住宅投资额类似，整体上升，并伴随季节变化。
随机选出一个地区，计算他的商品住宅投资规模对房价的相关系数并作出他们的散点图。随机取出的地区是安徽。
anhui=data2%>%
filter(region=="anhui")%>%
filter(!is.na(investment))%>%
select(year,investment:price)
ah=cor(anhui[complete.cases(anhui),-1])
knitr::kable(ah,caption="相关系数" ,align="c",digits = 2)
相关系数
| investment |
1.00 |
0.94 |
0.81 |
0.95 |
0.68 |
| completed |
0.94 |
1.00 |
0.72 |
0.89 |
0.55 |
| construction |
0.81 |
0.72 |
1.00 |
0.72 |
0.96 |
| land_fee |
0.95 |
0.89 |
0.72 |
1.00 |
0.58 |
| price |
0.68 |
0.55 |
0.96 |
0.58 |
1.00 |
plots=function(x,y){
ggplot(anhui,aes(price))+
geom_point(aes(y=x,color=y),show.legend = F)}
plots(anhui$investment,1)
plots(anhui$completed,2)
plots(anhui$construction,3)
plots(anhui$land_fee,4)
这些结果表明，房价与商品住宅施工面积具有较强的依赖关系，而与其他变量的关系不是很明显。故选取商品住宅施工面积作为解释变量。
宏观经济政策与利率
从2005年以来，政府一直不遗余力地控制房价过高，各项住房政策层出不穷。房地产开发投资层面，政府在土地开发、银行贷款方面施加强有力的控制，消费层面，政府规定第一套住房的首付比例和第二第三套住房的首付比例，以及银行贷款利率。银行利率的调整反应了政府对房价的调控，从国家房产政策可以看出，对房价的宏观调控，利率是必用的手段。
数据集interes记录了从2005年3月到2015年8月半年期贷款利率的变化。从半年期贷款利率图中看出，2007年是国家密集提高利率的年份，半年期贷款利率共五次加息，依次提高0.18，0.18，0.18，0.27，0.09个百分点。到2008年初半年期贷款利率达到了6.57%，2008年4次降息，2011年4次加息，2012年下调利率个0.5百分点。2007年经济形势正在趋热，政府采取紧缩政策，频频提高贷款利率。时隔两年经济形势逆转为趋冷，而此时却迎来一次房价暴增。利率对房价的影响在于，高利率下房产投资开发成本加大，利润空间变小，房产投机者的融资成本增加，持币待购者支付能力下降，多重因素迫使房价下跌。利率是影响房产泡沫的重要因素，当利率上升，预期房产实际价格下降，当降低到一定程度房产泡沫破灭。低利率下的情形相反。利率对房价的作用存在一个滞后效应，2008年利率处于低位，未来两三年内的房价房价高速增长；2011年利率处于高位，未来两三年房价仍然很高但是房价增长速度明显下降。因此认为，利率提高对未来两年内的房价增速有反作用力，利率下降则使房价增速变大。
由于利率是全国数据与房价数据结构不同，故回归模型中不考虑利率因素。
interests=xts(interests$interest,as.Date(interests$date,format='%Y/%m/%d'))
plot(as.xts(interests))
区域差异
前文的分析中已经发现，房价具有明显的地区差异，呈现从东部到西部递减的趋势。用人均GDP、城镇化水平、工业化水平这三个变量来衡量一个地区人口结构和经济发展状况，并通过他们用k-means聚类，将全国分为若干个大的地区，以探讨不同类别的房价差异。
先生成一个从2005年到2017年按年平均的人均GDP、城镇化水平、工业化水平的数据集region，它是一个31行x3列的数据框，每一行是一个地区的人均GDP、城镇化水平、工业化水平。后两者的取值再0-1之间，因此也将作人均GDP作一个极差正规化变换，使之取值再0-1之间。
region=data1%>%
select(region,pgdp,indlevel,urblevel)%>%
group_by(region)%>%
summarise_all(funs(mean))%>%
mutate(pgdp=(pgdp-min(pgdp))/(max(pgdp)-min(pgdp)))
region1=region%>%
select(pgdp:urblevel)
region1=as.data.frame(region1)
row.names(region1)=region$region
dn=dist(region1)
hn=hclust(dn,'ward.D2')
plot(hn)
km=kmeans(region1,5)
sort(km$cluster)
## xizang fujian guangdong jiangsu liaoning
## 1 2 2 2 2
## neimenggu shandong zhejiang chongqing hainan
## 2 2 2 3 3
## hebei heilongjiang hubei hunan jilin
## 3 3 3 3 3
## ningxia qinghai shanxi shanxi2 anhui
## 3 3 3 3 4
## gansu guangxi guizhou henan jiangxi
## 4 4 4 4 4
## sichuan xinjiang yunnan beijing shanghai
## 4 4 4 5 5
## tianjin
## 5
clusplot(region1,km$cluster,color=T,shade=T,labels=2,lines=0)
首先建立31个三维点的距离，采用默认的欧氏距离，画出系谱图。结果显示分成5类比较合适，再使用k-means聚类，从k-means聚类的结果可以看出，综合来说第四类各项指标都很高，故将第四类作为发达地区；将第三类作为较发达地区；第二类作为普通地区；第五类作为不发达地区；第一类作为贫困地区。分类结果如下：
发达地区：北京、上海、天津
较发达地区：福建、广东、江苏、辽宁、内蒙古、山东、浙江
普通地区：重庆、海南、河北、湖南、吉林、宁夏、青海、陕西、陕西
不发达地区：安徽、甘肃、广西、贵州、河南、江西、私闯、新疆、云南
贫困地区：西藏
发达地区为京津沪，较发达地区主要是东部沿海，落后地区为西部地区，总体上看分类效果还不错。在数据中分别用数字1到5标记五个地区。
data1=data1%>%group_by(region)%>%nest()%>%arrange(region)
h=as.data.frame(km$cluster)
h=h%>%mutate(region=data1$region)
h=h[rep(1:nrow(h),each=13),]
data1=data1%>%unnest()%>%
select(region,year,pgdp:urblevel)%>%
mutate(rank=c(1:403))
data1=h%>%mutate(rank=c(1:403))%>%
right_join(data1,"rank")
colnames(data1)[1]="labels"
colnames(data1)[2]="region"
data2=data2%>%arrange(region)%>%
filter(month==12)
regions=data1%>%mutate(price=data2$price)%>%
select(region,year,pgdp:price,labels)%>%
group_by(labels)
av_price=regions%>%select(pgdp:labels)%>%
summarise_all(funs(mean))
knitr::kable(av_price ,caption="平均房价",align="c",digits = 2)
平均房价
| 1 |
21711.46 |
102.87 |
0.87 |
0.02 |
3336.98 |
| 2 |
51291.51 |
102.50 |
0.92 |
0.59 |
5791.29 |
| 3 |
31174.69 |
102.83 |
0.88 |
0.50 |
3967.82 |
| 4 |
23803.96 |
102.80 |
0.85 |
0.40 |
3538.70 |
| 5 |
82227.23 |
102.50 |
0.99 |
0.85 |
13189.59 |
房价区域差异分析
数据集data1中包含各地区从2005年到2017年房价年度数据以及地区分类信息。计算可知，发达地区的房价均值为13189.59元，较发达地区为5791.29元，普通地区为3988.80元，不发达地区为3728.65元，贫困地区为3275.73元。房价均值从发达地区到贫困地区递减。不发达地区和普通地区的数据很接近。由于房价总体分布未知，为比较五类地区的房价是否有差异，采用Willcoxon秩和检验。price1、price2、price3、price4、price5分别是从发达地区到贫困地区的房价数据。只需分别检验零假设：
H01:price1大于price2
H02:price2大于price3
H03:price3大于price4
H04:price4大于price5
price1=filter(regions,price,labels==km$cluster[["beijing"]])$price
price2=filter(regions,price,labels==km$cluster[["fujian"]])$price
price3=filter(regions,price,labels==km$cluster[["chongqing"]])$price
price4=filter(regions,price,labels==km$cluster[["guizhou"]])$price
price5=filter(regions,price,labels==km$cluster[["xizang"]])$price
wilcox.test(price1,price2,alternative = "less")
##
## Wilcoxon rank sum test with continuity correction
##
## data: price1 and price2
## W = 3056, p-value = 1
## alternative hypothesis: true location shift is less than 0
wilcox.test(price2,price3,alternative = "less")
##
## Wilcoxon rank sum test with continuity correction
##
## data: price2 and price3
## W = 9235, p-value = 1
## alternative hypothesis: true location shift is less than 0
wilcox.test(price3,price4,alternative = "less")
##
## Wilcoxon rank sum test with continuity correction
##
## data: price3 and price4
## W = 9131, p-value = 0.8979
## alternative hypothesis: true location shift is less than 0
wilcox.test(price4,price5,alternative = "less")
##
## Wilcoxon rank sum test with continuity correction
##
## data: price4 and price5
## W = 827, p-value = 0.6985
## alternative hypothesis: true location shift is less than 0
从以上结果可以看书，除假设H04=0.70，其他的检验p值都表明强烈接受零假设，这说明出了普通地区和不发达地区外，各地区之间的房价按照从发达地区到贫困地区递减排列，各地区房价存在着明显的差异。
房价回归模型
前文讨论了影响房价的因素以及它们的影响程度，本章用多元回归分析建立房价与这些因素的回归方程，并利用回归方程定量的探讨房价影响因素。回归方程的建立依赖于变量选取、回归诊断，本章的主要内容是如何选取或变换变量建立回归方程，以及如何判断一个回归模型是合适的。
选择变量建立数据集
根据前文，本文选取的解释变量为城镇化水平、工业化水平、人均GDP、施工面积、区域差异。由于区域差异是定性因素，故回归模型中先不考虑区域差异。首先生成回归模型需要使用的数据集test，然后做出房价与各解释变量的散点图，查看是否能发现他们的函数关系。
test=data1%>%select(region,year,pgdp:urblevel,labels)%>%
mutate(price=data2$price)%>%
mutate(construction=data2$construction)
plots=function(x){
ggplot(test,aes(x,price))+
geom_point()
}
plots(test$pgdp)
plots(test$indlevel)
plots(test$urblevel)
plots(test$construction)
由散点图大致可以看出，房价与人均GDP、工业化水平、城镇化水平可以是指数的关系，与商品住宅施工面积可能是线性相关关系。因此，首先考虑模型：
lnprice=β_0+β_1 lnpgdp+β_2 lnindlevel+β_3 lnurblevel+β_4construction
mod=test%>%group_by(region,labels)%>%
nest()%>%
mutate(model=map(data,~lm(log(price)~log(pgdp)+log(indlevel)+log(urblevel)+log(construction),data=.)))%>%
mutate(resids=map2(data,model,add_residuals))
mod1=unnest(mod,resids)
ggplot(mod1,aes(year,resid))+
geom_line(aes(group=region),alpha=1/3)+
geom_smooth(se=F)
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
根据该模型的残差可视化图形可以看出，该模型存在区域差异。因此，按照聚类分析的结果作出各地区分别的残差分布图，同时作出各地区的拟合优度分布图。
ggplot(mod1,aes(year,resid))+
geom_line(aes(group=region),alpha=1/3)+
geom_smooth(se=F)+
facet_wrap(~labels)
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
mod=mod%>%
mutate(glance=map(model,glance))
mod2=unnest(mod,glance,drop=T)
ggplot(mod2,aes(labels,r.squared))+
geom_jitter(width=0.5)
结合各地区的残差分布图和拟合优度分布可以看出，贫困地区出现了异常，故模型中去掉西藏省。除了西藏省之外其他各地区的模型拟合优度基本在0.95以上，可见总体上模型拟合较好。
画出分地区残差分布的箱型图，可以看到四个地区中，发达地区和普通地区的波动相对较大，而较发达地区和普通地区均存少量异常值，删除这些异常值，然后对整体建立一个模型。
ggplot(filter(mod1,region!="xizang"),aes(labels,resid))+
geom_boxplot(aes(group=labels))
mod=mod%>%unnest(resids)%>%
filter(abs(resid)<0.1)%>%
filter(region!="xizang")%>%
select(price,pgdp,indlevel,urblevel,construction)
模型诊断
模型拟合出来后，首先应当通过回归方程检验、回归系数检验、自变量多重共线性检验、残差独立性检验。回归方程显著性是指因变量与自变量线性回归关系是否显著，回归方程检验即检验回归系数是否全为零；回归系数检验是指检验某个变量是否与因变量存在线性关系；残差独立性是指残差对各个自变量是否相互独立；自变量多重共线性是指在模型中可能存在线性相关的自变量，这可能导致回归系数的符号与实际问题相违背。
回归方程和回归系数检验
model=lm(log(price)~log(pgdp)+log(indlevel)+log(urblevel)+log(construction),data=mod)
summary(model)
##
## Call:
## lm(formula = log(price) ~ log(pgdp) + log(indlevel) + log(urblevel) +
## log(construction), data = mod)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.71646 -0.19730 -0.00753 0.19782 0.92800
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.69883 0.42361 1.650 0.099857 .
## log(pgdp) 0.77916 0.03954 19.705 < 2e-16 ***
## log(indlevel) -0.96922 0.31192 -3.107 0.002035 **
## log(urblevel) 0.30396 0.08166 3.722 0.000228 ***
## log(construction) -0.03794 0.01294 -2.933 0.003566 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.271 on 368 degrees of freedom
## Multiple R-squared: 0.7717, Adjusted R-squared: 0.7693
## F-statistic: 311.1 on 4 and 368 DF, p-value: < 2.2e-16
从回归结果可以看到，回归方程F检验值为311.1，p值接近0，故该模型通过回归方程检验。同时，除了截距项外各变量的p值都显著的小于0.01，这说明解释变量的系数显著不为0，回归系数检验通过。
多重共线性检验
一般用条件数k=max(λ)/min(λ)来估计解释变量的共线性,λ为x'x的特征值。一般认为k<100的共线性较小，若k>1000则存在严重的多重共线性。
kappa(cor(mod[,2:5]))
## [1] 12.99204
k=12.99<100，所以认为四个变量之间不存在多重共线性。
残差分析
用最小二乘回归建模型如果效果显著，那么它的残差服从正态分布、等方差、独立性。下面分别检验该模型的残差正态性、等方差、独立性。
残差的正态性是很重要的，如果建立的模型不符合正态性假设，如何处理数据使之符合正态性假设又成为回归分析的另一个重要问题。绘制残差的Q-Q图，图中显示残差正态性效果不太好，一头一尾偏差较大。使用R语言中的shapiro.test函数来判断随机变量是否符合正态分布。
qqnorm(model$residuals)
qqline(model$residuals)
shapiro.test(rstandard(model))
##
## Shapiro-Wilk normality test
##
## data: rstandard(model)
## W = 0.98527, p-value = 0.0007558
结果显示p值非常接近0，强烈否定零假设，可见该模型的残差正态性效果并不好。
同方差是指方差与因变量相互独立，如果模型中残差同方差性比较好，那么残差就不随因变量的变化发生改变。使用car包中的ncvTest函数进行残差同方差检验。
ncvTest(model)
## Non-constant Variance Score Test
## Variance formula: ~ fitted.values
## Chisquare = 45.6315, Df = 1, p = 1.4273e-11
p值非常小拒绝误差项具有相同方差的零假设，因此模型未通过等方差检验。它的方差随因变量的改变而变化，表明除了己考虑的自变量因素，房价还受到别的因素因素。前文对房价影响因素的定性探索中已经发现房价还与人口总量、地区和利率等有一定关系。
残差独立性指残差序列内部是否相关。D-W统计量可以用来检验残差是否自相关。
durbinWatsonTest(model)
## lag Autocorrelation D-W Statistic p-value
## 1 0.8729439 0.2457334 0
## Alternative hypothesis: rho != 0
经过检验，在置信水平为0.05的条件下，未通过独立性检验。