The followings are MatLab code and Julia code each.

```
%First, set up given data%
alpha=0.261;beta=0.99;delta=0.0176;n_star=1;
z_grid=[0.9526;0.9760;1;1.0246;1.0497];
Pi_z=[0.9149 0.0823 0.0028 0 0;...
0.0206 0.9163 0.0618 0.0014 0;...
0.0005 0.0412 0.9167 0.0412 0.0005;...
0 0.0014 0.0618 0.9163 0.0206;...
0 0 0.0028 0.0823 0.9149];
%Next, get the k_star and k_L and k_H%
k_star=(1/alpha*(1/beta-1+delta))^(1/(alpha-1));
k_L=k_star*0.85;k_H=1.15*k_star;
disp(round(k_star,4));disp(round(k_L,4));disp(round(k_H,4));
%Setting up an 499 points of kgrid%
k_grid=(k_L:((k_H-k_L)/498):k_H).';
disp(median(k_grid));
%%(b)%%
%First, set up the length variables%
n_k=length(k_grid);n_z=length(z_grid);
%Now initialize a V array%
V=zeros(n_k,n_z);
%Make settings for the iteration%
%First, make an array for each z_i%
k_prime_grid_mat_proto=repmat(k_grid,[1 n_z n_k]);%build prototype k' array with 3-d block array%
k_grid_mat=permute(k_prime_grid_mat_proto,[3 2 1]);%build k array with 3-d block array%;
k_prime_grid_mat=k_prime_grid_mat_proto.*(k_grid_mat.^alpha.*z_grid.'+k_grid_mat.*(1-delta)-k_prime_grid_mat_proto>=0);
k_prime_grid_mat(k_prime_grid_mat==0)=nan;
logC=log(k_grid_mat.^alpha.*z_grid.'+k_grid_mat.*(1-delta)-k_prime_grid_mat);
%build k' array that assigns nan to k' that is out of bound%
precision=0.01;distance=100;sz=size(k_prime_grid_mat);%precision and initial distance%
i=1;%iteration number%
%Start iteration%
while (distance > precision)
[TV,Tk_index]=max(logC+repmat(V*Pi_z.',[1 1 n_k]).*beta,[],'omitnan');%max within column omitting nan%
TV=reshape(permute(TV,[1,3,2]),[],size(TV,2));%change 1x5x499 array into 499x5 matrix%
distance=max(max(TV-V));
i=i+1;disp(i);disp(distance);
if distance<precision
break
else
V=TV;k_index=Tk_index;
end
end
%reviving k-rules from the k_index array%
k_prime_sol=k_prime_grid_mat_proto(Tk_index);
k_prime_sol_mat=reshape(permute(k_prime_sol,[1,3,2]),[],size(k_prime_sol,2));%change array into matrix%
k_prime_sol_previous=k_prime_grid_mat_proto(k_index);
k_prime_sol_previous_mat=reshape(permute(k_prime_sol,[1,3,2]),[],size(k_prime_sol_previous,2));%change array into matrix%
%matrix of previous capital to compare with solutions%
distance_k=max(max(abs(k_prime_sol_mat-k_prime_sol_previous_mat)));
%answers%
disp(i);disp(distance);disp(distance_k);
```

```
using Statistics
##(a)##
#First, set up the given data#
const alpha=0.261;const beta=0.99;const delta=0.0176;const n_star=1;
const z_grid=[0.9526;0.9760;1;1.0246;1.0497];
const Pi_z=[0.9149 0.0823 0.0028 0 0;
0.0206 0.9163 0.0618 0.0014 0;
0.0005 0.0412 0.9167 0.0412 0.0005;
0 0.0014 0.0618 0.9163 0.0206;
0 0 0.0028 0.0823 0.9149];
#Next, get the k_star and k_L and k_H#
const k_star=(1/alpha*(1/beta-1+delta))^(1/(alpha-1));
const k_L=k_star*0.85;const k_H=1.15*k_star;
println(round(k_star,digits=4));println(round(k_L,digits=4));println(round(k_H,digits=4));
#Setting up an 499 points of kgrid#
const k_grid=[k_L:((k_H-k_L)/498):k_H;];
println(median(k_grid))
##(b)##
#First, set up the length variables#
const n_k=length(k_grid);const n_z=length(z_grid);
#Make settings for the iteration#
#First, make an array for each z_i#
const k_prime_grid_mat_proto=repeat(k_grid,1, n_z, n_k);
const k_grid_mat=permutedims(k_prime_grid_mat_proto,(3,2,1));
k_prime_grid_mat=k_prime_grid_mat_proto.*(k_grid_mat.^alpha.*transpose(z_grid)+k_grid_mat.*(1-delta)-k_prime_grid_mat_proto.>0);
k_prime_grid_mat[k_prime_grid_mat.==0].=NaN;#build k' array that is out of bound#
logC=log.(k_grid_mat.^alpha.*transpose(z_grid)+k_grid_mat.*(1-delta)-k_prime_grid_mat);
logC[isnan.(logC)].=-Inf;
const precision=0.01;
const startdistance=100;#precision and initial distance#
#iteration number#
#start iteration#
function loop_over_global(distance, precision)
#Now initialize a V array#
global V=zeros(n_k,n_z);global i=1;
while distance>precision
global (TV,Tk_index)=findmax(logC::Array{Float64, 3}+repeat(V*transpose(Pi_z),1, 1, n_k).*beta,dims=1);
global TV=transpose(reshape(TV,5,499));
distance=maximum(TV-V);
println(i);println(distance);
if distance≤precision
break
else
global V=TV; global k_index=Tk_index; i=i+1;
end
end
return (TV,V,Tk_index,k_index,distance);
end
global (TV,V,Tk_index,k_index,distance)=loop_over_global(startdistance,precision);
#Reviving k-rules from the k_index array#
k_prime_sol_mat=transpose(reshape(k_prime_grid_mat_proto[Tk_index],5,499));
k_prime_sol_mat_previous=transpose(reshape(k_prime_grid_mat_proto[k_index],5,499));#matrix of previous capital to compare with solutions#
distance_k_prime_sol=maximum(abs.(k_prime_sol_mat-k_prime_sol_mat_previous));
##for information; how to convert a cartesian index array into matrix with only necessary info(first entry of the cartesian index)##
Tk_index=first.(Tuple.(Tk_index));Tk_index=transpose(reshape(Tk_index,5,499));
#answers#
println(i);println(distance);println(distance_k_prime_sol);
```

Thanks for reading this cumbersome codes. I wish I could get fluent in this fantastic language.