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

# 五、动态最优（III）动态规划方法

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

# 0. 准备工作

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

# 1. 模型设置

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

# 2. 价值函数迭代

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

# 3. 作图

library(ggplot2)value=ggplot(data,aes(x=k)) +                                # 价值函数  geom_line(aes(y=v1,colour="V(0)：Initial guess")) +  geom_line(aes(y=v2,colour="V(1)"))+  geom_line(aes(y=v11,colour="V(10)")) +  geom_line(aes(y=v51,colour="V(50)")) +    geom_line(aes(y=v,colour="V_converged")) +  theme_minimal() +  scale_colour_discrete(name="") +  scale_y_continuous("v(k)") +  scale_x_continuous("k") if (gamma==1 & delta==1){                                    # 特殊设定，用理论值检验V(k)   theta0=1/(1-discount)*(alpha*discount*log(alpha*discount)/(1-alpha*discount)+log(1-alpha*discount))  theta1=alpha/(1-alpha*discount)  vcheck.grid=theta0+theta1*log(k.grid)  data=data.frame(data,vcheck=vcheck.grid)  head(data.frame(v.grid, vcheck.grid))  tail(data.frame(v.grid, vcheck.grid))  value <- value + geom_line(aes(y=data\$vcheck,colour="V_theory")) }value

# 4. 策略函数

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

# 在线课程资料

1. Thomas J. Sargent and John Stachurski: Quantitative Economics
2. Paul Schrimpf: Mathematics for economics

• 私有
• 公开
• 删除