I have the following script where I am trying to solve a system of equations symbolically using sympy:
import sympy as sp
from sympy import Eq, Function, Symbol, symbols, solve, simplify, collect, pprint
from sympy import Sum, latex
## TODO Create pprinter function
def qpprint(lhs, rhs=0, simplify = True):
expr = Eq(lhs, rhs)
if simplify:
expr = expr.simplify()
return pprint(expr)
pi, g_z = symbols('pi g_{z}', positive = True, nonzero = True)
v, un, gamma_RD, iota = symbols('nu mu gamma_{RD} iota', positive = True, nonzero = True)
i, t_p, t_w, Gamma, B_Y, G_Y = symbols('i tau_{pi} tau_{w} Gamma B_{Y} G_{Y}' , positive = True, nonzero = True)
a_w, a_p, wnom, m_h, b_h = symbols('alpha_w alpha_p W m_h beta_h', positive = True, nonzero = True)
eta_f, eta_b = symbols('eta_f eta_b', positive = True, nonzero = True)
r = symbols('r', positive = True, nonzero = True)
theta = symbols('theta', positive = True, nonzero = True)
j = symbols('j') ## For iteration
lambda_h, lambda_f = symbols('lambda_h lambda_f', positive = True, nonzero = True)
Nf = symbols('N', positive = True, nonzero = True)
levh, levf = symbols('Lambda_h Lambda_f', positive = True, nonzero = True)
D, g_ss = symbols('Delta g_ss', nonzero = True)
ra_h, ra_f, sumAm_h, sumAm_f = symbols('ra_h ra_f sumAm_h sumAm_f')
NWh, NWf, NWb, NWg = symbols('NW_h NW_f NW_b NW_g')
SAVh, SAVf, SAVb, SAVg = symbols('SAV_h SAV_f SAV_b SAV_g')
Bh, Bb, Bg = symbols('B_h B_b B_g')
Mh, Mf, Mb = symbols('M_h M_f M_b')
Lh, Lf, Lb = symbols('L_h L_f L_b')
K, INV = symbols('K INV')
AM_h, AM_f, AM_b = symbols('AM_h AM_f AM_b')
FT_f, FD_f, FU_f, FN_f = symbols('FT_f FD_f FU_f FN_f')
FT_b, FD_b, FU_b = symbols('FT_b FD_b FU_b')
GDPNom, GDP, Yc, C, If, G, Z = symbols('GDPNom GDP Yc C If G Z')
Y_FC = symbols('Y_FC')
Pc, Pk, ULC = symbols('Pc, Pk ULC')
Sc, Sk, RD = symbols('Sc Sk RD')
ib_h, ib_f, spread = symbols('ib_h ib_f spread')
Th, Tf, Tb, Tg = symbols('T_h T_f T_b T_g')
WB, YD_h, YG_h, YT_h, YDk_h, YK_h = symbols('WB YD_h YG_h YT_h YDk_h YK_h')
A = symbols('A')
RDNom = symbols('RDNom')
h, ht = symbols('h h_T')
DepNom, Dep, delta = symbols("DepNom Dep delta")
_D = g_ss/(1 + g_ss) - D
qpprint(_D)
_g_ss = (1 + g_z)*(1 + pi) - 1 - g_ss
qpprint(_g_ss)
_sumAm_h = Sum((1/(1 + g_ss))**(j), (j, 1, lambda_h + 1)).doit() - sumAm_h ## TODO Check necessity of +1
qpprint(_sumAm_h)
_ra_h = (g_ss * sumAm_h)/(lambda_h - sumAm_h) - ra_h
qpprint(_ra_h)
_sumAm_f = Sum((1/(1 + g_ss))**(j), (j, 1, lambda_f + 1)).doit() - sumAm_f ## TODO Check necessity of +1
qpprint(_sumAm_f)
_ra_f = (g_ss * sumAm_f)/(lambda_f - sumAm_f) - ra_f
qpprint(_ra_f)
_GDPNom = (C + G + If)*Pc + RD*ULC - GDPNom
qpprint(_GDPNom)
_Y_FC = (K/(1 + g_z))*(1 / v) - Y_FC
qpprint(_Y_FC)
_un = ((Yc)/Y_FC) - un
qpprint(_un)
_RD = RDNom/ULC - RD
qpprint(_RD)
_RDNom = gamma_RD*(Pc * Sc)/(1 + g_ss) - RDNom
qpprint(_RDNom)
_Pc = (1 + theta)*ULC - Pc
qpprint(_Pc)
_ULC = (wnom)/A - ULC
qpprint(_ULC)
_INV = iota * Sc - INV
qpprint(_INV)
__Yc = Sc + (g_z/(1 + g_z))*INV - Yc
qpprint(__Yc)
_Sc = C + If + G - Sc
qpprint(_Sc)
_If = ht*Yc - If
qpprint(_If)
_ht = h + (DepNom/Pc)/(K) - ht
qpprint(_ht)
_h = (g_z*v)/un - h
qpprint(_h)
_Nf = (Yc / A) + RD / A - Nf
qpprint(_Nf)
_G = G_Y * Yc/(1 + g_z) - G
qpprint(_G)
_Tg = Tf + Th + Tb - Tg
qpprint(_Tg)
_NWg = - Bg - NWg
qpprint(_NWg)
_SAVg = D*NWg - SAVg
qpprint(_SAVg)
__SAVg = -Pc * G + Tg - i * (Bg / (1 + g_ss)) - SAVg
qpprint(__SAVg)
_SFCg = _SAVg - __SAVg
qpprint(_SFCg)
_Bg = B_Y * GDPNom - Bg
qpprint(_Bg)
_NWb = Bb + Lb - Mb - NWb
qpprint(_NWb)
_Mb = Mh + Mf - Mb
qpprint(_Mb)
_Lb = Lh + Lf - Lb
qpprint(_Lb)
_Bb = NWb - Lb + Mb - Bb
# _Bb = (SAVb/D) - Lb + Mb - Bb
qpprint(_Bb)
_SAVb = FU_b - SAVb
qpprint(_SAVb)
__SAVb = D*NWb - SAVb
qpprint(__SAVb)
_FT_b = (ib_h * Lh + i * Bb + ib_f * Lf)*(1/(1 + g_ss)) - FT_b
qpprint(_FT_b)
_FD_b = eta_b * (1 - t_p)*FT_b - FD_b
qpprint(_FD_b)
_FU_b = FT_b - FD_b - Tb - FU_b
qpprint(_FU_b)
_Tb = t_p * FT_b - Tb
qpprint(_Tb)
_ib_f = i + spread - ib_f
qpprint(_ib_f)
_ib_h = i + spread - ib_h
qpprint(_ib_h)
_SFCb = _SAVb - __SAVb
pprint(_SFCb)
_Lh = levh*GDPNom - Lh
qpprint(_Lh)
_Mh = m_h * GDPNom - Mh
qpprint(_Mh)
_Bh = b_h * Bg - Bh
qpprint(_Bh)
# _b_h = (((__SAVh + SAVh)/D) - Mh + Lh) - b_h
# qpprint(_b_h)
_NWh = Mh + Bh - Lh - NWh
qpprint(_NWh)
_SAVh = D*NWh - SAVh
qpprint(_SAVh)
__SAVh = WB + (FD_f + FD_b) - (Pc*C) + (i*Bh)/(1 + g_ss) - (ib_h*Lh)/(1 + g_ss) - Th - SAVh
qpprint(__SAVh)
_C = (a_w * WB * (1 - t_w) + a_p * YDk_h + Z)/Pc - C
qpprint(_C)
_YG_h = YT_h + (i * Bh)/(1 + g_ss) - YG_h
qpprint(_YG_h)
_YT_h = WB + (FD_f + FD_b) - YT_h
qpprint(_YT_h)
_YD_h = (1 - t_w)* YG_h - (ib_h * Lh)/(1 + g_ss) - YD_h
qpprint(_YD_h)
_WB = Nf * wnom - WB
qpprint(_WB)
_YK_h = (FD_f + FD_b)/(1 + g_ss) + (i * Bh)/(1 + g_ss) - YK_h
qpprint(_YK_h)
_YDk_h = (1 - t_w)*YK_h - (ib_h * Lh)/(1 + g_ss) - YDk_h
qpprint(_YDk_h)
_Th = t_w*YG_h - Th
qpprint(_Th)
_Z = D*Lh + AM_h - Z
qpprint(_Z)
_AM_h = ((g_ss + ra_h)/(lambda_h * (1 + g_ss)))*(Lh * sumAm_h) - AM_h ## TODO Change to a more human readable equation: sp.Sum(NL/(((1 + g)**(j))*lada), (j, t - lada, t)).doit().simplify().collect(g + 1)?
qpprint(_AM_h)
_SFCh = _SAVh - __SAVh
pprint(_SFCh)
NFWf = Mf - Lf ## Aux. variable. Not to be solved
NLf = AM_f + D*Lf ## Aux. variable. Not to be solved
CF_f = Pc*Sc - WB ## Aux. variable. Not to be solved
_NLf = Pc*If + wnom*Nf - (Mf)/(1 + g_ss) - NLf
qpprint(_NLf)
_SAVf = FU_f - Pc*If - SAVf
qpprint(_SAVf)
__SAVf = D * (NFWf) - SAVf ## Only valid considering NFW (not NW)
qpprint(__SAVf)
_NWf = NFWf + r*Pc*K + ULC*INV - NWf
qpprint(_NWf)
_AM_f = ((g_ss + ra_f)/(lambda_f * (1 + g_ss)))*(Lh * sumAm_f) - AM_f ## TODO Change to a more human readable equation
qpprint(_AM_f)
_FT_f = CF_f - ib_f*(Lf/(1+g_ss)) - FT_f
qpprint(_FT_f)
_FN_f = FT_f + D*INV*ULC - DepNom - FN_f
qpprint(_FN_f)
_DepNom = (delta * Pc * K * r ) /(1 + g_ss) - DepNom
qpprint(_DepNom)
_FD_f = (1 - t_p)*(eta_f)*FN_f - FD_f
qpprint(_FD_f)
_FU_f = FT_f - FD_f - Tf - FU_f
qpprint(_FU_f)
_Tf = t_p*FN_f - Tf
qpprint(_Tf)
_Lf = levf*GDPNom - Lf
qpprint(_Lf)
_SFCf = _SAVf - __SAVf
pprint(_SFCf)
eqs = [
# _sumAm_h,
# _ra_h,
# _sumAm_f,
# _ra_f,
_GDPNom,
_Y_FC,
_un,
_RD,
_RDNom,
_Pc,
_ULC,
_INV,
__Yc,
_Sc,
_G,
_h,
_Nf,
_NWb,
_Mb,
_Lb,
_Bb,
_SAVb,
__SAVb,
_FT_b,
_FD_b,
_FU_b,
_Tb,
_ib_f,
_ib_h,
_Lh,
_Mh,
_Bh,
_NWh,
_SAVh,
__SAVh,
_C,
_YG_h,
_YT_h,
_YD_h,
_WB,
_YK_h,
_YDk_h,
_Th,
_Z,
_AM_h,
_NLf,
_SAVf,
__SAVf,
_NWf,
_AM_f,
_If,
_ht,
_FT_f,
_FN_f,
_DepNom,
_FD_f,
_FU_f,
_Tf,
### Externally imposed
_Bg,
_Lf,
_SFCb,
_SFCg,
_SFCf,
_SFCh,
NWb + NWf + NWg + NWh - r*Pc*K - ULC*INV,
]
outputs = [
a_p,
eta_f,
# wnom,
# A,
# b_h,
# Pc,
# t_w,
]
exogs = [
a_w,
delta,
eta_b,
B_Y,
g_z,
G_Y,
gamma_RD,
i,
iota,
K,
lambda_f,
lambda_h,
levf,
levh,
m_h,
Nf,
v,
pi,
r,
t_p,
theta,
un,
j,
spread,
g_ss,
D,
]
sol = sp.solve(
eqs,
outputs,
# exclude = exogs,
# force = True, # Make positive all symbols without assumptions regarding sign.
# manual = True
)
print(sol)
The importat part is: when I specify outputs as:
outputs = [
a_p,
eta_f,
# wnom,
# A,
# b_h,
# Pc,
# t_w,
]
It returns a non-empty solution. However, if I uncomment any other variables (e.g. wnom), it returns an EmptySet.
In this sense, it looks like that the system has some inconsistency. However, if this is the case, is it not weird that it worked with fewer variables to solve for? Is this a mathematical problem, or am I hitting some internal max iteration in sympy?
Any suggestion on how to address this issue?
I've tried a combination of variables in outputs object, but only if I uncoment the first two it gives me a non-empty result. I was expecting that if the model is inconsistent, it would always return an empty set.
I've also tried the focus function in this issue, but I've got nothing new.
I've also tried linsolve, but since my system is non-linear, I coudn't move foward with it.
