Dynare forums

Dear All,
I would like to use the Dynare 2nd & 3rd order solvers inside a MATLAB function that takes model parameters and data as inputs, in order to compute decision rules that are then used to generate the likelihood (or some other function of the parameters and data). The key issue is that, in order to permit maximization of the likelihood with respect to the model parameters, the decision rules have to be computed INSIDE the Matlab function. (NB I do not want to use the particle filter-based methods available in Dynare.)

For second-order accurate decision rules, this can be done very easily and efficiently, using set_param_value and resol (inside the MATLAB function). Please see simplified code I) below, where MODEL_RBC.mod is a model with two parameters (RISKAV & ELS).

However, this does NOT work for third-order approximations! It seems that, unfortunately, set_param_value cannot be used to pass new parameter values to k_order_pert.

What works for third-order approximations: run the .mod file inside the Matlab function, for each new parameter vector. But this is VERY slow (see simplified code II below).

Is there any other faster method that could be used to compute the 3rd order decision rule repeatedly inside a Matlab function?

I think that this issue is of great interest for all users who want to write code for estimators (or for other optimization problems w.r.t. model parameters) based on higher order Dynare/Dynare++ DSGE model solutions.

Thank you very much for any feedback and advice!
Best, Robert Kollmann; http://www.robertkollmann.com
= = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = =

I) CODE (simplified ) THAT WORKS ONLY FOR SECOND-ORDER MODELS
1) MATLAB script file
clear all;
dynare MODEL_RBC.mod
fminsearch(MODEL_LL,x0)

2) MODEL_RBC.mod (steady state is computed in closed form using steady_state_model)
RISKAV=10; ELS=4; //any values can be set here (these values are only used when MODEL_RBC.mod is run for the first time in the script file
stoch_simul(periods=0,order=2,pruning,noprint,nograph,nocorr,…
nofunctions,nomoments,IRF=0);

3)Matlab Function that computes likelihood (L) for parameter vector PARAM
function [L]=MODEL_L(PARAM,DATA)
global oo_ M_ options_
RISKAV1=PARAM(1); % risk aversion
ELS1=PARAM(2); % elasticity of labor supply
set_param_value('RISKAV',RISKAV1);
set_param_value('ELS',ELS1);

[dr,info,M2_,options2_,oo2_] = resol(0,M_,options_,oo_);

//UNFORTUNATELY using k_order_pert here does NOT work
//[dr,info] = k_order_pert(oo_.dr,M_,options_);


= = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = =
II. CODE (simplified ) THAT WORKS FOR (SECOND AND) THIRD ORDER MODELS, BUT IS MUCH SLOWER. IS THERE A MORE EFFICIENT SET-UP?
1) MATLAB script file
clear all;
fminsearch(MODEL_LL,x0)

2) MODEL_RBC.mod
load PARAM1
set_param_value('RISKAV',RISKAV1);
set_param_value('ELS',ELS1);
stoch_simul(periods=0,order=3,pruning,noprint,nograph,nocorr,…
nofunctions,nomoments,IRF=0);

3)
function [L]=MODEL_L(PARAM,DATA)
global oo_ M_ options_

RISKAV1=PARAM(1); % risk aversion
ELS1=PARAM(2); % elasticity of labor supply
save PARAM1 RISKAV1 ELS1;
dynare MODEL_RBC.mod noclearall

Dear Robert, please take a look at the replication files to Born/Pfeifer (2014): "Risk Matters: A comment" at https://sites.google.com/site/pfeiferecon/20140525_replication_codes.zip?attredirects=0. There, we use SMM to estimate a third-order approximated model.

Thanks for your feedback, Johannes!
I just studied the code of your 2014 AER paper (congratulations!).

Please correct me if I am wrong, but it seems to me that your smm_diff_function.m function (that you employ for estimation) uses a second-order model approximation, not a third-order approximation.
smm_diff_function.m uses 'resol' to compute the model solution in Dynare, for the 'xopt' parameter vector with respect to which smm_diff_function is miniminzed.
Yet, as far as I understand it, resol can only compute first and second-order model approximations (see below).

