MATLAB Answers

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)

119 views (last 30 days)
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 Comments
Celestine Siameh
Celestine Siameh on 19 Sep 2021
Tested code to go at the end of your current code:
Does that mean same m file. but the tested code will be at the bottom of it. Thanks.

Sign in to comment.

Accepted Answer

Walter Roberson
Walter Roberson on 19 Sep 2021
Consolidated code, as you are having difficulty understanding me:
I do not recommend running this with maximum function evaluations set to 1e8, which is what I am currently running on my system. The display of all those iterations slows the program down a lot.
1e7 function evaluations took about half an hour on my system. I estimate it would take closer to an hour on your system, based upon the times you reported earlier.
  51 Comments