Index in position 1 exceeds array bounds (must not exceed 1). Error in pij (line 14) pmat(k,n) = pis(k,1) .* tau(k,n); Error in pfij (line 18)

Question

0 votos

Please my code is not running when I run the main script GE_fsolve. I can attached my main codes upon request. At the moment I inserted copies of all functions that I am using. Please kindly help me with this. Thanks!!!!

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Create output price function
function [Pis] = pis(wis)
global Ais beta IOcoef Nsec Ncoun
% reformat IOcoef to (NS,S) instead of the initial (NS,NS)
IOmat_cs_s = nan(Ncoun*Nsec,Nsec); %cs-s
for t = 1:1:Nsec
    IOmat_cs_s(:,t) = sum(IOcoef(:,((t-1)*Ncoun+1):Ncoun*t),2); %cs-s
end
IOmat_cs = sum(IOmat_cs_s,2); %change IOcoef to (NS,1) instead of (NS,S)
% create matrix
pimat = zeros(Nsec*Ncoun,1);
% compute pis(NS,1)
for i = 1:1:Nsec*Ncoun
    pimat(i) = ((beta(i)).^(-beta(i)) .* (1-beta(i)).^(-1+beta(i)))...
        .*(wis(i).^beta(i) .* (IOmat_cs(i).*PMis(pij(pis(i)))).^(1-beta(i)))./Ais(i); 
end
Pis = pimat; %(NS,1)
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Create pirce of output function
function [Pij] = pij(pis)
global tau Nsec Ncoun
% create matrix
pmat = zeros(Nsec*Ncoun,Ncoun);
 
% solve for price of output(n,k)
for k = 1:1:Nsec*Ncoun
    for n = 1:1:Ncoun
        pmat(k,n) = pis(k) .* tau(k,n); 
    end
end
Pij = pmat; %(NS,N)
end 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Create CES final demand function
function [Fij] = fij(pis,wis,mij)
global alpha sigma Nsec Ncoun
 
% PFis = PFis(pij,sigma,Ncoun,Nsec)
% reformat sigma to (N*S,1) instead of initial (S,1)
for s = 1:1:Nsec
    idx  = 1+(s-1)*Ncoun:1:Ncoun*s; %for every s what is idx 
    SG(idx,1) = ones(Ncoun,1)*sigma(s); %for every s use the respective idx as row and multiply the actual s value by col vector 
end
 
% create matrices
fmat = zeros(Nsec*Ncoun,Ncoun);
for i = 1:1:Nsec*Ncoun
    for j = 1:1:Ncoun
    fmat(:,:) = alpha(i,j) .* (pij(pis(i,1))).^(-(SG*ones(1,Ncoun)))...
                .*(PFis(pij(pis(i,1)))).^((SG*ones(1,Ncoun))-1) .* Yis(wis,pis,mij);   
    end
end     
Fij = fmat; %(NS,N)
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Create CES intermediate demand function
function [Mij] = mij(pis,mij)
global gama rho Nsec Ncoun
 
% reformat rho to (N*S,1) instead of initial (S,1)
for s = 1:1:Nsec
    idx  = 1+(s-1)*Ncoun:1:Ncoun*s; %for every s what is idx 
    SG(idx,1) = ones(Ncoun,1)*rho(s); %for every s use the respective idx as row and multiply the actual s value by col vector 
end
 
% create matrices
mmat = zeros(Nsec*Ncoun,Nsec*Ncoun);
for i = 1:1:Nsec*Ncoun
    for j = 1:1:Ncoun
    mmat(i,j) = gama(i,j) .* (pij(pis(i,1))).^(-(SG*ones(1,Ncoun)))...
              .* (PMis(pij(pis(i,1)))).^((SG*ones(1,Ncoun))-1) .* IMis(pis,mij); %(NS,N)
    end
end
Mij = mmat; %(NS*N)
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Create final demand total value of imports
function [Pfij] = pfij(pis,wis,mij)
global alpha sigma Nsec Ncoun
 
% reformat sigma to (N*S,1) instead of initial (S,1)
for s = 1:1:Nsec
    idx  = 1+(s-1)*Ncoun:1:Ncoun*s; %for every s what is idx 
    SG(idx,1) = ones(Ncoun,1)*sigma(s); %for every s use the respective idx as row and multiply the actual s value by col vector 
end
 
