@fanxy
2016-09-11T14:44:55.000000Z
字数 6274
阅读 2312
樊潇彦
复旦大学经济学院
中级宏观
# 准备工作
install.packages("readxl") # 读取excel数据
install.packages("stringr") # 处理字符串
install.packages("corrplot") # 作相关关系图
install.packages("igraph") # 作网络图
install.packages("forecast") # 作季节调整
library(readxl) # 读取excel数据
library(stringr) # 字符串处理
library(corrplot)
library(igraph)
library(forecast)
library(stats) # 基础包,不用安装直接调用
library(dplyr)
library(tidyr)
library(data.table)
library(foreign)
library(readstata13)
library(haven)
library(ggplot2)
library(ggrepel)
library(dygraphs)
library(plotrix)
library(lubridate)
library(zoo)
library(mFilter)
setwd("D:\\...") # 设定工作目录
rm(list=ls())
gdp=read_excel("Ch02_Data.xlsx",col_names=T,sheet="gdp_idx") # 读取excel数据
gdp=gdp%>%
filter(year<=2014)%>%
mutate(gdp_idx52=ifelse(year==1952,100,gdp_idx))%>%
mutate(gdp_idx52=cumprod(gdp_idx52/100))%>% # 1952年为基期的实际GDP指数
mutate(cg_idx52=ifelse(year==1952,100,gdpcg_gr+100))%>%
mutate(cg_idx52=cumprod(cg_idx52/100))%>% # 实际最终消费指数
mutate(i_idx52=ifelse(year==1952,100,gdpi_gr+100))%>%
mutate(i_idx52=cumprod(i_idx52/100))%>% # 实际资本形成指数
mutate(GDP52=gdp_idx52*gdp$GDP[1])%>% # 实际GDP
mutate(deflator=GDP/GDP52) # GDP平减指数
# 以下用HP滤波提取趋势和周期,年度数据lambda=100,季度和月度分别为1600和14400
# HP滤波得到有10个元素的列表(list),把cycle提取出来
gdp$gdp=hpfilter(log(gdp$gdp_idx52),drift=F,freq=100,type="lambda")$cycle
gdp$cg=hpfilter(log(gdp$cg_idx52),drift=F,freq=100,type="lambda")$cycle
gdp$i=hpfilter(log(gdp$i_idx52),drift=F,freq=100,type="lambda")$cycle
gdp$def=hpfilter(log(gdp$deflator),drift=F,freq=100,type="lambda")$cycle
round(cor(gdp[,c("gdp","cg","i","def")]),2) # 查看相关系数
# 比较两种方法提取出来的周期因素
gdp$gr=gdp$gdp_idx/100-1 # 计算实际GDP增长率
twoord.plot(lx = gdp$year, ly = gdp$gdp,
rx = gdp$year, ry = gdp$gr,
main = '两种周期提取方法的结果比较', xlab = '',
ylab = 'HP滤波', rylab = '实际GDP增长率',
type = c('line','line'))
# 产出、最终消费和资本形成
g_gdp=gdp%>%select(year,gdp:i)%>%gather(var,cycle,-year) # 选择要作图的数据
ggplot(g_gdp,aes(year,cycle,color=var))+geom_line(size=1)+
labs(title="产出、最终消费和资本形成",x="",y="")+
scale_colour_manual(values= c("cg"="green","gdp"="black","i"="red"),
labels = c('实际最终消费','实际GDP','实际资本形成'))+
guides(color = guide_legend(title = NULL)) +
theme_bw()+ theme(legend.position = 'bottom')
# 产出和GDP平减指数
g_gdp=gdp%>%select(year,gdp,def)%>%gather(var,cycle,-year)
ggplot(g_gdp,aes(year,cycle,color=var))+geom_line(size=1)+
labs(title="产出和GDP平减指数",x="",y="")+
scale_colour_manual(values= c("def"="red","gdp"="black"),
labels = c('GDP平减指数','实际GDP'))+
guides(color = guide_legend(title = NULL)) +
theme_bw()+ theme(legend.position = 'bottom')
# 财政支出
G=read_excel("Ch02_Data.xlsx",col_names=T,sheet="G")
str(G)
G=G%>%
mutate(year=year(year))%>%
filter(year<=2014)%>%
left_join(gdp,by="year")%>%
mutate(all=`国家财政支出`/deflator)%>% # 计算实际值
mutate(cen=`中央财政支出`/deflator)%>%
mutate(loc=`地方财政支出`/deflator)%>%
select(year,gdp,all,cen,loc)
G$all=hpfilter(log(G$all),drift=F,freq=100,type="lambda")$cycle
G$cen=hpfilter(log(G$cen),drift=F,freq=100,type="lambda")$cycle
G$loc=hpfilter(log(G$loc),drift=F,freq=100,type="lambda")$cycle
round(cor(G[,2:5]),2) # 查看相关系数
G=G %>% select(year,gdp,all)%>%gather(var,cycle,-year)
ggplot(G,aes(year,cycle,color=var))+geom_line(size=1)+
labs(title="产出和财政支出",x="",y="")+
scale_colour_manual(values= c("all"="red","gdp"="black"),
labels = c('实际财政支出','实际GDP'))+
guides(color = guide_legend(title = NULL)) +
theme_bw()+ theme(legend.position = 'bottom')
# 货币供给
M=read_excel("Ch02_Data.xlsx",col_names=T,sheet="M")
str(M)
# 对GDP做季度调整后HP滤波
gdp_ss=M%>%
select(time,`GDP:累计值`)%>%
mutate(time=as.yearmon(time))%>%
arrange(time)%>%
filter(!is.na(`GDP:累计值`))%>%
mutate(gdp=ts(`GDP:累计值`,frequency=4, start=c(1992,3))) # 转换为ts对象
lngdp_dec=stats::decompose(log(gdp_ss$gdp)) # 分解为趋势、季节和随机因素
plot(lngdp_dec) # 分解结果作图
lngdp_adj=seasadj(lngdp_dec) # 生成经季节调整的数据
plot(log(gdp_ss$gdp)); lines(lngdp_adj, col="red") # 比较季节调整前后的数据
gdp_ss$hp_gdp=hpfilter(lngdp_adj,drift=F,freq=1600,type="lambda")$cycle # 对季节调整数据作图
plot(gdp_ss$hp_gdp) # HP滤波结果作图
# 对货币供应做月度调整后HP滤波
mon=M%>%
select(time,`M0(流通中的现金):期末值`,`M2(货币和准货币):期末值`)%>%
mutate(time=as.yearmon(time))%>%
arrange(time)%>%
filter(time>=1996)%>%
mutate(m0=ts(`M0(流通中的现金):期末值`,frequency=12, start=c(1996,1)))%>% # 转换为ts对象
mutate(m2=ts(`M2(货币和准货币):期末值`,frequency=12, start=c(1996,1)))
str(mon)
lnm0=stats::decompose(log(mon$m0)); plot(lnm0)
mon$hp_m0=hpfilter(seasadj(lnm0),drift=F,freq=14400,type="lambda")$cycle
lnm2=stats::decompose(log(mon$m2)); plot(lnm2)
mon$hp_m2=hpfilter(seasadj(lnm2),drift=F,freq=14400,type="lambda")$cycle
mon=mon%>%
left_join(gdp_ss,by="time")%>%
select(time, contains("hp_"))%>%
mutate(hp_gdp=ifelse(is.na(hp_gdp),0,hp_gdp))
ggplot(mon,aes(time,hp_gdp)) +
geom_line(aes(time,hp_m0,color="M0"),size=1)+geom_line(aes(time,hp_m2,color="M2"),size=1)+
scale_colour_manual(values= c("M0"="orange","M2"="red"),
labels = c('M0','M2'))+
geom_bar(stat="identity", position = "dodge",fill="gray40")+
labs(title="产出和货币供给",x="",y="")+
guides(color = guide_legend(title = NULL)) +
theme_bw()+ theme(legend.position = 'bottom')
# 房地产市场与宏观经济
Guofang=read_excel("Ch02_Data.xlsx",col_names=T,sheet="Guofang")
str(Guofang)
Guofang=Guofang%>%
mutate(time=sub("月","",time))%>%
mutate(time=sub("年","-",time))%>%
mutate(time=as.yearmon(time))%>%
gather(var,idx,-time)%>%
arrange(var,time)
ggplot(Guofang,aes(time,idx,color=var))+geom_line(size=1)+
labs(title="国房景气指数和宏观经济景气指数",x="",y="")+
guides(color = guide_legend(title = NULL)) +
theme_bw()+ theme(legend.position = 'bottom')
# 行业相关性
indcor=read_excel("Ch02_Data.xlsx",col_names=T,sheet="indgdp_idx")
indcor123=indcor%>%
select(GDP,`第一产业`,`第二产业`,`第三产业`)
round(cor(indcor123),2) # 查看相关性
corr=cor(select(indcor,-year,-GDP,-`第二产业`,-`第三产业`))
corrplot(corr = corr,type="upper",tl.pos="d",order="hclust",method="ellipse")
corrplot(corr = corr,add=TRUE, type="lower", method="number",order="hclust",diag=FALSE,tl.pos="n", cl.pos="n")
# 地区相关性
gdpreg=read_excel("Ch02_Data.xlsx",col_names=T,sheet="reggdp_idx")
gdpreg=gdpreg%>%filter(!is.na(reg))
reginfo=gdpreg%>%mutate(regno=1:31)%>%select(reg,regno)
regcor=gdpreg%>%
gather(year,idx,-reg)%>%
filter(year>1952)%>%
spread(reg,idx)
regcor=cor(select(regcor,-year))
corrplot(corr = regcor,type="lower",tl.pos="d", # 地区相关关系图
order="hclust",method="shade")
# 每个地区选取与之相关性最高的地区,网络作图
diag(regcor)=0
maxregcor=data.frame(from=colnames(regcor),
to=colnames(regcor)[apply(regcor,1,which.max)],
cor=apply(regcor,1,max))
rownames(maxregcor)=NULL
g_regcor <- graph.data.frame(maxregcor, directed = T)
summary(g_regcor)
plot.igraph(g_regcor,
vertex.size = g_regcor$weight,
edge.arrow.size=0.5,
edge.color = g_regcor$weight,
edge.width= g_regcor$weight)