"solve" cannot find analytic solution to a large multi-equation system, unless one is willing to break up the system and solve the parts sequentially

Greetings,

I had this issue when I tried to use the Symbolic Toolkit (MATLAB R2009b on Windows PC and R2010a on Linux) to find an analytic solution to a system of equations.

In short (more background info and code below), there are 72 equation in 72 variables in my system. Many of the equations are trivial (x=0), the rest of the system has polynomials in the variables of at most second-order. I know that an analytic solution exists (actually 2 solutions). When I write

solve(equations,variables);

MATLAB says that it cannot find explicit solution (I get a warning on line 81 of solve). Since the system is really simple, I can go on and try to solve a subset of the equations (40 out of the 72):

solve(some_equations,some_variables);

gives a solution, which I substitute back into the remaining 32 eqns. Then, I solve a subset (8) of these 32 eqns, plug this solution back into the remaining 24 eqns, solve a subset of these eqns, and so on until I end up with the correct analytical solutions to the original problem.

Any idea what I should do differently to have MATLAB solve the large system right away? This is a fairly simple problem that I wanted to use to learn how to solve these types of problems in MATLAB. In other words, I'm hoping that in the future I could solve more complicated problems, where "substituting back in manually" is a lot more time consuming.

Background: I have a system of linear dynamic equations: 9 equations in 6 control vars and 3 endogenous state vars plus 5 exogenous forcing processes.

I want to find the state-space representation of this system using the "method of undetermined coefficients". Since the system is linear I guess that the state-space representation will be linear as well. So I have a linear guess for the endogenous variables (control and state) of the form:

x = ETA*y, where x is a vector of endogenous variables [control vars in time t; state vars in time t+1] y is a vector of state variables [state vars in time t; forcing processes in time t].

I plug these back into the original 9 equations, and look for the ETAs that satisfy the resulting equations. This is the 72(=(3+5)x9) equations in the 72 ETAs.

Thanks for your time and effort!

Best, Tamas

Source code:

%--------------------------------------------------------------------------
% Devereux and Sutherland (2007 IMF Working Paper)
% "Solving for Country Portfolios in Open Economy Macro Models"
% Example 1
%--------------------------------------------------------------------------
clear all
clc
% tic
%--------------------------------------------------------------------------
% Define variables s_t, x_t, c_t (see middle of p.21)
%--------------------------------------------------------------------------
endog_states = ['LW   ';'LPE  ';'LPEF '] ; 
exog_states  = ['Y    ';'YF   ';'M    ';'MF   '] ;
control_vars = ['C ';'CF';'rE';'P ';'PF';'YD'];
state_space  = [endog_states; exog_states;'XI   '] ;

syms LW LPE LPEF Y YF M MF XI LY LYF LM LMF epsY epsYF epsM epsMF...
     B RO zeta_Y zeta_M sigma_Y sigma_M ALPHA_tilde real
VARS = [LW LPE LPEF Y YF M MF XI];

%--------------------------------------------------------------------------
% Obtain state-space representation of system (eq. 20 and 21 on pp. 14-15)
%--------------------------------------------------------------------------
% Shock equations
%Shocks = cell(size(exog_states,1),1) ;
y  = 'zeta_Y'*Y  + epsY  ; ey  = 'zeta_Y'*Y  ;
yf = 'zeta_Y'*YF + epsYF ; eyf = 'zeta_Y'*YF ;
m  = 'zeta_M'*M  + epsM  ; em  = 'zeta_M'*M  ;
mf = 'zeta_M'*MF + epsMF ; emf = 'zeta_M'*MF ;

% State space representation for endogenous state variables
Next_states = cell(size(endog_states,1),1);
etas = cell(1,length(control_vars)*length(state_space)) ;
for i = 1:size(endog_states,1)
f = 0;
coeff = cell(size(state_space,1),1);
    for j = 1:size(state_space,1)
        coeff{j} = ['eta_' deblank(endog_states(i,:)) '_' deblank(state_space(j,:))] ;
        etas{1,(i-1)*size(state_space,1) + j}  = coeff{j} ; 
        f            = f + coeff{j}*VARS(j) ;
        clear temp2
    end

Next_states{i}  = f; 
end

% State space representation for exogenous variables
Control_vars = cell(size(control_vars,1),1);
for i = 1:size(control_vars,1)
f = 0;
coeff = cell(size(state_space,1),1);
    for j = 1:size(state_space,1)
        coeff{j} = ['eta_' deblank(control_vars(i,:)) '_' deblank(state_space(j,:))] ;
        etas{1,size(endog_states,1)*size(state_space,1) + ...
            (i-1)*size(state_space,1) + j}  = coeff{j} ;
        f     = f + coeff{j}*VARS(j) ;
        clear temp2
    end

