[关闭]
@fanxy 2017-05-11T15:57:09.000000Z 字数 2941 阅读 2209

四、动态最优(II)最优控制方法

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


0. 准备工作

  1. setwd("D:\\...") # 设定工作目录
  2. rm(list=ls()) # 清内存
  3. # "Ctrl"+"L" 清除命令窗口
  4. library(ggplot2)
  5. library(nleqslv)

对于连续时间的最优增长问题:


写出汉密尔顿函数:

根据最大值原理有:

前两个条件合并得到:

1. 设定参数

  1. rho <- 0.03 # 贴现因子
  2. gamma <- 1 # 效用函数
  3. if (gamma==1){u.fmla <- expression(log(c))
  4. }else{
  5. u.fmla <- expression(c^(1-gamma)/(1-gamma))}
  6. alpha <- 0.5 # 生产函数
  7. f.fmla <- expression(k^alpha)
  8. delta <- 0.1 # 折旧率

2. 微分方程组与稳态

  1. d2u.fmla <- deriv3(u.fmla,"c","c") # 效用函数求导
  2. u <- function(c) { eval(u.fmla) } # u
  3. du <- function(c) { as.vector(attr(d2u.fmla(c),"gradient")) } # u'
  4. d2u <- function(c) { as.vector(attr(d2u.fmla(c),"hessian")) } # u"
  5. df.fmla <- deriv(f.fmla,"k",function(k) {}) # 生产函数求导
  6. f <- function(k) { eval(f.fmla) } # y
  7. df <- function(k) { as.vector(attr(df.fmla(k),"gradient")) } # y'
  8. dcdt <- function(c,k) { # 微分方程 dc/dt
  9. -du(c)/d2u(c)*(df(k)-delta-rho)
  10. }
  11. dkdt <- function(c,k) { # 微分方程 dk/dt
  12. f(k) -c -delta*k
  13. }
  14. obj <- function(x) { # 联立微分方程组
  15. c.old <- x[1]
  16. k.old <- x[2]
  17. c(c.old+dcdt(c.old,k.old)*dt - c,
  18. k.old+dkdt(c.old,k.old)*dt - k)
  19. }
  20. kss <- ((rho+delta)/alpha)^(1/(alpha-1)) # 稳态
  21. css <- f(kss) - delta*kss

3. 搜索鞍点路径

  1. c.grid <- seq(0.01,2*css,length.out=20) # 设置取值范围和格点
  2. k.grid <- seq(0.01,2*kss,length.out=20)
  3. c.stable =vector()
  4. k.stable =vector()
  5. s=1
  6. c.stable[s] <- css; k.stable[s] <- kss
  7. dt <- 0.1
  8. # 左侧鞍点路径
  9. c <- css; k <- kss # 从稳态出发
  10. while (k>min(k.grid)) {
  11. left <- nleqslv(x=c(c-dt,k-dt),fn=obj) # 解联立方程组
  12. c <- left$x[1]
  13. k <- left$x[2]
  14. c.stable[s] <- c
  15. k.stable[s] <- k
  16. s <- s+1
  17. }
  18. # 右侧鞍点路径
  19. c <- css; k <- kss
  20. while (k<max(k.grid)) {
  21. right <- nleqslv(x=c(c+dt,k+dt),fn=obj)
  22. c <- right$x[1]
  23. k <- right$x[2]
  24. c.stable[s] <- c
  25. k.stable[s] <- k
  26. s <- s+1
  27. }

4. 相图

  1. d <- data.frame(c=as.vector(outer(0*k.grid,c.grid,FUN="+")),
  2. k=as.vector(outer(k.grid,0*c.grid,FUN="+")))
  3. d$dc <- dcdt(d$c,d$k)
  4. d$dk <- dkdt(d$c,d$k)
  5. d$zerodk <- f(d$k)-delta*d$k # dk=0的稳态分界线
  6. stable <- data.frame(k=k.stable,c=c.stable) # 鞍点路径
  7. ggplot(data=d) +
  8. geom_segment(data=d,aes(x=k,y=c,yend=(c+dc),xend=(k+dk)),
  9. arrow=arrow(length=unit(0.1,"cm")),colour="gray") +
  10. geom_line(aes(x=k,y=zerodk),colour="blue") +
  11. geom_vline(xintercept=kss,colour="red") + # dc=0的稳态分界线
  12. geom_line(data=stable,aes(x=k,y=c),colour="black") +
  13. xlim(min(k.grid),max(k.grid)) +
  14. ylim(min(c.grid),max(c.grid)) + theme_minimal()

phase.jpg-93.7kB

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