stepan
/
dynare
1
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

Adds the cycle reduction algorithm to solve the polynomial equation for retrieving the coefficients 

associated to the endogenous variables in the decision rule.
time-shift
Ferhat Mihoubi 2012-07-01 15:15:52 +02:00
parent 87ad577347
commit 4488357f59
6 changed files with 203 additions and 80 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

23
doc/dynare.texi
View File

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

69
matlab/cycle_reduction.m Normal file
View File

@ -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 <http://www.gnu.org/licenses/>.
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;

168
matlab/dr_block.m
View File

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

9
matlab/global_initialization.m
View File

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

11
preprocessor/DynareBison.yy
View File

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

3
preprocessor/DynareFlex.ll
View File

@ -302,6 +302,7 @@ string eofbuff;
<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>cycle_reduction {return token::CYCLE_REDUCTION;}
<DYNARE_STATEMENT>default {return token::DEFAULT;}
<DYNARE_STATEMENT>alpha {
yylval->string_val = new string(yytext);
@ -505,9 +506,11 @@ string eofbuff;
<DYNARE_STATEMENT>order {return token::ORDER;}
<DYNARE_STATEMENT>sylvester {return token::SYLVESTER;}
<DYNARE_STATEMENT>lyapunov {return token::LYAPUNOV;}
<DYNARE_STATEMENT>dr {return token::DR;}
<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>dr_cycle_reduction_tol {return token::DR_CYCLE_REDUCTION_TOL;}
<DYNARE_STATEMENT>replic {return token::REPLIC;}
<DYNARE_STATEMENT>ar {return token::AR;}
<DYNARE_STATEMENT>nofunctions {return token::NOFUNCTIONS;}