Tip revision: 2878a3d
tok_Bellmans_TL.m
``````% tok_Bellmans_TL.m
% Function used to obtain normative thresholds in token lead space via
% dynamic programming for the tokens task from Barendregt et al., 2022.

function N_theta = tok_Bellmans_TL(Nt,t_d,R_c,R_i,c,tol)

% Define time and discretizations:
N = 0:Nt; t_i = 0.5+t_d*(Nt:-1:0);
dt = 0.2; t = 0:dt:(dt*Nt);

% Define overall likelihood discretization:
g = [];
for i = 0:Nt
for j = 0:i
g = [g prob_choice(i-j,j,Nt)];
end
end
g = sort(unique(g)); g_i = find(g==0.5);

% Initialize secant method:
rho = 0.1; k = 1;

% Pre-allocate value function V and maximal index of value function V_I:
V = NaN(length(g),length(N)); V_I = NaN(length(g),length(N));

% Calculate value function using backward induction:
[V(1,end),V_I(1,end)] = max([g(1)*R_c+(1-g(1))*R_i-t_i(end)*rho(k) (1-g(1))*R_c+g(1)*R_i-t_i(end)*rho(k)],[],2);
[V(end,end),V_I(end,end)] = max([g(end)*R_c+(1-g(end))*R_i-t_i(end)*rho(k) (1-g(end))*R_c+g(end)*R_i-t_i(end)*rho(k)],[],2);
for j = (length(N)-1):-1:1

% Define current likelihood discretization:
gt = [];
for i = 0:(j-1)
gt = [gt prob_choice((j-1)-i,i,Nt)];
end
gt = unique(gt);

% Compute future value assuming singular responses at likelihoods of 0
% and 1:
for i = 1:length(g)
if sum(gt==g(i))~=0
if g(i) == 0
[V(i,j),V_I(i,j)] = max([g(1)*R_c+(1-g(1))*R_i-t_i(j)*rho(k) (1-g(1))*R_c+g(1)*R_i-t_i(j)*rho(k) V(1,j+1)-(c(t(j))+rho(k))*dt],[],2);
elseif g(i) == 1
[V(i,j),V_I(i,j)] = max([g(end)*R_c+(1-g(end))*R_i-t_i(j)*rho(k) (1-g(end))*R_c+g(end)*R_i-t_i(j)*rho(k) V(end,j+1)-(c(t(j))+rho(k))*dt],[],2);
else
f = find(~isnan(V(:,j+1)));
ip = f(f>=i); im = f(f<=i); ip = ip(1); im = im(end);
[V(i,j),V_I(i,j)] = max([g(i)*R_c+(1-g(i))*R_i-t_i(j)*rho(k) (1-g(i))*R_c+g(i)*R_i-t_i(j)*rho(k) 0.5*(V(ip,j+1)+V(im,j+1))-(c(t(j))+rho(k))*dt],[],2);
end
end
end
end

% Store initial value to measure convergence and perform second
% initialization of secant method:
V_rho(k) = V(g_i,1); rho = [rho 0.9]; k = 2;

V = NaN(length(g),length(N)); V_I = NaN(length(g),length(N));
[V(1,end),V_I(1,end)] = max([g(1)*R_c+(1-g(1))*R_i-t_i(end)*rho(k) (1-g(1))*R_c+g(1)*R_i-t_i(end)*rho(k)],[],2);
[V(end,end),V_I(end,end)] = max([g(end)*R_c+(1-g(end))*R_i-t_i(end)*rho(k) (1-g(end))*R_c+g(end)*R_i-t_i(end)*rho(k)],[],2);
for j = (length(N)-1):-1:1
gt = [];
for i = 0:(j-1)
gt = [gt prob_choice((j-1)-i,i,Nt)];
end
gt = unique(gt);