Control_vars{i}  = f;  
end
clear coeff

% Summarizing the results (all variables are meant to be deviations from ss)
%------------------------
% Wealth at the beginning of period t (W_t)
Wt = Next_states{1};
% pretty(collect(Wt,[LW LPE LPEF Y YF M MF XI]))
% E_{t-1}P_t
Pet = Next_states{2};
% pretty(collect(Pet,[LW LPE LPEF Y YF M MF XI]))
% E_{t-1}P^*_t
Peft =Next_states{3};
% pretty(collect(Peft,[LW LPE LPEF Y YF M MF XI]))
% Home consumption in period t (c_t)
ct = Control_vars{1};
% pretty(collect(ct,[LW LPE LPEF Y YF M MF XI]))
% Foreign consumption in period t (c^f_t)
cft = Control_vars{2};
% pretty(collect(cft,[LW LPE LPEF Y YF M MF XI]))
% Expected returns in period t (r^E_t)
ret = Control_vars{3};
% pretty(collect(ret,[LW LPE LPEF Y YF M MF XI]))
% Price level in period t (P_t)
pt = Control_vars{4};
% pretty(collect(pt,[LW LPE LPEF Y YF M MF XI]))
% Foreign price level in period t (P^*_t)
pft = Control_vars{5};
% pretty(collect(pft,[LW LPE LPEF Y YF M MF XI]))
% Output differential between countries
ydt = Control_vars{6};
% pretty(collect(ydt,[LW LPE LPEF Y YF M MF XI]))
%--------------------------------------------------------------------------
% Substitute into the log-linearized equations (eqs. 46, 48, 49, 50, 53)
%--------------------------------------------------------------------------

% Home Euler-equation
%--------------------
% Expectation of c_{t+1} in period t (E_t c_{t+1}):
% There is no economic intuition in the next equation, it is just the
% definition of E_t c_{t+1} using the objects created above...
E_ct = subs(ct,{'LW','LPE','LPEF','Y','YF','M','MF','XI'},...
    {Wt,Pet,Peft,ey,eyf,em,emf,0});
f_Euler = collect(RO*(ct - E_ct) + ret,[LW, LPE, LPEF, Y, YF, M, MF, XI]) ;
% pretty(f_Euler);

% Collecting the coefficients on the the state variables and putting them
% into a structure.  Since the model is linearized, the first derivatives
% with respect to the state variables will give the coefficient on the
% respective state variable. 
% These coefficient along with the ones from the following equations will 
% have to be set to zero, to obtain the solution of the problem, that is 
% to get matrices F_1, F_2, F_3, P_1, P_2, P_3 in eq. 21 on page 15.

diff_now = [LW, LPE, LPEF, Y, YF, M, MF, XI];
for k = 1:size(state_space)
    coeff = diff(f_Euler,diff_now(k),1);
    undet_c{k} = [char(coeff) ' = 0'];
end

% Foreign Euler-equation
%-----------------------
E_cft = subs(cft,{'LW','LPE','LPEF','Y','YF','M','MF','XI'},...
    {Wt,Pet,Peft,ey,eyf,em,emf,0});

f_EulerF = collect(RO*(cft - E_cft) + ret,[LW, LPE, LPEF, Y, YF, M, MF, XI]) ;
% pretty(f_EulerF);

% Continuing with the 'undet_c' structure...
diff_now = [LW, LPE, LPEF, Y, YF, M, MF, XI];
temp = length(undet_c);
for k = 1:size(state_space)
    coeff = diff(f_EulerF,diff_now(k),1);
    undet_c{temp+k} = [char(coeff) ' = 0'];
end

% Resource constraint
%---------------------
f_resource = collect(Y + YF - ct - cft,[LW, LPE, LPEF, Y, YF, M, MF, XI]) ;
% pretty(f_resource)

diff_now = [LW, LPE, LPEF, Y, YF, M, MF, XI];
temp = length(undet_c);
for k = 1:size(state_space)
    coeff = diff(f_resource,diff_now(k),1);
    undet_c{temp+k} = [char(coeff) ' = 0'];
end

% Budget constraint
%-------------------
B = sym('B');
f_budget = collect(LW/B + Y - ct + XI - Wt,...
    [LW, LPE, LPEF, Y, YF, M, MF, XI]) ;
% pretty(f_budget)

