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

Merge remote-tracking branch 'github/master' into ecb-master 

Fixed conflicts:
	matlab/modules/dates
time-shift
Stéphane Adjemian(Charybdis) 2018-05-16 17:32:39 +02:00
parent eee499faf7 ac5fcbe231
commit e62c2272b4
51 changed files with 901 additions and 787 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

92
doc/dynare.texi
View File

@ -222,6 +222,10 @@ The Model file
* Verbatim inclusion::
* Misc commands::
Variable declarations
* On-the-fly Model Variable Declaration::
Expressions
* Parameters and variables::
@ -382,7 +386,7 @@ Adjemian (Université du Maine, Gains and Cepremap), Houtan Bastani
(Université du Maine, Gains and Cepremap), Junior Maih (Norges Bank),
Ferhat Mihoubi (Université Paris-Est Créteil, Epee and Cepremap), George
Perendia, Johannes Pfeifer (University of Cologne), Marco Ratto (European Commission, Joint Research Centre - JRC)
and Sébastien Villemot (OFCE – Sciences Po).
and Sébastien Villemot (CEPREMAP).
Increasingly, the developer base is expanding, as tools developed by
researchers outside of Cepremap are integrated into Dynare. Financial
support is provided by Cepremap, Banque de France and DSGE-net (an
@ -985,6 +989,7 @@ Allows Dynare to issue a warning and continue processing when
@enumerate
@item there are more endogenous variables than equations
@item an undeclared symbol is assigned in @code{initval} or @code{endval}
@item an undeclared symbol is found in the @code{model} block; in this case, it is automatically declared exogenous
@item exogenous variables were declared but not used in the @code{model} block
@end enumerate
@ -1621,6 +1626,50 @@ model_local_variable GDP_US $GDPUS$;
@end deffn
@menu
* On-the-fly Model Variable Declaration::
@end menu
@node On-the-fly Model Variable Declaration
@subsection On-the-fly Model Variable Declaration
Endogenous variables, exogenous variables, and parameters can also be declared
inside the model block. To do this, simply follow the symbol name with a
vertical line (@code{|}) and either an @code{e}, an @code{x}, or a
@code{p}. For example, to declare a parameter named @code{alphaa} in the model
block, you could write @code{alphaa|p} directly in an equation where it
appears. Similarly, to declare an endogenous variable @code{c} in the model
block you could write @code{c|e}. These on-the-fly variable declarations do not
have to appear in the first place where this variable is encountered. Note that
on-the-fly variable declarations must be made on contemporaneous variables.
As an example, the following two snippets are equivalent:
@emph{Using on-the-fly variable and parameter declaration}
@example
model;
k(+1) = i|e + (1-delta|p)*k;
y|e = k|e^alpha|p;
@dots{}
end;
delta = 0.025;
alpha = 0.36;
@end example
@emph{Using standard variable declaration}
@example
var k, i, y;
parameters delta, alpha;
delta = 0.025;
alpha = 0.36;
@dots{}
model;
k(1) = i|e + (1-delta|p)*k;
y|e = k|e^alpha|p;
@dots{}
end;
@end example
@node Expressions
@section Expressions
@ -2337,7 +2386,7 @@ For the required @LaTeX{} packages, @pxref{write_latex_original_model}.
@table @code
@item write_equation_tags
@xref{write_equation_tags}
@xref{write_equation_tags}.
@end table
@ -2378,7 +2427,7 @@ For the required @LaTeX{} packages, @pxref{write_latex_original_model}.
@table @code
@item write_equation_tags
@xref{write_equation_tags}
@xref{write_equation_tags}.
@end table
@ -3732,6 +3781,9 @@ the square submatrix of the right Schur vectors corresponding to the
forward looking variables (jumpers) and to the explosive eigenvalues
must have full rank.
Note that the outcome may be different from what would be suggested by
@code{sum(abs(oo_.dr.eigval))} when eigenvalues are very close to @ref{qz_criterium}.
@optionshead
@table @code
@ -4005,7 +4057,7 @@ options).
@item 7
Allows the user to solve the perfect foresight model with the solvers available
through option @code{solve_algo} (@xref{solve_algo} for a list of possible
through option @code{solve_algo} (@pxref{solve_algo} for a list of possible
values, note that values 5, 6, 7 and 8, which require @code{bytecode} and/or
@code{block} options, are not allowed). For instance, the following commands:
@example
@ -4351,6 +4403,7 @@ exogenous variables are made available in @code{oo_.exo_simul}
(@pxref{oo_.exo_simul}). Default: @code{0}.
@item qz_criterium = @var{DOUBLE}
@anchor{qz_criterium}
Value used to split stable from unstable eigenvalues in reordering the
Generalized Schur decomposition used for solving 1^st order
problems. Default: @code{1.000001} (except when estimating with
@ -5612,12 +5665,17 @@ variance if the acceptance ratio is too low. In some situations it may
help to consider parameter-specific values for this scale parameter.
This can be done in the @ref{estimated_params}- block.
Note that @code{mode_compute=6} will tune the scale parameter to achieve an acceptance rate of
@ref{AcceptanceRateTarget}. The resulting scale parameter will be saved into a file
named @file{@var{MODEL_FILENAME}_mh_scale.mat}. This file can be loaded in subsequent runs
via the @code{posterior_sampler_options}-option @ref{scale_file}. Both @code{mode_compute=6}
and @code{scale_file} will overwrite any value specified in @code{estimated_params} with the tuned value.
Default: @code{0.2}
Note that @code{mode_compute=6} will tune the scale parameter to achieve an
acceptance rate of @ref{AcceptanceRateTarget}. The resulting scale parameter
will be saved into a file named @file{@var{MODEL_FILENAME}_mh_scale.mat}. This
file can be loaded in subsequent runs via the
@code{posterior_sampler_options}-option @ref{scale_file}. Both
@code{mode_compute=6} and @code{scale_file} will overwrite any value specified
in @code{estimated_params} with the tuned value. Default: @code{0.2}
Note also that for the Random Walk Metropolis Hastings algorithm, it is
possible to use option @ref{mh_tune_jscale}, to automatically tune the value of
@code{mh_jscale}.
@item mh_init_scale = @var{DOUBLE}
The scale to be used for drawing the initial value of the
@ -5634,6 +5692,16 @@ If the @ref{nointeractive}-option has been invoked, the program will instead aut
@code{mh_init_scale} by 10 percent after 100 futile draws and try another 100 draws. This iterative
procedure will take place at most 10 times, at which point Dynare will abort with an error message.
@item mh_tune_jscale [= @var{DOUBLE}]
@anchor{mh_tune_jscale} Automatically tunes the scale parameter of the jumping
distribution's covariance matrix (Metropolis-Hastings), so that the overall
acceptance ratio is close to the desired level. Default value is
@code{0.33}. It is not possible to match exactly the desired
acceptance ratio because of the stochastic nature of the algoirithm (the
proposals and the initialial conditions of the markov chains if
@code{mh_nblocks>1}). This option is only available for the Random Walk
Metropolis Hastings algorithm.
@item mh_recover
@anchor{mh_recover} Attempts to recover a Metropolis-Hastings
simulation that crashed prematurely, starting with the last available saved
@ -7466,7 +7534,7 @@ graphs. See @code{colormap} in Matlab/Octave manual for valid arguments.
@item nograph
@xref{nograph}. Suppresses the display and creation only within the
@code{shock_decomposition}-command, but does not affect other commands.
@xref{plot_shock_decomposition} for plotting graphs.
@xref{plot_shock_decomposition}, for plotting graphs.
@item init_state = @var{BOOLEAN}
@anchor{init_state} If equal to @math{0}, the shock decomposition is computed conditional on the smoothed state
@ -7599,7 +7667,7 @@ model).
@item nograph
@xref{nograph}. Only shock decompositions are computed and stored in @code{oo_.realtime_shock_decomposition},
@code{oo_.conditional_shock_decomposition} and @code{oo_.realtime_forecast_shock_decomposition} but no plot is made
(@xref{plot_shock_decomposition}).
(@pxref{plot_shock_decomposition}).
@item presample = @var{INTEGER}
@anchor{presample_shock_decomposition} Data point above which recursive

1
doc/parallel/parallel.tex
View File

@ -2,6 +2,7 @@
% AMS-LaTeX Paper ************************************************
% **** -----------------------------------------------------------
\documentclass[12pt,a4paper,pdftex,nofootinbib]{article}
\usepackage[cp1252]{inputenc}
\usepackage{amssymb,amsmath}
\usepackage[pdftex]{graphicx}
\usepackage{epstopdf}

56
julia/DynareModel.jl
View File

@ -1,6 +1,6 @@
module DynareModel
##
# Copyright (C) 2015 Dynare Team
# Copyright (C) 2015-2018 Dynare Team
#
# This file is part of Dynare.
#
@ -20,7 +20,7 @@ module DynareModel
export Model, Endo, Exo, ExoDet, Param, dynare_model
abstract Atom
abstract type Atom end
immutable Endo <: Atom
name::String
@ -88,6 +88,7 @@ type Model
aux_vars::Vector{AuxVars}
pred_vars::Vector{Int}
obs_vars::Vector{Int}
state_var::Vector{Int}
orig_endo_nbr::Int
orig_eq_nbr::Int
eq_nbr::Int
@ -106,6 +107,14 @@ type Model
maximum_endo_lead::Int
maximum_exo_lag::Int
maximum_exo_lead::Int
orig_maximum_lag::Int
orig_maximum_lead::Int
orig_maximum_endo_lag::Int
orig_maximum_endo_lead::Int
orig_maximum_exo_lag::Int
orig_maximum_exo_lead::Int
orig_maximum_exo_det_lag::Int
orig_maximum_exo_det_lead::Int
lead_lag_incidence::Matrix{Int}
nnzderivatives::Vector{Int}
analytical_steady_state::Bool
@ -130,18 +139,19 @@ function dynare_model()
return Model("", # fname
"", # dname
"", # dynare_version
Array(Endo,0), # endo
Array(Exo,0), # exo
Array(ExoDet,0), # exo_det
Array(Param,0), # param
Array(AuxVars,0), # aux_vars
Array(Int,0), # pred_vars
Array(Int,0), # obs_vars
Vector{Endo}(), # endo
Vector{Exo}(), # exo
Vector{ExoDet}(), # exo_det
Vector{Param}(), # param
Vector{AuxVars}(), # aux_vars
Vector{Int}(), # pred_vars
Vector{Int}(), # obs_vars
Vector{Int}(), # state_var
0, # orig_endo_nbr
0, # orig_eq_nbr
0, # eq_nbr
0, # ramsey_eq_nbr
Array(DetShocks,0), # det_shocks
Vector{DetShocks}(), # det_shocks
0, # nstatic
0, # nfwrd
0, # npred
@ -155,19 +165,27 @@ function dynare_model()
0, # maximum_endo_lead
0, # maximum_exo_lag
0, # maximum_exo_lead
Array(Int, 3, 0), # lead_lag_incidence
zeros(Int, 3), # nnzderivatives
0, # orig_maximum_lag
0, # orig_maximum_lead
0, # orig_maximum_endo_lag
0, # orig_maximum_endo_lead
0, # orig_maximum_exo_lag
0, # orig_maximum_exo_lead
0, # orig_maximum_exo_det_lag
0, # orig_maximum_exo_det_lead
Matrix{Int}(0,0), # lead_lag_incidence
zeros(Int64,3), # nnzderivatives
false, # analytical_steady_state
false, # user_written_analytical_steady_state
false, # static_and_dynamic_models_differ
Array(String,0), # equation_tags
Array(Int64,1), # exo_names_orig_ord
Array(Float64, 0, 0), # sigma_e (Cov matrix of the structural innovations)
Array(Float64, 0, 0), # correlation_matrix (Corr matrix of the structural innovations)
Array(Float64, 0, 0), # h (Cov matrix of the measurement errors)
Array(Float64, 0, 0), # correlation_matrix_me (Cov matrix of the measurement errors)
Vector{String}(), # equation_tags
Vector{Int}(), # exo_names_orig_ord
Matrix{Float64}(0,0), # sigma_e (Cov matrix of the structural innovations)
Matrix{Float64}(0,0), # correlation_matrix (Corr matrix of the structural innovations)
Matrix{Float64}(0,0), # h (Cov matrix of the measurement errors)
Matrix{Float64}(0,0), # correlation_matrix_me (Cov matrix of the measurement errors)
true, # sigma_e_is_diagonal
Array(Float64, 0), # params
Vector{Float64}(), # params
function()end, # static
function()end, # static_params_derivs
function()end, # dynamic

6
julia/DynareOutput.jl
View File

@ -1,6 +1,6 @@
module DynareOutput
##
# Copyright (C) 2015 Dynare Team
# Copyright (C) 2015-2018 Dynare Team
#
# This file is part of Dynare.
#
@ -28,8 +28,8 @@ end
function dynare_output()
return Output("", # dynare_version
Array(Float64, 0), # steady_state
Array(Float64, 0) # exo_steady_state
Vector{Float64}(), # steady_state
Vector{Float64}() # exo_steady_state
)
end

4
julia/PerfectForesightModelSolver.jl
View File

@ -112,7 +112,7 @@ function simulate_perfect_foresight_model!(endogenousvariables::Matrix{Float64},
i_cols_A += ny
end
end
err = maximum(abs(rd))
err = maximum(abs.(rd))
println("Iter. ", iteration, "\t err. ", round(err, 12))
if err<options.pfmsolver.tolf
iteration -= 1
@ -121,7 +121,7 @@ function simulate_perfect_foresight_model!(endogenousvariables::Matrix{Float64},
A = sparse(iA[1:m], jA[1:m], vA[1:m])
dy = -A\rd
Y[i_upd] += dy
if maximum(abs(dy))<options.pfmsolver.tolx
if maximum(abs.(dy))<options.pfmsolver.tolx
convergence = true
end
end

16
julia/SteadyState.jl
View File

@ -47,22 +47,26 @@ function steady!(model::Model, oo::Output)
end
function steady(model::Model, oo::Output, yinit::Vector{Float64})
ojectivefunction!(y::Vector{Float64}, fval::Vector{Float64}, fjac::Array{Float64}) = model.static(y, oo.exo_steady_state, model.params, fval, fjac)
r = nlsolve(only_fg!(ojectivefunction!), yinit, show_trace=false)
f!(fval::Vector{Float64}, y::Vector{Float64}) = model.static(y, oo.exo_steady_state, model.params, fval)
j!(fjac::Matrix{Float64}, y::Vector{Float64}) = model.static(y, oo.exo_steady_state, model.params, fjac)
fj!(fval::Vector{Float64}, fjac::Matrix{Float64}, y::Vector{Float64}) = model.static(y, oo.exo_steady_state, model.params, fval, fjac)
r = nlsolve(f!, j!, fj!, yinit, show_trace=false)
if converged(r)
return r.zero
else
return fill(nan(Float64), length(yinit))
return fill(NaN, length(yinit))
end
end
function steady!(model::Model, oo::Output, yinit::Vector{Float64})
ojectivefunction!(y::Vector{Float64}, fval::Vector{Float64}, fjac::Array{Float64}) = model.static(y, oo.exo_steady_state, model.params, fval, fjac)
r = nlsolve(only_fg!(ojectivefunction!), yinit, show_trace=false)
f!(fval::Vector{Float64}, y::Vector{Float64}) = model.static(y, oo.exo_steady_state, model.params, fval)
j!(fjac::Matrix{Float64}, y::Vector{Float64}) = model.static(y, oo.exo_steady_state, model.params, fjac)
fj!(fval::Vector{Float64}, fjac::Matrix{Float64}, y::Vector{Float64}) = model.static(y, oo.exo_steady_state, model.params, fval, fjac)
r = nlsolve(f!, j!, fj!, yinit, show_trace=false)
if converged(r)
oo.steady_state = r.zero
else
oo.steady_state = fill(nan(Float64), length(yinit))
oo.steady_state = fill(NaN, length(yinit))
end
end

6
license.txt
View File

@ -9,7 +9,7 @@ Upstream-Contact: Dynare Team, whose members in 2017 are:
Ferhat Mihoubi <fmihoubi@univ-evry.fr>
Johannes Pfeifer <jpfeifer@gmx.de>
Marco Ratto <marco.ratto@ec.europa.eu>
Sébastien Villemot <sebastien.villemot@sciencespo.fr>
Sébastien Villemot <sebastien@dynare.org>
Source: http://www.dynare.org
Files: *
@ -41,8 +41,8 @@ License: public-domain-aim
Files: matlab/optimization/bfgsi1.m matlab/csolve.m matlab/optimization/csminit1.m matlab/optimization/numgrad2.m
matlab/optimization/numgrad3.m matlab/optimization/numgrad3_.m matlab/optimization/numgrad5.m
matlab/optimization/numgrad5_.m matlab/optimization/csminwel1.m matlab/bvar_density.m
matlab/bvar_toolbox.m matlab/partial_information/PI_gensys.m matlab/qzswitch.m
matlab/qzdiv.m
matlab/bvar_toolbox.m matlab/partial_information/PI_gensys.m matlab/partial_information/qzswitch.m
matlab/partial_information/qzdiv.m
Copyright: 1993-2009 Christopher Sims
2006-2016 Dynare Team
License: GPL-3+

3
m4/ax_matlab_version.m4
View File

@ -22,6 +22,9 @@ AC_REQUIRE([AX_MATLAB])
AC_MSG_CHECKING([for MATLAB version])
if test "x$MATLAB_VERSION" != "x"; then
case $MATLAB_VERSION in
*2018a | *2018A)
MATLAB_VERSION="9.4"
;;
*2017b | *2017B)
MATLAB_VERSION="9.3"
;;

10
matlab/add_path_to_mex_files.m
View File

@ -51,7 +51,7 @@ else
addpath(mexpath);
end
end
else
elseif matlab_ver_less_than('9.4')
tmp = [dynareroot '../mex/matlab/win64-7.8-9.3/'];
if exist(tmp, 'dir')
mexpath = tmp;
@ -59,6 +59,14 @@ else
addpath(mexpath);
end
end
else
tmp = [dynareroot '../mex/matlab/win64-9.4/'];
if exist(tmp, 'dir')
mexpath = tmp;
if modifypath
addpath(mexpath);
end
end
end
end
% Add OS X 64bits specific paths for Dynare Mac package

47
matlab/backward/simul_backward_model.m
View File

@ -38,25 +38,44 @@ function simulation = simul_backward_model(initialconditions, samplesize, innova
global options_ M_ oo_
if nargin<3
Innovations =[];
Innovations = [];
else
if isdseries(innovations)
if isequal(innovations.dates(1)-1, initialconditions.dates(end))
if innovations.nobs<samplesize
error('Time span in third argument is too short (should not be less than %s, the value of the second argument)', num2str(samplesize))
end
% Set array holding innovations values.
Innovations = zeros(samplesize, M_.exo_nbr);
exonames = M_.exo_names;
for i=1:M_.exo_nbr
if ismember(exonames{i}, innovations.name)
Innovations(:,i) = innovations{exonames{i}}.data(1:samplesize);
else
disp(sprintf('Exogenous variable %s is not available in third argument, default value is zero.', exonames{i}));
if isdseries(initialconditions)
if isequal(innovations.dates(1)-1, initialconditions.dates(end))
if innovations.nobs<samplesize
error('Time span in third argument is too short (should not be less than %s, the value of the second argument)', num2str(samplesize))
end
% Set array holding innovations values.
Innovations = zeros(samplesize, M_.exo_nbr);
exonames = M_.exo_names;
for i=1:M_.exo_nbr
if ismember(exonames{i}, innovations.name)
Innovations(:,i) = innovations{exonames{i}}.data(1:samplesize);
else
disp(sprintf('Exogenous variable %s is not available in third argument, default value is zero.', exonames{i}));
end
end
else
error('Time spans in first and third arguments should be contiguous!')
end
else
error('Time spans in first and third arguments should be contiguous!')
if isempty(initialconditions)
if innovations.nobs<samplesize
error('Time span in third argument is too short (should not be less than %s, the value of the second argument)', num2str(samplesize))
end
Innovations = zeros(samplesize, M_.exo_nbr);
exonames = M_.exo_names;
for i=1:M_.exo_nbr
if ismember(exonames{i}, innovations.name)
Innovations(:,i) = innovations{exonames{i}}.data(1:samplesize);
else
disp(sprintf('Exogenous variable %s is not available in third argument, default value is zero.', exonames{i}));
end
end
else
error('First input must be an empty array!')
end
end
else
error('Third argument must be a dseries object!')

111
matlab/calibrate_mh_scale_parameter.m Normal file
View File

@ -0,0 +1,111 @@
function Scale = calibrate_mh_scale_parameter(ObjectiveFunction, CovarianceMatrix, Parameters, MhBounds, options, varargin)
% Tune the MH scale parameter so that the overall acceptance ratio is close to AcceptanceTarget.
%
% INPUTS
% - ObjectiveFunction [fhandle] Function (posterior kernel).
% - CovarianceMatrix [double] n*n matrix, covariance matrix of the jumping distribution.
% - Parameters [double] n*1 vector, parameter values.
% - MhBounds [double] n*2 matrix, bounds on the possible values for the parameters.
% - options [structure] content of options_.tune_mh_jscale.
% - varargin [cell] Additional arguments to be passed to ObjectiveFunction.
%
% OUTPUTS
% - Scale [double] scalar, optimal scale parameter for teh jumping distribution.
% Fire up the wait bar
hh = dyn_waitbar(0,'Tuning of the scale parameter...');
set(hh,'Name','Tuning of the scale parameter.');
% Intilialize various counters.
j = 1; jj = 1; isux = 0; jsux = 0; i = 0;
% Evaluate the objective function.
logpo0 = - feval(ObjectiveFunction, Parameters, varargin{:});
logpo1 = logpo0;
% Get the dimension of the problem.
n = length(Parameters);
% Initialize the correction on the scale factor.
correction = 1.0;
% Set the initial value of the scale parameter
Scale = options.guess;
% Transposition of some arrays.
MhBounds = MhBounds';
Parameters = Parameters';
% Compute the Cholesky of the covariance matrix, return an error if the
% matrix is not positive definite.
try
dd = chol(CovarianceMatrix);
catch
error('The covariance matrix has to be a symetric positive definite matrix!')
end
% Set parameters related to the proposal distribution
if options.rwmh.proposal_distribution=='rand_multivariate_normal'
nu = n;
elseif options.rwmh.proposal_distribution=='rand_multivariate_student'
nu = options.rwmh.student_degrees_of_freedom;
end
% Random Walk Metropolis Hastings iterations...
while j<=options.maxiter
% Obtain a proposal (jump)
proposal = feval(options.rwmh.proposal_distribution, Parameters, Scale*dd, nu);
% If out of boundaries set the posterior kernel equal to minus infinity
% so that the proposal will be rejected with probability one.
if all(proposal > MhBounds(1,:)) && all(proposal < MhBounds(2,:))
logpo0 = -feval(ObjectiveFunction, proposal(:), varargin{:});
else
logpo0 = -inf;
end
% Move if the proposal is enough likely...
if logpo0>-inf && log(rand)<logpo0-logpo1
Parameters = proposal;
logpo1 = logpo0;
isux = isux + 1;
jsux = jsux + 1;
end% ... otherwise I don't move.
prtfrc = j/options.maxiter;
% Update the waitbar
if ~mod(j, 10)
dyn_waitbar(prtfrc, hh, sprintf('Acceptance ratio [during last 500]: %f [%f]', isux/j, jsux/jj));
end
% Adjust the value of the scale parameter.
if ~mod(j, options.stepsize)
r1 = jsux/jj; % Local acceptance ratio
r2 = isux/j; % Overall acceptance ratio
% Set correction for the scale factor
c1 = r1/options.target;
if abs(c1-1)>.05
correction = correction^options.rho*c1^(1-options.rho);
else
correction = c1;
end
% Apply correction
if c1>0
Scale = Scale*correction;
else
Scale = Scale/10;
end
% Update some counters.
jsux = 0; jj = 0;
if abs(r2-options.target)<options.c2 && abs(r1-options.target)<options.c1
i = i+1;
else
i = 0;
end
% Test convergence.
if i>options.c3
break
end
end
j = j+1;
jj = jj + 1;
end
dyn_waitbar_close(hh);

57
matlab/check.m
View File

@ -1,46 +1,18 @@
function [eigenvalues_,result,info] = check(M, options, oo)
% Checks determinacy conditions by computing the generalized eigenvalues.
%
% INPUTS
% - M [structure] Matlab's structure describing the model (M_).
% - options [structure] Matlab's structure describing the current options (options_).
% - oo [structure] Matlab's structure containing the results (oo_).
%
% OUTPUTS
% - eigenvalues_ [double] vector, eigenvalues.
% - 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.
%@info:
%! @deftypefn {Function File} {[result,info] =} check (@var{M},@var{options},@var{oo})
%! @anchor{check}
%! @sp 1
%! Checks determinacy conditions by computing the generalized eigenvalues.
%! @sp 2
%! @strong{Inputs}
%! @sp 1
%! @table @ @var
%! @item M
%! Matlab's structure describing the model (initialized by dynare).
%! @item options
%! Matlab's structure describing the options (initialized by dynare).
%! @item oo
%! Matlab's structure gathering the results (initialized by dynare).
%! @end table
%! @sp 2
%! @strong{Outputs}
%! @sp 1
%! @table @ @var
%! @item eigenvalues_
%! Eigenvalues of the model.
%! @item result
%! Integer scalar equal to one (BK conditions are satisfied) or zero (otherwise).
%! @item info
%! Integer scalar, error code as returned by @ref{resol}.
%! @end table
%! @sp 2
%! @strong{This function is called by:}
%! @sp 1
%! @ref{smm_objective}
%! @sp 2
%! @strong{This function calls:}
%! @sp 1
%! @ref{resol}
%! None.
%! @end deftypefn
%@eod:
% Copyright (C) 2001-2014 Dynare Team
% Copyright (C) 2001-2018 Dynare Team
%
% This file is part of Dynare.
%
@ -78,7 +50,6 @@ end
eigenvalues_ = dr.eigval;
[m_lambda,i]=sort(abs(eigenvalues_));
n_explod = nnz(abs(eigenvalues_) > options.qz_criterium);
% Count number of forward looking variables
if ~options.block
@ -91,7 +62,7 @@ else
end
result = 0;
if (nyf == n_explod) && (dr.full_rank)
if (nyf == dr.edim) && (dr.full_rank)
result = 1;
end
@ -101,7 +72,7 @@ if options.noprint == 0
disp(sprintf('%16s %16s %16s\n','Modulus','Real','Imaginary'))
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 ', n_explod));
disp(sprintf('\nThere are %d eigenvalue(s) larger than 1 in modulus ', dr.edim));
disp(sprintf('for %d forward-looking variable(s)',nyf));
skipline()
if result

4
matlab/dr_block.m
View File

@ -593,7 +593,7 @@ for i = 1:Size
else
[err, ghx_other] = gensylv(1, A_, B_, C_, -D_);
end
if options_.aim_solver ~= 1 && options_.use_qzdiv
if options_.aim_solver ~= 1
% Necessary when using Sims' routines for QZ
ghx_other = real(ghx_other);
end
@ -650,7 +650,7 @@ for i = 1:Size
end
if options_.aim_solver ~= 1 && options_.use_qzdiv
if options_.aim_solver ~= 1
% Necessary when using Sims' routines for QZ
ghx = real(ghx);
ghu = real(ghu);

83
matlab/dyn_first_order_solver.m
View File

@ -1,53 +1,27 @@
function [dr,info] = dyn_first_order_solver(jacobia,DynareModel,dr,DynareOptions,task)
function [dr, info] = dyn_first_order_solver(jacobia, DynareModel, dr, DynareOptions, task)
%@info:
%! @deftypefn {Function File} {[@var{dr},@var{info}] =} dyn_first_order_solver (@var{jacobia},@var{DynareModel},@var{dr},@var{DynareOptions},@var{task})
%! @anchor{dyn_first_order_solver}
%! @sp 1
%! Computes the first order reduced form of the DSGE model
%! @sp 2
%! @strong{Inputs}
%! @sp 1
%! @table @ @var
%! @item jacobia
%! Matrix containing the Jacobian of the model
%! @item DynareModel
%! Matlab's structure describing the model (initialized by @code{dynare}).
%! @item dr
%! Matlab's structure describing the reduced form solution of the model.
%! @item qz_criterium
%! Double containing the criterium to separate explosive from stable eigenvalues
%! @end table
%! @sp 2
%! @strong{Outputs}
%! @sp 1
%! @table @ @var
%! @item dr
%! Matlab's structure describing the reduced form solution of the model.
%! @item info
%! Integer scalar, error code.
%! @sp 1
%! @table @ @code
%! @item info==0
%! No error.
%! @item info==1
%! The model doesn't determine the current variables uniquely.
%! @item info==2
%! MJDGGES returned an error code.
%! @item info==3
%! Blanchard & Kahn conditions are not satisfied: no stable equilibrium.
%! @item info==4
%! Blanchard & Kahn conditions are not satisfied: indeterminacy.
%! @item info==5
%! Blanchard & Kahn conditions are not satisfied: indeterminacy due to rank failure.
%! @item info==7
%! One of the generalized eigenvalues is close to 0/0
%! @end table
%! @end table
%! @end deftypefn
%@eod:
% Computes the first order reduced form of a DSGE model.
%
% INPUTS
% - jacobia [double] matrix, the jacobian of the dynamic model.
% - DynareModel [struct] Matlab's structre describing the model, M_ global.
% - dr [struct] Matlab's structure describing the reduced form model.
% - DynareOptions [struct] Matlab's structure containing the current state of the options, oo_ global.
% - task [integer] scalar, if task = 0 then decision rules are computed and if task = 1 then only eigenvales are computed.
%
% OUTPUTS
% - dr [struct] Matlab's structure describing the reduced form model.
% - info [integer] scalar, error code. Possible values are:
%
% info=0 -> no error,
% info=1 -> the model doesn't determine the current variables uniquely,
% info=2 -> mjdgges dll returned an error,
% info=3 -> Blanchard and Kahn conditions are not satisfied: no stable equilibrium,
% info=4 -> Blanchard and Kahn conditions are not satisfied: indeterminacy,
% info=5 -> Blanchard and Kahn conditions are not satisfied: indeterminacy due to rank failure,
% info=7 -> One of the eigenvalues is close to 0/0 (infinity of complex solutions)
% Copyright (C) 2001-2017 Dynare Team
% Copyright (C) 2001-2018 Dynare Team
%
% This file is part of Dynare.
%
@ -228,15 +202,17 @@ else
return
end
dr.sdim = sdim; % Number of stable eigenvalues.
dr.edim = length(dr.eigval)-sdim; % Number of exposive eigenvalues.
nba = nd-sdim;
if task == 1
if task==1
if rcond(w(npred+nboth+1:end,npred+nboth+1:end)) < 1e-9
dr.full_rank = 0;
else
dr.full_rank = 1;
end
return
end
if nba ~= nsfwrd
@ -251,8 +227,11 @@ else
info(2) = temp'*temp;
return
end
if task==1, return, end
%First order approximation
indx_stable_root = 1: (nd - nsfwrd); %=> index of stable roots
indx_stable_root = 1: (nd - nsfwrd); %=> index of stable roots
indx_explosive_root = npred + nboth + 1:nd; %=> index of explosive roots
% derivatives with respect to dynamic state variables
% forward variables
@ -318,7 +297,7 @@ end
dr.ghx = ghx;
dr.ghu = ghu;
if DynareOptions.aim_solver ~= 1 && DynareOptions.use_qzdiv
if DynareOptions.aim_solver ~= 1
% Necessary when using Sims' routines for QZ
dr.ghx = real(ghx);
dr.ghu = real(ghu);

19
matlab/dynare.m
View File

@ -70,9 +70,20 @@ dynareroot = dynare_config('', preprocessoroutput);
warning_config()
if ~isoctave
if isoctave
if octave_ver_less_than(supported_octave_version)
skipline()
warning(['This version of Octave is not supported. Consider installing ' ...
'version %s+ of Octave,\notherwise m files will be used instead ' ...
'of precompiled mex files and some features, like solution\n' ...
'of models approximated at third order, will not be available.'], supported_octave_version())
skipline()
end
else
if matlab_ver_less_than('7.5')
skipline()
warning('This version of Dynare has only been tested on MATLAB 7.5 (R2007b) and above. Since your MATLAB version is older than that, Dynare may fail to run, or give unexpected results. Consider upgrading your MATLAB installation, or switch to Octave.');
skipline()
end
end
@ -81,11 +92,7 @@ more off
% sets default format for save() command
if isoctave
if octave_ver_less_than('3.8')
default_save_options('-mat')
else
save_default_options('-mat')
end
save_default_options('-mat')
end
if nargin < 1

21
matlab/dynare_config.m
View File

@ -96,20 +96,6 @@ if isoctave
p{end+1} = '/missing/ordeig';
end
if isoctave && ~isequal(supported_octave_version(), version())
skipline()
warning(['This version of Octave is not supported. Consider installing ' ...
'version %s of Octave,\notherwise m files will be used instead ' ...
'of precompiled mex files and some features, like solution\n' ...
'of models approximated at third order, will not be available.'], supported_octave_version())
skipline()
end
% ilu is missing in Octave < 4.0
if isoctave && octave_ver_less_than('4.0')
p{end+1} = '/missing/ilu';
end
% corrcoef with two outputs is missing in Octave (ticket #796)
if isoctave && ~user_has_octave_forge_package('nan')
p{end+1} = '/missing/corrcoef';
@ -132,6 +118,11 @@ if isoctave || matlab_ver_less_than('9.3')
p{end+1} = '/missing/isfile';
end
% strsplit is missing in Matlab<R2013a
if ~isoctave && matlab_ver_less_than('8.1')
p{end+1} = '/missing/strsplit';
end
P = cellfun(@(c)[dynareroot(1:end-1) c], p, 'uni',false);
% Get mex files folder(s)
@ -249,4 +240,4 @@ end
% Initialization of the dates and dseries classes (recursive).
initialize_dseries_toolbox;
cd(origin);
cd(origin);

17
matlab/dynare_estimation_1.m
View File

@ -228,7 +228,7 @@ if ~isequal(options_.mode_compute,0) && ~options_.mh_posterior_mode_estimation
elseif isnumeric(options_.mode_compute) && options_.mode_compute==6 %save scaling factor
save([M_.fname '_optimal_mh_scale_parameter.mat'],'Scale');
options_.mh_jscale = Scale;
bayestopt_.jscale = ones(length(xparam1),1)*Scale;
bayestopt_.jscale(:) = options_.mh_jscale;
end
if ~isnumeric(options_.mode_compute) || ~isequal(options_.mode_compute,6) %always already computes covariance matrix
if options_.cova_compute == 1 %user did not request covariance not to be computed
@ -438,6 +438,21 @@ if (any(bayestopt_.pshape >0 ) && options_.mh_replic) || ...
error(['Estimation:: Mode value(s) of ', disp_string ,' are outside parameter bounds.'])
end
end
% Tunes the jumping distribution's scale parameter
if options_.mh_tune_jscale.status
if options_.posterior_sampler_options.posterior_sampling_method=='random_walk_metropolis_hastings'
options = options_.mh_tune_jscale;
options.rwmh = options_.posterior_sampler_options.rwmh;
options_.mh_jscale = calibrate_mh_scale_parameter(objective_function, ...
invhess, xparam1, [bounds.lb,bounds.ub], ...
options, dataset_, dataset_info, options_, M_, estim_params_, bayestopt_, bounds, oo_);
bayestopt_.jscale(:) = options_.mh_jscale;
disp(sprintf('mh_jscale has been set equal to %s', num2str(options_.mh_jscale)))
skipline()
else
warning('mh_tune_jscale is only available with Random Walk Metropolis Hastings!')
end
end
% runs MCMC
if options_.mh_replic || options_.load_mh_file
posterior_sampler_options = options_.posterior_sampler_options.current_options;

6
matlab/gensylv/sylvester3.m
View File

@ -38,11 +38,7 @@ end
if isoctave
[aa,bb,qq,zz]=qz(full(a),full(b));
for j=1:p
if octave_ver_less_than('3.4.0')
d(:,:,j)=qq'*d(:,:,j)*u;
else
d(:,:,j)=qq*d(:,:,j)*u;
end
d(:,:,j)=qq*d(:,:,j)*u;
end
else
[aa,bb,qq,zz]=qz(full(a),full(b),'real'); % available in Matlab version 6.0

18
matlab/global_initialization.m
View File

@ -334,15 +334,6 @@ options_.linear = 0;
options_.replic = 50;
options_.simul_replic = 1;
options_.drop = 100;
% if mjdgges.dll (or .mexw32 or ....) doesn't exist, matlab/qz is added to the path.
% There exists now qz/mjdgges.m that contains the calls to the old Sims code
% Hence, if mjdgges.m is visible exist(...)==2,
% this means that the DLL isn't avaiable and use_qzdiv is set to 1
if exist('mjdgges','file')==2
options_.use_qzdiv = 1;
else
options_.use_qzdiv = 0;
end
options_.aim_solver = 0; % i.e. by default do not use G.Anderson's AIM solver, use mjdgges instead
options_.k_order_solver=0; % by default do not use k_order_perturbation but mjdgges
options_.partial_information = 0;
@ -442,6 +433,15 @@ options_.mh_conf_sig = 0.90;
options_.prior_interval = 0.90;
options_.mh_drop = 0.5;
options_.mh_jscale = 0.2;
options_.mh_tune_jscale.status = false;
options_.mh_tune_jscale.guess = .2;
options_.mh_tune_jscale.target = .33;
options_.mh_tune_jscale.maxiter = 200000;
options_.mh_tune_jscale.rho = .7;
options_.mh_tune_jscale.stepsize = 1000;
options_.mh_tune_jscale.c1 = .02;
options_.mh_tune_jscale.c2 = .05;
options_.mh_tune_jscale.c3 = 4;
options_.mh_init_scale = 2*options_.mh_jscale;
options_.mh_mode = 1;
options_.mh_nblck = 2;

27
matlab/missing/strsplit/private/ischarint.m Normal file
View File

@ -0,0 +1,27 @@
function l = ischarint(x)
% Returns true if and only if char x represents an integer.
% Copyright © 2018 DynareTeam
%
% 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/>.
s = warning;
warning off;
l = isint(str2double(x));
warning(s);

34
matlab/missing/ilu/ilu.m → matlab/missing/strsplit/private/ischarnum.m
View File

@ -1,7 +1,8 @@
function [L, U, P] = ilu(A, setup)
% Partially implement function ilu using the (now deprecated) luinc
function l = ischarnum(x)
% Copyright (C) 2013-2017 Dynare Team
% Returns true if and only if char x represents a real number.
% Copyright © 2018 DynareTeam
%
% This file is part of Dynare.
%
@ -18,20 +19,17 @@ function [L, U, P] = ilu(A, setup)
% You should have received a copy of the GNU General Public License
% along with Dynare. If not, see <http://www.gnu.org/licenses/>.
if nargout ~= 2
error('Only two output arguments supported')
end
if nargin ~= 2
error('Only two input arguments supported')
end
l = false;
if isfield(setup, 'milu')
if setup.milu == 'off'
setup.milu = 0;
else
error(['Unsupported milu: ' setup.milu ])
s = warning;
warning off;
number = str2double(x);
warning(s);
if ~isempty(number)
if isreal(number)
l = true;
end
end
[L, U] = luinc(A, setup);
end
end

90
matlab/missing/strsplit/strsplit.m Normal file
View File

@ -0,0 +1,90 @@
function tok = strsplit(str, delimiters)
% Splits a string into multiple terms.
%
% INPUTS
% - str [char] String to be splitted.
% - delimiters [char, cell(char)] Delimiters.
%
% OUTPUTS
% - tok [cell(char)] Terms.
% Copyright © 2018 DynareTeam
%
% 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/>.
remove_empty = true;
remove_numbers = false;
% Check first input arguments
assert(ischar(str) && ndims(str)==2 && size(str,1)<=1, 'The first arugment has to be a row char array!');
% Set default value for second input arguments
if nargin<2
delimiters = {' '};
end
% If second input argument is a char transform it into a sigleton cell of char
if nargin>1
if ischar(delimiters)
assert(ndims(delimiters)==2 && size(delimiters,1)==1, 'The second input argument has to be be a char string!');
delimiters = {delimiters};
end
end
% Check that `delimiters` is a one dimensional cell
assert(ndim(delimiters)<=1, 'The second input argument has to be a one dimensional cell array!')
% Check that `delimiters` is a cell of row char arrays
assert(all(cellfun(@ischar, delimiters)) && all(cellfun(@rows, delimiters)==1), 'The second input argument has to be a cell of row char arrays!')
% If space is one of the delimiters obtain the index in `delimiters`
idspace = strmatch(' ', delimiters);
% Get the number of delimiters
n = length(delimiters);
% Remove unnecessary spaces
delimiters(setdiff(1:n, idspace)) = strtrim(delimiters(setdiff(1:n, idspace)));
% Join all the delimiters (strjoin is not available with matlab version less than R2013a)
if n>1
delimiter = '';
for i=1:n
if isspace(delimiters{i})
delimiter = horzcat(delimiter, '\s');
else
delimiter = horzcat(delimiter, delimiters{i});
end
delimiter = horzcat(delimiter,'|');
end
delimiter = horzcat(delimiter, '\W');
else
delimiter = delimiters{1};
end
% Get tokens
tok = regexp(str, delimiter, 'split');
if remove_empty
% Remove empty tokens
tok = tok(find(~cellfun(@isempty, tok)));
end
if remove_numbers
% Remove numbers
tok = tok(find(~cellfun(@ischarnum, tok)));
end

2
matlab/modules/dates

@ -1 +1 @@
Subproject commit 03a409517a9dc373eec8a0b6c98480f0fc8c145c
Subproject commit 5db5eb444d923ac74a553899b8dbfcc1cb6faa7a

0
matlab/qzdiv.m → matlab/partial_information/qzdiv.m
View File

0
matlab/qzswitch.m → matlab/partial_information/qzswitch.m
View File

5
matlab/plot_icforecast.m
View File

@ -27,11 +27,6 @@ function plot_icforecast(Variables,periods,options_)
% You should have received a copy of the GNU General Public License
% along with Dynare. If not, see <http://www.gnu.org/licenses/>.
if isoctave && octave_ver_less_than('3.4.0')
% The set() command on the handle returned by area() crashes in Octave 3.2
error('plot_conditional_forecast: you need Octave >= 3.4 (because of a bug in older versions)')
end
load conditional_forecasts;
forecast_periods = length(forecasts.cond.Mean.(Variables{1}));

73
matlab/qz/mjdgges.m
View File

@ -1,28 +1,25 @@
function [err,ss,tt,w,sdim,eigval,info] = mjdgges(e,d,qz_criterium, fake)
%function [err,ss,tt,w,sdim,eigval,info] = mjdgges(e,d,qz_criterium)
% QZ decomposition, Sims' codes are used.
function [err, ss, tt, zz, sdim, eigval, info] = mjdgges(e, d, qz_criterium, zhreshold)
%
% INPUTS
% e [double] real square (n*n) matrix.
% d [double] real square (n*n) matrix.
% qz_criterium [double] scalar (1+epsilon).
% zhreshold [double] ignored (in the DLL, this parameter is used for
% detecting eigenvalues too close to 0/0)
%
% OUTPUTS
% err [double] scalar: 1 indicates failure, 0 indicates success
% ss [complex] (n*n) matrix.
% tt [complex] (n*n) matrix.
% w [complex] (n*n) matrix.
% sdim [integer] scalar.
% eigval [complex] (n*1) vector.
% err [integer] scalar: 1 indicates failure, 0 indicates success
% ss [double] (n*n) quasi-triangular matrix.
% tt [double] (n*n) quasi-triangular matrix.
% zz [double] (n*n) orthogonal matrix.
% sdim [integer] scalar: number of stable eigenvalues.
% eigval [complex] (n*1) vector of generalized eigenvalues.
% info [integer] scalar.
%
% ALGORITHM
% Sims's qzdiv routine is used.
%
% SPECIAL REQUIREMENTS
% none.
% Copyright (C) 1996-2017 Dynare Team
% Copyright (C) 1996-2018 Dynare Team
%
% This file is part of Dynare.
%
@ -39,49 +36,35 @@ function [err,ss,tt,w,sdim,eigval,info] = mjdgges(e,d,qz_criterium, fake)
% You should have received a copy of the GNU General Public License
% along with Dynare. If not, see <http://www.gnu.org/licenses/>.
% Check number of inputs and outputs.
if nargin>4 || nargin<2 || nargout>7 || nargout==0
error('MJDGGES: takes 2 or 3 input arguments and between 1 and 7 output arguments.')
if nargin > 5 || nargin < 2 || nargout > 7 || nargout == 0
error('MJDGGES: takes 2, 3 or 4 input arguments and between 1 and 7 output arguments.')
end
% Check the first two inputs.
[me,ne] = size(e);
[md,nd] = size(d);
if ( ~isreal(e) || ~isreal(d) || me~=ne || md~=nd || me~=nd)
% info should be negative in this case, see dgges.f.
if isoctave
error('Octave unsupported, since it does not have real qz, ordqz and ordeig')
end
[me, ne] = size(e);
[md, nd] = size(d);
if ~isreal(e) || ~isreal(d) || me ~= ne || md ~= nd || me ~= nd
error('MJDGGES requires two square real matrices of the same dimension.')
end
% Set default value of qz_criterium.
if nargin <3
if nargin < 3
qz_criterium = 1 + 1e-6;
end
info = 0;
% Initialization of the output arguments.
ss = zeros(ne,ne);
tt = zeros(ne,ne);
w = zeros(ne,ne);
sdim = 0;
eigval = zeros(ne,1);
% Computational part.
try
if isoctave()
% Force complex QZ factorization.
e = complex(e);
d = complex(d);
end
[ss,tt,qq,w,a,b,c] = qz(e, d);
warning_old_state = warning;
warning off;
[tt,ss,qq,w] = qzdiv(qz_criterium, tt, ss, qq, w);
eigval = diag(ss)./diag(tt);
warning(warning_old_state);
sdim = sum(abs(eigval) < qz_criterium);
[ss, tt, qq, zz] = qz(e, d, 'real');
eigval = ordeig(ss, tt);
select = abs(eigval) < qz_criterium;
sdim = sum(select);
[ss, tt, qq, zz] = ordqz(ss, tt, qq, zz, select);
eigval = ordeig(ss, tt);
catch
info = 1;% Not as precise as lapack's info!
info = 1; % Not as precise as lapack's info!
end
err = 0;
err = 0;

111
matlab/resol.m
View File

@ -1,80 +1,37 @@
function [dr,info,M,options,oo] = resol(check_flag,M,options,oo)
function [dr, info, M, options, oo] = resol(check_flag, M, options, oo)
%@info:
%! @deftypefn {Function File} {[@var{dr},@var{info},@var{M},@var{options},@var{oo}] =} resol (@var{check_flag},@var{M},@var{options},@var{oo})
%! @anchor{resol}
%! @sp 1
%! Computes the perturbation-based decisions rules of the DSGE model (orders 1 to 3).
%! @sp 2
%! @strong{Inputs}
%! @sp 1
%! @table @ @var
%! @item check_flag
%! Integer scalar, equal to 0 if all the approximation is required, positive if only the eigenvalues are to be computed.
%! @item M
%! Matlab's structure describing the model (initialized by @code{dynare}).
%! @item options
%! Matlab's structure describing the options (initialized by @code{dynare}).
%! @item oo
%! Matlab's structure gathering the results (initialized by @code{dynare}).
%! @end table
%! @sp 2
%! @strong{Outputs}
%! @sp 1
%! @table @ @var
%! @item dr
%! Matlab's structure describing the reduced form solution of the model.
%! @item info
%! Integer scalar, error code.
%! @sp 1
%! @table @ @code
%! @item info==0
%! No error.
%! @item info==1
%! The model doesn't determine the current variables uniquely.
%! @item info==2
%! MJDGGES returned an error code.
%! @item info==3
%! Blanchard & Kahn conditions are not satisfied: no stable equilibrium.
%! @item info==4
%! Blanchard & Kahn conditions are not satisfied: indeterminacy.
%! @item info==5
%! Blanchard & Kahn conditions are not satisfied: indeterminacy due to rank failure.
%! @item info==6
%! The jacobian evaluated at the deterministic steady state is complex.
%! @item info==19
%! The steadystate routine has thrown an exception (inconsistent deep parameters).
%! @item info==20
%! Cannot find the steady state, info(2) contains the sum of square residuals (of the static equations).
%! @item info==21
%! The steady state is complex, info(2) contains the sum of square of imaginary parts of the steady state.
%! @item info==22
%! The steady has NaNs.
%! @item info==23
%! M_.params has been updated in the steadystate routine and has complex valued scalars.
%! @item info==24
%! M_.params has been updated in the steadystate routine and has some NaNs.
%! @item info==30
%! Ergodic variance can't be computed.
%! @end table
%! @sp 1
%! @item M
%! Matlab's structure describing the model (initialized by @code{dynare}).
%! @item options
%! Matlab's structure describing the options (initialized by @code{dynare}).
%! @item oo
%! Matlab's structure gathering the results (initialized by @code{dynare}).
%! @end table
%! @sp 2
%! @strong{This function is called by:}
%! @sp 1
%! @ref{dynare_estimation_init}
%! @sp 2
%! @strong{This function calls:}
%! @sp 1
%! None.
%! @end deftypefn
%@eod:
% Computes the perturbation based decision rules of the DSGE model (orders 1 to 3)
%
% INPUTS
% - check_flag [integer] scalar, equal to 0 if all the approximation is required, equal to 1 if only the eigenvalues are to be computed.
% - M [structure] Matlab's structure describing the model (M_).
% - options [structure] Matlab's structure describing the current options (options_).
% - oo [structure] Matlab's structure containing the results (oo_).
%
% OUTPUTS
% - dr [structure] Reduced form model.
% - info [integer] scalar or vector, error code.
% - M [structure] Matlab's structure describing the model (M_).
% - options [structure] Matlab's structure describing the current options (options_).
% - oo [structure] Matlab's structure containing the results (oo_).
%
% REMARKS
% Possible values for the error codes are:
%
% info(1)=0 -> No error.
% info(1)=1 -> The model doesn't determine the current variables uniquely.
% info(1)=2 -> MJDGGES returned an error code.
% info(1)=3 -> Blanchard & Kahn conditions are not satisfied: no stable equilibrium.
% info(1)=4 -> Blanchard & Kahn conditions are not satisfied: indeterminacy.
% info(1)=5 -> Blanchard & Kahn conditions are not satisfied: indeterminacy due to rank failure.
% info(1)=6 -> The jacobian evaluated at the deterministic steady state is complex.
% info(1)=19 -> The steadystate routine has thrown an exception (inconsistent deep parameters).
% info(1)=20 -> Cannot find the steady state, info(2) contains the sum of square residuals (of the static equations).
% info(1)=21 -> The steady state is complex, info(2) contains the sum of square of imaginary parts of the steady state.
% info(1)=22 -> The steady has NaNs.
% info(1)=23 -> M_.params has been updated in the steadystate routine and has complex valued scalars.
% 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 (C) 2001-2018 Dynare Team
%
@ -140,6 +97,8 @@ end
if options.block
[dr,info,M,options,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

50
matlab/stochastic_solvers.m
View File

@ -1,33 +1,28 @@
function [dr,info] = stochastic_solvers(dr,task,M_,options_,oo_)
% function [dr,info,M_,options_,oo_] = stochastic_solvers(dr,task,M_,options_,oo_)
% computes the reduced form solution of a rational expectations model (first, second or third
function [dr, info] = stochastic_solvers(dr, task, M_, options_, oo_)
% Computes the reduced form solution of a rational expectations model (first, second or third
% order approximation of the stochastic model around the deterministic steady state).
%
% INPUTS
% dr [matlab structure] Decision rules for stochastic simulations.
% task [integer] if task = 0 then dr1 computes decision rules.
% if task = 1 then dr1 computes eigenvalues.
% M_ [matlab structure] Definition of the model.
% options_ [matlab structure] Global options.
% oo_ [matlab structure] Results
% - dr [struct] Decision rules for stochastic simulations.
% - task [integer] scalar, if task = 0 then decision rules are computed and if task = 1 then only eigenvales are computed.
% - M_ [struct] Definition of the model.
% - options_ [struct] Options.
% - oo_ [struct] Results
%
% 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.
% info=9: k_order_pert was unable to compute the solution
% ALGORITHM
% ...
%
% SPECIAL REQUIREMENTS
% none.
% - dr [struct] Decision rules for stochastic simulations.
% - info [integer] scalar, error code:
%
% 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.
% info=9 -> k_order_pert was unable to compute the solution
% Copyright (C) 1996-2018 Dynare Team
%
@ -237,10 +232,11 @@ if M_.maximum_endo_lead == 0
end
dr.eigval = eig(kalman_transition_matrix(dr,nstatic+(1:nspred),1:nspred,M_.exo_nbr));
dr.full_rank = 1;
if any(abs(dr.eigval) > options_.qz_criterium)
dr.edim = nnz(abs(dr.eigval) > options_.qz_criterium);
dr.sdim = nd-dr.edim;
if dr.edim
temp = sort(abs(dr.eigval));
nba = nnz(abs(dr.eigval) > options_.qz_criterium);
temp = temp(nd-nba+1:nd)-1-options_.qz_criterium;
temp = temp(dr.sdim+1:nd)-1-options_.qz_criterium;
info(1) = 3;
info(2) = temp'*temp;
end

6
matlab/warning_config.m
View File

@ -34,11 +34,7 @@ warning('on', 'backtrace');
if isoctave
warning('off', 'Octave:separator-insert');
if octave_ver_less_than('4.0')
warning('off', 'Octave:matlab-incompatible');
else
warning('off', 'Octave:language-extension');
end
warning('off', 'Octave:language-extension');
warning('off', 'Octave:single-quote-string');
warning('off', 'Octave:missing-semicolon');
warning('off', 'Octave:empty-list-elements');

5
mex/build/matlab/configure.ac
View File

@ -58,6 +58,7 @@ CXXFLAGS="$CXXFLAGS -Wall -Wno-parentheses"
AC_PROG_F77([gfortran g77 f77])
AC_PROG_CC
AC_PROG_CC_C99 # mjdgges DLL now uses C99 features (variable declared in for loop)
AC_PROG_CXX
AC_PROG_RANLIB
AX_PROG_LN_S
@ -68,6 +69,10 @@ case ${host_os} in
*mingw32*)
# Ensure that -lpthread is statically linked under MinGW
PTHREAD_LIBS="-Wl,-Bstatic -lpthread -Wl,-Bdynamic"
# Kludge for bug in MinGW, that defines __STDC_UTF_16__ but not char16_t
# This breaks the matrix.h of older MATLABs (e.g. R2009a)
# Also see <uchar.h>
CFLAGS="$CFLAGS -include stdint.h -Dchar16_t=uint_least16_t"
;;
esac
AX_PTHREAD

8
mex/build/octave/Makefile.am
View File

@ -4,14 +4,6 @@ ACLOCAL_AMFLAGS = -I ../../../m4
if DO_SOMETHING
SUBDIRS = mjdgges kronecker bytecode libdynare++ gensylv qzcomplex block_kalman_filter sobol local_state_space_iterations
if COMPILE_LINSOLVE
SUBDIRS += linsolve
endif
if COMPILE_ORDSCHUR
SUBDIRS += ordschur
endif
if HAVE_MATIO
SUBDIRS += k_order_perturbation dynare_simul_
endif

26
mex/build/octave/configure.ac
View File

@ -32,20 +32,16 @@ if test "x$MKOCTFILE" != "x"; then
CXXFLAGS=`$MKOCTFILE -p CXXFLAGS`
LDFLAGS="`$MKOCTFILE -p LFLAGS` `$MKOCTFILE -p LDFLAGS`"
OCTAVE_VERSION=`$MKOCTFILE -v 2>&1 | sed 's/mkoctfile, version //'`
AX_COMPARE_VERSION([$OCTAVE_VERSION], [lt], [3.6], [AC_MSG_ERROR([Your Octave is too old, please upgrade to version 3.6 at least.])])
AX_COMPARE_VERSION([$OCTAVE_VERSION], [ge], [3.8], [OCTAVE38=yes])
AX_COMPARE_VERSION([$OCTAVE_VERSION], [ge], [4.0], [OCTAVE40=yes])
AX_COMPARE_VERSION([$OCTAVE_VERSION], [lt], [4.2], [AC_MSG_ERROR([Your Octave is too old, please upgrade to version 4.2 at least.])])
fi
AM_CONDITIONAL([COMPILE_LINSOLVE], [test "$OCTAVE38" != "yes"])
AM_CONDITIONAL([COMPILE_ORDSCHUR], [test "$OCTAVE40" != "yes"])
CFLAGS="$CFLAGS -Wall -Wno-parentheses"
FFLAGS="$FFLAGS -Wall"
CXXFLAGS="$CXXFLAGS -Wall -Wno-parentheses"
AC_PROG_F77([gfortran g77 f77])
AC_PROG_CC
AC_PROG_CC_C99 # mjdgges DLL now uses C99 features (variable declared in for loop)
AC_PROG_CXX
AC_PROG_RANLIB
AX_PROG_LN_S
@ -108,18 +104,6 @@ else
BUILD_MS_SBVAR_MEX_OCTAVE="no (missing GSL or MatIO library)"
fi
if test -n "$MKOCTFILE" -a "$OCTAVE38" != "yes"; then
BUILD_LINSOLVE_OCTAVE="yes"
else
BUILD_LINSOLVE_OCTAVE="no (Octave >= 3.8)"
fi
if test -n "$MKOCTFILE" -a "$OCTAVE40" != "yes"; then
BUILD_ORDSCHUR_OCTAVE="yes"
else
BUILD_ORDSCHUR_OCTAVE="no (Octave >= 4.0)"
fi
AC_ARG_ENABLE([openmp], AS_HELP_STRING([--enable-openmp], [use OpenMP for parallelization of some MEX files]), [
if test "x$enable_openmp" = "xyes"; then
CPPFLAGS="$CPPFLAGS -DUSE_OMP"
@ -137,8 +121,6 @@ Binaries (with "make"):
MS-SBVAR MEX files for Octave: $BUILD_MS_SBVAR_MEX_OCTAVE
Kalman Steady State MEX file for Octave: $BUILD_KALMAN_STEADY_STATE_OCTAVE
K-order and dynare_simul MEX for Octave: $BUILD_KORDER_DYNSIMUL_MEX_OCTAVE
Linsolve for Octave: $BUILD_LINSOLVE_OCTAVE
Ordschur for Octave: $BUILD_ORDSCHUR_OCTAVE
])
@ -151,12 +133,10 @@ AC_CONFIG_FILES([Makefile
k_order_perturbation/Makefile
dynare_simul_/Makefile
qzcomplex/Makefile
ordschur/Makefile
kalman_steady_state/Makefile
ms_sbvar/Makefile
block_kalman_filter/Makefile
sobol/Makefile
local_state_space_iterations/Makefile
linsolve/Makefile])
local_state_space_iterations/Makefile])
AC_OUTPUT

6
mex/build/octave/linsolve/Makefile.am
View File

@ -1,6 +0,0 @@
EXEEXT = .oct
include ../mex.am
mex_PROGRAMS = linsolve
nodist_linsolve_SOURCES = $(top_srcdir)/../../sources/linsolve/linsolve.cc

6
mex/build/octave/ordschur/Makefile.am
View File

@ -1,6 +0,0 @@
EXEEXT = .oct
include ../mex.am
mex_PROGRAMS = ordschur
nodist_ordschur_SOURCES = $(top_srcdir)/../../sources/ordschur/ordschur.cc

4
mex/sources/Makefile.am
View File

@ -8,13 +8,11 @@ EXTRA_DIST = \
bytecode \
qzcomplex \
k_order_perturbation \
ordschur \
kalman_steady_state \
ms-sbvar \
block_kalman_filter \
sobol \
local_state_space_iterations \
linsolve
local_state_space_iterations
clean-local:
rm -rf `find mex/sources -name *.o`

123
mex/sources/linsolve/linsolve.cc
View File

@ -1,123 +0,0 @@
/*
* Oct-file for bringing MATLAB's linsolve function to Octave.
*
* The implementation is incomplete:
* - it only knows about the TRANSA, LT, UT and SYM options
* - it only works with square matrices
* - it only works on double matrices (no single precision or complex)
*
* Written by Sébastien Villemot <sebastien@dynare.org>.
*/
/*
* Copyright (C) 2012-2017 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/>.
*/
#include <octave/oct.h>
#include <octave/ov-struct.h>
DEFUN_DLD(linsolve, args, nargout, "-*- texinfo -*-\n\
@deftypefn {Loadable Function} @var{x} = linsolve (@var{a}, @var{b})\n\
@deftypefnx {Loadable Function} [ @var{x}, @var{r} ] = linsolve (@var{a}, @var{b})\n\
@deftypefnx {Loadable Function} @var{x} = linsolve (@var{a}, @var{b}, @var{options})\n\
@deftypefnx {Loadable Function} [ @var{x}, @var{r} ] = linsolve (@var{a}, @var{b}, @var{options})\n\
\n\
Solves the linear system @math{A*X = B} and returns @var{X}.\n\
\n\
Alternatively, if @var{options} is provided and has a field @code{TRANSA} equal \
to @code{true}, then it solves the system @math{A'*X = B}.\n\
\n\
Also, the @code{LT} field of @var{options} (resp. the @code{UT} field) can be set \
to @code{true} to indicate that the matrix @var{a} is lower (resp. upper) \
triangular; similarly, the @code{SYM} field can be set to @code{true} to \
indicate that the matrix is symmetric.\n\
\n\
If requested, @var{r} will contain the reciprocal condition number.\n\
@end deftypefn\n\
")
{
int nargin = args.length();
octave_value_list retval;
if (nargin > 3 || nargin < 2 || nargout > 2)
{
print_usage();
return retval;
}
Matrix A = args(0).matrix_value();
Matrix B = args(1).matrix_value();
if (error_state)
return retval;
dim_vector dimA = A.dims();
dim_vector dimB = B.dims();
if (dimA(0) != dimB(0))
{
error("linsolve: must have same number of lines in A and B");
return retval;
}
if (dimA(0) != dimA(1))
{
error("linsolve: rectangular A not yet supported");
return retval;
}
MatrixType typA;
typA.mark_as_full();
bool transa = false;
if (nargin == 3)
{
octave_scalar_map opts = args(2).scalar_map_value();
if (error_state)
return retval;
octave_value tmp = opts.contents("TRANSA");
transa = tmp.is_defined() && tmp.bool_matrix_value().elem(0);
tmp = opts.contents("UT");
if (tmp.is_defined() && tmp.bool_matrix_value().elem(0))
typA.mark_as_upper_triangular();
tmp = opts.contents("LT");
if (tmp.is_defined() && tmp.bool_matrix_value().elem(0))
typA.mark_as_lower_triangular();
tmp = opts.contents("SYM");
if (tmp.is_defined() && tmp.bool_matrix_value().elem(0))
typA.mark_as_symmetric();
}
double rcond;
octave_idx_type info;
retval(0) = A.solve(typA, B, info, rcond, NULL, true, transa ? blas_trans : blas_no_trans);
if (nargout == 2)
{
if (dimB(1) > 0)
retval(1) = rcond;
else // If B has zero columns, A.solve() apparently does not compute A's rcond
retval(1) = A.rcond();
}
return retval;
}

81
mex/sources/mjdgges/mjdgges.c
View File

@ -32,17 +32,14 @@ my_criteria(const double *alphar, const double *alphai, const double *beta)
}
void
mjdgges(double *a, double *b, double *z, double *n, double *sdim, double *eval_r, double *eval_i, double zhreshold, double *info)
mjdgges(double *a, double *b, double *z, unsigned int n, double *sdim, double *alphar, double *alphai, double *beta, double *info)
{
lapack_int i_n, i_info, i_sdim, lwork;
double *alphar, *alphai, *beta, *work, *par, *pai, *pb, *per, *pei;
double *work;
double *junk;
lapack_int *bwork;
i_n = (lapack_int)*n;
alphar = mxCalloc(i_n, sizeof(double));
alphai = mxCalloc(i_n, sizeof(double));
beta = mxCalloc(i_n, sizeof(double));
i_n = (lapack_int) n;
lwork = 16*i_n+16;
work = mxCalloc(lwork, sizeof(double));
bwork = mxCalloc(i_n, sizeof(lapack_int));
@ -53,31 +50,6 @@ mjdgges(double *a, double *b, double *z, double *n, double *sdim, double *eval_r
*sdim = i_sdim;
*info = i_info;
par = alphar;
pai = alphai;
pb = beta;
pei = eval_i;
for (per = eval_r; per <= &eval_r[i_n-1]; ++per)
{
if ((fabs(*par) > zhreshold) || (fabs(*pb) > zhreshold))
*per = *par / *pb;
else
{
/* the ratio is too close to 0/0;
returns specific error number only if no other error */
if (i_info == 0)
*info = -30;
}
if (*pai == 0.0 && *pb == 0.0)
*pei = 0.0;
else
*pei = *pai / *pb;
++par;
++pai;
++pb;
++pei;
}
}
/* MATLAB interface */
@ -87,8 +59,7 @@ mexFunction(int nlhs, mxArray *plhs[],
{
unsigned int m1, n1, m2, n2;
double *s, *t, *z, *sdim, *eval_r, *eval_i, *info, *a, *b;
double n;
double *s, *t, *z, *sdim, *info, *a, *b;
/* Check for proper number of arguments */
@ -119,8 +90,12 @@ mexFunction(int nlhs, mxArray *plhs[],
t = mxGetPr(plhs[2]);
z = mxGetPr(plhs[3]);
sdim = mxGetPr(plhs[4]);
eval_r = mxGetPr(plhs[5]);
eval_i = mxGetPi(plhs[5]);
#if defined(MATLAB_MEX_FILE) && MATLAB_VERSION >= 0x0904
mxComplexDouble *gev = mxGetComplexDoubles(plhs[5]);
#else
double *gev_r = mxGetPr(plhs[5]);
double *gev_i = mxGetPi(plhs[5]);
#endif
info = mxGetPr(plhs[6]);
a = mxGetPr(prhs[0]);
@ -151,10 +126,40 @@ mexFunction(int nlhs, mxArray *plhs[],
memcpy(s, a, sizeof(double)*n1*n1);
memcpy(t, b, sizeof(double)*n1*n1);
n = n1;
double *alpha_r = mxCalloc(n1, sizeof(double));
double *alpha_i = mxCalloc(n1, sizeof(double));
double *beta = mxCalloc(n1, sizeof(double));
/* Do the actual computations in a subroutine */
mjdgges(s, t, z, &n, sdim, eval_r, eval_i, zhreshold, info);
mjdgges(s, t, z, n1, sdim, alpha_r, alpha_i, beta, info);
for (int i = 0; i < n1; i++)
{
if ((fabs(alpha_r[i]) > zhreshold) || (fabs(beta[i]) > zhreshold))
#if defined(MATLAB_MEX_FILE) && MATLAB_VERSION >= 0x0904
gev[i].real = alpha_r[i] / beta[i];
#else
gev_r[i] = alpha_r[i] / beta[i];
#endif
else
{
/* the ratio is too close to 0/0;
returns specific error number only if no other error */
if (*info == 0)
*info = -30;
}
if (alpha_i[i] == 0.0 && beta[i] == 0.0)
#if defined(MATLAB_MEX_FILE) && MATLAB_VERSION >= 0x0904
gev[i].imag = 0.0;
#else
gev_i[i] = 0.0;
#endif
else
#if defined(MATLAB_MEX_FILE) && MATLAB_VERSION >= 0x0904
gev[i].imag = alpha_i[i] / beta[i];
#else
gev_i[i] = alpha_i[i] / beta[i];
#endif
}
plhs[0] = mxCreateDoubleScalar(0);
}

117
mex/sources/ordschur/ordschur.cc
View File

@ -1,117 +0,0 @@
/*
* Oct-file for bringing MATLAB's ordschur function to Octave.
* Simple wrapper around LAPACK's dtrsen.
* Only supports real (double precision) decomposition.
* Only selection of eigenvalues with a boolean vector is supported.
*
* Written by Sébastien Villemot <sebastien@dynare.org>.
*/
/*
* Copyright (C) 2010-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/>.
*/
#include <octave/oct.h>
#include <octave/f77-fcn.h>
extern "C"
{
F77_RET_T
F77_FUNC(dtrsen, DTRSEN) (F77_CONST_CHAR_ARG_DECL, F77_CONST_CHAR_ARG_DECL,
const octave_idx_type *, const octave_idx_type &,
double *, const octave_idx_type &, double *, const octave_idx_type &,
double *, double *, octave_idx_type &, double &, double &, double *,
const octave_idx_type &, octave_idx_type *,
const octave_idx_type &, octave_idx_type &);
}
DEFUN_DLD(ordschur, args, nargout, "-*- texinfo -*-\n\
@deftypefn {Loadable Function} [ @var{us}, @var{ts} ] = ordschur (@var{u}, @var{t}, @var{select})\n\
\n\
Reorders the real Schur factorization @math{X = U*T*U'} so that selected\n\
eigenvalues appear in the upper left diagonal blocks of the quasi triangular\n\
Schur matrix @math{T}. The logical vector @var{select} specifies the selected\n\
eigenvalues as they appear along @math{T}'s diagonal.\n\
@end deftypefn\n\
")
{
int nargin = args.length();
octave_value_list retval;
if (nargin != 3 || nargout != 2)
{
print_usage();
return retval;
}
Matrix U = args(0).matrix_value();
Matrix T = args(1).matrix_value();
boolNDArray S = args(2).bool_array_value();
if (error_state)
return retval;
dim_vector dimU = U.dims();
dim_vector dimT = T.dims();
octave_idx_type n = dimU(0);
if (n != dimU(1) || n != dimT(0) || n != dimT(1))
{
error("ordschur: input matrices must be square and of same size");
return retval;
}
if (S.nelem() != n)
{
error("ordschur: selection vector has wrong size");
return retval;
}
octave_idx_type lwork = n, liwork = n;
OCTAVE_LOCAL_BUFFER(double, wr, n);
OCTAVE_LOCAL_BUFFER(double, wi, n);
OCTAVE_LOCAL_BUFFER(double, work, lwork);
OCTAVE_LOCAL_BUFFER(octave_idx_type, iwork, liwork);
octave_idx_type m, info;
double cond1, cond2;
OCTAVE_LOCAL_BUFFER(octave_idx_type, S2, n);
for (int i = 0; i < n; i++)
S2[i] = S(i);
F77_XFCN(dtrsen, dtrsen, (F77_CONST_CHAR_ARG("N"), F77_CONST_CHAR_ARG("V"),
S2, n, T.fortran_vec(), n, U.fortran_vec(), n,
wr, wi, m, cond1, cond2, work, lwork,
iwork, liwork, info));
if (info != 0)
{
error("ordschur: dtrsen failed");
return retval;
}
retval(0) = octave_value(U);
retval(1) = octave_value(T);
return retval;
}
/*
%!test
%! A = [1 2 3 -2; 4 5 6 -5 ; 7 8 9 -5; 10 11 12 4 ];
%! [U, T] = schur(A);
%! [US, TS] = ordschur(U, T, [ 0 0 1 1 ]);
%! assert(US*TS*US', A, sqrt(eps))
*/

2
preprocessor

@ -1 +1 @@
Subproject commit 3d946ec8b74aa18139a19f3512488886e2659a36
Subproject commit 1ba023e37d6ee7fc8d393840e42a4542274efbb7

2
tests/Makefile.am
View File

@ -42,6 +42,7 @@ MODFILES = \
estimation/MH_recover/fs2000_recover_2.mod \
estimation/MH_recover/fs2000_recover_3.mod \
estimation/t_proposal/fs2000_student.mod \
estimation/tune_mh_jscale/fs2000.mod \
moments/example1_var_decomp.mod \
moments/example1_bp_test.mod \
moments/test_AR1_spectral_density.mod \
@ -53,6 +54,7 @@ MODFILES = \
ramst_a.mod \
ramst_static_tag.mod \
example1.mod \
example1_on_the_fly.mod \
example2.mod \
example1_use_dll.mod \
example1_with_tags.mod \

105
tests/estimation/tune_mh_jscale/fs2000.mod Normal file
View File

@ -0,0 +1,105 @@
/*
* Copyright (C) 2004-2018 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/>.
*/
var m P c e W R k d n l gy_obs gp_obs y dA;
varexo e_a e_m;
parameters alp bet gam mst rho psi del;
alp = 0.33;
bet = 0.99;
gam = 0.003;
mst = 1.011;
rho = 0.7;
psi = 0.787;
del = 0.02;
model;
dA = exp(gam+e_a);
log(m) = (1-rho)*log(mst) + rho*log(m(-1))+e_m;
-P/(c(+1)*P(+1)*m)+bet*P(+1)*(alp*exp(-alp*(gam+log(e(+1))))*k^(alp-1)*n(+1)^(1-alp)+(1-del)*exp(-(gam+log(e(+1)))))/(c(+2)*P(+2)*m(+1))=0;
W = l/n;
-(psi/(1-psi))*(c*P/(1-n))+l/n = 0;
R = P*(1-alp)*exp(-alp*(gam+e_a))*k(-1)^alp*n^(-alp)/W;
1/(c*P)-bet*P*(1-alp)*exp(-alp*(gam+e_a))*k(-1)^alp*n^(1-alp)/(m*l*c(+1)*P(+1)) = 0;
c+k = exp(-alp*(gam+e_a))*k(-1)^alp*n^(1-alp)+(1-del)*exp(-(gam+e_a))*k(-1);
P*c = m;
m-1+d = l;
e = exp(e_a);
y = k(-1)^alp*n^(1-alp)*exp(-alp*(gam+e_a));
gy_obs = dA*y/y(-1);
gp_obs = (P/P(-1))*m(-1)/dA;
end;
shocks;
var e_a; stderr 0.014;
var e_m; stderr 0.005;
end;
steady_state_model;
dA = exp(gam);
gst = 1/dA;
m = mst;
khst = ( (1-gst*bet*(1-del)) / (alp*gst^alp*bet) )^(1/(alp-1));
xist = ( ((khst*gst)^alp - (1-gst*(1-del))*khst)/mst )^(-1);
nust = psi*mst^2/( (1-alp)*(1-psi)*bet*gst^alp*khst^alp );
n = xist/(nust+xist);
P = xist + nust;
k = khst*n;
l = psi*mst*n/( (1-psi)*(1-n) );
c = mst/P;
d = l - mst + 1;
y = k^alp*n^(1-alp)*gst^alp;
R = mst/bet;
W = l/n;
ist = y-c;
q = 1 - d;
e = 1;
gp_obs = m/dA;
gy_obs = dA;
end;
steady;
check;
estimated_params;
alp, beta_pdf, 0.356, 0.02;
bet, beta_pdf, 0.993, 0.002;
gam, normal_pdf, 0.0085, 0.003;
mst, normal_pdf, 1.0002, 0.007;
rho, beta_pdf, 0.129, 0.100;
psi, beta_pdf, 0.65, 0.05;
del, beta_pdf, 0.01, 0.005;
stderr e_a, inv_gamma_pdf, 0.035449, inf;
stderr e_m, inv_gamma_pdf, 0.008862, inf;
end;
varobs gp_obs gy_obs;
estimation(order=1, datafile='../fsdat_simul', nobs=192, loglinear, mh_replic=20000, mh_nblocks=2, mh_tune_jscale=0.33);
mhdata = load('fs2000/metropolis/fs2000_mh_history_0.mat');
if any(abs(mhdata.record.AcceptanceRatio-options_.mh_tune_jscale.target)>options_.mh_tune_jscale.c1)
error('Automagic tuning of the MCMC proposal scale parameter did not work as expected!')
end

39
tests/example1_on_the_fly.mod Normal file
View File

@ -0,0 +1,39 @@
// Example 1 from Collard's guide to Dynare
model;
c|e*theta|p*h|e^(1+psi|p)=(1-alpha)*y;
k|e = beta|p*(((exp(b)*c)/(exp(b(+1))*c(+1)))
*(exp(b(+1))*alpha|p*y(+1)+(1-delta)*k));
y|e = exp(a)*(k(-1)^alpha)*(h^(1-alpha));
k = exp(b)*(y-c)+(1-delta|p)*k(-1);
a|e = rho|p*a(-1)+tau*b(-1) + e|x;
b|e = tau|p*a(-1)+rho*b(-1) + u|x;
end;
alpha = 0.36;
rho = 0.95;
tau = 0.025;
beta = 0.99;
delta = 0.025;
psi = 0;
theta = 2.95;
phi = 0.1;
initval;
y = 1.08068253095672;
c = 0.80359242014163;
h = 0.29175631001732;
k = 11.08360443260358;
a = 0;
b = 0;
e = 0;
u = 0;
end;
shocks;
var e; stderr 0.009;
var u; stderr 0.009;
var e, u = phi*0.009*0.009;
end;
stoch_simul;

90
tests/julia/rbc/test1.jl
View File

@ -3,90 +3,78 @@ if isempty(findin([abspath("../../../julia")], LOAD_PATH))
unshift!(LOAD_PATH, abspath("../../../julia"))
end
# Load Dynare package
importall Dynare
using Dynare
using SteadyState
# Compile the rbc.mod file -> produce a module with the model definition.
@dynare "rbc1.mod"
importall SteadyState
# First call to the steady state routine (analytical)
@time SteadyState.steady!(model_, oo_)
# First call to the steady state routine (analytical)
@time SteadyState.steady!(model_, oo_)
paramsinit = copy(model_.params);
steadyState = deepcopy(oo_.steady_state)
yinit = deepcopy(oo_.steady_state)
y_init = copy(yinit)
y_init[1] = 1.1*yinit[1]
y_init[2] = 0.9*yinit[2]
y_init = copy(steadyState)
y_init[1] = 1.1*steadyState[1]
y_init[2] = 0.9*steadyState[2]
# First call to the steady state routine (numerical)
println("First call to the numerical steady state routine")
@time SteadyState.steady!(model_, oo_, yinit)
@time SteadyState.steady!(model_, oo_, y_init)
# Check results
@assert maximum(abs(oo_.steady_state-yinit))<1e-9
@assert maximum(abs.(oo_.steady_state-steadyState))<1e-6
y_init = copy(yinit)
yinit[1] = 1.1*yinit[1]
yinit[2] = 0.9*yinit[2]
yinit = deepcopy(steadyState)
yinit[1] = 1.1*steadyState[1]
yinit[2] = 0.9*steadyState[2]
# Second call to the steady state routine (numerical)
println("Second call to the numerical steady state routine")
@time SteadyState.steady!(model_, oo_, yinit)
params = model_.params
# change alpha
println("Change α")
model_.params[4] = max(min(1.0, model_.params[4]*1.2), 0.0)
ys = SteadyState.steady(model_, oo_)
y_init = copy(yinit)
@time SteadyState.steady!(model_, oo_, y_init)
y_init = copy(yinit)
@time SteadyState.steady!(model_, oo_, y_init)
@assert maximum(abs(oo_.steady_state-ys))<1e-9
model_.params = deepcopy(params)
model_.params[4] = max(min(1.0, params[4]*1.1), 0.0)
@time ys = SteadyState.steady(model_, oo_) # Analytical steady state
@time SteadyState.steady!(model_, oo_, steadyState)
@assert maximum(abs.(oo_.steady_state-ys))<1e-6
# change delta
println("Change δ")
model_.params[6] = max(min(1.0, model_.params[6]*1.2), 0.0)
ys = SteadyState.steady(model_, oo_)
y_init = copy(yinit)
@time SteadyState.steady!(model_, oo_, y_init)
y_init = copy(yinit)
@time SteadyState.steady!(model_, oo_, y_init)
@assert maximum(abs(oo_.steady_state-ys))<1e-9
model_.params = deepcopy(params)
model_.params[6] = max(min(1.0, params[6]*1.1), 0.0)
@time ys = SteadyState.steady(model_, oo_) # Analytical steady state
@time SteadyState.steady!(model_, oo_, steadyState)
@assert maximum(abs.(oo_.steady_state-ys))<1e-6
# change beta
println("Change β")
model_.params[1] = max(min(1-1e-6, model_.params[1]*0.99), 0.0)
ys = SteadyState.steady(model_, oo_)
y_init = copy(yinit)
@time SteadyState.steady!(model_, oo_, y_init)
y_init = copy(yinit)
@time SteadyState.steady!(model_, oo_, y_init)
@assert maximum(abs(oo_.steady_state-ys))<1e-9
model_.params = deepcopy(params)
model_.params[1] = max(min(1-1e-6, params[1]*0.99), 0.0)
@time ys = SteadyState.steady(model_, oo_) # Analytical steady state
@time SteadyState.steady!(model_, oo_, steadyState)
@assert maximum(abs.(oo_.steady_state-ys))<1e-6
# change tau
println("Change τ")
model_.params[3] /= 1.5
ys = SteadyState.steady(model_, oo_)
y_init = copy(yinit)
@time SteadyState.steady!(model_, oo_, y_init)
y_init = copy(yinit)
@time SteadyState.steady!(model_, oo_, y_init)
@assert maximum(abs(oo_.steady_state-ys))<1e-9
model_.params = deepcopy(params)
model_.params[3] = params[3]/1.1
@time ys = SteadyState.steady(model_, oo_)
@time SteadyState.steady!(model_, oo_, steadyState)
@assert maximum(abs.(oo_.steady_state-ys))<1e-6
# change Epsilon
println("Change ϵ")
model_.params[5] *= 1.5
ys = SteadyState.steady(model_, oo_)
y_init = copy(yinit)
@time SteadyState.steady!(model_, oo_, y_init)
y_init = copy(yinit)
@time SteadyState.steady!(model_, oo_, y_init)
@assert maximum(abs(oo_.steady_state-ys))<1e-9
model_.params = deepcopy(params)
model_.params[5] = params[5]*1.1
@time ys = SteadyState.steady(model_, oo_)
@time SteadyState.steady!(model_, oo_, steadyState)
@assert maximum(abs.(oo_.steady_state-ys))<1e-6

49
tests/julia/rbc/test1pfm.jl
View File

@ -3,18 +3,20 @@ if isempty(findin([abspath("../../../julia")], LOAD_PATH))
unshift!(LOAD_PATH, abspath("../../../julia"))
end
# Load Dynare package
importall Dynare
using PyPlot
# Load packages
using Plots
using Dynare
using SteadyState
using PerfectForesightModelSolver
plotlyjs() # Choose a backend
# Compile the rbc.mod file -> produce a module with the model definition.
@dynare "rbc1.mod"
importall SteadyState
importall PerfectForesightModelSolver
# First call to the steady state routine (analytical)
@time SteadyState.steady!(model_, oo_)
println(oo_.steady_state)
@ -32,10 +34,9 @@ exogenousvariables = repmat(oo_.exo_steady_state', options_.pfmsolver.periods+2,
n = 200
dates = collect(0:n-1)
plt[:figure](1)
#plt[:figure](1)
plot(dates, vec(endogenousvariables[1,1:n]), color="black", linewidth=2.0, linestyle="-")
# Initialize paths for the endogenous variables
endogenousvariables = repmat(oo_.steady_state, 1, options_.pfmsolver.periods+2)
# Destroy part of the initial stock of physical capital...
@ -65,19 +66,19 @@ exogenousvariables[10+1, 1] = 2
# Simulate the paths for the endogenous variables, given the expected shock
@time simulate_perfect_foresight_model!(endogenousvariables, exogenousvariables, oo_.steady_state, model_, options_)
n = 200
dates = collect(0:n-1)
plt[:figure](2)
subplot(221)
title("Efficiency")
plot(dates, vec(endogenousvariables[5,1:n]), color="black", linewidth=2.0, linestyle="-")
subplot(223)
title("Output")
plot(dates, vec(endogenousvariables[2,1:n]), color="black", linewidth=2.0, linestyle="-")
subplot(222)
title("Consumption")
plot(dates, vec(endogenousvariables[4,1:n]), color="black", linewidth=2.0, linestyle="-")
subplot(224)
title("Labour")
plot(dates, vec(endogenousvariables[3,1:n]), color="black", linewidth=2.0, linestyle="-")
suptitle("Expected positive expected shock")
# dates = collect(0:n-1)
# n = 200
# plt[:figure](2)
# subplot(221)
# title("Efficiency")
# plot(dates, vec(endogenousvariables[5,1:n]), color="black", linewidth=2.0, linestyle="-")
# subplot(223)
# title("Output")
# plot(dates, vec(endogenousvariables[2,1:n]), color="black", linewidth=2.0, linestyle="-")
# subplot(222)
# title("Consumption")
# plot(dates, vec(endogenousvariables[4,1:n]), color="black", linewidth=2.0, linestyle="-")
# subplot(224)
# title("Labour")
# plot(dates, vec(endogenousvariables[3,1:n]), color="black", linewidth=2.0, linestyle="-")
# suptitle("Expected positive expected shock")

14
tests/julia/rbc/test2.jl
View File

@ -3,14 +3,12 @@ if isempty(findin([abspath("../../../julia")], LOAD_PATH))
unshift!(LOAD_PATH, abspath("../../../julia"))
end
# Load Dynare package
importall Dynare
# Load Dynare package(s)
using Dynare
using SteadyState
# Compile the rbc.mod file -> produce a module with the model definition.
@dynare "rbc1.mod"
@dynare "rbc2.mod"
importall SteadyState
steady!(model_, oo_)
# Compute the steady state
steady_state!(model_, oo_)

1
tests/julia/rbc/test3.jl
View File

@ -17,4 +17,5 @@ importall Dynare
# using rbc1
print(rbc1.model_.fname)
print(rbc2.model_.fname)

4
tests/stochastic-backward-models/solow_cd.mod
View File

@ -54,6 +54,4 @@ shocks;
var e_n = 0.001;
end;
initialconditions = dseries([1 1.02 1 1.02 1], 2000Q1, {'Efficiency'; 'EfficiencyGrowth'; 'Population'; 'PopulationGrowth'; 'PhysicalCapitalStock'});
simulations = simul_backward_model(initialconditions, 5000);
simulations = simul_backward_model([], 5000);

16
tests/stochastic-backward-models/solow_cd_with_steadystate.mod
View File

@ -35,6 +35,14 @@ end;
d = dseries([.5 1 1 1.04 15], 2000Q1, {'Efficiency'; 'EfficiencyGrowth'; 'Population'; 'PopulationGrowth'; 'PhysicalCapitalStock'});
histval;
Efficiency(0) = .5;
EfficiencyGrowth(0) = 1;
Population(0) = 1;
PopulationGrowth(0) = 1.04;
PhysicalCapitalStock(0) = 15;
end;
LongRunEfficiency = d.Efficiency*d.EfficiencyGrowth^(rho_x/(1-rho_x));
LongRunPopulation = d.Population*d.PopulationGrowth^(rho_n/(1-rho_n));
LongRunEfficiencyGrowth = EfficiencyGrowth_ss;
@ -46,18 +54,18 @@ T = 5*floor(log(precision)/log(max(rho_x, rho_n)));
e = dseries(zeros(T, 2), 2000Q2, {'e_x'; 'e_n'});
simulations = simul_backward_model(d, T, e);
simulations = simul_backward_model([], T, e);
if abs(simulations.Efficiency.data(end)-LongRunEfficiency.data)>1e-10
error('Wrong long run level!')
error('Wrong long run level (Efficiency, diff=%s)!', num2str(abs(simulations.Efficiency.data(end)-LongRunEfficiency.data)))
end
if abs(simulations.Population.data(end)-LongRunPopulation.data)>1e-10
error('Wrong long run level!')
error('Wrong long run level (Population, diff=%s)!', num2str(abs(simulations.Population.data(end)-LongRunPopulation.data)))
end
IntensiveCapitalStock = simulations.PhysicalCapitalStock/(simulations.Efficiency*simulations.Population);
if abs(IntensiveCapitalStock.data(end)-LongRunIntensiveCapitalStock)>1e-3
error('Wrong long run level!')
error('Wrong long run level (Intensive capital stock)!')
end

10
tests/stochastic-backward-models/solow_ces.mod
View File

@ -35,11 +35,17 @@ model;
PhysicalCapitalStock = (1-delta)*PhysicalCapitalStock(-1) + s*Output;
end;
d = dseries([1 1.02 1 1.02 1], 2000Q1, {'Efficiency'; 'EfficiencyGrowth'; 'Population'; 'PopulationGrowth'; 'PhysicalCapitalStock'});
histval;
Efficiency(0) = 1;
EfficiencyGrowth(0) = 1.02;
Population(0) = 1;
PopulationGrowth(0) = 1.02;
PhysicalCapitalStock(0) = 1;
end;
shocks;
var e_x = 0.005;
var e_n = 0.001;
end;
simulations = simul_backward_model(d, 5000);
simulations = simul_backward_model([], 5000);

5
windows/dynare.nsi
View File

@ -103,6 +103,11 @@ Section "MEX files for MATLAB 64-bit, version 7.8 to 9.3 (R2009a to R2017b)"
File ..\mex\matlab\win64-7.8-9.3\*.mexw64
SectionEnd
Section "MEX files for MATLAB 64-bit, version 9.4 (R2018a)"
SetOutPath $INSTDIR\mex\matlab\win64-9.4
File ..\mex\matlab\win64-9.4\*.mexw64
SectionEnd
SectionGroupEnd
SectionGroup "MEX files for OCTAVE"