I am currently getting the error: Index exceeds number of array elements.

The exact error I get is,

I am currently getting the error: Index exceeds number of array elements. Index must not exceed 747. Error in indexing (line 1075) R_tilde = builtin('subsref',L_tilde,Idx);

I haven't included the full code but what is leading to the error is (note that the w2, w3, x2, x3 are calculated in earlier code but this is what they are given to be when I use them here, they have purposefully been rounded them to make them integers and to be used here):

% Multiple Switching - Equilibrium strategy starting with  = 0 on [0,1).

% Input values,  and 
beta = 1000
alpha = 600
kappa = 1 - alpha/beta % Influences the gamma* graphs, not sure how yet look at.
sigma = 4 % Appears to change the gamma graphs gradient/where it is at t = 1. It
          % also changes the value of theta1 (first intersection). The
          % larger the value of sigma and thus the larger the value of a,
          % the closer in value theta1 is to -log(kepa)/lambda.
a = sigma/2
lambda = 3 % Appears to change the gamma graphs gradient/where it is at t = 1. It
           % also changes the value of theta1 (first intersection).
b = -log(kappa)/lambda
% Calculating the  values up to 1:
% Calculating :
syms theta
g = 2*(normcdf(a) + (theta - 1) .* normcdf(a*(1-theta)) + 1/(a*sqrt(2*pi)) * (exp(-(a.^2)/2) - ...
    exp(-((a.*(1-theta)).^2)/2))) - theta == b;
soltheta1 = vpasolve(g, theta)
w = round(soltheta1, 3)
w1 = w + 0.001
w2 = 1000 * w1
w3 = (w + 0.002) * 1000
% Calculating :
syms theta2
g = 2*(normcdf(a) + (theta2 - 1) .* normcdf(a*(1-theta2)) + 1/(a*sqrt(2*pi)) * (exp(-(a.^2)/2) - ...
    exp(-((a.*(1-theta2)).^2)/2))) - theta2 == b + 2*(normcdf(a) + ((theta2 + (1 - kappa)/(kappa * lambda * normcdf(sigma*(1-theta2)/2))) - 1) .* normcdf(a*(1-(theta2 + (1 - kappa)/(kappa * lambda * normcdf(sigma*(1-theta2)/2))))) + 1/(a*sqrt(2*pi)) * (exp(-(a.^2)/2) - ...
    exp(-((a.*(1-(theta2 + (1 - kappa)/(kappa * lambda * normcdf(sigma*(1-theta2)/2))))).^2)/2))) - (theta2 + (1 - kappa)/(kappa * lambda * normcdf(sigma*(1-theta2)/2)));
soltheta2 = vpasolve(g, theta2)
x = round(soltheta2, 3)
x1 = x + 0.001
x2 = 1000 * x1
x3 = (x + 0.002) * 1000
% Calculating :
syms theta3
g = theta3 == soltheta2 + (1 - kappa)/(kappa * lambda * normcdf(sigma*(1-soltheta2)/2));
soltheta3 = vpasolve(g, theta3)
% Calculate  and the different types of  used.
t = 0:0.001:1
t1 = 0:0.001:w
t2 = w1:0.001:x
t3 = x1:0.001:1
sigmahat1 = sigma*(1-t1)
sigmahat2 = sigma*(1-t2)
sigmahat3 = sigma*(1-t3)
phi1 = normcdf(sigmahat1/2)
h1 = (2*phi1 - 1)
phi2 = normcdf(sigmahat2/2)
h2 = (2*phi2 - 1)
phi3 = normcdf(sigmahat3/2)
h3 = (2*phi3 - 1)
phi4 = normcdf(a)
phi5 = normcdf(a*(1-t1))
phi6 = normcdf(a*(1-t2))
phi7 = normcdf(a*(1-t3))
phi8 = normcdf(a*(1-t))
phi3theta = normcdf(a*(1-soltheta1))
% Is the small the derivative of big  and if so is this therefore the derivative of the cdf which would be the pdf?
% Calculate g and .
gt = 2*(phi4 + (t - 1) .* phi8 + 1/(a*sqrt(2.*pi)) * (exp(-(a.^2)/2) - exp(-(( ...
    a.*(1-t)).^2)/2))) - t
