Hey everyone
I m trying to replicate a model in a working paper by Garcia Cicco et al. named the real effects of global risk shocks in small open economies. Instead of taking a third orden aproximation of the policy function i am just trying a log.lineal approach. But when I run the .mod file i run into an error at calculating the steady state.
%----------------------------------------------------------------
% 0. Housekeeping (close all graphic windows)
%----------------------------------------------------------------
close all;
warning off;
%----------------------------------------------------------------
% 1. Defining variables
%----------------------------------------------------------------
// Endogenous variables 35
var cc p w l r D cn ct cx cm pn pt px lx lnt yx kx ix qx pi yn kn qn i icn icm in tb pco yr y rer rw re la;
// Exogenous variables (4)
varexo u_rw u_re u_pco u_px;
//PARAMETERS
parameters ZETA NU THETA PHI_D D_ss FI EPSILON CHI ALPHA_X DELTA PHI BETA ALPHA_N GAMMA RHO_rw RHO_re RHO_pco RHO_px yco;
parameters tms_ss r_ss la_ss cc_ss p_ss w_ss l_ss cn_ss ct_ss cx_ss cm_ss pn_ss;
parameters pt_ss px_ss lx_ss lnt_ss yx_ss ix_ss kx_ss qx_ss pi_ss yn_ss kn_ss qn_ss i_ss icn_ss icm_ss in_ss tb_ss pco_ss;
parameters yr_ss y_ss rer_ss D_to_gdpn_ss yn_to_kn_ss yx_to_kx_ss gdpn_ss tb_to_gdpn_ss;
parameters SIGMA_rw SIGMA_re SIGMA_pco SIGMA_px;
%----------------------------------------------------------------
% 2. Calibration ZETA FI se fijan en estado estacionario
%----------------------------------------------------------------
THETA=5;
NU=1.6;
//ZETA=0.00001;
ALPHA_X=0.64;
ALPHA_N=0.36;
GAMMA=0.4;
//FI=0.61;
CHI=0.52;
EPSILON=0.9;
DELTA=0.0051;
PHI_D=0.001;
PHI=6;
yco=8;
//r_ss=0.005;
BETA=0.995;
RHO_rw=0.979;
RHO_re=0.954;
RHO_pco=0.998;
RHO_px=0.99;
SIGMA_rw=0.01;
SIGMA_re=0.01;
SIGMA_pco=0.01;
SIGMA_px=0.01;
// Targeted steady state values
l_ss=0.3;
lx_ss=0.15;
D_to_gdpn_ss=0.5;
pco_ss=1.5;
px_ss=0.89;
//D_ss=1600;
%----------------------------------------------------------------
% 3. Model (30 equations)
%----------------------------------------------------------------
model;
// p*la=(cc - ZETA*NU^(-1)*l^(NU))^(-THETA)
exp(p*la)=(exp(cc)-(ZETA*(exp(l)^(NU))/(NU)))^(-THETA);
// w*la=(cc - ZETA*NU^(-1)*l^(NU))^(-THETA)*-ZETA*l^(NU-1)
exp(w*la)=(exp(cc)-(ZETA*(exp(l)^(NU))/(NU)))^(-THETA)*(NU)*exp(l)^(NU-1);
//la/(1+r)+la*PHI_D*(D-D_ss)=BETA*la(+1)
exp(la)/(1+r)+exp(la)*PHI_D*(D-D_ss)=BETA*exp(la(+1));
//cc=(FI^(1/EPSILON)*(cn^(1-1/EPSILON))+((1-FI)^(1/EPSILON))*(ct^(1-1/EPSILON)))^(EPSILON/(EPSILON-1))
exp(cc)=(FI^(1/EPSILON)*(exp(cn))^(1-1/EPSILON))+((1-FI)^(1/EPSILON))*(exp(ct)^(1-1/EPSILON))^(EPSILON/(EPSILON-1));
//ct=(cx/CHI)^(CHI)*(cm/(1-CHI))^(1-CHI);
exp(ct)=(exp(cx)/CHI)^(CHI)*(exp(cm)/(1-CHI))^(1-CHI);
//cn=FI*(p/pn)^(EPSILON)*cc
exp(cn)=FI*(exp(p)/exp(pn))^(EPSILON)*exp(cc);
//ct=(1-FI)*(p/pt)^(EPSILON)*cc
exp(ct)=(1-FI)*(exp(p)/exp(pt))^(EPSILON)*exp(cc);
//cx=CHI*(pt/px)*ct
exp(cx)=CHI*(exp(pt)/exp(px))*exp(ct);
//cm=(1-CHI)*pt*ct
exp(cm)=(1-CHI)*exp(pt)*exp(ct);
//l=lx+lnt
exp(l)=exp(lx)+exp(lnt);
//PRODUCCION DE BIENES TRANSABLES
//yx=kx(-1)^(ALPHA_X)*lx^(1-ALPHA_X)
exp(yx)=exp(kx(-1))^(ALPHA_X)*exp(lx)^(1-ALPHA_X);
//kx=(1-DELTA)*kx(-1)+ix-ix*PHI/2*(ix/ix(-1)-1)^2
exp(kx)=(1-DELTA)*exp(kx(-1))+exp(ix)-exp(ix)*PHI/2*(exp(ix)/exp(ix(-1))-1)^2;
//px*(1-ALPHA_X)*yx/lx=w
exp(px)*(1-ALPHA_X)*exp(yx)/exp(lx)=exp(w);
//qx=BETA*((la(+1)/la)*(px(+1)*ALPHA_X*yx(+1)/kx+qx(+1)*(1-DELTA)))
exp(qx)=BETA*((exp(la(+1))/exp(la))*(exp(px(+1))*ALPHA_X*exp(yx(+1))/exp(kx)+exp(qx(+1))*(1-DELTA)));
//pi=qx-PHI/2*(ix/ix(-1)-1)^2-PHI*ix/ix(-1)*(ix/ix(-1)-1)+BETA*((la(+1)/la)*qx(+1)*(PHI*(ix(+1)/ix)^2*(ix(+1)/ix-1)))
exp(pi)=exp(qx)-PHI/2*(exp(ix)/exp(ix(-1))-1)^2-PHI*exp(ix)/exp(ix(-1))*(exp(ix)/exp(ix(-1))-1);
//PRODUCCION DE BIENES NO TRANSABLES
//yn=kn(-1)^(ALPHA_X)*lnt^(1-ALPHA_X)
exp(yn)=exp(kn(-1))^(ALPHA_X)*exp(lnt)^(1-ALPHA_X);
//kn=(1-DELTA)*kn(-1)+in-in*PHI/2*(in/in(-1)-1)^2
exp(kn)=(1-DELTA)*exp(kn(-1))+exp(in)-exp(in)*PHI/2*(exp(in)/exp(in(-1))-1)^2;
//pn*(1-ALPHA_X)*yn/lnt=w
exp(pn)*(1-ALPHA_X)*exp(yn)/exp(lnt)=exp(w);
//qn=BETA*((la(+1)/la)*(pn(+1)*ALPHA_X*|n(+1)/kn+qn(+1)*(1-DELTA)))
exp(qn)=BETA*((exp(la(+1))/exp(la))*(exp(pn(+1))*ALPHA_X*exp(yn(+1))/exp(kn)+exp(qn(+1))*(1-DELTA)));
//pi=qn-PHI/2*(in/in(-1)-1)^2-PHI*in/in(-1)*(in/in(-1)-1)+BETA*((la(+1)/la)*qn(+1)*(PHI*(in(+1)/in)^2*(in(+1)/in-1)))
exp(pi)=exp(qn)-PHI/2*(exp(in)/exp(in(-1))-1)^2-PHI*exp(in)/exp(in(-1))*(exp(in)/exp(in(-1))-1);
//CANASTA DE INVERSION
//i=(icn/GAMMA)^(GAMMA)*(icm/1-GAMMA)^(1-GAMMA)
exp(i)=(exp(icn)/GAMMA)^(GAMMA)*(exp(icm)/1-GAMMA)^(1-GAMMA);
//i=in+ix
exp(i)=exp(in)+exp(ix);
//icn=GAMMA^(pi/pn)*i
exp(icn)=GAMMA^(exp(pi)/exp(pn))*exp(i);
//icm=(1-GAMMA)*(pi)*i
exp(icm)=(1-GAMMA)*(exp(pi))*exp(i);
//AGREGACIÓN
//yn=cn+icn
exp(yn)=exp(cn)+exp(icn);
//tb=pco*yco+px*(yx-cx)-(cm+icm)
tb=exp(pco)*exp(yco)+exp(px)*(exp(yx)-exp(cx))-(exp(cm)+exp(icm));
//D(-1)+(PHI_D/2)*(D-D_ss)^2=D/(1+r)+tb;
D(-1)+(exp(PHI_D)/2)*(D-D_ss)^2=D/(1+r)+tb;
//yr=yx+yn
exp(yr)=exp(yx)+exp(yn);
//y=yr+yco
exp(y)=exp(yr)+exp(yco);
//rer=1/pt
exp(rer)=1/exp(pt);
// Exogenous
//r=r_ss+rw+re
r=r_ss+rw+re;
//log(rw)=RHO_rw*log(rw(-1))
log(rw)=RHO_rw*log(rw(-1))+u_rw;
//log(re)=RHO_re*log(re(-1))
log(re)=RHO_re*log(re(-1))+u_re;
//log(pco/pco_ss)=RHO_pco*log(pco(-1))/(pco_ss))
log(exp(pco)/pco_ss)=RHO_pco*log(exp(pco(-1))/pco_ss)+u_pco;
//log(px/px_ss)=RHO_px*log(px(-1))/(px_ss))
log(exp(px)/px_ss)=RHO_px*log(exp(px(-1))/px_ss)+u_px;
end;
%----------------------------------------------------------------
% 4. Steady State
%----------------------------------------------------------------
steady_state_model;
r_ss=1/BETA-1;
lnt_ss=l_ss-lx_ss;
pn_ss=((px_ss)^(1/(1-ALPHA_X))*(1-ALPHA_X)*((1/BETA-(1-DELTA))/ALPHA_N)^(ALPHA_N/(1-ALPHA_N))/(1-ALPHA_N)*((1/BETA-(1-DELTA))/ALPHA_X)^(ALPHA_X/(1-ALPHA_X)))^((1-ALPHA_N)*(1-ALPHA_X)/(1-ALPHA_X+GAMMA*(ALPHA_X-ALPHA_N)));
pi_ss=(pn_ss)^GAMMA;
qx_ss=pi_ss;
qn_ss=pi_ss;
pt_ss=px_ss^CHI;
yx_to_kx_ss=pi_ss/px_ss*(1/BETA-(1-DELTA))/ALPHA_X;
yn_to_kn_ss=pi_ss/pn_ss*(1/BETA-(1-DELTA))/ALPHA_N;
kx_ss=lx_ss*(yx_to_kx_ss)^(-1/(1-ALPHA_X));
kn_ss=lnt_ss*(yn_to_kn_ss)^(-1/(1-ALPHA_N));
yx_ss=kx_ss*yx_to_kx_ss;
yn_ss=kn_ss*yn_to_kn_ss;
w_ss=px_ss*(1-ALPHA_X)*yx_ss/lx_ss;
ix_ss=DELTA*kx_ss;
in_ss=DELTA*kn_ss;
i_ss=in_ss+ix_ss;
icm_ss=(1-GAMMA)*pi_ss*i_ss;
icn_ss=GAMMA*(pi_ss/pn_ss)*i_ss;
cn_ss=yn_ss-icn_ss;
yr_ss=yn_ss+yx_ss;
y_ss=yr_ss+yco;
gdpn_ss=pn_ss*yn_ss+px_ss*yx_ss+pco_ss*yco;
tb_to_gdpn_ss=D_to_gdpn_ss*r_ss/(1+r_ss);
FI=(1-(tb_to_gdpn_ss*gdpn_ss+icm_ss-pco_ss*yco-px_ss*yx_ss)/(((pt_ss)^(1-EPSILON))*((pn_ss)^(EPSILON))*cn_ss))^(-1);
p_ss=(FI*(pn_ss)^(1-EPSILON)+(1-FI)*(pt_ss)^(1-EPSILON))^(1/(1-EPSILON));
cc_ss=cn_ss/FI*(pn_ss/p_ss)^(EPSILON);
ct_ss=(1-FI)*(p_ss/pt_ss)^(EPSILON)*cc_ss;
cx_ss=CHI*(pt_ss/px_ss)*ct_ss;
cm_ss=(1-CHI)*(pt_ss)*ct_ss;
tms_ss=-p_ss/w_ss;
//ZETA=cc_ss^(-THETA)*w_ss/(l_ss^(NU)*p_ss);
ZETA=p_ss*(l_ss)^(NU-1)/w_ss;
la_ss=(cc_ss-ZETA*(NU)^(-1)*l_ss^(NU))^(-THETA)/p_ss;
rer_ss=1/p_ss;
tb_ss=tb_to_gdpn_ss*gdpn_ss;
D_ss=D_to_gdpn_ss*gdpn_ss;
// Initial values for solver
la=log(la_ss);
cc=log(cc_ss);
w=log(w_ss);
l=log(l_ss);
p=log(p_ss);
r=r_ss;
D=D_ss;
cn=log(cn_ss);
ct=log(ct_ss);
cx=log(cx_ss);
cm=log(cm_ss);
pn=log(pn_ss);
pt=log(pt_ss);
lx=log(lx_ss);
lnt=log(lnt_ss);
y=log(y_ss);
yx=log(yx_ss);
kx=log(kx_ss);
ix=log(ix_ss);
qx=log(qx_ss);
pi=log(pi_ss);
yn=log(yn_ss);
kn=log(kn_ss);
in=log(in_ss);
qn=log(qn_ss);
i=log(i_ss);
icn=log(icn_ss);
icm=log(icm_ss);
tb=tb_ss;
pco=log(pco_ss);
px=log(px_ss);
yr=log(yr_ss);
rer=log(rer_ss);
rw=0;
re=0;
end;
%----------------------------------------------------------------
% 4. Computation
%----------------------------------------------------------------
steady;
check;
shocks;
var u_rw = SIGMA_rw^2;
var u_re = SIGMA_re^2;
var u_pco = SIGMA_pco^2;
var u_px = SIGMA_px^2;
end;
stoch_simul(periods=0, irf = 1000, order = 1, nograph);
save replica.mat;
Any comments will be much appreciated
Thanks