How to use fsolve to solve matrix equations involving matrices of different dimensions?

I have the following equation to solve using fsolve: (all elements in matrices A0, A1, A2, A3, B1, B2 are given parameter values)

A0 = [1 0 sigma; -kappa 1 0; -phi_x -phi_pi 1];

A1 = [1 sigma 0; 0 beta 0; 0 0 0];

A2 = [1 0; 0 1; 0 0];

A3 = [0 0 0; 0 0 0; sigma_i 0 0];

B1 = [rho_g 0, 0 rho_u];

B2 = [0 sigma_g 0; 0 0 sigma_u];

F0 = fsolve(@(x) A1*x*B1 - A0*x + A2, zeros(3,2));

However, I get the following error:

Error using *

Incorrect dimensions for matrix multiplication. Check that the number of columns in the first matrix matches

the number of rows in the second matrix. To operate on each element of the matrix individually, use TIMES (.*)

for elementwise multiplication.

Error in Standard_NK_economy>@(x)A1(ir,:)*x*B1(:,ir)-A0(ir,:)*x+A2(ir,:) (line 29)

f0 = fsolve(@(x) A1(ir,:)*x*B1(:,ir) - A0(ir,:)*x + A2(ir,:), zeros(1,2)); % x is 1x2

Error in fsolve (line 270)

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

Error in Standard_NK_economy (line 29)

f0 = fsolve(@(x) A1(ir,:)*x*B1(:,ir) - A0(ir,:)*x + A2(ir,:), zeros(1,2)); % x is 1x2

Caused by:

Failure in initial objective function evaluation. FSOLVE cannot continue.

Related documentation

[Dyuman Joshi]

Please share the mathematical formulation of the equations you are trying to solve.

There are undefined variables in the code above, thus we can not run your code.

[Smriti Verma]

% Calibration

sigma = 1;

kappa = 0.2;

beta = 0.99;

phi_pi = 1.5;

phi_x = 0.5;

sigma_i = 0.0025;

sigma_g = 0.0025;

sigma_u = 0.0025;

rho_g = 0.5;

rho_u = 0.3;

% Define A0, A1, A2, A3 for the 3 DSGE equations, where Z = [x; pi; i], S =

% [g; u], and eps = [eps_i; eps_g; eps_u]

A0 = [1 0 sigma; -kappa 1 0; -phi_x -phi_pi 1];

A1 = [1 sigma 0; 0 beta 0; 0 0 0];

A2 = [1 0; 0 1; 0 0];

A3 = [0 0 0; 0 0 0; sigma_i 0 0];

% Define B1, B2 for the two state equations: S_t = B1*S_t-1 + B2*eps_t

B1 = [rho_g 0, 0 rho_u];

B2 = [0 sigma_g 0; 0 0 sigma_u];

% Guessing Z_t = F0*S_t + F1*eps_t, we calculate F0 and F1:

F0 = fsolve(@(x) A1*x*B1 - A0*x + A2, zeros(3,2));

F1 = A0\A3; % inv(A0)*A3

G1 = F0*B1;

G2 = F0*B2 + F1;

Answer 1

[Torsten]

Your equations are linear in the unknowns - so you should use a linear solver, not a nonlinear one like "fsolve":

% Calibration
sigma = 1;
kappa = 0.2;
beta = 0.99;
phi_pi = 1.5;
phi_x = 0.5;
sigma_i = 0.0025;
sigma_g = 0.0025;
sigma_u = 0.0025;
rho_g = 0.5;
rho_u = 0.3;
% Define A0, A1, A2, A3 for the 3 DSGE equations, where Z = [x; pi; i], S =
% [g; u], and eps = [eps_i; eps_g; eps_u]
A0 = [1 0 sigma; -kappa 1 0; -phi_x -phi_pi 1];
A1 = [1 sigma 0; 0 beta 0; 0 0 0];
A2 = [1 0; 0 1; 0 0];
A3 = [0 0 0; 0 0 0; sigma_i 0 0];
% Define B1, B2 for the two state equations: S_t = B1*S_t-1 + B2*eps_t
B1 = [rho_g 0; 0 rho_u];
B2 = [0 sigma_g 0; 0 0 sigma_u];
X = sym('X',[size(A1,2),size(B1,1)]);
size(X)

ans = 1×2

3 2

eqns = A1*X*B1 - A0*X + A2 == 0;
symvar(eqns)

ans = 

[A,b] = equationsToMatrix(eqns)

A = 

b = 

F0 = A\b

F0 = 

F0 = [F0(1),F0(2);F0(3),F0(4);F0(5),F0(6)];
F1 = A0\A3 % inv(A0)*A3

F1 = 3×3

