I copied over a set of files from my workstation to use on my laptop when working in the labratory. The script runs fine when using the workstation but comes up with an error when used on my laptop. I didn’t write this code a former student at my university did but I need to alter his plots using the file. I’m new to MATLAB, and added the Raw Tension Data to the file path, and tried restarting the program/computer. I can’t seem to figure this out on my own. I’ve used other files the former student created with minimal issues on my laptop, this one is hanging me up.

Incorrect number or types of inputs or outputs for function resample.
Error in TensionPlots (line 7)
TPS_1 = resample(readmatrix(‘Raw Tension DataPS-1.csv’, ‘NumHeaderLines’,1),10,100);I copied over a set of files from my workstation to use on my laptop when working in the labratory. The script runs fine when using the workstation but comes up with an error when used on my laptop. I didn’t write this code a former student at my university did but I need to alter his plots using the file. I’m new to MATLAB, and added the Raw Tension Data to the file path, and tried restarting the program/computer. I can’t seem to figure this out on my own. I’ve used other files the former student created with minimal issues on my laptop, this one is hanging me up.

Incorrect number or types of inputs or outputs for function resample.
Error in TensionPlots (line 7)
TPS_1 = resample(readmatrix(‘Raw Tension DataPS-1.csv’, ‘NumHeaderLines’,1),10,100); I copied over a set of files from my workstation to use on my laptop when working in the labratory. The script runs fine when using the workstation but comes up with an error when used on my laptop. I didn’t write this code a former student at my university did but I need to alter his plots using the file. I’m new to MATLAB, and added the Raw Tension Data to the file path, and tried restarting the program/computer. I can’t seem to figure this out on my own. I’ve used other files the former student created with minimal issues on my laptop, this one is hanging me up.

Incorrect number or types of inputs or outputs for function resample.
Error in TensionPlots (line 7)
TPS_1 = resample(readmatrix(‘Raw Tension DataPS-1.csv’, ‘NumHeaderLines’,1),10,100); resample error MATLAB Answers — New Questions

​

Hi everyone,
I am new to coding so please bare with me as I am trying my best.
I am trying to plot and iterate over different data sets, for example, I want to plot data points 1-10, skip over 11 and 12, and then plot the next 10, 13-23 and so on and so forth. This is what I have so far. Please let me know where I have gone wrong. Thank you.

filename = "Mydatasheet.csv";
opts = detectImportOptions(‘filepathMydatasheet.csv’,’.csv’);
T = readtable("Mydatasheet.csv",opts);
dummy1 = table2array(:,1);
dummy2 = table2array(:,2);

for ii= 1: length(dummy,2)
plot(dummy1(1:10:230),dummy2(1:10:230));
endHi everyone,
I am new to coding so please bare with me as I am trying my best.
I am trying to plot and iterate over different data sets, for example, I want to plot data points 1-10, skip over 11 and 12, and then plot the next 10, 13-23 and so on and so forth. This is what I have so far. Please let me know where I have gone wrong. Thank you.

filename = "Mydatasheet.csv";
opts = detectImportOptions(‘filepathMydatasheet.csv’,’.csv’);
T = readtable("Mydatasheet.csv",opts);
dummy1 = table2array(:,1);
dummy2 = table2array(:,2);

for ii= 1: length(dummy,2)
plot(dummy1(1:10:230),dummy2(1:10:230));
end Hi everyone,
I am new to coding so please bare with me as I am trying my best.
I am trying to plot and iterate over different data sets, for example, I want to plot data points 1-10, skip over 11 and 12, and then plot the next 10, 13-23 and so on and so forth. This is what I have so far. Please let me know where I have gone wrong. Thank you.

filename = "Mydatasheet.csv";
opts = detectImportOptions(‘filepathMydatasheet.csv’,’.csv’);
T = readtable("Mydatasheet.csv",opts);
dummy1 = table2array(:,1);
dummy2 = table2array(:,2);

for ii= 1: length(dummy,2)
plot(dummy1(1:10:230),dummy2(1:10:230));
end plotting, for loop MATLAB Answers — New Questions

​

For my research, I measure a 2×2 scattering (S) matrix for a number of frequency points. Each component of the S matrix is complex, but varies smoothly and continuously. I am interested in the inner product of the two eigenvectors of the S matrix, but it is discontinuous due to the components of the eigenvectors behaving oddly. The issue isn’t that the label of which eigenvector is which has swapped, or a phase ambiguity of the eigenvectors. I define the eigenvectors as R1 = (A;B), R2 = (C;D) and choose A and D to be purely real. When the discontinuities in the inner product occur, a few things simultaneously occur. These include the label of eigenvalue 1 and 2 switching, a linear kink in the minima of A and D, a kink and switching of the real parts of B and C, and a kink or discontinuity in the imaginary parts of B and C. The code below and attached data show this behavior. Is there any way to fix this such that both the real and imaginary parts of the eigenvectors stay continuous, and the eigenvalues also be continuous with no switching?
%Experimental Data
load(‘Representative_Data.mat’)

%Eigensolver
for k=1:length(x_pts)
[Right_evec(:,:,k),eval(:,:,k)] = eig(smat(:,:,k));
Eval_1(k) = eval(1,1,k);
Eval_2(k) = eval(2,2,k);
end

%Get rid of eigenvector phase ambiguities
Phase_Component_1 = angle(Right_evec(1,1,:)); %A
Phase_Component_2 = angle(Right_evec(2,2,:)); %D
Right_evec(:,1,:) = Right_evec(:,1,:).*exp(-1i*Phase_Component_1);
Right_evec(:,2,:) = Right_evec(:,2,:).*exp(-1i*Phase_Component_2);

for k=1:length(x_pts)
Inner_Product(k) = Right_evec(:,1,k)’*Right_evec(:,2,k);
end

Marker_size = 20;

%Eigenvalues
figure(‘Position’,[200,100,1500,730]);
scatter(x_pts,real(Eval_1),Marker_size,’filled’); hold on;
scatter(x_pts,imag(Eval_1),Marker_size,’filled’);
scatter(x_pts,real(Eval_2),Marker_size,’filled’);
scatter(x_pts,imag(Eval_2),Marker_size,’filled’);
xlabel(‘Frequency (GHz)’)
ylabel(‘Eigenvalues’)
set(gca,’Fontsize’,26)
legend(‘Re lambda_1′,’Im lambda_1′,’Re lambda_2′,’Im lambda_2’)
grid on; box on;
%Eigenvector Components A and D
figure(‘Position’,[200,100,1500,730]);
scatter(x_pts,real(squeeze(Right_evec(1,1,:))),Marker_size,’filled’); hold on;
scatter(x_pts,real(squeeze(Right_evec(2,2,:))),Marker_size,’filled’);
%scatter(x_pts,imag(squeeze(Right_evec(1,1,:))),Marker_size,’filled’); %Imaginary part specified to be 0
%scatter(x_pts,imag(squeeze(Right_evec(2,2,:))),Marker_size,’filled’); %Imaginary part specified to be 0
xlabel(‘Frequency (GHz)’)
ylabel(‘Eigenvector components’)
set(gca,’Fontsize’,26)
legend(‘Re A’,’Re D’)
grid on; box on;
%Eigenvector Components B and C
figure(‘Position’,[200,100,1500,730]);
scatter(x_pts,real(squeeze(Right_evec(2,1,:))),Marker_size,’filled’); hold on;
scatter(x_pts,imag(squeeze(Right_evec(2,1,:))),Marker_size,’filled’);
scatter(x_pts,real(squeeze(Right_evec(1,2,:))),Marker_size,’filled’);
scatter(x_pts,imag(squeeze(Right_evec(1,2,:))),Marker_size,’filled’);
xlabel(‘Frequency (GHz)’)
ylabel(‘Eigenvector components’)
set(gca,’Fontsize’,26)
legend(‘Re B’,’Im B’,’Re C’,’Im C’)
grid on; box on;
%Inner Product
figure(‘Position’,[200,100,1500,730]);
scatter(x_pts,real(Inner_Product),Marker_size,’filled’); hold on;
%scatter(x_pts,imag(Inner_Product),Marker_size,’filled’); %Imaginary part is properly continuous
xlabel(‘Frequency (GHz)’)
ylabel(‘Re <R_1|R_2>’)
set(gca,’Fontsize’,26)
grid on; box on;For my research, I measure a 2×2 scattering (S) matrix for a number of frequency points. Each component of the S matrix is complex, but varies smoothly and continuously. I am interested in the inner product of the two eigenvectors of the S matrix, but it is discontinuous due to the components of the eigenvectors behaving oddly. The issue isn’t that the label of which eigenvector is which has swapped, or a phase ambiguity of the eigenvectors. I define the eigenvectors as R1 = (A;B), R2 = (C;D) and choose A and D to be purely real. When the discontinuities in the inner product occur, a few things simultaneously occur. These include the label of eigenvalue 1 and 2 switching, a linear kink in the minima of A and D, a kink and switching of the real parts of B and C, and a kink or discontinuity in the imaginary parts of B and C. The code below and attached data show this behavior. Is there any way to fix this such that both the real and imaginary parts of the eigenvectors stay continuous, and the eigenvalues also be continuous with no switching?
%Experimental Data
load(‘Representative_Data.mat’)

