# Different results from using dynare initval-block and using external Matlab program

Hi, everyone

different results from using dynare initval-block and using external Matlab program, while same analytical solutions of all endogenous variables.

That means your entered equations and the computed steady state values are inconsistent. There must be a mistake somewhere.

When I using external Matlab program to compute the steady state values, it turn out that the steady state contain NaN.

But when I enter the analytical steady state( as same as external Matlab program ) into Dynare initval-block, it is all ok, no error!

Sorry ，I can’t understand your suggestion.

My problem:

In my external Matlab program,

``````function [ys,check] = NK_baseline_steadystate(ys,exo)
% computes the steady state for the NK_baseline.mod and uses a numerical
% solver to do so
% Inputs:
%   - ys        [vector] vector of initial values for the steady state of
%                   the endogenous variables
%   - exo       [vector] vector of values for the exogenous variables
%
% Output:
%   - ys        [vector] vector of steady state values for the the endogenous variables
%   - check     [scalar] set to 0 if steady state computation worked and to
%                    1 of not (allows to impos restriction on parameters)

global M_

% read out parameters to access them with their name
NumberOfParameters = M_.param_nbr;
for ii = 1:NumberOfParameters
paramname = deblank(M_.param_names(ii,:));
eval([ paramname ' = M_.params(' int2str(ii) ');']);
end
% initialize indicator
check = 0;

%% Enter model equations here

options=optimset(); % set options for numerical solver

% the steady state computation follows FV (2006), section 4.1
PI=PIbar;
u=1;
q=1;
d=1;
phi=1;
m=0;
zeta=1;
LambdaYd= (LambdaA+alppha*Lambdamu)/(1-alppha);
mu_z=exp(LambdaYd);
mu_I=exp(Lambdamu);
mu_A=exp(LambdaA);

%set the parameter Lambdax
Lambdax=mu_z;

%set the parameter gammma1
gammma1=mu_z*mu_I/betta-(1-delta);
if gammma1<0 % parameter violates restriction; Preventing this cannot be implemented via prior restriction as it is a composite of different parameters and the valid prior region has unknown form
check=1; %set failure indicator
return; %return without updating steady states
end

r=1*gammma1;
R=1+(PI*mu_z/betta-1);

%set Rbar
Rbar=R;

PIstar=((1-thetap*PI^(-(1-epsilon)*(1-chi)))/(1-thetap))^(1/(1-epsilon));
PIstarw=((1-thetaw*PI^(-(1-chiw)*(1-eta))*mu_z^(-(1-eta)))/(1-thetaw))^(1/(1-eta));

mc=(epsilon-1)/epsilon*(1-betta*thetap*PI^((1-chi)*epsilon))/(1-betta*thetap*PI^(-(1-epsilon)*(1-chi)))*PIstar;
w=(1-alppha)*(mc*(alppha/r)^alppha)^(1/(1-alppha));
wstar=w*PIstarw;
vp=(1-thetap)/(1-thetap*PI^((1-chi)*epsilon))*PIstar^(-epsilon);
vw=(1-thetaw)/(1-thetaw*PI^((1-chiw)*eta)*mu_z^eta)*PIstarw^(-eta);
tempvaromega=alppha/(1-alppha)*w/r*mu_z*mu_I;

[ld,fval,exitflag]=fzero(@(ld)(1-betta*thetaw*mu_z^(eta-1)*PI^(-(1-chiw)*(1-eta)))/(1-betta*thetaw*mu_z^(eta*(1+gammma))*PI^(eta*(1-chiw)*(1+gammma)))...
-(eta-1)/eta*wstar/(varpsi*PIstarw^(-eta*gammma)*ld^gammma)*((1-h*mu_z^(-1))^(-1)-betta*h*(mu_z-h)^(-1))*...
((mu_A*mu_z^(-1)*vp^(-1)*tempvaromega^alppha-tempvaromega*(1-(1-delta)*(mu_z*mu_I)^(-1)))*ld-vp^(-1)*Phi)^(-1),0.25,options);
if exitflag <1
%indicate the SS computation was not sucessful; this would also be detected by Dynare
%setting the indicator here shows how to use this functionality to
%filter out parameter draws
check=1; %set failure indicator
return; %return without updating steady states
end

l=vw*ld;
k=tempvaromega*ld;
x=(1-(1-delta)*(mu_z*mu_I)^(-1))*k;
yd=(mu_A/mu_z*k^alppha*ld^(1-alppha)-Phi)/vp;
c=(mu_A*mu_z^(-1)*vp^(-1)*tempvaromega^alppha-tempvaromega*(1-(1-delta)*(mu_z*mu_I)^(-1)))*ld-vp^(-1)*Phi;
lambda=(1-h*betta*mu_z^(-1))*(1-h/mu_z)^(-1)/c;
F=yd-1/(1-alppha)*w*ld;
f=(eta-1)/eta*wstar*PIstarw^(-eta)*lambda*ld/(1-betta*thetaw*mu_z^(eta-1)*PI^(-(1-chiw)*(1-eta)));
f2=varpsi*d*phi*PIstarw^(-eta*(1+gammma))*ld^(1+gammma)/(1-betta*thetaw*(PI^chiw/PI)^(-eta*(1+gammma))*(wstar/wstar*mu_z)^(eta*(1+gammma)));

