[关闭]
@fanxy 2017-05-16T01:01:40.000000Z 字数 2610 阅读 961

三、动态最优(I)拉格朗日方法

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


0. 准备工作

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

1. 生命周期储蓄问题:射击方法求解

  1. T=10;
  2. k0=kT=0
  3. w=1
  4. beta=0.95
  5. R=1/beta+c(0.1,0,-0.1) # 比较三种beta*R的情况,利率 R=1+r
  6. dim1=c("c","k") # 把结果存为3维数组x,行为变量名 c,k
  7. dim2=c(1:T) # 列为时间 1,2...T
  8. dim3=paste("beta*R=",round(beta*R,2)) # 高为 beta*R 的三种情况
  9. x=array(rep(NaN,2*T*length(R)),
  10. c(2,T,length(R)),
  11. dimnames=list(dim1,dim2,dim3))
  12. tol=0.00001; # 循环的参数设置
  13. maxloop=100;
  14. par(mfrow=c(length(R),1)) # 做图窗口,3张图
  15. for (case in 1:length(R)){ # 分情况搜索的循环
  16. Phi = matrix(c(beta*R[case],-1,0,R[case]),2,2)
  17. c_max=w/(1/beta-1) # 猜测区间。c_max > w 允许负债 k<0
  18. c_min=0.01 # c_min > 0
  19. error=2*tol; loop=0
  20. while (error>tol & loop<maxloop){
  21. c_guess=(c_min+c_max)/2 # c 的初始猜测
  22. c=rep(c_guess,T) # 初始化路径
  23. k=rep(k0,T)
  24. x[,,case]=rbind(c,k)
  25. for (t in 2:T){ # 更新路径
  26. x[,t,case]=Phi %*% matrix(x[,t-1,case]) +c(0,w)
  27. }
  28. if (x[2,T,case]<kT){ # 判断初始猜测是否过高
  29. c_max=c_guess # 调低上限
  30. } else
  31. {
  32. c_min=c_guess # 调高下限
  33. }
  34. error= abs(x[2,T,case]-kT) # 用于决断是否收敛
  35. loop=loop+1
  36. }
  37. matplot(t(x[,,case]),col=c("red","blue"),type="l", ylim=range(-2,2), main=dim3[case], ylab="")
  38. }

shooting.jpeg-61.2kB

2. 最优增长问题

2.1 设置参数

  1. T=20
  2. beta=0.95; gamma=1
  3. z=5; alpha=0.4; delta=1
  4. kstar=(alpha*z/(1/beta+delta-1))^(1/(1-alpha)) # 稳态
  5. cstar=((1/beta+delta-1)/alpha-delta)*kstar
  6. casetype=c(0.5,1.5)
  7. k0vector=kstar*casetype # 两种初始k0
  8. tol=0.0001; # 循环的参数设置
  9. maxloop=100;
  10. datax=data.frame() # 准备存放数据

2.2 搜索最优路径

  1. for (case in 1:length(k0vector)){
  2. k0=k0vector[case]
  3. c_max=z*k0^alpha+(1-delta)*k0-0.01 # 生产函数要求 k_t+1 > 0,构成c的上限
  4. c_min=0.01 # c > 0
  5. error=2*tol; loop=0
  6. while (error>tol & loop<maxloop){
  7. c_guess=(c_min+c_max)/2 # c 的初始猜测
  8. c=rep(c_guess,T) # 初始化路径
  9. k=rep(k0,T)
  10. for (t in 2:T){ # 更新路径
  11. c[t]=c[t-1]*(beta*(alpha*z*k[t-1]^(alpha-1)+1-delta))^(1/gamma)
  12. k[t]=z*k[t-1]^alpha+(1-delta)*k[t-1]-c[t]
  13. if (c[t]>c[t-1]& k[t]<k[t-1]){ # 判断初始猜测是否过高
  14. c_max=c_guess # 落入左上区域,调低上限
  15. break} # 中止循环
  16. else if (c[t]<c[t-1] & k[t]>k[t-1]){
  17. c_min=c_guess # 落入右下区域,调高下限
  18. break # 中止循环
  19. }
  20. }
  21. if (t==T & k[T]<kstar){ # 判断初始猜测是否过高
  22. c_max=c_guess # 调低上限
  23. error= abs(k[T]-kstar) # 用于决断是否收敛
  24. } else if (t==T & k[T]>kstar){
  25. c_min=c_guess # 调高下限
  26. error= abs(k[T]-kstar) # 用于决断是否收敛
  27. } # end of for()
  28. loop=loop+1
  29. } # end of while() for shooting
  30. x=data.frame(c,k,rep(cstar,T),rep(kstar,T))
  31. datax=rbind(datax,x)
  32. matplot(as.matrix(x), type="l",col=c("red","blue","pink","green"),main=paste("c_t and k_t with k0=",casetype[case],"*kstar"))
  33. legend("bottomright",c("c_t","k_t","cstar","kstar"),col=c("red","blue","pink","green"),lty=1)
  34. } # end of grid of k0vector

ck_0.5.jpg-62.2kB
ck_1.5.jpg-56.5kB

2.3 绘制相图

  1. datax$kss=z*datax$k^alpha-delta*datax$k
  2. library(ggplot2)
  3. ggplot(data=datax,legend=TRUE) +
  4. geom_vline(xintercept=kstar,colour="red") +
  5. geom_line(aes(x=k,y=kss),colour="blue") +
  6. geom_line(aes(x=k,y=c),colour="black") +
  7. ylab("c_t") + xlab("k_t") + xlim(k0vector) +
  8. theme(panel.background=element_rect(fill='transparent', color='gray'))

phase.jpg-44.4kB

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