%Eigensolver
for k=1:length(x_pts)
[Right_evec(:,:,k),eval(:,:,k)] = eig(smat(:,:,k));
Eval_1(k) = eval(1,1,k);
Eval_2(k) = eval(2,2,k);
end

%Get rid of eigenvector phase ambiguities
Phase_Component_1 = angle(Right_evec(1,1,:)); %A
Phase_Component_2 = angle(Right_evec(2,2,:)); %D
Right_evec(:,1,:) = Right_evec(:,1,:).*exp(-1i*Phase_Component_1);
Right_evec(:,2,:) = Right_evec(:,2,:).*exp(-1i*Phase_Component_2);

for k=1:length(x_pts)
Inner_Product(k) = Right_evec(:,1,k)’*Right_evec(:,2,k);
end

Marker_size = 20;

%Eigenvalues
figure(‘Position’,[200,100,1500,730]);
scatter(x_pts,real(Eval_1),Marker_size,’filled’); hold on;
scatter(x_pts,imag(Eval_1),Marker_size,’filled’);
scatter(x_pts,real(Eval_2),Marker_size,’filled’);
scatter(x_pts,imag(Eval_2),Marker_size,’filled’);
xlabel(‘Frequency (GHz)’)
ylabel(‘Eigenvalues’)
set(gca,’Fontsize’,26)
legend(‘Re lambda_1′,’Im lambda_1′,’Re lambda_2′,’Im lambda_2’)
grid on; box on;
%Eigenvector Components A and D
figure(‘Position’,[200,100,1500,730]);
scatter(x_pts,real(squeeze(Right_evec(1,1,:))),Marker_size,’filled’); hold on;
scatter(x_pts,real(squeeze(Right_evec(2,2,:))),Marker_size,’filled’);
%scatter(x_pts,imag(squeeze(Right_evec(1,1,:))),Marker_size,’filled’); %Imaginary part specified to be 0
%scatter(x_pts,imag(squeeze(Right_evec(2,2,:))),Marker_size,’filled’); %Imaginary part specified to be 0
xlabel(‘Frequency (GHz)’)
ylabel(‘Eigenvector components’)
set(gca,’Fontsize’,26)
legend(‘Re A’,’Re D’)
grid on; box on;
%Eigenvector Components B and C
figure(‘Position’,[200,100,1500,730]);
scatter(x_pts,real(squeeze(Right_evec(2,1,:))),Marker_size,’filled’); hold on;
scatter(x_pts,imag(squeeze(Right_evec(2,1,:))),Marker_size,’filled’);
scatter(x_pts,real(squeeze(Right_evec(1,2,:))),Marker_size,’filled’);
scatter(x_pts,imag(squeeze(Right_evec(1,2,:))),Marker_size,’filled’);
xlabel(‘Frequency (GHz)’)
ylabel(‘Eigenvector components’)
set(gca,’Fontsize’,26)
legend(‘Re B’,’Im B’,’Re C’,’Im C’)
grid on; box on;
%Inner Product
figure(‘Position’,[200,100,1500,730]);
scatter(x_pts,real(Inner_Product),Marker_size,’filled’); hold on;
%scatter(x_pts,imag(Inner_Product),Marker_size,’filled’); %Imaginary part is properly continuous
xlabel(‘Frequency (GHz)’)
ylabel(‘Re <R_1|R_2>’)
set(gca,’Fontsize’,26)
grid on; box on; For my research, I measure a 2×2 scattering (S) matrix for a number of frequency points. Each component of the S matrix is complex, but varies smoothly and continuously. I am interested in the inner product of the two eigenvectors of the S matrix, but it is discontinuous due to the components of the eigenvectors behaving oddly. The issue isn’t that the label of which eigenvector is which has swapped, or a phase ambiguity of the eigenvectors. I define the eigenvectors as R1 = (A;B), R2 = (C;D) and choose A and D to be purely real. When the discontinuities in the inner product occur, a few things simultaneously occur. These include the label of eigenvalue 1 and 2 switching, a linear kink in the minima of A and D, a kink and switching of the real parts of B and C, and a kink or discontinuity in the imaginary parts of B and C. The code below and attached data show this behavior. Is there any way to fix this such that both the real and imaginary parts of the eigenvectors stay continuous, and the eigenvalues also be continuous with no switching?
%Experimental Data
load(‘Representative_Data.mat’)

%Eigensolver
for k=1:length(x_pts)
[Right_evec(:,:,k),eval(:,:,k)] = eig(smat(:,:,k));
Eval_1(k) = eval(1,1,k);
Eval_2(k) = eval(2,2,k);
end

%Get rid of eigenvector phase ambiguities
Phase_Component_1 = angle(Right_evec(1,1,:)); %A
Phase_Component_2 = angle(Right_evec(2,2,:)); %D
Right_evec(:,1,:) = Right_evec(:,1,:).*exp(-1i*Phase_Component_1);
Right_evec(:,2,:) = Right_evec(:,2,:).*exp(-1i*Phase_Component_2);

for k=1:length(x_pts)
Inner_Product(k) = Right_evec(:,1,k)’*Right_evec(:,2,k);
end

Marker_size = 20;

%Eigenvalues
figure(‘Position’,[200,100,1500,730]);
scatter(x_pts,real(Eval_1),Marker_size,’filled’); hold on;
scatter(x_pts,imag(Eval_1),Marker_size,’filled’);
scatter(x_pts,real(Eval_2),Marker_size,’filled’);
scatter(x_pts,imag(Eval_2),Marker_size,’filled’);
xlabel(‘Frequency (GHz)’)
ylabel(‘Eigenvalues’)
set(gca,’Fontsize’,26)
legend(‘Re lambda_1′,’Im lambda_1′,’Re lambda_2′,’Im lambda_2’)
grid on; box on;
%Eigenvector Components A and D
figure(‘Position’,[200,100,1500,730]);
scatter(x_pts,real(squeeze(Right_evec(1,1,:))),Marker_size,’filled’); hold on;
scatter(x_pts,real(squeeze(Right_evec(2,2,:))),Marker_size,’filled’);
%scatter(x_pts,imag(squeeze(Right_evec(1,1,:))),Marker_size,’filled’); %Imaginary part specified to be 0
%scatter(x_pts,imag(squeeze(Right_evec(2,2,:))),Marker_size,’filled’); %Imaginary part specified to be 0
xlabel(‘Frequency (GHz)’)
ylabel(‘Eigenvector components’)
set(gca,’Fontsize’,26)
legend(‘Re A’,’Re D’)
grid on; box on;
%Eigenvector Components B and C
figure(‘Position’,[200,100,1500,730]);
scatter(x_pts,real(squeeze(Right_evec(2,1,:))),Marker_size,’filled’); hold on;
scatter(x_pts,imag(squeeze(Right_evec(2,1,:))),Marker_size,’filled’);
scatter(x_pts,real(squeeze(Right_evec(1,2,:))),Marker_size,’filled’);
scatter(x_pts,imag(squeeze(Right_evec(1,2,:))),Marker_size,’filled’);
xlabel(‘Frequency (GHz)’)
ylabel(‘Eigenvector components’)
set(gca,’Fontsize’,26)
legend(‘Re B’,’Im B’,’Re C’,’Im C’)
grid on; box on;
%Inner Product
figure(‘Position’,[200,100,1500,730]);
scatter(x_pts,real(Inner_Product),Marker_size,’filled’); hold on;
%scatter(x_pts,imag(Inner_Product),Marker_size,’filled’); %Imaginary part is properly continuous
xlabel(‘Frequency (GHz)’)
ylabel(‘Re <R_1|R_2>’)
set(gca,’Fontsize’,26)
grid on; box on; eigenvectors MATLAB Answers — New Questions

