simul returns complex numbers with tiny imaginary part

[dominik]

Hi

when i run a non-stochastic simulation i get no error message but complex numbers in oo_.endo_simul. The imaginary part is very small though (no bigger than 6e-7). rplot fails. Is that a reason to worry or can i just ignore it?
Turning of shock e1 everything works fine.

Many thanks
Dominik

Code: Select all

var z a y c cr cw hw hr dw dr epsilon pi psi nw nr gamma omega bigomega pk lamda m inv g l lw lr k r w mc rk  infl t e b i f ;
var uc v inflstar g_ h_ Vw  CP;
var cshare crshare cwshare invshare capshare bshare lwshare replrate;


varexo e1 e2 e3 e4 e5 e6;

parameters   beta sigma rho gs_ss bs_ss eta delta alpha v1 v2 v3 omega_ss nw_ss nr_ss gamma_ss psi_ss xi  x_ss  theta;
parameters gammamp gammay gammatax gammasmooth calvo x mu;
parameters ssbigomega ssepsilon sscw sscr sspi sslamda sshr sshw sslw sslr ssdr ssdw ssc ssm ssa sst ssl ssr ssi ssinv ssk ssrk ssw sse;


//Structural parameters
beta=0.99;
sigma=1/3;
rho=1-1/sigma;
eta=0.1;
delta=0.05;     
alpha=0.66;
v1=0.64;
v2=0.002;
v3=1-v1-v2;
theta=10;
xi=0.34;   
chi_ss=(1/xi)^v3;
calvo=0.2;


// Demographic parameters
nw_ss=0.003;
nr_ss=nw_ss;
x_ss=0.01;
Tw=50;
Tr=12.5;
T=Tw+Tr;
omega_ss=1- 1/Tw;
gamma_ss=1-1/Tr;
x=x_ss;

psi_ss=(1-omega_ss)/(1+nw_ss-gamma_ss);


// Policy parameters
gammamp=1.5;
gammay=0.0;
gammatax=0.08;
gammasmooth=0.3;


// calibrations
gs_ss=0.18;           
bs_ss=0.7;       
es_ss=0.1038;

param=[gs_ss,beta,sigma,rho,chi_ss,xi,bs_ss,eta,delta,alpha,v1,v2,v3,omega_ss,...
        nw_ss,x_ss,gamma_ss,theta,psi_ss,(theta-1)/theta,0,es_ss,];

unknown0 = [1; 1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1];           
options=optimset('Display','iter','Algorithm','levenberg-marquardt','TolFun',1e-10,'TolX',0);
[unknown] = fsolve(@NKgertlerfun_esfix_grossinvdo,unknown0,options,param);

ssbigomega=    unknown(1);
ssepsilon=     unknown(2);
sscw=          unknown(3);
sscr=          unknown(4);
sspi=          unknown(5);
sslamda=       unknown(6);
sshr=          unknown(7);
sshw=          unknown(8);
sslw=          unknown(9);
sslr=          unknown(10);
ssdr=          unknown(11);
ssdw=          unknown(12);
ssc=           unknown(13);
ssm=           unknown(14);
ssa=           unknown(15);
sst=           unknown(16);
ssl=           unknown(17);
ssr=           unknown(18);
ssi=           unknown(19);
ssinv=         unknown(20);
ssk=           unknown(21);
ssrk=          unknown(22);
ssw=           unknown(23);
sse=           unknown(24);
mu            =sse/(ssw*psi_ss);
ssy=ssl^alpha*(ssk/((1+x_ss)*(1+nw_ss)))^(1-alpha);

//model(bytecode);
 model;




i                = 0.7*i(-1)+.3*(steady_state(i)+gammamp*(infl/steady_state(infl)-1)+gammay*(y/steady_state(y)-1));
1                = calvo*      (1/infl)^(1-theta)+(1-calvo)*inflstar^(1-theta);
v                = calvo*v(-1)*(1/infl)^(-theta) +(1-calvo)*inflstar^(-theta);
h_               = uc*mc*y + calvo * beta *  h_(+1)*(1/infl(+1))^(-theta)*(1+x)*(1+nw(-1)) ;
g_               = uc*inflstar*y + calvo * beta * g_(+1)*(1/infl(+1))^(1-theta)*(inflstar/inflstar(+1))*(1+x)*(1+nw(-1));
theta*h_         = (theta-1)*g_;


