Merge remote-tracking branch 'local_master/master' into dr1break

time-shift
Michel Juillard 2011-12-31 10:10:41 +01:00
commit 80ca47d62a
78 changed files with 16763 additions and 900 deletions

View File

@ -4,7 +4,9 @@ info_TEXINFOS = dynare.texi
if HAVE_TEXI2HTML
if HAVE_LATEX2HTML
html-local: dynare.texi
html-local: dynare.html
dynare.html: dynare.texi
rm -rf dynare.html
mkdir -p dynare.html
cd dynare.html && $(TEXI2HTML) --l2h --split section --prefix index ../dynare.texi

View File

@ -1954,6 +1954,12 @@ of the simulation is period @code{0}, then period @code{-1}, and so on.
If your lagged variables are linked by identities, be careful to
satisfy these identities when you set historical initial values.
Variables not initialized in the @code{histval} block are assumed to
have a value of zero at period 0 and before. Note that this behavior
differs from the case where there is no @code{histval} block, where all
variables are initialized at their steady state value at period 0 and
before.
@examplehead
@example

View File

@ -340,6 +340,24 @@ License: GPL-2+
along with this program. If not, see
<http://www.gnu.org/licenses/>.
Files: mex/sources/sobol/sobol.hh mex/sources/sobol/initialize_v_array.hh
mex/sources/sobol/initialize_v_array.inc
Copyright: 2009 John Burkardt
2010-2011 Dynare Team
License: LGPL-3+
This program is free software: you can redistribute it and/or modify
it under the terms of the GNU Lesser General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
.
This program 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 Lesser General Public License for more details.
.
You should have received a copy of the GNU Lesser General Public License
along with this program. If not, see <http://www.gnu.org/licenses/>.
Files: ms-sbvar/utilities_dw/*
Copyright: 1996-2011 Daniel Waggoner
License: GPL-3+

View File

@ -1,9 +1,9 @@
function [fval,exit_flag,ys,trend_coeff,info,Model,DynareOptions,BayesInfo,DynareResults,DLIK,AHess] = DsgeLikelihood(xparam1,DynareDataset,DynareOptions,Model,EstimatedParameters,BayesInfo,DynareResults,derivatives_info)
function [fval,exit_flag,ys,trend_coeff,info,Model,DynareOptions,BayesInfo,DynareResults,DLIK,AHess] = dsge_likelihood(xparam1,DynareDataset,DynareOptions,Model,EstimatedParameters,BayesInfo,DynareResults,derivatives_info)
% Evaluates the posterior kernel of a dsge model.
%@info:
%! @deftypefn {Function File} {[@var{fval},@var{exit_flag},@var{ys},@var{trend_coeff},@var{info},@var{Model},@var{DynareOptions},@var{BayesInfo},@var{DynareResults},@var{DLIK},@var{AHess}] =} DsgeLikelihood (@var{xparam1},@var{DynareDataset},@var{DynareOptions},@var{Model},@var{EstimatedParameters},@var{BayesInfo},@var{DynareResults},@var{derivatives_flag})
%! @anchor{DsgeLikelihood}
%! @deftypefn {Function File} {[@var{fval},@var{exit_flag},@var{ys},@var{trend_coeff},@var{info},@var{Model},@var{DynareOptions},@var{BayesInfo},@var{DynareResults},@var{DLIK},@var{AHess}] =} dsge_likelihood (@var{xparam1},@var{DynareDataset},@var{DynareOptions},@var{Model},@var{EstimatedParameters},@var{BayesInfo},@var{DynareResults},@var{derivatives_flag})
%! @anchor{dsge_likelihood}
%! @sp 1
%! Evaluates the posterior kernel of a dsge model.
%! @sp 2
@ -229,7 +229,7 @@ if EstimatedParameters.ncn
% Try to compute the cholesky decomposition of H (possible iff H is positive definite)
[CholH,testH] = chol(H);
if testH
% The variance-covariance matrix of the structural innovations is not definite positive. We have to compute the eigenvalues of this matrix in order to build the endogenous penalty.
% The variance-covariance matrix of the measurement errors is not definite positive. We have to compute the eigenvalues of this matrix in order to build the endogenous penalty.
a = diag(eig(H));
k = find(a < 0);
if k > 0
@ -428,13 +428,13 @@ switch DynareOptions.lik_init
[err,Pstar] = kalman_steady_state(transpose(T),R*Q*transpose(R),transpose(build_selection_matrix(Z,np,length(Z))),H);
end
if err
disp(['DsgeLikelihood:: I am not able to solve the Riccati equation, so I switch to lik_init=1!']);
disp(['dsge_likelihood:: I am not able to solve the Riccati equation, so I switch to lik_init=1!']);
DynareOptions.lik_init = 1;
Pstar = lyapunov_symm(T,R*Q*R',DynareOptions.qz_criterium,DynareOptions.lyapunov_complex_threshold);
end
Pinf = [];
otherwise
error('DsgeLikelihood:: Unknown initialization approach for the Kalman filter!')
error('dsge_likelihood:: Unknown initialization approach for the Kalman filter!')
end
if analytic_derivation

View File

@ -61,10 +61,10 @@ switch task
if horizon == 0
horizon = 5;
end
if size(oo_.endo_simul,2) < maximum_lag
if isempty(M_.endo_histval)
y0 = repmat(oo_.steady_state,1,maximum_lag);
else
y0 = oo_.endo_simul(:,1:maximum_lag);
y0 = M_.endo_histval;
end
case 'smoother'
horizon = options_.forecast;

View File

@ -37,6 +37,9 @@ if nargin && ~isempty(path_to_dynare)
end
dynareroot = strrep(which('dynare'),'dynare.m','');
origin = pwd();
cd([dynareroot '/..'])
if ~nargin || nargin==1
verbose = 1;
end
@ -67,7 +70,8 @@ if ~exist('OCTAVE_VERSION')
addpath([dynareroot '/missing/rows_columns'])
% Replacement for vec() (inexistent under MATLAB)
addpath([dynareroot '/missing/vec'])
if isempty(license('inuse','statistics_toolbox'))
[has_statistics_toolbox junk] = license('checkout','statistics_toolbox');
if ~has_statistics_toolbox
% Replacements for functions of the stats toolbox
addpath([dynareroot '/missing/stats/'])
end
@ -97,7 +101,8 @@ if exist('OCTAVE_VERSION')
addpath([dynareroot '/missing/nanmean'])
end
else
if isempty(license('inuse','statistics_toolbox'))
[has_statistics_toolbox junk] = license('checkout','statistics_toolbox');
if ~has_statistics_toolbox
addpath([dynareroot '/missing/nanmean'])
end
end
@ -246,5 +251,17 @@ else
end
if verbose
disp([ message 'k-order solution simulation.' ])
end
% Test if qmc_sequence DLL is present
if exist('qmc_sequence', 'file') == 3
message = '[mex] ';
else
message = '[no] ';
end
if verbose
disp([ message 'Quasi Monte-Carlo sequence (Sobol).' ])
disp(' ')
end
cd(origin);

View File

@ -32,7 +32,7 @@ function dynare_estimation_1(var_list_,dname)
global M_ options_ oo_ estim_params_ bayestopt_ dataset_
if ~options_.dsge_var
objective_function = str2func('DsgeLikelihood');
objective_function = str2func('dsge_likelihood');
else
objective_function = str2func('DsgeVarLikelihood');
end
@ -248,7 +248,7 @@ if ~isequal(options_.mode_compute,0) && ~options_.mh_posterior_mode_estimation
fval = feval(objective_function,xparam1,dataset_,options_,M_,estim_params_,bayestopt_,oo_);
options_.mh_jscale = Scale;
mouvement = max(max(abs(PostVar-OldPostVar)));
fval = DsgeLikelihood(xparam1,dataset_,options_,M_,estim_params_,bayestopt_,oo_);
fval = dsge_likelihood(xparam1,dataset_,options_,M_,estim_params_,bayestopt_,oo_);
disp(['Change in the covariance matrix = ' num2str(mouvement) '.'])
disp(['Mode improvement = ' num2str(abs(OldMode-fval))])
OldMode = fval;

View File

@ -41,7 +41,7 @@ function [A,B,ys,info,Model,DynareOptions,DynareResults] = dynare_resolve(Model,
%! @sp 2
%! @strong{This function is called by:}
%! @sp 1
%! @ref{DsgeLikelihood}, @ref{DsgeLikelihood_hh}, @ref{DsgeVarLikelihood}, @ref{dsge_posterior_kernel}, @ref{DsgeSmoother}, @ref{dynare_sensitivity}, @ref{gsa/thet2tau}, @ref{gsa/stab_map}, @ref{identification_analysis}, @ref{imcforecast}, @ref{thet2tau}
%! @ref{dsge_likelihood}, @ref{DsgeLikelihood_hh}, @ref{DsgeVarLikelihood}, @ref{dsge_posterior_kernel}, @ref{DsgeSmoother}, @ref{dynare_sensitivity}, @ref{gsa/thet2tau}, @ref{gsa/stab_map}, @ref{identification_analysis}, @ref{imcforecast}, @ref{thet2tau}
%! @sp 2
%! @strong{This function calls:}
%! @sp 1

View File

@ -38,8 +38,11 @@ global options_
options_ = set_default_option(options_,'solve_algo',2);
info = 0;
if options_.solve_algo == 0
if ~exist('OCTAVE_VERSION') && isempty(license('inuse','optimization_toolbox'))
error('You can''t use solve_algo=0 since you don''t have MATLAB''s Optimization Toolbox')
if ~exist('OCTAVE_VERSION')
[has_optimization_toolbox junk] = license('checkout','optimization_toolbox');
if ~has_optimization_toolbox
error('You can''t use solve_algo=0 since you don''t have MATLAB''s Optimization Toolbox')
end
end
options=optimset('fsolve');
options.MaxFunEvals = 50000;

View File

@ -33,6 +33,7 @@ function time_series = extended_path(initial_conditions,sample_size)
global M_ options_ oo_
debug = 0;
options_.verbosity = options_.ep.verbosity;
verbosity = options_.ep.verbosity+debug;
% Test if bytecode and block options are used (these options are mandatory)
@ -102,16 +103,14 @@ switch options_.ep.innovation_distribution
end
% Set future shocks (Stochastic Extended Path approach)
if options_.ep.stochastic
[r,w] = gauss_hermite_weights_and_nodes(options_.ep.number_of_nodes);
switch options_.ep.stochastic
case 1
if M_.exo_nbr>1
rr = cell(1);
ww = cell(1);
for i=1:size(M_.Sigma_e,1)
rr = {r};
ww = {w};
if options_.ep.stochastic.status
switch options_.ep.stochastic.method
case 'tensor'
[r,w] = gauss_hermite_weights_and_nodes(options_.ep.stochastic.nodes);
if options_.ep.stochastic.order*M_.exo_nbr>1
for i=1:options_.ep.stochastic.order*M_.exo_nbr
rr(i) = {r};
ww(i) = {w};
end
rrr = cartesian_product_of_sets(rr{:});
www = cartesian_product_of_sets(ww{:});
@ -120,9 +119,25 @@ if options_.ep.stochastic
www = w;
end
www = prod(www,2);
nnn = length(www);
number_of_nodes = length(www);
relative_weights = www/max(www);
switch options_.ep.stochastic.pruned.status
case 1
jdx = find(relative_weights>options_.ep.stochastic.pruned.relative);
www = www(jdx);
www = www/sum(www);
rrr = rrr(jdx,:);
case 2
jdx = find(weights>options_.ep.stochastic.pruned.level);
www = www(jdx);
www = www/sum(www);
rrr = rrr(jdx,:);
otherwise
% Nothing to be done!
end
nnn = length(www);
otherwise
error(['Order ' int2str(options_.ep.stochastic) ' Stochastic Extended Path method is not implemented!'])
error('extended_path:: Unknown stochastic_method option!')
end
else
rrr = zeros(1,number_of_structural_innovations);
@ -133,7 +148,6 @@ end
% Initializes some variables.
t = 0;
% Set waitbar (graphic or text mode)
graphic_waitbar_flag = ~( options_.console_mode || exist('OCTAVE_VERSION') );
@ -169,8 +183,10 @@ while (t<sample_size)
% Put it in oo_.exo_simul (second line).
oo_.exo_simul(2,positive_var_indx) = shocks;
for s = 1:nnn
oo_.exo_simul(3,positive_var_indx) = rrr(s,:)*covariance_matrix_upper_cholesky;
if options_.ep.init && s==1% Compute first order solution. t==1 &&
for u=1:options_.ep.stochastic.order
oo_.exo_simul(2+u,positive_var_indx) = rrr(s,(((u-1)*M_.exo_nbr)+1):(u*M_.exo_nbr))*covariance_matrix_upper_cholesky;
end
if options_.ep.init% Compute first order solution...
initial_path = simult_(initial_conditions,oo_.dr,oo_.exo_simul(2:end,:),1);
if options_.ep.init==1
oo_.endo_simul(:,1:end-1) = initial_path(:,1:end-1);% Last column is the steady state.

View File

@ -85,6 +85,6 @@ end
pshape_original = bayestopt_.pshape;
bayestopt_.pshape = Inf(size(bayestopt_.pshape));
llik = -DsgeLikelihood(parameters,dataset,options_,M_,estim_params_,bayestopt_,oo_);
llik = -dsge_likelihood(parameters,dataset,options_,M_,estim_params_,bayestopt_,oo_);
bayestopt_.pshape = pshape_original;

View File

@ -28,13 +28,15 @@ function global_initialization()
% You should have received a copy of the GNU General Public License
% along with Dynare. If not, see <http://www.gnu.org/licenses/>.
global oo_ M_ options_ estim_params_ bayestopt_
global oo_ M_ options_ estim_params_ bayestopt_ estimation_info
estim_params_ = [];
bayestopt_ = [];
options_.console_mode = 0;
options_.verbosity = 1;
options_.terminal_condition = 0;
options_.rplottype = 0;
options_.smpl = 0;
@ -54,7 +56,6 @@ options_.lyapunov_complex_threshold = 1e-15;
options_.solve_tolf = eps^(1/3);
options_.solve_tolx = eps^(2/3);
options_.solve_maxit = 500;
options_.deterministic_simulation_initialization = 0;
% Default number of threads for parallelized mex files.
options_.threads.kronecker.A_times_B_kronecker_C = 1;
@ -109,34 +110,76 @@ options_.SpectralDensity = 0;
% Extended path options
%
% Set verbose mode
options_.ep.verbosity = 1;
ep.verbosity = 0;
% Initialization of the perfect foresight equilibrium paths
% * init=0, previous solution is used.
% * init=1, a path generated with the first order reduced form is used.
% * init=2, mix of cases 0 and 1.
options_.ep.init = 0;
ep.init = 0;
% Maximum number of iterations for the deterministic solver.
options_.ep.maxit = 500;
ep.maxit = 500;
% Number of periods for the perfect foresight model.
options_.ep.periods = 200;
ep.periods = 200;
% Default step for increasing the number of periods if needed
options_.ep.step = 50;
ep.step = 50;
% Define last periods used to test if the solution is stable with respect to an increase in the number of periods.
options_.ep.lp = 5;
ep.lp = 5;
% Define first periods used to test if the solution is stable with respect to an increase in the number of periods.
options_.ep.fp = 100;
ep.fp = 100;
% Define the distribution for the structural innovations.
options_.ep.innovation_distribution = 'gaussian';
ep.innovation_distribution = 'gaussian';
% Set flag for the seed
options_.ep.set_dynare_seed_to_default = 1;
ep.set_dynare_seed_to_default = 1;
% Set algorithm for the perfect foresight solver
options_.ep.stack_solve_algo = 4;
% Set method
options_.ep.stochastic = 0;
% Set number of nodes for future shocks
options_.ep.number_of_nodes = 5;
ep.stack_solve_algo = 4;
% Stochastic extended path related options.
ep.stochastic.status = 0;
ep.stochastic.method = 'tensor';
ep.stochastic.order = 1;
ep.stochastic.nodes = 5;
ep.stochastic.pruned.status = 0;
ep.stochastic.pruned.relative = 1e-5;
ep.stochastic.pruned.level = 1e-5;
% Copy ep structure in options_ global structure
options_.ep = ep;
% Particle filter
%
% Default is that we do not use the non linear kalman filter
particle.status = 0;
% How do we initialize the states?
particle.initialization = 1;
particle_filter.initial_state_prior_std = .0001;
% Set the default order of approximation of the model (perturbation).
particle_filter.perturbation = 2;
% Set the default number of particles.
particle_filter.number_of_particles = 500;
% Set the default approximation order (Smolyak)
particle_filter.smolyak_accuracy = 3;
% By default we don't use pruning
particle_filter.pruning = 0;
% Set default algorithm
particle_filter.algorithm = 'sequential_importance_particle_filter';
% Set the Gaussian approximation method
particle_filter.approximation_method = 'unscented';
% Set unscented parameters alpha, beta and kappa for gaussian approximation
particle_filter.unscented.alpha = 1 ;
particle_filter.unscented.beta = 2 ;
particle_filter.unscented.kappa = 1 ;
% Configuration of resampling in case of particles
particle_filter.resampling = 'systematic' ;
% Choice of the resampling method
particle_filter.resampling_method = 'traditional' ;
% Configuration of the mixture filters
particle_filter.mixture_method = 'particles' ;
% Size of the different mixtures
particle_filter.mixture_state_variables = 5 ;
particle_filter.mixture_structural_shocks = 1 ;
particle_filter.mixture_measurement_shocks = 1 ;
% Copy ep structure in options_ global structure
options_.particle = particle;
% TeX output
options_.TeX = 0;
@ -182,6 +225,35 @@ options_.ramsey_policy = 0;
options_.timeless = 0;
% estimation
estimation_info.prior = struct('name', {}, 'shape', {}, 'mean', {}, ...
'mode', {}, 'stdev', {}, 'date1', {}, ...
'date2', {}, 'shift', {}, 'variance', {});
estimation_info.structural_innovation.prior = struct('name', {}, 'shape', {}, 'mean', {}, ...
'mode', {}, 'stdev', {}, 'date1', {}, ...
'date2', {}, 'shift', {}, 'variance', {});
estimation_info.structural_innovation_corr.prior = struct('name', {}, 'shape', {}, 'mean', {}, ...
'mode', {}, 'stdev', {}, 'date1', {}, ...
'date2', {}, 'shift', {}, 'variance', {});
estimation_info.measurement_error.prior = struct('name', {}, 'shape', {}, 'mean', {}, ...
'mode', {}, 'stdev', {}, 'date1', {}, ...
'date2', {}, 'shift', {}, 'variance', {});
estimation_info.measurement_error_corr.prior = struct('name', {}, 'shape', {}, 'mean', {}, ...
'mode', {}, 'stdev', {}, 'date1', {}, ...
'date2', {}, 'shift', {}, 'variance', {});
estimation_info.measurement_error.prior_index = {};
estimation_info.structural_innovation.prior_index = {};
estimation_info.measurement_error_corr.prior_index = {};
estimation_info.structural_innovation_corr.prior_index = {};
estimation_info.measurement_error.options_index = {};
estimation_info.structural_innovation.options_index = {};
estimation_info.measurement_error_corr.options_index = {};
estimation_info.structural_innovation_corr.options_index = {};
options_.initial_period = dynDate(1);
options_.dataset.firstobs = options_.initial_period;
options_.dataset.lastobs = NaN;
options_.dataset.nobs = NaN;
options_.dataset.xls_sheet = NaN;
options_.dataset.xls_range = NaN;
options_.Harvey_scale_factor = 10;
options_.MaxNumberOfBytes = 1e6;
options_.MaximumNumberOfMegaBytes = 111;
@ -260,6 +332,7 @@ oo_.exo_det_steady_state = [];
oo_.exo_det_simul = [];
M_.params = [];
M_.endo_histval = [];
% BVAR
M_.bvar = [];

View File

@ -1,512 +0,0 @@
function VECTOR = LPTAU(I, N)
%
% I.M. SOBOL', V.I. TURCHANINOV, Yu.L. LEVITAN, B.V. SHUKHMAN
% KELDYSH INSTITUTE OF APPLIED MATHEMATICS
% RUSSIAN ACADEMY OF SCIENCES
%
% QUASIRANDOM SEQUENCE GENERATORS
% -------------------------------
%
% 28.11.1991
%
% NOTE TO THE USER BY the NEA Data Bank:
% This quasi random number generator has been made available to
% you on condition that its identity is preserved when used
% in computer programs. If its use leads to scientific publication
% of results you should cite it in the references, in addition
% no commercial use should be made unless agreed upon with the
% main author (Prof. I.M. Sobol')
%
% ABSTRACT
% ........
%
% POINTS BELONGING TO LP-TAU SEQUENCES UNIFORMLY DISTRIBUTED IN THE
% N-DIMENSIONAL UNIT CUBE ARE OFTEN USED IN NUMERICAL MATHEMATICS:
%
% - AS NODES FOR MULTIDIMENSIONAL INTEGRATION;
% - AS SEARCHING POINTS IN GLOBAL OPTIMIZATION;
% - AS TRIAL POINTS IN MULTI-CRITERIA DECISION MAKING;
% - AS QUASIRANDOM POINTS FOR QUASI-MONTECARLO ALGORITHMS;
% - ETC.
%
% THIS SUBROUTINE CONTAINS THE ALGORITHM FOR FAST GENERATION OF
% LP-TAU SEQUENCES THAT ARE SUITABLE FOR MULTI-PROCESSOR COMPUTATIONS.
% THE DIMENSIONS N.LE.51, THE NUMBER OF POINTS N.LT.2**30.
% THE PROGRAMMING LANGUAGE IS FORTRAN-77. THIS SUBROUTINE IS AVAILABLE
% ALSO IN %-LANGUAGE.
% THE REPORT DESCRIBING THE ALGORITHM CONTAINS THE DESCRIPTION OF THE
% ALGORITHM AND CERTAIN IMPORTANT PROPERTIES OF LP-TAU SEQUENCES AND
% THEIR GENERALIZATIONS ARE DISCUSSED.
%
% REFERENCE:
% I.M. SOBOL', V.I. TURCHANINOV, Yu.L. LEVITAN, B.V. SHUKHMAN
% KELDYSH INSTITUTE OF APPLIED MATHEMATICS
% RUSSIAN ACADEMY OF SCIENCES
%
% QUASIRANDOM SEQUENCE GENERATORS
% MOSCOW 1992, IPM ZAK. NO.30 (100 COPIES)
%
% ------------------------------------------------------------------------
%
% INPUT PARAMETERS:
%
% I - NUMBER OF THE POINT (I=(0,2**30-1)),
% N - DIMENSION OF THE POINT (0<N<52);
%
% OUTPUT PARAMETER:
%
% VECTOR(N) - N-VECTOR CONTAINING THE CARTESIAN CO-ORDINATES OF
% THE I-TH POINT.
%
%
% TO CALL THE SUBROUTINE WRITE:
%
% CALL LPTAU(I,N,VECTOR)
% WHERE I, N: INTEGER CAPABLE OF STORING 2**30 (INTEGER*4 ON IBM
% OR OTHER 32 BIT/WORD MACHINES)
% VECTOR: DOUBLE PRECISION ARRAY WHOSE LENGTH < 52.
%
% INTEGER QP
% QP = QUANTITY POWER
%
% PARAMETER (MAXDIM=51, QP=30, MAXNUM=2**30-1)
MAXDIM=51; QP=30; MAXNUM=2^30-1;
%
% THE DIMENSION OF THE POINT CANNOT EXCEED MAXDIM
% THE TOTAL NUMBER OF GENERATED POINTS CANNOT EXCEED 2**QP
% MAXNUM=2**30-1 // 1073741823
%
% DOUBLE PRECISION VECTOR(N)
% INTEGER I,N
%
% INTEGER PRVNUM,PRVDIM
% INTEGER DIRECT(MAXDIM,QP), MASKV(MAXDIM)
% DOUBLE PRECISION SCALE
%
% Translated into MATLAB by M. Ratto
VECTOR=zeros(1,N);
SCALE =9.31322574615478516E-10;
persistent PRVNUM PRVDIM MASKV DIRECT
if isempty(PRVNUM), PRVNUM=-2; end,
if isempty(PRVDIM), PRVDIM=0; end,
if isempty(DIRECT), ...
DIRECT(1:MAXDIM,1)=[
536870912, 536870912, 536870912, 536870912, 536870912, ...
536870912, 536870912, 536870912, 536870912, 536870912, ...
536870912, 536870912, 536870912, 536870912, 536870912, ...
536870912, 536870912, 536870912, 536870912, 536870912, ...
536870912, 536870912, 536870912, 536870912, 536870912, ...
536870912, 536870912, 536870912, 536870912, 536870912, ...
536870912, 536870912, 536870912, 536870912, 536870912, ...
536870912, 536870912, 536870912, 536870912, 536870912, ...
536870912, 536870912, 536870912, 536870912, 536870912, ...
536870912, 536870912, 536870912, 536870912, 536870912, ...
536870912]';
DIRECT(1:MAXDIM,2)=[
805306368, 268435456, 805306368, 268435456, 805306368, ...
268435456, 805306368, 268435456, 268435456, 805306368, ...
268435456, 805306368, 268435456, 805306368, 268435456, ...
805306368, 805306368, 268435456, 805306368, 268435456, ...
805306368, 268435456, 805306368, 268435456, 268435456, ...
805306368, 268435456, 805306368, 268435456, 805306368, ...
268435456, 805306368, 805306368, 268435456, 805306368, ...
268435456, 805306368, 268435456, 805306368, 268435456, ...
268435456, 805306368, 268435456, 805306368, 268435456, ...
805306368, 268435456, 805306368, 805306368, 268435456, ...
805306368]';
DIRECT(1:MAXDIM,3)=[
939524096, 939524096, 134217728, 671088640, 402653184, ...
402653184, 671088640, 134217728, 671088640, 402653184, ...
939524096, 134217728, 671088640, 939524096, 134217728, ...
671088640, 134217728, 939524096, 939524096, 402653184, ...
402653184, 134217728, 671088640, 402653184, 671088640, ...
402653184, 402653184, 671088640, 134217728, 134217728, ...
939524096, 939524096, 939524096, 939524096, 134217728, ...
671088640, 402653184, 402653184, 671088640, 134217728, ...
671088640, 402653184, 939524096, 134217728, 671088640, ...
939524096, 134217728, 671088640, 134217728, 939524096, ...
939524096]';
DIRECT(1:MAXDIM,4)=[
1006632960, 67108864, 603979776, 1006632960, 335544320, ...
469762048, 872415232, 738197504, 469762048, 872415232, ...
1006632960, 67108864, 872415232, 469762048, 469762048, ...
469762048, 1006632960, 335544320, 872415232, 603979776, ...
201326592, 67108864, 335544320, 67108864, 201326592, ...
738197504, 1006632960, 738197504, 603979776, 335544320, ...
738197504, 67108864, 1006632960, 67108864, 603979776, ...
1006632960, 335544320, 469762048, 872415232, 738197504, ...
469762048, 872415232, 1006632960, 67108864, 872415232, ...
469762048, 469762048, 469762048, 1006632960, 335544320, ...
872415232]';
DIRECT(1:MAXDIM,5)=[
1040187392, 637534208, 1040187392, 838860800, 167772160, ...
234881024, 167772160, 1040187392, 436207616, 33554432, ...
570425344, 1040187392, 503316480, 838860800, 973078528, ...
167772160, 234881024, 436207616, 771751936, 100663296, ...
234881024, 905969664, 771751936, 33554432, 838860800, ...
973078528, 905969664, 33554432, 301989888, 838860800, ...
100663296, 234881024, 1040187392, 637534208, 1040187392, ...
838860800, 167772160, 234881024, 167772160, 1040187392, ...
436207616, 33554432, 570425344, 1040187392, 503316480, ...
838860800, 973078528, 167772160, 234881024, 436207616, ...
771751936]';
DIRECT(1:MAXDIM,6)=[
1056964608, 352321536, 50331648, 419430400, 956301312, ...
889192448, 620756992, 117440512, 956301312, 922746880, ...
822083584, 218103808, 587202560, 385875968, 452984832, ...
218103808, 184549376, 285212672, 150994944, 956301312, ...
352321536, 654311424, 553648128, 352321536, 788529152, ...
956301312, 587202560, 251658240, 218103808, 721420288, ...
251658240, 1056964608, 520093696, 889192448, 587202560, ...
956301312, 419430400, 352321536, 83886080, 654311424, ...
419430400, 385875968, 285212672, 754974720, 50331648, ...
922746880, 989855744, 754974720, 721420288, 822083584, ...
687865856]';
DIRECT(1:MAXDIM,7)=[
1065353216, 1065353216, 578813952, 25165824, 125829120, ...
964689920, 327155712, 310378496, 360710144, 360710144, ...
159383552, 444596224, 947912704, 662700032, 444596224, ...
528482304, 578813952, 578813952, 75497472, 1065353216, ...
92274688, 511705088, 897581056, 847249408, 662700032, ...
931135488, 411041792, 713031680, 696254464, 528482304, ...
209715200, 578813952, 1065353216, 1065353216, 578813952, ...
25165824, 125829120, 964689920, 327155712, 310378496, ...
360710144, 360710144, 159383552, 444596224, 947912704, ...
662700032, 444596224, 528482304, 578813952, 578813952, ...
75497472]';
DIRECT(1:MAXDIM,8)=[
1069547520, 4194304, 826277888, 624951296, 616562688, ...
801112064, 952107008, 406847488, 398458880, 331350016, ...
356515840, 46137344, 1010827264, 1069547520, 734003200, ...
113246208, 490733568, 633339904, 910163968, 801112064, ...
893386752, 851443712, 801112064, 918552576, 742391808, ...
499122176, 62914560, 264241152, 708837376, 717225984, ...
759169024, 499122176, 272629760, 423624704, 155189248, ...
239075328, 205520896, 943718400, 373293056, 457179136, ...
406847488, 71303168, 473956352, 54525952, 624951296, ...
104857600, 1019215872, 901775360, 213909504, 281018368, ...
851443712]';
DIRECT(1:MAXDIM,9)=[
1071644672, 543162368, 190840832, 329252864, 853540864, ...
132120576, 778043392, 73400320, 178257920, 522190848, ...
639631360, 534773760, 991952896, 333447168, 48234496, ...
1059061760, 761266176, 673185792, 220200960, 396361728, ...
362807296, 815792128, 819986432, 346030080, 39845888, ...
752877568, 387973120, 643825664, 291504128, 274726912, ...
568328192, 526385152, 673185792, 98566144, 396361728, ...
727711744, 1042284544, 106954752, 299892736, 912261120, ...
44040192, 895483904, 333447168, 551550976, 467664896, ...
618659840, 606076928, 274726912, 245366784, 446693376, ...
421527552]';
DIRECT(1:MAXDIM,10)=[
1072693248, 273678336, 644874240, 753926144, 495976448, ...
869269504, 355467264, 57671680, 816840704, 961544192, ...
804257792, 495976448, 347078656, 426770432, 1066401792, ...
372244480, 84934656, 208666624, 313524224, 598736896, ...
487587840, 965738496, 1011875840, 296747008, 393216000, ...
523239424, 720371712, 823132160, 128974848, 407896064, ...
747634688, 850395136, 873463808, 504365056, 481296384, ...
686817280, 592445440, 995098624, 498073600, 969932800, ...
586153984, 1039138816, 814743552, 523239424, 294649856, ...
305135616, 506462208, 11534336, 449839104, 619708416, ...
479199232]';
DIRECT(1:MAXDIM,11)=[
1073217536, 947388416, 1070071808, 977797120, 365428736, ...
702021632, 461897728, 829947904, 425197568, 634912768, ...
437780480, 582483968, 792199168, 315097088, 611844096, ...
667418624, 166199296, 513277952, 187170816, 1036517376, ...
25690112, 201850880, 443023360, 990380032, 63438848, ...
211288064, 983040000, 1069023232, 421003264, 742916096, ...
487063552, 363331584, 973602816, 286785536, 171442176, ...
669515776, 110624768, 383254528, 289931264, 352845824, ...
878182400, 655884288, 836239360, 765984768, 549978112, ...
655884288, 85458944, 591921152, 563609600, 277348352, ...
919076864]';
DIRECT(1:MAXDIM,12)=[
1073479680, 71565312, 2359296, 891551744, 158597120, ...
383516672, 1019478016, 947126272, 621019136, 714866688, ...
738459648, 265027584, 468975616, 131858432, 504627200, ...
581173248, 266600448, 865861632, 658243584, 546045952, ...
521404416, 304873472, 1060896768, 163840000, 305922048, ...
257163264, 50069504, 773062656, 59506688, 779354112, ...
165937152, 587988992, 486801408, 160694272, 90439680, ...
423362560, 536608768, 614203392, 56885248, 999030784, ...
10747904, 764674048, 25952256, 989069312, 352583680, ...
799801344, 261357568, 873201664, 40108032, 769392640, ...
254541824]';
DIRECT(1:MAXDIM,13)=[
1073610752, 644218880, 538836992, 455475200, 1062600704, ...
139329536, 205651968, 905052160, 797048832, 452329472, ...
973471744, 627703808, 614072320, 803078144, 637403136, ...
835059712, 949878784, 662044672, 767950848, 426901504, ...
448659456, 23986176, 1016201216, 524943360, 525991936, ...
618790912, 781058048, 761659392, 458096640, 226361344, ...
950665216, 952500224, 516030464, 337510400, 496107520, ...
830865408, 944111616, 636354560, 978452480, 921567232, ...
533594112, 7471104, 678035456, 471203840, 1065746432, ...
575275008, 996540416, 909246464, 879362048, 637927424, ...
25821184]';
DIRECT(1:MAXDIM,14)=[
1073676288, 357892096, 808648704, 2424832, 2555904, ...
624230400, 69271552, 456851456, 1052966912, 600637440, ...
487260160, 794624000, 386727936, 467599360, 798031872, ...
630652928, 340983808, 493944832, 37945344, 264175616, ...
263520256, 833421312, 235077632, 464846848, 534839296, ...
992411648, 10813440, 367067136, 116457472, 115015680, ...
928710656, 619773952, 813760512, 1043398656, 967770112, ...
912850944, 72155136, 1009057792, 668532736, 462356480, ...
267321344, 795803648, 635764736, 574160896, 1003421696, ...
181075968, 56688640, 388562944, 190906368, 657915904, ...
474939392]';
DIRECT(1:MAXDIM,15)=[
1073709056, 1073709056, 137003008, 547782656, 545095680, ...
26836992, 34701312, 354385920, 925663232, 656965632, ...
327581696, 894795776, 110067712, 1038057472, 209354752, ...
596541440, 42631168, 471433216, 52527104, 666861568, ...
706707456, 674070528, 824410112, 305496064, 136282112, ...
847740928, 531464192, 222920704, 379289600, 507740160, ...
11894784, 1053392896, 129990656, 557547520, 666468352, ...
1061912576, 576684032, 1041334272, 380469248, 114196480, ...
133070848, 517046272, 129990656, 790396928, 563773440, ...
388333568, 661749760, 446791680, 737378304, 229998592, ...
348225536]';
DIRECT(1:MAXDIM,16)=[
1073725440, 16384, 605372416, 275234816, 817971200, ...
603963392, 555335680, 721534976, 997801984, 1028767744, ...
407060480, 375275520, 256688128, 1021165568, 303349760, ...
1022476288, 234143744, 106708992, 732971008, 733954048, ...
789889024, 879575040, 764657664, 762658816, 1010843648, ...
941080576, 827932672, 98942976, 1051738112, 624934912, ...
993280000, 134070272, 201375744, 567558144, 882163712, ...
649084928, 356564992, 489439232, 637091840, 60637184, ...
199278592, 815677440, 927678464, 94519296, 419184640, ...
933838848, 426655744, 911130624, 171393024, 561332224, ...
471613440]';
DIRECT(1:MAXDIM,17)=[
1073733632, 536895488, 1043685376, 679944192, 417505280, ...
301981696, 832561152, 210542592, 167501824, 1071341568, ...
229302272, 970661888, 732176384, 576659456, 402464768, ...
451584000, 368467968, 928260096, 933847040, 29319168, ...
582934528, 772612096, 330014720, 647323648, 174071808, ...
1008689152, 295919616, 353869824, 177774592, 580198400, ...
381837312, 638574592, 637558784, 679370752, 504012800, ...
747118592, 429973504, 1032609792, 932667392, 583360512, ...
969498624, 1056333824, 660955136, 247488512, 153509888, ...
242180096, 205840384, 797499392, 824565760, 234348544, ...
842326016]';
DIRECT(1:MAXDIM,18)=[
1073737728, 268455936, 52785152, 1020628992, 345018368, ...
452972544, 704442368, 255987712, 750759936, 697692160, ...
196677632, 764604416, 485625856, 522022912, 680620032, ...
362270720, 838103040, 83972096, 629133312, 46108672, ...
867561472, 725422080, 184504320, 751112192, 191918080, ...
306425856, 507310080, 30453760, 281858048, 604000256, ...
208662528, 319557632, 318779392, 476139520, 863719424, ...
567062528, 521179136, 712790016, 610299904, 293687296, ...
1023086592, 549089280, 1065242624, 707751936, 363024384, ...
16674816, 197136384, 1037561856, 195112960, 372707328, ...
992751616]';
DIRECT(1:MAXDIM,19)=[
1073739776, 939554816, 580732928, 854333440, 172619776, ...
511694848, 936142848, 518199296, 593348608, 225527808, ...
900982784, 180279296, 168904704, 62814208, 754485248, ...
730691584, 1005996032, 411174912, 249866240, 641669120, ...
1008719872, 749066240, 860993536, 94177280, 432564224, ...
226355200, 925784064, 995657728, 967731200, 436226048, ...
913799168, 549894144, 964696064, 843315200, 445863936, ...
1047422976, 548947968, 492066816, 953870336, 1002653696, ...
861440000, 385636352, 325253120, 187353088, 653584384, ...
1008269312, 748693504, 1013016576, 55814144, 255170560, ...
260708352]';
DIRECT(1:MAXDIM,20)=[
1073740800, 67126272, 829514752, 423777280, 968297472, ...
205511680, 147076096, 926669824, 202300416, 118395904, ...
381332480, 1002738688, 743042048, 292551680, 584567808, ...
284339200, 183936000, 616762368, 435221504, 159376384, ...
907322368, 595696640, 247497728, 553735168, 826051584, ...
564454400, 446024704, 214236160, 33661952, 251685888, ...
660327424, 284244992, 859868160, 722502656, 622844928, ...
324342784, 682374144, 400579584, 405353472, 605187072, ...
840682496, 212956160, 157891584, 193201152, 437990400, ...
573578240, 368053248, 580197376, 937905152, 565527552, ...
89064448]';
DIRECT(1:MAXDIM,21)=[
1073741312, 637560320, 189496832, 27474432, 129338880, ...
908054016, 641870336, 186375680, 302677504, 763663872, ...
103878144, 325187072, 858254848, 922041856, 261924352, ...
954978816, 292822528, 849512960, 210311680, 933232128, ...
691981824, 155417088, 627070464, 416795136, 182081024, ...
513433088, 848658944, 515770880, 627273216, 629169664, ...
414566912, 147450368, 698353152, 244844032, 226578944, ...
1020087808, 886978048, 389697024, 1007004160, 839646720, ...
621924864, 549962240, 609583616, 735976960, 87342592, ...
1058542080, 163066368, 307997184, 876471808, 794280448, ...
675386880]';
DIRECT(1:MAXDIM,22)=[
1073741568, 352343296, 644236032, 636735232, 615860480, ...
959444224, 287380736, 1007410432, 890187008, 399480576, ...
520092928, 643311360, 816901376, 695310080, 1019229440, ...
77034240, 733295872, 1035127552, 986582784, 332381952, ...
334852352, 364956416, 596672256, 800381696, 480316672, ...
574863104, 647347968, 702910208, 499965184, 364968704, ...
120862976, 1023256320, 995114240, 13951232, 32520448, ...
702127360, 45176064, 444945664, 237860096, 152839936, ...
530633984, 429135616, 267272448, 884808960, 933712640, ...
61605632, 174335744, 564911360, 302327552, 650589440, ...
450649344]';
DIRECT(1:MAXDIM,23)=[
1073741696, 1065385856, 1073734272, 331949184, 842310784, ...
799537536, 965852032, 369351808, 662886016, 86119808, ...
865109888, 299633792, 422735488, 181087360, 174252416, ...
1041212544, 840196224, 750314368, 391053440, 903306880, ...
742365312, 236995200, 42492800, 946000512, 771692416, ...
897405824, 613803136, 924258688, 808338304, 1038125440, ...
683814272, 177186176, 766008960, 704549248, 194555008, ...
306383744, 496592512, 416020864, 655186816, 1032204928, ...
694773632, 577910144, 45797760, 910332544, 536014976, ...
675946368, 987635840, 788223872, 353993856, 96313472, ...
85248640]';
DIRECT(1:MAXDIM,24)=[
1073741760, 4210752, 12608, 744788544, 494377792, ...
115601344, 769248448, 990895808, 851706304, 979326784, ...
692061120, 429015104, 217132864, 736067008, 55694400, ...
456152640, 631601984, 787264192, 898599104, 478383808, ...
507774272, 458270272, 392995136, 872482496, 124824768, ...
1034161344, 362141632, 1053833280, 943810496, 428920128, ...
795835200, 835462848, 961843520, 606198080, 785652672, ...
154954432, 491617344, 297351232, 580735552, 634587968, ...
185277760, 141505600, 417673280, 106907456, 395575616, ...
566958656, 352651200, 415242688, 26635200, 1030317504, ...
247212992]';
DIRECT(1:MAXDIM,25)=[
1073741792, 543187040, 536882016, 981766752, 357858144, ...
810665952, 369083488, 613015520, 86115232, 845149216, ...
606076960, 1018241504, 35245536, 635403744, 236633696, ...
407451232, 427671904, 460141792, 116344544, 990246688, ...
1024892576, 883429408, 150655392, 173476896, 197666464, ...
1016740448, 605376608, 970487008, 1006881248, 598265632, ...
1022861728, 489055392, 216680608, 371439072, 53140576, ...
965114400, 98534112, 209692256, 264598496, 321437280, ...
545303840, 324119904, 876587488, 778239712, 802033120, ...
44406496, 665829024, 803213920, 309742112, 735181344, ...
1036125600]';
DIRECT(1:MAXDIM,26)=[
1073741808, 273698896, 805312112, 903123536, 153504624, ...
689632208, 432848944, 859888752, 510853776, 240805744, ...
250610192, 852489168, 139460144, 283082224, 222702928, ...
342181488, 922872432, 187528656, 360637424, 814514736, ...
301037840, 800872624, 272595184, 158627600, 238839088, ...
927181136, 710248688, 788854608, 1048376560, 484709072, ...
85709008, 1065469744, 429237328, 490778384, 848024144, ...
443114288, 687907504, 474573648, 45706416, 681465296, ...
805303216, 14567920, 369839504, 440402480, 797652624, ...
115959088, 929784560, 91226544, 205943056, 288358704, ...
950217296]';
DIRECT(1:MAXDIM,27)=[
1073741816, 947419256, 134232008, 460100200, 1073685336, ...
398354904, 1069834248, 686054696, 108299224, 273898840, ...
731382552, 938170664, 86812552, 677346808, 602208712, ...
652635752, 209797064, 323968376, 384078968, 561178056, ...
923399256, 996597176, 942777416, 885167400, 89791048, ...
956122504, 87393464, 987661336, 993412952, 827749496, ...
903211272, 33649816, 594875176, 65052056, 822835240, ...
625704888, 1065374584, 612232920, 536299720, 1046858504, ...
939518840, 160843032, 654656088, 744496936, 872197704, ...
219111576, 829112632, 455614152, 1064192952, 313010856, ...
820754920]';
DIRECT(1:MAXDIM,28)=[
1073741820, 71582788, 603989028, 37500, 120660, ...
182195932, 213921796, 473570692, 41763004, 156352828, ...
88343188, 698012780, 666108868, 337608156, 1054329884, ...
799472220, 641885244, 392104980, 502284340, 1002405628, ...
249907140, 75586228, 206587660, 565275348, 1021426548, ...
425987996, 598058580, 401940420, 919427316, 822049324, ...
115655868, 759299588, 562394644, 990942164, 1003212268, ...
598750292, 675334852, 675293340, 445665316, 903023756, ...
872412692, 172640916, 1051615436, 91229620, 94586692, ...
344873452, 7029156, 683421900, 434258804, 955002092, ...
436424380]';
DIRECT(1:MAXDIM,29)=[
1073741822, 644245094, 1040195102, 537039474, 537089354, ...
642656334, 73402402, 664239174, 1014620426, 500917234, ...
1042677826, 347788318, 1069154738, 373259762, 362587674, ...
239395914, 132283950, 903121466, 445210638, 968851198, ...
230153254, 935549622, 308815630, 861891286, 496995094, ...
670740382, 657311830, 573565382, 248331478, 1064650798, ...
338312798, 669059078, 1065190902, 379125622, 111897166, ...
520650422, 740300390, 1046483950, 66254898, 992513134, ...
234871606, 610553330, 450553866, 566758134, 787800806, ...
1071709866, 971732606, 528102354, 790919122, 917384366, ...
656244162]';
DIRECT(1:MAXDIM,30)=[
1073741823, 357913941, 50344755, 268538457, 805540473, ...
3487029, 3146769, 592038967, 729591963, 781578389, ...
66781679, 113956361, 331153483, 967802327, 935343371, ...
324351309, 609305915, 818137857, 131059769, 918519549, ...
827341719, 922770101, 456972047, 363883221, 1057014661, ...
900956063, 580478293, 520383687, 315470533, 227601711, ...
169598765, 909227271, 796983815, 349496709, 393974911, ...
378320037, 777004343, 927999269, 616377385, 345930959, ...
989849787, 627400495, 892675125, 260314121, 1041964063, ...
531367163, 195776757, 237309785, 949187605, 719002587, ...
495745187]';
end
%
% TRAP FOR A WRONG SUBROUTINE CALL
%
if ((I<0) | (N<1) | (I>MAXNUM) | (N>MAXDIM)),
disp('LP-TAU CALL FAILED')
disp(' PRESS <ENTER> TO EXIT LPTAU')
pause
return
end
if ((PRVNUM+1==I) & (N<=PRVDIM)),
%
% RECURRENT GENERATION OF THE POINT
%
%
% SEARCH POSITION OF THE RIGHTMOST "1"
% IN THE BINARY REPRESENTATION OF I
%
L=0;
POS=0;
while L<QP & POS==0,
L=L+1;
POS=bitand(bitshift(I,-(L-1)),1);
end
%
% RIGHTMOST POSITION IS L
%
for J=1:N
MASKV(J)=bitxor(MASKV(J),DIRECT(J,L));
VECTOR(J)=MASKV(J)*SCALE;
end
else
%
% GENERATION OF THE POINT FROM "I" AND "N"
%
MASKV=zeros(1,N);
IM=bitxor(I,bitshift(I,-1));
M=0;
while IM~=0 & M<QP,
M=M+1;
if (bitand(IM,1)==1),
for J=1:N
MASKV(J)=bitxor(MASKV(J),DIRECT(J,M));
end
end
IM=bitshift(IM,-1);
end
for J=1:N
VECTOR(J)=MASKV(J)*SCALE;
end
end
%
PRVNUM=I;
PRVDIM=N;
return

View File

@ -1,4 +1,4 @@
function s=skewness(y),
function s=gsa_skewness(y),
% y=stand_(y);
% s=mean(y.^3);

View File

@ -1,5 +1,5 @@
function [tadj, iff] = speed(A,B,mf,p),
% [tadj, iff] = speed(A,B,mf,p),
function [tadj, iff] = gsa_speed(A,B,mf,p),
% [tadj, iff] = gsa_speed(A,B,mf,p),
%
% Part of the Sensitivity Analysis Toolbox for DYNARE
%

View File

@ -3,10 +3,10 @@ function [yy, xdir, isig, lam]=log_trans_(y0,xdir0)
if nargin==1,
xdir0='';
end
f=inline('skewness(log(y+lam))','lam','y');
f=inline('gsa_skewness(log(y+lam))','lam','y');
isig=1;
if ~(max(y0)<0 | min(y0)>0)
if skewness(y0)<0,
if gsa_skewness(y0)<0,
isig=-1;
y0=-y0;
end

View File

@ -1,25 +0,0 @@
function [lpmat] = lptauSEQ(Nsam,Nvar)
% [lpmat] = lptauSEQ(Nsam,Nvar)
%
% function call LPTAU and generates a sample of dimension Nsam for a
% number of parameters Nvar
%
% Copyright (C) 2005 Marco Ratto
% THIS PROGRAM WAS WRITTEN FOR MATLAB BY
% Marco Ratto,
% Unit of Econometrics and Statistics AF
% (http://www.jrc.cec.eu.int/uasa/),
% IPSC, Joint Research Centre
% The European Commission,
% TP 361, 21020 ISPRA(VA), ITALY
% marco.ratto@jrc.it
%
clear lptau
lpmat = zeros(Nsam, Nvar);
for j=1:Nsam,
lpmat(j,:)=LPTAU(j,Nvar);
end
return

View File

@ -127,7 +127,7 @@ if opt_gsa.load_ident_files==0,
% ino=find(~ismember([1:nr],io));
% T2=A(ino,1:nr,:);
R=A(:,nr+1:nc,:);
% [tadj, iff] = speed(A(1:nr,1:nr,:),R,io,0.5);
% [tadj, iff] = gsa_speed(A(1:nr,1:nr,:),R,io,0.5);
% [tadj, j0, ir_tadj, ic_tadj] = teff(tadj,Nsam,istable);
% [iff, j0, ir_if, ic_if] = teff(iff,Nsam,istable);

View File

@ -28,8 +28,7 @@ function x0 = stab_map_(OutputDirectoryName)
% 3) Bivariate plots of significant correlation patterns
% ( abs(corrcoef) > alpha2) under the stable and unacceptable subsets
%
% USES lptauSEQ,
% stab_map_1, stab_map_2
% USES qmc_sequence, stab_map_1, stab_map_2
%
% Part of the Sensitivity Analysis Toolbox for DYNARE
%
@ -118,7 +117,7 @@ if fload==0,
Nsam=size(lpmat,1);
else
if np<52 & ilptau>0,
[lpmat] = lptauSEQ(Nsam,np); % lptau
[lpmat] = qmc_sequence(np, int64(0), 0, Nsam)';
if np>30 | ilptau==2, % scrambled lptau
for j=1:np,
lpmat(:,j)=lpmat(randperm(Nsam),j);

View File

@ -135,7 +135,7 @@ if info(1)==0,
data_info.data=oo_.endo_simul(options_.varobs_id,100+1:end);
% datax=data;
derivatives_info.no_DLIK=1;
[fval,cost_flag,ys,trend_coeff,info,M_,options_,bayestopt_,oo_,DLIK,AHess] = DsgeLikelihood(params',data_info,options_,M_,estim_params_,bayestopt_,oo_,derivatives_info);
[fval,cost_flag,ys,trend_coeff,info,M_,options_,bayestopt_,oo_,DLIK,AHess] = dsge_likelihood(params',data_info,options_,M_,estim_params_,bayestopt_,oo_,derivatives_info);
% fval = DsgeLikelihood(xparam1,data_info,options_,M_,estim_params_,bayestopt_,oo_);
AHess=-AHess;
ide_hess.AHess= AHess;

View File

@ -41,7 +41,7 @@ end
if DynareOptions.dsge_var
[fval,cost_flag,info] = DsgeVarLikelihood(xparam1,DynareDataset,DynareOptions,Model,EstimatedParameters,BayesInfo,DynareResults);
else
[fval,cost_flag,ys,trend_coeff,info] = DsgeLikelihood(xparam1,DynareDataset,DynareOptions,Model,EstimatedParameters,BayesInfo,DynareResults);
[fval,cost_flag,ys,trend_coeff,info] = dsge_likelihood(xparam1,DynareDataset,DynareOptions,Model,EstimatedParameters,BayesInfo,DynareResults);
end
if info(1) > 0

View File

@ -79,10 +79,10 @@ more off
if strcmpi(flag,'--test')
if nargin>1
dynare_config([],0);
dynare_path = dynare_config([],0);
number_of_matlab_routines = length(varargin);
for i=1:number_of_matlab_routines
dynTest(varargin{i});
dynTest(varargin{i},dynare_path);
end
else
disp('You have to specify at least one matlab routine after --test flag!')

View File

@ -2,7 +2,7 @@ function [LIK, likk, a, P] = kalman_filter(Y,start,last,a,P,kalman_tol,riccati_t
% Computes the likelihood of a stationnary state space model.
%@info:
%! @deftypefn {Function File} {[@var{LIK},@var{likk},@var{a},@var{P} ] =} DsgeLikelihood (@var{Y}, @var{start}, @var{last}, @var{a}, @var{P}, @var{kalman_tol}, @var{riccati_tol},@var{presample},@var{T},@var{Q},@var{R},@var{H},@var{Z},@var{mm},@var{pp},@var{rr},@var{Zflag},@var{diffuse_periods})
%! @deftypefn {Function File} {[@var{LIK},@var{likk},@var{a},@var{P} ] =} kalman_filter (@var{Y}, @var{start}, @var{last}, @var{a}, @var{P}, @var{kalman_tol}, @var{riccati_tol},@var{presample},@var{T},@var{Q},@var{R},@var{H},@var{Z},@var{mm},@var{pp},@var{rr},@var{Zflag},@var{diffuse_periods})
%! @anchor{kalman_filter}
%! @sp 1
%! Computes the likelihood of a stationary state space model, given initial condition for the states (mean and variance).
@ -63,7 +63,7 @@ function [LIK, likk, a, P] = kalman_filter(Y,start,last,a,P,kalman_tol,riccati_t
%! @sp 2
%! @strong{This function is called by:}
%! @sp 1
%! @ref{DsgeLikelihood}
%! @ref{dsge_likelihood}
%! @sp 2
%! @strong{This function calls:}
%! @sp 1

View File

@ -2,7 +2,7 @@ function [LIK, likk, a] = kalman_filter_ss(Y,start,last,a,T,K,iF,dF,Z,pp,Zflag)
% Computes the likelihood of a stationnary state space model (steady state kalman filter).
%@info:
%! @deftypefn {Function File} {[@var{LIK},@var{likk},@var{a},@var{P} ] =} DsgeLikelihood (@var{Y}, @var{start}, @var{last}, @var{a}, @var{P}, @var{kalman_tol}, @var{riccati_tol},@var{presample},@var{T},@var{Q},@var{R},@var{H},@var{Z},@var{mm},@var{pp},@var{rr},@var{Zflag},@var{diffuse_periods})
%! @deftypefn {Function File} {[@var{LIK},@var{likk},@var{a},@var{P} ] =} kalman_filter_ss (@var{Y}, @var{start}, @var{last}, @var{a}, @var{P}, @var{kalman_tol}, @var{riccati_tol},@var{presample},@var{T},@var{Q},@var{R},@var{H},@var{Z},@var{mm},@var{pp},@var{rr},@var{Zflag},@var{diffuse_periods})
%! @anchor{kalman_filter}
%! @sp 1
%! Computes the likelihood of a stationary state space model, given initial condition for the states (mean), the steady state kalman gain and the steady state inveverted covariance matrix of the prediction errors.

View File

@ -69,7 +69,7 @@ function [LIK, likk,a,P] = univariate_kalman_filter(data_index,number_of_observa
%! @sp 2
%! @strong{This function is called by:}
%! @sp 1
%! @ref{DsgeLikelihood}
%! @ref{dsge_likelihood}
%! @sp 2
%! @strong{This function calls:}
%! @sp 1

View File

@ -73,7 +73,7 @@ function [dLIK, dlikk, a, Pstar, llik] = univariate_kalman_filter_d(data_index,
%! @sp 2
%! @strong{This function is called by:}
%! @sp 1
%! @ref{DsgeLikelihood}, @ref{DsgeLikelihood_hh}
%! @ref{dsge_likelihood}, @ref{DsgeLikelihood_hh}
%! @sp 2
%! @strong{This function calls:}
%! @sp 1

View File

@ -35,29 +35,16 @@ if isempty(oo_.steady_state)
oo_.steady_state = zeros(M_.endo_nbr,1);
end
if isempty(oo_.endo_simul)
if isempty(M_.endo_histval)
if isempty(ys0_)
oo_.endo_simul = [oo_.steady_state*ones(1,M_.maximum_lag+options_.periods+M_.maximum_lead)];
else
oo_.endo_simul = [ys0_*ones(1,M_.maximum_lag) oo_.steady_state*ones(1,options_.periods+M_.maximum_lead)];
end
elseif size(oo_.endo_simul,2) < M_.maximum_lag+M_.maximum_lead+options_.periods
switch options_.deterministic_simulation_initialization
case 0
oo_.endo_simul = [oo_.endo_simul ...
oo_.steady_state*ones(1,M_.maximum_lag+options_.periods+M_.maximum_lead-size(oo_.endo_simul,2),1)];
case 1% A linear approximation is used to initialize the solution.
oldopt = options_;
options_.order = 1;
dr = oo_.dr;
dr.ys = oo_.steady_state;
[dr,info,M_,options_,oo_]=dr1(dr,0,M_,options_,oo_);
exogenous_variables = zeros(M_.maximum_lag+options_.periods+M_.maximum_lead-size(oo_.endo_simul,2)+1,0);
y0 = oo_.endo_simul(:,1:M_.maximum_lag);
oo_.endo_simul=simult_(y0,dr,exogenous_variables,1);
options_ = oldopt;
case 2% Homotopic mod: Leave endo_simul as it is.
otherwise
error('Unknown method.')
else
if ~isempty(ys0_)
error('histval and endval cannot be used simultaneously')
end
end
oo_.endo_simul = [M_.endo_histval ...
oo_.steady_state*ones(1,options_.periods+M_.maximum_lead)];
end

View File

@ -52,7 +52,7 @@ function ms_write_markov_file(fname, options)
%//=====================================================//
fprintf(fh,'//== Number of states for state_variable[%d] ==//\n', ...
i_chain);
n_states = length(options.ms.ms_chain(i_chain).state);
n_states = length(options.ms.ms_chain(i_chain).regime);
fprintf(fh,'%d\n\n',n_states);
%//== 03/15/06: DW TVBVAR code reads the data below and overwrite the prior data read somewhere else if any.
@ -64,7 +64,7 @@ function ms_write_markov_file(fname, options)
i_chain);
Alpha = ones(n_states,n_states);
for i_state = 1:n_states
p = 1-1/options.ms.ms_chain(i_chain).state(i_state).duration;
p = 1-1/options.ms.ms_chain(i_chain).regime(i_state).duration;
Alpha(i_state,i_state) = p*(n_states-1)/(1-p);
fprintf(fh,'%22.16f',Alpha(i_state,:));
fprintf(fh,'\n');

View File

@ -31,7 +31,7 @@ function plot_ms_probabilities(computed_probabilities, options_)
[T,num_grand_regimes] = size(computed_probabilities);
num_chains = length(options_.ms.ms_chain);
for i=1:num_chains
chains(i).num_regimes = length(options_.ms.ms_chain(i).state);
chains(i).num_regimes = length(options_.ms.ms_chain(i).regime);
chains(i).probabilities = zeros([T,chains(i).num_regimes]);
end

View File

@ -0,0 +1,375 @@
function [fval,exit_flag,ys,trend_coeff,info,Model,DynareOptions,BayesInfo,DynareResults] = non_linear_dsge_likelihood(xparam1,DynareDataset,DynareOptions,Model,EstimatedParameters,BayesInfo,DynareResults)
% Evaluates the posterior kernel of a dsge model using a non linear filter.
%@info:
%! @deftypefn {Function File} {[@var{fval},@var{exit_flag},@var{ys},@var{trend_coeff},@var{info},@var{Model},@var{DynareOptions},@var{BayesInfo},@var{DynareResults}] =} non_linear_dsge_likelihood (@var{xparam1},@var{DynareDataset},@var{DynareOptions},@var{Model},@var{EstimatedParameters},@var{BayesInfo},@var{DynareResults})
%! @anchor{dsge_likelihood}
%! @sp 1
%! Evaluates the posterior kernel of a dsge model using a non linear filter.
%! @sp 2
%! @strong{Inputs}
%! @sp 1
%! @table @ @var
%! @item xparam1
%! Vector of doubles, current values for the estimated parameters.
%! @item DynareDataset
%! Matlab's structure describing the dataset (initialized by dynare, see @ref{dataset_}).
%! @item DynareOptions
%! Matlab's structure describing the options (initialized by dynare, see @ref{options_}).
%! @item Model
%! Matlab's structure describing the Model (initialized by dynare, see @ref{M_}).
%! @item EstimatedParamemeters
%! Matlab's structure describing the estimated_parameters (initialized by dynare, see @ref{estim_params_}).
%! @item BayesInfo
%! Matlab's structure describing the priors (initialized by dynare, see @ref{bayesopt_}).
%! @item DynareResults
%! Matlab's structure gathering the results (initialized by dynare, see @ref{oo_}).
%! @end table
%! @sp 2
%! @strong{Outputs}
%! @sp 1
%! @table @ @var
%! @item fval
%! Double scalar, value of (minus) the likelihood.
%! @item exit_flag
%! Integer scalar, equal to zero if the routine return with a penalty (one otherwise).
%! @item ys
%! Vector of doubles, steady state level for the endogenous variables.
%! @item trend_coeffs
%! Matrix of doubles, coefficients of the deterministic trend in the measurement equation.
%! @item info
%! Integer scalar, error code.
%! @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 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.
%! @item info==41
%! At least one parameter is violating a lower bound condition.
%! @item info==42
%! At least one parameter is violating an upper bound condition.
%! @item info==43
%! The covariance matrix of the structural innovations is not positive definite.
%! @item info==44
%! The covariance matrix of the measurement errors is not positive definite.
%! @item info==45
%! Likelihood is not a number (NaN).
%! @item info==45
%! Likelihood is a complex valued number.
%! @end table
%! @item Model
%! Matlab's structure describing the model (initialized by dynare, see @ref{M_}).
%! @item DynareOptions
%! Matlab's structure describing the options (initialized by dynare, see @ref{options_}).
%! @item BayesInfo
%! Matlab's structure describing the priors (initialized by dynare, see @ref{bayesopt_}).
%! @item DynareResults
%! Matlab's structure gathering the results (initialized by dynare, see @ref{oo_}).
%! @end table
%! @sp 2
%! @strong{This function is called by:}
%! @sp 1
%! @ref{dynare_estimation_1}, @ref{mode_check}
%! @sp 2
%! @strong{This function calls:}
%! @sp 1
%! @ref{dynare_resolve}, @ref{lyapunov_symm}, @ref{priordens}
%! @end deftypefn
%@eod:
% Copyright (C) 2010-2011 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/>.
% AUTHOR(S) stephane DOT adjemian AT univ DASH lemans DOT fr
% frederic DOT karame AT univ DASH lemans DOT fr
% Declaration of the penalty as a persistent variable.
persistent penalty
persistent init_flag
persistent restrict_variables_idx observed_variables_idx state_variables_idx mf0 mf1
persistent sample_size number_of_state_variables number_of_observed_variables number_of_structural_innovations
% Initialization of the persistent variable.
if ~nargin || isempty(penalty)
penalty = 1e8;
if ~nargin, return, end
end
if nargin==1
penalty = xparam1;
return
end
% Initialization of the returned arguments.
fval = [];
ys = [];
trend_coeff = [];
cost_flag = 1;
% Set the number of observed variables
nvobs = DynareDataset.info.vobs;
%------------------------------------------------------------------------------
% 1. Get the structural parameters & define penalties
%------------------------------------------------------------------------------
% Return, with endogenous penalty, if some parameters are smaller than the lower bound of the prior domain.
if (DynareOptions.mode_compute~=1) & any(xparam1<BayesInfo.lb)
k = find(xparam1 < BayesInfo.lb);
fval = penalty+sum((BayesInfo.lb(k)-xparam1(k)).^2);
cost_flag = 0;
info = 41;
return
end
% Return, with endogenous penalty, if some parameters are greater than the upper bound of the prior domain.
if (DynareOptions.mode_compute~=1) & any(xparam1>BayesInfo.ub)
k = find(xparam1>BayesInfo.ub);
fval = penalty+sum((xparam1(k)-BayesInfo.ub(k)).^2);
cost_flag = 0;
info = 42;
return
end
% Get the diagonal elements of the covariance matrices for the structural innovations (Q) and the measurement error (H).
Q = Model.Sigma_e;
H = Model.H;
for i=1:EstimatedParameters_.nvx
k =EstimatedParameters_.var_exo(i,1);
Q(k,k) = xparam1(i)*xparam1(i);
end
offset = EstimatedParameters_.nvx;
if EstimatedParameters_.nvn
for i=1:EstimatedParameters_.nvn
k = EstimatedParameters_.var_endo(i,1);
H(k,k) = xparam1(i+offset)*xparam1(i+offset);
end
offset = offset+EstimatedParameters_.nvn;
else
H = zeros(nvobs);
end
% Get the off-diagonal elements of the covariance matrix for the structural innovations. Test if Q is positive definite.
if EstimatedParameters_.ncx
for i=1:EstimatedParameters_.ncx
k1 =EstimatedParameters_.corrx(i,1);
k2 =EstimatedParameters_.corrx(i,2);
Q(k1,k2) = xparam1(i+offset)*sqrt(Q(k1,k1)*Q(k2,k2));
Q(k2,k1) = Q(k1,k2);
end
% Try to compute the cholesky decomposition of Q (possible iff Q is positive definite)
[CholQ,testQ] = chol(Q);
if testQ
% The variance-covariance matrix of the structural innovations is not definite positive. We have to compute the eigenvalues of this matrix in order to build the endogenous penalty.
a = diag(eig(Q));
k = find(a < 0);
if k > 0
fval = penalty+sum(-a(k));
cost_flag = 0;
info = 43;
return
end
end
offset = offset+EstimatedParameters_.ncx;
end
% Get the off-diagonal elements of the covariance matrix for the measurement errors. Test if H is positive definite.
if EstimatedParameters_.ncn
for i=1:EstimatedParameters_.ncn
k1 = DynareOptions.lgyidx2varobs(EstimatedParameters_.corrn(i,1));
k2 = DynareOptions.lgyidx2varobs(EstimatedParameters_.corrn(i,2));
H(k1,k2) = xparam1(i+offset)*sqrt(H(k1,k1)*H(k2,k2));
H(k2,k1) = H(k1,k2);
end
% Try to compute the cholesky decomposition of H (possible iff H is positive definite)
[CholH,testH] = chol(H);
if testH
% The variance-covariance matrix of the measurement errors is not definite positive. We have to compute the eigenvalues of this matrix in order to build the endogenous penalty.
a = diag(eig(H));
k = find(a < 0);
if k > 0
fval = penalty+sum(-a(k));
cost_flag = 0;
info = 44;
return
end
end
offset = offset+EstimatedParameters_.ncn;
end
% Update estimated structural parameters in Mode.params.
if EstimatedParameters_.np > 0
Model.params(EstimatedParameters_.param_vals(:,1)) = xparam1(offset+1:end);
end
% Update Model.Sigma_e and Model.H.
Model.Sigma_e = Q;
Model.H = H;
%------------------------------------------------------------------------------
% 2. call model setup & reduction program
%------------------------------------------------------------------------------
% Linearize the model around the deterministic sdteadystate and extract the matrices of the state equation (T and R).
[T,R,SteadyState,info,Model,DynareOptions,DynareResults] = dynare_resolve(Model,DynareOptions,DynareResults,'restrict');
if info(1) == 1 || info(1) == 2 || info(1) == 5
fval = penalty+1;
cost_flag = 0;
return
elseif info(1) == 3 || info(1) == 4 || info(1)==6 ||info(1) == 19 || info(1) == 20 || info(1) == 21
fval = penalty+info(2);
cost_flag = 0;
return
end
% Define a vector of indices for the observed variables. Is this really usefull?...
BayesInfo.mf = BayesInfo.mf1;
% Define the deterministic linear trend of the measurement equation.
if DynareOptions.noconstant
constant = zeros(nvobs,1);
else
if DynareOptions.loglinear
constant = log(SteadyState(BayesInfo.mfys));
else
constant = SteadyState(BayesInfo.mfys);
end
end
% Define the deterministic linear trend of the measurement equation.
if BayesInfo.with_trend
trend_coeff = zeros(DynareDataset.info.nvobs,1);
t = DynareOptions.trend_coeffs;
for i=1:length(t)
if ~isempty(t{i})
trend_coeff(i) = evalin('base',t{i});
end
end
trend = repmat(constant,1,DynareDataset.info.ntobs)+trend_coeff*[1:DynareDataset.info.ntobs];
else
trend = repmat(constant,1,DynareDataset.info.ntobs);
end
% Get needed informations for kalman filter routines.
start = DynareOptions.presample+1;
np = size(T,1);
mf = BayesInfo.mf;
Y = transpose(dataset_.rawdata);
%------------------------------------------------------------------------------
% 3. Initial condition of the Kalman filter
%------------------------------------------------------------------------------
% Get decision rules and transition equations.
dr = DynareResults.dr;
% Set persistent variables (first call).
if isempty(init_flag)
mf0 = BayesInfo.mf0;
mf1 = BayesInfo.mf1;
restrict_variables_idx = BayesInfo.restrict_var_list;
observed_variables_idx = restrict_variables_idx(mf1);
state_variables_idx = restrict_variables_idx(mf0);
sample_size = size(Y,2);
number_of_state_variables = length(mf0);
number_of_observed_variables = length(mf1);
number_of_structural_innovations = length(Q);
init_flag = 1;
end
ReducedForm.ghx = dr.ghx(restrict_variables_idx,:);
ReducedForm.ghu = dr.ghu(restrict_variables_idx,:);
ReducedForm.ghxx = dr.ghxx(restrict_variables_idx,:);
ReducedForm.ghuu = dr.ghuu(restrict_variables_idx,:);
ReducedForm.ghxu = dr.ghxu(restrict_variables_idx,:);
ReducedForm.steadystate = dr.ys(dr.order_var(restrict_variables_idx));
ReducedForm.constant = ReducedForm.steadystate + .5*dr.ghs2(restrict_variables_idx);
ReducedForm.state_variables_steady_state = dr.ys(dr.order_var(state_variables_idx));
ReducedForm.Q = Q;
ReducedForm.H = H;
ReducedForm.mf0 = mf0;
ReducedForm.mf1 = mf1;
% Set initial condition.
switch DynareOptions.particle.initialization
case 1% Initial state vector variance is the ergodic variance associated to the first order Taylor-approximation of the model.
StateVectorMean = ReducedForm.constant(mf0);
StateVectorVariance = lyapunov_symm(ReducedForm.ghx(mf0,:),ReducedForm.ghu(mf0,:)*ReducedForm.Q*ReducedForm.ghu(mf0,:)',1e-12,1e-12);
case 2% Initial state vector variance is a monte-carlo based estimate of the ergodic variance (consistent with a k-order Taylor-approximation of the model).
StateVectorMean = ReducedForm.constant(mf0);
old_DynareOptionsperiods = DynareOptions.periods;
DynareOptions.periods = 5000;
y_ = simult(oo_.steady_state, dr);
y_ = y_(state_variables_idx,2001:5000);
StateVectorVariance = cov(y_');
DynareOptions.periods = old_DynareOptionsperiods;
clear('old_DynareOptionsperiods','y_');
case 3
StateVectorMean = ReducedForm.constant(mf0);
StateVectorVariance = DynareOptions.particle.initial_state_prior_std*eye(number_of_state_variables);
otherwise
error('Unknown initialization option!')
end
ReducedForm.StateVectorMean = StateVectorMean;
ReducedForm.StateVectorVariance = StateVectorVariance;
%------------------------------------------------------------------------------
% 4. Likelihood evaluation
%------------------------------------------------------------------------------
DynareOptions.warning_for_steadystate = 0;
LIK = feval(DynareOptions.particle.algorithm,ReducedForm,Y,[]);
if imag(LIK)
likelihood = penalty;
cost_flag = 0;
elseif isnan(LIK)
likelihood = penalty;
cost_flag = 0;
else
likelihood = LIK;
end
DynareOptions.warning_for_steadystate = 1;
% ------------------------------------------------------------------------------
% Adds prior if necessary
% ------------------------------------------------------------------------------
lnprior = priordens(xparam1,BayesInfo.pshape,BayesInfo.p6,BayesInfo.p7,BayesInfo.p3,BayesInfo.p4);
fval = (likelihood-lnprior);

View File

@ -73,7 +73,7 @@ if ~noprint
error('one (many) parameter(s) do(es) not satisfy the upper bound');
case 43
error('Covariance matrix of shocks is not positive definite')
case 44 %DsgeLikelihood_hh / DsgeLikelihood
case 44 %DsgeLikelihood_hh / dsge_likelihood
error('');
case 51
error('You are estimating a DSGE-VAR model, but the value of the dsge prior weight is too low!')

162
matlab/qmc_sequence.m Normal file
View File

@ -0,0 +1,162 @@
%@info:
%! @deftypefn {Mex File} {[@var{a}, @var{s}, @var{info}] =} qmc_sequence (@var{d}, @var{s}, @var{t}, @var{n}, @var{lu})
%! @anchor{qmc_sequence}
%! @sp 1
%! Computes quasi Monte-Carlo sequence (Sobol numbers).
%! @sp 2
%! @strong{Inputs}
%! @sp 1
%! @table @ @var
%! @item d
%! Integer scalar, dimension.
%! @item s
%! Integer scalar (int64), seed.
%! @item t
%! Integer scalar, sequence type:
%! @sp 1
%! @table @ @samp
%! @item @var{t}=0
%! Uniform numbers in a n-dimensional (unit by default) hypercube
%! @item @var{t}=1
%! Gaussian numbers
%! @item @var{t}=2
%! Uniform numbers on a n-dimensional (unit by default) hypersphere
%! @end table
%! @item n
%! Integer scalar, number of elements in the sequence.
%! @item lu
%! Optional argument, the interpretation depends on its size:
%! @sp 1
%! @table @ @samp
%! @item @var{d}x2 array of doubles
%! Lower and upper bounds of the hypercube (default is 0-1 in all dimensions). @var{t} must be equal to zero.
%! @item @var{d}x@var{d} array of doubles
%! Lower cholesky of the covariance matrix of the Gaussian variates (default is the identity matrix). @var{t} must be equal to one.
%! @item scalar double
%! Radius of the hypershere (default is one). @var{t} must be equal to two.
%! @end table
%! @end table
%! @sp 2
%! @strong{Outputs}
%! @sp 1
%! @table @ @var
%! @item a
%! @var{n}x@var{d} matrix of doubles, the Sobol sequence.
%! @item s
%! Integer scalar (int64), current value of the seed.
%! @item info
%! Integer scalar, equal to 1 if mex routine fails, 0 otherwise.
%! @end table
%! @sp 2
%! @strong{This function is called by:}
%! @sp 2
%! @strong{This function calls:}
%! @sp 2
%! @end deftypefn
%@eod:
%@test:1
%$ t = ones(3,1);
%$
%$ d = 2;
%$ n = 100;
%$ s = int64(0);
%$
%$ try
%$ [draws, S] = qmc_sequence(d,s,0,n);
%$ catch
%$ t(1) = 0;
%$ end
%$
%$ try
%$ [draws, S] = qmc_sequence(d,s,1,n);
%$ catch
%$ t(2) = 0;
%$ end
%$
%$ try
%$ [draws, S] = qmc_sequence(d,s,2,n);
%$ catch
%$ t(3) = 0;
%$ end
%$
%$ T = all(t);
%@eof:1
%@test:2
%$ t = ones(3,1);
%$
%$ d = 2;
%$ n = 100;
%$ s = int64(0);
%$
%$ [draws1, S] = qmc_sequence(d,s,0,n);
%$ [draws2, Q] = qmc_sequence(d,S,0,n);
%$ [draws3, P] = qmc_sequence(d,s,0,2*n);
%$
%$ t(1) = dyn_assert(s,int64(0));
%$ t(2) = dyn_assert(P,Q);
%$ t(3) = dyn_assert([draws1,draws2],draws3);
%$ T = all(t);
%@eof:2
%@test:3
%$
%$ d = 2;
%$ n = 100;
%$ s = int64(0);
%$
%$ [draws1, S] = qmc_sequence(d,s,0,n,[0 , 2; -1, 2]);
%$ [draws2, Q] = qmc_sequence(d,s,0,n);
%$
%$ draws3 = draws2;
%$ draws3(1,:) = 2*draws2(1,:);
%$ draws3(2,:) = 3*draws2(2,:)-1;
%$ t(1) = dyn_assert(S,Q);
%$ t(2) = dyn_assert(draws1,draws3);
%$ T = all(t);
%@eof:3
%@test:4
%$
%$ d = 2;
%$ n = 100;
%$ s = int64(0);
%$ radius = pi;
%$
%$ [draws, S] = qmc_sequence(d,s,2,n,radius);
%$
%$ t(1) = dyn_assert(sqrt(draws(:,3)'*draws(:,3)),radius,1e-14);
%$ t(2) = dyn_assert(sqrt(draws(:,5)'*draws(:,5)),radius,1e-14);
%$ t(3) = dyn_assert(sqrt(draws(:,7)'*draws(:,7)),radius,1e-14);
%$ t(4) = dyn_assert(sqrt(draws(:,11)'*draws(:,11)),radius,1e-14);
%$ t(5) = dyn_assert(sqrt(draws(:,13)'*draws(:,13)),radius,1e-14);
%$ t(6) = dyn_assert(sqrt(draws(:,17)'*draws(:,17)),radius,1e-14);
%$ t(7) = dyn_assert(sqrt(draws(:,19)'*draws(:,19)),radius,1e-14);
%$ t(8) = dyn_assert(sqrt(draws(:,23)'*draws(:,23)),radius,1e-14);
%$ t(9) = dyn_assert(sqrt(draws(:,29)'*draws(:,29)),radius,1e-14);
%$ T = all(t);
%@eof:4
%@test:5
%$
%$ d = 2;
%$ n = 100000;
%$ b = 100;
%$ s = int64(5);
%$
%$ covariance = [.4 -.1; -.1 .2];
%$ chol_covariance = transpose(chol(covariance));
%$
%$ draws = [];
%$
%$ for i=1:b
%$ [tmp, s] = qmc_sequence(d,s,1,n,chol_covariance);
%$ draws = [draws, tmp];
%$ end
%$
%$ COVARIANCE = draws*draws'/(b*n);
%$
%$ t(1) = dyn_assert(covariance,COVARIANCE,1e-6);
%$ T = all(t);
%@eof:5

View File

@ -1,9 +1,9 @@
function y_=simult(ys, dr)
% function y_=simult(ys, dr)
function y_=simult(y0, dr)
% function y_=simult(y0, dr)
% Recursive Monte Carlo simulations
%
% INPUTS
% ys: vector of variables in steady state
% y0: vector of variables in initial period of the simulation
% dr: structure of decisions rules for stochastic simulations
%
% OUTPUTS
@ -54,7 +54,7 @@ for i=1:replic
if ~isempty(M_.Sigma_e)
oo_.exo_simul(:,i_exo_var) = randn(options_.periods,nxs)*chol_S;
end
y_ = simult_(ys,dr,oo_.exo_simul,order);
y_ = simult_(y0,dr,oo_.exo_simul,order);
% elimninating initial value
y_ = y_(:,2:end);
if replic > 1

View File

@ -224,8 +224,11 @@ for it_=start:incr:finish
if (verbose == 1)
disp('steady: fsolve');
end
if ~exist('OCTAVE_VERSION') && ~license('test', 'optimization_toolbox')
error('SOLVE_ONE_BOUNDARY: you can''t use solve_algo=0 since you don''t have MATLAB''s Optimization Toolbox')
if ~exist('OCTAVE_VERSION')
[has_optimization_toolbox junk] = license('checkout','optimization_toolbox');
if ~has_optimization_toolbox
error('SOLVE_ONE_BOUNDARY: you can''t use solve_algo=0 since you don''t have MATLAB''s Optimization Toolbox')
end
end
options=optimset('fsolve');
options.MaxFunEvals = 50000;

View File

@ -62,13 +62,3 @@ end
disp_steady_state(M_,oo_);
M_.Sigma_e = Sigma_e;
if isempty(ys0_)
oo_.endo_simul(:,1:M_.maximum_lag) = oo_.steady_state * ones(1, M_.maximum_lag);
%%% Unless I'm wrong, this is (should be?) done in make_y_.m
% $$$ else
% $$$ options_ =set_default_option(options_,'periods',1);
% $$$ oo_.endo_simul(:,M_.maximum_lag+1:M_.maximum_lag+options_.periods+ ...
% $$$ M_.maximum_lead) = oo_.steady_state * ones(1,options_.periods+M_.maximum_lead);
end

View File

@ -133,7 +133,12 @@ if options_.periods > 0 && ~PI_PCL_solver
options_ =options_old;
return
end
oo_.endo_simul = simult(oo_.dr.ys,oo_.dr);
if isempty(M_.endo_histval)
y0 = oo_.dr.ys;
else
y0 = M_.endo_histval;
end
oo_.endo_simul = simult(y0,oo_.dr);
dyn2vec;
end

View File

@ -82,7 +82,12 @@ elseif options_.periods ~= 0
options_ =options_old;
return
end
oo_.endo_simul = simult(repmat(oo_.dr.ys,1,M_.maximum_lag),oo_.dr);
if isempty(M_.endo_histval)
y0 = oo_.dr.ys;
else
y0 = M_.endo_histval;
end
oo_.endo_simul = simult(y0,oo_.dr);
dyn2vec;
if options_.nomoments == 0
disp_moments(oo_.endo_simul,var_list);

View File

@ -1,4 +1,4 @@
function dynTest(fun)
function dynTest(fun,dynare_path)
%@info:
%! @deftypefn {Function File} dynTest (@var{fun})
@ -50,25 +50,32 @@ function dynTest(fun)
% You should have received a copy of the GNU General Public License
% along with Dynare. If not, see <http://www.gnu.org/licenses/>.
if isempty(strfind(fun,'@')) & (~isempty(strfind(fun,'/')) || ~isempty(strfind(fun,'\')) )
[pathstr1, name, ext] = fileparts(fun);
addpath(pathstr1);
rm_path = 1;
else
rm_path = 0;
original_directory = pwd();
cd([dynare_path filesep '..' filesep 'tests']);
[pathstr1, name1, ext1] = fileparts(fun);
mex_flag = 0;
if exist(name1)==3
mex_flag = 1;
end
[pathstr2, name, ext] = fileparts(which(fun));
class_flag = 0;
if ~isempty(strfind(fun,'@')) || ~isempty(strfind(which(name1),'@'))
class_flag = 1;
end
if ~( isempty(pathstr2) || isempty(name) || isempty(ext) ) && strcmp(ext(2:end),'m')
check = mtest(name,pathstr2);
if check
disp(['Succesfull test(s) for ' fun ' routine!'])
check = mtest(name1,pathstr1);
if check
if mex_flag
disp(['Succesfull test(s) for ' name1 ' mex file!'])
elseif class_flag
disp(['Succesfull test(s) for ' name1 ' method!'])
else
disp(['Succesfull test(s) for ' name1 ' routine!'])
end
else
disp([fun 'is not a known matlab/octave routine!'])
end
if rm_path
rmpath(pathstr1)
end
cd(original_directory);

View File

@ -2,7 +2,7 @@ ACLOCAL_AMFLAGS = -I ../../../m4
# libdynare++ must come before gensylv, k_order_perturbation, dynare_simul_
if DO_SOMETHING
SUBDIRS = mjdgges kronecker bytecode libdynare++ gensylv k_order_perturbation dynare_simul_ estimation block_kalman_filter
SUBDIRS = mjdgges kronecker bytecode libdynare++ gensylv k_order_perturbation dynare_simul_ estimation block_kalman_filter sobol
if HAVE_GSL
SUBDIRS += ms_sbvar

View File

@ -143,6 +143,7 @@ AC_CONFIG_FILES([Makefile
libslicot/Makefile
kalman_steady_state/Makefile
ms_sbvar/Makefile
block_kalman_filter/Makefile])
block_kalman_filter/Makefile
sobol/Makefile])
AC_OUTPUT

View File

@ -0,0 +1,2 @@
include ../mex.am
include ../../qmc_sequence.am

View File

@ -2,7 +2,7 @@ ACLOCAL_AMFLAGS = -I ../../../m4
# libdynare++ must come before gensylv, k_order_perturbation, dynare_simul_
if DO_SOMETHING
SUBDIRS = mjdgges kronecker bytecode libdynare++ gensylv k_order_perturbation dynare_simul_ qzcomplex ordschur block_kalman_filter
SUBDIRS = mjdgges kronecker bytecode libdynare++ gensylv k_order_perturbation dynare_simul_ qzcomplex ordschur block_kalman_filter sobol
if HAVE_GSL
SUBDIRS += ms_sbvar

View File

@ -130,6 +130,7 @@ AC_CONFIG_FILES([Makefile
libslicot/Makefile
kalman_steady_state/Makefile
ms_sbvar/Makefile
block_kalman_filter/Makefile])
block_kalman_filter/Makefile
sobol/Makefile])
AC_OUTPUT

View File

@ -14,10 +14,10 @@ LDFLAGS += $(shell $(MKOCTFILE) -p LDFLAGS)
LIBS += $(shell $(MKOCTFILE) -p OCTAVE_LIBS)
LIBS += $(shell $(MKOCTFILE) -p BLAS_LIBS)
LIBS += $(shell $(MKOCTFILE) -p LAPACK_LIBS)
LIBS += $(shell $(MKOCTFILE) -p FFTW_LIBS)
LIBS += $(shell $(MKOCTFILE) -p LIBS)
LIBS += $(shell $(MKOCTFILE) -p FLIBS)
LIBS += $(shell $(MKOCTFILE) -p CXXLIBS) # Only used for Octave/MinGW
all-local:
$(MKDIR_P) $(top_srcdir)/../../octave

View File

@ -0,0 +1,3 @@
EXEEXT = .mex
include ../mex.am
include ../../qmc_sequence.am

View File

@ -0,0 +1,5 @@
vpath %.cc $(top_srcdir)/../../sources/sobol
noinst_PROGRAMS = qmc_sequence
nodist_qmc_sequence_SOURCES = qmc_sequence.cc

View File

@ -13,7 +13,8 @@ EXTRA_DIST = \
libslicot \
kalman_steady_state \
ms-sbvar \
block_kalman_filter
block_kalman_filter \
sobol
clean-local:
rm -rf `find mex/sources -name *.o`

View File

@ -286,8 +286,7 @@ main(int nrhs, const char *prhs[])
nb_row_xd = row_x;
}
}
mxArray *ep = mxGetFieldByNumber(options_, 0, mxGetFieldNumber(options_, "ep"));
int verbose= int(*mxGetPr((mxGetFieldByNumber(ep, 0, mxGetFieldNumber(ep, "verbosity")))));
int verbose= int(*mxGetPr((mxGetFieldByNumber(options_, 0, mxGetFieldNumber(options_, "verbosity")))));
if (verbose)
print_it = true;
int maxit_ = int (floor(*(mxGetPr(mxGetFieldByNumber(options_, 0, mxGetFieldNumber(options_, "maxit_"))))));

View File

@ -0,0 +1,195 @@
/* Generates gaussian random deviates from uniform random deviates.
**
** Pseudo code of the algorithm is given at http://home.online.no/~pjacklam/notes/invnorm
**
** Copyright (C) 2010 Dynare Team
**
** This file is part of Dynare.
**
** Dynare is free software: you can redistribute it and/or modify
** it under the terms of the GNU General Public License as published by
** the Free Software Foundation, either version 3 of the License, or
** (at your option) any later version.
**
** Dynare is distributed in the hope that it will be useful,
** but WITHOUT ANY WARRANTY; without even the implied warranty of
** MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
** GNU General Public License for more details.
**
** You should have received a copy of the GNU General Public License
** along with Dynare. If not, see <http://www.gnu.org/licenses/>.
**
** AUTHOR(S): stephane DOT adjemian AT univ DASH lemans DOT fr
*/
#include <cstdlib>
#include <iostream>
#include <iomanip>
#include <cmath>
#include <ctime>
#include <dynblas.h>
using namespace std;
#define lb .02425
#define ub .97575
#ifdef USE_OMP
# include <omp.h>
#endif
#define DEBUG_OMP 0
template<typename T> T icdf( const T uniform )
/*
** This function invert the gaussian cumulative distribution function.
**
*/
{
static T A[6] =
{
-3.969683028665376e+01,
2.209460984245205e+02,
-2.759285104469687e+02,
1.383577518672690e+02,
-3.066479806614716e+01,
2.506628277459239e+00
};
static T B[5] =
{
-5.447609879822406e+01,
1.615858368580409e+02,
-1.556989798598866e+02,
6.680131188771972e+01,
-1.328068155288572e+01
};
static T C[6] =
{
-7.784894002430293e-03,
-3.223964580411365e-01,
-2.400758277161838e+00,
-2.549732539343734e+00,
4.374664141464968e+00,
2.938163982698783e+00
};
static T D[4] =
{
7.784695709041462e-03,
3.224671290700398e-01,
2.445134137142996e+00,
3.754408661907416e+00
};
T gaussian = (T)0.0;
if ( (0<uniform) && (uniform<lb) )
{
T tmp;
tmp = sqrt(-2*log(uniform));
gaussian = (((((C[0]*tmp+C[1])*tmp+C[2])*tmp+C[3])*tmp+C[4])*tmp+C[5])/((((D[0]*tmp+D[1])*tmp+D[2])*tmp+D[3])*tmp+1);
}
else
{
if ((lb <= uniform) && (uniform <= ub))
{
T tmp, TMP;
tmp = uniform - .5;
TMP = tmp*tmp;
gaussian = (((((A[0]*TMP+A[1])*TMP+A[2])*TMP+A[3])*TMP+A[4])*TMP+A[5])*tmp/(((((B[0]*TMP+B[1])*TMP+B[2])*TMP+B[3])*TMP+B[4])*TMP+1);
}
else
{
if ((ub < uniform) && (uniform < 1))
{
T tmp;
tmp = sqrt(-2*log(1-uniform));
gaussian = -(((((C[0]*tmp+C[1])*tmp+C[2])*tmp+C[3])*tmp+C[4])*tmp+C[5])/((((D[0]*tmp+D[1])*tmp+D[2])*tmp+D[3])*tmp+1);
}
}
}
if ( (0<uniform) && (uniform<1) )
{
T tmp, tmp_;
tmp = .5*erfc(-gaussian/sqrt(2.0))-uniform;
tmp_ = tmp*sqrt(2*M_PI)*exp(.5*gaussian*gaussian);
gaussian = gaussian - tmp_/(1+.5*gaussian*tmp_);
}
if ( uniform==0)
{
gaussian = -INFINITY;
}
if ( uniform==1)
{
gaussian = INFINITY;
}
return(gaussian);
}
template<typename T> void icdfm( const int n, T *U)
{
#if USE_OMP
#pragma omp parallel for num_threads(omp_get_num_threads())
#endif
for(int i=0; i<n; i++)
{
U[i] = icdf(U[i]);
}
return;
}
template<typename T> void icdfmSigma( const int d, const int n, T *U, const double *LowerCholSigma)
{
double one = 1.0;
double zero = 0.0;
blas_int dd(d);
blas_int nn(n);
icdfm(n*d, U);
double tmp[n*d];
dgemm("N","N",&dd,&nn,&dd,&one,LowerCholSigma,&dd,U,&dd,&zero,tmp,&dd);
memcpy(U,tmp,d*n*sizeof(double));
return;
}
template<typename T> void usphere( const int d, const int n, T *U)
{
icdfm(n*d, U);
#if USE_OMP
#pragma omp parallel for num_threads(omp_get_num_threads())
#endif
for (int j=0; j<n; j++)// sequence index.
{
int k = j*d;
double norm = 0.0;
for(int i=0; i<d; i++)// dimension index.
{
norm = norm + U[k+i]*U[k+i];
}
norm = sqrt(norm);
for(int i=0; i<d; i++)// dimension index.
{
U[k+i] = U[k+i]/norm;
}
}
return;
}
template<typename T> void usphereRadius( const int d, const int n, double radius, T *U)
{
icdfm(n*d, U);
#if USE_OMP
#pragma omp parallel for num_threads(omp_get_num_threads())
#endif
for (int j=0; j<n; j++)// sequence index.
{
int k = j*d;
double norm = 0.0;
for(int i=0; i<d; i++)// dimension index.
{
norm = norm + U[k+i]*U[k+i];
}
norm = sqrt(norm);
for(int i=0; i<d; i++)// dimension index.
{
U[k+i] = radius*U[k+i]/norm;
}
}
return;
}

