- 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

doc/dynare.texi

@@ -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
+Stéphane Adjemian (Université du Maine, Gains and Cepremap), Houtan
-Bastani (Cepremap), Michel Juillard (Banque de France), Frédéric Karamé
+Bastani (Cepremap), Michel Juillard (Banque de France), Frédéric Karamé
-(Université d'Évry, Epee and Cepremap), Junior Maih (Norges Bank),
+(Université d'Évry, Epee and Cepremap), Junior Maih (Norges Bank),
-Ferhat Mihoubi (Université d'Évry, Epee and Cepremap), George Perendia,
+Ferhat Mihoubi (Université d'Évry, Epee and Cepremap), George Perendia,
-Johannes Pfeifer, Marco Ratto (JRC) and Sébastien Villemot (Cepremap and
+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,
+Stéphane Adjemian, Houtan Bastani, Michel Juillard, Ferhat Mihoubi,
-George Perendia, Marco Ratto and Sébastien Villemot (2011), ``Dynare:
+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
+Laffargue, Jean-Pierre (1990): ``Résolution d'un modèle
-macroéconomique avec anticipations rationnelles'', @i{Annales
+macroéconomique avec anticipations rationnelles'', @i{Annales
-d'Économie et Statistique}, 17, 97--119.
+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
 

matlab/dr_block.m

@@ -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;

matlab/dsge_likelihood.m

@@ -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;

matlab/gensylv_fp.m

@@ -1,4 +1,4 @@
-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) ';']);

matlab/global_initialization.m

@@ -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
+% if equal to 1 use a fixed point method to solve Sylvester equation (for large scale models)
-% models)
 options_.sylvester_fp = 0;
 
-% use a fixed point method to solve Lyapunov equation (for large scale
+% convergence criteria to solve iteratively a sylvester equations
-% models)
+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;

matlab/lyapunov_symm.m

@@ -1,11 +1,14 @@
-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

matlab/partial_information/PCL_Part_info_moments.m

@@ -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
 

matlab/partial_information/disclyap_fast.m

@@ -1,4 +1,4 @@
-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; 

preprocessor/DynareBison.yy

@@ -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 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 REGIMES REGIME
+%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); };

preprocessor/DynareFlex.ll

@@ -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;}