Commit eed54fb0 authored by Ferhat's avatar Ferhat
Browse files

- Adds new algorithms to solve Lyapunov equations: Doubling algorithm and...

- Adds new algorithms to solve Lyapunov equations: Doubling algorithm and Square root solver. Their respective names are "doubling" and "square_root_solver".
- Adds the tolerance criteria for the iterative solvers (sylvester_fixed_point_tol, lyapunov_fixed_point_tol and lyapunov_doubling_tol)
- Updates the reference manual
parent 7d4a7b1b
......@@ -74,14 +74,14 @@ A copy of the license can be found at @uref{http://www.gnu.org/licenses/fdl.txt}
@titlepage
@title Dynare
@subtitle Reference Manual, version @value{VERSION}
@author Stéphane Adjemian
@author Stéphane Adjemian
@author Houtan Bastani
@author Michel Juillard
@author Junior Maih
@author Ferhat Mihoubi
@author George Perendia
@author Marco Ratto
@author Sébastien Villemot
@author Sébastien Villemot
@page
@vskip 0pt plus 1filll
@insertcopying
......@@ -291,11 +291,11 @@ The development of Dynare is mainly done at
@uref{http://www.cepremap.ens.fr, Cepremap} by a core team of
researchers who devote part of their time to software
development. Currently the development team of Dynare is composed of
Stéphane Adjemian (Université du Maine, Gains and Cepremap), Houtan
Bastani (Cepremap), Michel Juillard (Banque de France), Frédéric Karamé
(Université dvry, Epee and Cepremap), Junior Maih (Norges Bank),
Ferhat Mihoubi (Université d'Évry, Epee and Cepremap), George Perendia,
Johannes Pfeifer, Marco Ratto (JRC) and Sébastien Villemot (Cepremap and
Stéphane Adjemian (Université du Maine, Gains and Cepremap), Houtan
Bastani (Cepremap), Michel Juillard (Banque de France), Frédéric Karamé
(Université d'Évry, Epee and Cepremap), Junior Maih (Norges Bank),
Ferhat Mihoubi (Université d'Évry, Epee and Cepremap), George Perendia,
Johannes Pfeifer, Marco Ratto (JRC) and Sébastien Villemot (Cepremap and
Paris School of Economics). Increasingly, the developer base is
expanding, as tools developed by researchers outside of Cepremap are
integrated into Dynare. Financial support is provided by Cepremap,
......@@ -335,8 +335,8 @@ If you would like to refer to Dynare in a research article, the
recommended way is to cite the present manual, as follows:
@quotation
Stéphane Adjemian, Houtan Bastani, Michel Juillard, Ferhat Mihoubi,
George Perendia, Marco Ratto and Sébastien Villemot (2011), ``Dynare:
Stéphane Adjemian, Houtan Bastani, Michel Juillard, Ferhat Mihoubi,
George Perendia, Marco Ratto and Sébastien Villemot (2011), ``Dynare:
Reference Manual, Version 4,'' @i{Dynare Working Papers}, 1, CEPREMAP
@end quotation
......@@ -2890,7 +2890,7 @@ In a stochastic context, Dynare computes one or several simulations
corresponding to a random draw of the shocks. Dynare uses a Taylor
approximation, up to third order, of the expectation functions (see
@cite{Judd (1996)}, @cite{Collard and Juillard (2001a)}, @cite{Collard
and Juillard (2001b)}, and @cite{Schmitt-Grohé and Uríbe (2004)}). The
and Juillard (2001b)}, and @cite{Schmitt-Grohé and Uríbe (2004)}). The
details of the Dynare implementation of the first order solution are
given in @cite{Villemot (2011)}.
......@@ -3090,6 +3090,28 @@ that if @code{varobs} is not present or contains all endogenous
variables, then this is the full information case and this option has
no effect. More references can be found at
@uref{http://www.dynare.org/DynareWiki/PartialInformation}.
@item sylvester = OPTION
@anchor{sylvester}
Determines the algorithm used to solve the Sylvester equation for block decomposed model. Possible values for @code{OPTION} are:
@table @code
@item DEFAULT
Uses the default solver for Sylvester equations (@code{gensylv}) based on Ondra Kamenik algorithm (see @uref{http://www.dynare.org/documentation-and-support/dynarepp/sylvester.pdf/at_download/file} for more informations).
@item FIXED_POINT
Uses a fixed point algorithm to solve the Sylvester equation (@code{gensylv_fp}). This method is faster than the @code{DEFAULT} one for large scale models.
@end table
@noindent
Default value is @code{DEFAULT}
@item sylvester_fixed_point_tol = @var{DOUBLE}
@anchor{sylvester_fixed_point_tol}
It is the convergence criterion used in the fixed point sylvester solver. Its default value is 1e-12.
@end table
@outputhead
......@@ -4098,6 +4120,46 @@ with @ref{dsge_var}.
@item aim_solver
@xref{aim_solver}.
@item sylvester = OPTION
@xref{sylvester}.
@item sylvester_fixed_point_tol = @var{DOUBLE}
@xref{sylvester_fixed_point_tol}.
@item lyapunov = OPTION
@anchor{lyapunov}
Determines the algorithm used to solve the Laypunov equation to initialized the variance-covariance matrix of the Kalman filter using the steady-state value of state variables. Possible values for @code{OPTION} are:
@table @code
@item DEFAULT
Uses the default solver for Lyapunov equations (@code{lyapunov_symm} with method<=2) based on Bartels-Stewart algorithm.
@item FIXED_POINT
Uses a fixed point algorithm to solve the Lyapunov equation (@code{lyapunov_symm} with method=3). This method is faster than the @code{DEFAULT} one for large scale models, but it could require a large amount of iterations.
@item DOUBLING
Uses a doubling algorithm to solve the Lyapunov equation (@code{disclyap_fast}). This method is faster than the two previous one for large scale models.
@item SQUARE_ROOT_SOLVER
Uses a square-root solver for Lyapunov equations (@code{dlyapchol}). This method is fast for large scale models, but it requires the Control System Toolbox with MATLAB and the Control Package with OCTAVE.
@end table
@noindent
Default value is @code{DEFAULT}
@item lyapunov_fixed_point_tol = @var{DOUBLE}
@anchor{lyapunov_fixed_point_tol}
This is the convergence criterion used in the fixed point lyapunov solver. Its default value is 1e-10.
@item lyapunov_doubling_tol = @var{DOUBLE}
@anchor{lyapunov_doubling_tol}
This is the convergence criterion used in the doubling algorithm to solve the lyapunov equation. Its default value is 1e-16.
@end table
@customhead{Note}
......@@ -7068,9 +7130,9 @@ State Vector for Diffuse State Space Models,'' @i{Journal of Time
Series Analysis}, 24(1), 85--98.
@item
Laffargue, Jean-Pierre (1990): ``Résolution d'un modèle
macroéconomique avec anticipations rationnelles'', @i{Annales
d'Économie et Statistique}, 17, 97--119.
Laffargue, Jean-Pierre (1990): ``Résolution d'un modèle
macroéconomique avec anticipations rationnelles'', @i{Annales
d'Économie et Statistique}, 17, 97--119.
@item
Lubik, Thomas and Frank Schorfheide (2007): ``Do Central Banks Respond
......@@ -7100,7 +7162,7 @@ Schorfheide, Frank (2000): ``Loss Function-based evaluation of DSGE
models,'' @i{Journal of Applied Econometrics}, 15(6), 645--670.
@item
Schmitt-Grohé, Stephanie and Martin Uríbe (2004): ``Solving Dynamic
Schmitt-Grohé, Stephanie and Martin Uríbe (2004): ``Solving Dynamic
General Equilibrium Models Using a Second-Order Approximation to the
Policy Function,'' @i{Journal of Economic Dynamics and Control},
28(4), 755--775.
......@@ -7116,7 +7178,7 @@ Stochastic General Equilibrium Model of the Euro Area,'' @i{Journal of
the European Economic Association}, 1(5), 1123--1175.
@item
Villemot, Sébastien (2011): ``Solving rational expectations models at
Villemot, Sébastien (2011): ``Solving rational expectations models at
first order: what Dynare does,'' @i{Dynare Working Papers}, 2,
CEPREMAP
......
......@@ -568,7 +568,7 @@ for i = 1:Size;
vghx_other = - inv(kron(eye(size(D_,2)), A_) + kron(C_', B_)) * vec(D_);
ghx_other = reshape(vghx_other, size(D_,1), size(D_,2));
elseif options_.sylvester_fp == 1
ghx_other = gensylv_fp(A_, B_, C_, D_, i);
ghx_other = gensylv_fp(A_, B_, C_, D_, i, options_.sylvester_fixed_point_tol);
else
[err, ghx_other] = gensylv(1, A_, B_, C_, -D_);
end;
......
......@@ -367,7 +367,11 @@ switch DynareOptions.lik_init
kalman_algo = 1;
end
if DynareOptions.lyapunov_fp == 1
Pstar = lyapunov_symm(T,Q,DynareOptions.qz_criterium,DynareOptions.lyapunov_complex_threshold, 4, R);
Pstar = lyapunov_symm(T,Q,DynareOptions.lyapunov_fixed_point_tol,DynareOptions.lyapunov_complex_threshold, 3, R);
elseif DynareOptions.lyapunov_db == 1
Pstar = disclyap_fast(T,R*Q*R',DynareOptions.lyapunov_doubling_tol);
elseif DynareOptions.lyapunov_srs == 1
Pstar = lyapunov_symm(T,Q,DynareOptions.lyapunov_fixed_point_tol,DynareOptions.lyapunov_complex_threshold, 4, R);
else
Pstar = lyapunov_symm(T,R*Q*R',DynareOptions.qz_criterium,DynareOptions.lyapunov_complex_threshold);
end;
......
function X = gensylv_fp(A, B, C, D, block)
function X = gensylv_fp(A, B, C, D, block, tol)
% function X = gensylv_fp(A, B, C, D)
% Solve the Sylvester equation:
% A * X + B * X * C + D = 0
......@@ -8,6 +8,7 @@ function X = gensylv_fp(A, B, C, D, block)
% C
% D
% block : block number (for storage purpose)
% tol : convergence criteria
% OUTPUTS
% X solution
%
......@@ -39,9 +40,8 @@ function X = gensylv_fp(A, B, C, D, block)
% You should have received a copy of the GNU General Public License
% along with Dynare. If not, see <http://www.gnu.org/licenses/>.
%tol = 1e-07;
%tol = 1e-13;
tol = 1e-12;
%tol = 1e-12;
evol = 100;
A1 = inv(A);
eval(['persistent hxo_' int2str(block) ';']);
......
......@@ -413,13 +413,22 @@ options_.use_dll = 0;
% model evaluated using bytecode.dll
options_.bytecode = 0;
% use a fixed point method to solve Sylvester equation (for large scale
% models)
% if equal to 1 use a fixed point method to solve Sylvester equation (for large scale models)
options_.sylvester_fp = 0;
% use a fixed point method to solve Lyapunov equation (for large scale
% models)
% convergence criteria to solve iteratively a sylvester equations
options_.sylvester_fixed_point_tol = 1e-12;
% if 1 use a fixed point method to solve Lyapunov equation (for large scale models)
options_.lyapunov_fp = 0;
% if 1 use a doubling algorithm to solve Lyapunov equation (for large scale models)
options_.lyapunov_db = 0;
% if 1 use a squre root solver to solve Lyapunov equation (for large scale models)
options_.lyapunov_srs = 0;
% convergence criterion for iteratives methods to solve lyapunov equations
options_.lyapunov_fixed_point_tol = 1e-10;
options_.lyapunov_doubling_tol = 1e-16;
% dates for historical time series
options_.initial_date.freq = 1;
......
function [x,u] = lyapunov_symm(a,b,qz_criterium,lyapunov_complex_threshold,method, R)
function [x,u] = lyapunov_symm(a,b,third_argument,lyapunov_complex_threshold,method, R)
% Solves the Lyapunov equation x-a*x*a' = b, for b and x symmetric matrices.
% If a has some unit roots, the function computes only the solution of the stable subsystem.
%
% INPUTS:
% a [double] n*n matrix.
% b [double] n*n matrix.
% qz_criterium [double] scalar, unit root threshold for eigenvalues in a.
% third_argument [double] scalar, if method <= 2 :
% qz_criterium = third_argument unit root threshold for eigenvalues in a,
% elseif method = 3 :
% tol =third_argument the convergence criteria for fixed_point algorithm.
% lyapunov_complex_threshold [double] scalar, complex block threshold for the upper triangular matrix T.
% method [integer] Scalar, if method=0 [default] then U, T, n and k are not persistent.
% method=1 then U, T, n and k are declared as persistent
......@@ -48,10 +51,10 @@ if method == 3
if ~isempty(method1)
method = method1;
end;
tol = third_argument;
fprintf(' [methode=%d] ',method);
if method == 3
%tol = 1e-8;
tol = 1e-10;
%tol = 1e-10;
it_fp = 0;
evol = 100;
if isempty(X)
......@@ -81,13 +84,15 @@ if method == 3
end;
end;
elseif method == 4
% works only with Matlab System Control toolbox
% works only with Matlab System Control toolbox or octave the control package,
% the presence of the toolbox or package has to be tested
chol_b = R*chol(b,'lower');
Rx = dlyapchol(a,chol_b);
x = Rx' * Rx;
return;
end;
qz_criterium = third_argument;
if method
persistent U T k n
else
......
......@@ -134,7 +134,7 @@ exo_names=M_.exo_names(M_.exo_names_orig_ord,:);
DPDR=DD*PP*DD'+RR;
I_DPDR=inv(DPDR);
PDIDPDRD=PP*DD'*I_DPDR*DD;
MSIG=disclyap_fast(CCCC, CCCC*PDIDPDRD*PP*CCCC');
MSIG=disclyap_fast(CCCC, CCCC*PDIDPDRD*PP*CCCC', options_.lyapunov_doubling_tol);
COV_P=[ PP, PP; PP, PP+MSIG]; % P0
......
function X=disclyap_fast(G,V,ch)
function X=disclyap_fast(G,V,tol,ch)
% function X=disclyap_fast(G,V,ch)
%
% Solve the discrete Lyapunov Equation
......@@ -28,14 +28,14 @@ function X=disclyap_fast(G,V,ch)
% You should have received a copy of the GNU General Public License
% along with Dynare. If not, see <http://www.gnu.org/licenses/>.
if nargin == 2 || isempty( ch ) == 1
if nargin <= 3 || isempty( ch ) == 1
flag_ch = 0;
else
flag_ch = 1;
end
s=size(G,1);
tol = 1e-16;
%tol = 1e-16;
P0=V;
A0=G;
......
......@@ -96,7 +96,7 @@ class ParsingDriver;
%token BVAR_PRIOR_MU BVAR_PRIOR_OMEGA BVAR_PRIOR_TAU BVAR_PRIOR_TRAIN
%token BVAR_REPLIC BYTECODE
%token CHANGE_TYPE CHECK CONDITIONAL_FORECAST CONDITIONAL_FORECAST_PATHS CONF_SIG CONSTANT CONTROLLED_VAREXO CORR COVAR CUTOFF
%token DATAFILE FILE DR_ALGO DROP DSAMPLE DYNASAVE DYNATYPE CALIBRATION
%token DATAFILE FILE DOUBLING DR_ALGO DROP DSAMPLE DYNASAVE DYNATYPE CALIBRATION
%token END ENDVAL EQUAL ESTIMATION ESTIMATED_PARAMS ESTIMATED_PARAMS_BOUNDS ESTIMATED_PARAMS_INIT
%token FILENAME FILTER_STEP_AHEAD FILTERED_VARS FIRST_OBS LAST_OBS SET_TIME
%token <string_val> FLOAT_NUMBER
......@@ -110,7 +110,7 @@ class ParsingDriver;
%token INV_GAMMA_PDF INV_GAMMA1_PDF INV_GAMMA2_PDF IRF IRF_SHOCKS
%token KALMAN_ALGO KALMAN_TOL SUBSAMPLES OPTIONS
%token LABELS LAPLACE LIK_ALGO LIK_INIT LINEAR LOAD_IDENT_FILES LOAD_MH_FILE LOAD_PARAMS_AND_STEADY_STATE LOGLINEAR LYAPUNOV
%token MARKOWITZ MARGINAL_DENSITY MAX MAXIT
%token LYAPUNOV_FIXED_POINT_TOL LYAPUNOV_DOUBLING_TOL LYAPUNOV_SQUARE_ROOT_SOLVER_TOL MARKOWITZ MARGINAL_DENSITY MAX MAXIT
%token MFS MH_DROP MH_INIT_SCALE MH_JSCALE MH_MODE MH_NBLOCKS MH_REPLIC MH_RECOVER MIN MINIMAL_SOLVING_PERIODS
%token MODE_CHECK MODE_COMPUTE MODE_FILE MODEL MODEL_COMPARISON MODEL_INFO MSHOCKS ABS SIGN
%token MODIFIEDHARMONICMEAN MOMENTS_VARENDO DIFFUSE_FILTER SUB_DRAWS
......@@ -123,8 +123,8 @@ class ParsingDriver;
%token <string_val> QUOTED_STRING
%token QZ_CRITERIUM FULL DSGE_VAR DSGE_VARLAG DSGE_PRIOR_WEIGHT TRUNCATE
%token RELATIVE_IRF REPLIC RPLOT SAVE_PARAMS_AND_STEADY_STATE PARAMETER_UNCERTAINTY
%token SHOCKS SHOCK_DECOMPOSITION SIGMA_E SIMUL SIMUL_ALGO SIMUL_SEED SMOOTHER STACK_SOLVE_ALGO STEADY_STATE_MODEL SOLVE_ALGO
%token STDERR STEADY STOCH_SIMUL SYLVESTER REGIMES REGIME
%token SHOCKS SHOCK_DECOMPOSITION SIGMA_E SIMUL SIMUL_ALGO SIMUL_SEED SMOOTHER SQUARE_ROOT_SOLVER STACK_SOLVE_ALGO STEADY_STATE_MODEL SOLVE_ALGO
%token STDERR STEADY STOCH_SIMUL SYLVESTER SYLVESTER_FIXED_POINT_TOL REGIMES REGIME
%token TEX RAMSEY_POLICY PLANNER_DISCOUNT DISCRETIONARY_POLICY
%token <string_val> TEX_NAME
%token UNIFORM_PDF UNIT_ROOT_VARS USE_DLL USEAUTOCORR GSA_SAMPLE_FILE
......@@ -930,6 +930,7 @@ stoch_simul_options : o_dr_algo
| o_k_order_solver
| o_pruning
| o_sylvester
| o_sylvester_fixed_point_tol
;
symbol_list : symbol_list symbol
......@@ -1491,7 +1492,10 @@ estimation_options : o_datafile
| o_irf_shocks
| o_sub_draws
| o_sylvester
| o_sylvester_fixed_point_tol
| o_lyapunov
| o_lyapunov_fixed_point_tol
| o_lyapunov_doubling_tol
;
list_optim_option : QUOTED_STRING COMMA QUOTED_STRING
......@@ -2232,8 +2236,13 @@ o_sub_draws: SUB_DRAWS EQUAL INT_NUMBER {driver.option_num("sub_draws",$3);};
o_planner_discount : PLANNER_DISCOUNT EQUAL expression { driver.declare_optimal_policy_discount_factor_parameter($3); };
o_sylvester : SYLVESTER EQUAL FIXED_POINT {driver.option_num("sylvester_fp", "1"); }
| SYLVESTER EQUAL DEFAULT {driver.option_num("sylvester_fp", "0"); };
o_sylvester_fixed_point_tol : SYLVESTER_FIXED_POINT_TOL EQUAL non_negative_number {driver.option_num("sylvester_fixed_point_tol",$3);};
o_lyapunov : LYAPUNOV EQUAL FIXED_POINT {driver.option_num("lyapunov_fp", "1"); }
| LYAPUNOV EQUAL DEFAULT {driver.option_num("lyapunov_fp", "0"); };
| LYAPUNOV EQUAL DOUBLING {driver.option_num("lyapunov_db", "1"); }
| LYAPUNOV EQUAL SQUARE_ROOT_SOLVER {driver.option_num("lyapunov_srs", "1"); }
| LYAPUNOV EQUAL DEFAULT {driver.option_num("lyapunov_fp", "0");driver.option_num("lyapunov_db", "0"); driver.option_num("lyapunov_srs", "0");};
o_lyapunov_fixed_point_tol : LYAPUNOV_FIXED_POINT_TOL EQUAL non_negative_number {driver.option_num("lyapunov_fixed_point_tol",$3);};
o_lyapunov_doubling_tol : LYAPUNOV_DOUBLING_TOL EQUAL non_negative_number {driver.option_num("lyapunov_doubling_tol",$3);};
o_bvar_prior_tau : BVAR_PRIOR_TAU EQUAL signed_number { driver.option_num("bvar_prior_tau", $3); };
o_bvar_prior_decay : BVAR_PRIOR_DECAY EQUAL non_negative_number { driver.option_num("bvar_prior_decay", $3); };
......
......@@ -293,6 +293,8 @@ string eofbuff;
<DYNARE_STATEMENT>nstates {return token::NSTATES;}
<DYNARE_STATEMENT>indxscalesstates {return token::INDXSCALESSTATES;}
<DYNARE_STATEMENT>fixed_point {return token::FIXED_POINT;}
<DYNARE_STATEMENT>doubling {return token::DOUBLING;}
<DYNARE_STATEMENT>square_root_solver {return token::SQUARE_ROOT_SOLVER;}
<DYNARE_STATEMENT>default {return token::DEFAULT;}
<DYNARE_STATEMENT>alpha {
yylval->string_val = new string(yytext);
......@@ -491,6 +493,9 @@ string eofbuff;
<DYNARE_STATEMENT>order {return token::ORDER;}
<DYNARE_STATEMENT>sylvester {return token::SYLVESTER;}
<DYNARE_STATEMENT>lyapunov {return token::LYAPUNOV;}
<DYNARE_STATEMENT>sylvester_fixed_point_tol {return token::SYLVESTER_FIXED_POINT_TOL;}
<DYNARE_STATEMENT>lyapunov_fixed_point_tol {return token::LYAPUNOV_FIXED_POINT_TOL;}
<DYNARE_STATEMENT>lyapunov_doubling_tol {return token::LYAPUNOV_DOUBLING_TOL;}
<DYNARE_STATEMENT>replic {return token::REPLIC;}
<DYNARE_STATEMENT>ar {return token::AR;}
<DYNARE_STATEMENT>nofunctions {return token::NOFUNCTIONS;}
......
Supports Markdown
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