System of equations with array inputs

Hi,

I am trying to solve a system of 4 equations with 4 unknowns, but two of the inputs are 499x1 arrays. Essentially, I am trying to find the variable values for all 499 input values, and record in a 499x4 matrix. I have the below code (on two separate scripts) but I keep getting an error message (see below). Would anyone be able to guide me on how to approach a problem like this and whether I need to use a loop?

global Yc Yd
function R = eqns_optimtax(X,input)

% Parameters and State Variables
for P=1:499
eps= input(1);
Adt= input(2);
Act = input(3);
omega = input(4);
Yc(P,1) = input(5);
Yd(P,1) = input(6);
end

% Equilibrium equations (note: here, X(1))=wt, X(2)=pct, X(3)=pdt,
% X(4)=tau
% Equation 1: pct*Act-wt=0
% Equation 2: pdt*Adt*Ldt-wt=0
% Equation 3: ((1-omega)*Yct^((eps-1)/eps)+
% omega*Ydt^((eps-1)/eps))^(1/(eps-1))*Yct^(-1/eps)*(1-omega) - pct=0
% Equation 4: ((1-omega)*Yct^((eps-1)/eps)+
% omega*Ydt^((eps-1)/eps))^(1/(eps-1))*Ydt^(-1/eps)*omega - (1+tau)pdt=0

for P=1:499
    R(1,1)= X(2).*Act - X(1);
    R(2,1)= X(3).*Adt- X(1);
    R(3,1)= ((1-omega).*Yc(P,1).^((eps-1)/eps)+ omega.*Yd(P,1).^((eps-1)/eps)).^(1/(eps-1)).*Yc(P,1).^(-1/eps).*(1-omega) -X(2);
    R(4,1)= ((1-omega).*Yc(P,1).^((eps-1)/eps)+ omega.*Yd(P,1).^((eps-1)/eps)).^(1/(eps-1)).*Yd(P,1).^(-1/eps).*omega -(1+X(4)).*X(3);
end

This now runs on a separate script

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%% New Script %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
close all
%This code computes the optimal tax for the baseline model
% 4 equations, 4 unknowns(wt, pct, pdt, tau):
global Yc Yd
% Equation 1: pct*Act-wt=0
% Equation 2: pdt*Adt*Ldt-wt=0
% Equation 3: ((1-omega)*Yct^((eps-1)/eps)+
% omega*Ydt^((eps-1)/eps))^(1/(eps-1))*Yct^(-1/eps)*(1-omega) - pct=0
% Equation 4: ((1-omega)*Yct^((eps-1)/eps)+
% omega*Ydt^((eps-1)/eps))^(1/(eps-1))*Ydt^(-1/eps)*omega - (1+tau)pdt=0

% This is a non-linear system of equations. The code will use fsolve to
% find the solution for wt, pct, pdt,and tau for given parameter values of eps, omega; the given state variables Adt, Act and will match the level of 
% production to the optimal path determined by the social planner's solution.
% The system of equations is written in the other matlab file
% eqns_optimtax.m

% Parameter values and state variables
eps = 3;
Adt = 1;
Act = 1;
omega = 0.6;



% Guess vector of variables
for P=1:499
param(P,1:6)=[eps; Adt; Act;omega; Yc(P,1); Yd(P,1)];
X0(P,1:4)=[1, 0.2903, 0.7097, 0.8];
end

% Using fsolve to compute values of interest. The vector X is the solution
[X(P,1:4),error_of_solution]=fsolve('eqns_optimtax', X0(P,1:4) ,[], param(P,1:6));


% Report the largerst departure from 0 in absolute value
max_error=max(abs(error_of_solution));

Error Message:

All functions in a script must be closed with an 'end'.

Error in fsolve (line 260)

fuser = feval(funfcn{3},x,varargin{:});

Error in optimtax (line 33)

[X(P,1:4),error_of_solution]=fsolve('eqns_optimtax', X0(P,1:4) ,[], param(P,1:6));

Answer 1

[Torsten]

