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{:}});
%p
Related Question
- Double integration plot error. messed up with vector and syms
- Vectorisation of a loop
- How can i solve this equation
- I want to convert matlab code to verilog for the image processing project using hdl coder, i have the code but i dont know how to divide the code into function and test bench, please help me. I m using matlab r2018a version.
Best Answer