pi              = 1-  ( ((1+i(+1))*i/(i(+1)*(1+i)))^(v2*rho) * (w/w(+1)*1/(1+x))^(v3*rho) )^sigma   *(beta^sigma)*(((1+r)*bigomega(+1))^(sigma-1)) * pi/(pi(+1));
epsilon*pi      = 1-  ( ((1+i(+1))*i/(i(+1)*(1+i)))^(v2*rho) * (w/w(+1)*1/(1+x))^(v3*rho) )^sigma   *(beta^sigma)*((1+r)               ^(sigma-1))*gamma * epsilon*pi/(epsilon(+1)*pi(+1));
lamda*a         = lamda(-1)*(a(-1)*omega*(1-epsilon*pi)*(1+r(-1))/((1+nw(-1))*(1+x))) + (1-omega)*a + omega*psi*(dr -epsilon*pi*hr);
k(-1)           = w*(1-alpha)/(alpha*rk)*((1+nw(-1))*(1+x)*l);
Vw              = ((cw^v1*((v2/v1)*(1+i)/i*cw)^v2*(1-l)^v3)^rho /(1-beta*gamma))^(1/rho);
uc              = cw^(rho-1)*Vw^(rho*(1/rho-1));//1;
1+r             = (rk(+1)+(1-delta)*pk(+1))/pk;
k               = (   inv/k(-1) *((1+nw(-1))*(1+x)) - eta/2*( (  inv/k(-1)- (1-(1-delta)/((1+nw_ss)*(1+x)))  ) *((1+nw(-1))*(1+x))  )^2) * k(-1)/((1+nw(-1))*(1+x)) + k(-1)*(1-delta)/((1+nw(-1))*(1+x));
1               = pk*(1-eta* (  inv/k(-1)- (1-(1-delta)/((1+nw_ss)*(1+x)))  )*(1+nw(-1))*(1+x) )  ;
CP              = 1+pk*(   inv/k(-1) *((1+nw(-1))*(1+x)) - eta/2*( (  inv/k(-1)- (1-(1-delta)/((1+nw_ss)*(1+x)))  ) *((1+nw(-1))*(1+x))  )^2)* k(-1)/((1+nw(-1))*(1+x))-inv;

t/y             = sst/steady_state(y)+gammatax*(bshare-bs_ss)+gammasmooth*(bshare-bshare(-1));
b               = (b(-1)+m(-1)/(1+i(-1)))*(  (1+r(-1))/((1+nw(-1))*(1+x))  ) +g+e-t-m;
y               = z*l^alpha*(k/((1+x_ss)*(1+nw)))^(1-alpha)/v;
f               = y*(1-mc);
e/y             = sse/steady_state(y);
y               = g+c+inv;
cr*(1+v2/v1)    = epsilon*pi*((1+r(-1))*lamda(-1)*a(-1)    /((1+nr(-1))*(1+x))+hr);
c               = cw+cr*psi;
cw*(1+v2/v1)    =         pi*((1+r(-1))*(1-lamda(-1))*a(-1)/((1+nw(-1))*(1+x))+hw );
1-lw            = (v3/v1)*(cw/w);
1-lr            = (v3/v1)*(cr/(xi*w));
l               = lw+lr*psi*xi;
m               = ((1+i)/i)*(v2/v1)*c;
dr              = w*lr*xi+e/psi;
dw              = w*lw+f+(CP-1)-t;
hr              = hr(+1)*               gamma*((1+x)/(1+r))+dr;
hw              = hw(+1)*(omega/bigomega(+1))*((1+x)/(1+r))+dw+(1-omega/bigomega(+1))*((1+x)/(1+r))*hr(+1); %%my version
bigomega        = omega(-1)+(1-omega(-1))*(epsilon^(1/(1-sigma)))*(1/xi)^v3;
mc              = 1/z*(w/alpha)^alpha*(rk/(1-alpha))^(1-alpha);
1+i             = (1+r)*infl(+1);
a               = (m/(1+i))+b+k*pk;
b               = bshare*y;