g34 = 2*(phi4 + (t1 - 1) .* phi5 + 1/(a*sqrt(2*pi)) * (exp(-(a.^2)/2) - exp(-(( ...
    a.*(1-t1)).^2)/2))) - t1
g35 = 2*(phi4 + (t2 - 1) .* phi6 + 1/(a*sqrt(2*pi)) * (exp(-(a.^2)/2) - exp(-(( ...
    a.*(1-t2)).^2)/2))) - t2
g36 = 2*(phi4 + (t3 - 1) .* phi7 + 1/(a*sqrt(2*pi)) * (exp(-(a.^2)/2) - exp(-(( ...
    a.*(1-t3)).^2)/2))) - t3
gtlarge3 = 2*(phi4 + (soltheta1 - 1) * phi3theta + 1/(a*sqrt(2*pi)) * (exp(-(a^2)/2) - exp(-(( ...
    a*(1-soltheta1))^2)/2))) - soltheta1
gamma1 = exp(-lambda .* gt)
gamma34 = exp(-lambda .* g34)
gamma35 = exp(-lambda .* g35)
gamma36 = exp(-lambda .* g36)
% Calculating the switching curves gammastar t <=  when = 0, for and gammastar  t >  when = 0.
gammastar1 = zeros(1,1001)
gammastar2 = zeros(1,1001)
for index = 1:1:1001
    if index < w3
        gammastar1(index) = ((gamma34(index))./(kappa * (1 + lambda * (soltheta1 - t1(index)) .* h1(index))))
        gammastar2(index) = 1
    elseif (index < x3) && (index > w2)  
        gammastar1(index) = exp(-lambda.*(g35(index - w2) ...
            -gtlarge3)) % one of these is wrong
        % gammastar2(index) = ((gamma35(index - w2))./(kappa* (1 + lambda * (soltheta2 - t2(index - w2)) .* h2(index - w2))))
    else 
        gammastar1(index) = ((gamma36(index - x2))./(kappa * (1 + lambda * (soltheta2 - t3(index - x2)) .* h3(index - x2))))
        % gammastar2(index) = ((gamma36(index - x2))./(kappa * (1 + lambda * (soltheta2 - t3(index - x2)) .* h3(index - x2))))

    end
end


% gammastar1 = ((gamma1)./(kappa * (1 + lambda * (soltheta1 - t1) .* h1)))
% gammastar2 = ((gamma1)./(kappa * (1 + lambda * (soltheta2 - t) .* h)))
% gammat = exp(-lambda.*(gt-gtlarge3))
plot(t, gammastar1, t, gamma1, t, gammastar2)

I am not sure why this doesn't work because the elseif to me is meaning 229 < index < 977 which would be 747 elements (which is the amount of elements in g35), so why would it be saying index must not exceed 747? It seems to me that it is trying index = 977 which is causing an error because there is no 748th element in g35 but uf this is the case why would it be putting index = 977 through when I have said index < 977.

PS: Apologise if the code takes a while to run, it may not be written very eifficiently.

[Image Analyst]

w2 = 229.0

w2 = 229

w3 = 230.0

w3 = 230

x2 = 976.0

x2 = 976

x3 = 977.0

x3 = 977

gammastar1 = zeros(1,1001)

gammastar1 = 1×1001

0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0

gammastar2 = zeros(1,1001)

gammastar2 = 1×1001

0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0

for index = 1:1:1001
    if index < w3
        gammastar1(index) = ((gamma34(index))./(kappa * (1 + lambda * (soltheta1 - t1(index)) .* h1(index))))
        gammastar2(index) = 1
    elseif (index < x3) && (index > w2)  
        gammastar1(index) = exp(-lambda.*(g35(index - w2) ...
            -gtlarge3(index - w2))) % one of these is wrong
        gammastar2(index) = ((gamma35(index - w2))./(kappa* (1 + lambda * (soltheta2 - t2(index - w2)) .* h2(index - w2))))
    else 
        gammastar1(index) = ((gamma36(index - x2))./(kappa * (1 + lambda * (soltheta2 - t3(index - x2)) .* h3(index - x2))))
        gammastar2(index) = ((gamma36(index - x2))./(kappa * (1 + lambda * (soltheta2 - t3(index - x2)) .* h3(index - x2))))
        
    end