View File

@ -0,0 +1,28 @@
template<typename T> int initialize_v_array (int dim_max, int log_max, T **v)
/*
** This function initializes the v array used in the sobol routine.
**
** Original files downloaded from http://people.sc.fsu.edu/~burkardt/cpp_src/sobol/ (version 17-Feb-2009 09:46)
**
** Copyright (C) 2009 John Burkardt
** Copyright (C) 2010-2011 Dynare Team
**
** This program is free software: you can redistribute it and/or modify
** it under the terms of the GNU Lesser General Public License as published by
** the Free Software Foundation, either version 3 of the License, or
** (at your option) any later version.
**
** This program 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 Lesser General Public License for more details.
**
** You should have received a copy of the GNU Lesser General Public License
** along with this program. If not, see <http://www.gnu.org/licenses/>.
**
** AUTHOR(S): stephane DOT adjemian AT univ DASH lemans DOT fr
*/
{
#include "initialize_v_array.inc"
return 1;
}

File diff suppressed because it is too large Load Diff

View File

@ -0,0 +1,212 @@
/*
** Computes Quasi Monte-Carlo sequence.
**
** Copyright (C) 2010-2011 Dynare Team
**
** This file is part of Dynare (can be used outside 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/>.
**
** AUTHOR(S): stephane DOT adjemian AT univ DASH lemans DOT fr
**/
#include <string.h>
#include <stdint.h>
#include <dynmex.h>
#include "sobol.hh"
#include "gaussian.hh"
void mexFunction( int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])
{
/*
** INPUTS:
** prhs[0] [integer] scalar, dimension.
** prhs[1] [integer] scalar, seed.
** prhs[2] [integer] scalar, sequence type:
** 0 ==> uniform,
** 1 ==> gaussian,
** 2 ==> uniform on an hypershere.
** prhs[3] [integer] scalar, sequence size.
** prhs[4] [double] dimension*2 array, lower and upper bounds of the hypercube (default is 0-1 in all dimensions) if prhs[2]==0,
** dimension*dimension array, covariance of the multivariate gaussian distribution of prhs[2]==1 (default is the identity matrix),
** scalar, radius of the hypershere if prhs[2]==2 (default is one).
**
** OUTPUTS:
** plhs[0] [double] sequence_size*dimension array, the Sobol sequence.
** plhs[1] [integer] scalar, seed.
** plhs[2] [integer] zero in case of success, one in case of error
**
*/
/*
** Check the number of input and output arguments.
*/
if ( !( (nrhs==5) | (nrhs==4) | (nrhs==3) ) )
{
DYN_MEX_FUNC_ERR_MSG_TXT("qmc_sequence:: Five, four or three input arguments are required!");
}
if (nlhs == 0)
{
DYN_MEX_FUNC_ERR_MSG_TXT("qmc_sequence:: At least one output argument is required!");
}
/*
** Test the first input argument and assign it to dimension.
*/
if ( !( mxIsNumeric(prhs[0]) ) )
{
DYN_MEX_FUNC_ERR_MSG_TXT("qmc_sequence:: First input (dimension) has to be a positive integer!");
}
int dimension = ( int ) mxGetScalar(prhs[0]);
/*
** Test the second input argument and assign it to seed.
*/
if ( !( mxIsNumeric(prhs[1]) && mxIsClass(prhs[1],"int64") ) )
{
DYN_MEX_FUNC_ERR_MSG_TXT("qmc_sequence:: Second input (seed) has to be an integer [int64]!");
}
int64_T seed = ( int64_T ) mxGetScalar( prhs[1] );
/*
** Test the third input argument and assign it to type (kind of QMC sequence).
*/
int error_flag_3 = 0;
if ( !(mxIsNumeric(prhs[2])) )
{
error_flag_3 = 1;
}
int type = ( int ) mxGetScalar(prhs[2]);
if ( !(type==0 || type==1 || type==2) )
{
error_flag_3 = 1;
}
if (error_flag_3==1)
{
DYN_MEX_FUNC_ERR_MSG_TXT("qmc_sequence:: Third input (type of QMC sequence) has to be an integer equal to 0, 1 or 2!");
}
/*
** Test dimension>=2 when type==2
*/
if ( (type==2) && (dimension<2) )
{
DYN_MEX_FUNC_ERR_MSG_TXT("qmc_sequence:: First input (dimension) has to be greater than 1 for a uniform QMC on an hypershere!");
}
/*
** Test the optional fourth input argument and assign it to sequence_size.
*/
if ( ( nrhs>3 ) && !mxIsNumeric(prhs[3]) )
{
DYN_MEX_FUNC_ERR_MSG_TXT("qmc_sequence:: Fourth input (qmc sequence size) has to be a positive integer!");
}
int sequence_size;
if ( nrhs>3)
{
sequence_size = ( int ) mxGetScalar( prhs[3] );
}
else
{
sequence_size = 1;
}
/*
** Test the optional fifth input argument and assign it to lower_and_upper_bounds.
*/
if ( ( nrhs>4 ) && (type==0) && ( !( mxGetN(prhs[4])==2) ) )// Sequence of uniformly distributed numbers in an hypercube.
{
DYN_MEX_FUNC_ERR_MSG_TXT("qmc_sequence:: The fifth input argument must be an array with two columns!");
}
if ( (nrhs>4) && (type==0) && ( !( (int)mxGetM(prhs[4])==dimension) ) )
{
DYN_MEX_FUNC_ERR_MSG_TXT("qmc_sequence:: The fourth input argument must be an array with a number of lines equal to dimension (first input argument)!");
}
if ( ( nrhs>4 ) && (type==1) && ( !( ((int)mxGetN(prhs[4])==dimension) && ((int)mxGetM(prhs[4])==dimension) ) ) )// Sequence of normally distributed numbers.
{
DYN_MEX_FUNC_ERR_MSG_TXT("qmc_sequence:: The fifth input argument must be a squared matrix (whose dimension is given by the first input argument)!");
}
if ( ( nrhs>4 ) && (type==2) && ( !( (mxGetN(prhs[4])==1) && (mxGetM(prhs[4])==1) ) ) )// Sequence of uniformly distributed numbers on an hypershere.
{
DYN_MEX_FUNC_ERR_MSG_TXT("qmc_sequence:: The fifth input argument must be a positive scalar!");
}
double *lower_bounds, *upper_bounds;
int unit_hypercube_flag = 1;
if ( (type==0) && (nrhs>4) )
{
lower_bounds = (double *) mxCalloc(dimension,sizeof(double));
upper_bounds = (double *) mxCalloc(dimension,sizeof(double));
double *tmp;
tmp = (double *) mxCalloc(dimension*2,sizeof(double));
memcpy(tmp,mxGetPr(prhs[4]),dimension*2*sizeof(double));
lower_bounds = &tmp[0];
upper_bounds = &tmp[dimension];
unit_hypercube_flag = 0;
}
double *cholcov;
int identity_covariance_matrix = 1;
if ( (type==1) && (nrhs>4) )
{
cholcov = (double *) mxCalloc(dimension*dimension,sizeof(double));
double *tmp;
tmp = (double *) mxCalloc(dimension*dimension,sizeof(double));
memcpy(tmp,mxGetPr(prhs[4]),dimension*dimension*sizeof(double));
cholcov = &tmp[0];
identity_covariance_matrix = 0;
}
double radius = 1.0;
int unit_radius = 1;
if ( (type==2) && (nrhs>4) )
{
double *tmp;
tmp = (double *) mxCalloc(1,sizeof(double));
memcpy(tmp,mxGetPr(prhs[4]),sizeof(double));
radius = tmp[0];
unit_radius = 0;
}
/*
** Initialize outputs of the mex file.
*/
double *qmc_draws;
plhs[0] = mxCreateDoubleMatrix(dimension,sequence_size,mxREAL);
qmc_draws = mxGetPr(plhs[0]);
int64_T seed_out;
if (sequence_size==1)
{
next_sobol ( dimension, &seed, qmc_draws );
seed_out = seed;
}
else
seed_out = sobol_block( dimension, sequence_size, seed, qmc_draws);
if (type==0 && unit_hypercube_flag==0)// Uniform QMC sequence in an hypercube.
expand_unit_hypercube( dimension, sequence_size, qmc_draws, lower_bounds, upper_bounds);
else if (type==1)// Normal QMC sequance in R^n.
{
if (identity_covariance_matrix==1)
icdfm(dimension*sequence_size, qmc_draws);
else
icdfmSigma(dimension,sequence_size, qmc_draws, cholcov);
}
else if (type==2)// Uniform QMC sequence on an hypershere.
{
if (unit_radius==1)
usphere(dimension, sequence_size, qmc_draws);
else
usphereRadius(dimension, sequence_size, radius, qmc_draws);
}
if (nlhs >= 2)
{
plhs[1] = mxCreateNumericMatrix(1, 1, mxINT64_CLASS, mxREAL);
*((int64_T *) mxGetData(plhs[1])) = seed_out;
}
if (nlhs >= 3)
plhs[2] = mxCreateDoubleScalar(0);
}