diff_now = [LW, LPE, LPEF, Y, YF, M, MF, XI];
temp = length(undet_c);
for k = 1:size(state_space)
    coeff = diff(f_budget,diff_now(k),1);
    undet_c{temp+k} = [char(coeff) ' = 0'];
end

% Quantity theory
%-----------------
f_quant = collect(Y - M + pt,[LW, LPE, LPEF, Y, YF, M, MF, XI]) ;
% pretty(f_quant)

diff_now = [LW, LPE, LPEF, Y, YF, M, MF, XI];
temp = length(undet_c);
for k = 1:size(state_space)
    coeff = diff(f_quant,diff_now(k),1);
    undet_c{temp+k} = [char(coeff) ' = 0'];
end

% Quantity theory foreign
%-------------------------
f_quantF = collect(YF - MF + pft,[LW, LPE, LPEF, Y, YF, M, MF, XI]) ;
% pretty(f_quantF)

diff_now = [LW, LPE, LPEF, Y, YF, M, MF, XI];
temp = length(undet_c);
for k = 1:size(state_space)
    coeff = diff(f_quantF,diff_now(k),1);
    undet_c{temp+k} = [char(coeff) ' = 0'];
end

% Equations 54
%-------------
%E_pt = subs(pt,{'Y','YF','M','MF','XI'},{ey,eyf,em,emf,0});
f_pexp = collect(pt - Pet,[LW, LPE, LPEF, LY, LYF, LM, LMF, XI]);
%subs(pe,{'Y','YF','M','MF','XI'},{y,yf,m,mf}) ;
% pretty(f_pexp)

diff_now = [LW, LPE, LPEF, Y, YF, M, MF, XI];
temp = length(undet_c);
for k = 1:size(state_space)
    coeff = diff(f_pexp,diff_now(k),1);
    undet_c{temp+k} = [char(coeff) ' = 0'];
end

f_pexpf = collect(pft - Peft,[LW, LPE, LPEF, Y, YF, M, MF, XI]) ;

diff_now = [LW, LPE, LPEF, Y, YF, M, MF, XI];
temp = length(undet_c);
for k = 1:size(state_space)
    coeff = diff(f_pexpf,diff_now(k),1);
    undet_c{temp+k} = [char(coeff) ' = 0'];
end

% Output differential: Y_t - Y^*_t - Y_{D,t} = 0
%------------------------------------------------
f_yd = collect(Y - YF - ydt);
% pretty(f_yd)

diff_now = [LW, LPE, LPEF, Y, YF, M, MF, XI];
temp = length(undet_c);
for k = 1:size(state_space)
    coeff = diff(f_yd,diff_now(k),1);
    undet_c{temp+k} = [char(coeff) ' = 0'];
end

%%Solving the resulting equations
%----------------------------------
% This should ideally solve it...
% eta_sol = solve(undet_c{:},etas{:});

% Some equations give very trivial solutions, these are collected below
vals = cell(1,length(etas));
sol = solve(undet_c{end-39:end},etas{9:24},etas{end-23:end});
for k = 1:16
    vals{8+k} = sol.(etas{8+k});
end
for k = 1:24
   vals{end-24+k} = sol.(etas{end-24+k}); 
end

% From budget constraint get eta_LW_ as a function of eta_C_
diff_now = [LW, LPE, LPEF, Y, YF, M, MF, XI];
for k = 1:size(diff_now,2)
    coeff = diff(f_budget,diff_now(k),1);
    temp_eq{k} = [char(coeff) ' = 0'];
end
sol = solve(temp_eq{1:8},etas{1:8});
for k = 1:8
    temp_LW{k} = sol.(etas{k});
end

% from resource constraint get eta_cf_ coeffs as a function of eta_c coeffs
sol = solve(undet_c{17:24},etas{33:40});
for k = 1:8
    temp_CF{k} = sol.(etas{32+k});
end

% Plug in temp_CF and temp_LW into foreign Euler to get eta_re_ as a
% function of eta_C_
f_EF2 = subs(f_EulerF,{'eta_LW_LW','eta_LW_LPE','eta_LW_LPEF','eta_LW_Y','eta_LW_YF',...
    'eta_LW_M','eta_LW_MF','eta_LW_XI','eta_CF_LW','eta_CF_LPE','eta_CF_LPEF',...
    'eta_CF_Y','eta_CF_YF','eta_CF_M','eta_CF_MF','eta_CF_XI',...
    'eta_LPE_LW','eta_LPE_LPE','eta_LPE_LPEF','eta_LPE_Y','eta_LPE_YF',...
    'eta_LPE_M','eta_LPE_MF','eta_LPE_XI',...
    'eta_LPEF_LW','eta_LPEF_LPE','eta_LPEF_LPEF','eta_LPEF_Y','eta_LPEF_YF',...
    'eta_LPEF_M','eta_LPEF_MF','eta_LPEF_XI'},{temp_LW{:} ...
    temp_CF{:} vals{9:24}}) ;