​

Hi, I want to convert a transfer function to controllable and observable canonical form for the
num = [4];
den = [1 0.8 4];
Gp = tf (num , den)

Gp =

4
—————
s^2 + 0.8 s + 4Hi, I want to convert a transfer function to controllable and observable canonical form for the
num = [4];
den = [1 0.8 4];
Gp = tf (num , den)

Gp =

4
—————
s^2 + 0.8 s + 4 Hi, I want to convert a transfer function to controllable and observable canonical form for the
num = [4];
den = [1 0.8 4];
Gp = tf (num , den)

Gp =

4
—————
s^2 + 0.8 s + 4 transfer function, observability and controlability MATLAB Answers — New Questions

​

Hello! I am writing a script that makes the user draw a freehand object, then prompts them for a number, and finally attaches that number onto the produced ROI as a visible label. The currently generated label is nice but I was hoping to increase the fontsize to make it more visible on a screen. Is there a good way to do this? Thank you very much!
Below, I have included an executable sample of what I’m working with at the moment.
% Initiate drawfreehand function
Drawing = drawfreehand("Color","r", "LabelAlpha",1);

% Prompt User
prompt = "What value?n";
Value = input(prompt);

% Set the freehand object’s label
Drawing.Label = num2str(Value);Hello! I am writing a script that makes the user draw a freehand object, then prompts them for a number, and finally attaches that number onto the produced ROI as a visible label. The currently generated label is nice but I was hoping to increase the fontsize to make it more visible on a screen. Is there a good way to do this? Thank you very much!
Below, I have included an executable sample of what I’m working with at the moment.
% Initiate drawfreehand function
Drawing = drawfreehand("Color","r", "LabelAlpha",1);

% Prompt User
prompt = "What value?n";
Value = input(prompt);

% Set the freehand object’s label
Drawing.Label = num2str(Value); Hello! I am writing a script that makes the user draw a freehand object, then prompts them for a number, and finally attaches that number onto the produced ROI as a visible label. The currently generated label is nice but I was hoping to increase the fontsize to make it more visible on a screen. Is there a good way to do this? Thank you very much!
Below, I have included an executable sample of what I’m working with at the moment.
% Initiate drawfreehand function
Drawing = drawfreehand("Color","r", "LabelAlpha",1);

% Prompt User
prompt = "What value?n";
Value = input(prompt);

% Set the freehand object’s label
Drawing.Label = num2str(Value); drawfreehand, roi, image analysis MATLAB Answers — New Questions

​

I am trying to solve the problem using stackelberg sequential equilibrium method –

Understanding the Approach:
In a Stackelberg game, we have two players: the leader and the follower. The leader chooses its strategy first, knowing that the follower will then optimize their strategy in response to the leader’s choice. This sequential decision-making process requires us to iteratively solve for the best responses of each player until we converge to a Nash equilibrium.

-Follower’s Best Response (t given k):

The follower’s strategy t is a function of the leader’s strategy k. We need to find the value of t that maximizes the follower’s utility given the leader’s choice of k. We start by initial guess of k and compute corresponding optimal t.

Leader’s Best Response (k inputing response of t):

Once we have the optimal t, the leader then chooses its strategy k to maximize its utility, knowing the follower will choose t optimally in response.

We will do this iteration till initial guess converge to optimal k. Following code encapsulate my idea –

% Parameters
m_values = linspace(0, 1, 100); % Range of m values
tolerance = 1e-6; % Convergence tolerance
max_iter = 1000; % Maximum number of iterations

% Initialize arrays to store results
k_solutions = zeros(size(m_values));
t_solutions = zeros(size(m_values));

% Function to compute follower’s best response t given k and m
compute_t = @(k, m) fminbnd(@(t) abs(sqrt(k) * ((m./t + 2) / 3) – sqrt(1 – k) * (1 / 3) * (2 + m./t + (2 * m + t) ./ ((m – t).^2 – 2 * t))), 0, 1);

% Function to compute leader’s objective function given k, t, and m
compute_obj = @(k, t, m) -(sqrt(1 – k) * (t^6 * (2 * m * (3 * m + 2) + 4) – t^5 * (m * (m * (5 * m + 6) + 11) + 12) + m^6 + 2 * t^7 – t^8 – t^4 * (m * (3 * m * (m * (3 * m + 2) – 3) – 16) – 9) – 2 * m^3 * t^2 * (3 * m^3 + 2 * m + 3) + m^2 * t^3 * (m * (3 * m * (5 * m + 4) + 1) – 4)) + sqrt(k) * (t^4 * (m * (9 * m * (m * (m + 4) + 5) + 28) + 9) – m^4 * (3 * m – 2) – t^6 * (6 * m * (m + 2) + 10) + t^5 * (4 * m^3 + 12 * m + 8) + t^8 + m * t^2 * (m * (m * (4 * m^3 + 3 * m + 9) + 6) – 8) + 2 * m^2 * t * (3 * m – 2) – 2 * m * t^3 * (m * (m * (6 * m * (m + 2) + 7) + 6) – 2))) / (18 * t^3 * (2 * t – (m – t)^2));

% Main loop to solve for each m
for i = 1:length(m_values)
m = m_values(i);
k_guess = 0.75; % Initial guess for k
k_opt = k_guess;
t_opt = compute_t(k_opt, m);

% Iterative optimization
for iter = 1:max_iter
% Compute optimal t given k
t_opt = compute_t(k_opt, m);

% Optimize k given t
k_opt_new = fminbnd(@(k) compute_obj(k, t_opt, m), 0.5, 1); % Use fminbnd for bounded optimization

% Check for convergence
if abs(k_opt_new – k_opt) < tolerance
k_opt = k_opt_new;
break;
end
k_opt = k_opt_new;
end

% Store solutions
k_solutions(i) = k_opt;
t_solutions(i) = t_opt;
end

% Plot results
figure;
plot(m_values, k_solutions, ‘b-‘, ‘LineWidth’, 1.5);
hold on;
plot(m_values, t_solutions, ‘r–‘, ‘LineWidth’, 1.5);
xlabel(‘m’);
ylabel(‘Value’);
legend(‘Optimal k’, ‘Optimal t’);
title(‘Stackelberg Equilibrium Solutions’);

% Display final values
fprintf(‘Final solutions:n’);
fprintf(‘m t k t tn’);
for i = 1:length(m_values)
fprintf(‘%.4f t %.4f t %.4fn’, m_values(i), k_solutions(i), t_solutions(i));
end

I could have used calculus… But since expression of t cannot be explicitly expressed as a function of k, substituting t into k and then differentiating wrt k (where, t is also function k) would give an untractably long expression.

Is the Approach I used correct? Can someone please suggest a better approach? Anyone please helpI am trying to solve the problem using stackelberg sequential equilibrium method –

Understanding the Approach:
In a Stackelberg game, we have two players: the leader and the follower. The leader chooses its strategy first, knowing that the follower will then optimize their strategy in response to the leader’s choice. This sequential decision-making process requires us to iteratively solve for the best responses of each player until we converge to a Nash equilibrium.

-Follower’s Best Response (t given k):

The follower’s strategy t is a function of the leader’s strategy k. We need to find the value of t that maximizes the follower’s utility given the leader’s choice of k. We start by initial guess of k and compute corresponding optimal t.

Leader’s Best Response (k inputing response of t):

Once we have the optimal t, the leader then chooses its strategy k to maximize its utility, knowing the follower will choose t optimally in response.

