From 4488357f599d116f8d3a4521abf9d04452bde006 Mon Sep 17 00:00:00 2001 From: Ferhat Mihoubi Date: Sun, 1 Jul 2012 15:15:52 +0200 Subject: [PATCH] Adds the cycle reduction algorithm to solve the polynomial equation for retrieving the coefficients associated to the endogenous variables in the decision rule. --- doc/dynare.texi | 23 +++++ matlab/cycle_reduction.m | 69 ++++++++++++++ matlab/dr_block.m | 168 ++++++++++++++++++--------------- matlab/global_initialization.m | 9 +- preprocessor/DynareBison.yy | 11 ++- preprocessor/DynareFlex.ll | 3 + 6 files changed, 203 insertions(+), 80 deletions(-) create mode 100644 matlab/cycle_reduction.m diff --git a/doc/dynare.texi b/doc/dynare.texi index aa5eb015e..693373189 100644 --- a/doc/dynare.texi +++ b/doc/dynare.texi @@ -3168,6 +3168,29 @@ Default value is @code{default} @anchor{sylvester_fixed_point_tol} It is the convergence criterion used in the fixed point sylvester solver. Its default value is 1e-12. +@item dr = @var{OPTION} +@anchor{dr} +Determines the method used to compute the decision rule. Possible values for @code{@var{OPTION}} are: + +@table @code + +@item default +Uses the default method to compute the decision rule based on the generailzed schur decompostion +(see @uref{http://www.dynare.org/wp-repo/dynarewp002.pdf,the Dynare Website} for more information). + +@item cycle_reduction +Uses the cycle reduction algorithm to solve the polynomial equation for retrieving the coefficients +associated to the endogenous variables in the decision rule. This method is faster than the @code{default} one for large scale models. + +@end table + +@noindent +Default value is @code{default} + +@item dr_cycle_reduction_tol = @var{DOUBLE} +@anchor{dr_cycle_reduction_tol} +It is the convergence criterion used in the cycle reduction algorithm. Its default value is 1e-7. + @end table @outputhead diff --git a/matlab/cycle_reduction.m b/matlab/cycle_reduction.m new file mode 100644 index 000000000..355d5956e --- /dev/null +++ b/matlab/cycle_reduction.m @@ -0,0 +1,69 @@ +function [X, info] = cycle_reduction(A0, A1, A2, cvg_tol, ch) +% function [X, info] = cycle_reduction(A0,A1,A2,A3, cvg_tolch) +% +% Solves Polynomial Equation: +% A0 + A1 * X + A2 * X² = 0 +% Using Cyclic Reduction algorithm +% - D.A. Bini, G. Latouche, B. Meini (2002), "Solving matrix polynomial equations arising in +% queueing problems", Linear Algebra and its Applications 340 (2002) 225–244 +% - D.A. Bini, B. Meini, On the solution of a nonlinear matrix equation arising in queueing problems, +% SIAM J. Matrix Anal. Appl. 17 (1996) 906–926. +% ================================================================= + +% Copyright (C) 2006-2012 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 . + + max_it = 300; + it = 0; + info = 0; + crit = 1+cvg_tol; + A_0 = A1; + A0_0 = A0; + A1_0 = A1; + A2_0 = A2; + while crit > cvg_tol && it < max_it; + i_A1_0 = inv(A1_0); + A2_0_i_A1_0 = A2_0 * i_A1_0; + A0_0_i_A1_0 = A0_0 * i_A1_0; + A1_INC = A2_0_i_A1_0 * A0_0; + A_1 = A_0 - A1_INC; + A0_1 = - A0_0_i_A1_0 * A0_0; + A1_1 = A1_0 - A0_0_i_A1_0 * A2_0 - A1_INC; + A2_1 = - A2_0_i_A1_0 * A2_0; + + + crit = sum(sum(abs(A0_0))); + + A_0 = A_1; + A0_0 = A0_1; + A1_0 = A1_1; + A2_0 = A2_1; + it = it + 1; + %disp(['it=' int2str(it) ' crit = ' num2str(crit)]); + end; + if it==max_it + disp(['convergence not achieved after ' int2str(it) ' iterations']); + info = 1; + end + X = - inv(A_0) * A0; + if (nargin == 5 && ~isempty( ch ) == 1 ) + %check the solution + res = A0 + A1 * X + A2 * X * X; + if (sum(sum(abs(res))) > cvg_tol) + disp(['the norm residual of the residu ' num2str(res) ' compare to the tolerance criterion ' num2str(cvg_tol)]); + end; + end; \ No newline at end of file diff --git a/matlab/dr_block.m b/matlab/dr_block.m index 4e6f86c12..7bcdadedc 100644 --- a/matlab/dr_block.m +++ b/matlab/dr_block.m @@ -422,88 +422,102 @@ for i = 1:Size; row_indx = n_static+1:n; - D = [[aa(row_indx,index_0m) zeros(n_dynamic,n_both) aa(row_indx,index_p)] ; [zeros(n_both, n_pred) eye(n_both) zeros(n_both, n_both + n_fwrd)]]; - E = [-aa(row_indx,[index_m index_0p]) ; [zeros(n_both, n_both + n_pred) eye(n_both, n_both + n_fwrd) ] ]; - - [err, ss, tt, w, sdim, data(i).eigval, info1] = mjdgges(E,D,options_.qz_criterium); - - if (verbose) - disp('eigval'); - disp(data(i).eigval); - end; - if info1 - info(1) = 2; - info(2) = info1; - return + if task ~= 1 && options_.dr_cycle_reduction == 1 + A1 = [aa(row_indx,index_m ) zeros(n_dynamic,n_fwrd)]; + B1 = [aa(row_indx,index_0m) aa(row_indx,index_0p) ]; + C1 = [zeros(n_dynamic,n_pred) aa(row_indx,index_p)]; + [ghx, info] = cycle_reduction(A1, B1, C1, options_.dr_cycle_reduction_tol); + %ghx + ghx = ghx(:,index_m); + hx = ghx(1:n_pred+n_both,:); + gx = ghx(1+n_pred:end,:); end - nba = nd-sdim; - if task == 1 - data(i).rank = rank(w(nd-nyf+1:end,nd-nyf+1:end)); - dr.rank = dr.rank + data(i).rank; - if ~exist('OCTAVE_VERSION','builtin') - data(i).eigval = eig(E,D); + + if (task ~= 1 && ((options_.dr_cycle_reduction == 1 && info ==1) || options_.dr_cycle_reduction == 0)) || task == 1 + D = [[aa(row_indx,index_0m) zeros(n_dynamic,n_both) aa(row_indx,index_p)] ; [zeros(n_both, n_pred) eye(n_both) zeros(n_both, n_both + n_fwrd)]]; + E = [-aa(row_indx,[index_m index_0p]) ; [zeros(n_both, n_both + n_pred) eye(n_both, n_both + n_fwrd) ] ]; + + [err, ss, tt, w, sdim, data(i).eigval, info1] = mjdgges(E,D,options_.qz_criterium); + + if (verbose) + disp('eigval'); + disp(data(i).eigval); + end; + if info1 + info(1) = 2; + info(2) = info1; + return end - dr.eigval = [dr.eigval ; data(i).eigval]; - end - if (verbose) - disp(['sum eigval > 1 = ' int2str(sum(abs(data(i).eigval) > 1.)) ' nyf=' int2str(nyf) ' and dr.rank=' int2str(data(i).rank)]); - disp(['data(' int2str(i) ').eigval']); - disp(data(i).eigval); + nba = nd-sdim; + if task == 1 + data(i).rank = rank(w(nd-nyf+1:end,nd-nyf+1:end)); + dr.rank = dr.rank + data(i).rank; + if ~exist('OCTAVE_VERSION','builtin') + data(i).eigval = eig(E,D); + end + dr.eigval = [dr.eigval ; data(i).eigval]; + end + if (verbose) + disp(['sum eigval > 1 = ' int2str(sum(abs(data(i).eigval) > 1.)) ' nyf=' int2str(nyf) ' and dr.rank=' int2str(data(i).rank)]); + disp(['data(' int2str(i) ').eigval']); + disp(data(i).eigval); + end; + + %First order approximation + if task ~= 1 + if nba ~= nyf + sorted_roots = sort(abs(dr.eigval)); + if isfield(options_,'indeterminacy_continuity') + if options_.indeterminacy_msv == 1 + [ss,tt,w,q] = qz(e',d'); + [ss,tt,w,junk] = reorder(ss,tt,w,q); + ss = ss'; + tt = tt'; + w = w'; + %nba = nyf; + end + else + if nba > nyf + temp = sorted_roots(nd-nba+1:nd-nyf)-1-options_.qz_criterium; + info(1) = 3; + elseif nba < nyf; + temp = sorted_roots(nd-nyf+1:nd-nba)-1-options_.qz_criterium; + info(1) = 4; + end + info(2) = temp'*temp; + return + end + end + indx_stable_root = 1: (nd - nyf); %=> index of stable roots + indx_explosive_root = n_pred + n_both + 1:nd; %=> index of explosive roots + % derivatives with respect to dynamic state variables + % forward variables + Z = w'; + Z11t = Z(indx_stable_root, indx_stable_root)'; + Z21 = Z(indx_explosive_root, indx_stable_root); + Z22 = Z(indx_explosive_root, indx_explosive_root); + if ~isfloat(Z21) && (condest(Z21) > 1e9) + % condest() fails on a scalar under Octave + info(1) = 5; + info(2) = condest(Z21); + return; + else + %gx = -inv(Z22) * Z21; + gx = - Z22 \ Z21; + end + + % predetermined variables + hx = Z11t * inv(tt(indx_stable_root, indx_stable_root)) * ss(indx_stable_root, indx_stable_root) * inv(Z11t); + + k1 = 1:(n_pred+n_both); + k2 = 1:(n_fwrd+n_both); + + ghx = [hx(k1,:); gx(k2(n_both+1:end),:)]; + end; end; - %First order approximation - if task ~= 1 - if nba ~= nyf - sorted_roots = sort(abs(dr.eigval)); - if isfield(options_,'indeterminacy_continuity') - if options_.indeterminacy_msv == 1 - [ss,tt,w,q] = qz(e',d'); - [ss,tt,w,junk] = reorder(ss,tt,w,q); - ss = ss'; - tt = tt'; - w = w'; - %nba = nyf; - end - else - if nba > nyf - temp = sorted_roots(nd-nba+1:nd-nyf)-1-options_.qz_criterium; - info(1) = 3; - elseif nba < nyf; - temp = sorted_roots(nd-nyf+1:nd-nba)-1-options_.qz_criterium; - info(1) = 4; - end - info(2) = temp'*temp; - return - end - end - indx_stable_root = 1: (nd - nyf); %=> index of stable roots - indx_explosive_root = n_pred + n_both + 1:nd; %=> index of explosive roots - % derivatives with respect to dynamic state variables - % forward variables - Z = w'; - Z11t = Z(indx_stable_root, indx_stable_root)'; - Z21 = Z(indx_explosive_root, indx_stable_root); - Z22 = Z(indx_explosive_root, indx_explosive_root); - if ~isfloat(Z21) && (condest(Z21) > 1e9) - % condest() fails on a scalar under Octave - info(1) = 5; - info(2) = condest(Z21); - return; - else - %gx = -inv(Z22) * Z21; - gx = - Z22 \ Z21; - end - - % predetermined variables - hx = Z11t * inv(tt(indx_stable_root, indx_stable_root)) * ss(indx_stable_root, indx_stable_root) * inv(Z11t); - - k1 = 1:(n_pred+n_both); - k2 = 1:(n_fwrd+n_both); - - ghx = [hx(k1,:); gx(k2(n_both+1:end),:)]; - + if task~= 1 %lead variables actually present in the model - j4 = n_static+n_pred+1:n_static+n_pred+n_both+n_fwrd; % Index on the forward and both variables j3 = nonzeros(lead_lag_incidence(2,j4)) - n_static - 2 * n_pred - n_both; % Index on the non-zeros forward and both variables j4 = find(lead_lag_incidence(2,j4)); diff --git a/matlab/global_initialization.m b/matlab/global_initialization.m index 1711db5b3..28f77df91 100644 --- a/matlab/global_initialization.m +++ b/matlab/global_initialization.m @@ -436,13 +436,20 @@ options_.sylvester_fixed_point_tol = 1e-12; 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) +% if 1 use a square 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; +% if equal to 1 use a cycle reduction method to compute the decision rule (for large scale models) +options_.dr_cycle_reduction = 0; + +% convergence criterion for iteratives methods to solve the decision rule +options_.dr_cycle_reduction_tol = 1e-7; + + % dates for historical time series options_.initial_date.freq = 1; options_.initial_date.period = 1; diff --git a/preprocessor/DynareBison.yy b/preprocessor/DynareBison.yy index 738e6fe5e..83ddd5563 100644 --- a/preprocessor/DynareBison.yy +++ b/preprocessor/DynareBison.yy @@ -95,8 +95,8 @@ class ParsingDriver; %token BVAR_PRIOR_DECAY BVAR_PRIOR_FLAT BVAR_PRIOR_LAMBDA %token BVAR_PRIOR_MU BVAR_PRIOR_OMEGA BVAR_PRIOR_TAU BVAR_PRIOR_TRAIN %token BVAR_REPLIC BYTECODE -%token CALIB_SMOOTHER CHANGE_TYPE CHECK CONDITIONAL_FORECAST CONDITIONAL_FORECAST_PATHS CONF_SIG CONSTANT CONTROLLED_VAREXO CORR COVAR CUTOFF -%token DATAFILE FILE DOUBLING DR_ALGO DROP DSAMPLE DYNASAVE DYNATYPE CALIBRATION +%token CALIB_SMOOTHER CHANGE_TYPE CHECK CONDITIONAL_FORECAST CONDITIONAL_FORECAST_PATHS CONF_SIG CONSTANT CONTROLLED_VAREXO CORR COVAR CUTOFF CYCLE_REDUCTION +%token DATAFILE FILE DOUBLING DR DR_CYCLE_REDUCTION_TOL DR_ALGO DROP DSAMPLE DYNASAVE DYNATYPE CALIBRATION %token END ENDVAL EQUAL ESTIMATION ESTIMATED_PARAMS ESTIMATED_PARAMS_BOUNDS ESTIMATED_PARAMS_INIT EXTENDED_PATH %token FILENAME FILTER_STEP_AHEAD FILTERED_VARS FIRST_OBS LAST_OBS SET_TIME %token FLOAT_NUMBER @@ -937,6 +937,8 @@ stoch_simul_options : o_dr_algo | o_pruning | o_sylvester | o_sylvester_fixed_point_tol + | o_dr + | o_dr_cycle_reduction_tol ; symbol_list : symbol_list symbol @@ -1504,6 +1506,8 @@ estimation_options : o_datafile | o_lyapunov | o_lyapunov_fixed_point_tol | o_lyapunov_doubling_tol + | o_dr + | o_dr_cycle_reduction_tol | o_analytic_derivation ; @@ -2313,6 +2317,9 @@ o_lyapunov : LYAPUNOV EQUAL FIXED_POINT {driver.option_num("lyapunov_fp", "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_dr : DR EQUAL CYCLE_REDUCTION {driver.option_num("dr_cycle_reduction", "1"); } + | DR EQUAL DEFAULT {driver.option_num("dr_cycle_reduction", "0"); }; +o_dr_cycle_reduction_tol : DR_CYCLE_REDUCTION_TOL EQUAL non_negative_number {driver.option_num("dr_cycle_reduction_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); }; diff --git a/preprocessor/DynareFlex.ll b/preprocessor/DynareFlex.ll index 5e96dc374..c5413092f 100644 --- a/preprocessor/DynareFlex.ll +++ b/preprocessor/DynareFlex.ll @@ -302,6 +302,7 @@ string eofbuff; fixed_point {return token::FIXED_POINT;} doubling {return token::DOUBLING;} square_root_solver {return token::SQUARE_ROOT_SOLVER;} +cycle_reduction {return token::CYCLE_REDUCTION;} default {return token::DEFAULT;} alpha { yylval->string_val = new string(yytext); @@ -505,9 +506,11 @@ string eofbuff; order {return token::ORDER;} sylvester {return token::SYLVESTER;} lyapunov {return token::LYAPUNOV;} +dr {return token::DR;} sylvester_fixed_point_tol {return token::SYLVESTER_FIXED_POINT_TOL;} lyapunov_fixed_point_tol {return token::LYAPUNOV_FIXED_POINT_TOL;} lyapunov_doubling_tol {return token::LYAPUNOV_DOUBLING_TOL;} +dr_cycle_reduction_tol {return token::DR_CYCLE_REDUCTION_TOL;} replic {return token::REPLIC;} ar {return token::AR;} nofunctions {return token::NOFUNCTIONS;}