574
mex/sources/sobol/sobol.hh Normal file
View File

@ -0,0 +1,574 @@
/* Quasi Monte Carlo sequences (à la Sobol).
**
** Original files downloaded from http://people.sc.fsu.edu/~burkardt/cpp_src/sobol/ (version 17-Feb-2009 09:46)
**
** Copyright (C) 2009 John Burkardt
** Copyright (C) 2010-2011 Dynare Team
**
** This program is free software: you can redistribute it and/or modify
** it under the terms of the GNU Lesser General Public License as published by
** the Free Software Foundation, either version 3 of the License, or
** (at your option) any later version.
**
** This program 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 Lesser General Public License for more details.
**
** You should have received a copy of the GNU Lesser General Public License
** along with this program. If not, see <http://www.gnu.org/licenses/>.
**
** AUTHOR(S): stephane DOT adjemian AT univ DASH lemans DOT fr
*/
#include <cstdlib>
#include <iostream>
#include <iomanip>
#include <cmath>
#include <ctime>
#include "initialize_v_array.hh"
using namespace std;
#define DIM_MAX 1111
template<typename T> int bit_hi1(T n)
/*
** This function returns the position of the high 1 bit base 2 in an integer.
**
** Example:
**
** N Binary Hi 1
** ---- -------- ----
** 0 0 0
** 1 1 1
** 2 10 2
** 3 11 2
** 4 100 3
** 5 101 3
** 6 110 3
** 7 111 3
** 8 1000 4
** 9 1001 4
** 10 1010 4
** 11 1011 4
** 12 1100 4
** 13 1101 4
** 14 1110 4
** 15 1111 4
** 16 10000 5
** 17 10001 5
** 1023 1111111111 10
** 1024 10000000000 11
** 1025 10000000001 11
**
**
** Original files downloaded from http://people.sc.fsu.edu/~burkardt/cpp_src/sobol/ (version 17-Feb-2009 09:46)
**
** Input, int or long long, the integer to be measured.
** N should be nonnegative. If N is nonpositive, BIT_HI1 will always be 0.
**
** Output: the location of the high order bit.
*/
{
int bit = 0 ;
while ( n > 0 )
{
bit++ ;
n = n/2 ;
}
return bit ;
}
template<typename T> int bit_lo0 ( T n )
/*
** This function returns the position of the low 0 bit base 2 in an integer.
**
** Example:
**
** N Binary Lo 0
** ---- -------- ----
** 0 0 1
** 1 1 2
** 2 10 1
** 3 11 3
** 4 100 1
** 5 101 2
** 6 110 1
** 7 111 4
** 8 1000 1
** 9 1001 2
** 10 1010 1
** 11 1011 3
** 12 1100 1
** 13 1101 2
** 14 1110 1
** 15 1111 5
** 16 10000 1
** 17 10001 2
** 1023 1111111111 1
** 1024 10000000000 1
** 1025 10000000001 1
**
**
** Original files downloaded from http://people.sc.fsu.edu/~burkardt/cpp_src/sobol/ (version 17-Feb-2009 09:46)
**
** INPUTS
**
** Input, int N, the integer to be measured.
** N should be nonnegative.
**
** OUTPUTS (int) the position of the low 0 bit.
*/
{
int bit = 0;
while ( true )
{
bit++;
T n2 = n/2;
if ( n == 2*n2 )
{
break;
}
n = n2;
}
return bit;
}
template<typename T> T ixor ( T i, T j )
/*
** This function calculates the exclusive OR of two integers.
**
** Original files downloaded from http://people.sc.fsu.edu/~burkardt/cpp_src/sobol/ (version 17-Feb-2009 09:46)
**
** INPUTS I, J, two integer to be exclusive OR-ed.
**
** OUTPUTS (integer) the exclusive OR of I and J.
*/
{
T k = 0;
T l = 1;
while ( i != 0 || j != 0 )
{
T i2 = i / 2;
T j2 = j / 2;
if (
( ( i == 2 * i2 ) && ( j != 2 * j2 ) ) ||
( ( i != 2 * i2 ) && ( j == 2 * j2 ) ) )
{
k = k + l;
}
i = i2;
j = j2;
l = 2 * l;
}
return k;
}
template<typename T1, typename T2> void next_sobol ( int dim_num, T1 *seed, T2 *quasi )
/*
** This function generates a new quasirandom Sobol vector with each call.
**
** Discussion:
**
** The routine adapts the ideas of Antonov and Saleev.
**
** This routine uses LONG LONG INT for integers and DOUBLE for real values or
** INT for integers and FLOAT for real values.
**
** Thanks to Steffan Berridge for supplying (twice) the properly
** formatted V data needed to extend the original routine's dimension
** limit from 40 to 1111, 05 June 2007.
**
** Thanks to Francis Dalaudier for pointing out that the range of allowed
** values of DIM_NUM should start at 1, not 2! 17 February 2009.
**
** Original files downloaded from http://people.sc.fsu.edu/~burkardt/cpp_src/sobol/ (version 17-Feb-2009 09:46)
**
** Reference:
**
** IA Antonov, VM Saleev,
** An Economic Method of Computing LP Tau-Sequences,
** USSR Computational Mathematics and Mathematical Physics,
** Volume 19, 1980, pages 252 - 256.
**
** Paul Bratley, Bennett Fox,
** Algorithm 659:
** Implementing Sobol's Quasirandom Sequence Generator,
** ACM Transactions on Mathematical Software,
** Volume 14, Number 1, pages 88-100, 1988.
**
** Bennett Fox,
** Algorithm 647:
** Implementation and Relative Efficiency of Quasirandom
** Sequence Generators,
** ACM Transactions on Mathematical Software,
** Volume 12, Number 4, pages 362-376, 1986.
**
** Stephen Joe, Frances Kuo
** Remark on Algorithm 659:
** Implementing Sobol's Quasirandom Sequence Generator,
** ACM Transactions on Mathematical Software,
** Volume 29, Number 1, pages 49-57, March 2003.
**
** Ilya Sobol,
** USSR Computational Mathematics and Mathematical Physics,
** Volume 16, pages 236-242, 1977.
**
** Ilya Sobol, YL Levitan,
** The Production of Points Uniformly Distributed in a Multidimensional
** Cube (in Russian),
** Preprint IPM Akad. Nauk SSSR,
** Number 40, Moscow 1976.
**
** Parameters:
**
** Input, int DIM_NUM, the number of spatial dimensions.
** DIM_NUM must satisfy 1 <= DIM_NUM <= 1111.
**
** Input/output, long long int *SEED, the "seed" for the sequence.
** This is essentially the index in the sequence of the quasirandom
** value to be generated. On output, SEED has been set to the
** appropriate next value, usually simply SEED+1.
** If SEED is less than 0 on input, it is treated as though it were 0.
** An input value of 0 requests the first (0-th) element of the sequence.
**
** Output, double QUASI[DIM_NUM], the next quasirandom vector.
*/
{
static T1 atmost ;
static int dim_num_save = 0 ;
int LOG_MAX = sizeof(T1)*8-2 ;
bool includ[LOG_MAX];
static bool initialized = false;
static T1 lastq[DIM_MAX];
static T1 maxcol;
T1 l = 0;
static T1 poly[DIM_MAX] =
{
1, 3, 7, 11, 13, 19, 25, 37, 59, 47,
61, 55, 41, 67, 97, 91, 109, 103, 115, 131,
193, 137, 145, 143, 241, 157, 185, 167, 229, 171,
213, 191, 253, 203, 211, 239, 247, 285, 369, 299,
301, 333, 351, 355, 357, 361, 391, 397, 425, 451,
463, 487, 501, 529, 539, 545, 557, 563, 601, 607,
617, 623, 631, 637, 647, 661, 675, 677, 687, 695,
701, 719, 721, 731, 757, 761, 787, 789, 799, 803,
817, 827, 847, 859, 865, 875, 877, 883, 895, 901,
911, 949, 953, 967, 971, 973, 981, 985, 995, 1001,
1019, 1033, 1051, 1063, 1069, 1125, 1135, 1153, 1163, 1221,
1239, 1255, 1267, 1279, 1293, 1305, 1315, 1329, 1341, 1347,
1367, 1387, 1413, 1423, 1431, 1441, 1479, 1509, 1527, 1531,
1555, 1557, 1573, 1591, 1603, 1615, 1627, 1657, 1663, 1673,
1717, 1729, 1747, 1759, 1789, 1815, 1821, 1825, 1849, 1863,
1869, 1877, 1881, 1891, 1917, 1933, 1939, 1969, 2011, 2035,
2041, 2053, 2071, 2091, 2093, 2119, 2147, 2149, 2161, 2171,
2189, 2197, 2207, 2217, 2225, 2255, 2257, 2273, 2279, 2283,
2293, 2317, 2323, 2341, 2345, 2363, 2365, 2373, 2377, 2385,
2395, 2419, 2421, 2431, 2435, 2447, 2475, 2477, 2489, 2503,
2521, 2533, 2551, 2561, 2567, 2579, 2581, 2601, 2633, 2657,
2669, 2681, 2687, 2693, 2705, 2717, 2727, 2731, 2739, 2741,
2773, 2783, 2793, 2799, 2801, 2811, 2819, 2825, 2833, 2867,
2879, 2881, 2891, 2905, 2911, 2917, 2927, 2941, 2951, 2955,
2963, 2965, 2991, 2999, 3005, 3017, 3035, 3037, 3047, 3053,
3083, 3085, 3097, 3103, 3159, 3169, 3179, 3187, 3205, 3209,
3223, 3227, 3229, 3251, 3263, 3271, 3277, 3283, 3285, 3299,
3305, 3319, 3331, 3343, 3357, 3367, 3373, 3393, 3399, 3413,
3417, 3427, 3439, 3441, 3475, 3487, 3497, 3515, 3517, 3529,
3543, 3547, 3553, 3559, 3573, 3589, 3613, 3617, 3623, 3627,
3635, 3641, 3655, 3659, 3669, 3679, 3697, 3707, 3709, 3713,
3731, 3743, 3747, 3771, 3791, 3805, 3827, 3833, 3851, 3865,
3889, 3895, 3933, 3947, 3949, 3957, 3971, 3985, 3991, 3995,
4007, 4013, 4021, 4045, 4051, 4069, 4073, 4179, 4201, 4219,
4221, 4249, 4305, 4331, 4359, 4383, 4387, 4411, 4431, 4439,
4449, 4459, 4485, 4531, 4569, 4575, 4621, 4663, 4669, 4711,
4723, 4735, 4793, 4801, 4811, 4879, 4893, 4897, 4921, 4927,
4941, 4977, 5017, 5027, 5033, 5127, 5169, 5175, 5199, 5213,
5223, 5237, 5287, 5293, 5331, 5391, 5405, 5453, 5523, 5573,
5591, 5597, 5611, 5641, 5703, 5717, 5721, 5797, 5821, 5909,
5913, 5955, 5957, 6005, 6025, 6061, 6067, 6079, 6081, 6231,
6237, 6289, 6295, 6329, 6383, 6427, 6453, 6465, 6501, 6523,
6539, 6577, 6589, 6601, 6607, 6631, 6683, 6699, 6707, 6761,
6795, 6865, 6881, 6901, 6923, 6931, 6943, 6999, 7057, 7079,
7103, 7105, 7123, 7173, 7185, 7191, 7207, 7245, 7303, 7327,
7333, 7355, 7365, 7369, 7375, 7411, 7431, 7459, 7491, 7505,
7515, 7541, 7557, 7561, 7701, 7705, 7727, 7749, 7761, 7783,
7795, 7823, 7907, 7953, 7963, 7975, 8049, 8089, 8123, 8125,
8137, 8219, 8231, 8245, 8275, 8293, 8303, 8331, 8333, 8351,
8357, 8367, 8379, 8381, 8387, 8393, 8417, 8435, 8461, 8469,
8489, 8495, 8507, 8515, 8551, 8555, 8569, 8585, 8599, 8605,
8639, 8641, 8647, 8653, 8671, 8675, 8689, 8699, 8729, 8741,
8759, 8765, 8771, 8795, 8797, 8825, 8831, 8841, 8855, 8859,
8883, 8895, 8909, 8943, 8951, 8955, 8965, 8999, 9003, 9031,
9045, 9049, 9071, 9073, 9085, 9095, 9101, 9109, 9123, 9129,
9137, 9143, 9147, 9185, 9197, 9209, 9227, 9235, 9247, 9253,
9257, 9277, 9297, 9303, 9313, 9325, 9343, 9347, 9371, 9373,
9397, 9407, 9409, 9415, 9419, 9443, 9481, 9495, 9501, 9505,
9517, 9529, 9555, 9557, 9571, 9585, 9591, 9607, 9611, 9621,
9625, 9631, 9647, 9661, 9669, 9679, 9687, 9707, 9731, 9733,
9745, 9773, 9791, 9803, 9811, 9817, 9833, 9847, 9851, 9863,
9875, 9881, 9905, 9911, 9917, 9923, 9963, 9973,10003,10025,
10043,10063,10071,10077,10091,10099,10105,10115,10129,10145,
10169,10183,10187,10207,10223,10225,10247,10265,10271,10275,
10289,10299,10301,10309,10343,10357,10373,10411,10413,10431,
10445,10453,10463,10467,10473,10491,10505,10511,10513,10523,
10539,10549,10559,10561,10571,10581,10615,10621,10625,10643,
10655,10671,10679,10685,10691,10711,10739,10741,10755,10767,
10781,10785,10803,10805,10829,10857,10863,10865,10875,10877,
10917,10921,10929,10949,10967,10971,10987,10995,11009,11029,
11043,11045,11055,11063,11075,11081,11117,11135,11141,11159,
11163,11181,11187,11225,11237,11261,11279,11297,11307,11309,
11327,11329,11341,11377,11403,11405,11413,11427,11439,11453,
11461,11473,11479,11489,11495,11499,11533,11545,11561,11567,
11575,11579,11589,11611,11623,11637,11657,11663,11687,11691,
11701,11747,11761,11773,11783,11795,11797,11817,11849,11855,
11867,11869,11873,11883,11919,11921,11927,11933,11947,11955,
11961,11999,12027,12029,12037,12041,12049,12055,12095,12097,
12107,12109,12121,12127,12133,12137,12181,12197,12207,12209,
12239,12253,12263,12269,12277,12287,12295,12309,12313,12335,
12361,12367,12391,12409,12415,12433,12449,12469,12479,12481,
12499,12505,12517,12527,12549,12559,12597,12615,12621,12639,
12643,12657,12667,12707,12713,12727,12741,12745,12763,12769,
12779,12781,12787,12799,12809,12815,12829,12839,12857,12875,
12883,12889,12901,12929,12947,12953,12959,12969,12983,12987,
12995,13015,13019,13031,13063,13077,13103,13137,13149,13173,
13207,13211,13227,13241,13249,13255,13269,13283,13285,13303,
13307,13321,13339,13351,13377,13389,13407,13417,13431,13435,
13447,13459,13465,13477,13501,13513,13531,13543,13561,13581,
13599,13605,13617,13623,13637,13647,13661,13677,13683,13695,
13725,13729,13753,13773,13781,13785,13795,13801,13807,13825,
13835,13855,13861,13871,13883,13897,13905,13915,13939,13941,
13969,13979,13981,13997,14027,14035,14037,14051,14063,14085,
14095,14107,14113,14125,14137,14145,14151,14163,14193,14199,
14219,14229,14233,14243,14277,14287,14289,14295,14301,14305,
14323,14339,14341,14359,14365,14375,14387,14411,14425,14441,
14449,14499,14513,14523,14537,14543,14561,14579,14585,14593,
14599,14603,14611,14641,14671,14695,14701,14723,14725,14743,
14753,14759,14765,14795,14797,14803,14831,14839,14845,14855,
14889,14895,14909,14929,14941,14945,14951,14963,14965,14985,
15033,15039,15053,15059,15061,15071,15077,15081,15099,15121,
15147,15149,15157,15167,15187,15193,15203,15205,15215,15217,
15223,15243,15257,15269,15273,15287,15291,15313,15335,15347,
15359,15373,15379,15381,15391,15395,15397,15419,15439,15453,
15469,15491,15503,15517,15527,15531,15545,15559,15593,15611,
15613,15619,15639,15643,15649,15661,15667,15669,15681,15693,
15717,15721,15741,15745,15765,15793,15799,15811,15825,15835,
15847,15851,15865,15877,15881,15887,15899,15915,15935,15937,
15955,15973,15977,16011,16035,16061,16069,16087,16093,16097,
16121,16141,16153,16159,16165,16183,16189,16195,16197,16201,
16209,16215,16225,16259,16265,16273,16299,16309,16355,16375,
16381 };
static T2 recipd;
static T1 seed_save = - 1;
static T1** v ;
if ( !initialized || dim_num != dim_num_save )
{
v = new T1 *[DIM_MAX] ;
for( int i = 0 ; i < DIM_MAX ; i++ )
v[i] = new T1[LOG_MAX];
initialized = true;
initialize_v_array(DIM_MAX, LOG_MAX, v);
/*
** Check parameters.
*/
if ( dim_num < 1 || DIM_MAX < dim_num )
{
cout << "\n";
cout << "NEXT_SOBOL - Fatal error!\n";
cout << " The spatial dimension DIM_NUM should satisfy:\n";
cout << " 1 <= DIM_NUM <= " << DIM_MAX << "\n";
cout << " But this input value is DIM_NUM = " << dim_num << "\n";
exit ( 1 );
}
dim_num_save = dim_num;
/*
** Set ATMOST = 2^LOG_MAX - 1.
*/
atmost = (T1) 0;
for ( int i = 1; i <= LOG_MAX; i++ )
atmost = 2 * atmost + 1;
/*
** Find the highest 1 bit in ATMOST (should be LOG_MAX).
*/
maxcol = bit_hi1 ( atmost );
/*
** Initialize row 1 of V.
*/
for ( T1 j = 0; j < maxcol; j++ )
{
v[0][j] = (T1) 1;
}
/*
** Initialize the remaining rows of V.
*/
for ( int i = 1; i < dim_num; i++ )
{
/*
** The bit pattern of the integer POLY(I) gives the form
** of polynomial I.
**
** Find the degree of polynomial I from binary encoding.
*/
T1 j = poly[i];
T1 m = 0;
while ( true )
{
j = j / 2;
if ( j <= 0 )
{
break;
}
m = m + 1;
}
/*
** We expand this bit pattern to separate components
** of the logical array INCLUD.
*/
j = poly[i];
for ( T1 k = m-1; 0 <= k; k-- )
{
T1 j2 = j / 2;
includ[k] = ( j != ( 2 * j2 ) );
j = j2;
}
/*
** Calculate the remaining elements of row I as explained
** in Bratley and Fox, section 2.
**
** Some tricky indexing here. Did I change it correctly?
*/
for ( j = m; j < maxcol; j++ )
{
T1 newv = v[i][j-m];
l = 1;
for ( T1 k = 0; k < m; k++ )
{
l = 2 * l;
if ( includ[k] )
{
newv = ( newv ^ ( l * v[i][j-k-1] ) );
}
}
v[i][j] = newv;
}
}
/*
** Multiply columns of V by appropriate power of 2.
*/
l = 1;
for ( T1 j = maxcol - 2; 0 <= j; j-- )
{
l = 2 * l;
for ( int i = 0; i < dim_num; i++ )
{
v[i][j] = v[i][j] * l;
}
}
/*
** RECIPD is 1/(common denominator of the elements in V).
*/
recipd = 1.0E+00 / ( ( T2 ) ( 2 * l ) );
}
if ( *seed < 0 )
*seed = 0 ;
if ( *seed == 0 )
{
l = 1;
for ( int i = 0; i < dim_num; i++ )
{
lastq[i] = 0;
}
}
else if ( *seed == seed_save + 1 )
{
l = bit_lo0 ( *seed );
}
else if ( *seed <= seed_save )
{
seed_save = 0;
l = 1;
for ( int i = 0; i < dim_num; i++ )
lastq[i] = 0;
for ( T1 seed_temp = seed_save; seed_temp <= (*seed)-1; seed_temp++ )
{
l = bit_lo0 ( seed_temp );
for ( int i = 0; i < dim_num; i++ )
{
lastq[i] = ( lastq[i] ^ v[i][l-1] );
}
}
l = bit_lo0 ( *seed );
}
else if ( seed_save+1 < *seed )
{
for ( T1 seed_temp = seed_save+1; seed_temp <= (*seed)-1; seed_temp++ )
{
l = bit_lo0 ( seed_temp );
for ( int i = 0; i < dim_num; i++ )
{
lastq[i] = ( lastq[i] ^ v[i][l-1] );
}
}
l = bit_lo0 ( *seed );
}
/*
** Check that the user is not calling too many times!
*/
if ( maxcol < l )
{
cout << "\n";
cout << "NEXT_SOBOL - Fatal error!\n";
cout << " The value of SEED seems to be too large!\n";
cout << " SEED = " << *seed << "\n";
cout << " MAXCOL = " << maxcol << "\n";
cout << " L = " << l << "\n";
exit ( 2 );
}
/*
** Calculate the new components of QUASI.
** The caret indicates the bitwise exclusive OR.
*/
for ( int i = 0; i < dim_num; i++ )
{
quasi[i] = ( ( T2 ) lastq[i] ) * recipd;
lastq[i] = ( lastq[i]^v[i][l-1] );
}
seed_save = *seed;
*seed = *seed + 1;
return;
}
template<typename T1, typename T2> T1 sobol_block( int dimension, int block_size, T1 seed, T2 *block )
{
for ( int iter = 0 ; iter < block_size ; iter++ )
{
next_sobol ( dimension, &seed, &block[iter*dimension] );
}
return seed;
}
template<typename T> void expand_unit_hypercube( int dimension, int block_size, T *block, T *lower_bound, T* upper_bound )
{
T *hypercube_length = new T[dimension];
for(int dim = 0 ; dim < dimension ; dim++ )
{
hypercube_length[dim] = upper_bound[dim]-lower_bound[dim] ;
}
int base = 0;
for ( int sim = 0 ; sim < block_size ; sim++ )
{
for (int dim = 0 ; dim < dimension ; dim++ )
{
block[base+dim] = lower_bound[dim] + hypercube_length[dim]*block[base+dim];
}
base+= dimension;
}
delete[] hypercube_length;
}
#undef DIM_MAX

View File

@ -140,14 +140,15 @@ enum BlockSimulationType
/*! Warning: do not to change existing values for 0 to 4: the values matter for homotopy_setup command */
enum SymbolType
{
eEndogenous = 0, //!< Endogenous
eExogenous = 1, //!< Exogenous
eExogenousDet = 2, //!< Exogenous deterministic
eParameter = 4, //!< Parameter
eModelLocalVariable = 10, //!< Local variable whose scope is model (pound expression)
eModFileLocalVariable = 11, //!< Local variable whose scope is mod file (model excluded)
eExternalFunction = 12, //!< External (user-defined) function
eTrend = 13 //!< Trend variable
eEndogenous = 0, //!< Endogenous
eExogenous = 1, //!< Exogenous
eExogenousDet = 2, //!< Exogenous deterministic
eParameter = 4, //!< Parameter
eModelLocalVariable = 10, //!< Local variable whose scope is model (pound expression)
eModFileLocalVariable = 11, //!< Local variable whose scope is mod file (model excluded)
eExternalFunction = 12, //!< External (user-defined) function
eTrend = 13, //!< Trend variable
eStatementDeclaredVariable = 14 //!< Local variable assigned within a Statement (see subsample statement for example)
};
enum ExpressionType
@ -233,6 +234,18 @@ enum external_function_type
ExternalFunctionSecondDerivative
};
enum PriorDistributions
{
eNoShape = 0,
eBeta = 1,
eGamma = 2,
eNormal = 3,
eInvGamma = 4,
eInvGamma1 = 4,
eUniform = 5,
eInvGamma2 = 6
};
struct Block_contain_type
{
int Equation, Variable, Own_Derivative;

View File

@ -27,6 +27,11 @@ using namespace std;
#include "ComputingTasks.hh"
#include "Statement.hh"
#include <boost/algorithm/string/trim.hpp>
#include <boost/algorithm/string/split.hpp>
#include <boost/lexical_cast.hpp>
#include <boost/tokenizer.hpp>
SteadyStatement::SteadyStatement(const OptionsList &options_list_arg) :
options_list(options_list_arg)
{
@ -135,15 +140,6 @@ StochSimulStatement::checkPass(ModFileStructure &mod_file_struct)
cerr << "ERROR: in 'stoch_simul', you cannot use option 'pruning' with 'k_order_solver' option or with 3rd order approximation" << endl;
exit(EXIT_FAILURE);
}
// Workaround for ticket #157
it = options_list.num_options.find("periods");
if (it != options_list.num_options.end() && atoi(it->second.c_str()) > 0
&& mod_file_struct.histval_present)
{
cerr << "ERROR: the 'periods' option of 'stoch_simul' is not compatible with a 'histval' block" << endl;
exit(EXIT_FAILURE);
}
}
void
@ -325,6 +321,30 @@ EstimationStatement::checkPass(ModFileStructure &mod_file_struct)
cerr << "ERROR: An estimation statement cannot take more than one dsge_var option." << endl;
exit(EXIT_FAILURE);
}
if (options_list.string_options.find("datafile") == options_list.string_options.end() &&
!mod_file_struct.estimation_data_statement_present)
{
cerr << "ERROR: The estimation statement requires a data file to be supplied "
<< "either from the data statement or from the deprecated option datafile." << endl;
exit(EXIT_FAILURE);
}
if (options_list.string_options.find("datafile") != options_list.string_options.end())
cerr << "WARNING: The datafile option of estimation has been deprecated. "
<< "Use the data command instead." << endl;
if (options_list.string_options.find("xls_sheet") != options_list.string_options.end())
cerr << "WARNING: The xls_sheet option of estimation has been deprecated. "
<< "Use the data command instead." << endl;
if (options_list.string_options.find("xls_range") != options_list.string_options.end())
cerr << "WARNING: The xls_range option of estimation has been deprecated. "
<< "Use the data command instead." << endl;
if (options_list.num_options.find("first_obs") != options_list.num_options.end())
cerr << "WARNING: The first_obs option of estimation has been deprecated. "
<< "Use the data command instead." << endl;
}
void
@ -426,7 +446,7 @@ EstimatedParamsStatement::checkPass(ModFileStructure &mod_file_struct)
mod_file_struct.dsge_prior_weight_in_estimated_params = true;
// Handle case of degenerate beta prior
if (it->prior == "1") //BETA_PDF is associated with "1" in DynareBison.yy
if (it->prior == eBeta)
try
{
if (it->mean->eval(eval_context_t()) == 0.5
@ -1330,41 +1350,195 @@ SvarIdentificationStatement::writeOutput(ostream &output, const string &basename
MarkovSwitchingStatement::MarkovSwitchingStatement(const OptionsList &options_list_arg) :
options_list(options_list_arg)
{
OptionsList::num_options_t::const_iterator it_num = options_list.num_options.find("ms.restrictions");
if (it_num != options_list.num_options.end())
{
using namespace boost;
OptionsList::num_options_t::const_iterator it_num_regimes =
options_list.num_options.find("ms.number_of_regimes");
if (it_num_regimes == options_list.num_options.end())
{
cerr << "ERROR: should not arrive here: MarkovSwitchingStatement constructor" << endl;
exit(EXIT_FAILURE);
}
int num_regimes = lexical_cast< int >(it_num_regimes->second);
vector<string> tokenizedRestrictions;
split(tokenizedRestrictions, it_num->second, is_any_of("["), token_compress_on);
for (vector<string>::iterator it = tokenizedRestrictions.begin();
it != tokenizedRestrictions.end(); it++ )
if (it->size() > 0)
{
vector<string> restriction;
split(restriction, *it, is_any_of("], "));
for (vector<string>::iterator it1 = restriction.begin();
it1 != restriction.end(); )
if (it1->empty())
restriction.erase(it1);
else
it1++;
if (restriction.size() != 3)
{
cerr << "ERROR: restrictions in the subsample statement must be specified in the form "
<< "[current_period_regime, next_period_regime, transition_probability]" << endl;
exit(EXIT_FAILURE);
}
try
{
int from_regime = lexical_cast< int >(restriction[0]);
int to_regime = lexical_cast< int >(restriction[1]);
if (from_regime > num_regimes || to_regime > num_regimes)
{
cerr << "ERROR: the regimes specified in the restrictions option must be "
<< "<= the number of regimes specified in the number_of_regimes option" << endl;
exit(EXIT_FAILURE);
}
if (restriction_map.find(make_pair(from_regime, to_regime)) !=
restriction_map.end())
{
cerr << "ERROR: two restrictions were given for: " << from_regime << ", "
<< to_regime << endl;
exit(EXIT_FAILURE);
}
double transition_probability = lexical_cast< double >(restriction[2]);
if (transition_probability > 1.0)
{
cerr << "ERROR: the transition probability, " << transition_probability
<< " must be less than 1" << endl;
exit(EXIT_FAILURE);
}
restriction_map[make_pair(from_regime, to_regime)] = transition_probability;
}
catch (const bad_lexical_cast &)
{
cerr << "ERROR: The first two arguments for a restriction must be integers "
<< "specifying the regime and the last must be a double specifying the "
<< "transition probability. You wrote [" << *it << endl;
exit(EXIT_FAILURE);
}
}
}
}
void
MarkovSwitchingStatement::checkPass(ModFileStructure &mod_file_struct)
{
OptionsList::num_options_t::const_iterator it_num = options_list.num_options.find("ms.restrictions");
if (it_num != options_list.num_options.end())
{
using namespace boost;
OptionsList::num_options_t::const_iterator it_num_regimes =
options_list.num_options.find("ms.number_of_regimes");
int num_regimes = lexical_cast< int >(it_num_regimes->second);
vector<double> col_trans_prob_sum (num_regimes, 0);
vector<double> row_trans_prob_sum (num_regimes, 0);
vector<bool> all_restrictions_in_row (num_regimes, true);
vector<bool> all_restrictions_in_col (num_regimes, true);
for (int row=0; row<num_regimes; row++)
for (int col=0; col<num_regimes; col++)
if (restriction_map.find(make_pair(row+1, col+1)) != restriction_map.end())
{
row_trans_prob_sum[row] += restriction_map[make_pair(row+1, col+1)];
col_trans_prob_sum[col] += restriction_map[make_pair(row+1, col+1)];
}
else
{
all_restrictions_in_row[row] = false;
all_restrictions_in_col[col] = false;
}
for (int i=0; i<num_regimes; i++)
{
if (all_restrictions_in_row[i])
{
if (row_trans_prob_sum[i] != 1.0)
{
cerr << "ERROR: When all transitions probabilities are specified for a certain "
<< "regime, they must sum to 1" << endl;
exit(EXIT_FAILURE);
}
}
else
if (row_trans_prob_sum[i] >= 1.0)
{
cerr << "ERROR: When transition probabilites are not specified for every regime, "
<< "their sum must be < 1" << endl;
exit(EXIT_FAILURE);
}
if (all_restrictions_in_col[i])
{
if (col_trans_prob_sum[i] != 1.0)
{
cerr << "ERROR: When all transitions probabilities are specified for a certain "
<< "regime, they must sum to 1" << endl;
exit(EXIT_FAILURE);
}
}
else
if (col_trans_prob_sum[i] >= 1.0)
{
cerr << "ERROR: When transition probabilites are not specified for every regime, "
<< "their sum must be < 1" << endl;
exit(EXIT_FAILURE);
}
}
}
}
void
MarkovSwitchingStatement::writeOutput(ostream &output, const string &basename) const
{
OptionsList::num_options_t::const_iterator itChain, itState, itNOS, itDuration;
bool isDurationAVec = true;
string infStr("Inf");
OptionsList::num_options_t::const_iterator itChain, itNOR, itDuration;
map<pair<int, int>, double >::const_iterator itR;
itChain = options_list.num_options.find("ms.chain");
if (itChain == options_list.num_options.end())
{
cerr << "MarkovSwitchingStatement::writeOutput() Should not arrive here (1). Please report this to the Dynare Team." << endl;
cerr << "MarkovSwitchingStatement::writeOutput() Should not arrive here (1). "
<< "Please report this to the Dynare Team." << endl;
exit(EXIT_FAILURE);
}
itDuration = options_list.num_options.find("ms.duration");
if (itDuration == options_list.num_options.end())
{
cerr << "MarkovSwitchingStatement::writeOutput() Should not arrive here (2). Please report this to the Dynare Team." << endl;
cerr << "MarkovSwitchingStatement::writeOutput() Should not arrive here (2). "
<< "Please report this to the Dynare Team." << endl;
exit(EXIT_FAILURE);
}
else if (atof(itDuration->second.c_str()) || infStr.compare(itDuration->second) == 0)
isDurationAVec = false;
output << "options_.ms.duration = " << itDuration->second << ";" << endl;
itNOR = options_list.num_options.find("ms.number_of_regimes");
if (itNOR != options_list.num_options.end())
for (int i = 0; i < atoi(itNOR->second.c_str()); i++)
{
output << "options_.ms.ms_chain(" << itChain->second << ").regime("
<< i+1 << ").duration = options_.ms.duration";
if (isDurationAVec)
output << "(" << i+1 << ")";
output << ";" << endl;
}
else
{
cerr << "MarkovSwitchingStatement::writeOutput() Should not arrive here (3). "
<< "Please report this to the Dynare Team." << endl;
exit(EXIT_FAILURE);
}
itState = options_list.num_options.find("ms.state");
itNOS = options_list.num_options.find("ms.number_of_states");
if (itState != options_list.num_options.end()
&& itNOS == options_list.num_options.end())
output << "options_.ms.ms_chain(" << itChain->second << ").state(" << itState->second << ").duration = " << itDuration->second << ";" << endl;
else if (itState == options_list.num_options.end()
&& itNOS != options_list.num_options.end())
for (int i = 0; i < atoi(itNOS->second.c_str()); i++)
output << "options_.ms.ms_chain(" << itChain->second << ").state(" << i+1 << ").duration = " << itDuration->second << ";" << endl;
else
{
cerr << "MarkovSwitchingStatement::writeOutput() Should not arrive here (3). Please report this to the Dynare Team." << endl;
exit(EXIT_FAILURE);
}
int restrictions_index = 0;
for (itR=restriction_map.begin(); itR != restriction_map.end(); itR++)
output << "options_.ms.ms_chain(" << itChain->second << ").restrictions("
<< ++restrictions_index << ") = {[" << itR->first.first << ", "
<< itR->first.second << ", " << itR->second << "]};" << endl;
}
SvarStatement::SvarStatement(const OptionsList &options_list_arg) :
@ -1431,3 +1605,429 @@ SvarStatement::writeOutput(ostream &output, const string &basename) const
else
output << "'ALL';" << endl;
}
SetTimeStatement::SetTimeStatement(const OptionsList &options_list_arg) :
options_list(options_list_arg)
{
}
void
SetTimeStatement::writeOutput(ostream &output, const string &basename) const
{
options_list.writeOutput(output);
}
EstimationDataStatement::EstimationDataStatement(const OptionsList &options_list_arg) :
options_list(options_list_arg)
{
}
void
EstimationDataStatement::checkPass(ModFileStructure &mod_file_struct)
{
mod_file_struct.estimation_data_statement_present = true;
OptionsList::num_options_t::const_iterator it = options_list.num_options.find("nobs");
if (it != options_list.num_options.end())
if (atoi(it->second.c_str()) <= 0)
{
cerr << "ERROR: The nobs option of the data statement only accepts positive integers." << endl;
exit(EXIT_FAILURE);
}
if (options_list.string_options.find("file") == options_list.string_options.end())
{
cerr << "ERROR: The file option must be passed to the data statement." << endl;
exit(EXIT_FAILURE);
}
}
void
EstimationDataStatement::writeOutput(ostream &output, const string &basename) const
{
options_list.writeOutput(output, "options_.dataset");
if (options_list.date_options.find("first_obs") == options_list.date_options.end())
output << "options_.dataset.firstobs = options_.initial_period;" << endl;
}
BasicPriorStatement::~BasicPriorStatement()
{
}
BasicPriorStatement::BasicPriorStatement(const string &name_arg,
const PriorDistributions &prior_shape_arg,
const expr_t &variance_arg,
const OptionsList &options_list_arg) :
name(name_arg),
prior_shape(prior_shape_arg),
variance(variance_arg),
options_list(options_list_arg),
first_statement_encountered(false)
{
}
void
BasicPriorStatement::checkPass(ModFileStructure &mod_file_struct)
{
if (prior_shape == eNoShape)
{
cerr << "ERROR: You must pass the shape option to the prior statement." << endl;
exit(EXIT_FAILURE);
}
if (options_list.num_options.find("date1") != options_list.num_options.end() ||
options_list.num_options.find("date2") != options_list.num_options.end())
if (options_list.num_options.find("date1") == options_list.num_options.end() ||
options_list.num_options.find("date2") == options_list.num_options.end())
{
cerr << "ERROR: PriorStatement::checkPass(1). Should not arrive here. "
<< "Please inform Dynare Team." << endl;
exit(EXIT_FAILURE);
}
}
void
BasicPriorStatement::get_base_name(const SymbolType symb_type, string &lhs_field) const
{
if (symb_type == eExogenous || symb_type == eExogenousDet)
lhs_field = "structural_innovation";
else
lhs_field = "measurement_error";
}
void
BasicPriorStatement::writePriorIndex(ostream &output, const string &lhs_field) const
{
if (first_statement_encountered)
output << "prior_indx = 1;" << endl;
else
output << "prior_indx = size(estimation_info" << lhs_field << "_index, 2) + 1;" << endl;
}
void
BasicPriorStatement::writeVarianceOption(ostream &output, const string &lhs_field) const
{
if (variance)
{
output << "estimation_info" << lhs_field << "(prior_indx).variance = ";
variance->writeOutput(output);
output << ";" << endl;
}
}
void
BasicPriorStatement::writeOutputHelper(ostream &output, const string &field, const string &lhs_field) const
{
OptionsList::num_options_t::const_iterator itn = options_list.num_options.find(field);
if (itn != options_list.num_options.end())
output << "estimation_info" << lhs_field << "(prior_indx)." << field
<< " = " << itn->second << ";" << endl;
OptionsList::date_options_t::const_iterator itd = options_list.date_options.find(field);
if (itd != options_list.date_options.end())
output << "estimation_info" << lhs_field << "(prior_indx)." << field
<< " = '" << itd->second << "';" << endl;
}
void
BasicPriorStatement::writeShape(ostream &output, const string &lhs_field) const
{
assert(prior_shape != eNoShape);
output << "estimation_info" << lhs_field << "(prior_indx).shape = " << prior_shape << ";" << endl;
}
PriorStatement::PriorStatement(const string &name_arg,
const PriorDistributions &prior_shape_arg,
const expr_t &variance_arg,
const OptionsList &options_list_arg) :
BasicPriorStatement(name_arg, prior_shape_arg, variance_arg, options_list_arg)
{
}
void
PriorStatement::checkPass(ModFileStructure &mod_file_struct)
{
BasicPriorStatement::checkPass(mod_file_struct);
if (!mod_file_struct.prior_statement_present)
first_statement_encountered = true;
mod_file_struct.prior_statement_present = true;
}
void
PriorStatement::writeOutput(ostream &output, const string &basename) const
{
string lhs_field = ".prior";
writePriorIndex(output, lhs_field);
output << "estimation_info" << lhs_field << "_index(prior_indx) = {'" << name << "'};" << endl
<< "estimation_info" << lhs_field <<"(prior_indx).name = '" << name << "';" << endl;
writeShape(output, lhs_field);
writeOutputHelper(output, "mean", lhs_field);
writeOutputHelper(output, "mode", lhs_field);
writeOutputHelper(output, "stdev", lhs_field);
writeOutputHelper(output, "shape", lhs_field);
writeOutputHelper(output, "shift", lhs_field);
writeOutputHelper(output, "date1", lhs_field);
writeOutputHelper(output, "date2", lhs_field);
writeOutputHelper(output, "domain", lhs_field);
writeOutputHelper(output, "interval", lhs_field);
writeVarianceOption(output, lhs_field);
}
StdPriorStatement::StdPriorStatement(const string &name_arg,
const PriorDistributions &prior_shape_arg,
const expr_t &variance_arg,
const OptionsList &options_list_arg,
const SymbolTable &symbol_table_arg ) :
BasicPriorStatement(name_arg, prior_shape_arg, variance_arg, options_list_arg),
symbol_table(symbol_table_arg)
{
}
void
StdPriorStatement::checkPass(ModFileStructure &mod_file_struct)
{
BasicPriorStatement::checkPass(mod_file_struct);
if (!mod_file_struct.std_prior_statement_present)
first_statement_encountered = true;
mod_file_struct.std_prior_statement_present = true;
}
void
StdPriorStatement::writeOutput(ostream &output, const string &basename) const
{
string lhs_field;
get_base_name(symbol_table.getType(name), lhs_field);
lhs_field = "." + lhs_field + ".prior";
writePriorIndex(output, lhs_field);
output << "estimation_info" << lhs_field << "_index(prior_indx) = {'" << name << "'};" << endl;
output << "estimation_info" << lhs_field << "(prior_indx).name = '" << name << "';" << endl;
writeShape(output, lhs_field);
writeOutputHelper(output, "mean", lhs_field);
writeOutputHelper(output, "mode", lhs_field);
writeOutputHelper(output, "stdev", lhs_field);
writeOutputHelper(output, "shape", lhs_field);
writeOutputHelper(output, "shift", lhs_field);
writeOutputHelper(output, "domain", lhs_field);
writeOutputHelper(output, "interval", lhs_field);
writeVarianceOption(output, lhs_field);
}
CorrPriorStatement::CorrPriorStatement(const string &name_arg1, const string &name_arg2,
const PriorDistributions &prior_shape_arg,
const expr_t &variance_arg,
const OptionsList &options_list_arg,
const SymbolTable &symbol_table_arg ) :
BasicPriorStatement(name_arg1, prior_shape_arg, variance_arg, options_list_arg),
name1(name_arg2),
symbol_table(symbol_table_arg)
{
}
void
CorrPriorStatement::checkPass(ModFileStructure &mod_file_struct)
{
BasicPriorStatement::checkPass(mod_file_struct);
if (symbol_table.getType(name) != symbol_table.getType(name1))
{
cerr << "ERROR: In the corr(A,B).prior statement, A and B must be of the same type. "
<< "In your case, " << name << " and " << name1 << " are of different "
<< "types." << endl;
exit(EXIT_FAILURE);
}
if (!mod_file_struct.corr_prior_statement_present)
first_statement_encountered = true;
mod_file_struct.corr_prior_statement_present = true;
}
void
CorrPriorStatement::writeOutput(ostream &output, const string &basename) const
{
string lhs_field;
get_base_name(symbol_table.getType(name), lhs_field);
lhs_field = "." + lhs_field + "_corr.prior";
writePriorIndex(output, lhs_field);
output << "estimation_info" << lhs_field << "_index(prior_indx) = {'" << name << "_" << name1 << "'};" << endl;
output << "estimation_info" << lhs_field << "(prior_indx).name1 = '" << name << "';" << endl;
output << "estimation_info" << lhs_field << "(prior_indx).name2 = '" << name1 << "';" << endl;
writeShape(output, lhs_field);
writeOutputHelper(output, "mean", lhs_field);
writeOutputHelper(output, "mode", lhs_field);
writeOutputHelper(output, "stdev", lhs_field);
writeOutputHelper(output, "shape", lhs_field);
writeOutputHelper(output, "shift", lhs_field);
writeOutputHelper(output, "domain", lhs_field);
writeOutputHelper(output, "interval", lhs_field);
writeVarianceOption(output, lhs_field);
}
BasicOptionsStatement::~BasicOptionsStatement()
{
}
BasicOptionsStatement::BasicOptionsStatement(const string &name_arg,
const OptionsList &options_list_arg) :
name(name_arg),
options_list(options_list_arg),
first_statement_encountered(false)
{
}
void
BasicOptionsStatement::checkPass(ModFileStructure &mod_file_struct)
{
if (options_list.num_options.find("date1") != options_list.num_options.end() ||
options_list.num_options.find("date2") != options_list.num_options.end())
if (options_list.num_options.find("date1") == options_list.num_options.end() ||
options_list.num_options.find("date2") == options_list.num_options.end())
{
cerr << "ERROR: OptionsStatement::checkPass(1). Should not arrive here. "
<< "Please inform Dynare Team." << endl;
exit(EXIT_FAILURE);
}
}
void
BasicOptionsStatement::writeOptionsIndex(ostream &output, const string &lhs_field) const
{
if (first_statement_encountered)
output << "options_indx = 1;" << endl;
else
output << "options_indx = size(estimation_info" << lhs_field << "_index, 2) + 1;" << endl;
}
void
BasicOptionsStatement::get_base_name(const SymbolType symb_type, string &lhs_field) const
{
if (symb_type == eExogenous || symb_type == eExogenousDet)
lhs_field = "structural_innovation";
else
lhs_field = "measurement_error";
}
void
BasicOptionsStatement::writeOutputHelper(ostream &output, const string &field, const string &lhs_field) const
{
OptionsList::num_options_t::const_iterator itn = options_list.num_options.find(field);
if (itn != options_list.num_options.end())
output << "estimation_info" << lhs_field << "(options_indx)." << field
<< " = " << itn->second << ";" << endl;
OptionsList::date_options_t::const_iterator itd = options_list.date_options.find(field);
if (itd != options_list.date_options.end())
output << "estimation_info" << lhs_field << "(options_indx)." << field
<< " = '" << itd->second << "';" << endl;
}
OptionsStatement::OptionsStatement(const string &name_arg,
const OptionsList &options_list_arg) :
BasicOptionsStatement(name_arg, options_list_arg)
{
}
void
OptionsStatement::checkPass(ModFileStructure &mod_file_struct)
{
BasicOptionsStatement::checkPass(mod_file_struct);
if (!mod_file_struct.options_statement_present)
first_statement_encountered = true;
mod_file_struct.options_statement_present = true;
}
void
OptionsStatement::writeOutput(ostream &output, const string &basename) const
{
string lhs_field = ".options";
writeOptionsIndex(output, lhs_field);
output << "estimation_info" << lhs_field <<"_index(options_indx) = {'" << name << "'};" << endl
<< "estimation_info" << lhs_field << "(options_indx).name = '" << name << "';" << endl;
writeOutputHelper(output, "init", lhs_field);
writeOutputHelper(output, "bounds", lhs_field);
writeOutputHelper(output, "jscale", lhs_field);
writeOutputHelper(output, "date1", lhs_field);
writeOutputHelper(output, "date2", lhs_field);
}
StdOptionsStatement::StdOptionsStatement(const string &name_arg,
const OptionsList &options_list_arg,
const SymbolTable &symbol_table_arg ) :
BasicOptionsStatement(name_arg, options_list_arg),
symbol_table(symbol_table_arg)
{
}
void
StdOptionsStatement::checkPass(ModFileStructure &mod_file_struct)
{
BasicOptionsStatement::checkPass(mod_file_struct);
if (!mod_file_struct.std_options_statement_present)
first_statement_encountered = true;
mod_file_struct.std_options_statement_present = true;
}
void
StdOptionsStatement::writeOutput(ostream &output, const string &basename) const
{
string lhs_field;
get_base_name(symbol_table.getType(name), lhs_field);
lhs_field = "." + lhs_field + ".options";
writeOptionsIndex(output, lhs_field);
output << "estimation_info" << lhs_field << "_index(options_indx) = {'" << name << "'};" << endl;
output << "estimation_info" << lhs_field << "(options_indx).name = '" << name << "';" << endl;
writeOutputHelper(output, "init", lhs_field);
writeOutputHelper(output, "bounds", lhs_field);
writeOutputHelper(output, "jscale", lhs_field);
writeOutputHelper(output, "date1", lhs_field);
writeOutputHelper(output, "date2", lhs_field);
}
CorrOptionsStatement::CorrOptionsStatement(const string &name_arg1, const string &name_arg2,
const OptionsList &options_list_arg,
const SymbolTable &symbol_table_arg ) :
BasicOptionsStatement(name_arg1, options_list_arg),
name1(name_arg2),
symbol_table(symbol_table_arg)
{
}
void
CorrOptionsStatement::checkPass(ModFileStructure &mod_file_struct)
{
if (symbol_table.getType(name) != symbol_table.getType(name1))
{
cerr << "ERROR: In the corr(A,B).options statement, A and B must be of the same type. "
<< "In your case, " << name << " and " << name1 << " are of different "
<< "types." << endl;
exit(EXIT_FAILURE);
}
if (!mod_file_struct.corr_options_statement_present)
first_statement_encountered = true;
mod_file_struct.corr_prior_statement_present = true;
}
void
CorrOptionsStatement::writeOutput(ostream &output, const string &basename) const
{
string lhs_field;
get_base_name(symbol_table.getType(name), lhs_field);
lhs_field = "." + lhs_field + "_corr.options";
writeOptionsIndex(output, lhs_field);
output << "estimation_info" << lhs_field << "_index(options_indx) = {'" << name << "_" << name1 << "'};" << endl;
lhs_field += ".";
output << "estimation_info" << lhs_field << "(options_indx).name1 = '" << name << "';" << endl;
output << "estimation_info" << lhs_field << "(options_indx).name2 = '" << name1 << "';" << endl;
writeOutputHelper(output, "init", lhs_field);
writeOutputHelper(output, "bounds", lhs_field);
writeOutputHelper(output, "jscale", lhs_field);
writeOutputHelper(output, "date1", lhs_field);
writeOutputHelper(output, "date2", lhs_field);
}

View File

@ -251,7 +251,8 @@ class EstimationParams
{
public:
int type;
string name, name2, prior;
string name, name2;
PriorDistributions prior;
expr_t init_val, low_bound, up_bound, mean, std, p3, p4, jscale;
void
@ -260,7 +261,7 @@ public:
type = 0;
name = "";
name2 = "";
prior = "NaN";
prior = eNoShape;
init_val = datatree.NaN;
low_bound = datatree.MinusInfinity;
up_bound = datatree.Infinity;
@ -542,8 +543,10 @@ class MarkovSwitchingStatement : public Statement
{
private:
const OptionsList options_list;
map <pair<int, int >, double > restriction_map;
public:
MarkovSwitchingStatement(const OptionsList &options_list_arg);
virtual void checkPass(ModFileStructure &mod_file_struct);
virtual void writeOutput(ostream &output, const string &basename) const;
};
@ -556,4 +559,133 @@ public:
virtual void writeOutput(ostream &output, const string &basename) const;
};
class SetTimeStatement : public Statement
{
private:
const OptionsList options_list;
public:
SetTimeStatement(const OptionsList &options_list_arg);
virtual void writeOutput(ostream &output, const string &basename) const;
};
class EstimationDataStatement : public Statement
{
private:
const OptionsList options_list;
public:
EstimationDataStatement(const OptionsList &options_list_arg);
virtual void checkPass(ModFileStructure &mod_file_struct);
virtual void writeOutput(ostream &output, const string &basename) const;
};
class BasicPriorStatement : public Statement
{
public:
virtual ~BasicPriorStatement();
protected:
const string name;
const PriorDistributions prior_shape;
const expr_t variance;
const OptionsList options_list;
bool first_statement_encountered;
BasicPriorStatement(const string &name_arg,
const PriorDistributions &prior_shape_arg,
const expr_t &variance_arg,
const OptionsList &options_list_arg);
virtual void checkPass(ModFileStructure &mod_file_struct);
void get_base_name(const SymbolType symb_type, string &lhs_field) const;
void writePriorIndex(ostream &output, const string &lhs_field) const;
void writeVarianceOption(ostream &output, const string &lhs_field) const;
void writeOutputHelper(ostream &output, const string &field, const string &lhs_field) const;
void writeShape(ostream &output, const string &lhs_field) const;
};
class PriorStatement : public BasicPriorStatement
{
public:
PriorStatement(const string &name_arg,
const PriorDistributions &prior_shape_arg,
const expr_t &variance_arg,
const OptionsList &options_list_arg);
virtual void checkPass(ModFileStructure &mod_file_struct);
virtual void writeOutput(ostream &output, const string &basename) const;
};
class StdPriorStatement : public BasicPriorStatement
{
private:
const SymbolTable symbol_table;
public:
StdPriorStatement(const string &name_arg,
const PriorDistributions &prior_shape_arg,
const expr_t &variance_arg,
const OptionsList &options_list_arg,
const SymbolTable &symbol_table_arg);
virtual void checkPass(ModFileStructure &mod_file_struct);
virtual void writeOutput(ostream &output, const string &basename) const;
};
class CorrPriorStatement : public BasicPriorStatement
{
private:
const string name1;
const SymbolTable symbol_table;
public:
CorrPriorStatement(const string &name_arg1,
const string &name_arg2,
const PriorDistributions &prior_shape_arg,
const expr_t &variance_arg,
const OptionsList &options_list_arg,
const SymbolTable &symbol_table_arg);
virtual void checkPass(ModFileStructure &mod_file_struct);
virtual void writeOutput(ostream &output, const string &basename) const;
};
class BasicOptionsStatement : public Statement
{
public:
virtual ~BasicOptionsStatement();
protected:
const string name;
const OptionsList options_list;
bool first_statement_encountered;
BasicOptionsStatement(const string &name_arg,
const OptionsList &options_list_arg);
void get_base_name(const SymbolType symb_type, string &lhs_field) const;
virtual void checkPass(ModFileStructure &mod_file_struct);
void writeOptionsIndex(ostream &output, const string &lhs_field) const;
void writeOutputHelper(ostream &output, const string &field, const string &lhs_field) const;
};
class OptionsStatement : public BasicOptionsStatement
{
public:
OptionsStatement(const string &name_arg, const OptionsList &options_list_arg);
virtual void checkPass(ModFileStructure &mod_file_struct);
virtual void writeOutput(ostream &output, const string &basename) const;
};
class StdOptionsStatement : public BasicOptionsStatement
{
private:
const SymbolTable symbol_table;
public:
StdOptionsStatement(const string &name_arg, const OptionsList &options_list_arg,
const SymbolTable &symbol_table_arg);
virtual void checkPass(ModFileStructure &mod_file_struct);
virtual void writeOutput(ostream &output, const string &basename) const;
};
class CorrOptionsStatement : public BasicOptionsStatement
{
private:
const string name1;
const SymbolTable symbol_table;
public:
CorrOptionsStatement(const string &name_arg1, const string &name_arg2,
const OptionsList &options_list_arg, const SymbolTable &symbol_table_arg);
virtual void checkPass(ModFileStructure &mod_file_struct);
virtual void writeOutput(ostream &output, const string &basename) const;
};
#endif

View File

@ -1518,7 +1518,8 @@ void
DynamicModel::writeDynamicCFile(const string &dynamic_basename, const int order) const
{
string filename = dynamic_basename + ".c";
ofstream mDynamicModelFile;
string filename_mex = dynamic_basename + "_mex.c";
ofstream mDynamicModelFile, mDynamicMexFile;
mDynamicModelFile.open(filename.c_str(), ios::out | ios::binary);
if (!mDynamicModelFile.is_open())
@ -1531,12 +1532,16 @@ DynamicModel::writeDynamicCFile(const string &dynamic_basename, const int order)
<< " *" << endl
<< " * Warning : this file is generated automatically by Dynare" << endl
<< " * from model file (.mod)" << endl
<< endl
<< " */" << endl
<< "#include <math.h>" << endl
<< "#include \"mex.h\"" << endl
<< endl
<< "#define max(a, b) (((a) > (b)) ? (a) : (b))" << endl
<< "#include <math.h>" << endl;
if (external_functions_table.get_total_number_of_unique_model_block_external_functions())
// External Matlab function, implies Dynamic function will call mex
mDynamicModelFile << "#include \"mex.h\"" << endl;
else
mDynamicModelFile << "#include <stdlib.h>" << endl;
mDynamicModelFile << "#define max(a, b) (((a) > (b)) ? (a) : (b))" << endl
<< "#define min(a, b) (((a) > (b)) ? (b) : (a))" << endl;
// Write function definition if oPowerDeriv is used
@ -1545,76 +1550,93 @@ DynamicModel::writeDynamicCFile(const string &dynamic_basename, const int order)
// Writing the function body
writeDynamicModel(mDynamicModelFile, true);
// Writing the gateway routine
mDynamicModelFile << "/* The gateway routine */" << endl
<< "void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])" << endl
<< "{" << endl
<< " double *y, *x, *params, *steady_state;" << endl
<< " double *residual, *g1, *v2, *v3;" << endl
<< " int nb_row_x, it_;" << endl
<< endl
<< " /* Check that no derivatives of higher order than computed are being requested */ " << endl
<< " if (nlhs > " << order + 1 << ") " << endl
<< " mexErrMsgTxt(\"Derivatives of higher order than computed have been requested\"); " << endl
<< " /* Create a pointer to the input matrix y. */" << endl
<< " y = mxGetPr(prhs[0]);" << endl
<< endl
<< " /* Create a pointer to the input matrix x. */" << endl
<< " x = mxGetPr(prhs[1]);" << endl
<< endl
<< " /* Create a pointer to the input matrix params. */" << endl
<< " params = mxGetPr(prhs[2]);" << endl
<< endl
<< " /* Create a pointer to the input matrix steady_state. */" << endl
<< " steady_state = mxGetPr(prhs[3]);" << endl
<< endl
<< " /* Fetch time index */" << endl
<< " it_ = (int) mxGetScalar(prhs[4]) - 1;" << endl
<< endl
<< " /* Gets number of rows of matrix x. */" << endl
<< " nb_row_x = mxGetM(prhs[1]);" << endl
<< endl
<< " residual = NULL;" << endl
<< " if (nlhs >= 1)" << endl
<< " {" << endl
<< " /* Set the output pointer to the output matrix residual. */" << endl
<< " plhs[0] = mxCreateDoubleMatrix(" << equations.size() << ",1, mxREAL);" << endl
<< " /* Create a C pointer to a copy of the output matrix residual. */" << endl
<< " residual = mxGetPr(plhs[0]);" << endl
<< " }" << endl
<< endl
<< " g1 = NULL;" << endl
<< " if (nlhs >= 2)" << endl
<< " {" << endl
<< " /* Set the output pointer to the output matrix g1. */" << endl
<< " plhs[1] = mxCreateDoubleMatrix(" << equations.size() << ", " << dynJacobianColsNbr << ", mxREAL);" << endl
<< " /* Create a C pointer to a copy of the output matrix g1. */" << endl
<< " g1 = mxGetPr(plhs[1]);" << endl
<< " }" << endl
<< endl
<< " v2 = NULL;" << endl
<< " if (nlhs >= 3)" << endl
<< " {" << endl
<< " /* Set the output pointer to the output matrix v2. */" << endl
<< " plhs[2] = mxCreateDoubleMatrix(" << NNZDerivatives[1] << ", " << 3
<< ", mxREAL);" << endl
<< " v2 = mxGetPr(plhs[2]);" << endl
<< " }" << endl
<< endl
<< " v3 = NULL;" << endl
<< " if (nlhs >= 4)" << endl
<< " {" << endl
<< " /* Set the output pointer to the output matrix v3. */" << endl
<< " plhs[3] = mxCreateDoubleMatrix(" << NNZDerivatives[2] << ", " << 3 << ", mxREAL);" << endl
<< " v3 = mxGetPr(plhs[3]);" << endl
<< " }" << endl
<< endl
<< " /* Call the C subroutines. */" << endl
<< " Dynamic(y, x, nb_row_x, params, steady_state, it_, residual, g1, v2, v3);" << endl
<< "}" << endl;
writePowerDeriv(mDynamicModelFile, true);
mDynamicModelFile.close();
mDynamicMexFile.open(filename_mex.c_str(), ios::out | ios::binary);
if (!mDynamicMexFile.is_open())
{
cerr << "Error: Can't open file " << filename_mex << " for writing" << endl;
exit(EXIT_FAILURE);
}
// Writing the gateway routine
mDynamicMexFile << "/*" << endl
<< " * " << filename_mex << " : The gateway routine used to call the Dynamic function "
<< "located in " << filename << endl
<< " *" << endl
<< " * Warning : this file is generated automatically by Dynare" << endl
<< " * from model file (.mod)" << endl
<< endl
<< " */" << endl << endl
<< "#include \"mex.h\"" << endl << endl
<< "void Dynamic(double *y, double *x, int nb_row_x, double *params, double *steady_state, int it_, double *residual, double *g1, double *v2, double *v3);" << endl
<< "void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])" << endl
<< "{" << endl
<< " double *y, *x, *params, *steady_state;" << endl
<< " double *residual, *g1, *v2, *v3;" << endl
<< " int nb_row_x, it_;" << endl
<< endl
<< " /* Check that no derivatives of higher order than computed are being requested */" << endl
<< " if (nlhs > " << order + 1 << ")" << endl
<< " mexErrMsgTxt(\"Derivatives of higher order than computed have been requested\");" << endl
<< " /* Create a pointer to the input matrix y. */" << endl
<< " y = mxGetPr(prhs[0]);" << endl
<< endl
<< " /* Create a pointer to the input matrix x. */" << endl
<< " x = mxGetPr(prhs[1]);" << endl
<< endl
<< " /* Create a pointer to the input matrix params. */" << endl
<< " params = mxGetPr(prhs[2]);" << endl
<< endl
<< " /* Create a pointer to the input matrix steady_state. */" << endl
<< " steady_state = mxGetPr(prhs[3]);" << endl
<< endl
<< " /* Fetch time index */" << endl
<< " it_ = (int) mxGetScalar(prhs[4]) - 1;" << endl
<< endl
<< " /* Gets number of rows of matrix x. */" << endl
<< " nb_row_x = mxGetM(prhs[1]);" << endl
<< endl
<< " residual = NULL;" << endl
<< " if (nlhs >= 1)" << endl
<< " {" << endl
<< " /* Set the output pointer to the output matrix residual. */" << endl
<< " plhs[0] = mxCreateDoubleMatrix(" << equations.size() << ",1, mxREAL);" << endl
<< " /* Create a C pointer to a copy of the output matrix residual. */" << endl
<< " residual = mxGetPr(plhs[0]);" << endl
<< " }" << endl
<< endl
<< " g1 = NULL;" << endl
<< " if (nlhs >= 2)" << endl
<< " {" << endl
<< " /* Set the output pointer to the output matrix g1. */" << endl
<< " plhs[1] = mxCreateDoubleMatrix(" << equations.size() << ", " << dynJacobianColsNbr << ", mxREAL);" << endl
<< " /* Create a C pointer to a copy of the output matrix g1. */" << endl
<< " g1 = mxGetPr(plhs[1]);" << endl
<< " }" << endl
<< endl
<< " v2 = NULL;" << endl
<< " if (nlhs >= 3)" << endl
<< " {" << endl
<< " /* Set the output pointer to the output matrix v2. */" << endl
<< " plhs[2] = mxCreateDoubleMatrix(" << NNZDerivatives[1] << ", " << 3
<< ", mxREAL);" << endl
<< " v2 = mxGetPr(plhs[2]);" << endl
<< " }" << endl
<< endl
<< " v3 = NULL;" << endl
<< " if (nlhs >= 4)" << endl
<< " {" << endl
<< " /* Set the output pointer to the output matrix v3. */" << endl
<< " plhs[3] = mxCreateDoubleMatrix(" << NNZDerivatives[2] << ", " << 3 << ", mxREAL);" << endl
<< " v3 = mxGetPr(plhs[3]);" << endl
<< " }" << endl
<< endl
<< " /* Call the C subroutines. */" << endl
<< " Dynamic(y, x, nb_row_x, params, steady_state, it_, residual, g1, v2, v3);" << endl
<< "}" << endl;
mDynamicMexFile.close();
}
string

View File

@ -33,6 +33,7 @@ using namespace std;
class ParsingDriver;
#include "ExprNode.hh"
#include "CodeInterpreter.hh"
/* Little hack: we redefine the macro which computes the locations, because
we need to access the location from within the parsing driver for error
@ -75,6 +76,7 @@ class ParsingDriver;
SymbolType symbol_type_val;
vector<string *> *vector_string_val;
vector<int> *vector_int_val;
PriorDistributions prior_distributions_val;
};
%{
@ -94,17 +96,18 @@ class ParsingDriver;
%token BVAR_PRIOR_MU BVAR_PRIOR_OMEGA BVAR_PRIOR_TAU BVAR_PRIOR_TRAIN
%token BVAR_REPLIC BYTECODE
%token CHANGE_TYPE CHECK CONDITIONAL_FORECAST CONDITIONAL_FORECAST_PATHS CONF_SIG CONSTANT CONTROLLED_VAREXO CORR COVAR CUTOFF
%token DATAFILE DR_ALGO DROP DSAMPLE DYNASAVE DYNATYPE CALIBRATION
%token DATAFILE FILE DR_ALGO DROP DSAMPLE DYNASAVE DYNATYPE CALIBRATION
%token END ENDVAL EQUAL ESTIMATION ESTIMATED_PARAMS ESTIMATED_PARAMS_BOUNDS ESTIMATED_PARAMS_INIT
%token FILENAME FILTER_STEP_AHEAD FILTERED_VARS FIRST_OBS
%token FILENAME FILTER_STEP_AHEAD FILTERED_VARS FIRST_OBS LAST_OBS SET_TIME
%token <string_val> FLOAT_NUMBER
%token FORECAST K_ORDER_SOLVER INSTRUMENTS
%token GAMMA_PDF GRAPH CONDITIONAL_VARIANCE_DECOMPOSITION NOCHECK
%token FORECAST K_ORDER_SOLVER INSTRUMENTS PRIOR SHIFT MEAN STDEV VARIANCE MODE INTERVAL SHAPE DOMAINN
%token GAMMA_PDF GRAPH CONDITIONAL_VARIANCE_DECOMPOSITION NOCHECK STD
%token HISTVAL HOMOTOPY_SETUP HOMOTOPY_MODE HOMOTOPY_STEPS HP_FILTER HP_NGRID
%token IDENTIFICATION INF_CONSTANT INITVAL INITVAL_FILE
%token IDENTIFICATION INF_CONSTANT INITVAL INITVAL_FILE BOUNDS JSCALE INIT
%token <string_val> INT_NUMBER
%token <string_val> DATE_NUMBER
%token INV_GAMMA_PDF INV_GAMMA1_PDF INV_GAMMA2_PDF IRF IRF_SHOCKS
%token KALMAN_ALGO KALMAN_TOL
%token KALMAN_ALGO KALMAN_TOL SUBSAMPLES OPTIONS
%token LABELS LAPLACE LIK_ALGO LIK_INIT LINEAR LOAD_IDENT_FILES LOAD_MH_FILE LOAD_PARAMS_AND_STEADY_STATE LOGLINEAR
%token MARKOWITZ MARGINAL_DENSITY MAX MAXIT
%token MFS MH_DROP MH_INIT_SCALE MH_JSCALE MH_MODE MH_NBLOCKS MH_REPLIC MH_RECOVER MIN MINIMAL_SOLVING_PERIODS
@ -144,18 +147,18 @@ class ParsingDriver;
/* end of GSA analysis*/
%token FREQ INITIAL_YEAR INITIAL_SUBPERIOD FINAL_YEAR FINAL_SUBPERIOD DATA VLIST LOG_VAR PERCENT_VAR
%token VLISTLOG VLISTPER
%token RESTRICTION RESTRICTIONS RESTRICTION_FNAME CROSS_RESTRICTIONS NLAGS CONTEMP_REDUCED_FORM REAL_PSEUDO_FORECAST
%token RESTRICTION RESTRICTION_FNAME CROSS_RESTRICTIONS NLAGS CONTEMP_REDUCED_FORM REAL_PSEUDO_FORECAST
%token NONE DUMMY_OBS NSTATES INDXSCALESSTATES NO_BAYESIAN_PRIOR SPECIFICATION SIMS_ZHA
%token <string_val> ALPHA BETA ABAND NINV CMS NCMS CNUM
%token <string_val> ALPHA BETA ABAND NINV CMS NCMS CNUM GAMMA INV_GAMMA INV_GAMMA1 INV_GAMMA2 NORMAL UNIFORM
%token GSIG2_LMDM Q_DIAG FLAT_PRIOR NCSK NSTD
%token INDXPARR INDXOVR INDXAP APBAND INDXIMF IMFBAND INDXFORE FOREBAND INDXGFOREHAT INDXGIMFHAT
%token INDXESTIMA INDXGDLS EQ_MS FILTER_COVARIANCE FILTER_DECOMPOSITION
%token EQ_CMS TLINDX TLNUMBER BANACT
%token EQ_CMS TLINDX TLNUMBER BANACT RESTRICTIONS
%token OUTPUT_FILE_TAG DRAWS_NBR_BURN_IN_1 DRAWS_NBR_BURN_IN_2 HORIZON
%token SBVAR TREND_VAR DEFLATOR GROWTH_FACTOR MS_IRF MS_VARIANCE_DECOMPOSITION
%token MS_ESTIMATION MS_SIMULATION MS_COMPUTE_MDD MS_COMPUTE_PROBABILITIES MS_FORECAST
%token SVAR_IDENTIFICATION EQUATION EXCLUSION LAG UPPER_CHOLESKY LOWER_CHOLESKY MONTHLY QUARTERLY
%token MARKOV_SWITCHING CHAIN STATE DURATION NUMBER_OF_STATES
%token MARKOV_SWITCHING CHAIN DURATION NUMBER_OF_REGIMES
%token SVAR COEFF COEFFICIENTS VARIANCES CONSTANTS EQUATIONS
%token EXTERNAL_FUNCTION EXT_FUNC_NAME EXT_FUNC_NARGS FIRST_DERIV_PROVIDED SECOND_DERIV_PROVIDED
%token SELECTED_VARIABLES_ONLY COVA_COMPUTE SIMULATION_FILE_TAG FILE_TAG
@ -173,14 +176,14 @@ class ParsingDriver;
%type <node_val> expression expression_or_empty
%type <node_val> equation hand_side
%type <string_val> non_negative_number signed_number signed_integer
%type <string_val> filename symbol
%type <string_val> vec_value_1 vec_value
%type <string_val> range prior
%type <string_val> non_negative_number signed_number signed_integer date_number
%type <string_val> filename symbol vec_of_vec_value vec_value_list
%type <string_val> vec_value_1 vec_value signed_inf signed_number_w_inf
%type <string_val> range vec_value_w_inf vec_value_1_w_inf
%type <symbol_type_val> change_type_arg
%type <vector_string_val> change_type_var_list
%type <vector_int_val> vec_int_elem vec_int_1 vec_int vec_int_number
%type <prior_distributions_val> prior_pdf prior_distribution
%%
%start statement_list;
@ -213,6 +216,12 @@ statement : parameters
| estimated_params
| estimated_params_bounds
| estimated_params_init
| set_time
| data
| prior
| subsamples
| subsamples_eq
| options
| varobs
| observation_trends
| unit_root_vars
@ -733,9 +742,9 @@ ms_options_list : ms_options_list COMMA ms_options
;
ms_options : o_chain
| o_state
| o_duration
| o_number_of_states
| o_restrictions
| o_number_of_regimes
;
svar : SVAR '(' svar_options_list ')' ';'
@ -963,6 +972,10 @@ non_negative_number : INT_NUMBER
| FLOAT_NUMBER
;
date_number : DATE_NUMBER
| INT_NUMBER
;
signed_number : PLUS non_negative_number
{ $$ = $2; }
| MINUS non_negative_number
@ -970,6 +983,18 @@ signed_number : PLUS non_negative_number
| non_negative_number
;
signed_inf : PLUS INF_CONSTANT
{ $$ = new string ("Inf"); }
| MINUS INF_CONSTANT
{ $$ = new string ("-Inf"); }
| INF_CONSTANT
{ $$ = new string ("Inf"); }
;
signed_number_w_inf : signed_inf
| signed_number
;
estimated_params : ESTIMATED_PARAMS ';' estimated_list END ';' { driver.estimated_params(); };
estimated_list : estimated_list estimated_elem
@ -1007,24 +1032,21 @@ estimated_elem1 : STDERR symbol
}
;
estimated_elem2 : prior COMMA estimated_elem3
estimated_elem2 : prior_pdf COMMA estimated_elem3
{
driver.estim_params.prior = *$1;
delete $1;
driver.estim_params.prior = $1;
}
| expression_or_empty COMMA prior COMMA estimated_elem3
| expression_or_empty COMMA prior_pdf COMMA estimated_elem3
{
driver.estim_params.init_val = $1;
driver.estim_params.prior = *$3;
delete $3;
driver.estim_params.prior = $3;
}
| expression_or_empty COMMA expression_or_empty COMMA expression_or_empty COMMA prior COMMA estimated_elem3
| expression_or_empty COMMA expression_or_empty COMMA expression_or_empty COMMA prior_pdf COMMA estimated_elem3
{
driver.estim_params.init_val = $1;
driver.estim_params.low_bound = $3;
driver.estim_params.up_bound = $5;
driver.estim_params.prior = *$7;
delete $7;
driver.estim_params.prior = $7;
}
| expression
{
@ -1137,22 +1159,119 @@ estimated_bounds_elem : STDERR symbol COMMA expression COMMA expression ';'
}
;
prior : BETA_PDF
{ $$ = new string("1"); }
| GAMMA_PDF
{ $$ = new string("2"); }
| NORMAL_PDF
{ $$ = new string("3"); }
| INV_GAMMA_PDF
{ $$ = new string("4"); }
| INV_GAMMA1_PDF
{ $$ = new string("4"); }
| UNIFORM_PDF
{ $$ = new string("5"); }
| INV_GAMMA2_PDF
{ $$ = new string("6"); }
prior_distribution : BETA
{ $$ = eBeta; }
| GAMMA
{ $$ = eGamma; }
| NORMAL
{ $$ = eNormal; }
| INV_GAMMA
{ $$ = eInvGamma; }
| INV_GAMMA1
{ $$ = eInvGamma1; }
| UNIFORM
{ $$ = eUniform; }
| INV_GAMMA2
{ $$ = eInvGamma2; }
;
prior_pdf : BETA_PDF
{ $$ = eBeta; }
| GAMMA_PDF
{ $$ = eGamma; }
| NORMAL_PDF
{ $$ = eNormal; }
| INV_GAMMA_PDF
{ $$ = eInvGamma; }
| INV_GAMMA1_PDF
{ $$ = eInvGamma1; }
| UNIFORM_PDF
{ $$ = eUniform; }
| INV_GAMMA2_PDF
{ $$ = eInvGamma2; }
;
set_time : SET_TIME '(' date_number ')' ';'
{ driver.set_time($3); }
;
data : DATA '(' data_options_list ')'';'
{ driver.estimation_data(); }
;
data_options_list : data_options_list COMMA data_options
| data_options
;
data_options : o_file
| o_new_estimation_data_first_obs
| o_last_obs
| o_new_estimation_data_nobs
| o_xls_sheet
| o_xls_range
;
subsamples : symbol '.' SUBSAMPLES '(' subsamples_name_list ')' ';'
{ driver.set_subsamples($1); }
;
subsamples_eq : symbol '.' SUBSAMPLES EQUAL symbol '.' SUBSAMPLES ';'
{ driver.copy_subsamples($1, $5); }
;
subsamples_name_list : subsamples_name_list COMMA o_subsample_name
| o_subsample_name
;
prior : symbol '.' PRIOR { driver.set_prior_variance(); driver.prior_shape = eNoShape; } '(' prior_options_list ')' ';'
{ driver.set_prior($1); }
| symbol '.' symbol '.' PRIOR { driver.set_prior_variance(); driver.prior_shape = eNoShape; } '(' prior_options_list ')' ';'
{
driver.add_subsample_range(new string (*$1), $3);
driver.set_prior($1);
}
| STD '(' symbol ')' '.' PRIOR { driver.set_prior_variance(); driver.prior_shape = eNoShape; } '(' prior_options_list ')' ';'
{ driver.set_std_prior($3); }
| CORR '(' symbol COMMA symbol')' '.' PRIOR { driver.set_prior_variance(); driver.prior_shape = eNoShape; } '(' prior_options_list ')' ';'
{ driver.set_corr_prior($3, $5); }
;
prior_options_list : prior_options_list COMMA prior_options
| prior_options
;
prior_options : o_shift
| o_mean
| o_stdev
| o_variance
| o_mode
| o_interval
| o_shape
| o_domain
;
options : symbol '.' OPTIONS '(' options_options_list ')' ';'
{ driver.set_options($1); }
| symbol '.' symbol '.' OPTIONS '(' options_options_list ')' ';'
{
driver.add_subsample_range(new string (*$1), $3);
driver.set_options($1);
}
| STD '(' symbol ')' '.' OPTIONS '(' options_options_list ')' ';'
{ driver.set_std_options($3); }
| CORR '(' symbol COMMA symbol')' '.' OPTIONS '(' options_options_list ')' ';'
{ driver.set_corr_options($3, $5); }
;
options_options_list : options_options_list COMMA options_options
| options_options
;
options_options : o_jscale
| o_init
| o_bounds
;
estimation : ESTIMATION ';'
{ driver.run_estimation(); }
| ESTIMATION '(' estimation_options_list ')' ';'
@ -1852,6 +1971,7 @@ o_mfs : MFS EQUAL INT_NUMBER { driver.mfs($3); };
o_simul : SIMUL; // Do nothing, only here for backward compatibility
o_simul_seed : SIMUL_SEED EQUAL INT_NUMBER { driver.error("'simul_seed' option is no longer supported; use 'set_dynare_seed' command instead"); } ;
o_qz_criterium : QZ_CRITERIUM EQUAL non_negative_number { driver.option_num("qz_criterium", $3); };
o_file : FILE EQUAL filename { driver.option_str("file", $3); };
o_datafile : DATAFILE EQUAL filename { driver.option_str("datafile", $3); };
o_nobs : NOBS EQUAL vec_int
{ driver.option_vec_int("nobs", $3); }
@ -1864,6 +1984,24 @@ o_conditional_variance_decomposition : CONDITIONAL_VARIANCE_DECOMPOSITION EQUAL
{ driver.option_vec_int("conditional_variance_decomposition", $3); }
;
o_first_obs : FIRST_OBS EQUAL INT_NUMBER { driver.option_num("first_obs", $3); };
o_new_estimation_data_first_obs : FIRST_OBS EQUAL date_number
{ driver.option_date("first_obs", $3); }
;
o_last_obs : LAST_OBS EQUAL date_number
{ driver.option_date("last_obs", $3); }
;
o_shift : SHIFT EQUAL signed_number { driver.option_num("shift", $3); };
o_shape : SHAPE EQUAL prior_distribution { driver.prior_shape = $3; };
o_mode : MODE EQUAL signed_number { driver.option_num("mode", $3); };
o_mean : MEAN EQUAL signed_number { driver.option_num("mean", $3); };
o_stdev : STDEV EQUAL non_negative_number { driver.option_num("stdev", $3); };
o_jscale : JSCALE EQUAL non_negative_number { driver.option_num("jscale", $3); };
o_init : INIT EQUAL signed_number { driver.option_num("init", $3); };
o_bounds : BOUNDS EQUAL vec_value_w_inf { driver.option_num("bounds", $3); };
o_domain : DOMAINN EQUAL vec_value { driver.option_num("domain", $3); };
o_interval : INTERVAL EQUAL vec_value { driver.option_num("interval", $3); };
o_variance : VARIANCE EQUAL expression { driver.set_prior_variance($3); }
o_new_estimation_data_nobs : NOBS EQUAL INT_NUMBER { driver.option_num("nobs", $3); };
o_prefilter : PREFILTER EQUAL INT_NUMBER { driver.option_num("prefilter", $3); };
o_presample : PRESAMPLE EQUAL INT_NUMBER { driver.option_num("presample", $3); };
o_lik_algo : LIK_ALGO EQUAL INT_NUMBER { driver.option_num("lik_algo", $3); };
@ -1873,6 +2011,12 @@ o_nograph : NOGRAPH
| GRAPH
{ driver.option_num("nograph", "0"); }
;
o_subsample_name : symbol EQUAL date_number ':' date_number
{
driver.declare_statement_local_variable(new string (*$1));
driver.set_subsample_name_equal_to_date_range($1, $3, $5);
}
;
o_conf_sig : CONF_SIG EQUAL non_negative_number { driver.option_num("conf_sig", $3); };
o_mh_replic : MH_REPLIC EQUAL INT_NUMBER { driver.option_num("mh_replic", $3); };
o_mh_drop : MH_DROP EQUAL non_negative_number { driver.option_num("mh_drop", $3); };
@ -2075,13 +2219,15 @@ o_cnum : CNUM EQUAL INT_NUMBER {driver.option_num("ms.cnum",$3); };
o_k_order_solver : K_ORDER_SOLVER {driver.option_num("k_order_solver","1"); };
o_pruning : PRUNING { driver.option_num("pruning", "1"); };
o_chain : CHAIN EQUAL INT_NUMBER { driver.option_num("ms.chain",$3); };
o_state : STATE EQUAL INT_NUMBER { driver.option_num("ms.state",$3); };
o_restrictions : RESTRICTIONS EQUAL vec_of_vec_value
{ driver.option_num("ms.restrictions",$3); }
;
o_duration : DURATION EQUAL non_negative_number
{ driver.option_num("ms.duration",$3); }
| DURATION EQUAL INF_CONSTANT
{ driver.option_num("ms.duration","Inf"); }
| DURATION EQUAL vec_value_w_inf
{ driver.option_num("ms.duration",$3); }
;
o_number_of_states : NUMBER_OF_STATES EQUAL INT_NUMBER { driver.option_num("ms.number_of_states",$3); };
o_number_of_regimes : NUMBER_OF_REGIMES EQUAL INT_NUMBER { driver.option_num("ms.number_of_regimes",$3); };
o_coefficients : COEFFICIENTS { driver.option_str("ms.coefficients","svar_coefficients"); };
o_variances : VARIANCES { driver.option_str("ms.variances","svar_variances"); };
o_equations : EQUATIONS EQUAL vec_int
@ -2219,18 +2365,54 @@ vec_int : vec_int_1 ']'
{ $$ = $1; }
;
vec_value_1 : '[' signed_number
{ $2->insert(0, "["); $$ = $2;}
| vec_value_1 signed_number
{
$1->append(" ");
$1->append(*$2);
delete $2;
$$ = $1;
}
vec_value_1 : '[' signed_number { $2->insert(0,"["); $$ = $2; }
| '[' COMMA signed_number { $3->insert(0,"["); $$ = $3; }
| vec_value_1 signed_number
{
$1->append(" ");
$1->append(*$2);
delete $2;
$$ = $1;
}
| vec_value_1 COMMA signed_number
{
$1->append(" ");
$1->append(*$3);
delete $3;
$$ = $1;
}
;
vec_value : vec_value_1 ']' { $1->append("]"); $$ = $1; }
| vec_value_1 COMMA ']' { $1->append("]"); $$ = $1; }
;
vec_value : vec_value_1 ']' { $1->append("]"); $$ = $1; };
vec_value_list : vec_value_list COMMA vec_value
{
$1->append(",");
$1->append(*$3);
delete $3;
$$ = $1;
}
| vec_value
{ $$ = $1; }
;
vec_of_vec_value : '[' vec_value_list ']' { $$ = $2; }
| vec_value { $$ = $1; };
vec_value_1_w_inf : '[' signed_number_w_inf
{ $2->insert(0, "["); $$ = $2;}
| vec_value_1_w_inf signed_number_w_inf
{
$1->append(" ");
$1->append(*$2);
delete $2;
$$ = $1;
}
;
vec_value_w_inf : vec_value_1_w_inf ']' { $1->append("]"); $$ = $1; };
symbol : NAME
| ALPHA
@ -2240,6 +2422,12 @@ symbol : NAME
| CMS
| NCMS
| CNUM
| GAMMA
| INV_GAMMA
| INV_GAMMA1
| INV_GAMMA2
| NORMAL
| UNIFORM
;
%%

View File

@ -108,6 +108,8 @@ string eofbuff;
<INITIAL>periods {BEGIN DYNARE_STATEMENT; return token::PERIODS;}
<INITIAL>model_info {BEGIN DYNARE_STATEMENT; return token::MODEL_INFO;}
<INITIAL>estimation {BEGIN DYNARE_STATEMENT; return token::ESTIMATION;}
<INITIAL>set_time {BEGIN DYNARE_STATEMENT; return token::SET_TIME;}
<INITIAL>data {BEGIN DYNARE_STATEMENT; return token::DATA;}
<INITIAL>varobs {BEGIN DYNARE_STATEMENT; return token::VAROBS;}
<INITIAL>unit_root_vars {BEGIN DYNARE_STATEMENT; return token::UNIT_ROOT_VARS;}
<INITIAL>rplot {BEGIN DYNARE_STATEMENT; return token::RPLOT;}
@ -187,10 +189,29 @@ string eofbuff;
/* End of a Dynare block */
<DYNARE_BLOCK>end {BEGIN INITIAL; return token::END;}
<DYNARE_STATEMENT>subsamples {return token::SUBSAMPLES;}
<DYNARE_STATEMENT>options {return token::OPTIONS;}
<DYNARE_STATEMENT>prior {return token::PRIOR;}
<INITIAL>std {BEGIN DYNARE_STATEMENT; return token::STD;}
<INITIAL>corr {BEGIN DYNARE_STATEMENT; return token::CORR;}
/* Inside of a Dynare statement */
<DYNARE_STATEMENT>file {return token::FILE;}
<DYNARE_STATEMENT>datafile {return token::DATAFILE;}
<DYNARE_STATEMENT>nobs {return token::NOBS;}
<DYNARE_STATEMENT>last_obs {return token::LAST_OBS;}
<DYNARE_STATEMENT>first_obs {return token::FIRST_OBS;}
<DYNARE_STATEMENT>mean {return token::MEAN;}
<DYNARE_STATEMENT>stdev {return token::STDEV;}
<DYNARE_STATEMENT>domain {return token::DOMAINN;}
<DYNARE_STATEMENT>variance {return token::VARIANCE;}
<DYNARE_STATEMENT>mode {return token::MODE;}
<DYNARE_STATEMENT>interval {return token::INTERVAL;}
<DYNARE_STATEMENT>shape {return token::SHAPE;}
<DYNARE_STATEMENT>shift {return token::SHIFT;}
<DYNARE_STATEMENT>bounds {return token::BOUNDS;}
<DYNARE_STATEMENT>init {return token::INIT;}
<DYNARE_STATEMENT>jscale {return token::JSCALE;}
<DYNARE_STATEMENT>prefilter {return token::PREFILTER;}
<DYNARE_STATEMENT>presample {return token::PRESAMPLE;}
<DYNARE_STATEMENT>lik_algo {return token::LIK_ALGO;}
@ -249,7 +270,6 @@ string eofbuff;
<DYNARE_STATEMENT>nargs {return token::EXT_FUNC_NARGS;}
<DYNARE_STATEMENT>first_deriv_provided {return token::FIRST_DERIV_PROVIDED;}
<DYNARE_STATEMENT>second_deriv_provided {return token::SECOND_DERIV_PROVIDED;}
<DYNARE_STATEMENT>freq {return token::FREQ;}
<DYNARE_STATEMENT>monthly {return token::MONTHLY; }
<DYNARE_STATEMENT>quarterly {return token::QUARTERLY; }
@ -262,6 +282,7 @@ string eofbuff;
<DYNARE_STATEMENT>vlistper {return token::VLISTPER;}
<DYNARE_STATEMENT>restriction_fname {return token::RESTRICTION_FNAME;}
<DYNARE_STATEMENT>nlags {return token::NLAGS;}
<DYNARE_STATEMENT>restrictions {return token::RESTRICTIONS;}
<DYNARE_STATEMENT>cross_restrictions {return token::CROSS_RESTRICTIONS;}
<DYNARE_STATEMENT>contemp_reduced_form {return token::CONTEMP_REDUCED_FORM;}
<DYNARE_STATEMENT>real_pseudo_forecast {return token::REAL_PSEUDO_FORECAST;}
@ -277,6 +298,30 @@ string eofbuff;
yylval->string_val = new string(yytext);
return token::BETA;
}
<DYNARE_STATEMENT>gamma {
yylval->string_val = new string(yytext);
return token::GAMMA;
}
<DYNARE_STATEMENT>inv_gamma {
yylval->string_val = new string(yytext);
return token::INV_GAMMA;
}
<DYNARE_STATEMENT>inv_gamma1 {
yylval->string_val = new string(yytext);
return token::INV_GAMMA1;
}
<DYNARE_STATEMENT>inv_gamma2 {
yylval->string_val = new string(yytext);
return token::INV_GAMMA2;
}
<DYNARE_STATEMENT>normal {
yylval->string_val = new string(yytext);
return token::NORMAL;
}
<DYNARE_STATEMENT>uniform {
yylval->string_val = new string(yytext);
return token::UNIFORM;
}
<DYNARE_STATEMENT>gsig2_lmdm {return token::GSIG2_LMDM;}
<DYNARE_STATEMENT>specification {return token::SPECIFICATION;}
<DYNARE_STATEMENT>sims_zha {return token::SIMS_ZHA;}
@ -427,7 +472,6 @@ string eofbuff;
<DYNARE_BLOCK># {return Dynare::parser::token_type (yytext[0]);}
<DYNARE_BLOCK>autocorr {return token::AUTOCORR;}
<DYNARE_BLOCK>restrictions {return token::RESTRICTIONS;}
<DYNARE_BLOCK>restriction {return token::RESTRICTION;}
/* Inside Dynare statement */
@ -461,8 +505,7 @@ string eofbuff;
<DYNARE_STATEMENT,DYNARE_BLOCK>upper_cholesky {return token::UPPER_CHOLESKY;}
<DYNARE_STATEMENT,DYNARE_BLOCK>lower_cholesky {return token::LOWER_CHOLESKY;}
<DYNARE_STATEMENT>chain {return token::CHAIN;}
<DYNARE_STATEMENT>state {return token::STATE;}
<DYNARE_STATEMENT>number_of_states {return token::NUMBER_OF_STATES;}
<DYNARE_STATEMENT>number_of_regimes {return token::NUMBER_OF_REGIMES;}
<DYNARE_STATEMENT>duration {return token::DURATION;}
<DYNARE_STATEMENT>coefficients {return token::COEFFICIENTS;}
<DYNARE_STATEMENT>variances {return token::VARIANCES;}
@ -591,6 +634,11 @@ string eofbuff;
return token::INT_NUMBER;
}
<DYNARE_STATEMENT,DYNARE_BLOCK>([1-2][0-9]{3}[Mm](([1-9])|(1[0-2])))|([1-2][0-9]{3}[Qq][1-4])|([1-2][0-9]{3}[Ww](([1-9]{1})|([1-5][0-9]))) {
yylval->string_val = new string(yytext);
return token::DATE_NUMBER;
}
<DYNARE_STATEMENT,DYNARE_BLOCK>\'[^\']+\' {
yylval->string_val = new string(yytext + 1);
yylval->string_val->resize(yylval->string_val->length() - 1);