We will do this iteration till initial guess converge to optimal k. Following code encapsulate my idea –

% Parameters
m_values = linspace(0, 1, 100); % Range of m values
tolerance = 1e-6; % Convergence tolerance
max_iter = 1000; % Maximum number of iterations

% Initialize arrays to store results
k_solutions = zeros(size(m_values));
t_solutions = zeros(size(m_values));

% Function to compute follower’s best response t given k and m
compute_t = @(k, m) fminbnd(@(t) abs(sqrt(k) * ((m./t + 2) / 3) – sqrt(1 – k) * (1 / 3) * (2 + m./t + (2 * m + t) ./ ((m – t).^2 – 2 * t))), 0, 1);

% Function to compute leader’s objective function given k, t, and m
compute_obj = @(k, t, m) -(sqrt(1 – k) * (t^6 * (2 * m * (3 * m + 2) + 4) – t^5 * (m * (m * (5 * m + 6) + 11) + 12) + m^6 + 2 * t^7 – t^8 – t^4 * (m * (3 * m * (m * (3 * m + 2) – 3) – 16) – 9) – 2 * m^3 * t^2 * (3 * m^3 + 2 * m + 3) + m^2 * t^3 * (m * (3 * m * (5 * m + 4) + 1) – 4)) + sqrt(k) * (t^4 * (m * (9 * m * (m * (m + 4) + 5) + 28) + 9) – m^4 * (3 * m – 2) – t^6 * (6 * m * (m + 2) + 10) + t^5 * (4 * m^3 + 12 * m + 8) + t^8 + m * t^2 * (m * (m * (4 * m^3 + 3 * m + 9) + 6) – 8) + 2 * m^2 * t * (3 * m – 2) – 2 * m * t^3 * (m * (m * (6 * m * (m + 2) + 7) + 6) – 2))) / (18 * t^3 * (2 * t – (m – t)^2));

% Main loop to solve for each m
for i = 1:length(m_values)
m = m_values(i);
k_guess = 0.75; % Initial guess for k
k_opt = k_guess;
t_opt = compute_t(k_opt, m);

% Iterative optimization
for iter = 1:max_iter
% Compute optimal t given k
t_opt = compute_t(k_opt, m);

% Optimize k given t
k_opt_new = fminbnd(@(k) compute_obj(k, t_opt, m), 0.5, 1); % Use fminbnd for bounded optimization

% Check for convergence
if abs(k_opt_new – k_opt) < tolerance
k_opt = k_opt_new;
break;
end
k_opt = k_opt_new;
end

% Store solutions
k_solutions(i) = k_opt;
t_solutions(i) = t_opt;
end

% Plot results
figure;
plot(m_values, k_solutions, ‘b-‘, ‘LineWidth’, 1.5);
hold on;
plot(m_values, t_solutions, ‘r–‘, ‘LineWidth’, 1.5);
xlabel(‘m’);
ylabel(‘Value’);
legend(‘Optimal k’, ‘Optimal t’);
title(‘Stackelberg Equilibrium Solutions’);

% Display final values
fprintf(‘Final solutions:n’);
fprintf(‘m t k t tn’);
for i = 1:length(m_values)
fprintf(‘%.4f t %.4f t %.4fn’, m_values(i), k_solutions(i), t_solutions(i));
end

I could have used calculus… But since expression of t cannot be explicitly expressed as a function of k, substituting t into k and then differentiating wrt k (where, t is also function k) would give an untractably long expression.

Is the Approach I used correct? Can someone please suggest a better approach? Anyone please help I am trying to solve the problem using stackelberg sequential equilibrium method –

Understanding the Approach:
In a Stackelberg game, we have two players: the leader and the follower. The leader chooses its strategy first, knowing that the follower will then optimize their strategy in response to the leader’s choice. This sequential decision-making process requires us to iteratively solve for the best responses of each player until we converge to a Nash equilibrium.

-Follower’s Best Response (t given k):

The follower’s strategy t is a function of the leader’s strategy k. We need to find the value of t that maximizes the follower’s utility given the leader’s choice of k. We start by initial guess of k and compute corresponding optimal t.

Leader’s Best Response (k inputing response of t):

Once we have the optimal t, the leader then chooses its strategy k to maximize its utility, knowing the follower will choose t optimally in response.

We will do this iteration till initial guess converge to optimal k. Following code encapsulate my idea –

% Parameters
m_values = linspace(0, 1, 100); % Range of m values
tolerance = 1e-6; % Convergence tolerance
max_iter = 1000; % Maximum number of iterations

% Initialize arrays to store results
k_solutions = zeros(size(m_values));
t_solutions = zeros(size(m_values));

% Function to compute follower’s best response t given k and m
compute_t = @(k, m) fminbnd(@(t) abs(sqrt(k) * ((m./t + 2) / 3) – sqrt(1 – k) * (1 / 3) * (2 + m./t + (2 * m + t) ./ ((m – t).^2 – 2 * t))), 0, 1);

% Function to compute leader’s objective function given k, t, and m
compute_obj = @(k, t, m) -(sqrt(1 – k) * (t^6 * (2 * m * (3 * m + 2) + 4) – t^5 * (m * (m * (5 * m + 6) + 11) + 12) + m^6 + 2 * t^7 – t^8 – t^4 * (m * (3 * m * (m * (3 * m + 2) – 3) – 16) – 9) – 2 * m^3 * t^2 * (3 * m^3 + 2 * m + 3) + m^2 * t^3 * (m * (3 * m * (5 * m + 4) + 1) – 4)) + sqrt(k) * (t^4 * (m * (9 * m * (m * (m + 4) + 5) + 28) + 9) – m^4 * (3 * m – 2) – t^6 * (6 * m * (m + 2) + 10) + t^5 * (4 * m^3 + 12 * m + 8) + t^8 + m * t^2 * (m * (m * (4 * m^3 + 3 * m + 9) + 6) – 8) + 2 * m^2 * t * (3 * m – 2) – 2 * m * t^3 * (m * (m * (6 * m * (m + 2) + 7) + 6) – 2))) / (18 * t^3 * (2 * t – (m – t)^2));

% Main loop to solve for each m
for i = 1:length(m_values)
m = m_values(i);
k_guess = 0.75; % Initial guess for k
k_opt = k_guess;
t_opt = compute_t(k_opt, m);

% Iterative optimization
for iter = 1:max_iter
% Compute optimal t given k
t_opt = compute_t(k_opt, m);

% Optimize k given t
k_opt_new = fminbnd(@(k) compute_obj(k, t_opt, m), 0.5, 1); % Use fminbnd for bounded optimization

% Check for convergence
if abs(k_opt_new – k_opt) < tolerance
k_opt = k_opt_new;
break;
end
k_opt = k_opt_new;
end

% Store solutions
k_solutions(i) = k_opt;
t_solutions(i) = t_opt;
end

% Plot results
figure;
plot(m_values, k_solutions, ‘b-‘, ‘LineWidth’, 1.5);
hold on;
plot(m_values, t_solutions, ‘r–‘, ‘LineWidth’, 1.5);
xlabel(‘m’);
ylabel(‘Value’);
legend(‘Optimal k’, ‘Optimal t’);
title(‘Stackelberg Equilibrium Solutions’);

% Display final values
fprintf(‘Final solutions:n’);
fprintf(‘m t k t tn’);
for i = 1:length(m_values)
fprintf(‘%.4f t %.4f t %.4fn’, m_values(i), k_solutions(i), t_solutions(i));
end

I could have used calculus… But since expression of t cannot be explicitly expressed as a function of k, substituting t into k and then differentiating wrt k (where, t is also function k) would give an untractably long expression.

Is the Approach I used correct? Can someone please suggest a better approach? Anyone please help iteration, optimization, fmincon MATLAB Answers — New Questions

​

hi, i am new to MATLAB and i need help to use it for my assignment
so i was given a set of data and asked to plot it onto a graph. but i am not sure how to.

do i enter the values into a table or array? can anyone help me? any help would be deeply appreciated.hi, i am new to MATLAB and i need help to use it for my assignment
so i was given a set of data and asked to plot it onto a graph. but i am not sure how to.

do i enter the values into a table or array? can anyone help me? any help would be deeply appreciated. hi, i am new to MATLAB and i need help to use it for my assignment
so i was given a set of data and asked to plot it onto a graph. but i am not sure how to.

do i enter the values into a table or array? can anyone help me? any help would be deeply appreciated. homework, wavelength MATLAB Answers — New Questions

​

Hi there,
I think I have used this procedure for the following loop many times, but now it isn’t working:

h = 3
z = 4
y = 0
x = 0

for i = 1:3
z(i+1) = z(i) + (2*y(i) + 8*x(i)*(9 – x(i)))*h
y(i+1) = y(i) + z(i)*h
z(i) = z(i+1);
y(i) = y(i+1);
x(i) = x(i) + h
end

I keep getting an error saying Index must not exceed 1.

I don’t know why this is happening.

Can someone help please?Hi there,
I think I have used this procedure for the following loop many times, but now it isn’t working:

h = 3
z = 4
y = 0
x = 0

for i = 1:3
z(i+1) = z(i) + (2*y(i) + 8*x(i)*(9 – x(i)))*h
y(i+1) = y(i) + z(i)*h
z(i) = z(i+1);
y(i) = y(i+1);
x(i) = x(i) + h
end

I keep getting an error saying Index must not exceed 1.

I don’t know why this is happening.

Can someone help please? Hi there,
I think I have used this procedure for the following loop many times, but now it isn’t working:

h = 3
z = 4
y = 0
x = 0

for i = 1:3
z(i+1) = z(i) + (2*y(i) + 8*x(i)*(9 – x(i)))*h
y(i+1) = y(i) + z(i)*h
z(i) = z(i+1);
y(i) = y(i+1);
x(i) = x(i) + h
end

I keep getting an error saying Index must not exceed 1.

I don’t know why this is happening.

Can someone help please? for loops MATLAB Answers — New Questions

​

For the past month or so, when I update a file online, saving it in ForDriveConnector as I was doing before, I need to log out and back in for Matlab Online to see the file. Why is this so?For the past month or so, when I update a file online, saving it in ForDriveConnector as I was doing before, I need to log out and back in for Matlab Online to see the file. Why is this so? For the past month or so, when I update a file online, saving it in ForDriveConnector as I was doing before, I need to log out and back in for Matlab Online to see the file. Why is this so? file-saved-online MATLAB Answers — New Questions

​

I have to solve a dynamic programming problem with a finite horizon and I am trying to vectorize as much as possible for speed. I attach here a MWE so that everybody can run it
I have run the code with the Matlab profiler and indetified two bottlenecks, marked with the comment line % THIS IS SLOW ACCORDING TO PROFILER
The first bottleneck is the call to the function ReturnFn that builds the n_a*n_a matrix RetMat. I have done this in a vectorized way (inside the function RetMat, I do not use loops)
The second bottleneck is the maximization of RetMat along the first dimension.
I’d be very grateful for any comment/suggestion!

P.S. I have an if condition inside the loops to check for bugs, I removed it now that I am confident about the code but the speed improvement is marginal.
clear,clc,close all
% STATE VARIABLES
% V(a,h,z,j)
% a: asset holdings
% h: human capital (female)
% z: labor productivity shocks (eta_m, eta_f are shocks, theta is permanent)
% j: age (from 1 to N_j)
% CHOICE VARIABLES
% d: Female labor supply (only extensive margin: either 0 or 1)
% a’: Next-period assets
% h’ is implied by (d,a) and consumption is implied by the budget
% constraint
% DYNAMIC PROGRAMMING PROBLEM
% V(a,h,z,j) = max_{d,a’} F(d,a’,a,h,z,j)+beta*s_j*E[V(a’,h’,z’,j+1)|z]
% subject to
% h’=G(d,h), law of motion for human capital
verbose = 1;
%% Define grids and grid sizes
N_j = 80;
n_a = 51;
n_h = 11;
n_z = 50;
n_d = 2;
a_grid = linspace(0,450,n_a)’;
h_grid = linspace(0,0.72,n_h)’;
z_grid = linspace(0.9,1.1,n_z)’;
d_grid = [0,1]’;
z_grid = repmat(z_grid,[1,3]);
pi_z = rand(n_z,n_z);
pi_z = pi_z./sum(pi_z,2);
aprime_val = a_grid; %(a’,1)
a_val = a_grid’; %(1,a)
%% Set parameters that do not depend on age
beta = 0.98;
r = 0.04;
w_m = 1;
w_f = 0.75;
crra = 2;
nu = 0.12;
Jr = 45;
xi_1 = 0.05312;
xi_2 = -0.00188;
del_h = 0.074;
h_l = 0;
p.eff_j = ones(N_j,1);
p.pchild_j = ones(N_j,1);
p.pen_j = ones(N_j,1);
p.nchild_j = ones(N_j,1);
p.s_j = ones(N_j,1);
p.age_j = (1:1:N_j)’;
% Initialize output arrays
V = zeros(n_a,n_h,n_z,N_j);
Policy = zeros(2,n_a,n_h,n_z,N_j);
tic
%% Solve problem in the last period
% Set age-dependent parameters
eff_j = p.eff_j(N_j);
pchild_j = p.pchild_j(N_j);
pen_j = p.pen_j(N_j);
nchild_j = p.nchild_j(N_j);
age_j = p.age_j(N_j);
s_j = p.s_j(N_j);
% V(a,h,z,N_j) = max_{d,a’}
V_d = zeros(n_a,n_h,n_z,n_d);
Pol_aprime_d = zeros(n_a,n_h,n_z,n_d);
for d_c = 1:n_d
d_val = d_grid(d_c);
for z_c = 1:n_z
eta_m_val = z_grid(z_c,1);
eta_f_val = z_grid(z_c,2);
theta_val = z_grid(z_c,3);
for h_c = 1:n_h
h_val = h_grid(h_c);
% RetMat is (a’,a)
RetMat = ReturnFn(d_val,aprime_val,a_val,h_val,eta_m_val,eta_f_val,theta_val,…
w_m,eff_j,w_f,pchild_j,pen_j,r,nchild_j,crra,nu,age_j,Jr);
[max_val,max_ind] = max(RetMat,[],1);
V_d(:,h_c,z_c,d_c) = max_val;
Pol_aprime_d(:,h_c,z_c,d_c) = max_ind;
end %end h
end %end z
end %end d
[V(:,:,:,N_j),d_max] = max(V_d,[],4);
Policy(1,:,:,:,N_j) = d_max; % Optimal d
for z_c=1:n_z
for h_c=1:n_h
for a_c = 1:n_a
d_star = d_max(a_c,h_c,z_c);
Policy(2,a_c,h_c,z_c,N_j) = Pol_aprime_d(a_c,h_c,z_c,d_star); % Optimal a’
end
end
end
%% Backward iteration over age
for j = N_j-1:-1:1
if verbose==1; fprintf(‘Age %d out of %d n’,j,N_j); end
V_next = V(:,:,:,j+1); %V(a’,h’,z’)
% Set age-dependent parameters
eff_j = p.eff_j(j);
pchild_j = p.pchild_j(j);
pen_j = p.pen_j(j);
nchild_j = p.nchild_j(j);
age_j = p.age_j(j);
s_j = p.s_j(j);
for z_c = 1:n_z
eta_m_val = z_grid(z_c,1);
eta_f_val = z_grid(z_c,2);
theta_val = z_grid(z_c,3);
% Compute EV(a’,h’), given z
% EV = zeros(n_a,n_h);
z_prob = pi_z(z_c,:)’;
% for zprime_c = 1:n_z
% EV = EV+V_next(:,:,zprime_c)*z_prob(zprime_c);
% end %end z’
EV = V_next.*shiftdim(z_prob,-2); %V(a’,h’,z’)*Prob(1,1,z’)
EV = sum(EV,3); %EV(a’,h’)
for d_c=1:n_d
d_val = d_grid(d_c);
for h_c = 1:n_h
h_val = h_grid(h_c);
hprime_val = f_HC_accum(d_val,h_val,age_j,xi_1,xi_2,del_h,h_l);
% Ret_mat is (a’,a)
% THIS IS SLOW ACCORDING TO PROFILER
Ret_mat = ReturnFn(d_val,aprime_val,a_val,h_val,eta_m_val,eta_f_val,theta_val,…
w_m,eff_j,w_f,pchild_j,pen_j,r,nchild_j,crra,nu,age_j,Jr);
%[ind_l,weight_l] = interp_toolkit(hprime_val,h_grid);
[ind_l,weight_l] = find_loc(h_grid,hprime_val);
if ind_l>length(h_grid)-1
error(‘ind_l out of bounds’)
end
EV_interp = EV(:,ind_l)*weight_l+EV(:,ind_l+1)*(1-weight_l);
RHS_mat = Ret_mat+beta*s_j*EV_interp;
% THIS IS SLOW ACCORDING TO PROFILER
[max_val,max_ind] = max(RHS_mat,[],1);
% max_val and max_ind are (1,a)
V_d(:,h_c,z_c,d_c) = max_val; %best V given d
Pol_aprime_d(:,h_c,z_c,d_c) = max_ind; %best a’ given d
end %end h
end %end d
end %end z
[V(:,:,:,j),d_max] = max(V_d,[],4);
Policy(1,:,:,:,j) = d_max; % Optimal d
for z_c=1:n_z
for h_c=1:n_h
for a_c = 1:n_a
d_star = d_max(a_c,h_c,z_c);
Policy(2,a_c,h_c,z_c,j) = Pol_aprime_d(a_c,h_c,z_c,d_star); % Optimal a’
end
end
end
end %end j
toc
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function F = ReturnFn(l_f,aprime,a,h_f,eta_m,eta_f,theta,w_m,eff_j,w_f,…
pchild_j,pen_j,r,nchild_j,crra,nu,agej,Jr)
% Calculate earnings (incl. child care costs) of men and women
y_m = w_m*eff_j*theta*eta_m;
y_f = w_f*l_f*(exp(h_f)*theta*eta_f – pchild_j);
% l_f can be either 0 or 1
% calculate available resources
cash = (1+r)*a + pen_j*(agej>=Jr) + (y_m + y_f)*(agej<Jr);
cons = cash-aprime;
%pos = cons>0;
F = (cons/(sqrt(2+nchild_j))).^(1-crra)/(1-crra) – nu*l_f;
F(cons<=0) = -inf;
end %end function "f_ReturnFn"
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [h_f_prime] = f_HC_accum(l_f,h_f,age_j,xi_1,xi_2,del_h,h_l)
% l_f: d variable that affects h’
% h_f: current-period value of h
h_f_prime = h_f + (xi_1 + xi_2*age_j)*l_f – del_h*(1-l_f);
h_f_prime = max(h_f_prime, h_l);
%h_f_prime = h_f;
%h_f_prime = max(h_f_prime, h_l);
end %end function f_HC_accum
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [jl,omega] = find_loc(x_grid,xi)
%————————————————————————-%
% DESCRIPTION
% Find jl s.t. x_grid(jl)<=xi<x_grid(jl+1)
% for jl=1,..,N-1
% omega is the weight on x_grid(jl) so that
% omega*x_grid(jl)+(1-omega)*x_grid(jl+1)=xi
% INPUTS
% x_grid must be a strictly increasing column vector (nx,1)
% xi must be a scalar
% OUTPUTS
% jl: Left point (scalar)
% omega: weight on the left point (scalar)
% NOTES
% See find_loc_vec.m for a vectorized version.
%————————————————————————-%
nx = size(x_grid,1);
jl = max(min(locate(x_grid,xi),nx-1),1);
%Weight on x_grid(j)
omega = (x_grid(jl+1)-xi)/(x_grid(jl+1)-x_grid(jl));
omega = max(min(omega,1),0);
end %end function "find_loc"
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function jl = locate(xx,x)
%function jl = locate(xx,x)
%
% x is between xx(jl) and xx(jl+1)
%
% jl = 0 and jl = n means x is out of range
%
% xx is assumed to be monotone increasing
n = length(xx);
if x<xx(1)
jl = 0;
elseif x>xx(n)
jl = n;
else
jl = 1;
ju = n;
while (ju-jl>1)
jm = floor((ju+jl)/2);
if x>=xx(jm)
jl = jm;
else
ju=jm;
end
end
end
end %end function locateI have to solve a dynamic programming problem with a finite horizon and I am trying to vectorize as much as possible for speed. I attach here a MWE so that everybody can run it
I have run the code with the Matlab profiler and indetified two bottlenecks, marked with the comment line % THIS IS SLOW ACCORDING TO PROFILER
The first bottleneck is the call to the function ReturnFn that builds the n_a*n_a matrix RetMat. I have done this in a vectorized way (inside the function RetMat, I do not use loops)
The second bottleneck is the maximization of RetMat along the first dimension.
I’d be very grateful for any comment/suggestion!

