[关闭]
@fanxy 2016-09-11T06:44:55.000000Z 字数 6274 阅读 2674

第二讲 经济周期

樊潇彦 复旦大学经济学院 中级宏观

  1. # 准备工作
  2. install.packages("readxl")    # 读取excel数据
  3. install.packages("stringr")   # 处理字符串
  4. install.packages("corrplot")  # 作相关关系图
  5. install.packages("igraph")    # 作网络图
  6. install.packages("forecast")  # 作季节调整
  8. library(readxl)               # 读取excel数据
  9. library(stringr)              # 字符串处理
  10. library(corrplot)
  11. library(igraph)
  12. library(forecast) 
  13. library(stats)                # 基础包，不用安装直接调用
  15. library(dplyr)
  16. library(tidyr)
  17. library(data.table)  
  18. library(foreign)
  19. library(readstata13)
  20. library(haven)   
  21. library(ggplot2)
  22. library(ggrepel)
  23. library(dygraphs)     
  24. library(plotrix)      
  25. library(lubridate)    
  26. library(zoo)          
  27. library(mFilter)      
  29. setwd("D:\\...")              # 设定工作目录
  30. rm(list=ls())

1. 趋势和周期的分解

  1. gdp=read_excel("Ch02_Data.xlsx",col_names=T,sheet="gdp_idx") # 读取excel数据
  2. gdp=gdp%>%
  3.   filter(year<=2014)%>%
  4.   mutate(gdp_idx52=ifelse(year==1952,100,gdp_idx))%>%  
  5.   mutate(gdp_idx52=cumprod(gdp_idx52/100))%>%                # 1952年为基期的实际GDP指数
  6.   mutate(cg_idx52=ifelse(year==1952,100,gdpcg_gr+100))%>%
  7.   mutate(cg_idx52=cumprod(cg_idx52/100))%>%                  # 实际最终消费指数
  8.   mutate(i_idx52=ifelse(year==1952,100,gdpi_gr+100))%>%
  9.   mutate(i_idx52=cumprod(i_idx52/100))%>%                    # 实际资本形成指数
  10.   mutate(GDP52=gdp_idx52*gdp$GDP[1])%>%                      # 实际GDP
  11.   mutate(deflator=GDP/GDP52)                                 # GDP平减指数
  13. # 以下用HP滤波提取趋势和周期，年度数据lambda=100，季度和月度分别为1600和14400
  14. # HP滤波得到有10个元素的列表（list），把cycle提取出来
  15. gdp$gdp=hpfilter(log(gdp$gdp_idx52),drift=F,freq=100,type="lambda")$cycle
  16. gdp$cg=hpfilter(log(gdp$cg_idx52),drift=F,freq=100,type="lambda")$cycle
  17. gdp$i=hpfilter(log(gdp$i_idx52),drift=F,freq=100,type="lambda")$cycle
  18. gdp$def=hpfilter(log(gdp$deflator),drift=F,freq=100,type="lambda")$cycle
  19. round(cor(gdp[,c("gdp","cg","i","def")]),2)                  # 查看相关系数
  21. # 比较两种方法提取出来的周期因素
  22. gdp$gr=gdp$gdp_idx/100-1                                     # 计算实际GDP增长率
  23. twoord.plot(lx = gdp$year, ly = gdp$gdp, 
  24.             rx = gdp$year,  ry = gdp$gr, 
  25.             main = '两种周期提取方法的结果比较', xlab = '', 
  26.             ylab = 'HP滤波', rylab = '实际GDP增长率', 
  27.             type = c('line','line'))

