From 8f1a4cb3637fe4936ade877c88e6a120af2b5f24 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?S=C3=A9bastien=20Villemot?= Date: Fri, 13 Jan 2023 16:19:53 +0100 Subject: [PATCH] Drop block-decomposed first-order perturbation solution and kalman filter --- doc/manual/source/the-model-file.rst | 20 +- matlab/DsgeSmoother.m | 8 +- matlab/check.m | 16 +- matlab/disp_dr.m | 13 +- matlab/dr_block.m | 717 --------------- matlab/dsge_likelihood.m | 35 +- matlab/dynare_estimation_init.m | 72 +- matlab/dynare_resolve.m | 8 +- matlab/dynvars_from_endo_simul.m | 26 - matlab/evaluate_steady_state.m | 5 - matlab/model_diagnostics.m | 75 +- .../solve_block_decomposed_problem.m | 2 +- matlab/resol.m | 10 +- matlab/set_state_space.m | 8 +- matlab/simult_.m | 18 +- matlab/solve_two_boundaries.m | 4 +- matlab/th_autocovariances.m | 60 +- mex/build/block_kalman_filter.am | 15 - mex/build/matlab/Makefile.am | 2 +- .../matlab/block_kalman_filter/Makefile.am | 2 - mex/build/matlab/configure.ac | 1 - mex/build/octave/Makefile.am | 2 +- .../octave/block_kalman_filter/Makefile.am | 3 - mex/build/octave/configure.ac | 1 - mex/sources/Makefile.am | 1 - .../block_kalman_filter.cc | 850 ------------------ .../block_kalman_filter.hh | 68 -- preprocessor | 2 +- tests/Makefile.am | 4 - tests/block_bytecode/lola_common.inc | 7 - .../lola_solve_one_boundary.mod | 1 - .../lola_solve_one_boundary_mfs1.mod | 1 - .../lola_solve_one_boundary_mfs2.mod | 1 - .../lola_solve_one_boundary_mfs3.mod | 1 - tests/block_bytecode/lola_stochastic.mod | 5 - .../block_bytecode/lola_stochastic_block.mod | 26 - 36 files changed, 121 insertions(+), 1969 deletions(-) delete mode 100644 matlab/dr_block.m delete mode 100644 matlab/dynvars_from_endo_simul.m delete mode 100644 mex/build/block_kalman_filter.am delete mode 100644 mex/build/matlab/block_kalman_filter/Makefile.am delete mode 100644 mex/build/octave/block_kalman_filter/Makefile.am delete mode 100644 mex/sources/block_kalman_filter/block_kalman_filter.cc delete mode 100644 mex/sources/block_kalman_filter/block_kalman_filter.hh delete mode 100644 tests/block_bytecode/lola_stochastic.mod delete mode 100644 tests/block_bytecode/lola_stochastic_block.mod diff --git a/doc/manual/source/the-model-file.rst b/doc/manual/source/the-model-file.rst index 6594f21c9..740ab0197 100644 --- a/doc/manual/source/the-model-file.rst +++ b/doc/manual/source/the-model-file.rst @@ -4310,10 +4310,6 @@ Computing the stochastic solution When a list of ``VARIABLE_NAME`` is specified, results are displayed only for these variables. - The ``stoch_simul`` command with a first order approximation can - benefit from the block decomposition of the model (see - :opt:`block`). - *Options* .. option:: ar = INTEGER @@ -4982,8 +4978,8 @@ which is described below. the endogenous variables are generated by assuming that the agents believe that there will no more shocks after period :math:`t+S`. This is an experimental feature and can be quite - slow. A non-zero value is not compatible with either the - ``bytecode`` or the ``block`` option of the ``model`` block. + slow. A non-zero value is not compatible with the ``bytecode`` + option of the ``model`` block. Default: ``0``. .. option:: hybrid @@ -5110,15 +5106,10 @@ The coefficients of the decision rules are stored as follows: to all endogenous in the declaration order. * :math:`A` is stored in ``oo_.dr.ghx``. The matrix rows correspond to all endogenous in DR-order. The matrix columns correspond to state - variables in DR-order, as given by ``oo_.dr.state_var``. (N.B.: if the - ``block`` option to the ``model`` block has been specified, then rows - are in declaration order, and columns are ordered - according to ``oo_.dr.state_var`` which may differ from DR-order.) + variables in DR-order, as given by ``oo_.dr.state_var``. * :math:`B` is stored ``oo_.dr.ghu``. The matrix rows correspond to all endogenous in DR-order. The matrix columns correspond to - exogenous variables in declaration order. (N.B.: if the - ``block`` option to the ``model`` block has been specified, then rows - are in declaration order.) + exogenous variables in declaration order. Of course, the shown form of the approximation is only stylized, because it neglects the required different ordering in :math:`y^s` and @@ -5707,9 +5698,6 @@ Note that in order to avoid stochastic singularity, you must have at least as many shocks or measurement errors in your model as you have observed variables. -The estimation using a first order approximation can benefit from the -block decomposition of the model (see :opt:`block`). - .. _varobs: .. command:: varobs VARIABLE_NAME...; diff --git a/matlab/DsgeSmoother.m b/matlab/DsgeSmoother.m index 7d3ab0284..7ea0cd462 100644 --- a/matlab/DsgeSmoother.m +++ b/matlab/DsgeSmoother.m @@ -57,7 +57,7 @@ function [alphahat,etahat,epsilonhat,ahat,SteadyState,trend_coeff,aK,T,R,P,PK,de % SPECIAL REQUIREMENTS % None -% Copyright © 2006-2020 Dynare Team +% Copyright © 2006-2023 Dynare Team % % This file is part of Dynare. % @@ -359,11 +359,7 @@ if ~options_.smoother_redux oo_.dr.restrict_var_list = oldoo.restrict_var_list; oo_.dr.restrict_columns = oldoo.restrict_columns; else - if options_.block == 0 - ic = [ M_.nstatic+(1:M_.nspred) M_.endo_nbr+(1:size(oo_.dr.ghx,2)-M_.nspred) ]'; - else - ic = oo_.dr.restrict_columns; - end + ic = [ M_.nstatic+(1:M_.nspred) M_.endo_nbr+(1:size(oo_.dr.ghx,2)-M_.nspred) ]'; if isempty(options_.nk) nk=1; diff --git a/matlab/check.m b/matlab/check.m index 684325014..ec9ec6513 100644 --- a/matlab/check.m +++ b/matlab/check.m @@ -12,7 +12,7 @@ function [eigenvalues_,result,info] = check(M, options, oo) % - result [integer] scalar, equal to 1 if Blanchard and Kahn conditions are satisfied, zero otherwise. % - info [integer] scalar or vector, error code as returned by resol routine. -% Copyright © 2001-2019 Dynare Team +% Copyright © 2001-2023 Dynare Team % % This file is part of Dynare. % @@ -51,18 +51,8 @@ end eigenvalues_ = dr.eigval; [m_lambda,i]=sort(abs(eigenvalues_)); -% Count number of forward looking variables -if ~options.block - nyf = M.nsfwrd; -else - nyf = 0; - for j = 1:length(M.block_structure.block) - nyf = nyf + M.block_structure.block(j).n_forward + M.block_structure.block(j).n_mixed; - end -end - result = 0; -if (nyf == dr.edim) && (dr.full_rank) +if (M.nsfwrd == dr.edim) && (dr.full_rank) result = 1; end @@ -73,7 +63,7 @@ if ~options.noprint z=[m_lambda real(eigenvalues_(i)) imag(eigenvalues_(i))]'; disp(sprintf('%16.4g %16.4g %16.4g\n',z)) disp(sprintf('\nThere are %d eigenvalue(s) larger than 1 in modulus ', dr.edim)); - disp(sprintf('for %d forward-looking variable(s)',nyf)); + disp(sprintf('for %d forward-looking variable(s)', M.nsfwrd)); skipline() if result disp('The rank condition is verified.') diff --git a/matlab/disp_dr.m b/matlab/disp_dr.m index b959ad7b9..7c95714c5 100644 --- a/matlab/disp_dr.m +++ b/matlab/disp_dr.m @@ -10,7 +10,7 @@ function disp_dr(dr,order,var_list) % OUTPUTS % none -% Copyright © 2001-2019 Dynare Team +% Copyright © 2001-2023 Dynare Team % % This file is part of Dynare. % @@ -36,14 +36,9 @@ end nx =size(dr.ghx,2); nu =size(dr.ghu,2); -if options_.block - k = M_.nspred; - k1 = 1:M_.endo_nbr; -else - k = find(dr.kstate(:,2) <= M_.maximum_lag+1); - klag = dr.kstate(k,[1 2]); - k1 = dr.order_var; -end +k = find(dr.kstate(:,2) <= M_.maximum_lag+1); +klag = dr.kstate(k,[1 2]); +k1 = dr.order_var; if isempty(var_list) var_list = M_.endo_names(1:M_.orig_endo_nbr); diff --git a/matlab/dr_block.m b/matlab/dr_block.m deleted file mode 100644 index 6444523ec..000000000 --- a/matlab/dr_block.m +++ /dev/null @@ -1,717 +0,0 @@ -function [dr,info,M_,oo_] = dr_block(dr,task,M_,options_,oo_,varargin) -% function [dr,info,M_,oo_] = dr_block(dr,task,M_,options_,oo_,varargin) -% computes the reduced form solution of a rational expectations block-decomposed model -% (first order approximation of the stochastic model around the deterministic steady state). -% -% NB: This code does not work with option mfs > 0. The preprocessor has a check to avoid this -% configuration. See also #1726. -% -% INPUTS -% dr [matlab structure] Decision rules for stochastic simulations. -% task [integer] if task = 0 then dr_block computes decision rules. -% if task = 1 then dr_block computes eigenvalues. -% M_ [matlab structure] Definition of the model. -% options_ [matlab structure] Global options. -% oo_ [matlab structure] Results -% oo_ [matlab cell] Other input arguments -% -% OUTPUTS -% dr [matlab structure] Decision rules for stochastic simulations. -% info [integer] info=1: the model doesn't define current variables uniquely -% info=2: problem in mjdgges.dll info(2) contains error code. -% info=3: BK order condition not satisfied info(2) contains "distance" -% absence of stable trajectory. -% info=4: BK order condition not satisfied info(2) contains "distance" -% indeterminacy. -% info=5: BK rank condition not satisfied. -% info=6: The jacobian matrix evaluated at the steady state is complex. -% M_ [matlab structure] -% oo_ [matlab structure] -% -% ALGORITHM -% first order block relaxation method applied to the model -% E[A Yt-1 + B Yt + C Yt+1 + ut] = 0 -% -% SPECIAL REQUIREMENTS -% none. -% - -% Copyright © 2010-2023 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 . - -info = 0; -verbose = 0; -if nargin > 5 - verbose = varargin{1}; -end -%verbose = options_.verbosity; -if options_.order > 1 - error('2nd and 3rd order approximation not implemented with block option') -else - if options_.loglinear - error('The loglinear option is not yet supported in first order approximation for a block decomposed model'); - end -end - -z = repmat(dr.ys,1,M_.maximum_lead + M_.maximum_lag + 1); -zx = repmat([oo_.exo_simul oo_.exo_det_simul],M_.maximum_lead + M_.maximum_lag + 1, 1); -if isempty(zx) - zx = [repmat(oo_.exo_steady_state',M_.maximum_lead + M_.maximum_lag + 1,1) repmat(oo_.exo_det_steady_state',M_.maximum_lead + M_.maximum_lag + 1,1)]; -end - -if ~isfield(M_,'block_structure') - error('Option ''block'' has not been specified') -end -data = M_.block_structure.block; - -if options_.bytecode - [~, data]= bytecode('dynamic', 'evaluate', 'block_decomposed', z, zx, M_.params, dr.ys, 1, data); -else - T=NaN(M_.block_structure.dyn_tmp_nbr, 1); - it_=M_.maximum_lag+1; - y=dynvars_from_endo_simul(z, it_, M_); - for blk = 1:length(M_.block_structure.block) - [~, ~, T, data(blk).g1, data(blk).g1_x, data(blk).g1_xd, data(blk).g1_o]=feval([M_.fname '.dynamic'], blk, y, zx, M_.params, dr.ys, T, it_, true); - end -end -dr.full_rank = 1; -dr.eigval = []; -dr.nd = 0; - -dr.ghx = []; -dr.ghu = []; -%Determine the global list of state variables: -dr.state_var = M_.state_var; -M_.block_structure.state_var = dr.state_var; -n_sv = size(dr.state_var, 2); -dr.ghx = zeros(M_.endo_nbr, length(dr.state_var)); -dr.exo_var = 1:M_.exo_nbr; -dr.ghu = zeros(M_.endo_nbr, M_.exo_nbr); -for i = 1:length(data) - ghx = []; - indexi_0 = 0; - if (verbose) - disp('======================================================================'); - disp(['Block ' int2str(i)]); - disp('-----------'); - data(i) - end - n_pred = data(i).n_backward; - n_fwrd = data(i).n_forward; - n_both = data(i).n_mixed; - n_static = data(i).n_static; - nd = n_pred + n_fwrd + 2*n_both; - dr.nd = dr.nd + nd; - n_dynamic = n_pred + n_fwrd + n_both; - exo_nbr = M_.block_structure.block(i).exo_nbr; - exo_det_nbr = M_.block_structure.block(i).exo_det_nbr; - other_endo_nbr = M_.block_structure.block(i).other_endo_nbr; - jacob = full(data(i).g1); - lead_lag_incidence = data(i).lead_lag_incidence; - endo = data(i).variable; - exo = data(i).exogenous; - if (verbose) - disp('jacob'); - disp(jacob); - disp('lead_lag_incidence'); - disp(lead_lag_incidence); - end - maximum_lag = data(i).maximum_endo_lag; - maximum_lead = data(i).maximum_endo_lead; - n = n_dynamic + n_static; - - block_type = M_.block_structure.block(i).Simulation_Type; - if task ~= 1 - if block_type == 2 || block_type == 4 || block_type == 7 - block_type = 8; - end - end - if maximum_lag > 0 && (n_pred > 0 || n_both > 0) && block_type ~= 1 - indexi_0 = min(lead_lag_incidence(2,:)); - end - switch block_type - case 1 - %% ------------------------------------------------------------------ - %Evaluate Forward - if maximum_lag > 0 && n_pred > 0 - indx_r = find(M_.block_structure.block(i).lead_lag_incidence(1,:)); - indx_c = M_.block_structure.block(i).lead_lag_incidence(1,indx_r); - data(i).eigval = diag(jacob(indx_r, indx_c)); - data(i).rank = 0; - else - data(i).eigval = []; - data(i).rank = 0; - end - dr.eigval = [dr.eigval ; data(i).eigval]; - %First order approximation - if task ~= 1 - [tmp1, tmp2, indx_c] = find(M_.block_structure.block(i).lead_lag_incidence(2,:)); - B = jacob(:,indx_c); - if (maximum_lag > 0 && n_pred > 0) - [indx_r, tmp1, indx_r_v] = find(M_.block_structure.block(i).lead_lag_incidence(1,:)); - ghx = - B \ jacob(:,indx_r_v); - end - if other_endo_nbr - fx = data(i).g1_o; - % retrieves the derivatives with respect to endogenous - % variable belonging to previous blocks - fx_tm1 = zeros(n,other_endo_nbr); - fx_t = zeros(n,other_endo_nbr); - fx_tp1 = zeros(n,other_endo_nbr); - % stores in fx_tm1 the lagged values of fx - [r, c, lag] = find(data(i).lead_lag_incidence_other(1,:)); - fx_tm1(:,c) = fx(:,lag); - % stores in fx the current values of fx - [r, c, curr] = find(data(i).lead_lag_incidence_other(2,:)); - fx_t(:,c) = fx(:,curr); - % stores in fx_tp1 the leaded values of fx - [r, c, lead] = find(data(i).lead_lag_incidence_other(3,:)); - fx_tp1(:,c) = fx(:,lead); - - l_x = dr.ghx(data(i).other_endogenous,:); - l_x_sv = dr.ghx(dr.state_var, 1:n_sv); - - selector_tm1 = M_.block_structure.block(i).tm1; - - ghx_other = - B \ (fx_t * l_x + (fx_tp1 * l_x * l_x_sv) + fx_tm1 * selector_tm1); - dr.ghx(endo, :) = dr.ghx(endo, :) + ghx_other; - end - - if exo_nbr - fu = data(i).g1_x; - exo = dr.exo_var; - if other_endo_nbr > 0 - l_u_sv = dr.ghu(dr.state_var,:); - l_x = dr.ghx(data(i).other_endogenous,:); - l_u = dr.ghu(data(i).other_endogenous,:); - fu_complet = zeros(n, M_.exo_nbr); - fu_complet(:,data(i).exogenous) = fu; - ghu = - B \ (fu_complet + fx_tp1 * l_x * l_u_sv + (fx_t) * l_u ); - else - fu_complet = zeros(n, M_.exo_nbr); - fu_complet(:,data(i).exogenous) = fu; - ghu = - B \ fu_complet; - end - else - exo = dr.exo_var; - if other_endo_nbr > 0 - l_u_sv = dr.ghu(dr.state_var,:); - l_x = dr.ghx(data(i).other_endogenous,:); - l_u = dr.ghu(data(i).other_endogenous,:); - ghu = -B \ (fx_tp1 * l_x * l_u_sv + (fx_t) * l_u ); - else - ghu = []; - end - end - end - case 2 - %% ------------------------------------------------------------------ - %Evaluate Backward - if maximum_lead > 0 && n_fwrd > 0 - indx_r = find(M_.block_structure.block(i).lead_lag_incidence(3,:)); - indx_c = M_.block_structure.block(i).lead_lag_incidence(3,indx_r); - data(i).eigval = 1 ./ diag(jacob(indx_r, indx_c)); - data(i).rank = sum(abs(data(i).eigval) > 0); - full_rank = (rcond(jacob(indx_r, indx_c)) > 1e-9); - else - data(i).eigval = []; - data(i).rank = 0; - full_rank = 1; - end - dr.eigval = [dr.eigval ; data(i).eigval]; - dr.full_rank = dr.full_rank && full_rank; - %First order approximation - if task ~= 1 - if (maximum_lag > 0) - indx_r = find(M_.block_structure.block(i).lead_lag_incidence(3,:)); - indx_c = M_.block_structure.block(i).lead_lag_incidence(3,indx_r); - ghx = - inv(jacob(indx_r, indx_c)); - end - ghu = - inv(jacob(indx_r, indx_c)) * data(i).g1_x; - end - case 3 - %% ------------------------------------------------------------------ - %Solve Forward single equation - if maximum_lag > 0 && n_pred > 0 - data(i).eigval = - jacob(1 , 1 : n_pred) / jacob(1 , n_pred + n_static + 1 : n_pred + n_static + n_pred + n_both); - data(i).rank = 0; - else - data(i).eigval = []; - data(i).rank = 0; - end - dr.eigval = [dr.eigval ; data(i).eigval]; - %First order approximation - if task ~= 1 - if (maximum_lag > 0) - ghx = - jacob(1 , 1 : n_pred) / jacob(1 , n_pred + n_static + 1 : n_pred + n_static + n_pred + n_both); - else - ghx = 0; - end - if other_endo_nbr - fx = data(i).g1_o; - % retrieves the derivatives with respect to endogenous - % variable belonging to previous blocks - fx_tm1 = zeros(n,other_endo_nbr); - fx_t = zeros(n,other_endo_nbr); - fx_tp1 = zeros(n,other_endo_nbr); - % stores in fx_tm1 the lagged values of fx - [r, c, lag] = find(data(i).lead_lag_incidence_other(1,:)); - fx_tm1(:,c) = fx(:,lag); - % stores in fx the current values of fx - [r, c, curr] = find(data(i).lead_lag_incidence_other(2,:)); - fx_t(:,c) = fx(:,curr); - % stores in fx_tm1 the leaded values of fx - [r, c, lead] = find(data(i).lead_lag_incidence_other(3,:)); - fx_tp1(:,c) = fx(:,lead); - - l_x = dr.ghx(data(i).other_endogenous,:); - l_x_sv = dr.ghx(dr.state_var, 1:n_sv); - - selector_tm1 = M_.block_structure.block(i).tm1; - ghx_other = - (fx_t * l_x + (fx_tp1 * l_x * l_x_sv) + fx_tm1 * selector_tm1) / jacob(1 , n_pred + 1 : n_pred + n_static + n_pred + n_both); - dr.ghx(endo, :) = dr.ghx(endo, :) + ghx_other; - - end - if exo_nbr - fu = data(i).g1_x; - if other_endo_nbr > 0 - l_u_sv = dr.ghu(dr.state_var,:); - l_x = dr.ghx(data(i).other_endogenous,:); - l_u = dr.ghu(data(i).other_endogenous,:); - fu_complet = zeros(n, M_.exo_nbr); - fu_complet(:,data(i).exogenous) = fu; - ghu = -(fu_complet + fx_tp1 * l_x * l_u_sv + (fx_t) * l_u ) / jacob(1 , n_pred + 1 : n_pred + n_static + n_pred + n_both); - exo = dr.exo_var; - else - ghu = - fu / jacob(1 , n_pred + 1 : n_pred + n_static + n_pred + n_both); - end - else - if other_endo_nbr > 0 - l_u_sv = dr.ghu(dr.state_var,:); - l_x = dr.ghx(data(i).other_endogenous,:); - l_u = dr.ghu(data(i).other_endogenous,:); - ghu = -(fx_tp1 * l_x * l_u_sv + (fx_t) * l_u ) / jacob(1 , n_pred + 1 : n_pred + n_static + n_pred + n_both); - exo = dr.exo_var; - else - ghu = []; - end - end - end - case 4 - %% ------------------------------------------------------------------ - %Solve Backward single equation - if maximum_lead > 0 && n_fwrd > 0 - data(i).eigval = - jacob(1 , n_pred + n - n_fwrd + 1 : n_pred + n) / jacob(1 , n_pred + n + 1 : n_pred + n + n_fwrd) ; - data(i).rank = sum(abs(data(i).eigval) > 0); - full_rank = (abs(jacob(1,n_pred+n+1: n_pred+n+n_fwrd)) > 1e-9); - else - data(i).eigval = []; - data(i).rank = 0; - full_rank = 1; - end - dr.full_rank = dr.full_rank && full_rank; - dr.eigval = [dr.eigval ; data(i).eigval]; - case 6 - %% ------------------------------------------------------------------ - %Solve Forward complete - if (maximum_lag > 0) - ghx = - jacob(: , n_pred + 1 : n_pred + n_static ... - + n_pred + n_both) \ jacob(: , 1 : n_pred); - else - ghx = 0; - end - if maximum_lag > 0 && n_pred > 0 - data(i).eigval = -eig(ghx(n_static+1:end,:)); - data(i).rank = 0; - full_rank = (rcond(ghx(n_static+1:end,:)) > 1e-9); - else - data(i).eigval = []; - data(i).rank = 0; - full_rank = 1; - end - dr.eigval = [dr.eigval ; data(i).eigval]; - dr.full_rank = dr.full_rank && full_rank; - if task ~= 1 - if other_endo_nbr - fx = data(i).g1_o; - % retrieves the derivatives with respect to endogenous - % variable belonging to previous blocks - fx_tm1 = zeros(n,other_endo_nbr); - fx_t = zeros(n,other_endo_nbr); - fx_tp1 = zeros(n,other_endo_nbr); - % stores in fx_tm1 the lagged values of fx - [r, c, lag] = find(data(i).lead_lag_incidence_other(1,:)); - fx_tm1(:,c) = fx(:,lag); - % stores in fx the current values of fx - [r, c, curr] = find(data(i).lead_lag_incidence_other(2,:)); - fx_t(:,c) = fx(:,curr); - % stores in fx_tm1 the leaded values of fx - [r, c, lead] = find(data(i).lead_lag_incidence_other(3,:)); - fx_tp1(:,c) = fx(:,lead); - - l_x = dr.ghx(data(i).other_endogenous,:); - l_x_sv = dr.ghx(dr.state_var, 1:n_sv); - - selector_tm1 = M_.block_structure.block(i).tm1; - ghx_other = - (fx_t * l_x + (fx_tp1 * l_x * l_x_sv) + fx_tm1 * selector_tm1) / jacob(: , n_pred + 1 : n_pred + n_static + n_pred + n_both); - dr.ghx(endo, :) = dr.ghx(endo, :) + ghx_other; - end - if exo_nbr - fu = data(i).g1_x; - if other_endo_nbr > 0 - l_u_sv = dr.ghu(dr.state_var,:); - l_x = dr.ghx(data(i).other_endogenous,:); - l_u = dr.ghu(data(i).other_endogenous,:); - fu_complet = zeros(n, M_.exo_nbr); - fu_complet(:,data(i).exogenous) = fu; - ghu = -(fu_complet + fx_tp1 * l_x * l_u_sv + (fx_t) * l_u ) / jacob(: , n_pred + 1 : n_pred + n_static + n_pred + n_both); - exo = dr.exo_var; - else - ghu = - fu / jacob(: , n_pred + 1 : n_pred + n_static + n_pred + n_both); - end - else - if other_endo_nbr > 0 - l_u_sv = dr.ghu(dr.state_var,:); - l_x = dr.ghx(data(i).other_endogenous,:); - l_u = dr.ghu(data(i).other_endogenous,:); - ghu = -(fx_tp1 * l_x * l_u_sv + (fx_t) * l_u ) / jacob(1 , n_pred + 1 : n_pred + n_static + n_pred + n_both); - exo = dr.exo_var; - else - ghu = []; - end - end - end - case 7 - %% ------------------------------------------------------------------ - %Solve Backward complete - if maximum_lead > 0 && n_fwrd > 0 - tmp = -jacob(: , n_pred + 1 : n_pred + n) \ ... - jacob(: , n_pred + n + 1 : n_pred + n + n_fwrd); - data(i).eigval = 1 ./ eig(tmp(n_static+1:end, :)); - data(i).rank = sum(abs(data(i).eigval) > 0); - full_rank = (rcond(jacob(: , n_pred + 1 : n_pred + n)) > 1e-9); - else - data(i).eigval = []; - data(i).rank = 0; - full_rank = 1; - end - dr.full_rank = dr.full_rank && full_rank; - dr.eigval = [dr.eigval ; data(i).eigval]; - case {5,8} - %% ------------------------------------------------------------------ - %The lead_lag_incidence contains columns in the following order: - % static variables, backward variable, mixed variables and forward variables - % - %Proceeds to a QR decomposition on the jacobian matrix in order to reduce the problem size - index_c = lead_lag_incidence(2,:); % Index of all endogenous variables present at time=t - index_s = lead_lag_incidence(2,1:n_static); % Index of all static endogenous variables present at time=t - if n_static > 0 - [Q, ~] = qr(jacob(:,index_s)); - aa = Q'*jacob; - else - aa = jacob; - end - index_0m = (n_static+1:n_static+n_pred) + indexi_0 - 1; - index_0p = (n_static+n_pred+1:n) + indexi_0 - 1; - index_m = 1:(n_pred+n_both); - index_p = lead_lag_incidence(3,find(lead_lag_incidence(3,:))); - nyf = n_fwrd + n_both; - A = aa(:,index_m); % Jacobain matrix for lagged endogeneous variables - B = aa(:,index_c); % Jacobian matrix for contemporaneous endogeneous variables - C = aa(:,index_p); % Jacobain matrix for led endogeneous variables - - row_indx = n_static+1:n; - - if task ~= 1 && options_.dr_cycle_reduction - 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 - - if (task ~= 1 && ((options_.dr_cycle_reduction && info ==1) || ~options_.dr_cycle_reduction)) || 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) ] ]; - - [ss, tt, w, sdim, data(i).eigval, info1] = mjdgges(E,D,options_.qz_criterium,options_.qz_zero_threshold); - - if (verbose) - disp('eigval'); - disp(data(i).eigval); - end - if info1 - if info1 == -30 - % one eigenvalue is close to 0/0 - info(1) = 7; - else - info(1) = 2; - info(2) = info1; - info(3) = size(E,2); - end - return - end - nba = nd-sdim; - if task == 1 - data(i).rank = rank(w(nd-nyf+1:end,nd-nyf+1:end)); - dr.full_rank = dr.full_rank && (rcond(w(nd-nyf+1:end,nd- ... - nyf+1:end)) > 1e-9); - 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 - if isfield(options_,'indeterminacy_continuity') - if options_.indeterminacy_msv == 1 - [ss,tt,w,q] = qz(E',D'); - [ss,tt,w,~] = reorder(ss,tt,w,q); - ss = ss'; - tt = tt'; - w = w'; - nba = nyf; - end - else - sorted_roots = sort(abs(data(i).eigval)); - 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 - - 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)); - - if n_static > 0 - B_static = B(:,1:n_static); % submatrix containing the derivatives w.r. to static variables - else - B_static = []; - end - %static variables, backward variable, mixed variables and forward variables - B_pred = B(:,n_static+1:n_static+n_pred+n_both); - B_fyd = B(:,n_static+n_pred+n_both+1:end); - - % static variables - if n_static > 0 - temp = - C(1:n_static,j3)*gx(j4,:)*hx; - j5 = index_m; - b = aa(:,index_c); - b10 = b(1:n_static, 1:n_static); - b11 = b(1:n_static, n_static+1:n); - temp(:,j5) = temp(:,j5)-A(1:n_static,:); - temp = b10\(temp-b11*ghx); - ghx = [temp; ghx]; - temp = []; - end - - A_ = real([B_static C(:,j3)*gx+B_pred B_fyd]); % The state_variable of the block are located at [B_pred B_both] - - if other_endo_nbr - if n_static > 0 - fx = Q' * data(i).g1_o; - else - fx = data(i).g1_o; - end - % retrieves the derivatives with respect to endogenous - % variable belonging to previous blocks - fx_tm1 = zeros(n,other_endo_nbr); - fx_t = zeros(n,other_endo_nbr); - fx_tp1 = zeros(n,other_endo_nbr); - % stores in fx_tm1 the lagged values of fx - [r, c, lag] = find(data(i).lead_lag_incidence_other(1,:)); - fx_tm1(:,c) = fx(:,lag); - % stores in fx the current values of fx - [r, c, curr] = find(data(i).lead_lag_incidence_other(2,:)); - fx_t(:,c) = fx(:,curr); - % stores in fx_tp1 the leaded values of fx - [r, c, lead] = find(data(i).lead_lag_incidence_other(3,:)); - fx_tp1(:,c) = fx(:,lead); - - l_x = dr.ghx(data(i).other_endogenous,:); - - l_x_sv = dr.ghx(dr.state_var, :); - - selector_tm1 = M_.block_structure.block(i).tm1; - - B_ = [zeros(size(B_static)) zeros(n,n_pred) C(:,j3) ]; - C_ = l_x_sv; - D_ = (fx_t * l_x + fx_tp1 * l_x * l_x_sv + fx_tm1 * selector_tm1 ); - % Solve the Sylvester equation: - % A_ * gx + B_ * gx * C_ + D_ = 0 - 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)); - elseif options_.sylvester_fp - ghx_other = gensylv_fp(A_, B_, C_, D_, i, options_.sylvester_fixed_point_tol); - else - ghx_other = gensylv(1, A_, B_, C_, -D_); - end - if options_.aim_solver ~= 1 - % Necessary when using Sims' routines for QZ - ghx_other = real(ghx_other); - end - - dr.ghx(endo, :) = dr.ghx(endo, :) + ghx_other; - end - - if exo_nbr - if n_static > 0 - fu = Q' * data(i).g1_x; - else - fu = data(i).g1_x; - end - - if other_endo_nbr > 0 - l_u_sv = dr.ghu(dr.state_var,:); - l_x = dr.ghx(data(i).other_endogenous,:); - l_u = dr.ghu(data(i).other_endogenous,:); - fu_complet = zeros(n, M_.exo_nbr); - fu_complet(:,data(i).exogenous) = fu; - % Solve the equation in ghu: - % A_ * ghu + (fu_complet + fx_tp1 * l_x * l_u_sv + (fx_t + B_ * ghx_other) * l_u ) = 0 - - ghu = -A_\ (fu_complet + fx_tp1 * l_x * l_u_sv + fx_t * l_u + B_ * ghx_other * l_u_sv ); - exo = dr.exo_var; - else - ghu = - A_ \ fu; - end - else - if other_endo_nbr > 0 - l_u_sv = dr.ghu(dr.state_var,:); - l_x = dr.ghx(data(i).other_endogenous,:); - l_u = dr.ghu(data(i).other_endogenous,:); - % Solve the equation in ghu: - % A_ * ghu + (fx_tp1 * l_x * l_u_sv + (fx_t + B_ * ghx_other) * l_u ) = 0 - ghu = -real(A_)\ (fx_tp1 * l_x * l_u_sv + (fx_t * l_u + B_ * ghx_other * l_u_sv) ); - exo = dr.exo_var; - else - ghu = []; - end - end - - - - if options_.loglinear - error('The loglinear option is not yet supported in first order approximation for a block decomposed model'); - % k = find(dr.kstate(:,2) <= M_.maximum_endo_lag+1); - % klag = dr.kstate(k,[1 2]); - % k1 = dr.order_var; - % - % ghx = repmat(1./dr.ys(k1),1,size(ghx,2)).*ghx.* ... - % repmat(dr.ys(k1(klag(:,1)))',size(ghx,1),1); - % ghu = repmat(1./dr.ys(k1),1,size(ghu,2)).*ghu; - end - - - if options_.aim_solver ~= 1 - % Necessary when using Sims' routines for QZ - ghx = real(ghx); - ghu = real(ghu); - end - - %exogenous deterministic variables - if exo_det_nbr > 0 - error('Deterministic exogenous variables are not yet implemented in first order approximation for a block decomposed model'); - % f1 = sparse(jacobia_(:,nonzeros(M_.lead_lag_incidence(M_.maximum_endo_lag+2:end,order_var)))); - % f0 = sparse(jacobia_(:,nonzeros(M_.lead_lag_incidence(M_.maximum_endo_lag+1,order_var)))); - % fudet = data(i).g1_xd; - % M1 = inv(f0+[zeros(n,n_static) f1*gx zeros(n,nyf-n_both)]); - % M2 = M1*f1; - % dr.ghud = cell(M_.exo_det_length,1); - % dr.ghud{1} = -M1*fudet; - % for i = 2:M_.exo_det_length - % dr.ghud{i} = -M2*dr.ghud{i-1}(end-nyf+1:end,:); - % end - end - end - end - if task ~=1 - if (maximum_lag > 0 && (n_pred > 0 || n_both > 0)) - sorted_col_dr_ghx = M_.block_structure.block(i).sorted_col_dr_ghx; - dr.ghx(endo, sorted_col_dr_ghx) = dr.ghx(endo, sorted_col_dr_ghx) + ghx; - data(i).ghx = ghx; - data(i).pol.i_ghx = sorted_col_dr_ghx; - else - data(i).pol.i_ghx = []; - end - data(i).ghu = ghu; - dr.ghu(endo, exo) = ghu; - data(i).pol.i_ghu = exo; - end - - if (verbose) - disp('dr.ghx'); - dr.ghx - disp('dr.ghu'); - dr.ghu - end - -end -M_.block_structure.block = data ; -if (verbose) - disp('dr.ghx'); - disp(real(dr.ghx)); - disp('dr.ghu'); - disp(real(dr.ghu)); -end -if (task == 1) - return -end diff --git a/matlab/dsge_likelihood.m b/matlab/dsge_likelihood.m index f5dec3529..7c5bbac7c 100644 --- a/matlab/dsge_likelihood.m +++ b/matlab/dsge_likelihood.m @@ -115,7 +115,7 @@ function [fval,info,exit_flag,DLIK,Hess,SteadyState,trend_coeff,Model,DynareOpti %! @end deftypefn %@eod: -% Copyright © 2004-2021 Dynare Team +% Copyright © 2004-2023 Dynare Team % % This file is part of Dynare. % @@ -680,9 +680,7 @@ singularity_has_been_detected = false; % First test multivariate filter if specified; potentially abort and use univariate filter instead if ((kalman_algo==1) || (kalman_algo==3))% Multivariate Kalman Filter if no_missing_data_flag && ~DynareOptions.occbin.likelihood.status - if DynareOptions.block - [LIK,lik] = block_kalman_filter(T,R,Q,H,Pstar,Y,start,Z,kalman_tol,riccati_tol, Model.nz_state_var, Model.n_diag, Model.nobs_non_statevar); - elseif DynareOptions.fast_kalman_filter + if DynareOptions.fast_kalman_filter if diffuse_periods %kalman_algo==3 requires no diffuse periods (stationary %observables) as otherwise FE matrix will not be positive @@ -709,23 +707,18 @@ if ((kalman_algo==1) || (kalman_algo==3))% Multivariate Kalman Filter analytic_deriv_info{:}); end else - if 0 %DynareOptions.block - [LIK,lik] = block_kalman_filter(DatasetInfo.missing.aindex,DatasetInfo.missing.number_of_observations,DatasetInfo.missing.no_more_missing_observations,... - T,R,Q,H,Pstar,Y,start,Z,kalman_tol,riccati_tol, Model.nz_state_var, Model.n_diag, Model.nobs_non_statevar); - else - [LIK,lik] = missing_observations_kalman_filter(DatasetInfo.missing.aindex,DatasetInfo.missing.number_of_observations,DatasetInfo.missing.no_more_missing_observations,Y,diffuse_periods+1,size(Y,2), ... - a_0_given_tm1, Pstar, ... - kalman_tol, DynareOptions.riccati_tol, ... - DynareOptions.rescale_prediction_error_covariance, ... - DynareOptions.presample, ... - T,Q,R,H,Z,mm,pp,rr,Zflag,diffuse_periods, occbin_); - if occbin_.status && isinf(LIK) - fval = Inf; - info(1) = 320; - info(4) = 0.1; - exit_flag = 0; - return - end + [LIK,lik] = missing_observations_kalman_filter(DatasetInfo.missing.aindex,DatasetInfo.missing.number_of_observations,DatasetInfo.missing.no_more_missing_observations,Y,diffuse_periods+1,size(Y,2), ... + a_0_given_tm1, Pstar, ... + kalman_tol, DynareOptions.riccati_tol, ... + DynareOptions.rescale_prediction_error_covariance, ... + DynareOptions.presample, ... + T,Q,R,H,Z,mm,pp,rr,Zflag,diffuse_periods, occbin_); + if occbin_.status && isinf(LIK) + fval = Inf; + info(1) = 320; + info(4) = 0.1; + exit_flag = 0; + return end end if analytic_derivation diff --git a/matlab/dynare_estimation_init.m b/matlab/dynare_estimation_init.m index b4b68d1da..2be5b45fd 100644 --- a/matlab/dynare_estimation_init.m +++ b/matlab/dynare_estimation_init.m @@ -32,7 +32,7 @@ function [dataset_, dataset_info, xparam1, hh, M_, options_, oo_, estim_params_, % SPECIAL REQUIREMENTS % none -% Copyright © 2003-2021 Dynare Team +% Copyright © 2003-2023 Dynare Team % % This file is part of Dynare. % @@ -120,14 +120,9 @@ if options_.analytic_derivation && options_.fast_kalman_filter end % fast kalman filter is only available with kalman_algo == 0,1,3 -if options_.fast_kalman_filter - if ~ismember(options_.kalman_algo, [0,1,3]) - error(['estimation option conflict: fast_kalman_filter is only available ' ... +if options_.fast_kalman_filter && ~ismember(options_.kalman_algo, [0,1,3]) + error(['estimation option conflict: fast_kalman_filter is only available ' ... 'with kalman_algo = 0, 1 or 3']) - elseif options_.block - error(['estimation option conflict: fast_kalman_filter is not available ' ... - 'with block']) - end end % Set options_.lik_init equal to 3 if diffuse filter is used or kalman_algo refers to a diffuse filter algorithm. @@ -456,45 +451,23 @@ else end % Define union of observed and state variables -if options_.block - k1 = k1'; - [k2, i_posA, i_posB] = union(k1', M_.state_var', 'rows'); - % Set restrict_state to postion of observed + state variables in expanded state vector. - oo_.dr.restrict_var_list = [k1(i_posA) M_.state_var(sort(i_posB))]; - % set mf0 to positions of state variables in restricted state vector for likelihood computation. - [~,bayestopt_.mf0] = ismember(M_.state_var',oo_.dr.restrict_var_list); - % Set mf1 to positions of observed variables in restricted state vector for likelihood computation. - [~,bayestopt_.mf1] = ismember(k1,oo_.dr.restrict_var_list); - % Set mf2 to positions of observed variables in expanded state vector for filtering and smoothing. - bayestopt_.mf2 = var_obs_index_dr; - bayestopt_.mfys = k1; - oo_.dr.restrict_columns = [size(i_posA,1)+(1:size(M_.state_var,2))]; - [k2, i_posA, i_posB] = union(k3p, M_.state_var', 'rows'); - bayestopt_.smoother_var_list = [k3p(i_posA); M_.state_var(sort(i_posB))']; - [~,~,bayestopt_.smoother_saved_var_list] = intersect(k3p,bayestopt_.smoother_var_list(:)); - [~,ic] = intersect(bayestopt_.smoother_var_list,M_.state_var); - bayestopt_.smoother_restrict_columns = ic; - [~,bayestopt_.smoother_mf] = ismember(k1, bayestopt_.smoother_var_list); -else - % Define union of observed and state variables - k2 = union(var_obs_index_dr,[M_.nstatic+1:M_.nstatic+M_.nspred]', 'rows'); - % Set restrict_state to postion of observed + state variables in expanded state vector. - oo_.dr.restrict_var_list = k2; - % set mf0 to positions of state variables in restricted state vector for likelihood computation. - [~,bayestopt_.mf0] = ismember([M_.nstatic+1:M_.nstatic+M_.nspred]',k2); - % Set mf1 to positions of observed variables in restricted state vector for likelihood computation. - [~,bayestopt_.mf1] = ismember(var_obs_index_dr,k2); - % Set mf2 to positions of observed variables in expanded state vector for filtering and smoothing. - bayestopt_.mf2 = var_obs_index_dr; - bayestopt_.mfys = k1; - [~,ic] = intersect(k2,nstatic+(1:npred)'); - oo_.dr.restrict_columns = [ic; length(k2)+(1:nspred-npred)']; - bayestopt_.smoother_var_list = union(k2,k3); - [~,~,bayestopt_.smoother_saved_var_list] = intersect(k3,bayestopt_.smoother_var_list(:)); - [~,ic] = intersect(bayestopt_.smoother_var_list,nstatic+(1:npred)'); - bayestopt_.smoother_restrict_columns = ic; - [~,bayestopt_.smoother_mf] = ismember(var_obs_index_dr, bayestopt_.smoother_var_list); -end +k2 = union(var_obs_index_dr,[M_.nstatic+1:M_.nstatic+M_.nspred]', 'rows'); +% Set restrict_state to postion of observed + state variables in expanded state vector. +oo_.dr.restrict_var_list = k2; +% set mf0 to positions of state variables in restricted state vector for likelihood computation. +[~,bayestopt_.mf0] = ismember([M_.nstatic+1:M_.nstatic+M_.nspred]',k2); +% Set mf1 to positions of observed variables in restricted state vector for likelihood computation. +[~,bayestopt_.mf1] = ismember(var_obs_index_dr,k2); +% Set mf2 to positions of observed variables in expanded state vector for filtering and smoothing. +bayestopt_.mf2 = var_obs_index_dr; +bayestopt_.mfys = k1; +[~,ic] = intersect(k2,nstatic+(1:npred)'); +oo_.dr.restrict_columns = [ic; length(k2)+(1:nspred-npred)']; +bayestopt_.smoother_var_list = union(k2,k3); +[~,~,bayestopt_.smoother_saved_var_list] = intersect(k3,bayestopt_.smoother_var_list(:)); +[~,ic] = intersect(bayestopt_.smoother_var_list,nstatic+(1:npred)'); +bayestopt_.smoother_restrict_columns = ic; +[~,bayestopt_.smoother_mf] = ismember(var_obs_index_dr, bayestopt_.smoother_var_list); if options_.analytic_derivation if options_.lik_init == 3 @@ -661,9 +634,6 @@ if options_.heteroskedastic_filter if options_.fast_kalman_filter error('estimation option conflict: "heteroskedastic_filter" incompatible with "fast_kalman_filter"') end - if options_.block - error('estimation option conflict: "heteroskedastic_filter" incompatible with block kalman filter') - end if options_.analytic_derivation error(['estimation option conflict: analytic_derivation isn''t available ' ... 'for heteroskedastic_filter']) @@ -717,4 +687,4 @@ if options_.occbin.smoother.status && options_.occbin.smoother.inversion_filter fprintf('dynare_estimation_init: the inversion filter does not support smoothed_state_uncertainty. Disabling the option.\n') options_.smoothed_state_uncertainty=false; end -end \ No newline at end of file +end diff --git a/matlab/dynare_resolve.m b/matlab/dynare_resolve.m index 727d9f38b..aca22c68a 100644 --- a/matlab/dynare_resolve.m +++ b/matlab/dynare_resolve.m @@ -16,7 +16,7 @@ function [A,B,ys,info,M_,oo_] = dynare_resolve(M_,options_,oo_,mode) % - M_ [structure] Matlab's structure describing the model % - oo_ [structure] Matlab's structure containing the results -% Copyright © 2001-2021 Dynare Team +% Copyright © 2001-2023 Dynare Team % % This file is part of Dynare. % @@ -52,11 +52,7 @@ switch nargin nstatic = M_.nstatic; nspred = M_.nspred; iv = (1:endo_nbr)'; - if ~options_.block - ic = [ nstatic+(1:nspred) endo_nbr+(1:size(oo_.dr.ghx,2)-nspred) ]'; - else - ic = oo_.dr.restrict_columns; - end + ic = [ nstatic+(1:nspred) endo_nbr+(1:size(oo_.dr.ghx,2)-nspred) ]'; case 4 iv = oo_.dr.restrict_var_list; ic = oo_.dr.restrict_columns; diff --git a/matlab/dynvars_from_endo_simul.m b/matlab/dynvars_from_endo_simul.m deleted file mode 100644 index ceddbefe0..000000000 --- a/matlab/dynvars_from_endo_simul.m +++ /dev/null @@ -1,26 +0,0 @@ -function y2 = dynvars_from_endo_simul(y, it_, M_) -% Given the matrix y of paths for all endogenous (same format as -% oo_.endo_simul), and an iteration number (first simulation period corresponds -% to it_=M_.maximum_lag+1), return a vector of endogenous values in the format -% expected by the dynamic.m file (i.e. whose indices are described by -% M_.lead_lag_incidence) - -% Copyright © 2020 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 . - -y2 = y(:, it_+(-M_.maximum_endo_lag:M_.maximum_endo_lead)); -y2 = y2(find(M_.lead_lag_incidence')); diff --git a/matlab/evaluate_steady_state.m b/matlab/evaluate_steady_state.m index eae8bf782..2dba4e288 100644 --- a/matlab/evaluate_steady_state.m +++ b/matlab/evaluate_steady_state.m @@ -69,11 +69,6 @@ if length(M.aux_vars) > 0 && ~steadystate_flag && M.set_auxiliary_variables end if options.ramsey_policy - if options.block - % The current implementation needs the Jacobian of the full model, which is not - % provided by the block-decomposed routines. - error('The ''block'' option is not compatible with ''ramsey_model''/''ramsey_policy'''); - end if ~isfinite(M.params(strmatch('optimal_policy_discount_factor',M.param_names,'exact'))) fprintf('\nevaluate_steady_state: the planner_discount is NaN/Inf. That will cause problems.\n') end diff --git a/matlab/model_diagnostics.m b/matlab/model_diagnostics.m index 3abbb06c0..8336767ef 100644 --- a/matlab/model_diagnostics.m +++ b/matlab/model_diagnostics.m @@ -256,51 +256,46 @@ iyr0 = find(iyv) ; it_ = M.maximum_lag + 1; z = repmat(dr.ys,1,klen); -if ~options.block - if options.order == 1 - if (options.bytecode) - [~, loc_dr] = bytecode('dynamic','evaluate', z,exo_simul, ... - M.params, dr.ys, 1); - jacobia_ = [loc_dr.g1 loc_dr.g1_x loc_dr.g1_xd]; - else - [~,jacobia_] = feval([M.fname '.dynamic'],z(iyr0),exo_simul, ... - M.params, dr.ys, it_); - end - elseif options.order >= 2 - if (options.bytecode) - [~, loc_dr] = bytecode('dynamic','evaluate', z,exo_simul, ... - M.params, dr.ys, 1); - jacobia_ = [loc_dr.g1 loc_dr.g1_x]; - else - [~,jacobia_,hessian1] = feval([M.fname '.dynamic'],z(iyr0),... - exo_simul, ... - M.params, dr.ys, it_); - end +if options.order == 1 + if (options.bytecode) + [~, loc_dr] = bytecode('dynamic','evaluate', z,exo_simul, ... + M.params, dr.ys, 1); + jacobia_ = [loc_dr.g1 loc_dr.g1_x loc_dr.g1_xd]; + else + [~,jacobia_] = feval([M.fname '.dynamic'],z(iyr0),exo_simul, ... + M.params, dr.ys, it_); end +elseif options.order >= 2 + if (options.bytecode) + [~, loc_dr] = bytecode('dynamic','evaluate', z,exo_simul, ... + M.params, dr.ys, 1); + jacobia_ = [loc_dr.g1 loc_dr.g1_x]; + else + [~,jacobia_,hessian1] = feval([M.fname '.dynamic'],z(iyr0),... + exo_simul, ... + M.params, dr.ys, it_); + end +end - if any(any(isinf(jacobia_) | isnan(jacobia_))) +if any(any(isinf(jacobia_) | isnan(jacobia_))) + problem_dummy=1; + [infrow,infcol]=find(isinf(jacobia_) | isnan(jacobia_)); + fprintf('\nMODEL_DIAGNOSTICS: The Jacobian of the dynamic model contains Inf or NaN. The problem arises from: \n\n') + display_problematic_vars_Jacobian(infrow,infcol,M,dr.ys,'dynamic','MODEL_DIAGNOSTICS: ') +end +if any(any(~isreal(jacobia_))) + [imagrow,imagcol]=find(abs(imag(jacobia_))>1e-15); + if ~isempty(imagrow) problem_dummy=1; - [infrow,infcol]=find(isinf(jacobia_) | isnan(jacobia_)); - fprintf('\nMODEL_DIAGNOSTICS: The Jacobian of the dynamic model contains Inf or NaN. The problem arises from: \n\n') - display_problematic_vars_Jacobian(infrow,infcol,M,dr.ys,'dynamic','MODEL_DIAGNOSTICS: ') + fprintf('\nMODEL_DIAGNOSTICS: The Jacobian of the dynamic model contains imaginary parts. The problem arises from: \n\n') + display_problematic_vars_Jacobian(imagrow,imagcol,M,dr.ys,'dynamic','MODEL_DIAGNOSTICS: ') end - if any(any(~isreal(jacobia_))) - [imagrow,imagcol]=find(abs(imag(jacobia_))>1e-15); - if ~isempty(imagrow) - problem_dummy=1; - fprintf('\nMODEL_DIAGNOSTICS: The Jacobian of the dynamic model contains imaginary parts. The problem arises from: \n\n') - display_problematic_vars_Jacobian(imagrow,imagcol,M,dr.ys,'dynamic','MODEL_DIAGNOSTICS: ') - end +end +if exist('hessian1','var') + if any(any(isinf(hessian1) | isnan(hessian1))) + problem_dummy=1; + fprintf('\nMODEL_DIAGNOSTICS: The Hessian of the dynamic model contains Inf or NaN.\n') end - if exist('hessian1','var') - if any(any(isinf(hessian1) | isnan(hessian1))) - problem_dummy=1; - fprintf('\nMODEL_DIAGNOSTICS: The Hessian of the dynamic model contains Inf or NaN.\n') - end - end -else - fprintf('\nMODEL_DIAGNOSTICS: This command currently does not support the block option for checking.\n') - fprintf('\nMODEL_DIAGNOSTICS: the dynamic model. You may want to disable it for doing model_diagnostics. Skipping this part.\n') end if problem_dummy==0 diff --git a/matlab/perfect-foresight-models/solve_block_decomposed_problem.m b/matlab/perfect-foresight-models/solve_block_decomposed_problem.m index 79a74aa5b..85989f652 100644 --- a/matlab/perfect-foresight-models/solve_block_decomposed_problem.m +++ b/matlab/perfect-foresight-models/solve_block_decomposed_problem.m @@ -87,7 +87,7 @@ for blk = 1:length(M_.block_structure.block) [y, T, oo_] = solve_one_boundary(fh_dynamic, y, oo_.exo_simul, M_.params, oo_.steady_state, T, y_index, M_.block_structure.block(blk).NNZDerivatives, options_.periods, M_.block_structure.block(blk).is_linear, blk, M_.maximum_lag, options_.simul.maxit, options_.dynatol.f, cutoff, options_.stack_solve_algo, is_forward, true, false, M_, options_, oo_); elseif M_.block_structure.block(blk).Simulation_Type == 5 || ... % solveTwoBoundariesSimple M_.block_structure.block(blk).Simulation_Type == 8 % solveTwoBoundariesComplete - [y, T, oo_] = solve_two_boundaries(fh_dynamic, y, oo_.exo_simul, M_.params, oo_.steady_state, T, y_index, M_.block_structure.block(blk).NNZDerivatives, options_.periods, M_.block_structure.block(blk).maximum_lag, M_.block_structure.block(blk).maximum_lead, M_.block_structure.block(blk).is_linear, blk, M_.maximum_lag, options_.simul.maxit, options_.dynatol.f, cutoff, options_.stack_solve_algo, options_, M_, oo_); + [y, T, oo_] = solve_two_boundaries(fh_dynamic, y, oo_.exo_simul, M_.params, oo_.steady_state, T, y_index, M_.block_structure.block(blk).NNZDerivatives, options_.periods, M_.block_structure.block(blk).is_linear, blk, M_.maximum_lag, options_.simul.maxit, options_.dynatol.f, cutoff, options_.stack_solve_algo, options_, M_, oo_); end tmp = y(M_.block_structure.block(blk).variable, :); diff --git a/matlab/resol.m b/matlab/resol.m index 805ac3d46..6c3fd30e5 100644 --- a/matlab/resol.m +++ b/matlab/resol.m @@ -32,7 +32,7 @@ function [dr, info, M, oo] = resol(check_flag, M, options, oo) % info(1)=24 -> M_.params has been updated in the steadystate routine and has some NaNs. % info(1)=30 -> Ergodic variance can't be computed. -% Copyright © 2001-2018 Dynare Team +% Copyright © 2001-2023 Dynare Team % % This file is part of Dynare. % @@ -111,11 +111,5 @@ if options.loglinear end end -if options.block - [dr,info,M,oo] = dr_block(dr,check_flag,M,options,oo); - dr.edim = nnz(abs(dr.eigval) > options.qz_criterium); - dr.sdim = dr.nd-dr.edim; -else - [dr,info] = stochastic_solvers(dr,check_flag,M,options,oo); -end +[dr, info] = stochastic_solvers(dr, check_flag, M, options, oo); oo.dr = dr; diff --git a/matlab/set_state_space.m b/matlab/set_state_space.m index f3fe17dae..e442a2c16 100644 --- a/matlab/set_state_space.m +++ b/matlab/set_state_space.m @@ -34,7 +34,7 @@ function dr=set_state_space(dr,DynareModel,DynareOptions) %! @end deftypefn %@eod: -% Copyright © 1996-2017 Dynare Team +% Copyright © 1996-2023 Dynare Team % % This file is part of Dynare. % @@ -69,11 +69,7 @@ else both_var = []; stat_var = setdiff([1:endo_nbr]',fwrd_var); end -if DynareOptions.block - order_var = DynareModel.block_structure.variable_reordered; -else - order_var = [ stat_var(:); pred_var(:); both_var(:); fwrd_var(:)]; -end +order_var = [ stat_var(:); pred_var(:); both_var(:); fwrd_var(:)]; inv_order_var(order_var) = (1:endo_nbr); % building kmask for z state vector in t+1 diff --git a/matlab/simult_.m b/matlab/simult_.m index df98de1f3..37ac55d2e 100644 --- a/matlab/simult_.m +++ b/matlab/simult_.m @@ -17,7 +17,7 @@ function y_=simult_(M_,options_,y0,dr,ex_,iorder) % SPECIAL REQUIREMENTS % none -% Copyright © 2001-2021 Dynare Team +% Copyright © 2001-2023 Dynare Team % % This file is part of Dynare. % @@ -63,19 +63,9 @@ if options_.k_order_solver && ~options_.pruning % Call dynare++ routines. y_(dr.order_var,:) = y_; y_=[y_start y_]; else - if options_.block - if M_.maximum_lag > 0 - k2 = dr.state_var; - else - k2 = []; - end - order_var = 1:endo_nbr; - dr.order_var = order_var; - else - k2 = dr.kstate(find(dr.kstate(:,2) <= M_.maximum_lag+1),[1 2]); - k2 = k2(:,1)+(M_.maximum_lag+1-k2(:,2))*endo_nbr; - order_var = dr.order_var; - end + k2 = dr.kstate(find(dr.kstate(:,2) <= M_.maximum_lag+1),[1 2]); + k2 = k2(:,1)+(M_.maximum_lag+1-k2(:,2))*endo_nbr; + order_var = dr.order_var; switch iorder case 1 diff --git a/matlab/solve_two_boundaries.m b/matlab/solve_two_boundaries.m index 946e34e6f..7ad7a4f40 100644 --- a/matlab/solve_two_boundaries.m +++ b/matlab/solve_two_boundaries.m @@ -1,4 +1,4 @@ -function [y, T, oo]= solve_two_boundaries(fh, y, x, params, steady_state, T, y_index, nze, periods, y_kmin_l, y_kmax_l, is_linear, Block_Num, y_kmin, maxit_, solve_tolf, cutoff, stack_solve_algo,options,M, oo) +function [y, T, oo]= solve_two_boundaries(fh, y, x, params, steady_state, T, y_index, nze, periods, is_linear, Block_Num, y_kmin, maxit_, solve_tolf, cutoff, stack_solve_algo,options,M, oo) % Computes the deterministic simulation of a block of equation containing % both lead and lag variables using relaxation methods % @@ -14,8 +14,6 @@ function [y, T, oo]= solve_two_boundaries(fh, y, x, params, steady_state, T, y_i % nze [integer] number of non-zero elements in the % jacobian matrix % periods [integer] number of simulation periods -% y_kmin_l [integer] maximum number of lag in the block -% y_kmax_l [integer] maximum number of lead in the block % is_linear [logical] Whether the block is linear % Block_Num [integer] block number % y_kmin [integer] maximum number of lag in the model diff --git a/matlab/th_autocovariances.m b/matlab/th_autocovariances.m index e18f66982..3a01d9be4 100644 --- a/matlab/th_autocovariances.m +++ b/matlab/th_autocovariances.m @@ -42,7 +42,7 @@ function [Gamma_y,stationary_vars] = th_autocovariances(dr,ivar,M_,options_,node % E(x_t) = (I - {g_x}\right)^{- 1} 0.5\left( g_{\sigma\sigma} \sigma^2 + g_{xx} Var(\hat x_t) + g_{uu} Var(u_t) \right) % \] % -% Copyright © 2001-2020 Dynare Team +% Copyright © 2001-2023 Dynare Team % % This file is part of Dynare. % @@ -92,42 +92,32 @@ nspred = M_.nspred; nstatic = M_.nstatic; nx = size(ghx,2); -if ~options_.block - %order_var = dr.order_var; - inv_order_var = dr.inv_order_var; - kstate = dr.kstate; - ikx = [nstatic+1:nstatic+nspred]; - k0 = kstate(find(kstate(:,2) <= M_.maximum_lag+1),:); - i0 = find(k0(:,2) == M_.maximum_lag+1); - i00 = i0; - n0 = length(i0); - AS = ghx(:,i0); - ghu1 = zeros(nx,M_.exo_nbr); - ghu1(i0,:) = ghu(ikx,:); - for i=M_.maximum_lag:-1:2 - i1 = find(k0(:,2) == i); - n1 = size(i1,1); - j1 = zeros(n1,1); - for k1 = 1:n1 - j1(k1) = find(k0(i00,1)==k0(i1(k1),1)); - end - AS(:,j1) = AS(:,j1)+ghx(:,i1); - i0 = i1; + +inv_order_var = dr.inv_order_var; +kstate = dr.kstate; +ikx = [nstatic+1:nstatic+nspred]; +k0 = kstate(find(kstate(:,2) <= M_.maximum_lag+1),:); +i0 = find(k0(:,2) == M_.maximum_lag+1); +i00 = i0; +n0 = length(i0); +AS = ghx(:,i0); +ghu1 = zeros(nx,M_.exo_nbr); +ghu1(i0,:) = ghu(ikx,:); +for i=M_.maximum_lag:-1:2 + i1 = find(k0(:,2) == i); + n1 = size(i1,1); + j1 = zeros(n1,1); + for k1 = 1:n1 + j1(k1) = find(k0(i00,1)==k0(i1(k1),1)); end -else - ghu1 = zeros(nx,M_.exo_nbr); - trend = 1:M_.endo_nbr; - inv_order_var = trend(M_.block_structure.variable_reordered); - ghu1(1:length(dr.state_var),:) = ghu(dr.state_var,:); + AS(:,j1) = AS(:,j1)+ghx(:,i1); + i0 = i1; end b = ghu1*M_.Sigma_e*ghu1'; -if ~options_.block - ipred = nstatic+(1:nspred)'; -else - ipred = dr.state_var; -end +ipred = nstatic+(1:nspred)'; + % state space representation for state variables only [A,B] = kalman_transition_matrix(dr,ipred,1:nx,M_.exo_nbr); % Compute stationary variables (before HP filtering), @@ -135,11 +125,7 @@ end % HP filtering, this mean correction is computed *before* filtering) if local_order == 2 || options_.hp_filter == 0 [vx, u] = lyapunov_symm(A,B*M_.Sigma_e*B',options_.lyapunov_fixed_point_tol,options_.qz_criterium,options_.lyapunov_complex_threshold,[],options_.debug); - if ~options_.block - iky = inv_order_var(ivar); - else - iky = ivar; - end + iky = inv_order_var(ivar); stationary_vars = (1:length(ivar))'; if ~isempty(u) x = abs(ghx*u); diff --git a/mex/build/block_kalman_filter.am b/mex/build/block_kalman_filter.am deleted file mode 100644 index eaf08bde2..000000000 --- a/mex/build/block_kalman_filter.am +++ /dev/null @@ -1,15 +0,0 @@ -mex_PROGRAMS = block_kalman_filter - -TOPDIR = $(top_srcdir)/../../sources/block_kalman_filter - -block_kalman_filter_CPPFLAGS = $(AM_CPPFLAGS) -I$(TOPDIR) -block_kalman_filter_CXXFLAGS = $(AM_CXXFLAGS) -fopenmp -block_kalman_filter_LDFLAGS = $(AM_LDFLAGS) $(OPENMP_LDFLAGS) - -nodist_block_kalman_filter_SOURCES = block_kalman_filter.cc - -BUILT_SOURCES = $(nodist_block_kalman_filter_SOURCES) -CLEANFILES = $(nodist_block_kalman_filter_SOURCES) - -%.cc: $(TOPDIR)/%.cc - $(LN_S) -f $< $@ diff --git a/mex/build/matlab/Makefile.am b/mex/build/matlab/Makefile.am index 9ffb4102b..54a7dc1f6 100644 --- a/mex/build/matlab/Makefile.am +++ b/mex/build/matlab/Makefile.am @@ -1,6 +1,6 @@ ACLOCAL_AMFLAGS = -I ../../../m4 -SUBDIRS = mjdgges kronecker bytecode block_kalman_filter sobol perfect_foresight_problem num_procs block_trust_region disclyap_fast libkordersim local_state_space_iterations folded_to_unfolded_dr k_order_simul k_order_mean cycle_reduction logarithmic_reduction riccati_update +SUBDIRS = mjdgges kronecker bytecode sobol perfect_foresight_problem num_procs block_trust_region disclyap_fast libkordersim local_state_space_iterations folded_to_unfolded_dr k_order_simul k_order_mean cycle_reduction logarithmic_reduction riccati_update # libdynare++ must come before gensylv and k_order_perturbation if ENABLE_MEX_DYNAREPLUSPLUS diff --git a/mex/build/matlab/block_kalman_filter/Makefile.am b/mex/build/matlab/block_kalman_filter/Makefile.am deleted file mode 100644 index 1c603a9d8..000000000 --- a/mex/build/matlab/block_kalman_filter/Makefile.am +++ /dev/null @@ -1,2 +0,0 @@ -include ../mex.am -include ../../block_kalman_filter.am diff --git a/mex/build/matlab/configure.ac b/mex/build/matlab/configure.ac index 901576c47..67fa9913e 100644 --- a/mex/build/matlab/configure.ac +++ b/mex/build/matlab/configure.ac @@ -157,7 +157,6 @@ AC_CONFIG_FILES([Makefile k_order_welfare/Makefile kalman_steady_state/Makefile ms_sbvar/Makefile - block_kalman_filter/Makefile sobol/Makefile local_state_space_iterations/Makefile libkordersim/Makefile diff --git a/mex/build/octave/Makefile.am b/mex/build/octave/Makefile.am index 9f936fb2b..c89c9de0b 100644 --- a/mex/build/octave/Makefile.am +++ b/mex/build/octave/Makefile.am @@ -1,6 +1,6 @@ ACLOCAL_AMFLAGS = -I ../../../m4 -SUBDIRS = mjdgges kronecker bytecode block_kalman_filter sobol perfect_foresight_problem num_procs block_trust_region disclyap_fast libkordersim local_state_space_iterations folded_to_unfolded_dr k_order_simul k_order_mean cycle_reduction logarithmic_reduction riccati_update +SUBDIRS = mjdgges kronecker bytecode sobol perfect_foresight_problem num_procs block_trust_region disclyap_fast libkordersim local_state_space_iterations folded_to_unfolded_dr k_order_simul k_order_mean cycle_reduction logarithmic_reduction riccati_update # libdynare++ must come before gensylv and k_order_perturbation if ENABLE_MEX_DYNAREPLUSPLUS diff --git a/mex/build/octave/block_kalman_filter/Makefile.am b/mex/build/octave/block_kalman_filter/Makefile.am deleted file mode 100644 index aa99d4233..000000000 --- a/mex/build/octave/block_kalman_filter/Makefile.am +++ /dev/null @@ -1,3 +0,0 @@ -EXEEXT = .mex -include ../mex.am -include ../../block_kalman_filter.am diff --git a/mex/build/octave/configure.ac b/mex/build/octave/configure.ac index 815f1decc..10a582e08 100644 --- a/mex/build/octave/configure.ac +++ b/mex/build/octave/configure.ac @@ -160,7 +160,6 @@ AC_CONFIG_FILES([Makefile k_order_welfare/Makefile kalman_steady_state/Makefile ms_sbvar/Makefile - block_kalman_filter/Makefile sobol/Makefile local_state_space_iterations/Makefile libkordersim/Makefile diff --git a/mex/sources/Makefile.am b/mex/sources/Makefile.am index fc179768e..dc9ec6abd 100644 --- a/mex/sources/Makefile.am +++ b/mex/sources/Makefile.am @@ -14,7 +14,6 @@ EXTRA_DIST = \ k_order_welfare \ kalman_steady_state \ ms-sbvar \ - block_kalman_filter \ sobol \ local_state_space_iterations \ libkordersim \ diff --git a/mex/sources/block_kalman_filter/block_kalman_filter.cc b/mex/sources/block_kalman_filter/block_kalman_filter.cc deleted file mode 100644 index 536d8c140..000000000 --- a/mex/sources/block_kalman_filter/block_kalman_filter.cc +++ /dev/null @@ -1,850 +0,0 @@ -/* - * Copyright © 2007-2022 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 . - */ - -#include -#include -#include - -#include "block_kalman_filter.hh" - -#define BLAS -//#define CUBLAS - -#ifdef CUBLAS -# include -# include -#endif - -void -mexDisp(const mxArray *P) -{ - size_t n = mxGetN(P); - size_t m = mxGetM(P); - const double *M = mxGetPr(P); - mexPrintf("%d x %d\n", m, n); - mexEvalString("drawnow;"); - for (size_t i = 0; i < m; i++) - { - for (size_t j = 0; j < n; j++) - mexPrintf(" %9.4f", M[i+ j * m]); - mexPrintf("\n"); - } - mexEvalString("drawnow;"); -} - -void -mexDisp(const double *M, int m, int n) -{ - mexPrintf("%d x %d\n", m, n); - mexEvalString("drawnow;"); - for (int i = 0; i < m; i++) - { - for (int j = 0; j < n; j++) - mexPrintf(" %9.4f", M[i+ j * m]); - mexPrintf("\n"); - } - mexEvalString("drawnow;"); -} - -/*if block - %nz_state_var = M_.nz_state_var; - while notsteady && t riccati_tol; - oldK = K(:); - end - end; - else - while notsteady && t riccati_tol; - oldK = K(:); - - end - end - end -*/ - -bool -not_all_abs_F_bellow_crit(const double *F, int size, double crit) -{ - int i = 0; - while (i < size && abs(F[i]) < crit) - i++; - - if (i < size) - return false; - else - return true; -} - -double -det(const double *F, int dim, const lapack_int *ipiv) -{ - double det = 1.0; - for (int i = 0; i < dim; i++) - if (ipiv[i] - 1 == i) - det *= F[i + i * dim]; - else - det *= -F[i + i * dim]; - return det; -} - -BlockKalmanFilter::BlockKalmanFilter(int nrhs, const mxArray *prhs[]) -{ - if (nrhs != 13 && nrhs != 16) - mexErrMsgTxt("block_kalman_filter requires exactly \n 13 input arguments for standard Kalman filter \nor\n 16 input arguments for missing observations Kalman filter."); - if (nrhs == 16) - missing_observations = true; - else - missing_observations = false; - if (missing_observations) - { - if (!mxIsCell(prhs[0])) - mexErrMsgTxt("the first input argument of block_missing_observations_kalman_filter must be a Cell Array."); - pdata_index = prhs[0]; - if (!mxIsDouble(prhs[1])) - mexErrMsgTxt("the second input argument of block_missing_observations_kalman_filter must be a scalar."); - number_of_observations = ceil(mxGetScalar(prhs[1])); - if (!mxIsDouble(prhs[2])) - mexErrMsgTxt("the third input argument of block_missing_observations_kalman_filter must be a scalar."); - no_more_missing_observations = ceil(mxGetScalar(prhs[2])); - pT = mxDuplicateArray(prhs[3]); - pR = mxDuplicateArray(prhs[4]); - pQ = mxDuplicateArray(prhs[5]); - pH = mxDuplicateArray(prhs[6]); - pP = mxDuplicateArray(prhs[7]); - pY = mxDuplicateArray(prhs[8]); - start = mxGetScalar(prhs[9]); - mfd = mxGetPr(prhs[10]); - kalman_tol = mxGetScalar(prhs[11]); - riccati_tol = mxGetScalar(prhs[12]); - nz_state_var = mxGetPr(prhs[13]); - n_diag = mxGetScalar(prhs[14]); - pure_obs = mxGetScalar(prhs[15]); - } - else - { - no_more_missing_observations = 0; - pT = mxDuplicateArray(prhs[0]); - pR = mxDuplicateArray(prhs[1]); - pQ = mxDuplicateArray(prhs[2]); - pH = mxDuplicateArray(prhs[3]); - pP = mxDuplicateArray(prhs[4]); - pY = mxDuplicateArray(prhs[5]); - start = mxGetScalar(prhs[6]); - /*Defining the initials values*/ - n = mxGetN(pT); // Number of state variables. - pp = mxGetM(pY); // Maximum number of observed variables. - smpl = mxGetN(pY); // Sample size. ; - mfd = mxGetPr(prhs[7]); - kalman_tol = mxGetScalar(prhs[8]); - riccati_tol = mxGetScalar(prhs[9]); - nz_state_var = mxGetPr(prhs[10]); - n_diag = mxGetScalar(prhs[11]); - pure_obs = mxGetScalar(prhs[12]); - } - T = mxGetPr(pT); - R = mxGetPr(pR); - Q = mxGetPr(pQ); - H = mxGetPr(pH); - P = mxGetPr(pP); - Y = mxGetPr(pY); - - n = mxGetN(pT); // Number of state variables. - pp = mxGetM(pY); // Maximum number of observed variables. - smpl = mxGetN(pY); // Sample size. ; - n_state = n - pure_obs; - - /*mexPrintf("T\n"); - mexDisp(pT);*/ - - H_size = mxGetN(pH) * mxGetM(pH); - - n_shocks = mxGetM(pQ); - - if (missing_observations) - if (mxGetNumberOfElements(pdata_index) != static_cast(smpl)) - mexErrMsgTxt("the number of element in the cell array passed to block_missing_observation_kalman_filter as first argument has to be equal to the smpl size"); - - i_nz_state_var = std::make_unique(n); - for (int i = 0; i < n; i++) - i_nz_state_var[i] = nz_state_var[i]; - - pa = mxCreateDoubleMatrix(n, 1, mxREAL); // State vector. - a = mxGetPr(pa); - tmp_a = std::make_unique(n); - dF = 0.0; // det(F). - - p_tmp1 = mxCreateDoubleMatrix(n, n_shocks, mxREAL); - tmp1 = mxGetPr(p_tmp1); - t = 0; // Initialization of the time index. - plik = mxCreateDoubleMatrix(smpl, 1, mxREAL); - lik = mxGetPr(plik); - Inf = mxGetInf(); - LIK = 0.0; // Default value of the log likelihood. - notsteady = true; // Steady state flag. - F_singular = true; - v_pp = std::make_unique(pp); - v_n = std::make_unique(n); - mf = std::make_unique(pp); - for (int i = 0; i < pp; i++) - mf[i] = mfd[i] - 1; - - /*compute QQ = R*Q*transpose(R)*/ // Variance of R times the vector of structural innovations.; - // tmp = R * Q; - for (int i = 0; i < n; i++) - for (int j = 0; j < n_shocks; j++) - { - double res = 0.0; - for (int k = 0; k < n_shocks; k++) - res += R[i + k * n] * Q[j * n_shocks + k]; - tmp1[i + j * n] = res; - } - - // QQ = tmp * transpose(R) - pQQ = mxCreateDoubleMatrix(n, n, mxREAL); - QQ = mxGetPr(pQQ); - for (int i = 0; i < n; i++) - for (int j = i; j < n; j++) - { - double res = 0.0; - for (int k = 0; k < n_shocks; k++) - res += tmp1[i + k * n] * R[k * n + j]; - QQ[i + j * n] = QQ[j + i * n] = res; - } - mxDestroyArray(p_tmp1); - - pv = mxCreateDoubleMatrix(pp, 1, mxREAL); - v = mxGetPr(pv); - pF = mxCreateDoubleMatrix(pp, pp, mxREAL); - F = mxGetPr(pF); - piF = mxCreateDoubleMatrix(pp, pp, mxREAL); - iF = mxGetPr(piF); - lw = pp * 4; - w = std::make_unique(lw); - iw = std::make_unique(pp); - ipiv = std::make_unique(pp); - info = 0; -#if defined(BLAS) || defined(CUBLAS) - p_tmp = mxCreateDoubleMatrix(n, n, mxREAL); - tmp = mxGetPr(p_tmp); - p_P_t_t1 = mxCreateDoubleMatrix(n, n, mxREAL); - P_t_t1 = mxGetPr(p_P_t_t1); - pK = mxCreateDoubleMatrix(n, n, mxREAL); - K = mxGetPr(pK); - p_K_P = mxCreateDoubleMatrix(n, n, mxREAL); - K_P = mxGetPr(p_K_P); - oldK = std::make_unique(n * n); - P_mf = std::make_unique(n * n); - for (int i = 0; i < n * n; i++) - oldK[i] = Inf; -#else - p_tmp = mxCreateDoubleMatrix(n, n_state, mxREAL); - tmp = mxGetPr(p_tmp); - p_P_t_t1 = mxCreateDoubleMatrix(n_state, n_state, mxREAL); - P_t_t1 = mxGetPr(p_P_t_t1); - pK = mxCreateDoubleMatrix(n, pp, mxREAL); - K = mxGetPr(pK); - p_K_P = mxCreateDoubleMatrix(n_state, n_state, mxREAL); - K_P = mxGetPr(p_K_P); - oldK = std::make_unique(n * pp); - P_mf = std::make_unique(n * pp); - for (int i = 0; i < n * pp; i++) - oldK[i] = Inf; -#endif -} - -void -BlockKalmanFilter::block_kalman_filter_ss() -{ - if (t+1 < smpl) - while (t < smpl) - { - //v = Y(:,t)-a(mf); - for (int i = 0; i < pp; i++) - v[i] = Y[i + t * pp] - a[mf[i]]; - - //a = T*(a+K*v); - for (int i = pure_obs; i < n; i++) - { - double res = 0.0; - for (int j = 0; j < pp; j++) - res += K[j * n + i] * v[j]; - v_n[i] = res + a[i]; - } - for (int i = 0; i < n; i++) - { - double res = 0.0; - for (int j = pure_obs; j < n; j++) - res += T[j * n + i] * v_n[j]; - a[i] = res; - } - - //lik(t) = transpose(v)*iF*v; - for (int i = 0; i < pp; i++) - { - double res = 0.0; - for (int j = 0; j < pp; j++) - res += v[j] * iF[j * pp + i]; - v_pp[i] = res; - } - double res = 0.0; - for (int i = 0; i < pp; i++) - res += v_pp[i] * v[i]; - - lik[t] = (log(dF) + res + pp * log(2.0*M_PI))/2; - if (t + 1 > start) - LIK += lik[t]; - - t++; - } -} - -bool -BlockKalmanFilter::block_kalman_filter(int nlhs, mxArray *plhs[]) -{ - while (notsteady && t < smpl) - { - if (missing_observations) - { - // retrieve the d_index - pd_index = mxGetCell(pdata_index, t); - dd_index = mxGetPr(pd_index); - size_d_index = mxGetM(pd_index); - d_index.resize(size_d_index); - for (int i = 0; i < size_d_index; i++) - d_index[i] = ceil(dd_index[i]) - 1; - - //v = Y(:,t) - a(mf) - int i_i = 0; - //#pragma omp parallel for shared(v, i_i, d_index) - for (auto i = d_index.begin(); i != d_index.end(); i++) - { - //mexPrintf("i_i=%d, omp_get_max_threads()=%d\n",i_i,omp_get_max_threads()); - v[i_i] = Y[*i + t * pp] - a[mf[*i]]; - i_i++; - } - - //F = P(mf,mf) + H; - i_i = 0; - if (H_size == 1) - //#pragma omp parallel for shared(iF, F, i_i) - for (auto i = d_index.begin(); i != d_index.end(); i++, i_i++) - { - int j_j = 0; - for (auto j = d_index.begin(); j != d_index.end(); j++, j_j++) - iF[i_i + j_j * size_d_index] = F[i_i + j_j * size_d_index] = P[mf[*i] + mf[*j] * n] + H[0]; - } - else - //#pragma omp parallel for shared(iF, F, P, H, mf, i_i) - for (auto i = d_index.begin(); i != d_index.end(); i++, i_i++) - { - int j_j = 0; - for (auto j = d_index.begin(); j != d_index.end(); j++, j_j++) - iF[i_i + j_j * size_d_index] = F[i_i + j_j * size_d_index] = P[mf[*i] + mf[*j] * n] + H[*i + *j * pp]; - } - } - else - { - size_d_index = pp; - - //v = Y(:,t) - a(mf) - for (int i = 0; i < pp; i++) - v[i] = Y[i + t * pp] - a[mf[i]]; - - //F = P(mf,mf) + H; - if (H_size == 1) - for (int i = 0; i < pp; i++) - for (int j = 0; j < pp; j++) - iF[i + j * pp] = F[i + j * pp] = P[mf[i] + mf[j] * n] + H[0]; - else - for (int i = 0; i < pp; i++) - for (int j = 0; j < pp; j++) - iF[i + j * pp] = F[i + j * pp] = P[mf[i] + mf[j] * n] + H[i + j * pp]; - } - - /* Computes the norm of iF */ - double anorm = dlange("1", &size_d_index, &size_d_index, iF, &size_d_index, w.get()); - //mexPrintf("anorm = %f\n",anorm); - - /* Modifies F in place with a LU decomposition */ - dgetrf(&size_d_index, &size_d_index, iF, &size_d_index, ipiv.get(), &info); - if (info != 0) - mexPrintf("dgetrf failure with error %d\n", static_cast(info)); - - /* Computes the reciprocal norm */ - dgecon("1", &size_d_index, iF, &size_d_index, &anorm, &rcond, w.get(), iw.get(), &info); - if (info != 0) - mexPrintf("dgecon failure with error %d\n", static_cast(info)); - - if (rcond < kalman_tol) - if (not_all_abs_F_bellow_crit(F, size_d_index * size_d_index, kalman_tol)) //~all(abs(F(:))(info)); - - //lik(t) = log(dF)+transpose(v)*iF*v; -#pragma omp parallel for shared(v_pp) - for (int i = 0; i < size_d_index; i++) - { - double res = 0.0; - for (int j = 0; j < size_d_index; j++) - res += v[j] * iF[j * size_d_index + i]; - v_pp[i] = res; - } - double res = 0.0; - for (int i = 0; i < size_d_index; i++) - res += v_pp[i] * v[i]; - - lik[t] = (log(dF) + res + size_d_index * log(2.0*M_PI))/2; - if (t + 1 >= start) - LIK += lik[t]; - - if (missing_observations) - //K = P(:,mf)*iF; -#pragma omp parallel for shared(P_mf) - for (int i = 0; i < n; i++) - { - int j_j = 0; - //for (int j = 0; j < pp; j++) - for (auto j = d_index.begin(); j != d_index.end(); j++, j_j++) - P_mf[i + j_j * n] = P[i + mf[*j] * n]; - } - else - //K = P(:,mf)*iF; - for (int i = 0; i < n; i++) - for (int j = 0; j < pp; j++) - P_mf[i + j * n] = P[i + mf[j] * n]; - -#pragma omp parallel for shared(K) - for (int i = 0; i < n; i++) - for (int j = 0; j < size_d_index; j++) - { - double res = 0.0; - int j_pp = j * size_d_index; - for (int k = 0; k < size_d_index; k++) - res += P_mf[i + k * n] * iF[j_pp + k]; - K[i + j * n] = res; - } - - //a = T*(a+K*v); -#pragma omp parallel for shared(v_n) - for (int i = pure_obs; i < n; i++) - { - double res = 0.0; - for (int j = 0; j < size_d_index; j++) - res += K[j * n + i] * v[j]; - v_n[i] = res + a[i]; - } - -#pragma omp parallel for shared(a) - for (int i = 0; i < n; i++) - { - double res = 0.0; - for (int j = pure_obs; j < n; j++) - res += T[j * n + i] * v_n[j]; - a[i] = res; - } - - if (missing_observations) - { - //P = T*(P-K*P(mf,:))*transpose(T)+QQ; - int i_i = 0; - //#pragma omp parallel for shared(P_mf) - for (auto i = d_index.begin(); i != d_index.end(); i++, i_i++) - for (int j = pure_obs; j < n; j++) - P_mf[i_i + j * size_d_index] = P[mf[*i] + j * n]; - } - else - //P = T*(P-K*P(mf,:))*transpose(T)+QQ; -#pragma omp parallel for shared(P_mf) - for (int i = 0; i < pp; i++) - for (int j = pure_obs; j < n; j++) - P_mf[i + j * pp] = P[mf[i] + j * n]; - -#ifdef BLAS -# pragma omp parallel for shared(K_P) - for (int i = 0; i < n; i++) - for (int j = i; j < n; j++) - { - double res = 0.0; - //int j_pp = j * pp; - for (int k = 0; k < size_d_index; k++) - res += K[i + k * n] * P_mf[k + j * size_d_index]; - K_P[i * n + j] = K_P[j * n + i] = res; - } - //#pragma omp parallel for shared(P, K_P, P_t_t1) - for (int i = size_d_index; i < n; i++) - for (int j = i; j < n; j++) - { - unsigned int k = i * n + j; - P_t_t1[j * n + i] = P_t_t1[k] = P[k] - K_P[k]; - } - double one = 1.0; - double zero = 0.0; - std::copy_n(QQ, n * n, P); - blas_int n_b = n; - /*mexPrintf("sizeof(n_b)=%d, n_b=%d, sizeof(n)=%d, n=%d\n",sizeof(n_b),n_b,sizeof(n),n); - mexEvalString("drawnow;");*/ - dsymm("R", "U", &n_b, &n_b, - &one, P_t_t1, &n_b, - T, &n_b, &zero, - tmp, &n_b); - dgemm("N", "T", &n_b, &n_b, - &n_b, &one, tmp, &n_b, - T, &n_b, &one, - P, &n_b); -#else -# ifdef CUBLAS - for (int i = 0; i < n; i++) - for (int j = i; j < n; j++) - { - double res = 0.0; - //int j_pp = j * pp; - for (int k = 0; k < size_d_index; k++) - res += K[i + k * n] * P_mf[k + j * size_d_index]; - K_P[i * n + j] = K_P[j * n + i] = res; - } - //#pragma omp parallel for shared(P, K_P, P_t_t1) - for (int i = size_d_index; i < n; i++) - for (int j = i; j < n; j++) - { - unsigned int k = i * n + j; - P_t_t1[j * n + i] = P_t_t1[k] = P[k] - K_P[k]; - } - mexPrintf("CudaBLAS\n"); - mexEvalString("drawnow;"); - double one = 1.0; - double zero = 0.0; - cublasStatus_t status; - cublasHandle_t handle; - status = cublasCreate(&handle); - if (status != CUBLAS_STATUS_SUCCESS) - { - mexPrintf("!!!! CUBLAS initialization error\n"); - return false; - } - /*int device; - cudaGetDevice(&device);*/ - int n2 = n * n; - double *d_A = nullptr, *d_B = nullptr, *d_C = nullptr, *d_D = nullptr; - // Allocate device memory for the matrices - if (cudaMalloc(static_cast(&d_A), n2 * sizeof(double)) != cudaSuccess) - { - mexPrintf("!!!! device memory allocation error (allocate A)\n"); - return false; - } - if (cudaMalloc(static_cast(&d_B), n2 * sizeof(d_B[0])) != cudaSuccess) - { - mexPrintf("!!!! device memory allocation error (allocate B)\n"); - return false; - } - if (cudaMalloc(static_cast(&d_C), n2 * sizeof(d_C[0])) != cudaSuccess) - { - mexPrintf("!!!! device memory allocation error (allocate C)\n"); - return false; - } - if (cudaMalloc(static_cast(&d_D), n2 * sizeof(d_D[0])) != cudaSuccess) - { - mexPrintf("!!!! device memory allocation error (allocate D)\n"); - return false; - } - // Initialize the device matrices with the host matrices - status = cublasSetVector(n2, sizeof(P_t_t1[0]), P_t_t1, 1, d_A, 1); - if (status != CUBLAS_STATUS_SUCCESS) - { - mexPrintf("!!!! device access error (write A)\n"); - return false; - } - status = cublasSetVector(n2, sizeof(T[0]), T, 1, d_B, 1); - if (status != CUBLAS_STATUS_SUCCESS) - { - mexPrintf("!!!! device access error (write B)\n"); - return false; - } - status = cublasSetVector(n2, sizeof(tmp[0]), tmp, 1, d_C, 1); - if (status != CUBLAS_STATUS_SUCCESS) - { - mexPrintf("!!!! device access error (write C)\n"); - return false; - } - mexPrintf("just before calling\n"); - mexEvalString("drawnow;"); - status = cublasSetVector(n2, sizeof(QQ[0]), QQ, 1, d_D, 1); - if (status != CUBLAS_STATUS_SUCCESS) - { - mexPrintf("!!!! device access error (write D)\n"); - return false; - } - - // Performs operation using plain C code - - cublasDsymm(handle, CUBLAS_SIDE_RIGHT, CUBLAS_FILL_MODE_UPPER, n, n, - &one, d_A, n, - d_B, n, &zero, - d_C, n); - /*dgemm("N", "T", &n_b, &n_b, - &n_b, &one, tmp, &n_b, - T, &n_b, &one, - P, &n_b);*/ - cublasDgemm(handle, CUBLAS_OP_N, CUBLAS_OP_T, n, n, - n, &one, d_C, n, - d_B, n, &one, - d_D, n); - //double_symm(n, &one, h_A, h_B, &zero, h_C); - - status = cublasGetVector(n2, sizeof(P[0]), d_D, 1, P, 1); - if (status != CUBLAS_STATUS_SUCCESS) - { - mexPrintf("!!!! device access error (read P)\n"); - return false; - } - -# else -# pragma omp parallel for shared(K_P) - for (int i = pure_obs; i < n; i++) - { - unsigned int i1 = i - pure_obs; - for (int j = i; j < n; j++) - { - unsigned int j1 = j - pure_obs; - double res = 0.0; - int j_pp = j * size_d_index; - for (int k = 0; k < size_d_index; k++) - res += K[i + k * n] * P_mf[k + j_pp]; - K_P[i1 * n_state + j1] = K_P[j1 * n_state + i1] = res; - } - } - -# pragma omp parallel for shared(P_t_t1) - for (int i = pure_obs; i < n; i++) - { - unsigned int i1 = i - pure_obs; - for (int j = i; j < n; j++) - { - unsigned int j1 = j - pure_obs; - unsigned int k1 = i1 * n_state + j1; - P_t_t1[j1 * n_state + i1] = P_t_t1[k1] = P[i * n + j] - K_P[k1]; - } - } - - fill_n(tmp, 0, n * n_state); - -# pragma omp parallel for shared(tmp) - for (int i = 0; i < n; i++) - { - int max_k = i_nz_state_var[i]; - for (int j = pure_obs; j < n; j++) - { - int j1 = j - pure_obs; - int j1_n_state = j1 * n_state - pure_obs; - int indx_tmp = i + j1 * n; - for (int k = pure_obs; k < max_k; k++) - tmp[indx_tmp] += T[i + k * n] * P_t_t1[k + j1_n_state]; - } - } - - fill_n(P, 0, n * n); - - int n_n_obs = -n * pure_obs; -# pragma omp parallel for shared(P) - for (int i = 0; i < n; i++) - { - for (int j = i; j < n; j++) - { - int max_k = i_nz_state_var[j]; - int P_indx = i * n + j; - for (int k = pure_obs; k < max_k; k++) - { - int k_n = k * n; - P[P_indx] += tmp[i + k_n + n_n_obs] * T[j + k_n]; - } - } - } - -# pragma omp parallel for shared(P) - for (int i = 0; i < n; i++) - { - for (int j = i; j < n; j++) - P[j + i * n] += QQ[j + i * n]; - for (int j = i + 1; j < n; j++) - P[i + j * n] = P[j + i * n]; - } -# endif -#endif - if (t >= no_more_missing_observations) - { - double max_abs = 0.0; - for (int i = 0; i < n * size_d_index; i++) - { - double res = abs(K[i] - oldK[i]); - max_abs = std::max(res, max_abs); - } - notsteady = max_abs > riccati_tol; - - //oldK = K(:); - - std::copy_n(K, n * pp, oldK.get()); - } - } - t++; - } - - if (F_singular) - mexErrMsgTxt("The variance of the forecast error remains singular until the end of the sample\n"); - if (t < smpl) - block_kalman_filter_ss(); - return true; -} - -void -BlockKalmanFilter::return_results_and_clean(int nlhs, mxArray *plhs[]) -{ - if (nlhs > 2) - mexErrMsgTxt("block_kalman_filter provides at most 2 output argument."); - - if (nlhs >= 1) - { - plhs[0] = mxCreateDoubleMatrix(1, 1, mxREAL); - double *pind = mxGetPr(plhs[0]); - pind[0] = LIK; - } - - if (nlhs == 2) - plhs[1] = plik; - else - mxDestroyArray(plik); - - mxDestroyArray(pa); - mxDestroyArray(p_tmp); - mxDestroyArray(pQQ); - mxDestroyArray(pv); - mxDestroyArray(pF); - mxDestroyArray(piF); - mxDestroyArray(p_P_t_t1); - mxDestroyArray(pK); - mxDestroyArray(p_K_P); -} - -void -mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) -{ - BlockKalmanFilter block_kalman_filter(nrhs, prhs); - if (block_kalman_filter.block_kalman_filter(nlhs, plhs)) - block_kalman_filter.return_results_and_clean(nlhs, plhs); -} diff --git a/mex/sources/block_kalman_filter/block_kalman_filter.hh b/mex/sources/block_kalman_filter/block_kalman_filter.hh deleted file mode 100644 index 8687e92a7..000000000 --- a/mex/sources/block_kalman_filter/block_kalman_filter.hh +++ /dev/null @@ -1,68 +0,0 @@ -/* - * Copyright © 2007-2022 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 . - */ - -#ifndef BLOCK_KALMAN_FILTER -#define BLOCK_KALMAN_FILTER - -#ifndef DEBUG_EX -# include -#else -# include "mex_interface.hh" -#endif - -#include -#include - -#include -#include - -class BlockKalmanFilter -{ -public: - mxArray *pT, *pR, *pQ, *pH, *pP, *pY, *pQQ, *pv, *pF, *piF, *p_P_t_t1, *pK, *p_K_P; - double *T, *R, *Q, *H, *Y, *mfd, *QQ, *v, *F, *iF; - int start, pure_obs, smpl, n, n_state, n_shocks, H_size; - double kalman_tol, riccati_tol, dF, LIK, Inf, pi; - lapack_int pp, lw, info; - - double *nz_state_var; - std::unique_ptr i_nz_state_var, mf; - int n_diag, t; - mxArray *M_; - mxArray *pa, *p_tmp, *p_tmp1, *plik; - std::unique_ptr tmp_a; - double *tmp1, *lik; - std::unique_ptr v_n, v_pp, w, oldK, P_mf; - bool notsteady, F_singular, missing_observations; - std::unique_ptr iw, ipiv; - double anorm, rcond; - lapack_int size_d_index; - int no_more_missing_observations, number_of_observations; - const mxArray *pdata_index; - std::vector d_index; - const mxArray *pd_index; - double *dd_index; - double *K, *a, *K_P, *P_t_t1, *tmp, *P; -public: - BlockKalmanFilter(int nrhs, const mxArray *prhs[]); - bool block_kalman_filter(int nlhs, mxArray *plhs[]); - void block_kalman_filter_ss(); - void return_results_and_clean(int nlhs, mxArray *plhs[]); -}; -#endif diff --git a/preprocessor b/preprocessor index 35ac73fad..93b9ed695 160000 --- a/preprocessor +++ b/preprocessor @@ -1 +1 @@ -Subproject commit 35ac73fad8645302c10d0cd0068710ad3847366a +Subproject commit 93b9ed69577246cb13eef1890d52708f2aea4a9d diff --git a/tests/Makefile.am b/tests/Makefile.am index 478a2e4a4..b6ed5b9ae 100644 --- a/tests/Makefile.am +++ b/tests/Makefile.am @@ -198,8 +198,6 @@ MODFILES = \ block_bytecode/lola_solve_one_boundary_mfs1.mod \ block_bytecode/lola_solve_one_boundary_mfs2.mod \ block_bytecode/lola_solve_one_boundary_mfs3.mod \ - block_bytecode/lola_stochastic.mod \ - block_bytecode/lola_stochastic_block.mod \ k_order_perturbation/fs2000k2a.mod \ k_order_perturbation/fs2000k2_use_dll.mod \ k_order_perturbation/fs2000k_1_use_dll.mod \ @@ -837,8 +835,6 @@ block_bytecode/lola_solve_one_boundary_mfs2.m.trs: block_bytecode/lola_solve_one block_bytecode/lola_solve_one_boundary_mfs2.o.trs: block_bytecode/lola_solve_one_boundary.o.trs block_bytecode/lola_solve_one_boundary_mfs3.m.trs: block_bytecode/lola_solve_one_boundary.m.trs block_bytecode/lola_solve_one_boundary_mfs3.o.trs: block_bytecode/lola_solve_one_boundary.o.trs -block_bytecode/lola_stochastic_block.m.trs: block_bytecode/lola_stochastic.m.trs -block_bytecode/lola_stochastic_block.o.trs: block_bytecode/lola_stochastic.o.trs histval_initval_file/ramst_initval_file.m.trs: histval_initval_file/ramst_initval_file_data.m.tls histval_initval_file/ramst_data_generate.m.trs histval_initval_file/ramst_initval_file.o.trs: histval_initval_file/ramst_initval_file_data.o.tls histval_initval_file/ramst_data_generate.o.trs diff --git a/tests/block_bytecode/lola_common.inc b/tests/block_bytecode/lola_common.inc index ae582cea5..2c3bab87e 100644 --- a/tests/block_bytecode/lola_common.inc +++ b/tests/block_bytecode/lola_common.inc @@ -6,7 +6,6 @@ Macro-variables that can be defined to tune the computations: - block - mfs - - deterministic */ load lola_data.mat @@ -723,8 +722,6 @@ end; resid; steady(solve_algo=3); -@#if deterministic - endval; @#for i in wg n@{i}=n@{i}_fss; @@ -983,7 +980,3 @@ perfect_foresight_solver(maxit=100); if ~oo_.deterministic_simulation.status error('Perfect foresight simulation failed') end - -@#else // stochastic case, used by files under tests/block_bytecode/lola_* -stoch_simul(order=1); -@#endif diff --git a/tests/block_bytecode/lola_solve_one_boundary.mod b/tests/block_bytecode/lola_solve_one_boundary.mod index f72b8afca..01378dc2d 100644 --- a/tests/block_bytecode/lola_solve_one_boundary.mod +++ b/tests/block_bytecode/lola_solve_one_boundary.mod @@ -1,6 +1,5 @@ // Involves a call to solve_one_boundary.m that is tested here -@#define deterministic = true @#define block = true @#define mfs = 0 @#include "lola_common.inc" diff --git a/tests/block_bytecode/lola_solve_one_boundary_mfs1.mod b/tests/block_bytecode/lola_solve_one_boundary_mfs1.mod index 64ff378bf..8cd8a2e5c 100644 --- a/tests/block_bytecode/lola_solve_one_boundary_mfs1.mod +++ b/tests/block_bytecode/lola_solve_one_boundary_mfs1.mod @@ -1,6 +1,5 @@ // Tests option mfs=1 with block -@#define deterministic = true @#define block = true @#define mfs = 1 @#include "lola_common.inc" diff --git a/tests/block_bytecode/lola_solve_one_boundary_mfs2.mod b/tests/block_bytecode/lola_solve_one_boundary_mfs2.mod index 92a9c0f07..7a76efb28 100644 --- a/tests/block_bytecode/lola_solve_one_boundary_mfs2.mod +++ b/tests/block_bytecode/lola_solve_one_boundary_mfs2.mod @@ -1,6 +1,5 @@ // Tests option mfs=2 with block -@#define deterministic = true @#define block = true @#define mfs = 2 @#include "lola_common.inc" diff --git a/tests/block_bytecode/lola_solve_one_boundary_mfs3.mod b/tests/block_bytecode/lola_solve_one_boundary_mfs3.mod index 16b227463..fa549b68b 100644 --- a/tests/block_bytecode/lola_solve_one_boundary_mfs3.mod +++ b/tests/block_bytecode/lola_solve_one_boundary_mfs3.mod @@ -1,6 +1,5 @@ // Tests option mfs=3 with block -@#define deterministic = true @#define block = true @#define mfs = 3 @#include "lola_common.inc" diff --git a/tests/block_bytecode/lola_stochastic.mod b/tests/block_bytecode/lola_stochastic.mod deleted file mode 100644 index e6c0071fc..000000000 --- a/tests/block_bytecode/lola_stochastic.mod +++ /dev/null @@ -1,5 +0,0 @@ -// Stochastic version of LOLA model, as benchmark for lola_stochastic_block.mod - -@#define deterministic = false -@#define block = false -@#include "lola_common.inc" diff --git a/tests/block_bytecode/lola_stochastic_block.mod b/tests/block_bytecode/lola_stochastic_block.mod deleted file mode 100644 index c4db5454f..000000000 --- a/tests/block_bytecode/lola_stochastic_block.mod +++ /dev/null @@ -1,26 +0,0 @@ -/* Stochastic version of block decomposed LOLA model. - Check that policy functions are the same as in non-block version. */ - -@#define deterministic = false -@#define block = true -@#define mfs = 0 -@#include "lola_common.inc" - -[~, state_reorder] = sort(oo_.dr.state_var); - -ref = load(['lola_stochastic' filesep 'Output' filesep 'lola_stochastic_results.mat']); - -[~, ref_state_reorder] = sort(ref.oo_.dr.state_var); - -/* NB: With block, the rows of ghx and ghu are in declaration order (and not in - DR-order as in non-block mode) */ - -if max(max(abs(oo_.dr.ghx(:, state_reorder) - ref.oo_.dr.ghx(ref.oo_.dr.inv_order_var, ref_state_reorder)))) > 3e-9 - error('Error in ghx') -end - -if max(max(abs(oo_.dr.ghu - ref.oo_.dr.ghu(ref.oo_.dr.inv_order_var, :)))) > 5e-8 - error('Error in ghu') -end - -