What method would you use if you wanted to use a third-order model approximation for your SMM estimation ?
(Other than running the .mod file inside smm_diff_function.m ?)

Thank you for your kind attention.
Best regards, Robert

==================================
NB The resol.m file (and the stochastic_solvers.m file called by resol.m) explicitly state that these procs only handle first and second-order approximations.

To check that resol cannot handle third-order approximations, I ran a *.mod file with these commands:

"stoch_simul(order=3,pruning,periods=0,irf=0,nofunctions);
set_param_value('ALP', 0.8);
[oo_.dr,info,M_,options_,oo_] = resol(0,M_,options_,oo_); "

This produces the following error message:

" ??? Undefined function or variable "jacobia_".

Error in ==> stochastic_solvers at 118
[infrow,infcol]=find(isinf(jacobia_));

Error in ==> resol at 137
[dr,info] = stochastic_solvers(dr,check_flag,M,options,oo);

Error in ==> MODEL_RBC at 196
[oo_.dr,info,M_,options_,oo_] = resol(0,M_,options_,oo_);

Error in ==> dynare at 180
evalin('base',fname) ; "

Dear Robert,
we used a third order approximation (as we also did in our JME paper where we also did SMM at third order). Otherwise, uncertainty shocks would not matter at all.
Which Dynare version are you using?
resol calls k_order_pert, which computes the third-order decision rules. I need to change the headers of the functions you reference. They are wrong.
The problem when calling resol with third order is that in some versions you need to manually set

Code: Select all

options_.qz_criterium = 1+1e-6;

as it is locally set in stoch_simul and the absence of this option leads to crashes when directly calling resol.

If it does not work, please send me your code via email.

Dear Johannes,
thank you very much for your rapid answer!
I use Dynare version 4.4.3. and MATLAB Version 7.12.0.635 (R2011a)

Please find below (and in attachment) a simple .mod file that illustrates the problem in using resol with third order.

See last lines of the file:
"stoch_simul(periods=0,order=3,pruning,noprint,nograph,nocorr,nofunctions,nomoments,IRF=0);
set_param_value('ALP',.5);
[oo2_.dr,info2,M2_,options2_,oo2_] = resol(0,M_,options_,oo_); "

When I set order=2 in stoch_simul, then resol works correctly!

But for order=3 I get the following error message

---
??? Undefined function or variable "jacobia_".

Error in ==> stochastic_solvers at 118
[infrow,infcol]=find(isinf(jacobia_));

Error in ==> resol at 137
[dr,info] = stochastic_solvers(dr,check_flag,M,options,oo);

Error in ==> MODEL_RBC_TEST at 180
[oo2_.dr,info2,M2_,options2_,oo2_] = resol(0,M_,options_,oo_);

Error in ==> dynare at 180
evalin('base',fname) ;
----

I it tried putting "options_.qz_criterium = 1+1e-6;" before stoch_simul, but this still led to the same error.

So, I guess that "options_.qz_criterium = 1+1e-6;" has to be inserted somewhere in resol.m or in k_order_pert, or somewhere els?
Could you please tell me exactly where the change has to be made?

I am reluctant to changing these files without your advice, as this might lead to other problems.

If explainig/making the changes is too complicated, perhaps you could simply post (or send me) your correct version of resol?

It would also be great if this issue could be fixed in the next Dynare release.

Thank you very much to you and your colleagues in the Dynare team for the great work!
And thank you very much for your further feedback!

Best regards, Robert

========================================================================================================
ATTACHED FILE: MODEL_RBC_TEST.mod
// robert_kollmann@yahoo.com
// c:\dynare\workEZW\MODEL_RBC_TEST.mod Jan21 2015

// Simple RBC model with shocks to: TFP (th), Gov't purchases (g) and time preference (lam)
// u=(1/(1-RISKAV)) * C^(1-RISKAV) - MU * (1/(1+1/ELS))*lab^(1+1/ELS)
// ELS: labor supply elasticity. y=th*(kap^(1-ALP))*(lab^ALP);
// Normalization: set THss & MU such that in steady state y=lab=1
// Express variables in logs: Ly is log(GDP) etc.

