[关闭]
@fanxy 2018-04-23T15:52:22.000000Z 字数 3991 阅读 2731

二、动态最优数学基础

樊潇彦 复旦大学经济学院 数量经济学


0. 准备工作

数据下载:
CME_Qgdp.txt6.7kB

  1. setwd("D:\\...")
  2. rm(list=ls())
  3. install.packages(c("nleqslv","deSolve","scatterplot3d"))
  4. install.packages(c("zoo","quantmod","tseries","fBasics","forecast")) # 时间序列常用包
  5. library(dplyr)
  6. library(tidyr)
  7. library(ggplot2)
  8. library(nleqslv)
  9. library(deSolve)
  10. library(scatterplot3d)

1. 差(微)分方程

1.1 Solow-Swan模型

  1. delta=0.05; s=0.4; A=1; alpha=0.3 # 设定参数
  2. K=1; times = seq(0, 100, by = 1) # 初始值和时间轴
  3. # 求稳态资本 Kstar
  4. FunKstar<-function(Kstar) {
  5. error=s*A*Kstar^alpha - delta*Kstar # K_t+1 - K_t =0
  6. return(error)
  7. }
  8. library(nleqslv) # 调用程序包,内含命令nleqslv求解非线性方程(组)
  9. Kstar=nleqslv(5,FunKstar)$x
  10. # 解差分方程
  11. disK=rep(K,length(times))
  12. for (t in 2:length(times)) disK[t]=s*A*disK[t-1]^alpha + (1-delta)*disK[t-1]
  13. disK=data.frame(times=times,`差分方程`=disK)
  14. # 解微分方程
  15. Kfun=function(t,K,par){
  16. with(as.list(c(K,par)),{
  17. dK= s*A*K^alpha - delta*K
  18. list(dK)
  19. })
  20. }
  21. conK= ode(y=c(K=1), times=times, func=Kfun, # ode:程序包deSolve中命令,解常微分方程
  22. parms=c(delta=delta,s=s,A=A,alpha=alpha))
  23. # 合并数据,作图
  24. Capital=data.frame(times=times, `微分方程`=conK[,-1])%>%
  25. left_join(disK,by="times")
  26. tail(Capital,5) # 最后5期结果比较
  27. ggplot(Capital%>%gather(var,K,-times), aes(times,K, color=var))+geom_line()+
  28. geom_hline(yintercept=Kstar,linetype = "dotdash",col="red")+
  29. labs(title="Solow-Swan模型",x="时期",y="人均资本")+
  30. guides(color=guide_legend(title=NULL))+theme_bw()+theme(legend.position="bottom")

Solow_Swan.bmp-1912.6kB

1.2 洛伦兹方程与混沌

  1. parameters <- c(a = -8/3, b = -10, c = 28)
  2. state <- c(X = 1, Y = 1, Z = 1)
  3. times <- seq(0, 100, by = 0.01)
  4. Lorenz <- function(t, state, parameters) {
  5. with(as.list(c(state, parameters)), {
  6. dX <- a * X + Y * Z
  7. dY <- b * (Y - Z)
  8. dZ <- -X * Y + c * Y - Z
  9. list(c(dX, dY, dZ))
  10. })
  11. }
  12. out <- ode(y = state, times = times, func = Lorenz, parms = parameters)
  13. scatterplot3d(out[,-1], type = "l", highlight.3d=TRUE, col.axis="blue", col.grid="lightblue")

Lorenz.bmp-2626.2kB

2. 马尔可夫链与稳态人口分布

  1. # 自已写一个可以计算矩阵的幂的函数
  2. mat_power = function(A, n){
  3. Apower=A
  4. for (i in 2:n) Apower= Apower %*% A
  5. return(Apower)
  6. }
  7. P=matrix(c(0.95,0.05,0.03,0.97),nrow=2,byrow=T) # 人口迁移的转移矩阵
  8. P_star=mat_power(P,300) # 计算稳态转移矩阵
  9. round(eigen(t(P))$values,2) # 求特征值
  10. ev=eigen(t(P))$vectors # 求左特征向量
  11. x=ev[,1]/sum(ev[,1]) # 将特征值为1的特征向量正规化
  12. # 比较三种稳态人口比例
  13. x
  14. x%*%P
  15. P_star[1,]

