Replication Issues -- SS and Initval

[fabioit]

Hello,

I am a Master student with basic knowledge of dynare, i.e. I can run simple NK models and variants.
'Nuff said of me. I am trying to vaguely replicate "Fiscal Consolidation with Tax Evasion and Corruption" by Pappa et al (2015) with no success thus far.
My issues are:

1- Some variables are grouped together with the parameters in the "calibration"" section and are associated with some values: this puzzles me in general.

1.1- I have turned some of these variables into parameters to get #variables=#equations but I always get the following error when compiling:

Code: Select all

STEADY:  The Jacobian contains Inf or NaN. The problem arises from:

STEADY:  Derivative of Equation 12 with respect to Variable cc  (initial value of cc: 0)
STEADY:  Derivative of Equation 12 with respect to Variable c  (initial value of c: 0)

STEADY:  The problem most often occurs, because a variable with
STEADY:  exponent smaller than 1 has been initialized to 0. Taking the derivative
STEADY:  and evaluating it at the steady state then results in a division by 0.
Error using dynare_solve (line 60)
An element of the Jacobian is not finite or NaN
Error in evaluate_steady_state (line 66)
            [ys,check] = dynare_solve([M.fname '_static'],...
Error in resol (line 104)
[dr.ys,M.params,info] = evaluate_steady_state(oo.steady_state,M,options,oo,0);
Error in check (line 73)
[dr,info,M,options,oo] = resol(1,M,options,oo);
Error in rep_pappa4 (line 387)
oo_.dr.eigval = check(M_,options_,oo_);
Error in dynare (line 180)
evalin('base',fname) ;

The point is I have changed the values of cc and c to >0 with no avail. The same message pops up.

2- Steady state equations appear in the Appendix, as the standard r = 1/BETA, where and how do I specify them without turning R into a parameter? I simple models this isn't necessary as dynare computes them...

Hoping that someone will help a student in trouble,

Best,

Fabio

P.S.: my code at the moment is (don't laugh):

Code: Select all

var
//%controls
cc c nf ni u uf ui upsf upsi l k i pi r rr tr bigt lambdac lambdanf
lambdani wi wf px xi xf y mf mi surp tr
b
//%instruments
g taus;

varexo
eps_c eps_s eps_n eps_g;

parameters
ALPHA1 ALPHA2 ALPHAF ALPHAI BETA GAMMA DELTA ZETAPI ETA KAPPAF KAPPAI RHO SIGMAF
SIGMAI THETAF THETAI BIGPHI SMALLPHI OMEGA DIGAMMA XITR CHI PT ZSS MU1F MU1I
MU2 AF AI RRSS S psihf psihi psiff psifi tauc taun;

ALPHA1 = 0.85;
ALPHA2 = -0.25;
ALPHAF = 0.36;
ALPHAI = 0.4;
BETA = 0.96;
GAMMA = 0.3;
DELTA = 0.088;
ZETAPI = 1.5;
ETA = 2;
KAPPAF = 0.14;
KAPPAI = 0.13;
RHO = 0.02;
SIGMAF = 0.07;
SIGMAI = 0545;
THETAF = 0.22;
THETAI = 0.80;
BIGPHI = 0.7;
SMALLPHI = 2;
OMEGA = 0.5;
DIGAMMA = 0.2;      //% (VALUE) NOT FROM PAPPA
XITR = 0.2;
CHI = 0.25;
//%ZSS = 1.2;             //% NOT FROM PAPPA
PT = 0.8;      //% NOT FROM PAPPA
MU1F = 0.85;
MU1I = 0.12;
MU2 = 0.7;
AF = 1;
AI = 0.6;
RRSS = 0.3;      //%(VALUE) NOT FROM PAPPA
S = 0.10;
psihf = 0.63;
psihi = 0.91;
psiff = 0.96;
psifi = 0.05;
tauc = 0.16;
taun = 0.4;
//%------------------------------------------------------------
//% Model equations - Fabio
//%------------------------------------------------------------

model ;

//# S = s
//# PSIHF = psihf
//# PSIHI = psihi
//# PSIFF = psiff
//# PSIFI = psifi

//% ***LABOUR MARKETS***

//% Probability of filling formal vacancy
psiff = exp(mf)/exp(upsf) ;      //%d.f

//% Formal hiring probability
psihf = exp(mf)/exp(uf) ;      //%c.f

//% Probability of filling informal vacancy
psifi = exp(mi)/exp(upsi) ;      //%d.i

//% Informal hiring probability
psihi = exp(mi)/exp(ui) ;       //%c.i   

//% Law of motion of formal employment
exp(nf(+1)) = (1-SIGMAF)* exp(nf) + exp(mf)   ;   //% 2

//% Law of motion of informal employment
exp(ni(+1)) = (1-RHO-SIGMAI)*exp(ni) + exp(mi)   ;   //% 3

//% Formal matches
exp(mf) = MU1F*(exp(upsf)^MU2)*(exp(uf)^(1-MU2)) ;

//% Informal matches
exp(mi) = MU1I*(exp(upsi)^MU2)*(exp(ui)^(1-MU2)) ;

//% Total unmeployment
exp(u) = exp(uf) + exp(ui);

//%exp(ui) = exp(s)*exp(u)   ;   //%e
//%exp(uf) = (1-exp(s))*exp(u)   ;   //%f

//%==============================================================================

//% ***HOUSEHOLDS***

//% Final consumption Bundle
exp(cc) = ( ALPHA1*(exp(c))^(ALPHA2) + (1-ALPHA1)*(exp(g))^(ALPHA2) )^(1/ALPHA2);      //% b

//% Household composition
//% exp(nf) + exp(ni) + exp(u) + exp(l) = 1;      //% 4

//% Law of motion of (formal) capital
exp(k(+1)) = exp(i) + (1-DELTA)*exp(k) - (OMEGA/2)*( (exp(k(+1))/exp(k)-1)^2 ) * exp(k);      //% 5

//% [ct]
( (exp(cc))^(1-ETA-ALPHA2) ) * ALPHA1 * exp(((c))^(ALPHA2-1)) - exp(lambdac)*(1+tauc) = 0 ;
 
//% [k(+1)]
exp(lambdac) * (1 + OMEGA * ( exp(k(+1))/exp(k) -1) ) =
BETA * exp(lambdac(+1)) * ( (  1 - DELTA - ( (OMEGA/2*(exp(k(+2))/exp(k(+1)) - 1)^2) +
OMEGA*exp(k(+2))/exp(k(+1))*(exp(k(+2))/exp(k(+1)) - 1))) + r(+1) ) ;

//% [b(+1)]
exp(lambdac) * 1/exp(rr) = - BETA*exp(lambdac(+1)) * ( 1 / exp( pi(+1) ) ) ;

//% [nf(+1)]
exp(lambdanf) = - BETA*( BIGPHI*(exp(l(+1))^(-SMALLPHI)) -
exp(lambdac(+1)) * (1 - taun )*exp(wf(+1))*exp(lambdac(+1)) -
exp(lambdanf(+1))*(1-SIGMAF) ) ;

//% [ni(+1)]
exp(lambdani) = - BETA * ( BIGPHI*(exp(l(+1))^(-SMALLPHI)) -
exp(lambdac(+1)) * exp(wi(+1)) - exp(lambdani(+1))*(1-RHO-SIGMAI)) ;

//% [u]
//%BIGPHI*(exp(l(+1))^(-SMALLPHI)) = - lambdac*DIGAMMA*(1 - exp(s)) -
//%exp(lambdanf)*psihf*(1-exp(s)) - exp(lambdani)*psihi*exp(s) ;

//% [s]
//%exp(lambdani)*psihi = exp(lambdanf)*psihf + exp(lambdac)*DIGAMMA ;

//% [uf] [using uf==(1-s)*u] Household FOC wrt formal jobseekers
exp(lambdanf)*psihf = -DIGAMMA*exp(lambdac)+BIGPHI*(exp(l(+1))^(-SMALLPHI));
 
//% [ui] [using ui==s*u] Household FOC wrt informal jobseekers
exp(lambdani)*psihi = BIGPHI*(exp(l(+1))^(-SMALLPHI));

//%==============================================================================

//% ***PRODUCTION***

//% Production function in the formal sector
exp(xf) = (((AF*exp(nf))^(1-ALPHAF))*((exp(k))^ALPHAF));
//%xf = (1-ALPHAF)*(AF + nf) + ALPHAF*k;

//% Informal production function
exp(xi) = (AI*exp(ni))^(1-ALPHAI);
//%xi = (1-ALPHAI)*(AF + ni);

//% INTERMEDIATE FIRMS (FOCs)

//% [k]
exp(r) = (1-RHO*GAMMA)*exp(px)*ALPHAF*(exp(xf)/exp(k)) ;

//% [upsf]
KAPPAF/(psiff) =
BETA*( ( (exp(c(+1))/exp(c))^(ALPHA2-1)*((exp(cc(+1))/exp(cc))^(1-ALPHA2-ETA))*(1+tauc)/ (1+tauc) )*
( (1-RHO*GAMMA)*exp(px(+1))*(1-ALPHAF)*( exp(xf(+1)) / exp(nf(+1)) ) -
( 1+exp(taus(+1)) )*exp(wf(+1)) + (1-SIGMAF)*
KAPPAF / (exp(psiff(+1) )) ) ) ;

//% [upsi]
KAPPAI*(psifi) =
BETA*( (exp(c(+1))/exp(c))^(ALPHA2-1)*((exp(cc(+1))/exp(cc))^(1-ALPHA2-ETA))*(1+tauc)/ (1+tauc) )
*( (1-RHO*GAMMA)*exp(px(+1))*(1-ALPHAI)*( exp(xi(+1)) / exp(ni(+1)) ) -
exp(wi(+1)) + (1-RHO-SIGMAI)*KAPPAI / (exp(psifi(+1) )) ) ;

//% [Phillips Curve]
//%pi = BETA*pi(+1) - (1-CHI) * (z -log(ZSS)) ; //%FROM LECTURE'S ASSIGNMENT


//% DEFINITIONS OF MARGINAL UTILITY OF CONSUMPTION
//%vhnf = exp(lambdac) * (1 - taun ) * exp(wf) - BIGPHI*exp(l^(-SMALLPHI)) + exp((1-SIGMAF)* BETA * exp( vhnf(+1) )  ;
//%vhni = exp(lambdac) * exp(wi) - BIGPHI*exp(l^(-SMALLPHI)) + exp(lambdani)*(1-RHO-SIGMAI)*BETA* exp( vhni(+1) ) ;
//%vfnf = (1-RHO*GAMMA)*exp(px)*(1-ALPHAF)*( exp(xf) / exp(nf) ) - ( 1+exp(tauf) )*exp(wf) + (1-SIGMAF)*KAPPAF / (psiff) ;
//%vfni = (1-RHO*GAMMA)*exp(px)*(1-ALPHAI)*( exp(xi) / exp(ni) ) - exp(wi) + (1-RHO-SIGMAI)*KAPPAI / (psifi) ;


//% ***GOVERNMENT***

//% Govt Constraint
//% exp(df) = ( exp(b(+1))*exp(pi(+1)) ) / exp(r) - b ;      

//% Definition of Deficit
exp(surp) + exp(g) + DIGAMMA*exp(uf) - (1-XITR)*exp(tr) - RHO*GAMMA*exp(px)*exp(y) = 0 ;      

//% Total output
exp(y) = exp(xf) + exp(xi);

//% Resource constraint
//% exp (y) = exp(c) * exp(i) + exp(g) + KAPPAF*exp(upsf) + KAPPAI*exp(upsi) +
//% XITR*( (taun+exp(taus)*exp(wf)*exp(nf) + tauc*exp(c) + exp(bigt) ) ;   //%resource constraint


//% ***GOODS MARKET***

//% Aggregate resource constraint
exp(y) = exp(c) + exp(i) + exp(g) + KAPPAF*exp(upsf) + XITR * exp(tr) ;

//% Definition of tax revenue
exp(tr) = (taun+exp(taus))*exp(wf)*exp(nf) + tauc*exp(c) + exp(bigt);

//%     Definition of growth rate of y  (?????????)
//%     exp(gy) = exp(y)/exp(ylag);

//%     definition of lag of output  (?????????)
//%     exp(ylag1) = exp(y);

//%     definition of surp_Y
//%     exp(surp_Y) = exp(surp)/exp(y);

//%     Law of motion of debt-to-gdp
//%     exp(b)-exp(surp_Y)-(exp(b(+1))*exp(pi(+1))*exp(gy(+1)))/exp(rr) = 0 ;


//% ***MONETARY POLICY***

exp(rr) = RRSS * ZETAPI * ( exp(pi) - 1 ) ;


//% ***BARGAINING OVER WAGES***

//% THETAF * ( 1 + exp(taus) ) * exp(vhnf) = (1 - THETAF)*exp(lambdac)*(1 - taun)*exp(vfnf) ;
//% THETAI * exp(vhni) = (1 - THETAI)*exp(lambdac)*exp(vfni) ;

//% Formal sector wage
exp(wf) = THETAF / ( (1-taun)*exp(lambdac) ) * ( BIGPHI*(exp(l)^(-SMALLPHI)) - (1-SIGMAF)*exp(lambdanf) ) +
(1-THETAF)/(1+exp(taus))*( (1-RHO*GAMMA)*exp(px)*(1-ALPHAF) * ( exp(xf) / exp(nf) ) + ( (1-SIGMAF)*KAPPAF/psiff ) ) ;

//% Informal sector wage
exp(wi) = THETAI / exp(lambdac) * ( BIGPHI*(exp(l)^(-SMALLPHI)) - (1-RHO-SIGMAI)*exp(lambdani) ) +
(1-THETAI)*( (1-RHO*GAMMA)*exp(px)*(1-ALPHAI) * ( exp(xi) / exp(ni) ) + ( (1-RHO-SIGMAI)*KAPPAI/psifi ) ) ;



//% ***STOCHASTIC PROCESSES FOR THE SHOCKS***
//external_function(name = tauc);
//%tauc = PT * tauc(-1) + eps_c ;
//external_function(name = taus);
taus = PT * taus(-1) + eps_s ;
//external_function(name = taun);
//%taun = PT * taun(-1) + eps_n ;
g = PT * g(-1) + eps_g ;



end ;
check;
steady;


initval;
cc = 0.1;
c = 0.1;
g = 0.1;
nf = 0.3;
ni = 0.3;
uf = 0;
ui = 0;
u = 0.2;
l = 0.2;
k = 0;
i = 0;
pi = 0;
r = 0;
rr = 0;
lambdac = 1;
lambdanf = 1;
lambdani = 1;
wi = 0;
wf = 0;
px = 0;
y = 0;
xf = 0;
xi = 0;
bigt = 0;
tr = 0;
upsf = 0;
upsi = 0;
y = 0;
mf = 0;
mi = 0;
surp = 0;
taus = 0.16;
end;

write_latex_static_model ;
write_latex_dynamic_model ;


shocks;
//var eps_c ; stderr 0.002;
var eps_s ; stderr 0.002;
//var eps_n ; stderr 0.002;
var eps_g ; stderr 0.002;
end;
 
stoch_simul(order=1,nocorr,nomoments,irf=0);


[jpfeifer]

In the model, there should be as many equations as unknowns.

The calibration strategy is a way of pinning down the parameters. This does not affect which objects are parameters and which are variables. Everything that has a time index is a variable. The problem in their appendix is assigning values to parameters. The standard way of computing a steady state is computing the value of the endogenous variables, given the parameter values. However, you can turn this around. If you fix the value of a variable in steady state, you must have one degree of freedom left for this to be possible. This is achieved by setting one parameter to a value consistent with the restriction on the endogenous variable. The typical example is hours worked and the labor disutility parameter. If you fix labor to 1/3 in steady state, you need to compute the corresponding value of the labor disutility parameter that makes labor 1/3 in steady state given the value of all parameters.

Regarding the error message you get: if you do not use an analytically computed steady state (which you most probably need to do given the calibration strategy the authors use), you need to provide feasible initial conditions. Consumption cannot be 0 in steady state, which gives rise to the error.

[fabioit]

Thank you for the thorough reply.

After realising that I needed to take a few steps back, I am trying to build my first steadystate.m file using a much simpler model. I think I am pretty close from making it work, but there's still a bug in the (steady state) code that I cannot find. I used a file from Matteo Iacoviello website as a template, also writing the parameters on a separate matlab file.
It might have to do with the equation that gets shocked (#10, technology): ln(z) = PSI*ln(z(-1)) + eps;

I get this error message:

Code: Select all

Residuals of the static equations:

Equation number 1 : 27.8043
Equation number 2 : -0.67262
Equation number 3 : 27.9095
Equation number 4 : 0.71626
Equation number 5 : 26.8854
Equation number 6 : 7.028
Equation number 7 : -9.0572
Equation number 8 : 2.1896
Equation number 9 : 2.1481
Equation number 10 : NaN


Error using print_info (line 72)
The steadystate file did not compute the steady state
Error in resid (line 112)
    print_info(info,options_.noprint, options_)
Error in steady (line 90)
        resid;
Error in lse_dmp (line 132)
steady;
Error in dynare (line 180)
evalin('base',fname) ;

The steadystate file is:

Code: Select all


function [ys,check]=lse_dmp_steadystate(junk,ys)

global M_

check=0;

paramfile

nparams = size(M_.param_names,1);
for icount = 1:nparams
    eval(['M_.params(icount) = ',M_.param_names(icount,:),';'])
end

r_ss = 1/BETA - 1 - DELTA;
kn_ss = ( 1/ALPHA * r_ss )^(1/(ALPHA - 1));
p_ss = (1 - ALPHA) * (kn_ss)^ALPHA - W;
g_ss = BETA / ROX * p_ss;
uv_ss = (ZETA/(PHI0*g_ss))^(1/PHI1);
nv_ss = PHI0/ROX*(uv_ss)^PHI1;
v_ss = (uv_ss + nv_ss)^(-1);
n_ss = PHI0/ROX*(uv_ss)^PHI1*v_ss;
y_ss = ((kn_ss)^ALPHA)*n_ss;
u_ss = 1 - n_ss;
k_ss = n_ss*(1/ALPHA*r_ss)^(1/(ALPHA - 1));
c_ss = DELTA*k_ss + ZETA*v_ss - y_ss;

r = r_ss;
p = p_ss;
g = g_ss;
v = v_ss;
u = uv_ss*v;
n = nv_ss*v;
n = n_ss;
k = kn_ss*n;
y = y_ss;
u = u_ss;
k = k_ss;
c = c_ss;

z =  1;


xxx = [ ...
n
k
p
g
c
u
v
r
y
z  ];



logxxx = log(xxx) ;


ys = logxxx ;

I have also tried to use the dog_steadystate.m file that was posted years ago on this forum, but it seemed to work even less than this, which apparently only fails to find the last steady state.

[jpfeifer]

Please use the steady state file from the NK_baseline.mod in the Dynare examples folder as your starting point. Alternatively, use a steady_state_model-block.

[fabioit]

The standard way of computing a steady state is computing the value of the endogenous variables, given the parameter values. However, you can turn this around. If you fix the value of a variable in steady state, you must have one degree of freedom left for this to be possible. This is achieved by setting one parameter to a value consistent with the restriction on the endogenous variable.

How is this implemented in practice, considering I am also using a steadystate.m file?
I need to exogenously set the values of 4 endogenous variables and have 4 parameters to be determined exogenously.
Do I simply change their place in the .mod file(i.e. form "var" to "parameters" and viceversa)? Or is there a more correct way to do it, perhaps working on the steadystate.m file? I couldn't find any answer neither in the forum nor in the user guide.

[jpfeifer]

You simply do this in the steady_state file or steady_state_model-block. Please have a look at https://github.com/JohannesPfeifer/DSGE_mod/blob/master/RBC_news_shock_model/RBC_news_shock_model.mod. This simple RBC model for example calibrates labor to 0.33 by setting the labor disutility parameter psi to the corresponding value