diff --git a/.gitignore b/.gitignore
index bafc55989ccb6a3adb983ccc83912a7b628f5fc2..bd5ffa8c5faed4485bd9dc4640afc00d800be2ac 100644
--- a/.gitignore
+++ b/.gitignore
@@ -87,6 +87,7 @@ mex/build/matlab/run_m2html.m
/preprocessor/dynare_m.exe
/preprocessor/DynareBison.cc
/preprocessor/DynareBison.hh
+/preprocessor/FlexLexer.h
/preprocessor/DynareFlex.cc
/preprocessor/location.hh
/preprocessor/position.hh
diff --git a/doc/bvar-a-la-sims.tex b/doc/bvar-a-la-sims.tex
index ca1621a41b17fa9db8f3996570476a5a0ddd035f..79c7dddb77075d4903780a1e80ef2844420d3baf 100644
--- a/doc/bvar-a-la-sims.tex
+++ b/doc/bvar-a-la-sims.tex
@@ -27,7 +27,7 @@
\author{S\'ebastien Villemot\thanks{Paris School of Economics and
CEPREMAP. E-mail:
\href{mailto:sebastien.villemot@ens.fr}{\texttt{sebastien.villemot@ens.fr}}.}}
-\date{First version: September 2007 \hspace{1cm} This version: September 2011}
+\date{First version: September 2007 \hspace{1cm} This version: August 2012}
\maketitle
@@ -503,7 +503,7 @@ The command will actually compute the marginal density for several models: first
\subsection{Forecasting}
-The syntax for computing forecasts is:
+The syntax for computing (out-of-sample) forecasts is:
\medskip
\texttt{bvar\_forecast(}\textit{options\_list}\texttt{) }\textit{max\_number\_of\_lags}\texttt{;}
diff --git a/doc/dynare.texi b/doc/dynare.texi
index 014b0dd412d3c02c179144584e03b5d256d138f6..9e710937875831967c6c702a12bf9ffa145e2529 100644
--- a/doc/dynare.texi
+++ b/doc/dynare.texi
@@ -4044,7 +4044,7 @@ simulation that crashed prematurely. Shouldn't be used together with
@anchor{mode_file}
Name of the file containing previous value for the mode. When
computing the mode, Dynare stores the mode (@code{xparam1}) and the
-hessian (@code{hh}) in a file called
+hessian (@code{hh}, only if @code{cova_compute=1}) in a file called
@file{@var{MODEL_FILENAME}_mode.mat}
@item mode_compute = @var{INTEGER} | @var{FUNCTION_NAME}
@@ -4284,7 +4284,9 @@ computed after the computation of posterior mode (or maximum
likelihood). This increases speed of computation in large models
during development, when this information is not always necessary. Of
course, it will break all successive computations that would require
-this covariance matrix. Default is @code{1}.
+this covariance matrix. Otherwise, if this option is equal to
+@code{1}, the covariance matrix is computed and stored in variable
+@code{hh} of @file{@var{MODEL_FILENAME}_mode.mat}. Default is @code{1}.
@item solve_algo = @var{INTEGER}
@xref{solve_algo}.
@@ -4685,6 +4687,8 @@ See @file{bvar-a-la-sims.pdf}, which comes with Dynare distribution,
for more information on this command.
@end deffn
+Dynare can also run the smoother on a calibrated model:
+
@deffn Command calib_smoother [@var{VARIABLE_NAME}]@dots{};
@deffnx Command calib_smoother (@var{OPTIONS}@dots{}) [@var{VARIABLE_NAME}]@dots{};
@@ -4985,8 +4989,8 @@ declared in @code{conditional_forecast}.
@end deffn
@deffn Command bvar_forecast ;
-This command computes in-sample or out-sample forecasts for an
-estimated BVAR model, using Minnesota priors.
+This command computes (out-of-sample) forecasts for an estimated BVAR
+model, using Minnesota priors.
See @file{bvar-a-la-sims.pdf}, which comes with Dynare distribution,
for more information on this command.
@@ -5160,6 +5164,10 @@ This command computes an approximation of the optimal policy under
discretion. The algorithm implemented is essentially an LQ solver, and
is described by @cite{Dennis (2007)}.
+You should ensure that your model is linear and your objective is
+quadratic. Also, you should set the @code{linear} option of the
+@code{model} block.
+
@optionshead
This command accepts the same options than @code{ramsey_policy}, plus:
@@ -5188,8 +5196,10 @@ You need to give the one-period objective, not the discounted lifetime
objective. The discount factor is given by the @code{planner_discount}
option of @code{ramsey_policy} and @code{discretionary_policy}.
-Note that with this command you are not limited to quadratic
+With @code{ramsey_policy}, you are not limited to quadratic
objectives: you can give any arbitrary nonlinear expression.
+
+With @code{discretionary_policy}, the objective function must be quadratic.
@end deffn
@node Sensitivity and identification analysis
diff --git a/m4/ax_mexopts.m4 b/m4/ax_mexopts.m4
index a676a79c0113f4fbbc803757e5bc1f56ad57b1f5..69567864838b5abdceb0e32bc0c79e8651f7e777 100644
--- a/m4/ax_mexopts.m4
+++ b/m4/ax_mexopts.m4
@@ -1,4 +1,4 @@
-dnl Copyright (C) 2009-2011 Dynare Team
+dnl Copyright (C) 2009-2012 Dynare Team
dnl
dnl This file is part of Dynare.
dnl
@@ -62,7 +62,7 @@ case ${MATLAB_ARCH} in
ax_mexopts_ok="yes"
;;
maci | maci64)
- SDKROOT='/Developer/SDKs/MacOSX10.6.sdk'
+ SDKROOT='/Applications/Xcode.app/Contents/Developer/Platforms/MacOSX.platform/Developer/SDKs/MacOSX10.6.sdk'
MACOSX_DEPLOYMENT_TARGET='10.6'
if test "${MATLAB_ARCH}" = "maci"; then
ARCHS='i386'
diff --git a/matlab/DsgeVarLikelihood.m b/matlab/DsgeVarLikelihood.m
index 840842cef89a46fcd4534994a5248167574fe75f..c0a0196bcd2efbdb98fbac2d0e6acb2d69f3dc37 100644
--- a/matlab/DsgeVarLikelihood.m
+++ b/matlab/DsgeVarLikelihood.m
@@ -35,7 +35,9 @@ function [fval,grad,hess,exit_flag,info,PHI,SIGMAu,iXX,prior] = DsgeVarLikelihoo
% along with Dynare. If not, see <http://www.gnu.org/licenses/>.
% Declaration of the persistent variables.
-persistent penalty dsge_prior_weight_idx
+persistent dsge_prior_weight_idx
+
+penalty = BayesInfo.penalty;
grad=[];
hess=[];
@@ -46,16 +48,6 @@ SIGMAu = [];
iXX = [];
prior = [];
-% Initialization of the penalty
-if ~nargin || isempty(penalty)
- penalty = 1e8;
- if ~nargin, return, end
-end
-if nargin==1
- penalty = xparam1;
- return
-end
-
% Initialization of of the index for parameter dsge_prior_weight in Model.params.
if isempty(dsge_prior_weight_idx)
dsge_prior_weight_idx = strmatch('dsge_prior_weight',Model.param_names);
@@ -154,7 +146,8 @@ end
[T,R,SteadyState,info,Model,DynareOptions,DynareResults] = dynare_resolve(Model,DynareOptions,DynareResults,'restrict');
% Return, with endogenous penalty when possible, if dynare_resolve issues an error code (defined in resol).
-if info(1) == 1 || info(1) == 2 || info(1) == 5
+if info(1) == 1 || info(1) == 2 || info(1) == 5 || info(1) == 7 || info(1) ...
+ == 8 || info(1) == 22 || info(1) == 24
fval = penalty+1;
info = info(1);
exit_flag = 0;
diff --git a/matlab/cmaes.m b/matlab/cmaes.m
index a2ef613c07588427fa6615fa7bf38c3461d1ee41..c0b6d215c6da3312b283408d981bfa78fc69ff9b 100644
--- a/matlab/cmaes.m
+++ b/matlab/cmaes.m
@@ -644,11 +644,11 @@ else % flgresume
end
% initialize random number generator
- if ischar(opts.Seed)
- randn('state', eval(opts.Seed)); % random number generator state
- else
- randn('state', opts.Seed);
- end
+% $$$ if ischar(opts.Seed)
+% $$$ randn('state', eval(opts.Seed)); % random number generator state
+% $$$ else
+% $$$ randn('state', opts.Seed);
+% $$$ end
%qqq
% load(opts.SaveFilename, 'startseed');
% randn('state', startseed);
diff --git a/matlab/dsge_likelihood.m b/matlab/dsge_likelihood.m
index 0b770a0d4738b1d982c0c7a0194f1d966429e5ca..2a1a9c5a75223e64f69185c4f592e0ffedc08aa0 100644
--- a/matlab/dsge_likelihood.m
+++ b/matlab/dsge_likelihood.m
@@ -130,22 +130,7 @@ function [fval,DLIK,Hess,exit_flag,ys,trend_coeff,info,Model,DynareOptions,Bayes
% AUTHOR(S) stephane DOT adjemian AT univ DASH lemans DOT FR
-% Declaration of the penalty as a persistent variable.
-
-% Persistent variable 'penalty' is used to compute an endogenous penalty to
-% the value 'fval' when various conditions are encountered. These conditions
-% set also 'exit_flag' equal to 0 instead of 1. It is only when
-% dsge_likelihood() is called by an optimizer called by
-% dynare_estimation_1() that 'exit_flag' is ignored and penalized 'fval' is
-% actually used.
-% In that case, 'penalty' is properly initialized, at the very end of the
-% present function, by a call to dsge_likelihood() made in
-% initial_estimation_checks(). If a condition triggers exit_flag ==
-% 0, initial_estimation_checks() triggers an error.
-% In summary, an initial call to the present function, without triggering
-% any condition, guarantees that 'penalty' is properly initialized when needed.
-
-persistent penalty
+penalty = BayesInfo.penalty;
% Initialization of the returned variables and others...
fval = [];
@@ -262,7 +247,8 @@ end
[T,R,SteadyState,info,Model,DynareOptions,DynareResults] = dynare_resolve(Model,DynareOptions,DynareResults,'restrict');
% Return, with endogenous penalty when possible, if dynare_resolve issues an error code (defined in resol).
-if info(1) == 1 || info(1) == 2 || info(1) == 5 || info(1) == 7 || info(1) == 22 || info(1) == 24
+if info(1) == 1 || info(1) == 2 || info(1) == 5 || info(1) == 7 || info(1) ...
+ == 8 || info(1) == 22 || info(1) == 24
fval = penalty+1;
info = info(1);
exit_flag = 0;
@@ -568,16 +554,16 @@ if analytic_derivation,
if full_Hess
DTj = DT(:,:,j+offset);
DPj = dum;
- for i=1:j,
- DTi = DT(:,:,i+offset);
- DPi = DP(:,:,i+offset);
- D2Tij = D2T(:,:,i,j);
- D2Omij = D2Om(:,:,i,j);
+ for i=1:j+offset,
+ DTi = DT(:,:,i);
+ DPi = DP(:,:,i);
+ D2Tij = D2T(:,:,i,j+offset);
+ D2Omij = D2Om(:,:,i,j+offset);
tmp = D2Tij*Pstar*T' + T*Pstar*D2Tij' + DTi*DPj*T' + DTj*DPi*T' + T*DPj*DTi' + T*DPi*DTj' + DTi*Pstar*DTj' + DTj*Pstar*DTi' + D2Omij;
dum = lyapunov_symm(T,tmp,DynareOptions.qz_criterium,DynareOptions.lyapunov_complex_threshold);
dum(abs(dum)<1.e-12) = 0;
- D2P(:,:,i+offset,j+offset) = dum;
- D2P(:,:,j+offset,i+offset) = dum;
+ D2P(:,:,i,j+offset) = dum;
+ D2P(:,:,j+offset,i) = dum;
end
end
end
@@ -771,9 +757,6 @@ end
% Update DynareOptions.kalman_algo.
DynareOptions.kalman_algo = kalman_algo;
-% Update the penalty.
-penalty = fval;
-
if analytic_derivation==0 && nargout==2,
lik=lik(start:end,:);
DLIK=[-lnprior; lik(:)];
diff --git a/matlab/dyn_risky_steadystate_solver.m b/matlab/dyn_risky_steadystate_solver.m
index 6067aa19ebdbab2ac72cac45cdce962a70caf050..d739a8b19927290617c0d2a9afc6baa99a5170a9 100644
--- a/matlab/dyn_risky_steadystate_solver.m
+++ b/matlab/dyn_risky_steadystate_solver.m
@@ -1,92 +1,93 @@
function [dr,info] = dyn_risky_steadystate_solver(ys0,M, ...
dr,options,oo)
-%@info:
-%! @deftypefn {Function File} {[@var{dr},@var{info}] =} dyn_risky_steadystate_solver (@var{ys0},@var{M},@var{dr},@var{options},@var{oo})
-%! @anchor{dyn_risky_steadystate_solver}
-%! @sp 1
-%! Computes the second order risky steady state and first and second order reduced form of the DSGE model.
-%! @sp 2
-%! @strong{Inputs}
-%! @sp 1
-%! @table @ @var
-%! @item ys0
-%! Vector containing a guess value for the risky steady state
-%! @item M
-%! Matlab's structure describing the model (initialized by @code{dynare}).
-%! @item dr
-%! Matlab's structure describing the reduced form solution of the model.
-%! @item options
-%! Matlab's structure describing the options (initialized by @code{dynare}).
-%! @item oo
-%! Matlab's structure gathering the results (initialized by @code{dynare}).
-%! @end table
-%! @sp 2
-%! @strong{Outputs}
-%! @sp 1
-%! @table @ @var
-%! @item dr
-%! Matlab's structure describing the reduced form solution of the model.
-%! @item info
-%! Integer scalar, error code.
-%! @sp 1
-%! @table @ @code
-%! @item info==0
-%! No error.
-%! @item info==1
-%! The model doesn't determine the current variables uniquely.
-%! @item info==2
-%! MJDGGES returned an error code.
-%! @item info==3
-%! Blanchard & Kahn conditions are not satisfied: no stable equilibrium.
-%! @item info==4
-%! Blanchard & Kahn conditions are not satisfied: indeterminacy.
-%! @item info==5
-%! Blanchard & Kahn conditions are not satisfied: indeterminacy due to rank failure.
-%! @item info==6
-%! The jacobian evaluated at the deterministic steady state is complex.
-%! @item info==19
-%! The steadystate routine 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.
-%! @end table
-%! @end table
-%! @end deftypefn
-%@eod:
-
-% Copyright (C) 2001-2012 Dynare Team
-%
-% This file is part of Dynare.
-%
-% Dynare is free software: you can redistribute it and/or modify
-% it under the terms of the GNU General Public License as published by
-% the Free Software Foundation, either version 3 of the License, or
-% (at your option) any later version.
-%
-% Dynare is distributed in the hope that it will be useful,
-% but WITHOUT ANY WARRANTY; without even the implied warranty of
-% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
-% GNU General Public License for more details.
-%
-% You should have received a copy of the GNU General Public License
-% along with Dynare. If not, see <http://www.gnu.org/licenses/>.
+ %@info:
+ %! @deftypefn {Function File} {[@var{dr},@var{info}] =} dyn_risky_steadystate_solver (@var{ys0},@var{M},@var{dr},@var{options},@var{oo})
+ %! @anchor{dyn_risky_steadystate_solver}
+ %! @sp 1
+ %! Computes the second order risky steady state and first and second order reduced form of the DSGE model.
+ %! @sp 2
+ %! @strong{Inputs}
+ %! @sp 1
+ %! @table @ @var
+ %! @item ys0
+ %! Vector containing a guess value for the risky steady state
+ %! @item M
+ %! Matlab's structure describing the model (initialized by @code{dynare}).
+ %! @item dr
+ %! Matlab's structure describing the reduced form solution of the model.
+ %! @item options
+ %! Matlab's structure describing the options (initialized by @code{dynare}).
+ %! @item oo
+ %! Matlab's structure gathering the results (initialized by @code{dynare}).
+ %! @end table
+ %! @sp 2
+ %! @strong{Outputs}
+ %! @sp 1
+ %! @table @ @var
+ %! @item dr
+ %! Matlab's structure describing the reduced form solution of the model.
+ %! @item info
+ %! Integer scalar, error code.
+ %! @sp 1
+ %! @table @ @code
+ %! @item info==0
+ %! No error.
+ %! @item info==1
+ %! The model doesn't determine the current variables uniquely.
+ %! @item info==2
+ %! MJDGGES returned an error code.
+ %! @item info==3
+ %! Blanchard & Kahn conditions are not satisfied: no stable equilibrium.
+ %! @item info==4
+ %! Blanchard & Kahn conditions are not satisfied: indeterminacy.
+ %! @item info==5
+ %! Blanchard & Kahn conditions are not satisfied: indeterminacy due to rank failure.
+ %! @item info==6
+ %! The jacobian evaluated at the deterministic steady state is complex.
+ %! @item info==19
+ %! The steadystate routine 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.
+ %! @end table
+ %! @end table
+ %! @end deftypefn
+ %@eod:
+
+ % Copyright (C) 2001-2012 Dynare Team
+ %
+ % This file is part of Dynare.
+ %
+ % Dynare is free software: you can redistribute it and/or modify
+ % it under the terms of the GNU General Public License as published by
+ % the Free Software Foundation, either version 3 of the License, or
+ % (at your option) any later version.
+ %
+ % Dynare is distributed in the hope that it will be useful,
+ % but WITHOUT ANY WARRANTY; without even the implied warranty of
+ % MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ % GNU General Public License for more details.
+ %
+ % You should have received a copy of the GNU General Public License
+ % along with Dynare. If not, see <http://www.gnu.org/licenses/>.
info = 0;
lead_lag_incidence = M.lead_lag_incidence;
order_var = dr.order_var;
+ endo_nbr = M.endo_nbr;
exo_nbr = M.exo_nbr;
M.var_order_endo_names = M.endo_names(dr.order_var,:);
-
+
[junk,dr.i_fwrd_g,i_fwrd_f] = find(lead_lag_incidence(3,order_var));
dr.i_fwrd_f = i_fwrd_f;
nd = nnz(lead_lag_incidence) + M.exo_nbr;
@@ -104,54 +105,41 @@ function [dr,info] = dyn_risky_steadystate_solver(ys0,M, ...
end
if isfield(options,'portfolio') && options.portfolio == 1
- eq_tags = M.equations_tags;
- n_tags = size(eq_tags,1);
- l_var = zeros(n_tags,1);
- for i=1:n_tags
- l_var(i) = find(strncmp(eq_tags(i,3),M.endo_names, ...
- length(cell2mat(eq_tags(i,3)))));
+ pm = portfolio_model_structure(M,options);
+
+ x0 = ys0(pm.v_p);
+ n = length(x0);
+ [x, info] = solve1(@risky_residuals_ds,x0,1:n,1:n,0,1, options.gstep, ...
+ options.solve_tolf,options.solve_tolx, ...
+ options.solve_maxit,options.debug,pm,M,dr, ...
+ options,oo);
+ if info
+ error('DS approach can''t be computed')
end
- dr.ys = ys0;
- x0 = ys0(l_var);
- % dr = first_step_ds(x0,M,dr,options,oo);
- n = size(ys0);
- %x0 = ys0;
- [x, info] = solve1(@risky_residuals_ds,x0,1:n_tags,1:n_tags,0,1,M,dr,options,oo);
- %[x, info] = solve1(@risky_residuals,x0,1:n,1:n,0,1,M,dr,options,oo);
+ %[x, info] = csolve(@risky_residuals_ds,x0,[],1e-10,100,M,dr,options,oo);
% ys0(l_var) = x;
- ys0(l_var) = x;
- dr.ys = ys0;
- oo.dr = dr;
- oo.steady_state = ys0;
- disp_steady_state(M,oo);
+ [resids,dr1] = risky_residuals_ds(x,pm,M,dr,options,oo);
+ ys1 = dr1.ys;
end
-
- [ys, info] = csolve(func,ys0,[],1e-10,100,M,dr,options,oo);
- if options.k_order_solver
- [resid,dr] = risky_residuals_k_order(ys,M,dr,options,oo);
- else
- [resid,dr] = risky_residuals(ys,M,dr,options,oo);
+ [ys, info] = solve1(func,ys0,1:endo_nbr,1:endo_nbr,0,1, options.gstep, ...
+ options.solve_tolf,options.solve_tolx, ...
+ options.solve_maxit,options.debug,pm,M,dr,options,oo);
+ % [ys, info] = csolve(func,ys0,[],1e-10,100,M,dr,options,oo);
+ if info
+ error('RSS approach can''t be computed')
end
-
dr.ys = ys;
+
+ [resid,dr] = func(ys,pm,M,dr,options,oo);
+
for i=1:M.endo_nbr
- disp(sprintf('%16s %12.6f %12.6f',M.endo_names(i,:),ys0(i), ys(i)))
+ disp(sprintf('%16s %12.6f %12.6f',M.endo_names(i,:),ys1(i), ys(i)))
end
-
- dr.ghs2 = zeros(size(dr.ghs2));
-
- k_var = setdiff(1:M.endo_nbr,l_var);
- dr.ghx(k_var,:) = dr.ghx;
- dr.ghu(k_var,:) = dr.ghu;
- dr.ghxx(k_var,:) = dr.ghxx;
- dr.ghxu(k_var,:) = dr.ghxu;
- dr.ghuu(k_var,:) = dr.ghuu;
- dr.ghs2(k_var,:) = dr.ghs2;
+
end
-function [resid,dr] = risky_residuals(ys,M,dr,options,oo)
- persistent old_ys old_resid
+function [resid,dr] = risky_residuals(ys,pm,M,dr,options,oo)
lead_lag_incidence = M.lead_lag_incidence;
iyv = lead_lag_incidence';
@@ -165,7 +153,7 @@ function [resid,dr] = risky_residuals(ys,M,dr,options,oo)
z = repmat(ys,1,3);
z = z(iyr0) ;
[resid1,d1,d2] = feval([M.fname '_dynamic'],z,...
- [oo.exo_simul ...
+ [oo.exo_simul ...
oo.exo_det_simul], M.params, dr.ys, 2);
if ~isreal(d1) || ~isreal(d2)
pause
@@ -174,150 +162,53 @@ function [resid,dr] = risky_residuals(ys,M,dr,options,oo)
if options.use_dll
% In USE_DLL mode, the hessian is in the 3-column sparse representation
d2 = sparse(d2(:,1), d2(:,2), d2(:,3), ...
- size(d1, 1), size(d1, 2)*size(d1, 2));
+ size(d1, 1), size(d1, 2)*size(d1, 2));
end
if isfield(options,'portfolio') && options.portfolio == 1
- eq_tags = M.equations_tags;
- n_tags = size(eq_tags,1);
- portfolios_eq = cell2mat(eq_tags(strcmp(eq_tags(:,2), ...
- 'portfolio'),1));
- eq = setdiff(1:M.endo_nbr,portfolios_eq);
- l_var = zeros(n_tags,1);
- for i=1:n_tags
- l_var(i) = find(strncmp(eq_tags(i,3),M.endo_names, ...
- length(cell2mat(eq_tags(i,3)))));
- end
- k_var = setdiff(1:M.endo_nbr,l_var);
- lli1 = lead_lag_incidence(:,k_var);
- lead_incidence = lli1(3,:)';
- k = find(lli1');
- lli2 = lli1';
- lli2(k) = 1:nnz(lli1);
- lead_lag_incidence = lli2';
- x = ys(l_var);
- dr = first_step_ds(x,M,dr,options,oo);
-
-
- M.lead_lag_incidence = lead_lag_incidence;
- lli1a = [nonzeros(lli1'); size(d1,2)+(-M.exo_nbr+1:0)'];
- d1a = d1(eq,lli1a);
- ih = 1:size(d2,2);
- ih = reshape(ih,size(d1,2),size(d1,2));
- ih1 = ih(lli1a,lli1a);
- d2a = d2(eq,ih1);
-
- M.endo_nbr = M.endo_nbr-n_tags;
- dr = set_state_space(dr,M,options);
-
- [junk,dr.i_fwrd_g] = find(lead_lag_incidence(3,dr.order_var));
- i_fwrd_f = nonzeros(lead_incidence(dr.order_var));
- i_fwrd2_f = ih(i_fwrd_f,i_fwrd_f);
- dr.i_fwrd_f = i_fwrd_f;
- dr.i_fwrd2_f = i_fwrd2_f;
- nd = nnz(lead_lag_incidence) + M.exo_nbr;
- dr.nd = nd;
- kk = reshape(1:nd^2,nd,nd);
- kkk = reshape(1:nd^3,nd^2,nd);
- dr.i_fwrd2a_f = kk(i_fwrd_f,:);
- % dr.i_fwrd3_f = kkk(i_fwrd2_f,:);
- dr.i_uu = kk(end-M.exo_nbr+1:end,end-M.exo_nbr+1:end);
+ pm = portfolio_model_structure(M,options);
+ x = ys(pm.v_p);
+ dr = first_step_ds(x,pm,M,dr,options,oo);
+ dr.ys = ys;
else
d1a = d1;
d2a = d2;
end
-% $$$ [junk,cols_b,cols_j] = find(lead_lag_incidence(2,dr.order_var));
-% $$$ b = zeros(M.endo_nbr,M.endo_nbr);
-% $$$ b(:,cols_b) = d1a(:,cols_j);
-% $$$
-% $$$ [dr,info] = dyn_first_order_solver(d1a,b,M,dr,options,0);
-% $$$ if info
-% $$$ [m1,m2]=max(abs(ys-old_ys));
-% $$$ disp([m1 m2])
-% $$$ % print_info(info,options.noprint);
-% $$$ resid = old_resid+info(2)/40;
-% $$$ return
-% $$$ end
-% $$$
-% $$$ dr = dyn_second_order_solver(d1a,d2a,dr,M);
-
- gu1 = dr.ghu(dr.i_fwrd_g,:);
-
- resid = resid1+0.5*(d1(:,dr.i_fwrd_f)*dr.ghuu(dr.i_fwrd_g,:)+ ...
- d2(:,dr.i_fwrd2_f)*kron(gu1,gu1))*vec(M.Sigma_e);
- disp(d1(:,dr.i_fwrd_f)*dr.ghuu(dr.i_fwrd_g,:)*vec(M.Sigma_e));
- old_ys = ys;
- disp(max(abs(resid)))
- old_resid = resid;
+ gu1 = dr.ghu(pm.i_fwrd_g,:);
+
+ resid = resid1+0.5*(d1(:,pm.i_fwrd_f1)*dr.ghuu(pm.i_fwrd_g,:)+ ...
+ d2(:,pm.i_fwrd_f2)*kron(gu1,gu1))*vec(M.Sigma_e);
end
-function [resid,dr] = risky_residuals_ds(x,M,dr,options,oo)
- persistent old_ys old_resid old_resid1 old_d1 old_d2
+function [resid,dr] = risky_residuals_ds(x,pm,M,dr,options,oo)
+
+ v_p = pm.v_p;
+ v_np = pm.v_np;
+
+ % computing steady state of non-portfolio variables consistent with
+ % assumed portfolio
+ dr.ys(v_p) = x;
+ ys0 = dr.ys(v_np);
+ f_h =str2func([M.fname '_static']);
+ [dr.ys(v_np),info] = csolve(@ds_static_model,ys0,[],1e-10,100,f_h,x,pm.eq_np,v_np,v_p, ...
+ M.endo_nbr,M.exo_nbr,M.params);
+ if info
+ error('can''t compute non-portfolio steady state')
+ end
- dr = first_step_ds(x,M,dr,options,oo);
+ dr_np = first_step_ds(x,pm,M,dr,options,oo);
lead_lag_incidence = M.lead_lag_incidence;
iyv = lead_lag_incidence';
iyv = iyv(:);
iyr0 = find(iyv) ;
- if M.exo_nbr == 0
- oo.exo_steady_state = [] ;
- end
-
- eq_tags = M.equations_tags;
- n_tags = size(eq_tags,1);
- portfolios_eq = cell2mat(eq_tags(strcmp(eq_tags(:,2), ...
- 'portfolio'),1));
- eq = setdiff(1:M.endo_nbr,portfolios_eq);
- l_var = zeros(n_tags,1);
- for i=1:n_tags
- l_var(i) = find(strncmp(eq_tags(i,3),M.endo_names, ...
- length(cell2mat(eq_tags(i,3)))));
- end
- k_var = setdiff(1:M.endo_nbr,l_var);
- lli1 = lead_lag_incidence(:,k_var);
- k = find(lli1');
- lli2 = lli1';
- lli2(k) = 1:nnz(lli1);
- lead_lag_incidence = lli2';
-
- ys = dr.ys;
- ys(l_var) = x;
-
- z = repmat(ys,1,3);
+ z = repmat(dr.ys,1,3);
z = z(iyr0) ;
[resid1,d1,d2] = feval([M.fname '_dynamic'],z,...
- [oo.exo_simul ...
+ [oo.exo_simul ...
oo.exo_det_simul], M.params, dr.ys, 2);
-% $$$ if isempty(old_resid)
-% $$$ old_resid1 = resid1;
-% $$$ old_d1 = d1;
-% $$$ old_d2 = d2;
-% $$$ old_ys = ys;
-% $$$ else
-% $$$ if ~isequal(resid1,old_resid)
-% $$$ disp('ys')
-% $$$ disp((ys-old_ys)');
-% $$$ disp('resids1')
-% $$$ disp((resid1-old_resid1)')
-% $$$ old_resid1 = resid1;
-% $$$ pause
-% $$$ end
-% $$$ if ~isequal(d1,old_d1)
-% $$$ disp('d1')
-% $$$ disp(d1-old_d1);
-% $$$ old_d1 = d1;
-% $$$ pause
-% $$$ end
-% $$$ if ~isequal(d2,old_d2)
-% $$$ disp('d2')
-% $$$ disp(d2-old_d2);
-% $$$ old_d2 = d2;
-% $$$ pause
-% $$$ end
-% $$$ end
if ~isreal(d1) || ~isreal(d2)
pause
end
@@ -325,95 +216,31 @@ function [resid,dr] = risky_residuals_ds(x,M,dr,options,oo)
if options.use_dll
% In USE_DLL mode, the hessian is in the 3-column sparse representation
d2 = sparse(d2(:,1), d2(:,2), d2(:,3), ...
- size(d1, 1), size(d1, 2)*size(d1, 2));
+ size(d1, 1), size(d1, 2)*size(d1, 2));
end
-% $$$ if isfield(options,'portfolio') && options.portfolio == 1
-% $$$ lli1a = [nonzeros(lli1'); size(d1,2)+(-M.exo_nbr+1:0)'];
-% $$$ d1a = d1(eq,lli1a);
-% $$$ ih = 1:size(d2,2);
-% $$$ ih = reshape(ih,size(d1,2),size(d1,2));
-% $$$ ih1 = ih(lli1a,lli1a);
-% $$$ d2a = d2(eq,ih1);
-% $$$
-% $$$ M.endo_nbr = M.endo_nbr-n_tags;
-% $$$ dr = set_state_space(dr,M);
-% $$$
-% $$$ dr.i_fwrd_g = find(lead_lag_incidence(3,dr.order_var)');
-% $$$ else
-% $$$ d1a = d1;
-% $$$ d2a = d2;
-% $$$ end
-% $$$
-% $$$ [junk,cols_b,cols_j] = find(lead_lag_incidence(2,dr.order_var));
-% $$$ b = zeros(M.endo_nbr,M.endo_nbr);
-% $$$ b(:,cols_b) = d1a(:,cols_j);
-% $$$
-% $$$ [dr,info] = dyn_first_order_solver(d1a,b,M,dr,options,0);
-% $$$ if info
-% $$$ [m1,m2]=max(abs(ys-old_ys));
-% $$$ disp([m1 m2])
-% $$$ % print_info(info,options.noprint);
-% $$$ resid = old_resid+info(2)/40;
-% $$$ return
-% $$$ end
-% $$$
-% $$$ dr = dyn_second_order_solver(d1a,d2a,dr,M);
- gu1 = dr.ghu(dr.i_fwrd_g,:);
-
- % resid = resid1+0.5*(d1(:,dr.i_fwrd_f)*dr.ghuu(dr.i_fwrd_g,:)+ ...
- % d2(:,dr.i_fwrd2_f)*kron(gu1,gu1))*vec(M.Sigma_e);
- resid = resid1+0.5*(d2(:,dr.i_fwrd2_f)*kron(gu1,gu1))*vec(M.Sigma_e);
-
-% $$$ if isempty(old_resid)
-% $$$ old_resid = resid;
-% $$$ else
-% $$$ disp('resid')
-% $$$ dr = (resid-old_resid)';
-% $$$ % disp(dr)
-% $$$ % disp(dr(portfolios_eq))
-% $$$ old_resid = resid;
-% $$$ end
- resid = resid(portfolios_eq)
+ gu1 = dr_np.ghu(pm.i_fwrd_g,:);
+
+ resid = resid1+0.5*(d2(:,pm.i_fwrd_f2)*kron(gu1,gu1))*vec(M.Sigma_e);
+
+ resid = resid(pm.eq_p)
end
-function [dr] = first_step_ds(x,M,dr,options,oo)
-
+function dr_np = first_step_ds(x,pm,M,dr,options,oo)
+
lead_lag_incidence = M.lead_lag_incidence;
iyv = lead_lag_incidence';
iyv = iyv(:);
iyr0 = find(iyv) ;
-
- if M.exo_nbr == 0
- oo.exo_steady_state = [] ;
- end
-
- eq_tags = M.equations_tags;
- n_tags = size(eq_tags,1);
- portfolios_eq = cell2mat(eq_tags(strcmp(eq_tags(:,2), ...
- 'portfolio'),1));
- eq = setdiff(1:M.endo_nbr,portfolios_eq);
- l_var = zeros(n_tags,1);
- for i=1:n_tags
- l_var(i) = find(strncmp(eq_tags(i,3),M.endo_names, ...
- length(cell2mat(eq_tags(i,3)))));
- end
- k_var = setdiff(1:M.endo_nbr,l_var);
- lli1 = lead_lag_incidence(:,k_var);
- k = find(lli1');
- lli2 = lli1';
- lli2(k) = 1:nnz(lli1);
- lead_lag_incidence = lli2';
- M.lead_lag_incidence = lead_lag_incidence;
ys = dr.ys;
- ys(l_var) = x;
+ ys(pm.v_p) = x;
z = repmat(ys,1,3);
z = z(iyr0) ;
[resid1,d1,d2] = feval([M.fname '_dynamic'],z,...
- [oo.exo_simul ...
+ [oo.exo_simul ...
oo.exo_det_simul], M.params, dr.ys, 2);
if ~isreal(d1) || ~isreal(d2)
pause
@@ -422,102 +249,139 @@ function [dr] = first_step_ds(x,M,dr,options,oo)
if options.use_dll
% In USE_DLL mode, the hessian is in the 3-column sparse representation
d2 = sparse(d2(:,1), d2(:,2), d2(:,3), ...
- size(d1, 1), size(d1, 2)*size(d1, 2));
+ size(d1, 1), size(d1, 2)*size(d1, 2));
end
- if isfield(options,'portfolio') && options.portfolio == 1
- lli1a = [nonzeros(lli1'); size(d1,2)+(-M.exo_nbr+1:0)'];
- d1a = d1(eq,lli1a);
- ih = 1:size(d2,2);
- ih = reshape(ih,size(d1,2),size(d1,2));
- ih1 = ih(lli1a,lli1a);
- d2a = d2(eq,ih1);
-
- M.endo_nbr = M.endo_nbr-n_tags;
- dr = set_state_space(dr,M,options);
+ d1_np = d1(pm.eq_np,pm.i_d1_np);
+ d2_np = d2(pm.eq_np,pm.i_d2_np);
- dr.i_fwrd_g = find(lead_lag_incidence(3,dr.order_var)');
- else
- d1a = d1;
- d2a = d2;
- end
-
- [junk,cols_b,cols_j] = find(lead_lag_incidence(2,dr.order_var));
- b = zeros(M.endo_nbr,M.endo_nbr);
- b(:,cols_b) = d1a(:,cols_j);
-
- [dr,info] = dyn_first_order_solver(d1a,M,dr,options,0);
+ [dr_np,info] = dyn_first_order_solver(d1_np,pm.M_np,pm.dr_np,options,0);
if info
- [m1,m2]=max(abs(ys-old_ys));
- disp([m1 m2])
- % print_info(info,options.noprint);
- resid = old_resid+info(2)/40;
+ print_info(info,0);
return
end
- dr = dyn_second_order_solver(d1a,d2a,dr,M,...
- options.threads.kronecker.A_times_B_kronecker_C,...
- options.threads.kronecker.sparse_hessian_times_B_kronecker_C);
+ dr_np = dyn_second_order_solver(d1_np,d2_np,dr_np,pm.M_np,...
+ options.threads.kronecker.A_times_B_kronecker_C,...
+ options.threads.kronecker.sparse_hessian_times_B_kronecker_C);
end
-function [resid,dr] = risky_residuals_k_order(ys,M,dr,options,oo)
-
- lead_lag_incidence = M.lead_lag_incidence;
- npred = dr.npred;
+function [resid,dr] = risky_residuals_k_order(ys,pm,M,dr,options,oo)
exo_nbr = M.exo_nbr;
- vSigma_e = vec(M.Sigma_e);
+ endo_nbr = M.endo_nbr;
- iyv = lead_lag_incidence';
+ iyv = M.lead_lag_incidence';
iyv = iyv(:);
iyr0 = find(iyv) ;
- if M.exo_nbr == 0
+ if exo_nbr == 0
oo.exo_steady_state = [] ;
end
z = repmat(ys,1,3);
z = z(iyr0) ;
- [resid1,d1,d2,d3] = feval([M.fname '_dynamic'],z,...
- [oo.exo_simul ...
+ [resid1,d1,d2] = feval([M.fname '_dynamic'],z,...
+ [oo.exo_simul ...
oo.exo_det_simul], M.params, dr.ys, 2);
-
-% $$$ hessian = sparse(d2(:,1), d2(:,2), d2(:,3), ...
-% $$$ size(d1, 1), size(d1, 2)*size(d1, 2));
-% $$$ fy3 = sparse(d2(:,1), d2(:,2), d2(:,3), ...
-% $$$ size(d1, 1), size(d1, 2)^3);
-
- options.order = 3;
- nu2 = exo_nbr*(exo_nbr+1)/2;
-% $$$ d1_0 = d1;
-% $$$ gu1 = dr.ghu(dr.i_fwrd_g,:);
-% $$$ guu = dr.ghuu;
-% $$$ for i=1:2
-% $$$ d1 = d1_0 + 0.5*(hessian(:,dr.i_fwrd2a_f)*kron(eye(dr.nd),guu(dr.i_fwrd_g,:)*vSigma_e)+ ...
-% $$$ fy3(:,dr.i_fwrd3_f)*kron(eye(dr.nd),kron(gu1,gu1)*vSigma_e));
-% $$$ [junk,cols_b,cols_j] = find(lead_lag_incidence(2,dr.order_var));
-% $$$ b = zeros(M.endo_nbr,M.endo_nbr);
-% $$$ b(:,cols_b) = d1(:,cols_j);
-
-% $$$ [dr,info] = dyn_first_order_solver(d1,b,M,dr,options,0);
- [err,g_0, g_1, g_2, g_3] = k_order_perturbation(dr,M,options);
- mexErrCheck('k_order_perturbation', err);
- gu1 = g_1(dr.i_fwrd_g,end-exo_nbr+1:end);
- guu = unfold(g_2(:,end-nu2+1:end),exo_nbr);
- d1old = d1;
- % disp(max(max(abs(d1-d1old))));
- % end
+ if isfield(options,'portfolio') && options.portfolio == 1
+ eq_np = pm.eq_np;
+
+ d1_np = d1(eq_np,pm.i_d1_np);
+ d2_np = d2(eq_np,pm.i_d2_np);
+
+ M_np = pm.M_np;
+ dr_np = pm.dr_np;
+
+ [dr_np,info] = dyn_first_order_solver(d1_np,pm.M_np,pm.dr_np,options,0);
+ if info
+ print_info(info,0);
+ return
+ end
+
+ dr_np = dyn_second_order_solver(d1_np,d2_np,dr_np,pm.M_np,...
+ options.threads.kronecker.A_times_B_kronecker_C,...
+ options.threads.kronecker.sparse_hessian_times_B_kronecker_C);
+ end
- [junk,cols_b,cols_j] = find(lead_lag_incidence(2,dr.order_var));
+ i_fwrd_f1 = pm.i_fwrd_f1;
+ i_fwrd_f2 = pm.i_fwrd_f2;
+ i_fwrd_f3 = pm.i_fwrd_f3;
+ i_fwrd_g = pm.i_fwrd_g;
+ gu1 = dr_np.ghu(i_fwrd_g,:);
+ ghuu = dr_np.ghuu;
- resid = resid1+0.5*(d1(:,dr.i_fwrd_f)*guu(dr.i_fwrd_g,:)+hessian(:,dr.i_fwrd2_f)*kron(gu1,gu1))*vec(M.Sigma_e);
+ resid = resid1+0.5*(d1(:,i_fwrd_f1)*ghuu(i_fwrd_g,:)+d2(:,i_fwrd_f2)* ...
+ kron(gu1,gu1))*vec(M.Sigma_e);
if nargout > 1
- [dr,info] = k_order_pert(dr,M,options,oo);
+ [resid1,d1,d2,d3] = feval([M.fname '_dynamic'],z,...
+ [oo.exo_simul ...
+ oo.exo_det_simul], M.params, dr.ys, 2);
+
+
+ [a,b,c] = find(d2(eq_np,pm.i_d2_np));
+ d2_np = [a b c];
+
+ [a,b,c] = find(d3(eq_np,pm.i_d3_np));
+ d3_np = [a b c];
+
+ options.order = 3;
+ % space holder, unused by k_order_pertrubation
+ dr_np.ys = dr.ys(pm.v_np);
+ nu2 = exo_nbr*(exo_nbr+1)/2;
+ nu3 = exo_nbr*(exo_nbr+1)*(exo_nbr+2)/3;
+ M_np.NZZDerivatives = [nnz(d1_np); nnz(d2_np); nnz(d3_np)];
+ [err,g_0, g_1, g_2, g_3] = k_order_perturbation(dr_np,M_np,options,d1_np,d2_np,d3_np);
+ mexErrCheck('k_order_perturbation', err);
+
+ gu1 = g_1(i_fwrd_g,end-exo_nbr+1:end);
+ ghuu = unfold2(g_2(:,end-nu2+1:end),exo_nbr);
+ ghsuu = get_ghsuu(g_3,size(g_1,2),exo_nbr);
+
+ i_fwrd1_f2 = pm.i_fwrd1_f2;
+ i_fwrd1_f3 = pm.i_fwrd1_f3;
+ n = size(d1,2);
+ d1b = d1 + 0.5*( ...
+ d1(:,i_fwrd_f1)*...
+ d2(:,i_fwrd1_f2)*kron(eye(n),dr_np.ghuu(i_fwrd_g,:)*vec(M.Sigma_e))...
+ + 0.5*d3(:,i_fwrd1_f3)*kron(eye(n),kron(gu1,gu1)*vec(M.Sigma_e)));
+ format short
+ kk1 = [nonzeros(M.lead_lag_incidence(:,1:6)'); ...
+ nnz(M.lead_lag_incidence)+[1; 2]]
+ kk2 = [nonzeros(M.lead_lag_incidence(:,1:6)'); ...
+ nnz(M.lead_lag_incidence)+[3; 4]]
+ format short
+ gu1
+ kron(gu1,gu1)*vec(M.Sigma_e)
+ disp(d1(:,:))
+ disp(d1b(:,:))
+ aa2=d2(:,i_fwrd1_f2)*kron(eye(n),dr_np.ghuu(i_fwrd_g,:)*vec(M.Sigma_e));
+ aa3=d3(:,i_fwrd1_f3)*kron(eye(n),kron(gu1,gu1)*vec(M.Sigma_e));
+ disp(d3(4,7+6*n+6*n*n))
+ disp(d3(4,8+16*n+17*n*n)) %8,17,18
+ disp(d3(4,8+17*n+16*n*n)) %8,17,18
+ disp(d3(4,7*n+17+17*n*n)) %8,17,18
+ disp(d3(4,7*n+18+16*n*n)) %8,17,18
+ disp(d3(4,7*n*n+16*n+18)) %8,17,18
+ disp(d3(4,7*n*n+17+17*n)) %8,17,18
+ pause
+ disp(aa2(:,kk1))
+ disp(aa2(:,kk2))
+ disp(aa3(:,kk1))
+ disp(aa3(:,kk2))
+ [dr,info] = dyn_first_order_solver(d1b,M,dr,options,0);
+ if info
+ print_info(info,0);
+ return
+ end
+
+ disp_dr(dr,dr.order_var,[]);
+
end
end
-function y=unfold(x,n)
+function y=unfold2(x,n)
y = zeros(size(x,1),n*n);
k = 1;
for i=1:n
@@ -526,6 +390,122 @@ function y=unfold(x,n)
if i ~= j
y(:,(j-1)*n+i) = x(:,k);
end
+ k = k+1;
end
end
end
+
+function y=unfold3(x,n)
+ y = zeros(size(x,1),n*n*n);
+ k = 1;
+ for i=1:n
+ for j=i:n
+ for m=j:n
+ y(:,(i-1)*n*n+(j-1)*n+m) = x(:,k);
+ y(:,(i-1)*n*n+(m-1)*n+j) = x(:,k);
+ y(:,(j-1)*n*n+(i-1)*n+m) = x(:,k);
+ y(:,(j-1)*n*n+(m-1)*n+i) = x(:,k);
+ y(:,(m-1)*n*n+(i-1)*n+j) = x(:,k);
+ y(:,(m-1)*n*n+(j-1)*n+i) = x(:,k);
+
+ k = k+1;
+ end
+ end
+ end
+end
+
+function pm = portfolio_model_structure(M,options)
+
+ i_d3_np = [];
+ i_d3_p = [];
+
+ lead_index = M.maximum_endo_lag+2;
+ lead_lag_incidence = M.lead_lag_incidence;
+ eq_tags = M.equations_tags;
+ n_tags = size(eq_tags,1);
+ eq_p = cell2mat(eq_tags(strcmp(eq_tags(:,2), ...
+ 'portfolio'),1));
+ pm.eq_p = eq_p;
+ pm.eq_np = setdiff(1:M.endo_nbr,eq_p);
+ v_p = zeros(n_tags,1);
+ for i=1:n_tags
+ v_p(i) = find(strncmp(eq_tags(i,3),M.endo_names, ...
+ length(cell2mat(eq_tags(i,3)))));
+ end
+ if any(lead_lag_incidence(lead_index,v_p))
+ error(['portfolio variables appear in the model as forward ' ...
+ 'variable'])
+ end
+ pm.v_p = v_p;
+ v_np = setdiff(1:M.endo_nbr,v_p);
+ pm.v_np = v_np;
+ lli_np = lead_lag_incidence(:,v_np)';
+ k = find(lli_np);
+ lead_lag_incidence_np = lli_np;
+ lead_lag_incidence_np(k) = 1:nnz(lli_np);
+ lead_lag_incidence_np = lead_lag_incidence_np';
+ pm.lead_lag_incidence_np = lead_lag_incidence_np;
+ i_d1_np = [nonzeros(lli_np); nnz(lead_lag_incidence)+(1:M.exo_nbr)'];
+ pm.i_d1_np = i_d1_np;
+
+ n = nnz(lead_lag_incidence)+M.exo_nbr;
+ ih = reshape(1:n*n,n,n);
+ i_d2_np = ih(i_d1_np,i_d1_np);
+ pm.i_d2_np = i_d2_np(:);
+
+ ih = reshape(1:n*n*n,n,n,n);
+ i_d3_np = ih(i_d1_np,i_d1_np,i_d1_np);
+ pm.i_d3_np = i_d3_np(:);
+
+ M_np = M;
+ M_np.lead_lag_incidence = lead_lag_incidence_np;
+ M_np.lead_lag_incidence = lead_lag_incidence_np;
+ M_np.endo_nbr = length(v_np);
+ M_np.endo_names = M.endo_names(v_np,:);
+ dr_np = struct();
+ dr_np = set_state_space(dr_np,M_np,options);
+ pm.dr_np = dr_np;
+ M_np.var_order_endo_names = M_np.endo_names(dr_np.order_var,:);
+ pm.M_np = M_np;
+ pm.i_fwrd_g = find(lead_lag_incidence_np(lead_index,dr_np.order_var)');
+
+ i_fwrd_f1 = nonzeros(lead_lag_incidence(lead_index,:));
+ pm.i_fwrd_f1 = i_fwrd_f1;
+ n = nnz(lead_lag_incidence)+M.exo_nbr;
+ ih = reshape(1:n*n,n,n);
+ i_fwrd_f2 = ih(i_fwrd_f1,i_fwrd_f1);
+ pm.i_fwrd_f2 = i_fwrd_f2(:);
+ i_fwrd1_f2 = ih(i_fwrd_f1,:);
+ pm.i_fwrd1_f2 = i_fwrd1_f2(:);
+
+ ih = reshape(1:n*n*n,n,n,n);
+ i_fwrd_f3 = ih(i_fwrd_f1,i_fwrd_f1,i_fwrd_f1);
+ pm.i_fwrd_f3 = i_fwrd_f3(:);
+ i_fwrd1_f3 = ih(i_fwrd_f1,i_fwrd_f1,:);
+ pm.i_fwrd1_f3 = i_fwrd1_f3(:);
+end
+
+function r=ds_static_model(y0,f_h,p0,eq_np,v_np,v_p,endo_nbr,exo_nbr,params)
+ ys = zeros(endo_nbr,1);
+ ys(v_p) = p0;
+ ys(v_np) = y0;
+ r = f_h(ys,zeros(exo_nbr,1),params);
+ r = r(eq_np);
+end
+
+function ghsuu = get_ghsuu(g,ns,nx)
+ nxx = nx*(nx+1)/2;
+ m1 = 0;
+ m2 = ns*(ns+1)/2;
+ kk = 1:(nx*nx);
+ ghsuu = zeros(size(g,1),(ns*nx*nx));
+
+ for i=1:n
+ j = m1+(1:m2);
+ k = j(end-nxx+1:end);
+ ghsuu(:,kk) = unfold2(g(:,k),nx);
+ m1 = m1+m2;
+ m2 = m2 - (n-i+1);
+ kk = kk + nx*nx;
+ end
+end
\ No newline at end of file
diff --git a/matlab/dynare_config.m b/matlab/dynare_config.m
index 0644f8e7f7afefb437902c9eb232c1b8066bbe3e..c5963cc36300355a670bb5c2e4fc70131d487e5b 100644
--- a/matlab/dynare_config.m
+++ b/matlab/dynare_config.m
@@ -66,8 +66,9 @@ addpath([dynareroot '/utilities/general/'])
% or some MATLAB versions, and for which we provide some replacement functions
if ~exist('OCTAVE_VERSION')
- % Replacements for rows() and columns() (inexistent under MATLAB)
+ % Replacements for rows(), columns() and issquare() (inexistent under MATLAB)
addpath([dynareroot '/missing/rows_columns'])
+ addpath([dynareroot '/missing/issquare'])
% Replacement for vec() (inexistent under MATLAB)
addpath([dynareroot '/missing/vec'])
if ~user_has_matlab_license('statistics_toolbox')
@@ -93,11 +94,6 @@ if (exist('OCTAVE_VERSION') && ~user_has_octave_forge_package('statistics')) ...
addpath([dynareroot '/missing/nanmean'])
end
-% linsolve is missing in Octave
-if (exist('OCTAVE_VERSION'))
- addpath([dynareroot '/missing/linsolve'])
-end
-
% Add path to MEX files
if exist('OCTAVE_VERSION')
addpath([dynareroot '../mex/octave/']);
diff --git a/matlab/dynare_solve.m b/matlab/dynare_solve.m
index 72555ed859c44dc0d1741e6ef574244fd600a00a..da0e58921ada2e5c4372ffdfd133ca70cfff2669 100644
--- a/matlab/dynare_solve.m
+++ b/matlab/dynare_solve.m
@@ -74,7 +74,9 @@ if options_.solve_algo == 0
end
elseif options_.solve_algo == 1
nn = size(x,1);
- [x,info]=solve1(func,x,1:nn,1:nn,jacobian_flag,1,varargin{:});
+ [x,info]=solve1(func,x,1:nn,1:nn,jacobian_flag,1,options_.gstep, ...
+ options_.solve_tolf,options_.solve_tolx, ...
+ options_.solve_maxit,options_.debug,varargin{:});
elseif options_.solve_algo == 2 || options_.solve_algo == 4
nn = size(x,1) ;
tolf = options_.solve_tolf ;
@@ -125,14 +127,19 @@ elseif options_.solve_algo == 2 || options_.solve_algo == 4
if options_.debug
disp(['DYNARE_SOLVE (solve_algo=2|4): solving block ' num2str(i) ', of size ' num2str(r(i+1)-r(i)) ]);
end
- [x,info]=solve1(func,x,j1(r(i):r(i+1)-1),j2(r(i):r(i+1)-1),jacobian_flag, bad_cond_flag, varargin{:});
+ [x,info]=solve1(func,x,j1(r(i):r(i+1)-1),j2(r(i):r(i+1)-1),jacobian_flag, ...
+ bad_cond_flag, options_.gstep, ...
+ options_.solve_tolf,options_.solve_tolx, ...
+ options_.solve_maxit,options_.debug,varargin{:});
if info
return
end
end
fvec = feval(func,x,varargin{:});
if max(abs(fvec)) > tolf
- [x,info]=solve1(func,x,1:nn,1:nn,jacobian_flag, bad_cond_flag, varargin{:});
+ [x,info]=solve1(func,x,1:nn,1:nn,jacobian_flag, bad_cond_flag, ...
+ options_.gstep, options_.solve_tolf,options_.solve_tolx, ...
+ options_.solve_maxit,options_.debug,varargin{:});
end
elseif options_.solve_algo == 3
if jacobian_flag
diff --git a/matlab/metropolis_hastings_initialization.m b/matlab/metropolis_hastings_initialization.m
index 80c2f3532a809cd39a3c376ab98e134a5775a37c..5b5f255626ac06ebc73eaf592f2eef52b14e36cd 100644
--- a/matlab/metropolis_hastings_initialization.m
+++ b/matlab/metropolis_hastings_initialization.m
@@ -1,5 +1,7 @@
function [ ix2, ilogpo2, ModelName, MhDirectoryName, fblck, fline, npar, nblck, nruns, NewFile, MAX_nruns, d ] = ...
- metropolis_hastings_initialization(TargetFun, xparam1, vv, mh_bounds, varargin)
+ metropolis_hastings_initialization(TargetFun, xparam1, vv, mh_bounds,dataset_,options_,M_,estim_params_,bayestopt_,oo_)
+%function [ ix2, ilogpo2, ModelName, MhDirectoryName, fblck, fline, npar, nblck, nruns, NewFile, MAX_nruns, d ] =
+% metropolis_hastings_initialization(TargetFun, xparam1, vv, mh_bounds, dataset_,options_,M_,estim_params_,bayestopt_,oo_)
% Metropolis-Hastings initialization.
%
% INPUTS
@@ -8,7 +10,12 @@ function [ ix2, ilogpo2, ModelName, MhDirectoryName, fblck, fline, npar, nblck,
% o xparam1 [double] (p*1) vector of parameters to be estimated (initial values).
% o vv [double] (p*p) matrix, posterior covariance matrix (at the mode).
% o mh_bounds [double] (p*2) matrix defining lower and upper bounds for the parameters.
-% o varargin list of argument following mh_bounds
+% o dataset_ data structure
+% o options_ options structure
+% o M_ model structure
+% o estim_params_ estimated parameters structure
+% o bayestopt_ estimation options structure
+% o oo_ outputs structure
%
% OUTPUTS
% None
@@ -33,8 +40,6 @@ function [ ix2, ilogpo2, ModelName, MhDirectoryName, fblck, fline, npar, nblck,
% You should have received a copy of the GNU General Public License
% along with Dynare. If not, see <http://www.gnu.org/licenses/>.
-global M_ options_ bayestopt_
-
ix2 = [];
ilogpo2 = [];
ModelName = [];
@@ -54,8 +59,6 @@ if ~isempty(M_.bvar)
ModelName = [M_.fname '_bvar'];
end
-bayestopt_.penalty = 1e8;
-
MhDirectoryName = CheckPath('metropolis',M_.dname);
nblck = options_.mh_nblck;
@@ -107,8 +110,9 @@ if ~options_.load_mh_file && ~options_.mh_recover
candidate = rand_multivariate_normal( transpose(xparam1), d * options_.mh_init_scale, npar);
if all(candidate(:) > mh_bounds(:,1)) && all(candidate(:) < mh_bounds(:,2))
ix2(j,:) = candidate;
- ilogpo2(j) = - feval(TargetFun,ix2(j,:)',varargin{:});
- if ilogpo2(j) <= - bayestopt_.penalty+1e-6
+ ilogpo2(j) = - feval(TargetFun,ix2(j,:)',dataset_,options_,M_,estim_params_,bayestopt_,oo_);
+ if ~isfinite(ilogpo2(j)) % if returned log-density is
+ % Inf or Nan (penalized value)
validate = 0;
else
fprintf(fidlog,[' Blck ' int2str(j) ':\n']);
@@ -142,7 +146,7 @@ if ~options_.load_mh_file && ~options_.mh_recover
candidate = transpose(xparam1);
if all(candidate' > mh_bounds(:,1)) && all(candidate' < mh_bounds(:,2))
ix2 = candidate;
- ilogpo2 = - feval(TargetFun,ix2',varargin{:});
+ ilogpo2 = - feval(TargetFun,ix2',dataset_,options_,M_,estim_params_,bayestopt_,oo_);
disp('MH: Initialization at the posterior mode.')
disp(' ')
fprintf(fidlog,[' Blck ' int2str(1) 'params:\n']);
@@ -180,9 +184,6 @@ if ~options_.load_mh_file && ~options_.mh_recover
% separate initializaton for each chain
JSUM = 0;
for j=1:nblck,
- JSUM = JSUM + sum(100*clock);
- randn('state',JSUM);
- rand('state',JSUM);
record.Seeds(j).Normal = randn('state');
record.Seeds(j).Unifor = rand('state');
end
diff --git a/matlab/issquare.m b/matlab/missing/issquare/issquare.m
similarity index 87%
rename from matlab/issquare.m
rename to matlab/missing/issquare/issquare.m
index 71233cbf6a4d93fc1b6190efebbf585ebfdfec07..4d4eb2b572c8a771080d57d35f8ffdc0904131af 100644
--- a/matlab/issquare.m
+++ b/matlab/missing/issquare/issquare.m
@@ -4,7 +4,7 @@ function i = issquare(A)
%! @deftypefn {Function File} {@var{i} =} issquare (@var{A})
%! @anchor{issquare}
%! @sp 1
-%! Returns 1 if @var{A} is a square matrix, 0 otherwise.
+%! If @var{A} is a square matrix, returns its dimension; otherwise return 0.
%! @sp 2
%! @strong{Inputs}
%! @sp 1
@@ -17,7 +17,7 @@ function i = issquare(A)
%! @sp 1
%! @table @ @var
%! @item i
-%! Integer scalar (0 or 1).
+%! Integer scalar.
%! @end table
%! @sp 2
%! @strong{This function is called by:}
@@ -45,4 +45,8 @@ function i = issquare(A)
% along with Dynare. If not, see <http://www.gnu.org/licenses/>.
d = size(A);
-i = (length(d)==2) && (d(1)==d(2));
\ No newline at end of file
+if (length(d)==2) && (d(1)==d(2))
+ i = d(1);
+else
+ i = 0;
+end
diff --git a/matlab/missing/linsolve/linsolve.m b/matlab/missing/linsolve/linsolve.m
deleted file mode 100644
index 2ae167bd4455b97020a3ae46ee3d113f5ea0e238..0000000000000000000000000000000000000000
--- a/matlab/missing/linsolve/linsolve.m
+++ /dev/null
@@ -1,33 +0,0 @@
-function [x,c] = linsolve(A,B,opts)
-% (very imperfect) Clone of Matlab's linsolve.
-
-% 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/>.
-
- c = [];
- x = [];
- if nargin == 3
- if isfield(opts,'TRANSA')
- A = A';
- end
- end
- if nargout == 2
- c = rcond(A);
- end
-
- x = A\B;
-
diff --git a/matlab/print_info.m b/matlab/print_info.m
index b406e964a004111787713a13f76ca84c7c827925..15a793614fe22d181594d8ed8cc4380b175e69c6 100644
--- a/matlab/print_info.m
+++ b/matlab/print_info.m
@@ -51,7 +51,20 @@ if ~noprint
case 7
error(['One of the eigenvalues is close to 0/0 (the absolute ' ...
'value of numerator and denominator is smaller than 1e-6)'])
- case 19
+ case 8
+ if ~isempty(info(2))
+ global M_;
+ disp_string=deblank(M_.param_names(info(2),:));
+ for ii=1:length(info)-2
+ disp_string=[disp_string,', ',deblank(M_.param_names(info(2+ii),:))];
+ end
+ error(['The Jacobian contains NaNs because the following parameters are NaN: '...
+ disp_string])
+ else
+ error(['The Jacobian contains NaNs'])
+ end
+
+ case 19
error('The steadystate file did not compute the steady state')
case 20
error(['Impossible to find the steady state. Either the model' ...
diff --git a/matlab/random_walk_metropolis_hastings.m b/matlab/random_walk_metropolis_hastings.m
index 48c20fa947aa9eb1c158859b980ffc93a32de0d3..6e5d782d4428e6983d73bd36524dc99f4a6df484 100644
--- a/matlab/random_walk_metropolis_hastings.m
+++ b/matlab/random_walk_metropolis_hastings.m
@@ -1,5 +1,5 @@
-function record=random_walk_metropolis_hastings(TargetFun,ProposalFun,xparam1,vv,mh_bounds,varargin)
-%function record=random_walk_metropolis_hastings(TargetFun,ProposalFun,xparam1,vv,mh_bounds,varargin)
+function record=random_walk_metropolis_hastings(TargetFun,ProposalFun,xparam1,vv,mh_bounds,dataset_,options_,M_,estim_params_,bayestopt_,oo_)
+%function record=random_walk_metropolis_hastings(TargetFun,ProposalFun,xparam1,vv,mh_bounds,dataset_,options_,M_,estim_params_,bayestopt_,oo_)
% Random walk Metropolis-Hastings algorithm.
%
% INPUTS
@@ -8,7 +8,12 @@ function record=random_walk_metropolis_hastings(TargetFun,ProposalFun,xparam1,vv
% o xparam1 [double] (p*1) vector of parameters to be estimated (initial values).
% o vv [double] (p*p) matrix, posterior covariance matrix (at the mode).
% o mh_bounds [double] (p*2) matrix defining lower and upper bounds for the parameters.
-% o varargin list of argument following mh_bounds
+% o dataset_ data structure
+% o options_ options structure
+% o M_ model structure
+% o estim_params_ estimated parameters structure
+% o bayestopt_ estimation options structure
+% o oo_ outputs structure
%
% OUTPUTS
% o record [struct] structure describing the iterations
@@ -51,12 +56,17 @@ function record=random_walk_metropolis_hastings(TargetFun,ProposalFun,xparam1,vv
% You should have received a copy of the GNU General Public License
% along with Dynare. If not, see <http://www.gnu.org/licenses/>.
-global M_ options_ bayestopt_ estim_params_ oo_
+
+% In Metropolis, we set penalty to Inf to as to reject all parameter sets
+% triggering error in target density computation
+
+bayestopt_.penalty = Inf;
+
%%%%
%%%% Initialization of the random walk metropolis-hastings chains.
%%%%
[ ix2, ilogpo2, ModelName, MhDirectoryName, fblck, fline, npar, nblck, nruns, NewFile, MAX_nruns, d ] = ...
- metropolis_hastings_initialization(TargetFun, xparam1, vv, mh_bounds, varargin{:});
+ metropolis_hastings_initialization(TargetFun, xparam1, vv, mh_bounds,dataset_,options_,M_,estim_params_,bayestopt_,oo_);
InitSizeArray = min([repmat(MAX_nruns,nblck,1) fline+nruns-1],[],2);
@@ -93,10 +103,16 @@ localVars = struct('TargetFun', TargetFun, ...
'nruns', nruns, ...
'NewFile', NewFile, ...
'MAX_nruns', MAX_nruns, ...
- 'd', d);
-localVars.InitSizeArray=InitSizeArray;
-localVars.record=record;
-localVars.varargin=varargin;
+ 'd', d, ...
+ 'InitSizeArray',InitSizeArray, ...
+ 'record', record, ...
+ 'dataset_', dataset_, ...
+ 'options_', options_, ...
+ 'M_',M_, ...
+ 'bayestopt_', bayestopt_, ...
+ 'estim_params_', estim_params_, ...
+ 'oo_', oo_,...
+ 'varargin',[]);
% The user don't want to use parallel computing, or want to compute a
@@ -110,11 +126,7 @@ if isnumeric(options_.parallel) || (nblck-fblck)==0,
% Parallel in Local or remote machine.
else
% Global variables for parallel routines.
- globalVars = struct('M_',M_, ...
- 'options_', options_, ...
- 'bayestopt_', bayestopt_, ...
- 'estim_params_', estim_params_, ...
- 'oo_', oo_);
+ globalVars = struct();
% which files have to be copied to run remotely
NamFileInput(1,:) = {'',[ModelName '_static.m']};
diff --git a/matlab/random_walk_metropolis_hastings_core.m b/matlab/random_walk_metropolis_hastings_core.m
index 1898d478ab48dca812c639a50c11bc3ba09a7de1..84f180c97d2472782fd05a37f4755c49155a0ef0 100644
--- a/matlab/random_walk_metropolis_hastings_core.m
+++ b/matlab/random_walk_metropolis_hastings_core.m
@@ -68,9 +68,6 @@ if nargin<4,
whoiam=0;
end
-
-global bayestopt_ estim_params_ options_ M_ oo_
-
% reshape 'myinputs' for local computation.
% In order to avoid confusion in the name space, the instruction struct2local(myinputs) is replaced by:
@@ -90,6 +87,12 @@ MAX_nruns=myinputs.MAX_nruns;
d=myinputs.d;
InitSizeArray=myinputs.InitSizeArray;
record=myinputs.record;
+dataset_ = myinputs.dataset_;
+bayestopt_ = myinputs.bayestopt_;
+estim_params_ = myinputs.estim_params_;
+options_ = myinputs.options_;
+M_ = myinputs.M_;
+oo_ = myinputs.oo_;
varargin=myinputs.varargin;
% Necessary only for remote computing!
@@ -100,9 +103,6 @@ if whoiam
bayestopt_.p3,bayestopt_.p4,1);
end
-% (re)Set the penalty
-bayestopt_.penalty = Inf;
-
MhDirectoryName = CheckPath('metropolis',M_.dname);
options_.lik_algo = 1;
@@ -168,7 +168,7 @@ for b = fblck:nblck,
par = feval(ProposalFun, ix2(b,:), proposal_covariance_Cholesky_decomposition, n);
if all( par(:) > mh_bounds(:,1) ) && all( par(:) < mh_bounds(:,2) )
try
- logpost = - feval(TargetFun, par(:),varargin{:});
+ logpost = - feval(TargetFun, par(:),dataset_,options_,M_,estim_params_,bayestopt_,oo_);
catch
logpost = -inf;
end
diff --git a/matlab/sim1_purely_backward.m b/matlab/sim1_purely_backward.m
index b432c4789f2cfb9d970cbf31805c367299747752..beb699880d253338e69620dc824ddef18b6ddf81 100644
--- a/matlab/sim1_purely_backward.m
+++ b/matlab/sim1_purely_backward.m
@@ -34,7 +34,11 @@ function sim1_purely_backward()
yb = oo_.endo_simul(:,it-1); % Values at previous period, also used as guess value for current period
yb1 = yb(iyb);
- tmp = solve1(model_dynamic, [yb1; yb], 1:M_.endo_nbr, nyb+1:nyb+M_.endo_nbr, 1, 1, oo_.exo_simul, M_.params, oo_.steady_state, it);
+ tmp = solve1(model_dynamic, [yb1; yb], 1:M_.endo_nbr, nyb+1:nyb+ ...
+ M_.endo_nbr, 1, 1, options_.gstep, ...
+ options_.solve_tolf,options_.solve_tolx, ...
+ options_.solve_maxit,options_.debug,oo_.exo_simul, ...
+ M_.params, oo_.steady_state, it);
oo_.endo_simul(:,it) = tmp(nyb+1:nyb+M_.endo_nbr);
end
\ No newline at end of file
diff --git a/matlab/sim1_purely_forward.m b/matlab/sim1_purely_forward.m
index 65a874b5055f035e8c0657712cb135bb994dd64c..795cb25e281c305f242ea56cf21ffc3467bc795c 100644
--- a/matlab/sim1_purely_forward.m
+++ b/matlab/sim1_purely_forward.m
@@ -33,7 +33,11 @@ function sim1_purely_forward()
yf = oo_.endo_simul(:,it+1); % Values at next period, also used as guess value for current period
yf1 = yf(iyf);
- tmp = solve1(model_dynamic, [yf; yf1], 1:M_.endo_nbr, 1:M_.endo_nbr, 1, 1, oo_.exo_simul, M_.params, oo_.steady_state, it);
+ tmp = solve1(model_dynamic, [yf; yf1], 1:M_.endo_nbr, 1:M_.endo_nbr, ...
+ 1, 1, options_.gstep, options_.solve_tolf, ...
+ options_.solve_tolx, options_.solve_maxit, ...
+ options_.debug,oo_.exo_simul, M_.params, oo_.steady_state, ...
+ it);
oo_.endo_simul(:,it) = tmp(1:M_.endo_nbr);
end
\ No newline at end of file
diff --git a/matlab/simul_backward_nonlinear_model.m b/matlab/simul_backward_nonlinear_model.m
index 546990cb5eb656cd581f22117b77937563febcfc..80dec87a2720502e510271791c5947db58c411a7 100644
--- a/matlab/simul_backward_nonlinear_model.m
+++ b/matlab/simul_backward_nonlinear_model.m
@@ -104,6 +104,10 @@ DynareOutput.endo_simul(:,1) = DynareOutput.steady_state;
for it = 2:sample_size+1
y(jdx) = DynareOutput.endo_simul(:,it-1); % A good guess for the initial conditions is the previous values for the endogenous variables.
y(hdx) = y(jdx(iy1)); % Set lagged variables.
- y(jdx) = solve1(model_dynamic, y, idx, jdx, 1, 1, DynareOutput.exo_simul, DynareModel.params, DynareOutput.steady_state, it);
+ y(jdx) = solve1(model_dynamic, y, idx, jdx, 1, 1, DynareOptions.gstep, ...
+ DynareOptions.solve_tolf,DynareOptions.solve_tolx, ...
+ DynareOptions.solve_maxit,DynareOptions.debug, ...
+ DynareOutput.exo_simul, DynareModel.params, ...
+ DynareOutput.steady_state, it);
DynareOutput.endo_simul(:,it) = y(jdx);
end
\ No newline at end of file
diff --git a/matlab/solve1.m b/matlab/solve1.m
index 74111276f6468499888ea48c1f86f46792d75f02..97cb5e24e3fab5514f7c29f3010d5aa0c5f06bd7 100644
--- a/matlab/solve1.m
+++ b/matlab/solve1.m
@@ -1,4 +1,4 @@
-function [x,check] = solve1(func,x,j1,j2,jacobian_flag,bad_cond_flag,varargin)
+function [x,check] = solve1(func,x,j1,j2,jacobian_flag,bad_cond_flag,gstep,tolf,tolx,maxit,debug,varargin)
% function [x,check] = solve1(func,x,j1,j2,jacobian_flag,bad_cond_flag,varargin)
% Solves systems of non linear equations of several variables
%
@@ -11,6 +11,12 @@ function [x,check] = solve1(func,x,j1,j2,jacobian_flag,bad_cond_flag,varargin)
% jacobian_flag=0: jacobian obtained numerically
% bad_cond_flag=1: when Jacobian is badly conditionned, use an
% alternative formula to Newton step
+% gstep increment multiplier in numercial derivative
+% computation
+% tolf tolerance for residuals
+% tolx tolerance for solution variation
+% maxit maximum number of iterations
+% debug debug flag
% varargin: list of arguments following bad_cond_flag
%
% OUTPUTS
@@ -37,18 +43,13 @@ function [x,check] = solve1(func,x,j1,j2,jacobian_flag,bad_cond_flag,varargin)
% You should have received a copy of the GNU General Public License
% along with Dynare. If not, see <http://www.gnu.org/licenses/>.
-global M_ options_ fjac
-
nn = length(j1);
fjac = zeros(nn,nn) ;
g = zeros(nn,1) ;
-tolf = options_.solve_tolf ;
-tolx = options_.solve_tolx;
tolmin = tolx ;
stpmx = 100 ;
-maxit = options_.solve_maxit ;
check = 0 ;
@@ -77,7 +78,7 @@ for its = 1:maxit
fvec = fvec(j1);
fjac = fjac(j1,j2);
else
- dh = max(abs(x(j2)),options_.gstep(1)*ones(nn,1))*eps^(1/3);
+ dh = max(abs(x(j2)),gstep(1)*ones(nn,1))*eps^(1/3);
for j = 1:nn
xdh = x ;
@@ -89,33 +90,10 @@ for its = 1:maxit
end
g = (fvec'*fjac)';
- if options_.debug
+ if debug
disp(['cond(fjac) ' num2str(cond(fjac))])
end
- M_.unit_root = 0;
- if M_.unit_root
- if first_time
- first_time = 0;
- [q,r,e]=qr(fjac);
- n = sum(abs(diag(r)) < 1e-12);
- fvec = q'*fvec;
- p = e*[-r(1:end-n,1:end-n)\fvec(1:end-n);zeros(n,1)];
- disp(' ')
- disp('STEADY with unit roots:')
- disp(' ')
- if n > 0
- disp([' The following variable(s) kept their value given in INITVAL' ...
- ' or ENDVAL'])
- disp(char(e(:,end-n+1:end)'*M_.endo_names))
- else
- disp(' STEADY can''t find any unit root!')
- end
- else
- [q,r]=qr(fjac*e);
- fvec = q'*fvec;
- p = e*[-r(1:end-n,1:end-n)\fvec(1:end-n);zeros(n,1)];
- end
- elseif bad_cond_flag && cond(fjac) > 1/sqrt(eps)
+ if bad_cond_flag && rcond(fjac) < sqrt(eps)
fjac2=fjac'*fjac;
p=-(fjac2+1e6*sqrt(nn*eps)*max(sum(abs(fjac2)))*eye(nn))\(fjac'*fvec);
else
@@ -126,7 +104,7 @@ for its = 1:maxit
[x,f,fvec,check]=lnsrch1(xold,fold,g,p,stpmax,func,j1,j2,varargin{:});
- if options_.debug
+ if debug
disp([its f])
disp([xold x])
end
diff --git a/matlab/stochastic_solvers.m b/matlab/stochastic_solvers.m
index 84f2c312fe7cb9f1082ef18db45e79baac3ee49b..8a688a6269b6f65e517761502040d6392c0be740 100644
--- a/matlab/stochastic_solvers.m
+++ b/matlab/stochastic_solvers.m
@@ -120,6 +120,13 @@ if ~isreal(jacobia_)
end
end
+if any(any(isnan(jacobia_)))
+ info(1) = 8;
+ NaN_params=find(isnan(M_.params));
+ info(2:length(NaN_params)+1) = NaN_params;
+ return
+end
+
kstate = dr.kstate;
kad = dr.kad;
kae = dr.kae;
diff --git a/mex/build/matlab/configure.ac b/mex/build/matlab/configure.ac
index 4a3f115d79d463506c0a6457f20042f1971fd699..ca948b4dba099f07da4484df7401bc7c7209d785 100644
--- a/mex/build/matlab/configure.ac
+++ b/mex/build/matlab/configure.ac
@@ -1,6 +1,6 @@
dnl Process this file with autoconf to produce a configure script.
-dnl Copyright (C) 2009-2011 Dynare Team
+dnl Copyright (C) 2009-2012 Dynare Team
dnl
dnl This file is part of Dynare.
dnl
@@ -17,7 +17,7 @@ dnl
dnl You should have received a copy of the GNU General Public License
dnl along with Dynare. If not, see <http://www.gnu.org/licenses/>.
-AC_PREREQ([2.61])
+AC_PREREQ([2.62])
AC_INIT([dynare], [4.4-unstable])
AC_CONFIG_SRCDIR([configure.ac])
AM_INIT_AUTOMAKE([-Wall -Werror foreign])
diff --git a/mex/build/octave/Makefile.am b/mex/build/octave/Makefile.am
index 132ddbe15889f4c84ed13e1e636454fdd66da581..2689afab4a680d4fec1573075bd34ff89cd2ca75 100644
--- a/mex/build/octave/Makefile.am
+++ b/mex/build/octave/Makefile.am
@@ -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 sobol local_state_space_iterations
+SUBDIRS = mjdgges kronecker bytecode libdynare++ gensylv k_order_perturbation dynare_simul_ qzcomplex ordschur block_kalman_filter sobol local_state_space_iterations linsolve
if HAVE_GSL
if HAVE_MATIO
diff --git a/mex/build/octave/configure.ac b/mex/build/octave/configure.ac
index 34855798155da65422b3d43aaead3a97fbdbeaed..0a136279976f1667c1d903037faaacc6635b9ef0 100644
--- a/mex/build/octave/configure.ac
+++ b/mex/build/octave/configure.ac
@@ -1,6 +1,6 @@
dnl Process this file with autoconf to produce a configure script.
-dnl Copyright (C) 2009-2011 Dynare Team
+dnl Copyright (C) 2009-2012 Dynare Team
dnl
dnl This file is part of Dynare.
dnl
@@ -17,7 +17,7 @@ dnl
dnl You should have received a copy of the GNU General Public License
dnl along with Dynare. If not, see <http://www.gnu.org/licenses/>.
-AC_PREREQ([2.61])
+AC_PREREQ([2.62])
AC_INIT([dynare], [4.4-unstable])
AC_CONFIG_SRCDIR([configure.ac])
AM_INIT_AUTOMAKE([-Wall -Werror foreign])
@@ -132,6 +132,7 @@ AC_CONFIG_FILES([Makefile
ms_sbvar/Makefile
block_kalman_filter/Makefile
sobol/Makefile
- local_state_space_iterations/Makefile])
+ local_state_space_iterations/Makefile
+ linsolve/Makefile])
AC_OUTPUT
diff --git a/mex/build/octave/linsolve/Makefile.am b/mex/build/octave/linsolve/Makefile.am
new file mode 100644
index 0000000000000000000000000000000000000000..4f3a5aaf195c136db20d93a63286d9e666e172a5
--- /dev/null
+++ b/mex/build/octave/linsolve/Makefile.am
@@ -0,0 +1,6 @@
+EXEEXT = .oct
+include ../mex.am
+
+noinst_PROGRAMS = linsolve
+
+nodist_linsolve_SOURCES = $(top_srcdir)/../../sources/linsolve/linsolve.cc
diff --git a/mex/sources/Makefile.am b/mex/sources/Makefile.am
index 41b94ffa67d9bd6a879dbca859f322daabaacc53..57934b16407763c5867c47ed8cf96eddd53ab7ab 100644
--- a/mex/sources/Makefile.am
+++ b/mex/sources/Makefile.am
@@ -15,7 +15,8 @@ EXTRA_DIST = \
ms-sbvar \
block_kalman_filter \
sobol \
- local_state_space_iterations
+ local_state_space_iterations \
+ linsolve
clean-local:
rm -rf `find mex/sources -name *.o`
diff --git a/mex/sources/estimation/logMHMCMCposterior.cc b/mex/sources/estimation/logMHMCMCposterior.cc
index c388767af8d8114f7dca05089337eb5c8d71cfcc..a9308c0d47584e7012773be5a2d9c4016fb970f8 100644
--- a/mex/sources/estimation/logMHMCMCposterior.cc
+++ b/mex/sources/estimation/logMHMCMCposterior.cc
@@ -130,7 +130,13 @@ sampleMHMC(LogPosteriorDensity &lpd, RandomWalkMetropolisHastings &rwmh,
MATFile *drawmat; // MCMC draws output file pointer
int matfStatus;
#else // OCTAVE_MEX_FILE e.t.c.
+# if MATIO_MAJOR_VERSION > 1 || (MATIO_MAJOR_VERSION == 1 && MATIO_MINOR_VERSION >= 5)
+ size_t dims[2];
+ const matio_compression compression = MAT_COMPRESSION_NONE;
+# else
int dims[2];
+ const int compression = COMPRESSION_NONE;
+# endif
mat_t *drawmat;
matvar_t *matvar;
int matfStatus;
@@ -304,7 +310,7 @@ sampleMHMC(LogPosteriorDensity &lpd, RandomWalkMetropolisHastings &rwmh,
}
else
{
- int start[2] = {0, 0}, edge[2] = {2, 2}, stride[2] = {1, 1}, err = 0;
+ int start[2] = {0, 0}, edge[2], stride[2] = {1, 1}, err = 0;
mexPrintf("MHMCMC: Using interim partial draws file %s \n", mhFName.c_str());
// matvar = Mat_VarReadInfo(drawmat, "x2");
matvar = Mat_VarReadInfo(drawmat, (char *) "x2");
@@ -317,9 +323,9 @@ sampleMHMC(LogPosteriorDensity &lpd, RandomWalkMetropolisHastings &rwmh,
else
{
// GetVariable(drawmat, "x2");
- dims[0] = matvar->dims[0]-1;
- dims[1] = matvar->dims[1]-1;
- err = Mat_VarReadData(drawmat, matvar, mxGetPr(mxMhParamDrawsPtr), start, stride, matvar->dims);
+ edge[0] = matvar->dims[0];
+ edge[1] = matvar->dims[1];
+ err = Mat_VarReadData(drawmat, matvar, mxGetPr(mxMhParamDrawsPtr), start, stride, edge);
if (err)
{
fline(b) = 1;
@@ -339,9 +345,9 @@ sampleMHMC(LogPosteriorDensity &lpd, RandomWalkMetropolisHastings &rwmh,
else
{
// GetVariable(drawmat, "x2");
- dims[0] = matvar->dims[0]-1;
- dims[1] = matvar->dims[1]-1;
- err = Mat_VarReadData(drawmat, matvar, mxGetPr(mxMhLogPostDensPtr), start, stride, matvar->dims);
+ edge[0] = matvar->dims[0];
+ edge[1] = matvar->dims[1];
+ err = Mat_VarReadData(drawmat, matvar, mxGetPr(mxMhLogPostDensPtr), start, stride, edge);
if (err)
{
fline(b) = 1;
@@ -513,7 +519,7 @@ sampleMHMC(LogPosteriorDensity &lpd, RandomWalkMetropolisHastings &rwmh,
dims[0] = currInitSizeArray;
dims[1] = npar;
matvar = Mat_VarCreate("x2", MAT_C_DOUBLE, MAT_T_DOUBLE, 2, dims, mxGetPr(mxMhParamDrawsPtr), 0);
- matfStatus = Mat_VarWrite(drawmat, matvar, 0);
+ matfStatus = Mat_VarWrite(drawmat, matvar, compression);
Mat_VarFree(matvar);
if (matfStatus)
{
@@ -523,7 +529,7 @@ sampleMHMC(LogPosteriorDensity &lpd, RandomWalkMetropolisHastings &rwmh,
//matfStatus = matPutVariable(drawmat, "logpo2", mxMhLogPostDensPtr);
dims[1] = 1;
matvar = Mat_VarCreate("logpo2", MAT_C_DOUBLE, MAT_T_DOUBLE, 2, dims, mxGetPr(mxMhLogPostDensPtr), 0);
- matfStatus = Mat_VarWrite(drawmat, matvar, 0);
+ matfStatus = Mat_VarWrite(drawmat, matvar, compression);
Mat_VarFree(matvar);
if (matfStatus)
{
diff --git a/mex/sources/linsolve/linsolve.cc b/mex/sources/linsolve/linsolve.cc
new file mode 100644
index 0000000000000000000000000000000000000000..6c43f3f195ea6cf25bcb55ecb3e1dd44abc032e2
--- /dev/null
+++ b/mex/sources/linsolve/linsolve.cc
@@ -0,0 +1,120 @@
+/*
+ * Oct-file for bringing MATLAB's linsolve function to Octave.
+ *
+ * The implementation is incomplete:
+ * - it only knows about the TRANSA, LT, UT and SYM options
+ * - it only works with square matrices
+ * - it only works on double matrices (no single precision or complex)
+ *
+ * Written by Sebastien Villemot <sebastien.villemot@ens.fr>.
+ */
+
+/*
+ * Copyright (C) 2012 Dynare Team
+ *
+ * This file is part of Dynare.
+ *
+ * Dynare is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * Dynare is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with Dynare. If not, see <http://www.gnu.org/licenses/>.
+ */
+
+
+#include <octave/oct.h>
+#include <octave/ov-struct.h>
+
+DEFUN_DLD(linsolve, args, nargout, "-*- texinfo -*-\n\
+@deftypefn {Loadable Function} @var{x} = linsolve (@var{a}, @var{b})\n\
+@deftypefnx {Loadable Function} [ @var{x}, @var{r} ] = linsolve (@var{a}, @var{b})\n\
+@deftypefnx {Loadable Function} @var{x} = linsolve (@var{a}, @var{b}, @var{options})\n\
+@deftypefnx {Loadable Function} [ @var{x}, @var{r} ] = linsolve (@var{a}, @var{b}, @var{options})\n\
+\n\
+Solves the linear system @math{A*X = B} and returns @var{X}.\n\
+\n\
+Alternatively, if @var{options} is provided and has a field @code{TRANSA} equal \
+to @code{true}, then it solves the system @math{A'*X = B}.\n\
+\n\
+Also, the @code{LT} field of @var{options} (resp. the @code{UT} field) can be set \
+to @code{true} to indicate that the matrix @var{a} is lower (resp. upper) \
+triangular; similarly, the @code{SYM} field can be set to @code{true} to \
+indicate that the matrix is symmetric.\n\
+\n\
+If requested, @var{r} will contain the reciprocal condition number.\n\
+@end deftypefn\n\
+")
+{
+ int nargin = args.length();
+ octave_value_list retval;
+
+ if (nargin > 3 || nargin < 2 || nargout > 2)
+ {
+ print_usage();
+ return retval;
+ }
+
+ Matrix A = args(0).matrix_value();
+ Matrix B = args(1).matrix_value();
+ if (error_state)
+ return retval;
+
+ dim_vector dimA = A.dims();
+ dim_vector dimB = B.dims();
+
+ if (dimA(0) != dimB(0))
+ {
+ error("linsolve: must have same number of lines in A and B");
+ return retval;
+ }
+
+ if (dimA(0) != dimA(1))
+ {
+ error("linsolve: rectangular A not yet supported");
+ return retval;
+ }
+
+
+ MatrixType typA;
+ typA.mark_as_full();
+
+ bool transa = false;
+ if (nargin == 3)
+ {
+ octave_scalar_map opts = args(2).scalar_map_value();
+ if (error_state)
+ return retval;
+
+ octave_value tmp = opts.contents("TRANSA");
+ transa = tmp.is_defined() && tmp.bool_matrix_value().elem(0);
+
+ tmp = opts.contents("UT");
+ if (tmp.is_defined() && tmp.bool_matrix_value().elem(0))
+ typA.mark_as_upper_triangular();
+
+ tmp = opts.contents("LT");
+ if (tmp.is_defined() && tmp.bool_matrix_value().elem(0))
+ typA.mark_as_lower_triangular();
+
+ tmp = opts.contents("SYM");
+ if (tmp.is_defined() && tmp.bool_matrix_value().elem(0))
+ typA.mark_as_symmetric();
+ }
+
+ double rcond;
+ octave_idx_type info;
+
+ retval(0) = A.solve(typA, B, info, rcond, NULL, true, transa ? blas_trans : blas_no_trans);
+
+ if (nargout == 2)
+ retval(1) = rcond;
+
+ return retval;
+}
diff --git a/osx/README.txt b/osx/README.txt
index ef1331f87d5681e3c1fd66119934ed9df3d60972..0760d5051ba185758734dc63942014f3ac75aed1 100644
--- a/osx/README.txt
+++ b/osx/README.txt
@@ -53,13 +53,14 @@ specify a MOD file.
Using Dynare with Octave
------------------------
-Dynare is now available for GNU Octave, a free clone of MATLAB (R) (see
+Dynare also works on top of GNU Octave, a free clone of MATLAB (R) (see
<http://www.octave.org>).
-The recommended Octave distribution is the Octave 3.4
-precompiled binaries from Octave Forge, available at:
+This version of Dynare is compiled for Octave 3.6.2, and may not work
+with other versions of Octave. You can download Octave through the Homebrew
+package manager:
- http://sourceforge.net/projects/octave/files/Octave%20MacOSX%20Binary/
+ http://mxcl.github.com/homebrew
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'
@@ -75,10 +76,6 @@ You can test your installation by typing 'dynare' at the Octave prompt. This
should give you an error message complaining that you did not specify a MOD
file.
-For more information about Dynare for Octave, go to:
-
- http://www.dynare.org/DynareWiki/DynareOctave
-
Dynamic Loadable Libraries
--------------------------
@@ -91,13 +88,16 @@ If the libraries are correctly detected by MATLAB (R) or Octave, the following s
be displayed when you launch Dynare:
Configuring Dynare ...
+
[mex] Generalized QZ.
[mex] Sylvester equation solution.
[mex] Kronecker products.
[mex] Sparse kronecker products.
+ [mex] Local state space iteration (second order).
[mex] Bytecode evaluation.
[mex] k-order perturbation solver.
[mex] k-order solution simulation.
+ [mex] Quasi Monte-Carlo sequence (Sobol).
On the contrary, if the libraries are not detected, Dynare will fallback on
slower alternatives written in M-files (which exist for some of the libraries),
@@ -108,9 +108,11 @@ and display the following:
[m] Sylvester equation solution.
[m] Kronecker products.
[m] Sparse kronecker products.
+ [m] Local state space iteration (second order).
[no] Bytecode evaluation.
[no] k-order perturbation solver.
[no] k-order solution simulation.
+ [no] Quasi Monte-Carlo sequence (Sobol).
In this last case, Dynare will run correctly on the basic features,
but with suboptimal speed, and some features will be missing. There
diff --git a/osx/createOsxFolderForPkg.sh b/osx/createOsxFolderForPkg.sh
index d429886b636b00ff6dfab8f4b9897f276fdb392f..5aa234e40baafa653e4d062d32febde1064794b0 100755
--- a/osx/createOsxFolderForPkg.sh
+++ b/osx/createOsxFolderForPkg.sh
@@ -47,10 +47,19 @@ cp -r $TOP_DYN_DIR/examples $INSTALLDIR
##########################################################
# FIRST BUILD 32 BIT EVERYTHING, 32 BIT MATLAB < 7.5 MEX #
##########################################################
-./configure CFLAGS='-arch i386' CXXFLAGS='-arch i386' FFLAGS='-arch i386' LDFLAGS='-arch i386' --with-matlab=/Applications/MATLAB/R2007a MATLAB_VERSION=7.4 --with-gsl=/usr/local/Cellar/gsl/1.15_32bit
+./configure CFLAGS='-arch i386' CXXFLAGS='-arch i386' FFLAGS='-arch i386' LDFLAGS='-arch i386' --with-matlab=/Applications/MATLAB/R2007a MATLAB_VERSION=7.4 --with-gsl=/usr/local32
make pdf
+
+cd $TOP_DYN_DIR/preprocessor
+make
+
+cd $TOP_DYN_DIR/dynare++
+make
+
+cd $TOP_DYN_DIR/mex/build/matlab
make
+
# Matlab
# Must come after configure because matlab/dynare_version.m is created by configure script
cp -r $TOP_DYN_DIR/matlab $INSTALLDIR
@@ -74,21 +83,6 @@ cp $TOP_DYN_DIR/mex/build/matlab/mjdgges/*.mexmaci
cp $TOP_DYN_DIR/mex/build/matlab/ms_sbvar/*.mexmaci $INSTALLDIR/mex/matlab/osx32-7.4
cp $TOP_DYN_DIR/mex/build/matlab/sobol/*.mexmaci $INSTALLDIR/mex/matlab/osx32-7.4
-# Octave
-cp $TOP_DYN_DIR/mex/build/octave/block_kalman_filter/*.mex $INSTALLDIR/mex/octave
-cp $TOP_DYN_DIR/mex/build/octave/bytecode/*.mex $INSTALLDIR/mex/octave
-cp $TOP_DYN_DIR/mex/build/octave/dynare_simul_/*.mex $INSTALLDIR/mex/octave
-cp $TOP_DYN_DIR/mex/build/octave/gensylv/*.mex $INSTALLDIR/mex/octave
-cp $TOP_DYN_DIR/mex/build/octave/k_order_perturbation/*.mex $INSTALLDIR/mex/octave
-cp $TOP_DYN_DIR/mex/build/octave/kalman_steady_state/*.mex $INSTALLDIR/mex/octave
-cp $TOP_DYN_DIR/mex/build/octave/kronecker/*.mex $INSTALLDIR/mex/octave
-cp $TOP_DYN_DIR/mex/build/octave/local_state_space_iterations/*.mex $INSTALLDIR/mex/octave
-cp $TOP_DYN_DIR/mex/build/octave/mjdgges/*.mex $INSTALLDIR/mex/octave
-cp $TOP_DYN_DIR/mex/build/octave/ms_sbvar/*.mex $INSTALLDIR/mex/octave
-cp $TOP_DYN_DIR/mex/build/octave/ordschur/*.oct $INSTALLDIR/mex/octave
-cp $TOP_DYN_DIR/mex/build/octave/qzcomplex/*.oct $INSTALLDIR/mex/octave
-cp $TOP_DYN_DIR/mex/build/octave/sobol/*.mex $INSTALLDIR/mex/octave
-
# dynare++
cp $TOP_DYN_DIR/dynare++/src/dynare++ $INSTALLDIR/dynare++
cp $TOP_DYN_DIR/dynare++/extern/matlab/dynare_simul.m $INSTALLDIR/dynare++
@@ -118,7 +112,7 @@ cp $TOP_DYN_DIR/dynare++/kord/kord.pdf
##############################################
make clean
cd $TOP_DYN_DIR/mex/build/matlab
-./configure CFLAGS='-arch i386' CXXFLAGS='-arch i386' FFLAGS='-arch i386' LDFLAGS='-arch i386' --with-matlab=/Applications/MATLAB/MATLAB_R2009b_32bit/MATLAB_R2009b.app MATLAB_VERSION=7.9 MEXEXT='mexmaci' --with-gsl=/usr/local/Cellar/gsl/1.15_32bit
+./configure CFLAGS='-arch i386' CXXFLAGS='-arch i386' FFLAGS='-arch i386' LDFLAGS='-arch i386' --with-matlab=/Applications/MATLAB/MATLAB_R2009b_32bit/MATLAB_R2009b.app MATLAB_VERSION=7.9 MEXEXT='mexmaci' --with-gsl=/usr/local32
make
# Matlab
@@ -140,7 +134,7 @@ cp $TOP_DYN_DIR/mex/build/matlab/sobol/*.mexmaci
#####################################
make clean
cd $TOP_DYN_DIR/mex/build/matlab
-./configure --with-matlab=/Applications/MATLAB/MATLAB_R2009b.app MATLAB_VERSION=7.9 --with-gsl=/usr/local/Cellar/gsl/1.15
+./configure --with-matlab=/Applications/MATLAB/MATLAB_R2009b.app MATLAB_VERSION=7.9
make
# Matlab
@@ -156,6 +150,30 @@ cp $TOP_DYN_DIR/mex/build/matlab/mjdgges/*.mexmaci64
cp $TOP_DYN_DIR/mex/build/matlab/ms_sbvar/*.mexmaci64 $INSTALLDIR/mex/matlab/osx64
cp $TOP_DYN_DIR/mex/build/matlab/sobol/*.mexmaci64 $INSTALLDIR/mex/matlab/osx64
+#####################################
+# RETURN TO BUILD 64 BIT OCTAVE MEX #
+#####################################
+make clean
+cd $TOP_DYN_DIR/mex/build/octave
+./configure
+make
+
+# Octave
+cp $TOP_DYN_DIR/mex/build/octave/block_kalman_filter/*.mex $INSTALLDIR/mex/octave
+cp $TOP_DYN_DIR/mex/build/octave/bytecode/*.mex $INSTALLDIR/mex/octave
+cp $TOP_DYN_DIR/mex/build/octave/dynare_simul_/*.mex $INSTALLDIR/mex/octave
+cp $TOP_DYN_DIR/mex/build/octave/gensylv/*.mex $INSTALLDIR/mex/octave
+cp $TOP_DYN_DIR/mex/build/octave/k_order_perturbation/*.mex $INSTALLDIR/mex/octave
+cp $TOP_DYN_DIR/mex/build/octave/kalman_steady_state/*.mex $INSTALLDIR/mex/octave
+cp $TOP_DYN_DIR/mex/build/octave/kronecker/*.mex $INSTALLDIR/mex/octave
+cp $TOP_DYN_DIR/mex/build/octave/local_state_space_iterations/*.mex $INSTALLDIR/mex/octave
+cp $TOP_DYN_DIR/mex/build/octave/mjdgges/*.mex $INSTALLDIR/mex/octave
+cp $TOP_DYN_DIR/mex/build/octave/ms_sbvar/*.mex $INSTALLDIR/mex/octave
+cp $TOP_DYN_DIR/mex/build/octave/ordschur/*.oct $INSTALLDIR/mex/octave
+cp $TOP_DYN_DIR/mex/build/octave/qzcomplex/*.oct $INSTALLDIR/mex/octave
+cp $TOP_DYN_DIR/mex/build/octave/sobol/*.mex $INSTALLDIR/mex/octave
+
+
# clean everything
cd $TOP_DYN_DIR
make distclean
@@ -169,4 +187,4 @@ chmod -R g+w $INSTALLDIR
echo DONE
# NEED TO BUILD DYNARE.HTML DOCUMENTION ON DEBIAN
-# AND INCLUDE IN DISTRIBUTION BY HAND
\ No newline at end of file
+# AND INCLUDE IN DISTRIBUTION BY HAND
diff --git a/preprocessor/FlexLexer.h b/preprocessor/FlexLexer.h
deleted file mode 100644
index bad4ce03fb65bdcec983d91b7de3fa1265515170..0000000000000000000000000000000000000000
--- a/preprocessor/FlexLexer.h
+++ /dev/null
@@ -1,206 +0,0 @@
-// -*-C++-*-
-// FlexLexer.h -- define interfaces for lexical analyzer classes generated
-// by flex
-
-// Copyright (c) 1993 The Regents of the University of California.
-// All rights reserved.
-//
-// This code is derived from software contributed to Berkeley by
-// Kent Williams and Tom Epperly.
-//
-// Redistribution and use in source and binary forms, with or without
-// modification, are permitted provided that the following conditions
-// are met:
-
-// 1. Redistributions of source code must retain the above copyright
-// notice, this list of conditions and the following disclaimer.
-// 2. Redistributions in binary form must reproduce the above copyright
-// notice, this list of conditions and the following disclaimer in the
-// documentation and/or other materials provided with the distribution.
-
-// Neither the name of the University nor the names of its contributors
-// may be used to endorse or promote products derived from this software
-// without specific prior written permission.
-
-// THIS SOFTWARE IS PROVIDED ``AS IS'' AND WITHOUT ANY EXPRESS OR
-// IMPLIED WARRANTIES, INCLUDING, WITHOUT LIMITATION, THE IMPLIED
-// WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
-// PURPOSE.
-
-// This file defines FlexLexer, an abstract class which specifies the
-// external interface provided to flex C++ lexer objects, and yyFlexLexer,
-// which defines a particular lexer class.
-//
-// If you want to create multiple lexer classes, you use the -P flag
-// to rename each yyFlexLexer to some other xxFlexLexer. You then
-// include <FlexLexer.h> in your other sources once per lexer class:
-//
-// #undef yyFlexLexer
-// #define yyFlexLexer xxFlexLexer
-// #include <FlexLexer.h>
-//
-// #undef yyFlexLexer
-// #define yyFlexLexer zzFlexLexer
-// #include <FlexLexer.h>
-// ...
-
-#ifndef __FLEX_LEXER_H
-// Never included before - need to define base class.
-#define __FLEX_LEXER_H
-
-#include <iostream>
-# ifndef FLEX_STD
-# define FLEX_STD std::
-# endif
-
-extern "C++" {
-
-struct yy_buffer_state;
-typedef int yy_state_type;
-
-class FlexLexer {
-public:
- virtual ~FlexLexer() { }
-
- const char* YYText() const { return yytext; }
- int YYLeng() const { return yyleng; }
-
- virtual void
- yy_switch_to_buffer( struct yy_buffer_state* new_buffer ) = 0;
- virtual struct yy_buffer_state*
- yy_create_buffer( FLEX_STD istream* s, int size ) = 0;
- virtual void yy_delete_buffer( struct yy_buffer_state* b ) = 0;
- virtual void yyrestart( FLEX_STD istream* s ) = 0;
-
- virtual int yylex() = 0;
-
- // Call yylex with new input/output sources.
- int yylex( FLEX_STD istream* new_in, FLEX_STD ostream* new_out = 0 )
- {
- switch_streams( new_in, new_out );
- return yylex();
- }
-
- // Switch to new input/output streams. A nil stream pointer
- // indicates "keep the current one".
- virtual void switch_streams( FLEX_STD istream* new_in = 0,
- FLEX_STD ostream* new_out = 0 ) = 0;
-
- int lineno() const { return yylineno; }
-
- int debug() const { return yy_flex_debug; }
- void set_debug( int flag ) { yy_flex_debug = flag; }
-
-protected:
- char* yytext;
- int yyleng;
- int yylineno; // only maintained if you use %option yylineno
- int yy_flex_debug; // only has effect with -d or "%option debug"
-};
-
-}
-#endif // FLEXLEXER_H
-
-#if defined(yyFlexLexer) || ! defined(yyFlexLexerOnce)
-// Either this is the first time through (yyFlexLexerOnce not defined),
-// or this is a repeated include to define a different flavor of
-// yyFlexLexer, as discussed in the flex manual.
-#define yyFlexLexerOnce
-
-extern "C++" {
-
-class yyFlexLexer : public FlexLexer {
-public:
- // arg_yyin and arg_yyout default to the cin and cout, but we
- // only make that assignment when initializing in yylex().
- yyFlexLexer( FLEX_STD istream* arg_yyin = 0, FLEX_STD ostream* arg_yyout = 0 );
-
- virtual ~yyFlexLexer();
-
- void yy_switch_to_buffer( struct yy_buffer_state* new_buffer );
- struct yy_buffer_state* yy_create_buffer( FLEX_STD istream* s, int size );
- void yy_delete_buffer( struct yy_buffer_state* b );
- void yyrestart( FLEX_STD istream* s );
-
- void yypush_buffer_state( struct yy_buffer_state* new_buffer );
- void yypop_buffer_state();
-
- virtual int yylex();
- virtual void switch_streams( FLEX_STD istream* new_in, FLEX_STD ostream* new_out = 0 );
- virtual int yywrap();
-
-protected:
- virtual int LexerInput( char* buf, int max_size );
- virtual void LexerOutput( const char* buf, int size );
- virtual void LexerError( const char* msg );
-
- void yyunput( int c, char* buf_ptr );
- int yyinput();
-
- void yy_load_buffer_state();
- void yy_init_buffer( struct yy_buffer_state* b, FLEX_STD istream* s );
- void yy_flush_buffer( struct yy_buffer_state* b );
-
- int yy_start_stack_ptr;
- int yy_start_stack_depth;
- int* yy_start_stack;
-
- void yy_push_state( int new_state );
- void yy_pop_state();
- int yy_top_state();
-
- yy_state_type yy_get_previous_state();
- yy_state_type yy_try_NUL_trans( yy_state_type current_state );
- int yy_get_next_buffer();
-
- FLEX_STD istream* yyin; // input source for default LexerInput
- FLEX_STD ostream* yyout; // output sink for default LexerOutput
-
- // yy_hold_char holds the character lost when yytext is formed.
- char yy_hold_char;
-
- // Number of characters read into yy_ch_buf.
- int yy_n_chars;
-
- // Points to current character in buffer.
- char* yy_c_buf_p;
-
- int yy_init; // whether we need to initialize
- int yy_start; // start state number
-
- // Flag which is used to allow yywrap()'s to do buffer switches
- // instead of setting up a fresh yyin. A bit of a hack ...
- int yy_did_buffer_switch_on_eof;
-
-
- size_t yy_buffer_stack_top; /**< index of top of stack. */
- size_t yy_buffer_stack_max; /**< capacity of stack. */
- struct yy_buffer_state ** yy_buffer_stack; /**< Stack as an array. */
- void yyensure_buffer_stack(void);
-
- // The following are not always needed, but may be depending
- // on use of certain flex features (like REJECT or yymore()).
-
- yy_state_type yy_last_accepting_state;
- char* yy_last_accepting_cpos;
-
- yy_state_type* yy_state_buf;
- yy_state_type* yy_state_ptr;
-
- char* yy_full_match;
- int* yy_full_state;
- int yy_full_lp;
-
- int yy_lp;
- int yy_looking_for_trail_begin;
-
- int yy_more_flag;
- int yy_more_len;
- int yy_more_offset;
- int yy_prev_more_offset;
-};
-
-}
-
-#endif // yyFlexLexer || ! yyFlexLexerOnce
-
diff --git a/preprocessor/Makefile.am b/preprocessor/Makefile.am
index 38cf4082aa681e51c8f5f0c5d48836fab746eee7..6556fd445a4c7617d76e03865c1a5f2752ac87f0 100644
--- a/preprocessor/Makefile.am
+++ b/preprocessor/Makefile.am
@@ -1,6 +1,6 @@
SUBDIRS = macro
-BUILT_SOURCES = DynareBison.hh stack.hh position.hh location.hh DynareBison.cc DynareFlex.cc
+BUILT_SOURCES = DynareBison.hh stack.hh position.hh location.hh DynareBison.cc DynareFlex.cc FlexLexer.h
matlabdir = $(libdir)/matlab
@@ -47,7 +47,6 @@ dynare_m_SOURCES = \
DynareMain.cc \
DynareMain2.cc \
CodeInterpreter.hh \
- FlexLexer.h \
ExternalFunctionsTable.cc \
ExternalFunctionsTable.hh \
SteadyStateModel.hh \
@@ -60,8 +59,9 @@ dynare_m_CPPFLAGS = $(BOOST_CPPFLAGS) -I.
dynare_m_LDFLAGS = $(BOOST_LDFLAGS)
dynare_m_LDADD = macro/libmacro.a
-DynareFlex.cc: DynareFlex.ll
+DynareFlex.cc FlexLexer.h: DynareFlex.ll
$(LEX) -oDynareFlex.cc DynareFlex.ll
+ cp /usr/include/FlexLexer.h .
DynareBison.cc DynareBison.hh location.hh stack.hh position.hh: DynareBison.yy
$(YACC) -o DynareBison.cc DynareBison.yy
diff --git a/preprocessor/ModFile.cc b/preprocessor/ModFile.cc
index 4496b65bfa53068e9e0922aa500859b8ee2d3054..3f6b4a429937b63782876ee6dbe01bc5576df678 100644
--- a/preprocessor/ModFile.cc
+++ b/preprocessor/ModFile.cc
@@ -377,10 +377,11 @@ ModFile::computingPass(bool no_tmp_terms)
|| mod_file_struct.estimation_present || mod_file_struct.osr_present
|| mod_file_struct.ramsey_policy_present || mod_file_struct.identification_present)
static_model.set_cutoff_to_zero();
- if (mod_file_struct.identification_present)
- static_model.computingPass(global_eval_context, no_tmp_terms, true, block, byte_code);
- else
- static_model.computingPass(global_eval_context, no_tmp_terms, false, block, byte_code);
+
+ const bool static_hessian = mod_file_struct.identification_present
+ || mod_file_struct.estimation_analytic_derivation;
+ static_model.computingPass(global_eval_context, no_tmp_terms, static_hessian,
+ block, byte_code);
}
// Set things to compute for dynamic model
if (mod_file_struct.simul_present || mod_file_struct.check_present
@@ -402,7 +403,7 @@ ModFile::computingPass(bool no_tmp_terms)
exit(EXIT_FAILURE);
}
bool hessian = mod_file_struct.order_option >= 2 || mod_file_struct.identification_present || mod_file_struct.estimation_analytic_derivation;
- bool thirdDerivatives = mod_file_struct.order_option == 3;
+ bool thirdDerivatives = mod_file_struct.order_option == 3 || mod_file_struct.estimation_analytic_derivation;
bool paramsDerivatives = mod_file_struct.identification_present || mod_file_struct.estimation_analytic_derivation;
dynamic_model.computingPass(true, hessian, thirdDerivatives, paramsDerivatives, global_eval_context, no_tmp_terms, block, use_dll, byte_code);
}