二、动态最优数学基础

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

---

0. 准备工作

数据下载：
CME_Qgdp.txt6.7kB

setwd("D:\\...")
rm(list=ls())
install.packages(c("nleqslv","deSolve","scatterplot3d"))
install.packages(c("zoo","quantmod","tseries","fBasics","forecast")) # 时间序列常用包
library(dplyr)
library(tidyr)
library(ggplot2)
library(nleqslv)
library(deSolve)
library(scatterplot3d)

1. 差（微）分方程

1.1 Solow-Swan模型

$${K_{t + 1}} = \left( {1 - \delta } \right){K_t} + sY_t = \left( {1 - \delta } \right){K_t} + sAK_t^\alpha$$

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

1.2 洛伦兹方程与混沌

$$\left\{ {\begin{array}{*{20}{l}} {\frac{{dx}}{{dt}} = a x + yz}\\ {\frac{{dy}}{{dt}} = b (y - z)}\\ {\frac{{dz}}{{dt}} = - xy + c y - z} \end{array}} \right.$$

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

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

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

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

3.1 获取宏观和金融数据

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

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

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

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

3.2 数据平稳性检验

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

3.3 平稳时间序列的ARMA模型

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

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