log(z)          = 0.9*log(z(-1))+e2+e4+e5+e6;
psi             = (1-omega(-1))/(1+nw(-1))+gamma(-1)/(1+nw(-1))*psi(-1);
1+nr            = (1-omega)/psi+gamma;
nw          = nw_ss+e1;
gamma           = gamma_ss*(1+e3);
omega           = omega_ss;
g               = gs_ss*steady_state(y);

cshare=c/y;
crshare=cr/y;
cwshare=cw/y;
invshare=inv/y;
capshare=k/y;
lwshare=lw/l;
replrate=e/sigma/psi;
end;
 
// STEADY STATE BLOCK
steady_state_model;
CP=1;
bigomega=ssbigomega   ;
epsilon=ssepsilon    ;
cw=sscw         ;
cr=sscr        ;
pi=sspi         ;
lamda=sslamda      ;
hr=sshr         ;
hw=sshw         ;
lw=sslw         ;
lr=sslr         ;
dr=ssdr         ;
dw=ssdw         ;
c=ssc          ;
m=ssm          ;
a=ssa          ;
t=sst          ;
l=ssl       ;
r=ssr          ;
i=ssi          ;
inv=ssinv        ;
k=ssk          ;
rk=ssrk         ;
w=ssw          ;
e=sse          ;
v               = 1;
gamma           = gamma_ss;
Vw              = ((cw^v1*((v2/v1)*(1+i)/i*cw)^v2*(1-l)^v3)^rho /(1-beta*gamma))^(1/rho);
uc              = cw^(rho-1)*Vw^(rho*(1/rho-1));//1;
z               = 1;
psi             = psi_ss;
nw              = nw_ss;
omega           = omega_ss;
nr=nw_ss;
bshare=bs_ss;
infl=1;
inflstar=1;
pk=1;
mc=(theta-1)/theta;
y               = z*l^alpha*(k/((1+x_ss)*(1+nw)))^(1-alpha)/v;
g              = gs_ss*y;

h_=uc*mc*y / (1-calvo * beta *  (1/infl)^(-theta)*(1+x)*(1+nw));
g_=uc*inflstar*y /(1- calvo * beta *(1/infl)^(1-theta)*(inflstar/inflstar)*(1+x)*(1+nw));
f=(1-mc)*y;
b=bshare*y;
cshare=c/y;
crshare=cr/y;
cwshare=cw/y;
invshare=inv/y;
capshare=k/y;
lwshare=lw/l;
replrate=e/sigma/psi;

end;
 
steady;
check;

//endval;
//e3=0.025691079;
//e1=0.006316351;
//end;

shocks;
//var e1= 0;
var e2= 0.0;
//var e3= 0;
var e4= 0;
var e5= 0;
var e6= 0;


var e1;
periods   1   2   3   4   5   6   7   8   9   10   11   12   13   14   15   16   17   18   19   20   21   22   23   24   25   26   27   28   29   30   31   32   33   34   35   36   37   38   39   40   41   42   43   44   45   ;
values    -3.24862E-07   0.000238555   -0.000225492   -0.001277732   -0.001789681   -0.000800866   -0.003260379   -0.003922711   -0.004217081   -0.00416992   -0.004116136   -0.003850664   -0.004132674   -0.004033812   -0.004603801   -0.004942889   -0.005351577   -0.005752738   -0.006053941   -0.006869655   -0.007390552   -0.007873152   -0.008233119   -0.00895226   -0.009604756   -0.009556714   -0.009698426   -0.009435267   -0.009422349   -0.009447366   -0.009097359   -0.009081987   -0.008592505   -0.007995061   -0.007916783   -0.007542935   -0.007480244   -0.007528256   -0.00741608   -0.007355205   -0.007760983   -0.007464075   -0.007225171   -0.006510635   -0.006316351   ;