% pretty(collect(f_EF2,[LW, LPE, LPEF, Y, YF, M, MF, XI]))

diff_now = [LW, LPE, LPEF, Y, YF, M, MF, XI];
for k = 1:size(diff_now,2)
    coeff = diff(f_EF2,diff_now(k),1);
    temp_eq{k} = [char(coeff) ' = 0'];
end
sol = solve(temp_eq{1:8},etas{41:48});
for k = 1:8
    temp_RE{k} = sol.(etas{40+k});
end

% Plug all of the above into home Euler eq
f_E2 = subs(f_Euler,{'eta_LW_LW','eta_LW_LPE','eta_LW_LPEF','eta_LW_Y','eta_LW_YF',...
    'eta_LW_M','eta_LW_MF','eta_LW_XI','eta_CF_LW','eta_CF_LPE','eta_CF_LPEF',...
    'eta_CF_Y','eta_CF_YF','eta_CF_M','eta_CF_MF','eta_CF_XI',...
    'eta_LPE_LW','eta_LPE_LPE','eta_LPE_LPEF','eta_LPE_Y','eta_LPE_YF',...
    'eta_LPE_M','eta_LPE_MF','eta_LPE_XI',...
    'eta_LPEF_LW','eta_LPEF_LPE','eta_LPEF_LPEF','eta_LPEF_Y','eta_LPEF_YF',...
    'eta_LPEF_M','eta_LPEF_MF','eta_LPEF_XI'...
    'eta_rE_LW','eta_rE_LPE','eta_rE_LPEF','eta_rE_Y','eta_rE_YF',...
    'eta_rE_M','eta_rE_MF','eta_rE_XI',},...
    {temp_LW{:} temp_CF{:} vals{9:24} temp_RE{:}});
%pretty(collect(f_E2,[LW, LPE, LPEF, Y, YF, M, MF, XI]))

diff_now = [LW, LPE, LPEF, Y, YF, M, MF, XI];
for k = 1:size(diff_now,2)
    coeff = diff(f_E2,diff_now(k),1);
    temp_eq{k} = [char(coeff) ' = 0'];
end
sol = solve(temp_eq{1:8},etas{25:32});
for k = 1:8
    vals{24+k} = sol.(etas{24+k})(2);
end
% Substitute back solutions
% eta_RE_
for k = 1:8
vals{40+k} = simplify(subs(temp_RE{k},{'eta_C_LW','eta_C_LPE','eta_C_LPEF','eta_C_Y','eta_C_YF',...
    'eta_C_M','eta_C_MF','eta_C_XI'},{vals{25:32}}));
%display(vals{40+k})
end

% eta_CF_
for k = 1:8
vals{32+k} = simplify(subs(temp_CF{k},{'eta_C_LW','eta_C_LPE','eta_C_LPEF','eta_C_Y','eta_C_YF',...
    'eta_C_M','eta_C_MF','eta_C_XI'},{vals{25:32}}));
%display(vals{32+k})
end

% eta_LW_
for k = 1:8
vals{k} = simplify(subs(temp_LW{k},{'eta_C_LW','eta_C_LPE','eta_C_LPEF','eta_C_Y','eta_C_YF',...
    'eta_C_M','eta_C_MF','eta_C_XI'},{vals{25:32}}));
%display(vals{k})
end

Answer 1

[Stefan Wehmeier]

Hi Tamas, this works in newer versions of the symbolic toolbox (R2011a) as we have added some heuristics for systems containing many trivial equations. What can happen if you have many variables and many equations is an internal timeout because not enough "effort" is left. We prefer to give up soon, rather than letting the user wait for hours. If you want to experiment with your systems and are prepared to wait if necessary, just allow the solver to spend much more "effort". Example: evalin(symengine, 'MAXEFFORT:= 10^35')

In your example, we used to overestimate the necessary effort to solve the system in previous versions.

Best regards,

Stefan Wehmeier (symbolic toolbox developer)

[Tamas]

Thanks so much for the prompt answer, Stefan! I figured that there should be an easy fix to this issue. It works fine with the evalin(symengine, 'MAXEFFORT:= 10^35') command.

Best,

Tamas