2. 特征事实

  1. # 产出、最终消费和资本形成
  2. g_gdp=gdp%>%select(year,gdp:i)%>%gather(var,cycle,-year)    # 选择要作图的数据
  3. ggplot(g_gdp,aes(year,cycle,color=var))+geom_line(size=1)+ 
  4.   labs(title="产出、最终消费和资本形成",x="",y="")+ 
  5.   scale_colour_manual(values= c("cg"="green","gdp"="black","i"="red"),
  6.                         labels = c('实际最终消费','实际GDP','实际资本形成'))+
  7.   guides(color = guide_legend(title = NULL)) +  
  8.   theme_bw()+ theme(legend.position = 'bottom') 
  10. # 产出和GDP平减指数
  11. g_gdp=gdp%>%select(year,gdp,def)%>%gather(var,cycle,-year)
  12. ggplot(g_gdp,aes(year,cycle,color=var))+geom_line(size=1)+ 
  13.   labs(title="产出和GDP平减指数",x="",y="")+ 
  14.   scale_colour_manual(values= c("def"="red","gdp"="black"),
  15.                         labels = c('GDP平减指数','实际GDP'))+
  16.   guides(color = guide_legend(title = NULL)) +  
  17.   theme_bw()+ theme(legend.position = 'bottom')

3. 财政支出和货币供给

  1. # 财政支出
  2. G=read_excel("Ch02_Data.xlsx",col_names=T,sheet="G") 
  3. str(G)
  4. G=G%>%
  5.   mutate(year=year(year))%>%
  6.   filter(year<=2014)%>%
  7.   left_join(gdp,by="year")%>%
  8.   mutate(all=`国家财政支出`/deflator)%>%                          # 计算实际值
  9.   mutate(cen=`中央财政支出`/deflator)%>%
  10.   mutate(loc=`地方财政支出`/deflator)%>%
  11.   select(year,gdp,all,cen,loc)
  12. G$all=hpfilter(log(G$all),drift=F,freq=100,type="lambda")$cycle
  13. G$cen=hpfilter(log(G$cen),drift=F,freq=100,type="lambda")$cycle
  14. G$loc=hpfilter(log(G$loc),drift=F,freq=100,type="lambda")$cycle
  15. round(cor(G[,2:5]),2)                                             # 查看相关系数
  17. G=G %>% select(year,gdp,all)%>%gather(var,cycle,-year)
  18. ggplot(G,aes(year,cycle,color=var))+geom_line(size=1)+ 
  19.   labs(title="产出和财政支出",x="",y="")+ 
  20.   scale_colour_manual(values= c("all"="red","gdp"="black"),
  21.                         labels = c('实际财政支出','实际GDP'))+
  22.   guides(color = guide_legend(title = NULL)) +  
  23.   theme_bw()+ theme(legend.position = 'bottom') 
  25. # 货币供给
  26. M=read_excel("Ch02_Data.xlsx",col_names=T,sheet="M")
  27. str(M)
  28. # 对GDP做季度调整后HP滤波
  29. gdp_ss=M%>%
  30.   select(time,`GDP:累计值`)%>%
  31.   mutate(time=as.yearmon(time))%>%
  32.   arrange(time)%>%
  33.   filter(!is.na(`GDP:累计值`))%>%
  34.   mutate(gdp=ts(`GDP:累计值`,frequency=4, start=c(1992,3)))   # 转换为ts对象
  35. lngdp_dec=stats::decompose(log(gdp_ss$gdp))                   # 分解为趋势、季节和随机因素
  36. plot(lngdp_dec)                                               # 分解结果作图
  37. lngdp_adj=seasadj(lngdp_dec)                                  # 生成经季节调整的数据
  38. plot(log(gdp_ss$gdp)); lines(lngdp_adj, col="red")            # 比较季节调整前后的数据
  39. gdp_ss$hp_gdp=hpfilter(lngdp_adj,drift=F,freq=1600,type="lambda")$cycle  # 对季节调整数据作图
  40. plot(gdp_ss$hp_gdp)                                           # HP滤波结果作图
  42. # 对货币供应做月度调整后HP滤波
  43. mon=M%>%
  44.   select(time,`M0(流通中的现金):期末值`,`M2(货币和准货币):期末值`)%>%
  45.   mutate(time=as.yearmon(time))%>%
  46.   arrange(time)%>%
  47.   filter(time>=1996)%>%
  48.   mutate(m0=ts(`M0(流通中的现金):期末值`,frequency=12, start=c(1996,1)))%>%  # 转换为ts对象
  49.   mutate(m2=ts(`M2(货币和准货币):期末值`,frequency=12, start=c(1996,1))) 
  50. str(mon)  
  51. lnm0=stats::decompose(log(mon$m0)); plot(lnm0)
  52. mon$hp_m0=hpfilter(seasadj(lnm0),drift=F,freq=14400,type="lambda")$cycle
  53. lnm2=stats::decompose(log(mon$m2)); plot(lnm2)
  54. mon$hp_m2=hpfilter(seasadj(lnm2),drift=F,freq=14400,type="lambda")$cycle
  56. mon=mon%>%
  57.   left_join(gdp_ss,by="time")%>%
  58.   select(time, contains("hp_"))%>%
  59.   mutate(hp_gdp=ifelse(is.na(hp_gdp),0,hp_gdp)) 
  60. ggplot(mon,aes(time,hp_gdp)) + 
  61.   geom_line(aes(time,hp_m0,color="M0"),size=1)+geom_line(aes(time,hp_m2,color="M2"),size=1)+
  62.   scale_colour_manual(values= c("M0"="orange","M2"="red"),
  63.                         labels = c('M0','M2'))+
  64.   geom_bar(stat="identity", position = "dodge",fill="gray40")+
  65.   labs(title="产出和货币供给",x="",y="")+ 
  66.   guides(color = guide_legend(title = NULL)) +  
  67.   theme_bw()+ theme(legend.position = 'bottom')

