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

第二讲 经济周期

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


  1. # 准备工作
  2. install.packages("readxl") # 读取excel数据
  3. install.packages("stringr") # 处理字符串
  4. install.packages("corrplot") # 作相关关系图
  5. install.packages("igraph") # 作网络图
  6. install.packages("forecast") # 作季节调整
  7. library(readxl) # 读取excel数据
  8. library(stringr) # 字符串处理
  9. library(corrplot)
  10. library(igraph)
  11. library(forecast)
  12. library(stats) # 基础包,不用安装直接调用
  13. library(dplyr)
  14. library(tidyr)
  15. library(data.table)
  16. library(foreign)
  17. library(readstata13)
  18. library(haven)
  19. library(ggplot2)
  20. library(ggrepel)
  21. library(dygraphs)
  22. library(plotrix)
  23. library(lubridate)
  24. library(zoo)
  25. library(mFilter)
  26. setwd("D:\\...") # 设定工作目录
  27. 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平减指数
  12. # 以下用HP滤波提取趋势和周期,年度数据lambda=100,季度和月度分别为1600和14400
  13. # HP滤波得到有10个元素的列表(list),把cycle提取出来
  14. gdp$gdp=hpfilter(log(gdp$gdp_idx52),drift=F,freq=100,type="lambda")$cycle
  15. gdp$cg=hpfilter(log(gdp$cg_idx52),drift=F,freq=100,type="lambda")$cycle
  16. gdp$i=hpfilter(log(gdp$i_idx52),drift=F,freq=100,type="lambda")$cycle
  17. gdp$def=hpfilter(log(gdp$deflator),drift=F,freq=100,type="lambda")$cycle
  18. round(cor(gdp[,c("gdp","cg","i","def")]),2) # 查看相关系数
  19. # 比较两种方法提取出来的周期因素
  20. gdp$gr=gdp$gdp_idx/100-1 # 计算实际GDP增长率
  21. twoord.plot(lx = gdp$year, ly = gdp$gdp,
  22. rx = gdp$year, ry = gdp$gr,
  23. main = '两种周期提取方法的结果比较', xlab = '',
  24. ylab = 'HP滤波', rylab = '实际GDP增长率',
  25. 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')
  9. # 产出和GDP平减指数
  10. g_gdp=gdp%>%select(year,gdp,def)%>%gather(var,cycle,-year)
  11. ggplot(g_gdp,aes(year,cycle,color=var))+geom_line(size=1)+
  12. labs(title="产出和GDP平减指数",x="",y="")+
  13. scale_colour_manual(values= c("def"="red","gdp"="black"),
  14. labels = c('GDP平减指数','实际GDP'))+
  15. guides(color = guide_legend(title = NULL)) +
  16. 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) # 查看相关系数
  16. G=G %>% select(year,gdp,all)%>%gather(var,cycle,-year)
  17. ggplot(G,aes(year,cycle,color=var))+geom_line(size=1)+
  18. labs(title="产出和财政支出",x="",y="")+
  19. scale_colour_manual(values= c("all"="red","gdp"="black"),
  20. labels = c('实际财政支出','实际GDP'))+
  21. guides(color = guide_legend(title = NULL)) +
  22. theme_bw()+ theme(legend.position = 'bottom')
  23. # 货币供给
  24. M=read_excel("Ch02_Data.xlsx",col_names=T,sheet="M")
  25. str(M)
  26. # 对GDP做季度调整后HP滤波
  27. gdp_ss=M%>%
  28. select(time,`GDP:累计值`)%>%
  29. mutate(time=as.yearmon(time))%>%
  30. arrange(time)%>%
  31. filter(!is.na(`GDP:累计值`))%>%
  32. mutate(gdp=ts(`GDP:累计值`,frequency=4, start=c(1992,3))) # 转换为ts对象
  33. lngdp_dec=stats::decompose(log(gdp_ss$gdp)) # 分解为趋势、季节和随机因素
  34. plot(lngdp_dec) # 分解结果作图
  35. lngdp_adj=seasadj(lngdp_dec) # 生成经季节调整的数据
  36. plot(log(gdp_ss$gdp)); lines(lngdp_adj, col="red") # 比较季节调整前后的数据
  37. gdp_ss$hp_gdp=hpfilter(lngdp_adj,drift=F,freq=1600,type="lambda")$cycle # 对季节调整数据作图
  38. plot(gdp_ss$hp_gdp) # HP滤波结果作图
  39. # 对货币供应做月度调整后HP滤波
  40. mon=M%>%
  41. select(time,`M0(流通中的现金):期末值`,`M2(货币和准货币):期末值`)%>%
  42. mutate(time=as.yearmon(time))%>%
  43. arrange(time)%>%
  44. filter(time>=1996)%>%
  45. mutate(m0=ts(`M0(流通中的现金):期末值`,frequency=12, start=c(1996,1)))%>% # 转换为ts对象
  46. mutate(m2=ts(`M2(货币和准货币):期末值`,frequency=12, start=c(1996,1)))
  47. str(mon)
  48. lnm0=stats::decompose(log(mon$m0)); plot(lnm0)
  49. mon$hp_m0=hpfilter(seasadj(lnm0),drift=F,freq=14400,type="lambda")$cycle
  50. lnm2=stats::decompose(log(mon$m2)); plot(lnm2)
  51. mon$hp_m2=hpfilter(seasadj(lnm2),drift=F,freq=14400,type="lambda")$cycle
  52. mon=mon%>%
  53. left_join(gdp_ss,by="time")%>%
  54. select(time, contains("hp_"))%>%
  55. mutate(hp_gdp=ifelse(is.na(hp_gdp),0,hp_gdp))
  56. ggplot(mon,aes(time,hp_gdp)) +
  57. geom_line(aes(time,hp_m0,color="M0"),size=1)+geom_line(aes(time,hp_m2,color="M2"),size=1)+
  58. scale_colour_manual(values= c("M0"="orange","M2"="red"),
  59. labels = c('M0','M2'))+
  60. geom_bar(stat="identity", position = "dodge",fill="gray40")+
  61. labs(title="产出和货币供给",x="",y="")+
  62. guides(color = guide_legend(title = NULL)) +
  63. 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')
  14. # 行业相关性
  15. indcor=read_excel("Ch02_Data.xlsx",col_names=T,sheet="indgdp_idx")
  16. indcor123=indcor%>%
  17. select(GDP,`第一产业`,`第二产业`,`第三产业`)
  18. round(cor(indcor123),2) # 查看相关性
  19. corr=cor(select(indcor,-year,-GDP,-`第二产业`,-`第三产业`))
  20. corrplot(corr = corr,type="upper",tl.pos="d",order="hclust",method="ellipse")
  21. corrplot(corr = corr,add=TRUE, type="lower", method="number",order="hclust",diag=FALSE,tl.pos="n", cl.pos="n")
  22. # 地区相关性
  23. gdpreg=read_excel("Ch02_Data.xlsx",col_names=T,sheet="reggdp_idx")
  24. gdpreg=gdpreg%>%filter(!is.na(reg))
  25. reginfo=gdpreg%>%mutate(regno=1:31)%>%select(reg,regno)
  26. regcor=gdpreg%>%
  27. gather(year,idx,-reg)%>%
  28. filter(year>1952)%>%
  29. spread(reg,idx)
  30. regcor=cor(select(regcor,-year))
  31. corrplot(corr = regcor,type="lower",tl.pos="d", # 地区相关关系图
  32. order="hclust",method="shade")
  33. # 每个地区选取与之相关性最高的地区,网络作图
  34. diag(regcor)=0
  35. maxregcor=data.frame(from=colnames(regcor),
  36. to=colnames(regcor)[apply(regcor,1,which.max)],
  37. cor=apply(regcor,1,max))
  38. rownames(maxregcor)=NULL
  39. g_regcor <- graph.data.frame(maxregcor, directed = T)
  40. summary(g_regcor)
  41. plot.igraph(g_regcor,
  42. vertex.size = g_regcor$weight,
  43. edge.arrow.size=0.5,
  44. edge.color = g_regcor$weight,
  45. edge.width= g_regcor$weight)
添加新批注
在作者公开此批注前,只有你和作者可见。
回复批注