% Your code:
a_min = 1e-10;
a_max = 1;
na = 200;
a_grid = linspace(a_min,a_max,na)';
Nsim = 10000; % Desired value is larger, around 50000
r = 0.04;
nz = 14; % no. of points 2nd dim of pol_ap
TT = 80; % no. of points 3rd dim of pol_ap
pol_ap = rand(na,nz,TT); % fake data
z_ind_sim = randi(nz,Nsim,TT); % fake data
income_sim = rand(Nsim,TT); % fake data
%% Simulate Nsim agents for TT periods
a_sim = zeros(Nsim,TT+1);
c_sim = zeros(Nsim,TT);
% Set initial condition for assets
a_sim(:,1) = a_grid(1); % All agents start their life with a_grid(1)
tic
for j = 1:TT % outer loop over j
z_c = z_ind_sim(:,j); % dim: (Nsim,1)
for i=1:Nsim
a_sim(i,j+1) = interp1q(a_grid,pol_ap(:,z_c(i),j), a_sim(i,j));
end
% Here it is easy to vectorize the loop over i
c_sim(:,j) = (1+r)*a_sim(:,j)+income_sim(:,j)-a_sim(:,j+1);
end %END j
toc
% Interpolate once only to get the indices:
a_sim2 = zeros(Nsim,TT+1);
c_sim2 = zeros(Nsim,TT);
% Set initial condition for assets
a_sim2(:,1) = a_grid(1); % All agents start their life with a_grid(1)
tic
for j = 1:TT % outer loop over j
z_c = z_ind_sim(:, j); % dim: (Nsim,1)
index = interp1(a_grid, 1:size(pol_ap, 1), a_sim(:,j));
F = floor(index);
S = index - F;
for i = 1:Nsim
a_sim2(i, j+1) = pol_ap(F(i), z_c(i), j) * (1 - S(i)) + ...
pol_ap(F(i) + 1, z_c(i), j) * S(i);
end
% Here it is easy to vectorize the loop over i
c_sim2(:,j) = (1+r)*a_sim2(:,j)+income_sim(:,j)-a_sim(:,j+1);
end %END j
toc
max(abs(a_sim(:) - a_sim2(:)))
max(abs(c_sim(:) - c_sim2(:)))