var Ly Lc Linvest Llab Lkap Lth Lg Llam;
varexo Eth Eg Elam;
parameters ALP DEPREC BETTA RISKAV ELS Gss RHOth RHOg RHOlam SSth SSg SSlam;

ALP=.7; DEPREC=.025; BETTA=0.99; Gss=.2; RISKAV=10; ELS=4; RHOth=.99;
RHOg=.99; RHOlam=.99; SSth=.10; SSg=.05; SSlam=.0025;
//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
model;
exp(Ly)=exp(Lth)*(exp(Llab)^ALP)*(exp(Lkap(-1))^(1-ALP));
exp(Lkap)=exp(Lkap(-1))*(1-DEPREC)+exp(Linvest);
# Kss=(1-ALP)/(DEPREC+(1-BETTA)/BETTA);
# THss=Kss^(ALP-1);
# Css=1-Gss-DEPREC*Kss;
# MU=ALP*Css^-RISKAV;
(ALP*exp(Ly)/exp(Llab))*(exp(Lc)^-RISKAV)=MU*exp(Llab)^(1/ELS);
exp(Ly)=exp(Lc)+exp(Linvest)+exp(Lg);
1=BETTA*exp(Llam)*((exp(Lc(+1))/exp(Lc))^-RISKAV)*((1-ALP)*exp(Ly(+1))/exp(Lkap) + 1-DEPREC);
Lth-log(THss)=RHOth*(Lth(-1)-log(THss)) + SSth*Eth;
Lg-log(Gss)=RHOg*(Lg(-1)-log(Gss)) + SSg*Eg;
Llam=RHOlam*Llam(-1) + SSlam*Elam;
end;
//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

steady_state_model;
Ly=0;
Lkap=log((1-ALP)/(DEPREC+(1-BETTA)/BETTA));
Lg=log(Gss);
Lc=log(1-exp(Lg)-DEPREC*exp(Lkap));
Linvest=log(DEPREC*exp(Lkap));
Llab=0;
Lth=log(exp(Lkap)^(ALP-1));
Llam=0;
end;

//steady;
//check;
//model_diagnostics;

shocks;
var Eth = 1;
var Eg = 1;
var Elam = 1;
end;
//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
//options_.qz_criterium = 1+1e-6
stoch_simul(periods=0,order=3,pruning,noprint,nograph,nocorr,nofunctions,nomoments,IRF=0);
//options_.qz_criterium = 1+1e-6
set_param_value('ALP',.5);
[oo2_.dr,info2,M2_,options2_,oo2_] = resol(0,M_,options_,oo_);

This seems to be a bug. Although you are using order=3, the k_order_solver option is not enabled. If you add

Code: Select all

k_order_solver

to the options of stoch_simul, it will run.

Dear Johannes,
thank you very much!
resol works for the third-order approximation when the k_order_solver option is included in stoch_simul !

I find it mysterious that in your AER code, resol could handle the third-order approximation, although you did not use the k_order_solver option. It would be great if you and your colleagues could fix thisin the next version of dynare (or please mention it in the manual).

I had a terrible 'headache' during the whole of last week about not being able to use resol for the third-order. But fortunately that problem is solved now.

Congratulations to you for knowing the code so well! And thanks for your generosity, and for responding so fast.

Best wishes, Robert

We will deal with the issue. If order==3, there are two places where the k_order_solver option is set. In stoch_simul and simult_. However, stoch_simul does not globally set this option, which explains the problems in your case. In our replication code, we directly call simult_ so that the option is correctly set afterward.

First, thanks for this post. I was also having problem with the jacobia_ issue.

Further, I am also trying to implement an optimization routine and liked the idea that I don't have to run Dynare prepocessor again.
I am running into a strange issue, which is related more to using resol.m (less to the optimization).