3. 随机过程(选读,更多材料参见《金融数据处理与分析》)

3.1 获取宏观和金融数据

学校电子图书馆有很多中文数据库,下面以国泰安中季度国内生产总值为例,说明宏观经济数据的特征。

  1. Qgdp=read.table("CME_Qgdp.txt", header = T, stringsAsFactors =F)
  2. # summary(Qgdp); str(Qgdp) # 查看数据结构
  3. library(zoo)
  4. Qgdp=Qgdp%>%mutate(time=as.yearmon(Staper))%>%
  5. rename(`季度GDP(累计值)`=Gdpq0101)%>%rename(`同比增长率(名义)`=Gdpq0105)%>%
  6. select(time,`季度GDP(累计值)`,`同比增长率(名义)`)
  7. ggplot(Qgdp%>%gather(var,value,-time),aes(time,value))+
  8. geom_line()+facet_wrap(~var,scales="free")+
  9. labs(title="中国宏观经济数据",x="",y="")+theme_bw()

Gdp_Q.bmp-2260.1kB

quantmod包中的getSymbols()函数从yahoo财经下载上证综指,更多指标可以查看相关网站。

  1. library(quantmod)
  2. ssec<-getSymbols("^SSEC",from = "1991-07-15",to = Sys.Date(),src = "yahoo",auto.assign=FALSE)
  3. dim(ssec)
  4. tail(ssec)
  5. plot(ssec$SSEC.Adjusted, main = "上证综指(SSEC)")
  6. ssec.m<-to.monthly(ssec) # 转换成月度数据
  7. chartSeries(ssec.m,type="candlesticks",theme="white", name = "上证综指月线图",time.scale = "",up.col="red",dn.col="green")

SSEC_M.bmp-2260.1kB

3.2 数据平稳性检验

  1. ln_ssec=log(ssec.m$ssec.Adjusted)
  2. # 计算月收益率,由于月收益率可能较大,不用diff(log(x))
  3. ssec.m$ret=diff(ssec.m$ssec.Adjusted) / lag(ssec.m$ssec.Adjusted) * 100
  4. ret=ssec.m$ret[-1,] # 去掉第一行缺失值
  5. plot(ret, main="上证综指月收益率(%)")
  6. library(tseries)
  7. pp.test(ln_ssec) # DF检验,H0:有单位根
  8. pp.test(ret)
  9. adf.test(ln_ssec) # ADF检验,H0:有单位根
  10. adf.test(ret)
  11. kpss.test(ln_ssec) # KPSS检验,H0:平稳
  12. kpss.test(ret)

3.3 平稳时间序列的ARMA模型

实际应用中,forecast 包中的 auto.arima() 命令,可以自动识别 ARIMA(p,d,q) 中各个阶数,非常方便。

  1. library(fBasics)
  2. basicStats(ret) # 对月度收益率的基本统计描述
  3. t.test(ret) # 均值检验 H0: x=0
  4. normalTest(ret,method='jb',na.rm=T) # 正态分布检验,JB-test
  5. par(mfcol=c(1,1))
  6. acf(ret, lag=12, plot=T) # 自相关系数作图
  7. Box.test(ret,lag=6,type='Ljung') # m阶自相关系数为零的检验
  8. library(forecast)
  9. stock=auto.arima(ret,stationary = TRUE, seasonal = FALSE,
  10. ic="aic")
  11. # 拟合与预测
  12. par(mfcol=c(2,1))
  13. plot(stock$x, lty = 1, main = "上证综指收益率(%)")
  14. lines(fitted(stock), lty = 2,lwd = 2, col = "red")
  15. predict(stock, n.ahead=3)
  16. plot(forecast(stock))

est_SSEC_M.bmp-2260.1kB

添加新批注
在作者公开此批注前,只有你和作者可见。
回复批注