% create matrices
pfmat = zeros(Nsec*Ncoun,Ncoun);
for i = 1:1:Nsec*Ncoun
    for j = 1:1:Ncoun
    pfmat(:,:) = alpha(i,j) .* (pij(pis(i,1))).^(1-(SG*ones(1,Ncoun)))...
                .*(PFis(pij(pis(i,1)))).^(1-(SG*ones(1,Ncoun))) .* Yis(wis,pis,mij);  
    end
end
Pfij = pfmat; %(NS,N)
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Create intermediate demand total value of imports
function [Pmij] = pmij(pis,mij)
global gama rho Nsec Ncoun
 
% reformat rho to (N*S,1) instead of initial (S,1)
for s = 1:1:Nsec
    idx  = 1+(s-1)*Ncoun:1:Ncoun*s; %for every s what is idx 
    SG(idx,1) = ones(Ncoun,1)*rho(s); %for every s use the respective idx as row and multiply the actual s value by col vector 
end
 
% create matrices
pmmat = zeros(Nsec*Ncoun,Nsec*Ncoun);
for i = 1:1:Nsec*Ncoun
    for j = 1:1:Ncoun
    pmmat(i,j) = gama(i,j) .* (pij(pis(i,1))).^(1-(SG*ones(1,Ncoun)))...
              .* (PMis(pij(pis(i,1)))).^(1-(SG*ones(1,Ncoun))) .* IMis(pis,mij); %(NS,N)
    end
end
Pmij = pmmat; %(NS,N)
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Solve GE Model
function syseqn = GEsol(x)
global Ais Lis beta IOmat_cs Ncoun   
    % Define unknown variables x
    % Break x into 4 separate variables
    % System x=(   mij,     fij,     pis,      wi )
    %           (N*S,N*S)  (N*S,N)  (N*S,1)  (NS,1)    
    mij = x(:,1:Ncoun);
    fij = x(:,Ncoun+1:end-2);
    pis = x(:,end-1);
    wis = x(:,end);
 
    % Safeguard restrictions
    tol  = 1e-10; %Tolerance, use to keep certain variables away from zero
    mij  = max(mij, tol);
    fij  = max(fij, tol);
    pis  = max(pis, tol);
    wis  = max(wis, tol);
   
    % Equations to be solve simultaneously
    eq1 = pij(pis).*fij - pfij(pis,wis,mij);  %final demand import
    eq2 = pij(pis).*mij - pmij(pis,wis,mij);  %intermeditate demand import
    eq3 = pis - (beta.^(-beta).*(1-beta).^(-1+beta))... 
          .* (wis.^beta .* (IOmat_cs.*PMis(pij(pis))).^(1-beta))./Ais;  %output price equation
    eq4 = Lis -  (((sum(fij,2) + sum(mij,2)))./Ais) * ((beta.*PMis(pij(pis)))./...
          ((1-beta)*wis)).^(1-beta);  %labor market clearing condition
 
syseqn = [eq1, eq2, eq3, eq4];
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Main Script PrograM
% ***********************************
close all;
clear;
clc; 
tic;
 
global alpha gama sigma rho tau Ais Lis beta IOcoef IOmat_cs Nsec Ncoun
 
%******************************************
%    Define country-sector inputs         %
%******************************************
 
Nsec = 2;
Ncoun = 4;
 
%******************************************
%         Define parameters               %
%******************************************
 
sigma  = 5*ones(Nsec,1); %elasticity of substitution for households
rho    = 5*ones(Nsec,1); %elasticity of substitution for firm
beta   = 0.6*ones(Nsec*Ncoun,1); %value-added share
Lis     = 1*ones(Nsec*Ncoun,1); %labor (supply) endowment
Ais    = 1*ones(Nsec*Ncoun,1); %TFP
tau    = 1.5*ones(Nsec*Ncoun,Ncoun); %tariff
alpha  = 0.6*ones(Nsec*Ncoun,Ncoun); %sector share of expenditure for final demand
gama   = 0.4*ones(Nsec*Ncoun,Ncoun); %sector share of expenditure for intermediate demand
IOcoef = importdata('IOcoef.txt');
 
% reformat IOcoef to (NS,S) instead of the initial (NS,NS)
IOmat_cs_s = nan(Ncoun*Nsec,Nsec); %cs-s
for t = 1:1:Ncoun
    IOmat_cs_s(:,t) = sum(IOcoef(:,((t-1)*Nsec+1):Nsec*t),2); %cs-s
end
IOmat_cs = sum(IOmat_cs_s,2); %change IOcoef to (NS,1) instead of (NS,S)
 