View File

@ -496,6 +496,7 @@ VariableNode::prepareForDerivation()
non_null_derivatives = datatree.local_variables_table[symb_id]->non_null_derivatives;
break;
case eModFileLocalVariable:
case eStatementDeclaredVariable:
// Such a variable is never derived
break;
case eExternalFunction:
@ -523,6 +524,9 @@ VariableNode::computeDerivative(int deriv_id)
case eModFileLocalVariable:
cerr << "ModFileLocalVariable is not derivable" << endl;
exit(EXIT_FAILURE);
case eStatementDeclaredVariable:
cerr << "eStatementDeclaredVariable is not derivable" << endl;
exit(EXIT_FAILURE);
case eExternalFunction:
cerr << "Impossible case!" << endl;
exit(EXIT_FAILURE);
@ -735,6 +739,7 @@ VariableNode::writeOutput(ostream &output, ExprNodeOutputType output_type,
case eExternalFunction:
case eTrend:
case eStatementDeclaredVariable:
cerr << "Impossible case" << endl;
exit(EXIT_FAILURE);
}
@ -917,6 +922,9 @@ VariableNode::getChainRuleDerivative(int deriv_id, const map<int, expr_t> &recur
case eModFileLocalVariable:
cerr << "ModFileLocalVariable is not derivable" << endl;
exit(EXIT_FAILURE);
case eStatementDeclaredVariable:
cerr << "eStatementDeclaredVariable is not derivable" << endl;
exit(EXIT_FAILURE);
case eExternalFunction:
cerr << "Impossible case!" << endl;
exit(EXIT_FAILURE);

View File

@ -449,7 +449,7 @@ ModFile::writeOutputFiles(const string &basename, bool clear_all, bool console,
}
mOutputFile << "tic;" << endl
<< "global M_ oo_ options_ ys0_ ex0_" << endl
<< "global M_ oo_ options_ ys0_ ex0_ estimation_info" << endl
<< "options_ = [];" << endl
<< "M_.fname = '" << basename << "';" << endl
<< "%" << endl
@ -545,31 +545,33 @@ ModFile::writeOutputFiles(const string &basename, bool clear_all, bool console,
#if defined(_WIN32) || defined(__CYGWIN32__)
if (msvc)
// MATLAB/Windows + Microsoft Visual C++
mOutputFile << " eval('mex -O LINKFLAGS=\"$LINKFLAGS /export:Dynamic\" " << basename << "_dynamic.c')" << endl
<< " eval('mex -O LINKFLAGS=\"$LINKFLAGS /export:Dynamic\" " << basename << "_static.c')" << endl;
mOutputFile << " eval('mex -O LINKFLAGS=\"$LINKFLAGS /export:Dynamic\" " << basename << "_dynamic.c " << basename << "_dynamic_mex.c')" << endl
<< " eval('mex -O LINKFLAGS=\"$LINKFLAGS /export:Dynamic\" " << basename << "_static.c "<< basename << "_static_mex.c')" << endl;
else if (cygwin)
// MATLAB/Windows + Cygwin g++
mOutputFile << " eval('mex -O PRELINK_CMDS1=\"echo EXPORTS > mex.def & echo mexFunction >> mex.def & echo Dynamic >> mex.def\" " << basename << "_dynamic.c')" << endl
<< " eval('mex -O PRELINK_CMDS1=\"echo EXPORTS > mex.def & echo mexFunction >> mex.def & echo Dynamic >> mex.def\" " << basename << "_static.c')" << endl;
mOutputFile << " eval('mex -O PRELINK_CMDS1=\"echo EXPORTS > mex.def & echo mexFunction >> mex.def & echo Dynamic >> mex.def\" " << basename << "_dynamic.c " << basename << "_dynamic_mex.c')" << endl
<< " eval('mex -O PRELINK_CMDS1=\"echo EXPORTS > mex.def & echo mexFunction >> mex.def & echo Dynamic >> mex.def\" " << basename << "_static.c "<< basename << "_static_mex.c')" << endl;
else
mOutputFile << " error('When using the USE_DLL option, you must give either ''cygwin'' or ''msvc'' option to the ''dynare'' command')" << endl;
#else
# ifdef __linux__
// MATLAB/Linux
mOutputFile << " eval('mex -O LDFLAGS=''-pthread -shared -Wl,--no-undefined'' " << basename << "_dynamic.c')" << endl
<< " eval('mex -O LDFLAGS=''-pthread -shared -Wl,--no-undefined'' " << basename << "_static.c')" << endl;
mOutputFile << " eval('mex -O LDFLAGS=''-pthread -shared -Wl,--no-undefined'' " << basename << "_dynamic.c " << basename << "_dynamic_mex.c')" << endl
<< " eval('mex -O LDFLAGS=''-pthread -shared -Wl,--no-undefined'' " << basename << "_static.c "<< basename << "_static_mex.c')" << endl;
# else // MacOS
// MATLAB/MacOS
mOutputFile << " eval('mex -O LDFLAGS=''-Wl,-twolevel_namespace -undefined error -arch \\$ARCHS -Wl,-syslibroot,\\$SDKROOT -mmacosx-version-min=\\$MACOSX_DEPLOYMENT_TARGET -bundle'' " << basename << "_dynamic.c')" << endl
<< " eval('mex -O LDFLAGS=''-Wl,-twolevel_namespace -undefined error -arch \\$ARCHS -Wl,-syslibroot,\\$SDKROOT -mmacosx-version-min=\\$MACOSX_DEPLOYMENT_TARGET -bundle'' " << basename << "_static.c')" << endl;
mOutputFile << " eval('mex -O LDFLAGS=''-Wl,-twolevel_namespace -undefined error -arch \\$ARCHS -Wl,-syslibroot,\\$SDKROOT -mmacosx-version-min=\\$MACOSX_DEPLOYMENT_TARGET -bundle'' "
<< basename << "_dynamic.c " << basename << "_dynamic_mex.c')" << endl
<< " eval('mex -O LDFLAGS=''-Wl,-twolevel_namespace -undefined error -arch \\$ARCHS -Wl,-syslibroot,\\$SDKROOT -mmacosx-version-min=\\$MACOSX_DEPLOYMENT_TARGET -bundle'' "
<< basename << "_static.c " << basename << "_static_mex.c')" << endl;
# endif
#endif
mOutputFile << "else" << endl // Octave
<< " if ~octave_ver_less_than('3.2.0')" << endl // Workaround for bug in Octave >= 3.2, see http://bugs.debian.org/cgi-bin/bugreport.cgi?bug=550823
<< " sleep(2)" << endl
<< " end" << endl
<< " mex " << basename << "_dynamic.c" << endl
<< " mex " << basename << "_static.c" << endl
<< " mex " << basename << "_dynamic.c " << basename << "_dynamic_mex.c" << endl
<< " mex " << basename << "_static.c " << basename << "_static_mex.c" << endl
<< "end" << endl;
}

View File

@ -132,8 +132,7 @@ InitValStatement::writeOutput(ostream &output, const string &basename) const
void
InitValStatement::writeOutputPostInit(ostream &output) const
{
output << "oo_.endo_simul=[oo_.steady_state*ones(1,M_.maximum_lag)];" << endl
<< "if M_.exo_nbr > 0;" << endl
output << "if M_.exo_nbr > 0;" << endl
<< "\too_.exo_simul = [ones(M_.maximum_lag,1)*oo_.exo_steady_state'];" << endl
<<"end;" << endl
<< "if M_.exo_det_nbr > 0;" << endl
@ -189,7 +188,7 @@ HistValStatement::writeOutput(ostream &output, const string &basename) const
output << "%" << endl
<< "% HISTVAL instructions" << endl
<< "%" << endl
<< "oo_.endo_simul = zeros(M_.endo_nbr,M_.maximum_lag);" << endl;
<< "M_.endo_histval = zeros(M_.endo_nbr,M_.maximum_lag);" << endl;
for (hist_values_t::const_iterator it = hist_values.begin();
it != hist_values.end(); it++)
@ -225,7 +224,7 @@ HistValStatement::writeOutput(ostream &output, const string &basename) const
int tsid = symbol_table.getTypeSpecificID(symb_id) + 1;
if (type == eEndogenous)
output << "oo_.endo_simul( " << tsid << ", M_.maximum_lag + " << lag << ") = ";
output << "M_.endo_histval( " << tsid << ", M_.maximum_lag + " << lag << ") = ";
else if (type == eExogenous)
output << "oo_.exo_simul( M_.maximum_lag + " << lag << ", " << tsid << " ) = ";
else if (type != eExogenousDet)

View File

@ -71,7 +71,7 @@ public:
InitValStatement(const init_values_t &init_values_arg,
const SymbolTable &symbol_table_arg);
virtual void writeOutput(ostream &output, const string &basename) const;
//! Writes initializations for oo_.endo_simul, oo_.exo_simul and oo_.exo_det_simul
//! Writes initializations for oo_.exo_simul and oo_.exo_det_simul
void writeOutputPostInit(ostream &output) const;
};

View File

@ -46,6 +46,15 @@ ParsingDriver::check_symbol_existence(const string &name)
error("Unknown symbol: " + name);
}
void
ParsingDriver::check_symbol_is_parameter(string *name)
{
check_symbol_existence(*name);
int symb_id = mod_file->symbol_table.getID(*name);
if (mod_file->symbol_table.getType(symb_id) != eParameter)
error(*name + " is not a parameter");
}
void
ParsingDriver::set_current_data_tree(DataTree *data_tree_arg)
{
@ -166,6 +175,16 @@ ParsingDriver::declare_parameter(string *name, string *tex_name)
delete tex_name;
}
void
ParsingDriver::declare_statement_local_variable(string *name)
{
if (mod_file->symbol_table.exists(*name))
error("Symbol " + *name + " cannot be assigned within a statement " +
"while being assigned elsewhere in the modfile");
declare_symbol(name, eStatementDeclaredVariable, NULL);
delete name;
}
void
ParsingDriver::declare_optimal_policy_discount_factor_parameter(expr_t exprnode)
{
@ -372,13 +391,9 @@ ParsingDriver::dsample(string *arg1, string *arg2)
void
ParsingDriver::init_param(string *name, expr_t rhs)
{
check_symbol_existence(*name);
check_symbol_is_parameter(name);
int symb_id = mod_file->symbol_table.getID(*name);
if (mod_file->symbol_table.getType(symb_id) != eParameter)
error(*name + " is not a parameter");
mod_file->addStatement(new InitParamStatement(symb_id, rhs, mod_file->symbol_table));
delete name;
}
@ -1056,6 +1071,23 @@ ParsingDriver::option_str(const string &name_option, const string &opt)
options_list.string_options[name_option] = opt;
}
void
ParsingDriver::option_date(const string &name_option, string *opt)
{
option_date(name_option, *opt);
delete opt;
}
void
ParsingDriver::option_date(const string &name_option, const string &opt)
{
if (options_list.date_options.find(name_option)
!= options_list.date_options.end())
error("option " + name_option + " declared twice");
options_list.date_options[name_option] = opt;
}
void
ParsingDriver::option_symbol_list(const string &name_option)
{
@ -1184,6 +1216,177 @@ ParsingDriver::set_unit_root_vars()
symbol_list.clear();
}
void
ParsingDriver::set_time(string *arg)
{
string arg1 = *arg;
for (size_t i=0; i<arg1.length(); i++)
arg1[i]= toupper(arg1[i]);
option_date("initial_period", arg1);
mod_file->addStatement(new SetTimeStatement(options_list));
options_list.clear();
}
void
ParsingDriver::estimation_data()
{
mod_file->addStatement(new EstimationDataStatement(options_list));
options_list.clear();
}
void
ParsingDriver::copy_subsamples(string *to_parameter, string *from_parameter)
{
check_symbol_is_parameter(to_parameter);
check_symbol_is_parameter(from_parameter);
if (subsample_declarations.find(*to_parameter) != subsample_declarations.end())
error("Parameter " + *to_parameter + " has more than one subsample statement." +
"You may only have one subsample statement per parameter.");
if (subsample_declarations.find(*from_parameter) == subsample_declarations.end())
error("Parameter " + *from_parameter + " does not have an associated subsample statement.");
subsample_declarations[*to_parameter] = subsample_declarations[*from_parameter];
delete to_parameter;
delete from_parameter;
}
void
ParsingDriver::set_subsamples(string *name)
{
check_symbol_is_parameter(name);
if (subsample_declarations.find(*name) != subsample_declarations.end())
error("Parameter " + *name + " has more than one subsample statement." +
"You may only have one subsample statement per parameter.");
subsample_declarations[*name] = subsample_declaration_map;
subsample_declaration_map.clear();
delete name;
}
void
ParsingDriver::check_symbol_is_statement_variable(string *name)
{
check_symbol_existence(*name);
int symb_id = mod_file->symbol_table.getID(*name);
if (mod_file->symbol_table.getType(symb_id) != eStatementDeclaredVariable)
error(*name + " is not a variable assigned in a statement");
}
void
ParsingDriver::set_subsample_name_equal_to_date_range(string *name, string *date1, string *date2)
{
check_symbol_is_statement_variable(name);
if (subsample_declaration_map.find(*name) != subsample_declaration_map.end())
error("Symbol " + *name + " may only be assigned once in a SUBSAMPLE statement");
subsample_declaration_map[*name] = make_pair(*date1, *date2);
delete name;
delete date1;
delete date2;
}
void
ParsingDriver::add_subsample_range(string *parameter, string *subsample_name)
{
check_symbol_is_parameter(parameter);
check_symbol_is_statement_variable(subsample_name);
subsample_declarations_t::const_iterator it = subsample_declarations.find(*parameter);
if (it == subsample_declarations.end())
error("A subsample statement has not been issued for " + *parameter);
subsample_declaration_map_t tmp_map = it->second;
if (tmp_map.find(*subsample_name) == tmp_map.end())
error("The subsample name " + *subsample_name + " was not previously declared in a subsample statement.");
option_date("date1", tmp_map[*subsample_name].first);
option_date("date2", tmp_map[*subsample_name].second);
delete parameter;
delete subsample_name;
}
void
ParsingDriver::set_prior(string *name)
{
check_symbol_is_parameter(name);
mod_file->addStatement(new PriorStatement(*name, prior_shape, prior_variance, options_list));
options_list.clear();
set_prior_variance();
prior_shape = eNoShape;
delete name;
}
void
ParsingDriver::set_prior_variance(expr_t variance)
{
prior_variance = variance;
}
void
ParsingDriver::set_options(string *name)
{
check_symbol_is_parameter(name);
mod_file->addStatement(new OptionsStatement(*name, options_list));
options_list.clear();
delete name;
}
void
ParsingDriver::check_symbol_is_endogenous_or_exogenous(string *name)
{
check_symbol_existence(*name);
int symb_id = mod_file->symbol_table.getID(*name);
switch(mod_file->symbol_table.getType(symb_id))
{
case eEndogenous:
case eExogenous:
case eExogenousDet:
break;
default:
error(*name + " is neither endogenous or exogenous.");
}
}
void
ParsingDriver::set_std_prior(string *name)
{
check_symbol_is_endogenous_or_exogenous(name);
mod_file->addStatement(new StdPriorStatement(*name, prior_shape, prior_variance,
options_list, mod_file->symbol_table));
options_list.clear();
set_prior_variance();
prior_shape = eNoShape;
delete name;
}
void
ParsingDriver::set_std_options(string *name)
{
check_symbol_is_endogenous_or_exogenous(name);
// mod_file->addStatement(new StdOptionsStatement(*name, options_list, mod_file->symbol_table));
options_list.clear();
delete name;
}
void
ParsingDriver::set_corr_prior(string *name1, string *name2)
{
check_symbol_is_endogenous_or_exogenous(name1);
check_symbol_is_endogenous_or_exogenous(name2);
mod_file->addStatement(new CorrPriorStatement(*name1, *name2, prior_shape, prior_variance,
options_list, mod_file->symbol_table));
options_list.clear();
set_prior_variance();
prior_shape = eNoShape;
delete name1;
delete name2;
}
void
ParsingDriver::set_corr_options(string *name1, string *name2)
{
check_symbol_is_endogenous_or_exogenous(name1);
check_symbol_is_endogenous_or_exogenous(name2);
// mod_file->addStatement(new CorrOptionsStatement(*name1, *name2, options_list, mod_file->symbol_table));
options_list.clear();
delete name1;
delete name2;
}
void
ParsingDriver::run_estimation()
{
@ -1540,7 +1743,7 @@ ParsingDriver::svar()
void
ParsingDriver::markov_switching()
{
OptionsList::num_options_t::const_iterator it0, it1;
OptionsList::num_options_t::const_iterator it0;
it0 = options_list.num_options.find("ms.chain");
if (it0 == options_list.num_options.end())
@ -1548,31 +1751,15 @@ ParsingDriver::markov_switching()
else if (atoi(it0->second.c_str()) <= 0)
error("The value passed to the chain option must be greater than zero.");
it0 = options_list.num_options.find("ms.state");
it1 = options_list.num_options.find("ms.number_of_states");
if ((it0 == options_list.num_options.end())
&& (it1 == options_list.num_options.end()))
error("Either a state option or a number_of_states option must be passed to the markov_switching statement.");
it0 = options_list.num_options.find("ms.number_of_regimes");
if (it0 == options_list.num_options.end())
error("A number_of_regimes option must be passed to the markov_switching statement.");
else if (atoi(it0->second.c_str()) <= 0)
error("The value passed to the number_of_regimes option must be greater than zero.");
if ((it0 != options_list.num_options.end())
&& (it1 != options_list.num_options.end()))
error("You cannot pass both a state option and a number_of_states option to the markov_switching statement.");
if (it0 != options_list.num_options.end())
if (atoi(it0->second.c_str()) <= 0)
error("The value passed to the state option must be greater than zero.");
if (it1 != options_list.num_options.end())
if (atoi(it1->second.c_str()) <= 0)
error("The value passed to the number_of_states option must be greater than zero.");
string infStr("Inf");
it0 = options_list.num_options.find("ms.duration");
if (it0 == options_list.num_options.end())
error("A duration option must be passed to the markov_switching statement.");
else if (infStr.compare(it0->second) != 0)
if (atof(it0->second.c_str()) <= 0.0)
error("The value passed to the duration option must be greater than zero.");
mod_file->addStatement(new MarkovSwitchingStatement(options_list));
options_list.clear();

View File

@ -81,6 +81,15 @@ private:
//! Checks that a given symbol exists, and stops with an error message if it doesn't
void check_symbol_existence(const string &name);
//! Checks that a given symbol exists and is a parameter, and stops with an error message if it isn't
void check_symbol_is_parameter(string *name);
//! Checks that a given symbol was assigned within a Statement
void check_symbol_is_statement_variable(string *name);
//! Checks that a given symbol exists and is a endogenous or exogenous, and stops with an error message if it isn't
void check_symbol_is_endogenous_or_exogenous(string *name);
//! Helper to add a symbol declaration
void declare_symbol(const string *name, SymbolType type, const string *tex_name);
@ -178,6 +187,14 @@ private:
vector<int> declared_trend_vars;
//! Temporary storage for declaring nonstationary variables
vector<int> declared_nonstationary_vars;
//! Temporary storage for a variance declared in the prior statement
expr_t prior_variance;
//! Temporary storage for declaring subsamples: map<statement_local_var, <date1, date2 >
typedef map<string, pair<string, string> > subsample_declaration_map_t;
subsample_declaration_map_t subsample_declaration_map;
//! Temporary storage for subsample statement: map<parameter, subsample_declaration_map >
typedef map<string, subsample_declaration_map_t > subsample_declarations_t;
subsample_declarations_t subsample_declarations;
//! reset the values for temporary storage
void reset_current_external_function_options();
//! Adds a model lagged variable to ModelTree and VariableTable
@ -200,6 +217,9 @@ public:
//! Estimation parameters
EstimationParams estim_params;
//! Temporary storage for the prior shape
PriorDistributions prior_shape;
//! Error handler with explicit location
void error(const Dynare::parser::location_type &l, const string &m) __attribute__ ((noreturn));
//! Error handler using saved location
@ -231,6 +251,16 @@ public:
void declare_exogenous_det(string *name, string *tex_name = NULL);
//! Declares a parameter
void declare_parameter(string *name, string *tex_name = NULL);
//! Declares a statement local variable
void declare_statement_local_variable(string *name);
//! Completes a subsample statement
void set_subsamples(string *name);
//! Declares a subsample, assigning the value to name
void set_subsample_name_equal_to_date_range(string *name, string *date1, string *date2);
//! Adds a subsample range to the list of options for the prior statement
void add_subsample_range(string *parameter, string *subsample_name);
//! Copies the set of subsamples from_parameter to_parameter
void copy_subsamples(string *to_parameter, string *from_parameter);
//! Declares declare_optimal_policy_discount_factor as a parameter and initializes it to exprnode
void declare_optimal_policy_discount_factor_parameter(expr_t exprnode);
//! Adds a predetermined_variable
@ -319,6 +349,10 @@ public:
void option_str(const string &name_option, string *opt);
//! Sets an option to a string value
void option_str(const string &name_option, const string &opt);
//! Sets an option to a date value
void option_date(const string &name_option, string *opt);
//! Sets an option to a date value
void option_date(const string &name_option, const string &opt);
//! Sets an option to a list of symbols (used in conjunction with add_in_symbol_list())
void option_symbol_list(const string &name_option);
//! Sets an option to a vector of integers
@ -351,6 +385,24 @@ public:
void external_function_option(const string &name_option, const string &opt);
//! Add a line in an estimated params block
void add_estimated_params_element();
//! Sets the frequency of the data
void set_time(string *arg);
//! Estimation Data
void estimation_data();
//! Sets the prior for a parameter
void set_prior(string *arg);
//! Adds the variance option to its temporary holding place
void set_prior_variance(expr_t variance=NULL);
//! Sets the options for a parameter
void set_options(string *arg);
//! Sets the prior for estimated std dev
void set_std_prior(string *arg);
//! Sets the options for estimated std dev
void set_std_options(string *arg);
//! Sets the prior for estimated correlation
void set_corr_prior(string *arg1, string *arg2);
//! Sets the options for estimated correlation
void set_corr_options(string *arg1, string *arg2);
//! Runs estimation process
void run_estimation();
//! Runs dynare_sensitivy()

View File

@ -88,6 +88,10 @@ OptionsList::writeOutput(ostream &output) const
it != string_options.end(); it++)
output << "options_." << it->first << " = '" << it->second << "';" << endl;
for (date_options_t::const_iterator it = date_options.begin();
it != date_options.end(); it++)
output << "options_." << it->first << " = dynDate('" << it->second << "');" << endl;
for (symbol_list_options_t::const_iterator it = symbol_list_options.begin();
it != symbol_list_options.end(); it++)
it->second.writeOutput("options_." + it->first, output);
@ -127,6 +131,10 @@ OptionsList::writeOutput(ostream &output, const string &option_group) const
it != string_options.end(); it++)
output << option_group << "." << it->first << " = '" << it->second << "';" << endl;
for (date_options_t::const_iterator it = date_options.begin();
it != date_options.end(); it++)
output << option_group << "." << it->first << " = dynDate('" << it->second << "');" << endl;
for (symbol_list_options_t::const_iterator it = symbol_list_options.begin();
it != symbol_list_options.end(); it++)
it->second.writeOutput(option_group + "." + it->first, output);
@ -154,6 +162,7 @@ OptionsList::clear()
num_options.clear();
paired_num_options.clear();
string_options.clear();
date_options.clear();
symbol_list_options.clear();
vector_int_options.clear();
}