-0.0014 0 0 -0.0003 0 0 0.0014 0 0

G1 = F0*B1

G1 = 

G2 = F0*B2 + F1

G2 = 

[Dyuman Joshi]

Note - Requires Symbolic Math Toolbox.

[Torsten]

Seems it also works this way:

% Calibration
sigma = 1;
kappa = 0.2;
beta = 0.99;
phi_pi = 1.5;
phi_x = 0.5;
sigma_i = 0.0025;
sigma_g = 0.0025;
sigma_u = 0.0025;
rho_g = 0.5;
rho_u = 0.3;
% Define A0, A1, A2, A3 for the 3 DSGE equations, where Z = [x; pi; i], S =
% [g; u], and eps = [eps_i; eps_g; eps_u]
A0 = [1 0 sigma; -kappa 1 0; -phi_x -phi_pi 1];
A1 = [1 sigma 0; 0 beta 0; 0 0 0];
A2 = [1 0; 0 1; 0 0];
A3 = [0 0 0; 0 0 0; sigma_i 0 0];
% Define B1, B2 for the two state equations: S_t = B1*S_t-1 + B2*eps_t
B1 = [rho_g 0; 0 rho_u];
B2 = [0 sigma_g 0; 0 0 sigma_u];
F0 = fsolve(@(x) A1*x*B1 - A0*x + A2, zeros(3,2));

Equation solved.

fsolve completed because the vector of function values is near zero
as measured by the value of the function tolerance, and
the problem appears regular as measured by the gradient.

F1 = A0\A3 % inv(A0)*A3

F1 = 3×3

-0.0014 0 0 -0.0003 0 0 0.0014 0 0

G1 = F0*B1

G1 = 3×2

0.3582 -0.3322 0.1418 0.3322 0.3918 0.3322

G2 = F0*B2 + F1

G2 = 3×3

-0.0014 0.0018 -0.0028 -0.0003 0.0007 0.0028 0.0014 0.0020 0.0028

The error you made in your code was that you specified B1 as a 1x4 vector instead of a 2x2 matrix:

B1 = [rho_g 0, 0 rho_u];

instead of

B1 = [rho_g 0; 0 rho_u];

[Smriti Verma]

Right, of course, thanks so much!

[Smriti Verma]

I am trying to write code for solving the New Keynesian model steady state and transition dynamics, is it okay to ask about it in this thread?

[Dyuman Joshi]

Best to ask a new question.

[Smriti Verma]

I have a posted a question about modeling in Economics, would you have any idea by any chance? I could provide the link to the question.

Answer 2

[Dyuman Joshi]

% Calibration
sigma = 1;
kappa = 0.2;
beta = 0.99;
phi_pi = 1.5;
phi_x = 0.5;
sigma_i = 0.0025;
sigma_g = 0.0025;
sigma_u = 0.0025;
rho_g = 0.5;
rho_u = 0.3;
% Define A0, A1, A2, A3 for the 3 DSGE equations, where Z = [x; pi; i], S =
% [g; u], and eps = [eps_i; eps_g; eps_u]
A0 = [1 0 sigma; -kappa 1 0; -phi_x -phi_pi 1];
A1 = [1 sigma 0; 0 beta 0; 0 0 0];
A2 = [1 0; 0 1; 0 0];
A3 = [0 0 0; 0 0 0; sigma_i 0 0];
% Define B1, B2 for the two state equations: S_t = B1*S_t-1 + B2*eps_t
B1 = [rho_g 0; 0 rho_u];
B2 = [0 sigma_g 0; 0 0 sigma_u];

%size of x
r = 3;
c = 2;
F0 = fsolve(@(x) fun(x, A0, A1, A2, B1, r, c), zeros(r, c))

Equation solved.

fsolve completed because the vector of function values is near zero
as measured by the value of the function tolerance, and
the problem appears regular as measured by the gradient.

F0 = 3×2

0.7163 -1.1074 0.2837 1.1074 0.7837 1.1074

F1 = A0\A3 % inv(A0)*A3

F1 = 3×3

-0.0014 0 0 -0.0003 0 0 0.0014 0 0

G1 = F0*B1

G1 = 3×2

0.3582 -0.3322 0.1418 0.3322 0.3918 0.3322

G2 = F0*B2 + F1

G2 = 3×3

-0.0014 0.0018 -0.0028 -0.0003 0.0007 0.0028 0.0014 0.0020 0.0028

%% Function definition
%Define equations to be solved as a function
function y = fun(x, A0, A1, A2, B1, r, c)
x = reshape(x, r, c); 
y = A1*x*B1 - A0*x + A2;
end

[Smriti Verma]

Thanks!