Linearsolve issue

Python Help
1

Hello,

I am new to Python - I am trying to do some stochastic simulations of the model but my code won’t run - it reports invertibility condition violation, but I cannot figure out why. The code works fine in Dynare but not sure why it won’t run in Python.

I would appreciate any suggestion.

Thank you in advance for your help!

```python
# Input model parameters using a pandas Series
parameters = pd.Series(dtype=float)

# Parameter values for the model
parameters['alpha'] =0.5
parameters['beta'] = 0.97
parameters['alpham'] = 0.43
parameters['vartheta'] = 0.5
parameters['xi_bar'] = 0.233
parameters['iota'] = 1
parameters['phi'] = 1
parameters['eps_p'] = 6
parameters['eta'] = 1
parameters['theta'] = 0.2
parameters['psi'] = 2
parameters['eta_z'] = 5.565
parameters['zeta_bar'] = 1
parameters['i_bar'] = 0.003
parameters['B'] = 0.7
parameters['b_bar_ss'] = 0.31
parameters['rho_r'] = 0.3389  
parameters['rho_ss'] = 0.967
parameters['sigma_r'] = 0.0028
parameters['rho_b'] = 0.9448
parameters['sigma_b'] = 0.0038
parameters['rho_z'] = 0.9752
parameters['sigma_z'] = 0.0031
parameters['rho_i'] = 0.8527
parameters['sigma_i'] = 0.0013
parameters['eps_ss'] = 0
parameters['sigma_mw'] = 0.0216
parameters['sigma_tau'] = 0.001
parameters['z_ss'] = 1.00  

# model equations
def equations(variables_forward, variables_current, parameters):

    # Extract the parameters
    p = parameters

    # Extract the variables for the current period and the forward period
    fwd = variables_forward
    cur = variables_current
    
    # Euler
    euler_eq= p.beta*((1+fwd.i)/fwd.pi)*((fwd.e/fwd.c_e) + (fwd.u/fwd.c_u)+((1-fwd.l)/fwd.c_n))-(1/cur.c_e)

    # BC and liquidity constraints
    bc_eq=fwd.b+(1+cur.i)*cur.a-(1+cur.i)*cur.b-cur.x
    consU_eq=cur.x+cur.tau_u-cur.c_u;
    consN_eq=cur.x-cur.c_n;
    asset_eq=cur.x+(1-cur.tau)*(cur.e*cur.w+cur.D)+cur.u*cur.tau_u -(cur.e*cur.c_e+cur.u*cur.c_u+(1-cur.l)*cur.c_n)-fwd.a

   # employment
    emp_eq=fwd.rho*cur.e+fwd.fs*fwd.sigma*fwd.s-fwd.e

   # searchers
    s_eq=fwd.l-cur.e-fwd.s

   # participation condition
    part_eq=cur.fs*cur.sigma*((1-cur.tau)*cur.w -cur.Gamma)+cur.fs*cur.sigma*p.beta*(cur.c_e/fwd.c_e)*((fwd.rho-fwd.fs*fwd.sigma)/fwd.fs*fwd.sigma)*(fwd.MRS_NC-fwd.Omega)-cur.MRS_NC+cur.Omega

   #where

    homeut_eq=p.xi_bar*cur.c_e/(1-cur.l)-cur.MRS_NC;

    unoppc_eq=cur.tau_u + (cur.c_n-cur.c_u)+cur.c_e*(np.log(cur.c_u) - np.log(cur.c_n)-(p.zeta_bar/(1+p.eta_z))*cur.sigma**(1+p.eta_z))-cur.Omega
    eoppc_eq=cur.tau_u + (cur.c_e-cur.c_u)+cur.c_e*(np.log(cur.c_u) - np.log(cur.c_e)+(cur.e**p.iota))-p.beta*cur.c_e*(p.zeta_bar/(1+p.eta_z))*fwd.sigma**(1+p.eta_z) - cur.Gamma

   # optimal search intensity

    searchint_eq=(1/cur.sigma)*(cur.MRS_NC-cur.Omega)-cur.c_e*p.zeta_bar*cur.sigma**p.eta_z
 
   # aggregate employment and unemployment 
    l_eq=cur.e+cur.u-cur.l;
    n_eq=1-cur.l-cur.n;

   # discount factor 
    discf_eq=p.beta*((1-fwd.tau)/(1-cur.tau))*(cur.c_e/fwd.c_e)-fwd.Lambda

   # hiring costs
    kappa_eq=cur.z*p.B*(cur.fs)**p.eta-cur.k

   #asset market equilibrium
    am_eq=cur.b_bar-cur.a
    bm_eq=cur.b_bar-cur.b

# Firms:

   #optimal hiring
    opthiring_eq=cur.q*cur.z-cur.w+ fwd.Lambda*fwd.rho*(fwd.k/fwd.fv)-(cur.k/cur.fv)

   # dividends definition

    dw_eq=cur.q*cur.z*cur.e-cur.w*cur.e-cur.k*cur.v-cur.dw

   # price
 
    optprice_eq=(cur.p_a/cur.p_b)-cur.p_starP
    pa_eq=(p.eps_p/(p.eps_p-1))*cur.q*cur.Y + fwd.Lambda*(1-p.theta)*fwd.p_a*(fwd.pi**p.eps_p)-cur.p_a
    pb_eq=cur.Y+fwd.Lambda*(1-p.theta)*fwd.p_b*(fwd.pi**(p.eps_p-1))-cur.p_b

   # inflation

    infl_eq=((1-p.theta)/(1-p.theta*((cur.p_starP)**(1-p.eps_p))))**(1/(1-p.eps_p))-cur.pi


   #Output
    y_eq=cur.z*cur.e/cur.varsigma-cur.Y

   # price dispersion
    varsigma_eq=p.theta*(fwd.p_starP)**(-p.eps_p)+(1-p.theta)*(fwd.pi**(p.eps_p))*cur.varsigma-fwd.varsigma

    # dividends
    D_eq=cur.Y-cur.q*cur.z*cur.e+cur.dw-cur.D

   # government
    tau_eq=cur.u*cur.tau_u-cur.tau*(cur.e*cur.w+cur.D)

   # unemployment benefits
    tauU_eq=0.3*cur.w-fwd.tau_u; 

   # Taylor rule 
    taylor_eq = (1+p.i_bar)*(cur.pi**(p.psi))*np.exp(cur.eps)-1-i

   # job finding rates
    fs_eq=p.alpham*(cur.v/cur.s*cur.sigma)**(1-p.alpha)-cur.fs
    fv_eq=p.alpham*(cur.v/cur.s*cur.sigma)**(-p.alpha)-cur.fv

   # bargained wage
    w_eq=p.vartheta*(cur.q*cur.z+fwd.Lambda*fwd.fs*fwd.sigma*(fwd.k/fwd.fv))+((1-p.vartheta)*cur.Gamma/(1-cur.tau))-cur.w
   
   # Agg consumption 
    cons_eq=cur.e*cur.c_e+cur.u*cur.c_u+(1-cur.l)*cur.c_n-cur.C
    
    # Exogenous technology process
    z_eq=(1-p.rho_z)*np.log(p.z_ss)+p.rho_z*np.log(cur.z)-fwd.z
    rho_eq=(1-p.rho_r)*np.log(p.rho_ss)+p.rho_r*np.log(cur.rho)-np.log(fwd.rho)
    bbar_eq=(1-p.rho_b)*p.b_bar_ss+p.rho_b*cur.b_bar-fwd.b_bar
    eps_eq =(1-p.rho_i)*p.eps_ss+p.rho_i*cur.eps-fwd.eps

    # Stack equilibrium conditions into a numpy array
    return np.array([
        euler_eq,
        bc_eq,
        consU_eq,
        consN_eq,
        asset_eq,
        emp_eq,
        s_eq,
        part_eq,
        homeut_eq,
        unoppc_eq,
        eoppc_eq,
        searchint_eq,
        l_eq,
        n_eq,
        discf_eq,
        kappa_eq,
        am_eq,
        bm_eq,
        opthiring_eq,
        dw_eq,
        optprice_eq,
        pa_eq,
        pb_eq,
        infl_eq,
        y_eq,
        varsigma_eq,
        D_eq,
        tau_eq,
        tauU_eq,
        taylor_eq,
        fs_eq,
        fv_eq,
        w_eq,
        cons_eq,
        z_eq,
        rho_eq,
        bbar_eq,
        eps_eq
    ])

SM = ls.model(equations=equations,
               n_states=7,
               n_exo_states=4,
               var_names=['c_e', 'c_u', 'l', 'u', 'e', 's', 'k', 'v', 'fs', 'fv', 'C', 'tau_u', 'w', 'n', 'sigma', 'MRS_NC', 'Omega', 'Gamma', 'tau', 'Y', 'D', 'dw', 'p_a', 'p_b', 'z', 'rho', 'b_bar', 'eps', 'varsigma',  'Lambda', 'i', 'pi',  'q', 'p_starP', 'x',  'c_n', 'a', 'b'],
               parameters=parameters)



#initial guess for the steady-state values of the model variables
#guess = [0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5 , 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5]

# Call the 'compute_ss' method of the  model and pass the 'guess' as input
#SM.compute_ss(guess)

# Print the computed steady-state values
#print(SM.ss)

#SM.linear_approximation()

# Solve the model
#SM.solve_klein(SM.a,SM.b)

# Print the matrices F and P 
#print('The matrix F:\n\n', np.around(SM.f, 2), '\n\n')
#print('The matrix P:\n\n', np.around(SM.p, 2))

# II way: solving manually for the steady state instead of using linearsolve

beta= 0.97
alpham=0.43
alpha=0.5
vartheta=0.5
xi_bar=0.233
iota=1
phi=1
eps_p=6
eta=1
theta=0.2
psi=2
eta_z=5.565
zeta_bar=1
i_bar=0.003
B=0.7
z_ss=1
b_bar_ss=0.31
rho_ss=0.967
rho_r=0.3389
sigma_r=0.0028
rho_b=0.9448
sigma_b=0.0038
rho_z=0.9752
sigma_z=0.0031
rho_i =0.8527
sigma_i =0.0013
sigma_mw =0.0216
sigma_tau=0.001

z=z_ss
rho=rho_ss
b_bar=b_bar_ss
eps=0
varsigma=1
i=i_bar
pi=1
p_starP=1
q=(eps_p-1)/eps_p
Lambda=beta
x=b_bar_ss
c_n=b_bar_ss
a=b_bar_ss
b=b_bar_ss

def system(variables):
    # Unpack the variables from the input tuple
    c_e, c_u, l, u, e, s, k, v, fs, fv, C, tau_u, w, n, sigma, MRS_NC, Omega, Gamma, tau, Y, D, dw, p_a, p_b = variables

    # Equilibrium conditions
    # Euler equation
    x1 =  beta*((1+i)/pi)*((e/c_e) + (u/c_u)+((1-l)/c_n))-(1/c_e)

    # consumption of the unemployed
    x2 =x+ tau_u-c_u
    
    # end of period asset constraint 
    x3 = (1-tau)*e*w+(1-tau)*D+u*tau_u -(e*c_e+u*c_u+(1-l)*c_n)

    # employment
    x4 = rho*e+fs*sigma*s-e
    
    # searchers 
    x5=l-e-s
    
    # participation condition
    x6= MRS_NC-Omega - fs*sigma*((1-tau)*w -Gamma)-beta*(rho-fs*sigma)*(MRS_NC-Omega)
    
    # marginal utility from home production 
    x7 = MRS_NC-xi_bar*c_e/(1-l)
    
    # opportunity cost of employment
    x8 = tau_u + (c_e-c_u)+c_e*(np.log(c_u) - np.log(c_e)+(e**iota))-beta*c_e*(zeta_bar/(1+eta_z))*sigma**(1+eta_z)-Gamma

    # opportunity cost of nonparticipation
    x9 = tau_u + (c_n-c_u)+c_e*(np.log(c_u) - np.log(c_n)-(zeta_bar/(1+eta_z))*sigma**(1+eta_z)) - Omega
    
    # search intensity
    x10 = c_e*zeta_bar*sigma**eta_z-(1/sigma)*(MRS_NC-Omega)
    
    #employment and unemployment
    x11=e+u-l
    x12=1-n-l
    
    # hiring costs
    x13=z*B*(fs)**eta-k

    # optimal hiring
    x14 =q*z-w+beta*rho*(k/fv)-(k/fv)
    
    # dividends
    x15 = q*z*e-w*e-k*v-dw
    
    # prices
    x16=(eps_p/(eps_p-1))*q*Y + Lambda*(1-theta)*p_a*(pi**eps_p)-p_a
    x17=Y+Lambda*(1-theta)*p_b*(pi**(eps_p-1))-p_b
    
    # output
    x18 = Y-z*e
    
    # total dividends
    x19=Y-q*z*e+dw-D
    
    # total 
    x20=tau*(e*w+D)-u*tau_u
    
    # unemployment benefits 
    x21=tau_u-0.3*w
    
    # job finding/filling rates 
    x22=fs-alpham*(v/s*sigma)**(1-alpha)
    x23=fv-alpham*(v/s*sigma)**(-alpha)
    
    # wage 
    x24=w- vartheta*(q*z+Lambda*fs*sigma*(k/fv))-(1-vartheta)*Gamma/(1-tau)

        
    # Return the system of equilibrium conditions as a tuple
    return x1, x2, x3, x4, x5, x6, x7, x8, x9, x10, x11, x12, x13, x14, x15, x16, x17, x18, x19, x20, x21, x22, x23, x24


guess =[0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5]

# Find the root of the system of equations using scipy.optimize.root
sol = scipy.optimize.root(system, guess).x

c_e, c_u, l, u, e, s, k, v, fs, fv, C, tau_u, w, n, sigma, MRS_NC, Omega, Gamma, tau, Y, D, dw, p_a, p_b  = sol

sol


# Set the steady-state values for the model 'SM' using the 'set_ss' method
# The input list contains the steady-state values for the model variables in the following order:
# [c_e, c_u, l, u, e, s, k, v, fs, fv, C, tau_u, w, n, sigma, MRS_NC, Omega, Gamma, tau, Y, D, dw, p_a, p_b, z, rho, b_bar, eps, varsigma, Lambda, i, pi, q, p_starP, x, c_n, a, b]

SM.set_ss([3.69532502e-01, 5.45157372e-01, 6.52403573e-01, 1.13545498e-01, 5.37484294e-01, 1.09805537e-01, 3.67914880e-01, 4.73705833e-01, 5.21852766e-01, 3.49653704e-01, 5.42416027e+24, 2.46132090e-01, 7.57740932e-01, 3.47641784e-01, 3.42227646e-01, 2.35790147e-01, 2.32562219e-01, 4.26236441e-01, 1.14461118e-01, 5.40639579e-01, -3.19475457e-02, -1.28813194e-01,  2.41260087e+00, 2.41260087e+00, 1, 0.967, 0.31, 0, 1,0.97, 0.003,1, 0.8333, 1, 0.31, 0.31, 0.31, 0.31])


SM.linear_approximation()

# Print the matrices A and B after rounding the values to four decimal places
print('The matrix A:\n\n',np.around(SM.a,4),'\n\n')
print('The matrix B:\n\n',np.around(SM.b,4),'\n\n')

# Solve the model
SM.solve_klein(SM.a,SM.b)
2

Please include the exact error message, including traceback information.

The code you’ve posted is incomplete; it’s missing at least your imports, possibly more. This makes it difficult to tell what might be wrong with it.

3

Thank you, Alexander, for your reply.

It reports the following errors:

```python
SM.linear_approximation()
<ipython-input-11-dfc4f4c4360e>:105: RuntimeWarning: invalid value encountered in double_scalars
  infl_eq=((1-p.theta)/(1-p.theta*((cur.p_starP)**(1-p.eps_p))))**(1/(1-p.eps_p))-cur.pi