function main
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %%%%%%%%%%%%%%%%%%%%%%%%%% New Script %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%This code computes the optimal tax for the baseline model
    % 4 equations, 4 unknowns(wt, pct, pdt, tau):
    global Yc Yd
    % Equation 1: pct*Act-wt=0
    % Equation 2: pdt*Adt*Ldt-wt=0
    % Equation 3: ((1-omega)*Yct^((eps-1)/eps)+
    % omega*Ydt^((eps-1)/eps))^(1/(eps-1))*Yct^(-1/eps)*(1-omega) - pct=0
    % Equation 4: ((1-omega)*Yct^((eps-1)/eps)+
    % omega*Ydt^((eps-1)/eps))^(1/(eps-1))*Ydt^(-1/eps)*omega - (1+tau)pdt=0
    % This is a non-linear system of equations. The code will use fsolve to
    % find the solution for wt, pct, pdt,and tau for given parameter values of eps, omega; the given state variables Adt, Act and will match the level of 
    % production to the optimal path determined by the social planner's solution.
    % The system of equations is written in the other matlab file
    % eqns_optimtax.m
    % Parameter values and state variables
    eps = 3;
    Adt = 1;
    Act = 1;
    omega = 0.6;
    % Guess vector of variables
    x0 = [1, 0.2903, 0.7097, 0.8];
    X = zeros(4,499);
    % Using fsolve to compute values of interest. The vector X is the solution
    for P = 1:499 
        param =[eps; Adt; Act;omega; Yc(P,1); Yd(P,1)];
        [x,error_of_solution]=fsolve(@(x)eqns_optimtax(x,param), x0);
        X(:,P) = x;
        ERROR_OF_SOLUTION(P) = norm(error_of_solution);
        x0 = x;
    end
    % Report the largerst departure from 0 in absolute value
    max_error = max(ERROR_OF_SOLUTION);
    plot(1:499,X(1,:))
end

function R = eqns_optimtax(X,input)
    % Parameters and State Variables
        
    eps= input(1);
    Adt= input(2);
    Act = input(3);
    omega = input(4);
    Yc = input(5);
    Yd = input(6);
        
    % Equilibrium equations (note: here, X(1))=wt, X(2)=pct, X(3)=pdt,
    % X(4)=tau
    % Equation 1: pct*Act-wt=0
    % Equation 2: pdt*Adt*Ldt-wt=0
    % Equation 3: ((1-omega)*Yct^((eps-1)/eps)+
    % omega*Ydt^((eps-1)/eps))^(1/(eps-1))*Yct^(-1/eps)*(1-omega) - pct=0
    % Equation 4: ((1-omega)*Yct^((eps-1)/eps)+
    % omega*Ydt^((eps-1)/eps))^(1/(eps-1))*Ydt^(-1/eps)*omega - (1+tau)pdt=0
        
    R(1,1)= X(2)*Act - X(1);
    R(2,1)= X(3)*Adt- X(1);
    R(3,1)= ((1-omega)*Yc^((eps-1)/eps)+ omega*Yd^((eps-1)/eps))^(1/(eps-1))*Yc^(-1/eps)*(1-omega) -X(2);
    R(4,1)= ((1-omega)*Yc^((eps-1)/eps)+ omega*Yd^((eps-1)/eps))^(1/(eps-1))*Yd^(-1/eps)*omega -(1+X(4))*X(3);
end

[Irem Sara]

Thank you very much! When I run the code, I can see the equation is being solved in each step but unfortunately X (499x4 matrix with the solutions) isn't in the variable list. Do you know how to deal with this?

[Torsten]

Yes. Remove the line

function main

and the

end

at the end of the function.

[Irem Sara]

Great- thanks a lot!

Answer 2

[Jan]

Your file eqns_optimtax.m starts with the line:

global Yc Yd

An m file starting with code instead of the keyword "function" is a script, not a function. If you define function inside scripts, they must be closed by an "end", which is missing.

I guess, that the line "globak Yc Yd" is complete orphaned. Simply delete it because it does not make any sense here. Of course eqns_optimtax must be a function, not a script.

Prefer to use function handles in the call of fsolve. Using a char vector containing the name of the function is still supported, but outdated for 20 years now.

[X(P,1:4),error_of_solution]=fsolve(@eqns_optimtax, X0(P,1:4) ,[], param(P,1:6));

Hint: Clean up the loop.

for P=1:499
    eps= input(1);
    Adt= input(2);
    Act = input(3);
    omega = input(4);
    Yc(P,1) = input(5);
    Yd(P,1) = input(6);
end

% Smarter without a loop:
eps= input(1);
Adt= input(2);
Act = input(3);
omega = input(4);
Yc(1:499,1) = input(5);
Yd(1:499,1) = input(6);

Using "input" and "end" as names of variables causes unexpected behavior frequently, because this shadows important built-in functions.