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.

Thanks in advance!

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

Thanks for your reply!

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!

Your steady state file does have the recommended structure. Please see https://git.dynare.org/Dynare/dynare/blob/master/examples/NK_baseline_steadystate.m

Sorry ,I can’t understand your suggestion.

My problem:

In my external Matlab program,

function [ys,check] = NK_baseline_steadystate(ys,exo)
% 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

A proper steady state file dsge_steadystate.m (2.5 KB)
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!