There seem to be some NaN values in the SM.a and SM.b matrices, so that when I run:

```python
SM.solve_klein(SM.a,SM.b)

it reports: An exception has occurred, use %tb to see the full traceback.
SystemExit: Invertibility condition violated.

I am not sure what you mean with the imports, below is the code again, hope it runs now showing the above errors. Thanks again!

```python
parameters = pd.Series(dtype=float)

# Parameter values for the model
parameters['alpha'] =0.5
parameters['beta'] = 0.97
parameters['alpham'] = 0.43
parameters['vartheta'] = 0.5
parameters['xi_bar'] = 0.233
parameters['iota'] = 1
parameters['phi'] = 1
parameters['eps_p'] = 6
parameters['eta'] = 1
parameters['theta'] = 0.2
parameters['psi'] = 2
parameters['eta_z'] = 5.565
parameters['zeta_bar'] = 1
parameters['i_bar'] = 0.003
parameters['B'] = 0.7
parameters['b_bar_ss'] = 0.31
parameters['rho_r'] = 0.3389  
parameters['rho_ss'] = 0.967
parameters['sigma_r'] = 0.0028
parameters['rho_b'] = 0.9448
parameters['sigma_b'] = 0.0038
parameters['rho_z'] = 0.9752
parameters['sigma_z'] = 0.0031
parameters['rho_i'] = 0.8527
parameters['sigma_i'] = 0.0013
parameters['eps_ss'] = 0
parameters['sigma_mw'] = 0.0216
parameters['sigma_tau'] = 0.001
parameters['z_ss'] = 1.00  