% % parameters put in a vector to enable easy passing to the function
% P = [beta,Lis,Ais,IOcoef];
%******************************************
%   Solve non-linear equations(fsolve)    %
%******************************************
 
options = optimset('Display', 'iter', 'MaxFunEvals', 20000, 'MaxIter', 20000, 'TolFun', 1e-10, 'TolX', 1e-10);
 
% initial guess vakues for x matrix
x0 = 0.5.*ones(8,10);
 
% call output of equation function under initial guess
% syseqn = GEsol(x0,P);
% syseqn = GEsol(x0,P,Nsec,Ncoun);
syseqn = GEsol(x0);
 
% fsolve solves the function by looking for matrix x0 using initial guess
% x = fsolve('GEsol', x0, options,P);
x = fsolve('GEsol', x0, options);
 
% solution of the model in terms of x
mij = x(:,1:Ncoun);
fij = x(:,Ncoun+1:end-2);
pis = x(:,end-1);
wis = x(:,end);

55 comentarios
Mostrar 53 comentarios más antiguos Ocultar 53 comentarios más antiguos

Mathieu NOE el 6 de Sept. de 2021

Abrir en MATLAB Online

hello

your code is not very easy to debugg , but I could see already a couple of issues :

probably the most annoying one - as far as I understand how you constructed the whole thing :

so we start with the main code (GE_fsolve.m) , and first subfunction called is GEsol.m
in GEsol.m , there are these lines

    % Equations to be solve simultaneously
    eq1 = mypij(wis).*fij - mypfij(wis,mij);  %final demand import
    eq2 = mypij(wis).*mij - mypmij(wis,mij);  %intermeditate demand import
    eq3 = pis - (beta.^(-beta).*(1-beta).^(-1+beta))... 
          .* (wis.^beta .* (IOmat_cs.*myPMis(mypij(wis))).^(1-beta))./Ais;  %output price equation
    eq4 = Lis -  (((sum(fij,2) + sum(mij,2)))./Ais) * ((beta.*PMis(pij(wis)))./...
          ((1-beta)*wis)).^(1-beta);  %labor market clearing condition

where we call the subfunction mypij.m , (not only) and this function calls mypis.m

(BTW why the index k is no more used in the computation ? )

% solve for price of output(k,n)
for k = 1:1:Nsec*Ncoun
    for n = 1:1:Ncoun
        pmat(:,n) = mypis(wis(:,1)) .* tau(:,n); 
    end
end

and mypis.m then use again mypij.m .... this cannot work !

% compute pis(NS,1)
for i = 1:1:Nsec*Ncoun
    pimat(:,1) = ((beta(:,1)).^(-beta(:,1)) .* (1-beta(:,1)).^(-1+beta(:,1)))...
        .*(wis(:,1).^beta(:,1) .* (IOmat_cs(:,1).*myPMis(mypij(wis(:,1))))...
        .^(1-beta(:,1)))./Ais(:,1); 
end

so we have a problem here; a subfunction cannot call herself even through another (intermediate) subfunction

this could be the root cause of your infinite loop - or - matlab is stucked in this impossible scenario;

Walter Roberson el 6 de Sept. de 2021

Abrir en MATLAB Online

@Mathieu NOE's reply was not present when I wrote mine... not sure why it did not show up.

An example of something that does not work:

function C = A(x)
   C = B(x-1);
end
function D = B(x)
   D = A(x-1);
end

In order to execute A, you need to finish executing a call to B. In order to finish executing a call to B, you need to finish executing a call to A -- but in order to finish that call to A, you need to call B to call A to call B to ... it never ends until MATLAB says you are out of memory or you have made too many calls.

An example of something that does work

function C = A(x)
   C = B(x-1);
end
function D = B(x)
   if all(isinf(x) | x < 0)
     D = max(-2, x);
   else
     D = A(x-1);
   end
end

For example, call A([2 4]) then that calls B([1 3]) . In B, no values of [1 3] are infinite or negative, so you take the else and call A([0 2]) . A then calls B([-1 1]) . None of the elements of [-1 1] are infinite, and at least one of the values is not < 0, so the else is taken and A([-2 0]) is called. A then calls B([-3 -1]) . None of the elements of [-3 -1] are infinite, but all of the values of [-3 -1] are < 0, so the if is taken. max(-2, [-3 -1]) is [-2 -1] and that is the output, so that call to B with [-3 -1] finishes, and returns to A, which finishes and returns to the B instance processing [-1 1], which finishes and returns to A, which finishes and returns to the B that is processing [1 3], which returns to the A that is processing [2 4], which finishes -- and the final output is [-2 -1] .