var e3;
periods   1   2   3   4   5   6   7   8   9   10   11   12   13   14   15   16   17   18   19   20   21   22   23   24   25   26   27   28   29   30   31   32   33   34   35   36   37   38   39   40   41   42   43   44   45   ;
values    0.000208603   0.000416241   0.000622922   0.000828652   0.001033438   0.002200423   0.003337205   0.004444942   0.005524729   0.00657761   0.007560042   0.008519608   0.009457097   0.01037326   0.011268814   0.012093837   0.012901745   0.013693065   0.0144683   0.015227934   0.015884003   0.016528651   0.017162174   0.017784857   0.018396974   0.018896869   0.019389826   0.019875986   0.02035549   0.020828473   0.021216735   0.021600656   0.021980309   0.022355765   0.022727091   0.023038706   0.023347443   0.02365334   0.023956437   0.024256772   0.024526162   0.024821339   0.02511386   0.025403762   0.025691079   ;
//var e2=0.001;
end;



//stoch_simul(loglinear,order=1,nomoments);
simul(periods=400,maxit=100);

rplot i;rplot y;rplot c;rplot cw;rplot cr; rplot psi;rplot nw; rplot gamma;rplot k;


[jpfeifer]

Use model_diagnostics. There is a singularity problem in your model that should explain the results.

[dominik]

Hi Johannes

you are right yet I dont really know how to interpret the output of model diagnostics. Model diagnostics points out a collinearity in all (!) variables and equation 1, which is the monetary policy rule, which is fairly standard i should say:

i = 0.7*i(-1)+.3*(steady_state(i)+gammamp*(infl/steady_state(infl)-1)+gammay*(y/steady_state(y)-1));

(i is the net interest rate, infl is the gross inflation, y is output in levels, the rule should be the smoothed taylor rule)

Can you help me understand what the warning tells me?

[jpfeifer]

My guess is that usually in New Keynesian models the inflation rate is not endogenously determined and by this token also the nominal interest rate. Thus, you cannot use the steady_state-operator. See the end of http://www.dynare.org/phpBB3/viewtopic.php?f=1&t=5407. Basically, fix the steady state inflation rate. If that does not help, get back to me as there is a bug in simul regarding imaginary numbers: https://github.com/DynareTeam/dynare/pull/775

[dominik]

Hi Johannes
I did as you proposed. i replace steady_state(i) by ssi (which is a parameter to which i assign the value of the steady state nominal IR (given the zero inflation target)) in the monetary rule.
the error message of model diagnostics disappears. this in itself is a surprising feature of model diagnostics, because the steady state block defines steady_state(i) to be equal to iss. so that - i thought - would be the same. but anyway.
more importantly, while the error disappears, everything else remains as it was. i.e. the numbers remain complex.
Many thanks,
Dominik

[dominik]

By the way, another way i found to make the results real, besides as i said setting shock e1 to 0 and only keeping shock e3, is to use this policy rule instead which should be the same in a linear approximation.
1+i = (1+steady_state(i))*infl^gammamp*(y/steady_state(y))^gammay;

[jpfeifer]

The two things are not the same. What you are saying by the steady state operator is that for the dyamic model, steady_state(pi) must be replaced by pi in steady state. For the Jacobian, Dynare then takes the derivative w.r.t. to pi, which is equal to 1. It correctly detects that this steady state value cannot be endogenously computed from the model because for any value of pi there exists a steady state.
In contrast, when steady_state(pi) is replaced by a parameter for pi_ss, the derivative w.r.t. pi is 0. Dynare will then correctly detect that it can now solve for the other occurrences of pi in steady state.
I made your model run by using the current unstable version (November 5, 2014) where the attached file replaces the packaged one.

[dominik]

Dear Joahannes,

thanks for the answer. By just overwriting the old sim1 in the dynare/matlab folder with the new sim1 I got an error
(Reference to non-existent field 'endogenous_terminal_period'.
Error in sim1 (line 35)
endogenous_terminal_period = options_.endogenous_terminal_period;)
and unfortunatley i wasn't able to install dynare unstable on my institution's computer. But anyway, using model(bytecode) i also get a normal solution (which is essentially the same but in real numbers). I guess that will do the trick the same way as the unstable version, right?

I still don't understand the first point, since i thought the steady state command should just be a pointer to the steady state file block where in turn i provide a pointer to the parameter. I never ask dynare to determine a steady state, since i provide it in the steady state file in form of parameters. But anyway, I do understand what i should do to avoid the problem, so that's all that i need to know.

I guess the issue is solved. Once again, many thanks for your immediate and effective help!