# model equations
def equations(variables_forward, variables_current, parameters):

    # Extract the parameters
    p = parameters

    # Extract the variables for the current period and the forward period
    fwd = variables_forward
    cur = variables_current
    
    # Euler
    euler_eq= p.beta*((1+fwd.i)/fwd.pi)*((fwd.e/fwd.c_e) + (fwd.u/fwd.c_u)+((1-fwd.l)/fwd.c_n))-(1/cur.c_e)

    # BC and liquidity constraints
    bc_eq=fwd.b+(1+cur.i)*cur.a-(1+cur.i)*cur.b-cur.x
    consU_eq=cur.x+cur.tau_u-cur.c_u;
    consN_eq=cur.x-cur.c_n;
    asset_eq=cur.x+(1-cur.tau)*(cur.e*cur.w+cur.D)+cur.u*cur.tau_u -(cur.e*cur.c_e+cur.u*cur.c_u+(1-cur.l)*cur.c_n)-fwd.a

   # employment
    emp_eq=fwd.rho*cur.e+fwd.fs*fwd.sigma*fwd.s-fwd.e

   # searchers
    s_eq=fwd.l-cur.e-fwd.s

   # participation condition
    part_eq=cur.fs*cur.sigma*((1-cur.tau)*cur.w -cur.Gamma)+cur.fs*cur.sigma*p.beta*(cur.c_e/fwd.c_e)*((fwd.rho-fwd.fs*fwd.sigma)/fwd.fs*fwd.sigma)*(fwd.MRS_NC-fwd.Omega)-cur.MRS_NC+cur.Omega

   #where

    homeut_eq=p.xi_bar*cur.c_e/(1-cur.l)-cur.MRS_NC;

    unoppc_eq=cur.tau_u + (cur.c_n-cur.c_u)+cur.c_e*(np.log(cur.c_u) - np.log(cur.c_n)-(p.zeta_bar/(1+p.eta_z))*cur.sigma**(1+p.eta_z))-cur.Omega
    eoppc_eq=cur.tau_u + (cur.c_e-cur.c_u)+cur.c_e*(np.log(cur.c_u) - np.log(cur.c_e)+(cur.e**p.iota))-p.beta*cur.c_e*(p.zeta_bar/(1+p.eta_z))*fwd.sigma**(1+p.eta_z) - cur.Gamma

   # optimal search intensity

    searchint_eq=(1/cur.sigma)*(cur.MRS_NC-cur.Omega)-cur.c_e*p.zeta_bar*cur.sigma**p.eta_z
 
   # aggregate employment and unemployment 
    l_eq=cur.e+cur.u-cur.l;
    n_eq=1-cur.l-cur.n;

   # discount factor 
    discf_eq=p.beta*((1-fwd.tau)/(1-cur.tau))*(cur.c_e/fwd.c_e)-fwd.Lambda

   # hiring costs
    kappa_eq=cur.z*p.B*(cur.fs)**p.eta-cur.k

   #asset market equilibrium
    am_eq=cur.b_bar-cur.a
    bm_eq=cur.b_bar-cur.b