So the problem is not having one function calling another and the other calling the first: the problem is that you have to hae some condition in at least one of the function under which the calls terminate. Each cycle of calls must, in some sense, be reducing the amount of remaining work, and then at some point the calls have to detect that there is no remaining work and return back values.

However... as outsiders we have to ask whether mypis() really depends upon mypij() and mypij() really depends upon mypis() . Sorry, but at the moment as an outsider I get the impression that you are flailing around trying things without thinking through questions about which routine needs to process vectors and which routines need to process invididual elements.

We cannot make any recommendations for you because there is not enough documentation about what the code is intended to do, or why you have those multiple layers of calls.

Mathieu NOE el 6 de Sept. de 2021

hello again

I don't know what the code is supposed to do (as most of the time when you are not expert in the specific topic)

but if I were you I would first make a graph to show how should the data flow be organized - and avoid this situation.

Try to make the whole process simpler like :

one main function call one or several subfunctions , that can work in parallel but are not dependent one each other
if you really needs several hierarchical levels , a subfunction of level -1 can call a subfunction of level -2 , which can also call a subfunction of level -3 , etc... but the problem that we see in your code is that you create a kind of algebraic loop that cannot be solved : a first subfunction call a second subfunction that call itself the first subfunction. That is simply not good coding practice ...

If I had a lot of free time , I could try rewritte your program , but I need to see how the data flow is / should be organized, where we start, what we need to compute (some infos about the whole process would be needed), etc...

so this can be a lot of work and even if I am also keen to give a hand within the possibilities of my skills and agenda, I will not promise to solve all issues for tomorrow !

Celestine Siameh el 6 de Sept. de 2021

equations matlab.pdf

@Mathieu NOE @Walter Roberson

Thanks for the feedback, I will try and incorporate all the feedback. I agree my model is complicated to solve, that is why I need help from you all, as it is beyond my own coding skills. I realized you both have the same feedback on my code. Moreover, my background is in economics and I am not perfectly excellent at coding, I am still improving on that.

Summary of what I am trying to achieve and the details involved. I want to solve for the price (pis) of output (qis), wages (wis), final demand(fij), intermediate demand(mij).

*The price(pis) is a function of wage(wis) and a function of price index(PMis).

*Now if this output is purchased by another country, then that couuntry will be charged for the price of the output plus the import tariff, this new price then becomes pij. So pij is a function of pis and tariff(tau).

*Now this output purchased is either consumed or use as intermediate input by the country buying it.

*If is used as final consumption, then its demand function is fij, which is a function of pij, PFis, and Yis.

*PFis is similar to PMis, which are all price index, and they depend on pij not pis, but pij depends on pis, and pis depend on wis and PMis.

*Yis is the total income which is a function of wis, total labor endowment (Lis), and total revenue (TRis).

*TRis is the revenue generated from the import tariff imposed by that country purchasing it, so is a function of (pis*fij*tau)+(pis*tau*mij).

*Now if this output is used as intermediate input, then its demand function is mij, which is a function of pij, PMis, and IMis.

*IMis is the total value of the intermediate input used; which is sourced across all sectors of a country, and it is a function of mij and PMis.

*So to solve for the 4 unknown variables, I need to solve the 4 equations in the GEsol, where equation1 solves fij (NS,N) iteratively, equation 2 solves mij (NS,N) iteratively, equation 3 solves pis (NS,1) iteratively and equation 4 solves wis(NS,1) iteralively. But these variables are functions of each other and solves simultaneously.

N is the same Ncoun (countries) and S is the same Nsec (sectors).

So I want to solve one matrix for all these unknown variables then disintegrate them. That is why I have the initial guess as:

% initial guess vakues for x matrix

x0 = 0.5.*ones(8,10);

% solution of the model in terms of x

mij = x(:,1:Ncoun);

fij = x(:,Ncoun+1:end-2);

pis = x(:,end-1);

wis = x(:,end);

*the rows are the same, which is N*S; and the first N cols is for fij, then second N cols is for mij, the last but one col is for pis, and the final col is for wis.

Please I have attached the equations to be solved and I will be happy if you can check it again Maybe I will have to change everything, which I am ready to do.

Hope to hear from you all soon.

Thanks,

Celestine.

Celestine Siameh el 7 de Sept. de 2021

@Mathieu NOE @Walter Roberson

Hello Please can I get help on my coding issue below.