g1=lambda*mc*yd/(1-betta*thetap*PI^((1-chi)*epsilon));
g2=epsilon/(epsilon-1)*g1;

V=1/(1-betta)*(log((1-h)*c)-varpsi*l^(1+gammma)/(1+gammma));

%% end own model equations

for iter = 1:length(M_.params) %update parameters set in the file
eval([ 'M_.params(' num2str(iter) ') = ' M_.param_names(iter,:) ';' ])
end

NumberOfEndogenousVariables = M_.orig_endo_nbr; %auxiliary variables are set automatically
for ii = 1:NumberOfEndogenousVariables
varname = deblank(M_.endo_names(ii,:));
eval(['ys(' int2str(ii) ') = ' varname ';']);
end
``````

The results:

``````ans =
1 to 8 column
2.6112    3.0029    0.3369    0.0251    0.0200    0.0395    1.0000    1.0000
9 to 16 column
0.1192    1.0000    6.0000    1.0382    1.0923    1.3571    3.4113    0.3679
17 to 24 column
4.7675    4.7675    0.6654    0.0534    3.1816    4.6209    0.5701    0.1140
25 to 30 column
1.0000    1.0025    0.1500    0.0500    0.1250 -478.3098
``````

but when I use initval-block in dynare, and the same steady state solutions are entered

``````initval;
A=1;
Z=1;
v=1;
pssi=psis;
tauc=taucs;
taun=tauns;
tauk=tauks;
pi=pi_target;
i=1/betta*(1+pi)-1;
u=1;
R=(1+tauc)/(1-tauk)*(1/betta-(1-delta0));
pi_star=(((1+pi)^(1-epsilonp)-fip*(1+pi)^(etap*(1-epsilonp)))/(1-fip))^(1/(1-epsilonp))-1;
v_p=(1-fip)*((1+pi_star)/(1+pi))^(-epsilonp)/(1-fip*(1+pi)^(epsilonp*(1-etap)));
mc=(epsilonp-1)/epsilonp*(1+pi_star)/(1+pi)*(1-fip*betta*(1+pi)^(epsilonp*(1-etap)))/(1-fip*betta*(1+pi)^((1-epsilonp)*(etap-1)));
w=(1-alfa)*mc*(alfa*mc/R)^(alfa/(1-alfa));
w_star=w*((1-fiw*(1+pi)^((epsilonw-1)*(1-etaw)))/(1-fiw))^(1/(1-epsilonw));
N_d=((epsilonw-1)/epsilonw/psis*(1-taun)/(1+tauc)/(1-b)*(1-betta*b)*((1-omegga)*(alfa*mc/R)^(alfa/(1-alfa))/v_p-delta0*(alfa*mc/R)^(1/(1-alfa)))^(-1)*(w/w_star)^(-epsilonw*ka)*((1-fiw*betta*(1+pi)^(epsilonw*(1+ka)*(1-etaw)))/(1-fiw*betta*(1+pi)^((epsilonw-1)*(1-etaw)))))^(1/(1+ka));
K_hat=(alfa*mc/R)^(1/(1-alfa))*N_d;
K=K_hat;
I=delta0*K;
Y=w*N_d+R*K;
G=omegga*Y;
C=Y-G-I;
lamb=1/C*(1-betta*b)/(1-b)/(1+tauc);
mu=(1+tauc)*lamb;
x1=lamb*mc*Y/(1-fip*betta*(1+pi)^((1-etap)*epsilonp));
x2=lamb*Y/(1-fip*betta*(1+pi)^((etap-1)*(1-epsilonp)));
h1=psis*(w/w_star)^(epsilonw*(1+ka))*N_d^(1+ka)/(1-fiw*betta*(1+pi)^(epsilonw*(1+ka)*(1-etaw)));
h2=(1-taun)*lamb*(w/w_star)^epsilonw*N_d/(1-fiw*betta*(1+pi)^((epsilonw-1)*(1-etaw)));
Welfare=1/(1-betta)*(log(C-b*C)-pssi*N_d^(1+ka)/(1+ka));
end;
``````

In `oo_.steady_state` file, steadystate values:

``````1.51302016047266
1.73997318454356
0.581426577621927
0.0251256281407035
0.0200000000000000
0.0394615936826992
1.00000000000000
1
0.122748716814203
1
6
1.03824378411395
1.09232823900298
1.48241405128077
2.03567156553141
0.378889958696544
4.90994867256813
4.90994867256812
0.665393362613623
0.0534454955762316
2.84632657916216
4.13393923502041
0.880219117262940
0.176043822826811
1
1.00246304513304
0.150000000000000
0.0500000000000000
0.125000000000000
-371.005689992056
``````

results in

``````Equation number 8 : 0.49559 : ·½³Ì8£º×îÓÅ¹¤×Ê¶¨¼Û
Equation number 16 : -0.28528 : ·½³Ì16£ºÉú²úº¯Êý
``````

That means your steady state file is inconsistent with the entered equations.

It’s done well!
Thank you very very much!

all of Residuals is zero.

My steadystate file is not named correctly!

Thanks again!