end

Unrecognized function or variable 'gamma34'.

plot(t, gammastar1, t, gamma1, t, gammastar2)

Please give all missing variables so we can run your code.

[Charles]

Apologies, I have updated the question.

[Cris LaPierre]

Also, please share the entire error message (all the red text).

This error appears to be related to symbolic variables, so it would be helpful to have a working example that duplicates the error you are getting.

Answer 1

[David Hill]

Change symbolic to double for comparison. Below works.

for index = 1:1:1001
    if index < double(w3)
        gammastar1(index) = ((gamma34(index))./(kappa * (1 + lambda * (soltheta1 - t1(index)) .* h1(index))))
        gammastar2(index) = 1
    elseif (index < double(x3)) && (index > double(w2))  
        gammastar1(index) = exp(-lambda.*(g35(index - w2) ...
            -gtlarge3)) % one of these is wrong
        % gammastar2(index) = ((gamma35(index - w2))./(kappa* (1 + lambda * (soltheta2 - t2(index - w2)) .* h2(index - w2))))
    else 
        gammastar1(index) = ((gamma36(index - x2))./(kappa * (1 + lambda * (soltheta2 - t3(index - x2)) .* h3(index - x2))))
        % gammastar2(index) = ((gamma36(index - x2))./(kappa * (1 + lambda * (soltheta2 - t3(index - x2)) .* h3(index - x2))))
    end
end

Answer 2

[Walter Roberson]

Change some of the variables to double, and make sure to round() some of the computations.

% Multiple Switching - Equilibrium strategy starting with  = 0 on [0,1).

% Input values,  and 
beta = 1000;
alpha = 600;
kappa = 1 - alpha/beta; % Influences the gamma* graphs, not sure how yet look at.
sigma = 4; % Appears to change the gamma graphs gradient/where it is at t = 1. It
          % also changes the value of theta1 (first intersection). The
          % larger the value of sigma and thus the larger the value of a,
          % the closer in value theta1 is to -log(kepa)/lambda.
a = sigma/2;
lambda = 3; % Appears to change the gamma graphs gradient/where it is at t = 1. It
           % also changes the value of theta1 (first intersection).
b = -log(kappa)/lambda;
% Calculating the  values up to 1:
% Calculating :
syms theta
g = 2*(normcdf(a) + (theta - 1) .* normcdf(a*(1-theta)) + 1/(a*sqrt(2*pi)) * (exp(-(a.^2)/2) - ...
    exp(-((a.*(1-theta)).^2)/2))) - theta == b;
soltheta1 = vpasolve(g, theta);
disp(soltheta1)

0.34069559421320661731328878287072

w = round(double(soltheta1), 3);
w1 = w + 0.001;
w2 = round(1000 * w1);
w3 = round((w + 0.002) * 1000);
display(w); display(w1); display(w2); display(w3);

w = 0.3410

w1 = 0.3420

w2 = 342

w3 = 343

% Calculating :
syms theta2
g = 2*(normcdf(a) + (theta2 - 1) .* normcdf(a*(1-theta2)) + 1/(a*sqrt(2*pi)) * (exp(-(a.^2)/2) - ...
    exp(-((a.*(1-theta2)).^2)/2))) - theta2 == b + 2*(normcdf(a) + ((theta2 + (1 - kappa)/(kappa * lambda * normcdf(sigma*(1-theta2)/2))) - 1) .* normcdf(a*(1-(theta2 + (1 - kappa)/(kappa * lambda * normcdf(sigma*(1-theta2)/2))))) + 1/(a*sqrt(2*pi)) * (exp(-(a.^2)/2) - ...
    exp(-((a.*(1-(theta2 + (1 - kappa)/(kappa * lambda * normcdf(sigma*(1-theta2)/2))))).^2)/2))) - (theta2 + (1 - kappa)/(kappa * lambda * normcdf(sigma*(1-theta2)/2)));
soltheta2 = vpasolve(g, theta2);
x = round(double(soltheta2), 3);
x1 = x + 0.001;
x2 = 1000 * x1;
x3 = (x + 0.002) * 1000;
display(x); display(x1); display(x2); display(x3);