P.S. I have an if condition inside the loops to check for bugs, I removed it now that I am confident about the code but the speed improvement is marginal.
clear,clc,close all
% STATE VARIABLES
% V(a,h,z,j)
% a: asset holdings
% h: human capital (female)
% z: labor productivity shocks (eta_m, eta_f are shocks, theta is permanent)
% j: age (from 1 to N_j)
% CHOICE VARIABLES
% d: Female labor supply (only extensive margin: either 0 or 1)
% a’: Next-period assets
% h’ is implied by (d,a) and consumption is implied by the budget
% constraint
% DYNAMIC PROGRAMMING PROBLEM
% V(a,h,z,j) = max_{d,a’} F(d,a’,a,h,z,j)+beta*s_j*E[V(a’,h’,z’,j+1)|z]
% subject to
% h’=G(d,h), law of motion for human capital
verbose = 1;
%% Define grids and grid sizes
N_j = 80;
n_a = 51;
n_h = 11;
n_z = 50;
n_d = 2;
a_grid = linspace(0,450,n_a)’;
h_grid = linspace(0,0.72,n_h)’;
z_grid = linspace(0.9,1.1,n_z)’;
d_grid = [0,1]’;
z_grid = repmat(z_grid,[1,3]);
pi_z = rand(n_z,n_z);
pi_z = pi_z./sum(pi_z,2);
aprime_val = a_grid; %(a’,1)
a_val = a_grid’; %(1,a)
%% Set parameters that do not depend on age
beta = 0.98;
r = 0.04;
w_m = 1;
w_f = 0.75;
crra = 2;
nu = 0.12;
Jr = 45;
xi_1 = 0.05312;
xi_2 = -0.00188;
del_h = 0.074;
h_l = 0;
p.eff_j = ones(N_j,1);
p.pchild_j = ones(N_j,1);
p.pen_j = ones(N_j,1);
p.nchild_j = ones(N_j,1);
p.s_j = ones(N_j,1);
p.age_j = (1:1:N_j)’;
% Initialize output arrays
V = zeros(n_a,n_h,n_z,N_j);
Policy = zeros(2,n_a,n_h,n_z,N_j);
tic
%% Solve problem in the last period
% Set age-dependent parameters
eff_j = p.eff_j(N_j);
pchild_j = p.pchild_j(N_j);
pen_j = p.pen_j(N_j);
nchild_j = p.nchild_j(N_j);
age_j = p.age_j(N_j);
s_j = p.s_j(N_j);
% V(a,h,z,N_j) = max_{d,a’}
V_d = zeros(n_a,n_h,n_z,n_d);
Pol_aprime_d = zeros(n_a,n_h,n_z,n_d);
for d_c = 1:n_d
d_val = d_grid(d_c);
for z_c = 1:n_z
eta_m_val = z_grid(z_c,1);
eta_f_val = z_grid(z_c,2);
theta_val = z_grid(z_c,3);
for h_c = 1:n_h
h_val = h_grid(h_c);
% RetMat is (a’,a)
RetMat = ReturnFn(d_val,aprime_val,a_val,h_val,eta_m_val,eta_f_val,theta_val,…
w_m,eff_j,w_f,pchild_j,pen_j,r,nchild_j,crra,nu,age_j,Jr);
[max_val,max_ind] = max(RetMat,[],1);
V_d(:,h_c,z_c,d_c) = max_val;
Pol_aprime_d(:,h_c,z_c,d_c) = max_ind;
end %end h
end %end z
end %end d
[V(:,:,:,N_j),d_max] = max(V_d,[],4);
Policy(1,:,:,:,N_j) = d_max; % Optimal d
for z_c=1:n_z
for h_c=1:n_h
for a_c = 1:n_a
d_star = d_max(a_c,h_c,z_c);
Policy(2,a_c,h_c,z_c,N_j) = Pol_aprime_d(a_c,h_c,z_c,d_star); % Optimal a’
end
end
end
%% Backward iteration over age
for j = N_j-1:-1:1
if verbose==1; fprintf(‘Age %d out of %d n’,j,N_j); end
V_next = V(:,:,:,j+1); %V(a’,h’,z’)
% Set age-dependent parameters
eff_j = p.eff_j(j);
pchild_j = p.pchild_j(j);
pen_j = p.pen_j(j);
nchild_j = p.nchild_j(j);
age_j = p.age_j(j);
s_j = p.s_j(j);
for z_c = 1:n_z
eta_m_val = z_grid(z_c,1);
eta_f_val = z_grid(z_c,2);
theta_val = z_grid(z_c,3);
% Compute EV(a’,h’), given z
% EV = zeros(n_a,n_h);
z_prob = pi_z(z_c,:)’;
% for zprime_c = 1:n_z
% EV = EV+V_next(:,:,zprime_c)*z_prob(zprime_c);
% end %end z’
EV = V_next.*shiftdim(z_prob,-2); %V(a’,h’,z’)*Prob(1,1,z’)
EV = sum(EV,3); %EV(a’,h’)
for d_c=1:n_d
d_val = d_grid(d_c);
for h_c = 1:n_h
h_val = h_grid(h_c);
hprime_val = f_HC_accum(d_val,h_val,age_j,xi_1,xi_2,del_h,h_l);
% Ret_mat is (a’,a)
% THIS IS SLOW ACCORDING TO PROFILER
Ret_mat = ReturnFn(d_val,aprime_val,a_val,h_val,eta_m_val,eta_f_val,theta_val,…
w_m,eff_j,w_f,pchild_j,pen_j,r,nchild_j,crra,nu,age_j,Jr);
%[ind_l,weight_l] = interp_toolkit(hprime_val,h_grid);
[ind_l,weight_l] = find_loc(h_grid,hprime_val);
if ind_l>length(h_grid)-1
error(‘ind_l out of bounds’)
end
EV_interp = EV(:,ind_l)*weight_l+EV(:,ind_l+1)*(1-weight_l);
RHS_mat = Ret_mat+beta*s_j*EV_interp;
% THIS IS SLOW ACCORDING TO PROFILER
[max_val,max_ind] = max(RHS_mat,[],1);
% max_val and max_ind are (1,a)
V_d(:,h_c,z_c,d_c) = max_val; %best V given d
Pol_aprime_d(:,h_c,z_c,d_c) = max_ind; %best a’ given d
end %end h
end %end d
end %end z
[V(:,:,:,j),d_max] = max(V_d,[],4);
Policy(1,:,:,:,j) = d_max; % Optimal d
for z_c=1:n_z
for h_c=1:n_h
for a_c = 1:n_a
d_star = d_max(a_c,h_c,z_c);
Policy(2,a_c,h_c,z_c,j) = Pol_aprime_d(a_c,h_c,z_c,d_star); % Optimal a’
end
end
end
end %end j
toc
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function F = ReturnFn(l_f,aprime,a,h_f,eta_m,eta_f,theta,w_m,eff_j,w_f,…
pchild_j,pen_j,r,nchild_j,crra,nu,agej,Jr)
% Calculate earnings (incl. child care costs) of men and women
y_m = w_m*eff_j*theta*eta_m;
y_f = w_f*l_f*(exp(h_f)*theta*eta_f – pchild_j);
% l_f can be either 0 or 1
% calculate available resources
cash = (1+r)*a + pen_j*(agej>=Jr) + (y_m + y_f)*(agej<Jr);
cons = cash-aprime;
%pos = cons>0;
F = (cons/(sqrt(2+nchild_j))).^(1-crra)/(1-crra) – nu*l_f;
F(cons<=0) = -inf;
end %end function "f_ReturnFn"
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [h_f_prime] = f_HC_accum(l_f,h_f,age_j,xi_1,xi_2,del_h,h_l)
% l_f: d variable that affects h’
% h_f: current-period value of h
h_f_prime = h_f + (xi_1 + xi_2*age_j)*l_f – del_h*(1-l_f);
h_f_prime = max(h_f_prime, h_l);
%h_f_prime = h_f;
%h_f_prime = max(h_f_prime, h_l);
end %end function f_HC_accum
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [jl,omega] = find_loc(x_grid,xi)
%————————————————————————-%
% DESCRIPTION
% Find jl s.t. x_grid(jl)<=xi<x_grid(jl+1)
% for jl=1,..,N-1
% omega is the weight on x_grid(jl) so that
% omega*x_grid(jl)+(1-omega)*x_grid(jl+1)=xi
% INPUTS
% x_grid must be a strictly increasing column vector (nx,1)
% xi must be a scalar
% OUTPUTS
% jl: Left point (scalar)
% omega: weight on the left point (scalar)
% NOTES
% See find_loc_vec.m for a vectorized version.
%————————————————————————-%
nx = size(x_grid,1);
jl = max(min(locate(x_grid,xi),nx-1),1);
%Weight on x_grid(j)
omega = (x_grid(jl+1)-xi)/(x_grid(jl+1)-x_grid(jl));
omega = max(min(omega,1),0);
end %end function "find_loc"
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function jl = locate(xx,x)
%function jl = locate(xx,x)
%
% x is between xx(jl) and xx(jl+1)
%
% jl = 0 and jl = n means x is out of range
%
% xx is assumed to be monotone increasing
n = length(xx);
if x<xx(1)
jl = 0;
elseif x>xx(n)
jl = n;
else
jl = 1;
ju = n;
while (ju-jl>1)
jm = floor((ju+jl)/2);
if x>=xx(jm)
jl = jm;
else
ju=jm;
end
end
end
end %end function locate I have to solve a dynamic programming problem with a finite horizon and I am trying to vectorize as much as possible for speed. I attach here a MWE so that everybody can run it
I have run the code with the Matlab profiler and indetified two bottlenecks, marked with the comment line % THIS IS SLOW ACCORDING TO PROFILER
The first bottleneck is the call to the function ReturnFn that builds the n_a*n_a matrix RetMat. I have done this in a vectorized way (inside the function RetMat, I do not use loops)
The second bottleneck is the maximization of RetMat along the first dimension.
I’d be very grateful for any comment/suggestion!