# Firms:

   #optimal hiring
    opthiring_eq=cur.q*cur.z-cur.w+ fwd.Lambda*fwd.rho*(fwd.k/fwd.fv)-(cur.k/cur.fv)

   # dividends definition

    dw_eq=cur.q*cur.z*cur.e-cur.w*cur.e-cur.k*cur.v-cur.dw

   # price
 
    optprice_eq=(cur.p_a/cur.p_b)-cur.p_starP
    pa_eq=(p.eps_p/(p.eps_p-1))*cur.q*cur.Y + fwd.Lambda*(1-p.theta)*fwd.p_a*(fwd.pi**p.eps_p)-cur.p_a
    pb_eq=cur.Y+fwd.Lambda*(1-p.theta)*fwd.p_b*(fwd.pi**(p.eps_p-1))-cur.p_b

   # inflation

    infl_eq=((1-p.theta)/(1-p.theta*((cur.p_starP)**(1-p.eps_p))))**(1/(1-p.eps_p))-cur.pi


   #Output
    y_eq=cur.z*cur.e/cur.varsigma-cur.Y

   # price dispersion
    varsigma_eq=p.theta*(fwd.p_starP)**(-p.eps_p)+(1-p.theta)*(fwd.pi**(p.eps_p))*cur.varsigma-fwd.varsigma

    # dividends
    D_eq=cur.Y-cur.q*cur.z*cur.e+cur.dw-cur.D

   # government
    tau_eq=cur.u*cur.tau_u-cur.tau*(cur.e*cur.w+cur.D)

   # unemployment benefits
    tauU_eq=0.3*cur.w-fwd.tau_u; 

   # Taylor rule 
    taylor_eq = (1+p.i_bar)*(cur.pi**(p.psi))*np.exp(cur.eps)-1-i

   # job finding rates
    fs_eq=p.alpham*(cur.v/cur.s*cur.sigma)**(1-p.alpha)-cur.fs
    fv_eq=p.alpham*(cur.v/cur.s*cur.sigma)**(-p.alpha)-cur.fv

   # bargained wage
    w_eq=p.vartheta*(cur.q*cur.z+fwd.Lambda*fwd.fs*fwd.sigma*(fwd.k/fwd.fv))+((1-p.vartheta)*cur.Gamma/(1-cur.tau))-cur.w
   
   # Agg consumption 
    cons_eq=cur.e*cur.c_e+cur.u*cur.c_u+(1-cur.l)*cur.c_n-cur.C
    
    # Exogenous technology process
    z_eq=(1-p.rho_z)*np.log(p.z_ss)+p.rho_z*np.log(cur.z)-fwd.z
    rho_eq=(1-p.rho_r)*np.log(p.rho_ss)+p.rho_r*np.log(cur.rho)-np.log(fwd.rho)
    bbar_eq=(1-p.rho_b)*p.b_bar_ss+p.rho_b*cur.b_bar-fwd.b_bar
    eps_eq =(1-p.rho_i)*p.eps_ss+p.rho_i*cur.eps-fwd.eps

    # Stack equilibrium conditions into a numpy array
    return np.array([
        euler_eq,
        bc_eq,
        consU_eq,
        consN_eq,
        asset_eq,
        emp_eq,
        s_eq,
        part_eq,
        homeut_eq,
        unoppc_eq,
        eoppc_eq,
        searchint_eq,
        l_eq,
        n_eq,
        discf_eq,
        kappa_eq,
        am_eq,
        bm_eq,
        opthiring_eq,
        dw_eq,
        optprice_eq,
        pa_eq,
        pb_eq,
        infl_eq,
        y_eq,
        varsigma_eq,
        D_eq,
        tau_eq,
        tauU_eq,
        taylor_eq,
        fs_eq,
        fv_eq,
        w_eq,
        cons_eq,
        z_eq,
        rho_eq,
        bbar_eq,
        eps_eq
    ])