View File

@ -86,6 +86,20 @@ public:
bool dsge_var_estimated;
//! Whether there is a bayesian_irf option passed to the estimation statement
bool bayesian_irf_present;
//! Whether there is a data statement present
bool estimation_data_statement_present;
//! Whether there is a prior statement present
bool prior_statement_present;
//! Whether there is a std prior statement present
bool std_prior_statement_present;
//! Whether there is a corr prior statement present
bool corr_prior_statement_present;
//! Whether there is a options statement present
bool options_statement_present;
//! Whether there is a std options statement present
bool std_options_statement_present;
//! Whether there is a corr options statement present
bool corr_options_statement_present;
};
class Statement
@ -118,11 +132,13 @@ public:
typedef map<string, string> num_options_t;
typedef map<string, pair<string, string> > paired_num_options_t;
typedef map<string, string> string_options_t;
typedef map<string, string> date_options_t;
typedef map<string, SymbolList> symbol_list_options_t;
typedef map<string, vector<int> > vec_int_options_t;
num_options_t num_options;
paired_num_options_t paired_num_options;
string_options_t string_options;
date_options_t date_options;
symbol_list_options_t symbol_list_options;
vec_int_options_t vector_int_options;
void writeOutput(ostream &output) const;

View File

@ -1288,6 +1288,7 @@ StaticModel::writeStaticCFile(const string &func_name) const
{
// Writing comments and function definition command
string filename = func_name + "_static.c";
string filename_mex = func_name + "_static_mex.c";
ofstream output;
output.open(filename.c_str(), ios::out | ios::binary);
@ -1303,10 +1304,15 @@ StaticModel::writeStaticCFile(const string &func_name) const
<< " * Warning : this file is generated automatically by Dynare" << endl
<< " * from model file (.mod)" << endl << endl
<< " */" << endl
<< "#include <math.h>" << endl
<< "#include \"mex.h\"" << endl
<< endl
<< "#define max(a, b) (((a) > (b)) ? (a) : (b))" << endl
<< "#include <math.h>" << endl;
if (external_functions_table.get_total_number_of_unique_model_block_external_functions())
// External Matlab function, implies Static function will call mex
output << "#include \"mex.h\"" << endl;
else
output << "#include <stdlib.h>" << endl;
output << "#define max(a, b) (((a) > (b)) ? (a) : (b))" << endl
<< "#define min(a, b) (((a) > (b)) ? (b) : (a))" << endl;
@ -1317,8 +1323,26 @@ StaticModel::writeStaticCFile(const string &func_name) const
writeStaticModel(output, true);
output << "}" << endl << endl;
writePowerDeriv(output, true);
output.close();
output.open(filename_mex.c_str(), ios::out | ios::binary);
if (!output.is_open())
{
cerr << "ERROR: Can't open file " << filename_mex << " for writing" << endl;
exit(EXIT_FAILURE);
}
// Writing the gateway routine
output << "/* The gateway routine */" << endl
output << "/*" << endl
<< " * " << filename_mex << " : The gateway routine used to call the Static function "
<< "located in " << filename << endl
<< " *" << endl
<< " * Warning : this file is generated automatically by Dynare" << endl
<< " * from model file (.mod)" << endl << endl
<< " */" << endl << endl
<< "#include \"mex.h\"" << endl << endl
<< "void Static(double *y, double *x, int nb_row_x, double *params, double *residual, double *g1, double *v2);" << endl
<< "void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])" << endl
<< "{" << endl
<< " double *y, *x, *params;" << endl
@ -1367,7 +1391,6 @@ StaticModel::writeStaticCFile(const string &func_name) const
<< " /* Call the C subroutines. */" << endl
<< " Static(y, x, nb_row_x, params, residual, g1, v2);" << endl
<< "}" << endl << endl;
writePowerDeriv(output, true);
output.close();
}