4. 产业与地区关联

  1. # 房地产市场与宏观经济
  2. Guofang=read_excel("Ch02_Data.xlsx",col_names=T,sheet="Guofang") 
  3. str(Guofang)
  4. Guofang=Guofang%>%
  5.   mutate(time=sub("月","",time))%>%
  6.   mutate(time=sub("年","-",time))%>%
  7.   mutate(time=as.yearmon(time))%>%
  8.   gather(var,idx,-time)%>%
  9.   arrange(var,time)
  10. ggplot(Guofang,aes(time,idx,color=var))+geom_line(size=1)+
  11.   labs(title="国房景气指数和宏观经济景气指数",x="",y="")+ 
  12.   guides(color = guide_legend(title = NULL)) +  
  13.   theme_bw()+ theme(legend.position = 'bottom') 
  15. # 行业相关性
  16. indcor=read_excel("Ch02_Data.xlsx",col_names=T,sheet="indgdp_idx")
  17. indcor123=indcor%>%
  18.   select(GDP,`第一产业`,`第二产业`,`第三产业`)
  19. round(cor(indcor123),2)                                 # 查看相关性
  21. corr=cor(select(indcor,-year,-GDP,-`第二产业`,-`第三产业`))
  22. corrplot(corr = corr,type="upper",tl.pos="d",order="hclust",method="ellipse")
  23. corrplot(corr = corr,add=TRUE, type="lower", method="number",order="hclust",diag=FALSE,tl.pos="n", cl.pos="n")
  25. # 地区相关性
  26. gdpreg=read_excel("Ch02_Data.xlsx",col_names=T,sheet="reggdp_idx")
  27. gdpreg=gdpreg%>%filter(!is.na(reg))
  28. reginfo=gdpreg%>%mutate(regno=1:31)%>%select(reg,regno)
  29. regcor=gdpreg%>%
  30.   gather(year,idx,-reg)%>%
  31.   filter(year>1952)%>%
  32.   spread(reg,idx)
  34. regcor=cor(select(regcor,-year))
  36. corrplot(corr = regcor,type="lower",tl.pos="d",         # 地区相关关系图
  37.          order="hclust",method="shade")  
  38. # 每个地区选取与之相关性最高的地区，网络作图
  39. diag(regcor)=0
  40. maxregcor=data.frame(from=colnames(regcor),
  41.                to=colnames(regcor)[apply(regcor,1,which.max)], 
  42.                cor=apply(regcor,1,max))
  43. rownames(maxregcor)=NULL
  44. g_regcor <- graph.data.frame(maxregcor, directed = T)
  45. summary(g_regcor)
  46. plot.igraph(g_regcor, 
  47.             vertex.size = g_regcor$weight, 
  48.             edge.arrow.size=0.5, 
  49.             edge.color = g_regcor$weight, 
  50.             edge.width= g_regcor$weight)
添加新批注
在作者公开此批注前，只有你和作者可见。
回复批注