SM = ls.model(equations=equations,
               n_states=7,
               n_exo_states=4,
               var_names=['c_e', 'c_u', 'l', 'u', 'e', 's', 'k', 'v', 'fs', 'fv', 'C', 'tau_u', 'w', 'n', 'sigma', 'MRS_NC', 'Omega', 'Gamma', 'tau', 'Y', 'D', 'dw', 'p_a', 'p_b', 'z', 'rho', 'b_bar', 'eps', 'varsigma',  'Lambda', 'i', 'pi',  'q', 'p_starP', 'x',  'c_n', 'a', 'b'],
               parameters=parameters)


# Compute the steady-state values of the model using the 'linearsolve' package instead

# 'guess' is a list containing the initial guess for the steady-state values of the model variables
guess = [0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5 , 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5]

# Call the 'compute_ss' method of the model and pass the 'guess' as input
SM.compute_ss(guess)

# Print the computed steady-state values
print(SM.ss)

SM.linear_approximation()

print('The matrix A:\n\n',np.around(SM.a,4),'\n\n')
print('The matrix B:\n\n',np.around(SM.b,4),'\n\n')

# Solve the nk model
SM.solve_klein(SM.a,SM.b)

# Print the matrices F and P 
print('The matrix F:\n\n', np.around(SM.f, 2), '\n\n')
print('The matrix P:\n\n', np.around(SM.p, 2))
4