View File

@ -162,12 +162,9 @@ public:
//! Execute computations (variable sorting + derivation)
/*!
\param jacobianExo whether derivatives w.r. to exo and exo_det should be in the Jacobian (derivatives w.r. to endo are always computed)
\param hessian whether 2nd derivatives w.r. to exo, exo_det and endo should be computed (implies jacobianExo = true)
\param thirdDerivatives whether 3rd derivatives w.r. to endo/exo/exo_det should be computed (implies jacobianExo = true)
\param paramsDerivatives whether 2nd derivatives w.r. to a pair (endo/exo/exo_det, parameter) should be computed (implies jacobianExo = true)
\param eval_context evaluation context for normalization
\param no_tmp_terms if true, no temporary terms will be computed in the static files
\param hessian whether 2nd derivatives w.r. to exo, exo_det and endo should be computed
*/
void computingPass(const eval_context_t &eval_context, bool no_tmp_terms, bool hessian, bool block, bool bytecode);

View File

@ -53,7 +53,6 @@ end;
steady;
options_.ep.verbosity = 0;
options_.ep.stochastic = 0;
options_.console_mode = 0;
ts = extended_path([],1000);

104
tests/ep/rbcii.mod Normal file
View File

@ -0,0 +1,104 @@
@#define extended_path_version = 1
var Capital, Output, Labour, Consumption, Efficiency, efficiency, ExpectedTerm, LagrangeMultiplier;
varexo EfficiencyInnovation;
parameters beta, theta, tau, alpha, psi, delta, rho, effstar, sigma2;
/*
** Calibration
*/
beta = 0.990;
theta = 0.357;
tau = 2.000;
alpha = 0.450;
psi = -0.500;
delta = 0.020;
rho = 0.995;
effstar = 1.000;
sigma2 = 0.001;
@#if extended_path_version
rho = 0.800;
@#endif
external_function(name=mean_preserving_spread);
model(block,bytecode,cutoff=0);
// Eq. n°1:
efficiency = rho*efficiency(-1) + EfficiencyInnovation;
// Eq. n°2:
Efficiency = effstar*exp(efficiency-mean_preserving_spread(rho));
// Eq. n°3:
Output = Efficiency*(alpha*(Capital(-1)^psi)+(1-alpha)*(Labour^psi))^(1/psi);
// Eq. n°4:
Capital = max(Output-Consumption + (1-delta)*Capital(-1),(1-delta)*Capital(-1));
// Eq. n°5:
((1-theta)/theta)*(Consumption/(1-Labour)) - (1-alpha)*(Output/Labour)^(1-psi);
// Eq. n°6:
(((Consumption^theta)*((1-Labour)^(1-theta)))^(1-tau))/Consumption - LagrangeMultiplier - ExpectedTerm(1);
// Eq. n°7:
(Capital==(1-delta)*Capital(-1))*(Output-Consumption) + (1-(Capital==(1-delta)*Capital(-1)))*LagrangeMultiplier = 0;
// Eq. n°8:
ExpectedTerm = beta*(((((Consumption^theta)*((1-Labour)^(1-theta)))^(1-tau))/Consumption)*(alpha*((Output/Capital(-1))^(1-psi))+(1-delta))-(1-delta)*LagrangeMultiplier);
end;
@#if extended_path_version
shocks;
var EfficiencyInnovation = sigma2;
end;
steady;
options_.maxit_ = 100;
options_.ep.verbosity = 0;
options_.ep.stochastic.status = 0;
options_.console_mode = 0;
ts = extended_path([],1000);
options_.ep.stochastic.status = 1;
sts = extended_path([],1000);
figure(1)
plot(ts(2,:)-ts(4,:));
figure(2)
plot(sts(2,:)-sts(4,:));
figure(3)
plot(sts(2,:)-ts(2,:))
@#else
shocks;
var EfficiencyInnovation;
periods 1;
values -.4;
end;
steady;
options_.maxit_ = 100;
simul(periods=4000);
n = 100;
plot(Output(1:n)-Consumption(1:n),'-b','linewidth',2)
@#endif

