diff --git a/matlab/dr_block.m b/matlab/dr_block.m
index 7f7df97cc83b1bd3bf2f92050d6cf3fab2c7c920..e7d7d0398f8817467b7c0818e5dbc634a3dfa57d 100644
--- a/matlab/dr_block.m
+++ b/matlab/dr_block.m
@@ -567,7 +567,9 @@ for i = 1:Size;
if block_type == 5
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));
- else
+ elseif options_.sylvester_fp == 1
+ ghx_other = gensylv_fp(A_, B_, C_, D_, i);
+ else
[err, ghx_other] = gensylv(1, A_, B_, C_, -D_);
end;
if options_.aim_solver ~= 1 && options_.use_qzdiv
diff --git a/matlab/dsge_likelihood.m b/matlab/dsge_likelihood.m
index fbb2d1eab32148bd45b1813bcbe388a62591ca74..282a8abddf715eedda08f6f28e378379efd08c2e 100644
--- a/matlab/dsge_likelihood.m
+++ b/matlab/dsge_likelihood.m
@@ -366,7 +366,11 @@ switch DynareOptions.lik_init
% Use standard kalman filter except if the univariate filter is explicitely choosen.
kalman_algo = 1;
end
- Pstar = lyapunov_symm(T,R*Q*R',DynareOptions.qz_criterium,DynareOptions.lyapunov_complex_threshold);
+ if DynareOptions.lyapunov_fp == 1
+ Pstar = lyapunov_symm(T,Q,DynareOptions.qz_criterium,DynareOptions.lyapunov_complex_threshold, 4, R);
+ else
+ Pstar = lyapunov_symm(T,R*Q*R',DynareOptions.qz_criterium,DynareOptions.lyapunov_complex_threshold);
+ end;
Pinf = [];
a = zeros(mm,1);
Zflag = 0;
diff --git a/matlab/gensylv_fp.m b/matlab/gensylv_fp.m
new file mode 100644
index 0000000000000000000000000000000000000000..f56a65a9183b7939b5ef45d0ef67a529e3dc71bb
--- /dev/null
+++ b/matlab/gensylv_fp.m
@@ -0,0 +1,73 @@
+function X = gensylv_fp(A, B, C, D, block)
+% function X = gensylv_fp(A, B, C, D)
+% Solve the Sylvester equation:
+% A * X + B * X * C + D = 0
+% INPUTS
+% A
+% B
+% C
+% D
+% block : block number (for storage purpose)
+% OUTPUTS
+% X solution
+%
+% ALGORITHM
+% fixed point method
+% MARLLINY MONSALVE (2008): "Block linear method for large scale
+% Sylvester equations", Computational & Applied Mathematics, Vol 27, n�1,
+% p47-59
+% ||A^-1||.||B||.||C|| < 1 is a suffisant condition:
+% - to get a unique solution for the Sylvester equation
+% - to get a convergent fixed-point algorithm
+%
+% SPECIAL REQUIREMENTS
+% none.
+% Copyright (C) 1996-2010 Dynare Team
+%
+% This file is part of Dynare.
+%
+% Dynare is free software: you can redistribute it and/or modify
+% it under the terms of the GNU General Public License as published by
+% the Free Software Foundation, either version 3 of the License, or
+% (at your option) any later version.
+%
+% Dynare is distributed in the hope that it will be useful,
+% but WITHOUT ANY WARRANTY; without even the implied warranty of
+% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+% GNU General Public License for more details.
+%
+% 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;
+evol = 100;
+A1 = inv(A);
+eval(['persistent hxo_' int2str(block) ';']);
+hxo = eval(['hxo_' int2str(block) ';']);
+if isempty(hxo)
+ X = zeros(size(B, 2), size(C, 1));
+else
+ X = hxo;
+end;
+it_fp = 0;
+maxit_fp = 1000;
+Z = - (B * X * C + D);
+while it_fp < maxit_fp && evol > tol;
+ %X_old = X;
+ %X = - A1 * ( B * X * C + D);
+ %evol = max(max(abs(X - X_old)));
+ X = A1 * Z;
+ Z_old = Z;
+ Z = - (B * X * C + D);
+ evol = max(sum(abs(Z - Z_old))); %norm_1
+ %evol = max(sum(abs(Z - Z_old)')); %norm_inf
+ it_fp = it_fp + 1;
+end;
+%fprintf('sylvester it_fp=%d evol=%g | ',it_fp,evol);
+if evol < tol
+ eval(['hxo_' int2str(block) ' = X;']);
+else
+ error(['convergence not achieved in fixed point solution of Sylvester equation after ' int2str(it_fp) ' iterations']);
+end;
\ No newline at end of file
diff --git a/matlab/global_initialization.m b/matlab/global_initialization.m
index 2ee53d90f3d7a75d351085da5a8509ce5c1d7be0..0e5dd3eefc9c00d3585ca5ca60331254b2208de6 100644
--- a/matlab/global_initialization.m
+++ b/matlab/global_initialization.m
@@ -382,6 +382,14 @@ 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)
+options_.sylvester_fp = 0;
+
+% use a fixed point method to solve Lyapunov equation (for large scale
+% models)
+options_.lyapunov_fp = 0;
+
% dates for historical time series
options_.initial_date.freq = 1;
options_.initial_date.period = 1;
diff --git a/matlab/lyapunov_symm.m b/matlab/lyapunov_symm.m
index d32f644394bb8f537abe97ad2b64768e4669c169..afa763c82f43036ff0332c643a79609bf1d7438d 100644
--- a/matlab/lyapunov_symm.m
+++ b/matlab/lyapunov_symm.m
@@ -1,4 +1,4 @@
-function [x,u] = lyapunov_symm(a,b,qz_criterium,lyapunov_complex_threshold,method)
+function [x,u] = lyapunov_symm(a,b,qz_criterium,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.
%
@@ -12,6 +12,7 @@ function [x,u] = lyapunov_symm(a,b,qz_criterium,lyapunov_complex_threshold,metho
% variables and the schur decomposition is triggered.
% method=2 then U, T, n and k are declared as persistent
% variables and the schur decomposition is not performed.
+% method=3 fixed point method
% OUTPUTS
% x: [double] m*m solution matrix of the lyapunov equation, where m is the dimension of the stable subsystem.
% u: [double] Schur vectors associated with unit roots
@@ -38,10 +39,55 @@ function [x,u] = lyapunov_symm(a,b,qz_criterium,lyapunov_complex_threshold,metho
%
% 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<5
method = 0;
end
+
+if method == 3
+ persistent X method1;
+ if ~isempty(method1)
+ method = method1;
+ end;
+ fprintf(' [methode=%d] ',method);
+ if method == 3
+ %tol = 1e-8;
+ tol = 1e-10;
+ it_fp = 0;
+ evol = 100;
+ if isempty(X)
+ X = b;
+ max_it_fp = 2000;
+ else
+ max_it_fp = 300;
+ end;
+ at = a';
+ %fixed point iterations
+ while evol > tol && it_fp < max_it_fp;
+ X_old = X;
+ X = a * X * at + b;
+ evol = max(sum(abs(X - X_old))); %norm_1
+ %evol = max(sum(abs(X - X_old)')); %norm_inf
+ it_fp = it_fp + 1;
+ end;
+ fprintf('lyapunov it_fp=%d evol=%g\n',it_fp,evol);
+ if it_fp >= max_it_fp
+ disp(['convergence not achieved in solution of Lyapunov equation after ' int2str(it_fp) ' iterations, switching method from 3 to 0']);
+ method1 = 0;
+ method = 0;
+ else
+ method1 = 3;
+ x = X;
+ return;
+ end;
+ end;
+elseif method == 4
+ % works only with Matlab System Control toolbox
+ chol_b = R*chol(b,'lower');
+ Rx = dlyapchol(a,chol_b);
+ x = Rx' * Rx;
+ return;
+end;
+
if method
persistent U T k n
else
@@ -113,4 +159,4 @@ if i == 1
x(1,1) = (B(1,1)+c)/(1-T(1,1)*T(1,1));
end
x = U(:,k+1:end)*x*U(:,k+1:end)';
-u = U(:,1:k);
\ No newline at end of file
+u = U(:,1:k);
diff --git a/preprocessor/DynareBison.yy b/preprocessor/DynareBison.yy
index 6a5a2055b873cebe307af96a71ee7a4932a564c0..40857352d314936f48c81608ed7f7e783bc19946 100644
--- a/preprocessor/DynareBison.yy
+++ b/preprocessor/DynareBison.yy
@@ -100,6 +100,7 @@ class ParsingDriver;
%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
+%token DEFAULT FIXED_POINT
%token FORECAST K_ORDER_SOLVER INSTRUMENTS PRIOR SHIFT MEAN STDEV VARIANCE MODE INTERVAL SHAPE DOMAINN
%token GAMMA_PDF GRAPH CONDITIONAL_VARIANCE_DECOMPOSITION NOCHECK STD
%token HISTVAL HOMOTOPY_SETUP HOMOTOPY_MODE HOMOTOPY_STEPS HP_FILTER HP_NGRID
@@ -108,7 +109,7 @@ class ParsingDriver;
%token <string_val> DATE_NUMBER
%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
+%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 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
@@ -123,7 +124,7 @@ class ParsingDriver;
%token QZ_CRITERIUM FULL DSGE_VAR DSGE_VARLAG DSGE_PRIOR_WEIGHT
%token RELATIVE_IRF REPLIC RPLOT SAVE_PARAMS_AND_STEADY_STATE
%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
+%token STDERR STEADY STOCH_SIMUL SYLVESTER
%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
@@ -694,8 +695,8 @@ svar_identification_elem : EXCLUSION LAG INT_NUMBER ';' svar_equation_list
{ driver.add_constants_exclusion(); }
| RESTRICTION EQUATION INT_NUMBER COMMA
{ driver.add_restriction_equation_nbr($3);}
- restriction_expression EQUAL
- {driver.add_restriction_equal();}
+ restriction_expression EQUAL
+ {driver.add_restriction_equal();}
restriction_expression ';'
| UPPER_CHOLESKY ';'
{ driver.add_upper_cholesky(); }
@@ -925,6 +926,7 @@ stoch_simul_options : o_dr_algo
| o_conditional_variance_decomposition
| o_k_order_solver
| o_pruning
+ | o_sylvester
;
symbol_list : symbol_list symbol
@@ -1339,6 +1341,8 @@ estimation_options : o_datafile
| o_cova_compute
| o_irf_shocks
| o_sub_draws
+ | o_sylvester
+ | o_lyapunov
;
list_optim_option : QUOTED_STRING COMMA QUOTED_STRING
@@ -2072,6 +2076,10 @@ o_aim_solver: AIM_SOLVER {driver.option_num("aim_solver", "1"); };
o_partial_information : PARTIAL_INFORMATION {driver.option_num("partial_information", "1"); };
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_lyapunov : LYAPUNOV EQUAL FIXED_POINT {driver.option_num("lyapunov_fp", "1"); }
+ | LYAPUNOV EQUAL DEFAULT {driver.option_num("lyapunov_fp", "0"); };
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); };
diff --git a/preprocessor/DynareFlex.ll b/preprocessor/DynareFlex.ll
index 2d91d096e1330661937395ceac979b689a4cf6e0..a5010b5bdc198b1b411b8b09f1c68e324504268a 100644
--- a/preprocessor/DynareFlex.ll
+++ b/preprocessor/DynareFlex.ll
@@ -290,6 +290,8 @@ string eofbuff;
<DYNARE_STATEMENT>dummy_obs {return token::DUMMY_OBS;}
<DYNARE_STATEMENT>nstates {return token::NSTATES;}
<DYNARE_STATEMENT>indxscalesstates {return token::INDXSCALESSTATES;}
+<DYNARE_STATEMENT>fixed_point {return token::FIXED_POINT;}
+<DYNARE_STATEMENT>default {return token::DEFAULT;}
<DYNARE_STATEMENT>alpha {
yylval->string_val = new string(yytext);
return token::ALPHA;
@@ -481,6 +483,8 @@ string eofbuff;
<DYNARE_STATEMENT>stack_solve_algo {return token::STACK_SOLVE_ALGO;}
<DYNARE_STATEMENT>drop {return token::DROP;}
<DYNARE_STATEMENT>order {return token::ORDER;}
+<DYNARE_STATEMENT>sylvester {return token::SYLVESTER;}
+<DYNARE_STATEMENT>lyapunov {return token::LYAPUNOV;}
<DYNARE_STATEMENT>replic {return token::REPLIC;}
<DYNARE_STATEMENT>ar {return token::AR;}
<DYNARE_STATEMENT>nofunctions {return token::NOFUNCTIONS;}