That’s not an error, it’s a warning. It may or may not be related to the SystemExit error you see later.

Do what it says: “use %tb to see the full traceback”, and post it here.

It does not, it errors out on line 1 with

NameError: name 'pd' is not defined

Somewhere above this line you have done import pandas as pd. You also have import numpy as np, and probably import linearsolve as ls. These are “imports”; they make it so that you can use code from libraries in your own code.

You haven’t included these in the code you posted. I can guess what some of them are, but please just include them in your post instead. Since you omitted part of your code, I can’t know if you also omitted other relevant code.

5

Also, on this line:

the variable i is undefined.

6

Here is the full traceback:

```python
---------------------------------------------------------------------------
SystemExit                                Traceback (most recent call last)
<ipython-input-17-23186cfbe049> in <module>
    213 
    214 # Solve the model
--> 215 SM.solve_klein(SM.a,SM.b)
    216 
    217 # Print the matrices F and P

~\Anaconda3\lib\site-packages\linearsolve\__init__.py in solve_klein(self, a, b, eigenvalue_warnings)
    605             b = self.b
    606 
--> 607         self.f,n,self.p,l,self.stab,self.eig = klein(a=a,b=b,c=None,phi=None,n_states=self.n_states,eigenvalue_warnings=eigenvalue_warnings)
    608 
    609         if not np.iscomplexobj(self.parameters):

~\Anaconda3\lib\site-packages\linearsolve\__init__.py in klein(a, b, c, phi, n_states, eigenvalue_warnings)
    924     if n_states>0:
    925         if np.linalg.matrix_rank(z11)<n_states:
--> 926             sys.exit("Invertibility condition violated.")
    927 
    928     s11 = s[0:n_states,0:n_states];

SystemExit: Invertibility condition violated.

Oh, sorry about the imports. Thanks a lot! Here is the code again, corrected for i as well:

```python

#pip install linearsolve
import numpy as np  
np.set_printoptions(threshold=np.inf)
import pandas as pd  
import linearsolve as ls  
import matplotlib.pyplot as plt  
import scipy.optimize  

parameters = pd.Series(dtype=float)

# Parameter values for the model
parameters['alpha'] =0.5
parameters['beta'] = 0.97
parameters['alpham'] = 0.43
parameters['vartheta'] = 0.5
parameters['xi_bar'] = 0.233
parameters['iota'] = 1
parameters['phi'] = 1
parameters['eps_p'] = 6
parameters['eta'] = 1
parameters['theta'] = 0.2
parameters['psi'] = 2
parameters['eta_z'] = 5.565
parameters['zeta_bar'] = 1
parameters['i_bar'] = 0.003
parameters['B'] = 0.7
parameters['b_bar_ss'] = 0.31
parameters['rho_r'] = 0.3389  
parameters['rho_ss'] = 0.967
parameters['sigma_r'] = 0.0028
parameters['rho_b'] = 0.9448
parameters['sigma_b'] = 0.0038
parameters['rho_z'] = 0.9752
parameters['sigma_z'] = 0.0031
parameters['rho_i'] = 0.8527
parameters['sigma_i'] = 0.0013
parameters['eps_ss'] = 0
parameters['sigma_mw'] = 0.0216
parameters['sigma_tau'] = 0.001
parameters['z_ss'] = 1.00  

