%本部分旨在将写出的贝尔曼方程中的值函数进行不断迭代，从而求解消费与资本的政策函数c(k);k'(k)
clear all
%设置参数
alpha=0.3;
beta=0.6;
delta=0.3;
epsi=0.00001;
%创建grid
k_min=0.04;
k_max=0.20;
K=linspace(k_min,k_max,1000);
%初始化
v_last=zeros(1,length(K));
v_new=ones(1,length(K));
v_order=zeros(1,length(K)); %在k' j循环的时候记录rhs,便于后续选出最大的rhs赋给v_new
C_op=zeros(1,length(K));
K_op=zeros(1,length(K));
KI_op=zeros(1,length(K));
%迭代过程
while max(abs(v_last-v_new))>epsi %判断误差
v_last=v_new; %更新值函数
for i =1:length(K)
for j =1:length(K)
if K(i)^alpha+(1-delta)*K(i)-K(j)>=0
rhs=log(K(i)^alpha+(1-delta)*K(i)-K(j))+beta*v_last(j);
v_order(j)=rhs; %把每一次j循环的rhs录入v_oder，便于后续排序。
end
end
[M,I] = max(v_order); % M是rhs的最大值，I是rhs最大时j的索引，以便于求出K_op&C_op。
v_new(i)=M;%选出最大的rhs赋给v_new
C_op(i)=K(i)^alpha+(1-delta)*K(i)-K(I); %使用I:rhs最大时的j取值
K_op(i)=K(I); %使用I:rhs最大时的j取值
KI_op(i)=I;%记录不同K取值时的rhs最大时的j，便于画图。
end
end
subplot(2,2,1)
plot(K,v_new)
title('值函数迭代')
xlabel('K')
x = linspace(0.04,0.20);
ylabel('v new')
subplot(2,2,2)
plot(K,K_op);
title('资本变量的政策函数图像')
xlabel('K')
x = linspace(0.04,0.20);
ylabel('K’')
subplot(2,2,3)
plot(K,C_op);
title('消费变量的政策函数图像')
xlabel('K')
x = linspace(0.04,0.20);
ylabel('C')
subplot(2,2,4)
plot(K,KI_op);
title('给定K，K’的j取值')
xlabel('K')
x = linspace(0.04,0.20);
ylabel('K’的j取值')
结果：