View File

@ -0,0 +1,49 @@
function [ys, info] = rbcii_steadystate(ys, exogenous)
% Steady state routine for rbc.mod (Business Cycle model with endogenous labour and CES production function)
% AUTHOR(S)
% stephane DOT adjemian AT univ DASH lemans DOT fr
% frederic DOT karame AT univ DASH evry DOT fr
% Output_per_unit_of_Capital = (((1/beta)-1+delta)/alpha)^(1/(1-psi));
% Consumption_per_unit_of_Capital = Output_per_unit_of_Capital - delta;
% Labour_per_unit_of_Capital = (((Output_per_unit_of_Capital/effstar)^psi-alpha)/(1-alpha))^(1/psi);
% Output_per_unit_of_Labour = Output_per_unit_of_Capital/Labour_per_unit_of_Capital;
% Consumption_per_unit_of_Labour = Consumption_per_unit_of_Capital/Labour_per_unit_of_Capital;
% SteadyStateLabour = 1/(1 + Consumption_per_unit_of_Labour/((theta*(1-alpha)/(1-theta))*(Output_per_unit_of_Labour^(1-psi))));
% SteadyStateConsumption = Consumption_per_unit_of_Labour*SteadyStateLabour;
% SteadyStateCapital = SteadyStateLabour/Labour_per_unit_of_Capital;
% SteadyStateOutput = Output_per_unit_of_Capital*SteadyStateCapital;
% ShareOfCapital = alpha/(alpha+(1-alpha)*Labour_per_unit_of_Capital^psi);
global M_
info = 0;
% Compute steady state ratios.
Output_per_unit_of_Capital=((1/M_.params(1)-1+M_.params(6))/M_.params(4))^(1/(1-M_.params(5)));
Consumption_per_unit_of_Capital=Output_per_unit_of_Capital-M_.params(6);
Labour_per_unit_of_Capital=(((Output_per_unit_of_Capital/M_.params(8))^M_.params(5)-M_.params(4))/(1-M_.params(4)))^(1/M_.params(5));
Output_per_unit_of_Labour=Output_per_unit_of_Capital/Labour_per_unit_of_Capital;
Consumption_per_unit_of_Labour=Consumption_per_unit_of_Capital/Labour_per_unit_of_Capital;
% Compute steady state share of capital.
ShareOfCapital=M_.params(4)/(M_.params(4)+(1-M_.params(4))*Labour_per_unit_of_Capital^M_.params(5));
% Compute steady state of the endogenous variables.
SteadyStateLabour=1/(1+Consumption_per_unit_of_Labour/((1-M_.params(4))*M_.params(2)/(1-M_.params(2))*Output_per_unit_of_Labour^(1-M_.params(5))));
SteadyStateConsumption=Consumption_per_unit_of_Labour*SteadyStateLabour;
SteadyStateCapital=SteadyStateLabour/Labour_per_unit_of_Capital;
SteadyStateOutput=Output_per_unit_of_Capital*SteadyStateCapital;
% Fill returned argument ys with steady state values.
ys(2)=SteadyStateOutput;
ys(4)=SteadyStateConsumption;
ys(1)=SteadyStateCapital;
ys(3)=SteadyStateLabour;
ys(5)=M_.params(8);
ys(6)=0;
ys(7)=M_.params(1)*((((SteadyStateConsumption^M_.params(2))*((1-SteadyStateLabour)^(1-M_.params(2))))^(1-M_.params(3)))/SteadyStateConsumption)* ...
(M_.params(4)*((SteadyStateOutput/SteadyStateCapital)^(1-M_.params(5)))+1-M_.params(6));
ys(8)=0;

View File

@ -49,4 +49,7 @@ a(-1) = 0.3;
u(-1) = 0.1;
end;
stoch_simul(nograph);
stoch_simul(nograph, periods = 200);
forecast;

View File

@ -6,7 +6,7 @@ svar_identification;
lower_cholesky;
end;
markov_switching(chain=1,number_of_states=2,duration=2.5);
markov_switching(chain=1,number_of_regimes=2,duration=2.5);
svar(variances, chain=1);

View File

@ -6,7 +6,7 @@ svar_identification;
lower_cholesky;
end;
markov_switching(chain=1,number_of_states=2,duration=2.5);
markov_switching(chain=1,number_of_regimes=2,duration=2.5);
svar(variances, chain=1);

View File

@ -71,6 +71,7 @@ end
failedBlock = {};
num_block_tests = 0;
cd([top_test_dir '/block_bytecode']);
[has_optimization_toolbox junk] = license('checkout','optimization_toolbox');
for blockFlag = 0:1
for bytecodeFlag = 0:1
default_solve_algo = 2;
@ -85,7 +86,7 @@ for blockFlag = 0:1
solve_algos = 1:8;
stack_solve_algos = 0:5;
end
if license('test', 'optimization_toolbox')
if has_optimization_toolbox
solve_algos = [ solve_algos 0 ];
end

View File

@ -57,10 +57,10 @@ Using Dynare with Octave
Dynare is now available for GNU Octave, a free clone of MATLAB (R) (see
<http://www.octave.org>).
The recommended Octave distribution is the Octave/MinGW 3.2.4
The recommended Octave distribution is the Octave/MinGW 3.4.3
precompiled binaries from Octave Forge, available at:
http://sourceforge.net/projects/octave/files/Octave_Windows%20-%20MinGW/
http://sourceforge.net/projects/octave/files/Octave_Windows%20-%20MinGW/Octave%203.4.3%20for%20Windows%20MinGW%20Installer/
Every time you run Octave, you should type the two following commands (assuming
that you have installed Dynare at the standard location, and replacing '4.x.y'

View File

@ -108,7 +108,7 @@ SectionEnd
SectionGroupEnd
Section "MEX files for Octave 3.2.4 (MinGW build)"
Section "MEX files for Octave 3.4.3 (MinGW build)"
SetOutPath $INSTDIR\mex\octave
File ..\mex\octave\*.mex ..\mex\octave\*.oct
SectionEnd