# model equations
def equations(variables_forward, variables_current, parameters):

    # Extract the parameters
    p = parameters

    # Extract the variables for the current period and the forward period
    fwd = variables_forward
    cur = variables_current
    
    # Euler
    euler_eq= p.beta*((1+fwd.i)/fwd.pi)*((fwd.e/fwd.c_e) + (fwd.u/fwd.c_u)+((1-fwd.l)/fwd.c_n))-(1/cur.c_e)

    # BC and liquidity constraints
    bc_eq=fwd.b+(1+cur.i)*cur.a-(1+cur.i)*cur.b-cur.x
    consU_eq=cur.x+cur.tau_u-cur.c_u;
    consN_eq=cur.x-cur.c_n;
    asset_eq=cur.x+(1-cur.tau)*(cur.e*cur.w+cur.D)+cur.u*cur.tau_u -(cur.e*cur.c_e+cur.u*cur.c_u+(1-cur.l)*cur.c_n)-fwd.a

   # employment
    emp_eq=fwd.rho*cur.e+fwd.fs*fwd.sigma*fwd.s-fwd.e

   # searchers
    s_eq=fwd.l-cur.e-fwd.s

   # participation condition
    part_eq=cur.fs*cur.sigma*((1-cur.tau)*cur.w -cur.Gamma)+cur.fs*cur.sigma*p.beta*(cur.c_e/fwd.c_e)*((fwd.rho-fwd.fs*fwd.sigma)/fwd.fs*fwd.sigma)*(fwd.MRS_NC-fwd.Omega)-cur.MRS_NC+cur.Omega

   #where

    homeut_eq=p.xi_bar*cur.c_e/(1-cur.l)-cur.MRS_NC;

    unoppc_eq=cur.tau_u + (cur.c_n-cur.c_u)+cur.c_e*(np.log(cur.c_u) - np.log(cur.c_n)-(p.zeta_bar/(1+p.eta_z))*cur.sigma**(1+p.eta_z))-cur.Omega
    eoppc_eq=cur.tau_u + (cur.c_e-cur.c_u)+cur.c_e*(np.log(cur.c_u) - np.log(cur.c_e)+(cur.e**p.iota))-p.beta*cur.c_e*(p.zeta_bar/(1+p.eta_z))*fwd.sigma**(1+p.eta_z) - cur.Gamma

   # optimal search intensity

    searchint_eq=(1/cur.sigma)*(cur.MRS_NC-cur.Omega)-cur.c_e*p.zeta_bar*cur.sigma**p.eta_z
 
   # aggregate employment and unemployment 
    l_eq=cur.e+cur.u-cur.l;
    n_eq=1-cur.l-cur.n;

   # discount factor 
    discf_eq=p.beta*((1-fwd.tau)/(1-cur.tau))*(cur.c_e/fwd.c_e)-fwd.Lambda

   # hiring costs
    kappa_eq=cur.z*p.B*(cur.fs)**p.eta-cur.k

   #asset market equilibrium
    am_eq=cur.b_bar-cur.a
    bm_eq=cur.b_bar-cur.b