P.S. I have an if condition inside the loops to check for bugs, I removed it now that I am confident about the code but the speed improvement is marginal.
clear,clc,close all
% STATE VARIABLES
% V(a,h,z,j)
% a: asset holdings
% h: human capital (female)
% z: labor productivity shocks (eta_m, eta_f are shocks, theta is permanent)
% j: age (from 1 to N_j)
% CHOICE VARIABLES
% d: Female labor supply (only extensive margin: either 0 or 1)
% a’: Next-period assets
% h’ is implied by (d,a) and consumption is implied by the budget
% constraint
% DYNAMIC PROGRAMMING PROBLEM
% V(a,h,z,j) = max_{d,a’} F(d,a’,a,h,z,j)+beta*s_j*E[V(a’,h’,z’,j+1)|z]
% subject to
% h’=G(d,h), law of motion for human capital
verbose = 1;
%% Define grids and grid sizes
N_j = 80;
n_a = 51;
n_h = 11;
n_z = 50;
n_d = 2;
a_grid = linspace(0,450,n_a)’;
h_grid = linspace(0,0.72,n_h)’;
z_grid = linspace(0.9,1.1,n_z)’;
d_grid = [0,1]’;
z_grid = repmat(z_grid,[1,3]);
pi_z = rand(n_z,n_z);
pi_z = pi_z./sum(pi_z,2);
aprime_val = a_grid; %(a’,1)
a_val = a_grid’; %(1,a)
%% Set parameters that do not depend on age
beta = 0.98;
r = 0.04;
w_m = 1;
w_f = 0.75;
crra = 2;
nu = 0.12;
Jr = 45;
xi_1 = 0.05312;
xi_2 = -0.00188;
del_h = 0.074;
h_l = 0;
p.eff_j = ones(N_j,1);
p.pchild_j = ones(N_j,1);
p.pen_j = ones(N_j,1);
p.nchild_j = ones(N_j,1);
p.s_j = ones(N_j,1);
p.age_j = (1:1:N_j)’;
% Initialize output arrays
V = zeros(n_a,n_h,n_z,N_j);
Policy = zeros(2,n_a,n_h,n_z,N_j);
tic
%% Solve problem in the last period
% Set age-dependent parameters
eff_j = p.eff_j(N_j);
pchild_j = p.pchild_j(N_j);
pen_j = p.pen_j(N_j);
nchild_j = p.nchild_j(N_j);
age_j = p.age_j(N_j);
s_j = p.s_j(N_j);
% V(a,h,z,N_j) = max_{d,a’}
V_d = zeros(n_a,n_h,n_z,n_d);
Pol_aprime_d = zeros(n_a,n_h,n_z,n_d);
for d_c = 1:n_d
d_val = d_grid(d_c);
for z_c = 1:n_z
eta_m_val = z_grid(z_c,1);
eta_f_val = z_grid(z_c,2);
theta_val = z_grid(z_c,3);
for h_c = 1:n_h
h_val = h_grid(h_c);
% RetMat is (a’,a)
RetMat = ReturnFn(d_val,aprime_val,a_val,h_val,eta_m_val,eta_f_val,theta_val,…
w_m,eff_j,w_f,pchild_j,pen_j,r,nchild_j,crra,nu,age_j,Jr);
[max_val,max_ind] = max(RetMat,[],1);
V_d(:,h_c,z_c,d_c) = max_val;
Pol_aprime_d(:,h_c,z_c,d_c) = max_ind;
end %end h
end %end z
end %end d
[V(:,:,:,N_j),d_max] = max(V_d,[],4);
Policy(1,:,:,:,N_j) = d_max; % Optimal d
for z_c=1:n_z
for h_c=1:n_h
for a_c = 1:n_a
d_star = d_max(a_c,h_c,z_c);
Policy(2,a_c,h_c,z_c,N_j) = Pol_aprime_d(a_c,h_c,z_c,d_star); % Optimal a’
end
end
end
%% Backward iteration over age
for j = N_j-1:-1:1
if verbose==1; fprintf(‘Age %d out of %d n’,j,N_j); end
V_next = V(:,:,:,j+1); %V(a’,h’,z’)
% Set age-dependent parameters
eff_j = p.eff_j(j);
pchild_j = p.pchild_j(j);
pen_j = p.pen_j(j);
nchild_j = p.nchild_j(j);
age_j = p.age_j(j);
s_j = p.s_j(j);
for z_c = 1:n_z
eta_m_val = z_grid(z_c,1);
eta_f_val = z_grid(z_c,2);
theta_val = z_grid(z_c,3);
% Compute EV(a’,h’), given z
% EV = zeros(n_a,n_h);
z_prob = pi_z(z_c,:)’;
% for zprime_c = 1:n_z
% EV = EV+V_next(:,:,zprime_c)*z_prob(zprime_c);
% end %end z’
EV = V_next.*shiftdim(z_prob,-2); %V(a’,h’,z’)*Prob(1,1,z’)
EV = sum(EV,3); %EV(a’,h’)
for d_c=1:n_d
d_val = d_grid(d_c);
for h_c = 1:n_h
h_val = h_grid(h_c);
hprime_val = f_HC_accum(d_val,h_val,age_j,xi_1,xi_2,del_h,h_l);
% Ret_mat is (a’,a)
% THIS IS SLOW ACCORDING TO PROFILER
Ret_mat = ReturnFn(d_val,aprime_val,a_val,h_val,eta_m_val,eta_f_val,theta_val,…
w_m,eff_j,w_f,pchild_j,pen_j,r,nchild_j,crra,nu,age_j,Jr);
%[ind_l,weight_l] = interp_toolkit(hprime_val,h_grid);
[ind_l,weight_l] = find_loc(h_grid,hprime_val);
if ind_l>length(h_grid)-1
error(‘ind_l out of bounds’)
end
EV_interp = EV(:,ind_l)*weight_l+EV(:,ind_l+1)*(1-weight_l);
RHS_mat = Ret_mat+beta*s_j*EV_interp;
% THIS IS SLOW ACCORDING TO PROFILER
[max_val,max_ind] = max(RHS_mat,[],1);
% max_val and max_ind are (1,a)
V_d(:,h_c,z_c,d_c) = max_val; %best V given d
Pol_aprime_d(:,h_c,z_c,d_c) = max_ind; %best a’ given d
end %end h
end %end d
end %end z
[V(:,:,:,j),d_max] = max(V_d,[],4);
Policy(1,:,:,:,j) = d_max; % Optimal d
for z_c=1:n_z
for h_c=1:n_h
for a_c = 1:n_a
d_star = d_max(a_c,h_c,z_c);
Policy(2,a_c,h_c,z_c,j) = Pol_aprime_d(a_c,h_c,z_c,d_star); % Optimal a’
end
end
end
end %end j
toc
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function F = ReturnFn(l_f,aprime,a,h_f,eta_m,eta_f,theta,w_m,eff_j,w_f,…
pchild_j,pen_j,r,nchild_j,crra,nu,agej,Jr)
% Calculate earnings (incl. child care costs) of men and women
y_m = w_m*eff_j*theta*eta_m;
y_f = w_f*l_f*(exp(h_f)*theta*eta_f – pchild_j);
% l_f can be either 0 or 1
% calculate available resources
cash = (1+r)*a + pen_j*(agej>=Jr) + (y_m + y_f)*(agej<Jr);
cons = cash-aprime;
%pos = cons>0;
F = (cons/(sqrt(2+nchild_j))).^(1-crra)/(1-crra) – nu*l_f;
F(cons<=0) = -inf;
end %end function "f_ReturnFn"
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [h_f_prime] = f_HC_accum(l_f,h_f,age_j,xi_1,xi_2,del_h,h_l)
% l_f: d variable that affects h’
% h_f: current-period value of h
h_f_prime = h_f + (xi_1 + xi_2*age_j)*l_f – del_h*(1-l_f);
h_f_prime = max(h_f_prime, h_l);
%h_f_prime = h_f;
%h_f_prime = max(h_f_prime, h_l);
end %end function f_HC_accum
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [jl,omega] = find_loc(x_grid,xi)
%————————————————————————-%
% DESCRIPTION
% Find jl s.t. x_grid(jl)<=xi<x_grid(jl+1)
% for jl=1,..,N-1
% omega is the weight on x_grid(jl) so that
% omega*x_grid(jl)+(1-omega)*x_grid(jl+1)=xi
% INPUTS
% x_grid must be a strictly increasing column vector (nx,1)
% xi must be a scalar
% OUTPUTS
% jl: Left point (scalar)
% omega: weight on the left point (scalar)
% NOTES
% See find_loc_vec.m for a vectorized version.
%————————————————————————-%
nx = size(x_grid,1);
jl = max(min(locate(x_grid,xi),nx-1),1);
%Weight on x_grid(j)
omega = (x_grid(jl+1)-xi)/(x_grid(jl+1)-x_grid(jl));
omega = max(min(omega,1),0);
end %end function "find_loc"
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function jl = locate(xx,x)
%function jl = locate(xx,x)
%
% x is between xx(jl) and xx(jl+1)
%
% jl = 0 and jl = n means x is out of range
%
% xx is assumed to be monotone increasing
n = length(xx);
if x<xx(1)
jl = 0;
elseif x>xx(n)
jl = n;
else
jl = 1;
ju = n;
while (ju-jl>1)
jm = floor((ju+jl)/2);
if x>=xx(jm)
jl = jm;
else
ju=jm;
end
end
end
end %end function locate vectorization, performance MATLAB Answers — New Questions

​