I am explaining my problem using a toy example (the relevant files are also attached). Let's say I am estimating beta and gamma.
1. I run the Rec1.mod file with beta=0.95 and gamma=2. All results are fine.
2. I set the initial value beta=0.97 and gamma=1.5.
3. Then from smm.m I call objfn.m. Here I use resol.m using the M_, oo_ etc.
4. I get an error 'dynare:k_order_perturbation: Caught Kord exception: The model is not Blanchard-Kahn stable.
5. Surprisingly, if I go back and use beta=0.97 and gamma=1.5 in Rec1.mod file there are no errors.
Ideally if I am using same values for beta and gamma I shouldn't get any error/warning (i.e. in the first run of the optimizer).

I have been able to trace the error to k_order_pert.m where the cc code is throwing err=1.

I am not showing the results further on in the optimization, as they are not related to the problem. My concern is that getting this error will be fine if the model violates BK conditions - but in this case the model for same parameter value is solved through Rec1.mod file. So when optimizing chances might be there that obj function will not be evaluated for certain parameter set (which is actually solvable).

I referred to replication codes in Benjamin Born and Johannes Pfeifer (2014).

Please help me in understanding if this a bug or I am committing some error.

Thank you.

Just to update.
I am using Dynare version 4.4.3 and Matlab R2013b.

Infact when I used Dynare Version 4.3.3, mexErrCheck failed.
Error using mexErrCheck (line 42)
Error encountered in: k_order_perturbation.
This makes sense as in 4.4.3, dynare saves err as info and doesnt raise an error.

That is because your mod-file does not handle parameter dependence correctly. theta is a function of GAM. Use a model-local variable for theta and everything works:

Code: Select all

var a c k v vk;
varexo epsilon;

predetermined_variables k;

parameters    SIG DELTA ALFA BETTA RHO GAM SIGZ;
BETTA=0.95; %discount rate
DELTA=0.1; %depreciation rate
ALFA=0.3; %capital share
RHO=0.95; %persistence of technology shock
SIG=0.5; %intertemporal elasticity of substitution
GAM=2;
SIGZ = 0.005;

model;
#theta  = (1-GAM)/(1-(SIG));
0 = exp(c) + exp(k(+1)) - (1-DELTA) * exp(k) - exp(a) * exp(k)^ALFA;
0 = exp(c)^(-SIG) - BETTA * exp(c(+1))^(-SIG) *  (exp(v(+1))^(SIG-GAM) / exp(vk)^(1-1/theta)) * (exp(a(+1)) * ALFA * exp(k(+1))^(ALFA-1) + 1 - DELTA);
0 = a - RHO * a(-1) - SIGZ*epsilon;
0 = exp(v)^(1-SIG) - exp(c)^(1-SIG) - BETTA* exp(vk)^(1/theta);
0 = exp(vk) - exp(v(+1))^(1-GAM);
end;

steady_state_model;
k = log(((1/BETTA+DELTA-1)/ALFA)^(1/(ALFA-1)));
c = log(exp(k)^(ALFA)-DELTA*exp(k)); %steady-state value of consumption
v = (1/(1-SIG))*log(((1/(1-BETTA)) * (exp(c)^(1-SIG))));
vk = (1-GAM)*log(exp(v));
a = 0;
end;

shocks;
var epsilon; stderr 1;
end;

steady(solve_algo=2, maxit=1000, nocheck);
check;
options_.k_order_solver=1; %this is important as then for order=3 there will be jacobia_ not found error.
set_dynare_seed(201602)
stoch_simul(order=3, periods=20, drop=0, nodisplay,noprint,nograph,nocorr,nofunctions,nomoments);


Thank you Johannes!

Hi Johannes,

Quick question in regards with using model-local variables. I understand their scope is restricted to model block only. However, in some cases I need to use this model-local variable in steady state model block as well. But I am not able to use model-local variables there. Is there any way I can do that ?

Thanks.

If those model-local variables are parameters you can define them as actual parameters instead of model-local variables and then update them in the steady_state_model-block.

Thanks ! Gotta it.

For example, I was focusing on using 'theta' example I posted in this post. You mentioned I should use theta as model-local variable for handling parameter dependency correctly. So I was wondering if there is any way of using theta in steady_state_model block in this situation (I am not defining theta as a parameter.).