# Firms:

   #optimal hiring
    opthiring_eq=cur.q*cur.z-cur.w+ fwd.Lambda*fwd.rho*(fwd.k/fwd.fv)-(cur.k/cur.fv)

   # dividends definition

    dw_eq=cur.q*cur.z*cur.e-cur.w*cur.e-cur.k*cur.v-cur.dw

   # price
 
    optprice_eq=(cur.p_a/cur.p_b)-cur.p_starP
    pa_eq=(p.eps_p/(p.eps_p-1))*cur.q*cur.Y + fwd.Lambda*(1-p.theta)*fwd.p_a*(fwd.pi**p.eps_p)-cur.p_a
    pb_eq=cur.Y+fwd.Lambda*(1-p.theta)*fwd.p_b*(fwd.pi**(p.eps_p-1))-cur.p_b

   # inflation

    infl_eq=((1-p.theta)/(1-p.theta*((cur.p_starP)**(1-p.eps_p))))**(1/(1-p.eps_p))-cur.pi


   #Output
    y_eq=cur.z*cur.e/cur.varsigma-cur.Y

   # price dispersion
    varsigma_eq=p.theta*(fwd.p_starP)**(-p.eps_p)+(1-p.theta)*(fwd.pi**(p.eps_p))*cur.varsigma-fwd.varsigma

    # dividends
    D_eq=cur.Y-cur.q*cur.z*cur.e+cur.dw-cur.D

   # government
    tau_eq=cur.u*cur.tau_u-cur.tau*(cur.e*cur.w+cur.D)

   # unemployment benefits
    tauU_eq=0.3*cur.w-fwd.tau_u; 

   # Taylor rule 
    taylor_eq = (1+p.i_bar)*(cur.pi**(p.psi))*np.exp(cur.eps)-1-cur.i

   # job finding rates
    fs_eq=p.alpham*(cur.v/cur.s*cur.sigma)**(1-p.alpha)-cur.fs
    fv_eq=p.alpham*(cur.v/cur.s*cur.sigma)**(-p.alpha)-cur.fv

   # bargained wage
    w_eq=p.vartheta*(cur.q*cur.z+fwd.Lambda*fwd.fs*fwd.sigma*(fwd.k/fwd.fv))+((1-p.vartheta)*cur.Gamma/(1-cur.tau))-cur.w
   
   # Agg consumption 
    cons_eq=cur.e*cur.c_e+cur.u*cur.c_u+(1-cur.l)*cur.c_n-cur.C
    
    # Exogenous technology process
    z_eq=(1-p.rho_z)*np.log(p.z_ss)+p.rho_z*np.log(cur.z)-fwd.z
    rho_eq=(1-p.rho_r)*np.log(p.rho_ss)+p.rho_r*np.log(cur.rho)-np.log(fwd.rho)
    bbar_eq=(1-p.rho_b)*p.b_bar_ss+p.rho_b*cur.b_bar-fwd.b_bar
    eps_eq =(1-p.rho_i)*p.eps_ss+p.rho_i*cur.eps-fwd.eps

    # Stack equilibrium conditions into a numpy array
    return np.array([
        euler_eq,
        bc_eq,
        consU_eq,
        consN_eq,
        asset_eq,
        emp_eq,
        s_eq,
        part_eq,
        homeut_eq,
        unoppc_eq,
        eoppc_eq,
        searchint_eq,
        l_eq,
        n_eq,
        discf_eq,
        kappa_eq,
        am_eq,
        bm_eq,
        opthiring_eq,
        dw_eq,
        optprice_eq,
        pa_eq,
        pb_eq,
        infl_eq,
        y_eq,
        varsigma_eq,
        D_eq,
        tau_eq,
        tauU_eq,
        taylor_eq,
        fs_eq,
        fv_eq,
        w_eq,
        cons_eq,
        z_eq,
        rho_eq,
        bbar_eq,
        eps_eq
    ])

SM = ls.model(equations=equations,
               n_states=7,
               n_exo_states=4,
               var_names=['c_e', 'c_u', 'l', 'u', 'e', 's', 'k', 'v', 'fs', 'fv', 'C', 'tau_u', 'w', 'n', 'sigma', 'MRS_NC', 'Omega', 'Gamma', 'tau', 'Y', 'D', 'dw', 'p_a', 'p_b', 'z', 'rho', 'b_bar', 'eps', 'varsigma',  'Lambda', 'i', 'pi',  'q', 'p_starP', 'x',  'c_n', 'a', 'b'],
               parameters=parameters)


# Compute the steady-state values

# 'guess' is a list containing the initial guess for the steady-state values of the model variables
guess = [0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5 , 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5]

# Call the 'compute_ss' method of the model and pass the 'guess' as input
SM.compute_ss(guess)

# Print the computed steady-state values
print(SM.ss)

SM.linear_approximation()

print('The matrix A:\n\n',np.around(SM.a,4),'\n\n')
print('The matrix B:\n\n',np.around(SM.b,4),'\n\n')

# Solve the model
SM.solve_klein(SM.a,SM.b)

# Print the matrices F and P 
print('The matrix F:\n\n', np.around(SM.f, 2), '\n\n')
print('The matrix P:\n\n', np.around(SM.p, 2))
7

linearsolve does not like your initial guess. The RuntimeWarnings are raised first by SM.compute_ss(guess), then again by SM.linear_approximation().

Later, SystemExit is raised in SM.solve_klein(SM.a, SM.b) because the QZ decomposition of your matrix pair A & B has insufficient rank. I can’t tell whether this is a bug in linearsolve or a problem with your code, but perhaps it may give you some idea how to proceed.

8

Thank you very much for your help. I tried using the same steady state values for the initial guess as the ones for which the code runs in Dynare but it did not help. Would you perhaps have any other suggestion what else I could try?

9

Try using the solution from Dynare as the initial guess for linearsolve. If linearsolve is able to solve the problem then, it is simply more sensitive to a poor initial guess than Dynare. You will need to somehow refine your initial guess before passing it to linearsolve.

If linearsolve still cannot solve the problem even with an initial guess very close to the solution, you may have found a bug in linearsolve. In that case, your options are to

  1. File a bug report on linearsolve’s issue tracker
  2. Fix the bug yourself
  3. Pay someone to fix it for you
  4. Use some other software
10

I see. Thanks a lot once again for all the help!