Commit aa6ff983 authored by sebastien's avatar sebastien
Browse files

AIM solver:

* added interface in the MOD file, with a new option "aim_solver" to stoch_simul and estimation
* documented the option in the reference manual


git-svn-id: https://www.dynare.org/svn/dynare/trunk@3059 ac1d8469-bf42-47a9-8791-bf33cf982152
parent 217462af
......@@ -1692,6 +1692,11 @@ The simulated endogenous variables are available in global matrix <varname>oo_.e
<term><option>solve_algo</option> = <replaceable>INTEGER</replaceable></term>
<listitem><para>See <link linkend="solve_algo">there</link> for the possible values and their meaning</para></listitem>
</varlistentry>
<varlistentry id="aim_solver">
<term><option>aim_solver</option></term>
<listitem><para>Use the Anderson-Moore Algorithm (AIM) to compute the decision rules, instead of using Dynare's default method based on a generalized Schur decomposition. This option is only valid for first order approximation. See <ulink url="http://www.federalreserve.gov/Pubs/oss/oss4/aimindex.html">AIM website</ulink> for more details on the algorithm.</para>
</listitem>
</varlistentry>
</variablelist>
</refsect1>
......@@ -2425,6 +2430,10 @@ end;
<term><option>irf</option> = <replaceable>INTEGER</replaceable></term>
<listitem><para>See <link linkend="irf">there</link></para></listitem>
</varlistentry>
<varlistentry>
<term><option>aim_solver</option></term>
<listitem><para>See <link linkend="aim_solver">there</link></para></listitem>
</varlistentry>
</variablelist>
......@@ -3186,7 +3195,7 @@ This problem is solved using a numerical optimizer.
</refsynopsisdiv>
<refsect1><title>Description</title>
<para>This function is an interface to the global sensitivity analysis (GSA) toolbox developed by the Joint Research Center (JRC) of the European Commission. The GSA toolbox needs to be downloaded separately from the JRC web site (<ulink url="http://eemc.jrc.ec.europa.eu/Software-DYNARE.htm">http://eemc.jrc.ec.europa.eu/Software-DYNARE.htm</ulink>).</para>
<para>This function is an interface to the global sensitivity analysis (GSA) toolbox developed by the Joint Research Center (JRC) of the European Commission. The GSA toolbox needs to be downloaded separately from the <ulink url="http://eemc.jrc.ec.europa.eu/Software-DYNARE.htm">JRC web site</ulink>.</para>
</refsect1>
<refsect1><title>Options</title>
......
......@@ -7,12 +7,12 @@ A subset of routines from Gary Anderson AIM package starting with SP... needed t
The path to the AIM directory,if exists, is added by dynare_config.m using addpath.
DR1 tries to invoke AIM if options_.useAIM == 1 is set and, if not check only, and if 1st order only:
if (options_.useAIM == 1) && (task == 0) && (options_.order == 1)
DR1 tries to invoke AIM if options_.aim_solver == 1 is set and, if not check only, and if 1st order only:
if (options_.aim_solver == 1) && (task == 0) && (options_.order == 1)
for start, options_.useAIM = 0 is set by default in global_initialization.m so that system uses mjdgges by default.
for start, options_.aim_solver = 0 is set by default in global_initialization.m so that system uses mjdgges by default.
If AIM is to be used, options_.useAIM = 1 needs to be set either in the model <>.mod file, before invoking, estimate and/or stoch_simul, or by issuing appropriate command for estimate and/or stoch_simul.
If AIM is to be used, options_.aim_solver = 1 needs to be set either in the model <>.mod file, before invoking, estimate and/or stoch_simul, or by issuing appropriate command for estimate and/or stoch_simul.
NOTE: in the current implementation, as of July 2008, handling of exceptions is rather fundamental and, in particular, when Blanchard and Kahn conditions are not met, only a large penalty value 1.0e+8 is being set.
......@@ -22,7 +22,7 @@ or
Error in AIM: aimcode=3 : Aim: too many big roots
especially close to the point of convergence.
However, if other exceptions occur and aimcode (see codes below) is higher than 5, the system resets options_.useAIM = 0 and tries to use mjdgges instead.
However, if other exceptions occur and aimcode (see codes below) is higher than 5, the system resets options_.aim_solver = 0 and tries to use mjdgges instead.
APPENDIX
......
......@@ -230,9 +230,11 @@ function [dr,info,M_,options_,oo_] = dr1(dr,task,M_,options_,oo_)
b(:,cols_b) = jacobia_(:,cols_j);
if M_.maximum_endo_lead == 0; % backward models
% If required, try Gary Anderson and G Moore AIM solver if not
% check only and if 1st order (added by GP July'08)
if (options_.useAIM == 1) && (task == 0) && (options_.order == 1)
% If required, use AIM solver if not check only
if (options_.aim_solver == 1) && (task == 0)
if options_.order > 1
error('Option "aim_solver" is incompatible with order >= 2')
end
try
[dr,aimcode]=dynAIMsolver1(jacobia_,M_,dr);
if aimcode ~=1
......@@ -241,13 +243,10 @@ function [dr,info,M_,options_,oo_] = dr1(dr,task,M_,options_,oo_)
return
end
catch
%warning('Problem with using AIM solver - Using Dynare solver instead');
disp('Problem with using AIM solver - Using Dynare solver instead');
options_.useAIM = 0; % and try mjdgges instead
disp(lasterror.message)
error('Problem with AIM solver - Try to remove the "aim_solver" option');
end
end % end if useAIM and...
%else use original Dynare solver
if ~((options_.useAIM == 1) && (task == 0) && (options_.order == 1))
else % use original Dynare solver
[k1,junk,k2] = find(kstate(:,4));
dr.ghx(:,k1) = -b\jacobia_(:,k2);
if M_.exo_nbr
......@@ -274,9 +273,11 @@ function [dr,info,M_,options_,oo_] = dr1(dr,task,M_,options_,oo_)
aa = jacobia_;
end
% If required, try Gary Anderson and G Moore AIM solver, that is, if not check
% only and if 1st order (added by GP July'08)
if (options_.useAIM == 1) && (task == 0) && (options_.order == 1)
% If required, use AIM solver if not check only
if (options_.aim_solver == 1) && (task == 0)
if options_.order > 1
error('Option "aim_solver" is incompatible with order >= 2')
end
try
[dr,aimcode]=dynAIMsolver1(aa,M_,dr);
......@@ -289,14 +290,6 @@ function [dr,info,M_,options_,oo_] = dr1(dr,task,M_,options_,oo_)
end
[A,B] =transition_matrix(dr);
dr.eigval = eig(A);
% if any(abs(dr.eigval) > options_.qz_criterium)
% temp = sort(abs(dr.eigval));
% nba = nnz(abs(dr.eigval) > options_.qz_criterium);
% temp = temp(nd-nba+1:nd)-1-options_.qz_criterium;
% info(1) = 3;
% info(2) = temp'*temp;
% return
% end
sdim = sum( abs(dr.eigval) < options_.qz_criterium );
nba = nd-sdim;
......@@ -314,13 +307,10 @@ function [dr,info,M_,options_,oo_] = dr1(dr,task,M_,options_,oo_)
return
end
catch
%warning('Problem with using AIM solver - Using Dynare solver instead');
disp('Problem with using AIM solver - Using Dynare solver instead');
options_.useAIM = 0; % and then try mjdgges instead
disp(lasterror.message)
error('Problem with AIM solver - Try to remove the "aim_solver" option')
end
end % end if useAIM and...
%else % use original Dynare solver
if ~((options_.useAIM == 1)&& (task == 0) && (options_.order == 1)) % || isempty(options_.useAIM)
else % use original Dynare solver
k1 = M_.lead_lag_incidence(find([1:klen] ~= M_.maximum_endo_lag+1),:);
a = aa(:,nonzeros(k1'));
b(:,cols_b) = aa(:,cols_j);
......@@ -456,15 +446,13 @@ function [dr,info,M_,options_,oo_] = dr1(dr,task,M_,options_,oo_)
dr.ghu = repmat(1./dr.ys(k1),1,size(dr.ghu,2)).*dr.ghu;
end
if ~((options_.useAIM == 1) && (options_.order == 1)) % if not use AIM ...
if options_.aim_solver ~= 1 && options_.use_qzdiv
%% Necessary when using Sims' routines for QZ
if options_.use_qzdiv
gx = real(gx);
hx = real(hx);
dr.ghx = real(dr.ghx);
dr.ghu = real(dr.ghu);
end
end % if not use AIM
gx = real(gx);
hx = real(hx);
dr.ghx = real(dr.ghx);
dr.ghu = real(dr.ghu);
end
%exogenous deterministic variables
if M_.exo_det_nbr > 0
......
......@@ -127,7 +127,7 @@ function global_initialization()
else
options_.use_qzdiv = 0;
end
options_.useAIM = 0; % i.e. by default do not use G.Anderson's AIM solver, use mjdgges instead
options_.aim_solver = 0; % i.e. by default do not use G.Anderson's AIM solver, use mjdgges instead
options_.use_k_order=0; % by default do not use k_order_perturbation but mjdgges
% Ramsey policy
......
......@@ -86,7 +86,7 @@ class ParsingDriver;
#define yylex driver.lexer->lex
%}
%token AR AUTOCORR
%token AIM_SOLVER AR AUTOCORR
%token BAYESIAN_IRF BETA_PDF BLOCK
%token BVAR_DENSITY BVAR_FORECAST
%token BVAR_PRIOR_DECAY BVAR_PRIOR_FLAT BVAR_PRIOR_LAMBDA
......@@ -716,6 +716,7 @@ stoch_simul_options : o_dr_algo
| o_qz_criterium
| o_print
| o_noprint
| o_aim_solver
;
symbol_list : symbol_list symbol
......@@ -1016,6 +1017,7 @@ estimation_options : o_datafile
| o_diffuse_filter
| o_plot_priors
| o_order
| o_aim_solver
;
prior_analysis : PRIOR_ANALYSIS '(' prior_posterior_options_list ')' ';'
......@@ -1087,11 +1089,11 @@ osr_params : OSR_PARAMS symbol_list ';' { driver.set_osr_params(); };
osr : OSR ';'
{ driver.run_osr(); }
| OSR '(' stoch_simul_options ')' ';'
| OSR '(' stoch_simul_options_list ')' ';'
{ driver.run_osr(); }
| OSR symbol_list ';'
{ driver.run_osr(); }
| OSR '(' stoch_simul_options ')' symbol_list ';'
| OSR '(' stoch_simul_options_list ')' symbol_list ';'
{driver.run_osr(); }
;
......@@ -1624,6 +1626,7 @@ o_noconstant : NOCONSTANT { driver.option_num("noconstant", "1"); };
o_mh_recover : MH_RECOVER { driver.option_num("mh_recover", "1"); };
o_diffuse_filter: DIFFUSE_FILTER {driver.option_num("diffuse_filter", "1"); };
o_plot_priors: PLOT_PRIORS {driver.option_num("plot_priors", "1"); };
o_aim_solver: AIM_SOLVER {driver.option_num("aim_solver", "1"); };
o_planner_discount : PLANNER_DISCOUNT EQUAL number { driver.option_num("planner_discount",$3); };
......
......@@ -221,6 +221,7 @@ int sigma_e = 0;
<DYNARE_STATEMENT>filename {return token::FILENAME;}
<DYNARE_STATEMENT>diffuse_filter {return token::DIFFUSE_FILTER;}
<DYNARE_STATEMENT>plot_priors {return token::PLOT_PRIORS;}
<DYNARE_STATEMENT>aim_solver {return token::AIM_SOLVER;}
<DYNARE_STATEMENT>freq {return token::FREQ;}
<DYNARE_STATEMENT>initial_year {return token::INITIAL_YEAR;}
......
......@@ -78,6 +78,5 @@ unit_root_vars P_obs Y_obs;
steady;
options_.useAIM = 0;
stoch_simul(order=1,irf=0);
......@@ -78,13 +78,13 @@ unit_root_vars P_obs Y_obs;
steady;
options_.useAIM = 1;
stoch_simul(order=1,irf=0);
stoch_simul(aim_solver, order=1, irf=0);
benchmark = load('fs2000_b1L1L_results');
threshold = 1e-8;
if max(max(abs(benchmark.oo_.dr.ghx-oo_.dr.ghx) > 1e-10));
exit('error in ghx');
elseif max(max(abs(benchmark.oo_.dr.ghu-oo_.dr.ghu) > 1e-10));
exit('error in ghy');
if max(max(abs(benchmark.oo_.dr.ghx-oo_.dr.ghx) > threshold));
error('error in ghx');
elseif max(max(abs(benchmark.oo_.dr.ghu-oo_.dr.ghu) > threshold));
error('error in ghy');
end;
......@@ -74,6 +74,5 @@ steady;
check;
options_.useAIM = 0;
stoch_simul(order=1,irf=0);
......@@ -74,13 +74,13 @@ steady;
check;
options_.useAIM = 1;
stoch_simul(order=1,irf=0);
stoch_simul(aim_solver, order=1,irf=0);
benchmark = load('fs2000x10L9_L_results');
threshold = 1e-8;
if max(max(abs(benchmark.oo_.dr.ghx-oo_.dr.ghx) > 1e-10));
exit('error in ghx');
elseif max(max(abs(benchmark.oo_.dr.ghu-oo_.dr.ghu) > 1e-10));
exit('error in ghy');
if max(max(abs(benchmark.oo_.dr.ghx-oo_.dr.ghx) > threshold));
error('error in ghx');
elseif max(max(abs(benchmark.oo_.dr.ghu-oo_.dr.ghu) > threshold));
error('error in ghy');
end;
......@@ -55,6 +55,5 @@ end;
steady;
options_.useAIM = 0;
stoch_simul(order=1,irf=0);
......@@ -55,13 +55,13 @@ end;
steady;
options_.useAIM = 1;
stoch_simul(order=1,irf=0);
stoch_simul(aim_solver, order=1,irf=0);
benchmark = load('fs2000x10_L9_L_results');
if max(max(abs(benchmark.oo_.dr.ghx-oo_.dr.ghx) > 1e-10));
exit('error in ghx');
elseif max(max(abs(benchmark.oo_.dr.ghu-oo_.dr.ghu) > 1e-10));
exit('error in ghy');
threshold = 1e-8;
if max(max(abs(benchmark.oo_.dr.ghx-oo_.dr.ghx) > threshold));
error('error in ghx');
elseif max(max(abs(benchmark.oo_.dr.ghu-oo_.dr.ghu) > threshold));
error('error in ghy');
end;
......@@ -43,6 +43,5 @@ var e_ys = 1.89;
var e_pies = 1.89;
end;
options_.useAIM = 0;
stoch_simul(order=1,irf=0);
......@@ -43,13 +43,13 @@ var e_ys = 1.89;
var e_pies = 1.89;
end;
options_.useAIM = 1;
stoch_simul(order=1,irf=0);
stoch_simul(aim_solver, order=1,irf=0);
benchmark = load('ls2003_2L0L_results');
threshold = 1e-8;
if max(max(abs(benchmark.oo_.dr.ghx-oo_.dr.ghx) > 1e-10));
exit('error in ghx');
elseif max(max(abs(benchmark.oo_.dr.ghu-oo_.dr.ghu) > 1e-10));
exit('error in ghy');
if max(max(abs(benchmark.oo_.dr.ghx-oo_.dr.ghx) > threshold));
error('error in ghx');
elseif max(max(abs(benchmark.oo_.dr.ghu-oo_.dr.ghu) > threshold));
error('error in ghy');
end;
......@@ -41,6 +41,5 @@ var e_ys = 1.89;
var e_pies = 1.89;
end;
options_.useAIM = 0;
stoch_simul(order=1,irf=0);
......@@ -41,14 +41,14 @@ var e_ys = 1.89;
var e_pies = 1.89;
end;
options_.useAIM = 1;
stoch_simul(order=1,irf=0);
stoch_simul(aim_solver, order=1,irf=0);
benchmark = load('ls2003_2L2L_results');
threshold = 1e-8;
if max(max(abs(benchmark.oo_.dr.ghx-oo_.dr.ghx) > 1e-10));
exit('error in ghx');
elseif max(max(abs(benchmark.oo_.dr.ghu-oo_.dr.ghu) > 1e-10));
exit('error in ghy');
if max(max(abs(benchmark.oo_.dr.ghx-oo_.dr.ghx) > threshold));
error('error in ghx');
elseif max(max(abs(benchmark.oo_.dr.ghu-oo_.dr.ghu) > threshold));
error('error in ghy');
end;
......@@ -16,14 +16,24 @@ OCTAVE_MODS = \
homotopy/homotopy2_test.mod \
homotopy/homotopy3_test.mod \
bvar_a_la_sims/bvar_standalone.mod \
bvar_a_la_sims/bvar_and_dsge.mod
bvar_a_la_sims/bvar_and_dsge.mod \
AIM/fs2000x10L9_L.mod \
AIM/fs2000x10L9_L_AIM.mod \
AIM/fs2000x10_L9_L.mod \
AIM/fs2000x10_L9_L_AIM.mod
MODS = $(OCTAVE_MODS) \
arima/mod1b.mod \
arima/mod1c.mod \
arima/mod2a.mod \
arima/mod2b.mod \
fs2000/fs2000a.mod
fs2000/fs2000a.mod \
AIM/fs2000_b1L1L.mod \
AIM/fs2000_b1L1L_AIM.mod \
AIM/ls2003_2L0L.mod \
AIM/ls2003_2L0L_AIM.mod \
AIM/ls2003_2L2L.mod \
AIM/ls2003_2L2L_AIM.mod
EXTRA_DIST = $(MODS) \
run_test_octave.m \
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment