[关闭]
@fanxy 2017-05-23T02:28:56.000000Z 字数 2937 阅读 1560

五、动态最优(III)动态规划方法

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


0. 准备工作

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

1. 模型设置

  1. discount <- 0.95 # 贴现因子 beta^t
  2. gamma <- 1
  3. if (gamma==1){u.fmla <- expression(log(c))
  4. }else{
  5. u.fmla <- expression(c^(1-gamma)/(1-gamma))} # 效用函数表达式
  6. u <- deriv(u.fmla, "c", function(c) {}) # u(c)
  7. alpha <- 0.5
  8. f.fmla <- expression(z*k^alpha) # 生产函数表达式
  9. f <- deriv(f.fmla, c("k","z"), function(k,z) {}) # y=f(k,z)
  10. delta <- 1 # 折旧率
  11. kprime.fmla =expression(z*k^alpha + k*(1-delta) - c) # k_t+1表达式
  12. kprime <- deriv(kprime.fmla, c("k","z","c"), # k_t+1(k,z,c)
  13. function(k,z,c) {})
  14. z=zss=1
  15. kss <- ((1/discount+delta-1)/alpha)^(1/(alpha-1)) # 稳态
  16. css <- f(kss,zss) - delta*kss
  17. NK=100
  18. k.grid <- seq(0.01,3*kss,length.out=NK) # 设置k的格点

2. 价值函数迭代

  1. library(parallel)
  2. v.change <- 100; tol <- 1e-4 # 迭代设置
  3. iterations <- 0; maxiter=100
  4. c_min=kprime_min=0.01
  5. # 在每个格点上给出初始猜测,v0(k_i),比较好的初始猜测将大大提高迭代效率
  6. #v.grid <- rep(0,NK)
  7. v.grid <- rep(u(css)/(1-discount),NK)
  8. v0 <- approxfun(k.grid, v.grid, rule=2, method="linear") # 在格点之间用内插方法得到连续的v0(k)
  9. v.app <- list() # 用list存放v的迭代数据
  10. v.app[[1]] <- v0
  11. while(v.change > tol || iterations<maxiter) {
  12. bellman <- mcmapply(function(k) { # v(k_i)
  13. c_max=f(k,zss) + (1-delta)*k - 0.01
  14. optimize(function(c) { # = max_c (u + beta*v(kprime))
  15. u(c) + discount*v0(kprime(k,zss,c)) },
  16. interval=c(c_min, c_max),
  17. maximum=TRUE
  18. )
  19. }, k.grid)
  20. v.g.new <- unlist(bellman["objective",]) # v1(k_i)
  21. v.change <- max(abs(v.g.new-v.grid)) # 判断是否收敛
  22. iterations <- iterations+1
  23. v.grid <- v.g.new # 更新v0(k_i)=v1(k_i)
  24. v0 <- approxfun(k.grid,v.grid,rule=2,method="linear") # 拟合v0(k)
  25. v.app[[iterations+1]] <- v0 # 存放
  26. print(sprintf("After %d iterations v.change=%.2g\n",iterations,v.change))
  27. }
  28. c.grid=unlist(bellman["maximum",]) # 最优策略c(k_i)
  29. data<- data.frame(k=k.grid, c=c.grid, k.new=kprime(k.grid,zss,c.grid), f=f(k.grid,zss),
  30. v1=v.app[[1]](k.grid), v2=v.app[[2]](k.grid),
  31. v11=v.app[[11]](k.grid), v51=v.app[[51]](k.grid),
  32. v=v0(k.grid))

3. 作图

  1. library(ggplot2)
  2. value=ggplot(data,aes(x=k)) + # 价值函数
  3. geom_line(aes(y=v1,colour="V(0):Initial guess")) +
  4. geom_line(aes(y=v2,colour="V(1)"))+
  5. geom_line(aes(y=v11,colour="V(10)")) +
  6. geom_line(aes(y=v51,colour="V(50)")) +
  7. geom_line(aes(y=v,colour="V_converged")) +
  8. theme_minimal() +
  9. scale_colour_discrete(name="") +
  10. scale_y_continuous("v(k)") +
  11. scale_x_continuous("k")
  12. if (gamma==1 & delta==1){ # 特殊设定,用理论值检验V(k)
  13. theta0=1/(1-discount)*(alpha*discount*log(alpha*discount)/(1-alpha*discount)+log(1-alpha*discount))
  14. theta1=alpha/(1-alpha*discount)
  15. vcheck.grid=theta0+theta1*log(k.grid)
  16. data=data.frame(data,vcheck=vcheck.grid)
  17. head(data.frame(v.grid, vcheck.grid))
  18. tail(data.frame(v.grid, vcheck.grid))
  19. value <- value + geom_line(aes(y=data$vcheck,colour="V_theory"))
  20. }
  21. value

value.jpg-54.3kB

4. 策略函数

  1. policy <- ggplot(data,aes(x=k))+ # 策略函数
  2. geom_line(aes(y=c,colour="consumption")) +
  3. geom_line(aes(y=k.new,colour="next capital")) +
  4. geom_line(aes(y=f,colour="production")) +
  5. theme_minimal() +
  6. scale_colour_discrete(name="") +
  7. scale_y_continuous("consumption, next capital, or production") +
  8. scale_x_continuous("current capital")
  9. policy

policy.jpg-50.1kB

在线课程资料

  1. Thomas J. Sargent and John Stachurski: Quantitative Economics
  2. Paul Schrimpf: Mathematics for economics
添加新批注
在作者公开此批注前,只有你和作者可见。
回复批注