x = 0.8600

x1 = 0.8610

x2 = 861

x3 = 862

% Calculating :
syms theta3
g = theta3 == soltheta2 + (1 - kappa)/(kappa * lambda * normcdf(sigma*(1-soltheta2)/2));
soltheta3 = vpasolve(g, theta3);
% Calculate  and the different types of  used.
t = 0:0.001:1;
t1 = 0:0.001:w;
t2 = w1:0.001:x;
t3 = x1:0.001:1;
sigmahat1 = sigma*(1-t1);
sigmahat2 = sigma*(1-t2);
sigmahat3 = sigma*(1-t3);
phi1 = normcdf(sigmahat1/2);
h1 = (2*phi1 - 1);
phi2 = normcdf(sigmahat2/2);
h2 = (2*phi2 - 1);
phi3 = normcdf(sigmahat3/2);
h3 = (2*phi3 - 1);
phi4 = normcdf(a);
phi5 = normcdf(a*(1-t1));
phi6 = normcdf(a*(1-t2));
phi7 = normcdf(a*(1-t3));
phi8 = normcdf(a*(1-t));
phi3theta = normcdf(a*(1-soltheta1));
% Is the small the derivative of big  and if so is this therefore the derivative of the cdf which would be the pdf?
% Calculate g and .
gt = 2*(phi4 + (t - 1) .* phi8 + 1/(a*sqrt(2.*pi)) * (exp(-(a.^2)/2) - exp(-(( ...
    a.*(1-t)).^2)/2))) - t;
g34 = 2*(phi4 + (t1 - 1) .* phi5 + 1/(a*sqrt(2*pi)) * (exp(-(a.^2)/2) - exp(-(( ...
    a.*(1-t1)).^2)/2))) - t1;
g35 = 2*(phi4 + (t2 - 1) .* phi6 + 1/(a*sqrt(2*pi)) * (exp(-(a.^2)/2) - exp(-(( ...
    a.*(1-t2)).^2)/2))) - t2;
g36 = 2*(phi4 + (t3 - 1) .* phi7 + 1/(a*sqrt(2*pi)) * (exp(-(a.^2)/2) - exp(-(( ...
    a.*(1-t3)).^2)/2))) - t3;
gtlarge3 = 2*(phi4 + (soltheta1 - 1) * phi3theta + 1/(a*sqrt(2*pi)) * (exp(-(a^2)/2) - exp(-(( ...
    a*(1-soltheta1))^2)/2))) - soltheta1;
gamma1 = exp(-lambda .* gt);
gamma34 = exp(-lambda .* g34);
gamma35 = exp(-lambda .* g35);
gamma36 = exp(-lambda .* g36);
% Calculating the switching curves gammastar t <=  when = 0, for and gammastar  t >  when = 0.
gammastar1 = zeros(1,1001);
gammastar2 = zeros(1,1001);
for index = 1:1:1001
    if index < w3
        gammastar1(index) = ((gamma34(index))./(kappa * (1 + lambda * (soltheta1 - t1(index)) .* h1(index))));
        gammastar2(index) = 1;
    elseif (index < x3) && (index > w2)  
        gammastar1(index) = exp(-lambda.*(g35(index - w2) ...
            -gtlarge3)); % one of these is wrong
        % gammastar2(index) = ((gamma35(index - w2))./(kappa* (1 + lambda * (soltheta2 - t2(index - w2)) .* h2(index - w2))))
    else 
        gammastar1(index) = ((gamma36(index - x2))./(kappa * (1 + lambda * (soltheta2 - t3(index - x2)) .* h3(index - x2))));
        % gammastar2(index) = ((gamma36(index - x2))./(kappa * (1 + lambda * (soltheta2 - t3(index - x2)) .* h3(index - x2))))

    end
end


% gammastar1 = ((gamma1)./(kappa * (1 + lambda * (soltheta1 - t1) .* h1)))
% gammastar2 = ((gamma1)./(kappa * (1 + lambda * (soltheta2 - t) .* h)))
% gammat = exp(-lambda.*(gt-gtlarge3))
plot(t, gammastar1, t, gamma1, t, gammastar2)