for i = 1:length(g)
if sum(gt==g(i))~=0
if g(i) == 0
[V(i,j),V_I(i,j)] = max([g(1)*R_c+(1-g(1))*R_i-t_i(j)*rho(k) (1-g(1))*R_c+g(1)*R_i-t_i(j)*rho(k) V(1,j+1)-(c(t(j))+rho(k))*dt],[],2);
elseif g(i) == 1
[V(i,j),V_I(i,j)] = max([g(end)*R_c+(1-g(end))*R_i-t_i(j)*rho(k) (1-g(end))*R_c+g(end)*R_i-t_i(j)*rho(k) V(end,j+1)-(c(t(j))+rho(k))*dt],[],2);
else
f = find(~isnan(V(:,j+1)));
ip = f(f>=i); im = f(f<=i); ip = ip(1); im = im(end);
[V(i,j),V_I(i,j)] = max([g(i)*R_c+(1-g(i))*R_i-t_i(j)*rho(k) (1-g(i))*R_c+g(i)*R_i-t_i(j)*rho(k) 0.5*(V(ip,j+1)+V(im,j+1))-(c(t(j))+rho(k))*dt],[],2);
end
end
end
end

V_rho(k) = V(g_i,1);

% Continue interating using secant method until initial value has
% sufficiently converged:
while abs(V_rho(k)) > tol

% Update new reward rate using secant method:
k = k+1; rho(k) = rho(k-1)-V_rho(k-1)*(rho(k-1)-rho(k-2))/(V_rho(k-1)-V_rho(k-2));

V = NaN(length(g),length(N)); V_I = NaN(length(g),length(N));
[V(1,end),V_I(1,end)] = max([g(1)*R_c+(1-g(1))*R_i-t_i(end)*rho(k) (1-g(1))*R_c+g(1)*R_i-t_i(end)*rho(k)],[],2);
[V(end,end),V_I(end,end)] = max([g(end)*R_c+(1-g(end))*R_i-t_i(end)*rho(k) (1-g(end))*R_c+g(end)*R_i-t_i(end)*rho(k)],[],2);
for j = (length(N)-1):-1:1
gt = [];
for i = 0:(j-1)
gt = [gt prob_choice((j-1)-i,i,Nt)];
end
gt = unique(gt);

for i = 1:length(g)
if sum(gt==g(i))~=0
if g(i) == 0
[V(i,j),V_I(i,j)] = max([g(1)*R_c+(1-g(1))*R_i-t_i(j)*rho(k) (1-g(1))*R_c+g(1)*R_i-t_i(j)*rho(k) V(1,j+1)-(c(t(j))+rho(k))*dt],[],2);
elseif g(i) == 1
[V(i,j),V_I(i,j)] = max([g(end)*R_c+(1-g(end))*R_i-t_i(j)*rho(k) (1-g(end))*R_c+g(end)*R_i-t_i(j)*rho(k) V(end,j+1)-(c(t(j))+rho(k))*dt],[],2);
else
f = find(~isnan(V(:,j+1)));
ip = f(f>=i); im = f(f<=i); ip = ip(1); im = im(end);
[V(i,j),V_I(i,j)] = max([g(i)*R_c+(1-g(i))*R_i-t_i(j)*rho(k) (1-g(i))*R_c+g(i)*R_i-t_i(j)*rho(k) 0.5*(V(ip,j+1)+V(im,j+1))-(c(t(j))+rho(k))*dt],[],2);
end
end
end
end
V_rho(k) = V(g_i,1);
end

% Pre-allocate normative thresholds in token lead space:
N_theta = NaN(1,length(N));

% Construct normative thresholds in token lead space based off maximal
% index V_I:
for i =(length(N_theta)):-1:1
if ~isempty(find(V_I(:,i)==1,1))
if ~isempty(find(V_I(:,i)==3,1))
g_theta = g(find(V_I(:,i)==3,1,'last'));
for j = 0:(i-1)
if prob_choice(i-j-1,j,Nt) == g_theta
N_theta(i) = i-2*j-1;
end
end
else
N_theta(i) = 0;
end
else
N_theta(i) = 0.5*(Nt+1);
end
if N_theta(i) == 0
N_theta(i:end) = 0;
end
end``````