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

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

time-shift
Michel Juillard 2012-03-10 18:21:14 +01:00
parent 80ca47d62a 4ba016cc16
commit 1fcf708b80
208 changed files with 10878 additions and 4600 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

13
configure.ac
View File

@ -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
@ -30,7 +30,7 @@ case ${host_os} in
*mingw32*)
# On mingw32, we don't want dynamic libgcc
# Note that static-libstdc++ is only supported since GCC 4.5 (but generates no error on older versions)
LDFLAGS="$LDFLAGS -static-libgcc -static-libstdc++"
LDFLAGS="$LDFLAGS -static-libgcc -static-libstdc++ -static-libgfortran"
have_windows="yes"
;;
*cygwin*)
@ -155,6 +155,15 @@ AC_CHECK_PROG([CWEAVE], [cweave], [cweave])
AM_CONDITIONAL([HAVE_CWEAVE], [test "x$CWEAVE" != "x"])
AC_PROG_F77
AC_F77_LIBRARY_LDFLAGS
# Hack to get static linking of libgfortran on MinGW
# (-static-libgfortran doesn't act on gcc/g++)
case ${host_os} in
*mingw32*)
FLIBS=`echo $FLIBS | sed 's/-lgfortran/-Wl,-Bstatic -lgfortran -Wl,-Bdynamic/'`
;;
esac
if test "x$F77" != "x"; then
AX_BLAS
AX_LAPACK

124
doc/dynare.texi
View File

@ -295,15 +295,15 @@ Stéphane Adjemian (Université du Maine, Gains and Cepremap), Houtan
Bastani (Cepremap), Michel Juillard (Banque de France), Frédéric Karamé
(Université d'Évry, Epee and Cepremap), Junior Maih (Norges Bank),
Ferhat Mihoubi (Université d'Évry, Epee and Cepremap), George Perendia,
Marco Ratto (JRC) and Sébastien Villemot (Cepremap and Paris School of
Economics). Increasingly, the developer base is expanding, as tools
developed by researchers outside of Cepremap are integrated into
Dynare. Financial support is provided by Cepremap, Banque de France and
DSGE-net (an international research network for DSGE modeling). The
Dynare project also received funding through the Seventh Framework
Programme for Research (FP7) of the European Commission's Socio-economic
Sciences and Humanities (SSH) Program from October 2008 to September
2011 under grant agreement SSH-CT-2009-225149.
Johannes Pfeifer, Marco Ratto (JRC) and Sébastien Villemot (Cepremap and
Paris School of Economics). Increasingly, the developer base is
expanding, as tools developed by researchers outside of Cepremap are
integrated into Dynare. Financial support is provided by Cepremap,
Banque de France and DSGE-net (an international research network for
DSGE modeling). The Dynare project also received funding through the
Seventh Framework Programme for Research (FP7) of the European
Commission's Socio-economic Sciences and Humanities (SSH) Program from
October 2008 to September 2011 under grant agreement SSH-CT-2009-225149.
Interaction between developers and users of Dynare is central to the
project. A @uref{http://www.dynare.org/phpBB3, web forum} is available
@ -3888,6 +3888,9 @@ Uses Dynare implementation of the Nelder-Mead simplex based optimization
routine (generally more efficient than the MATLAB or Octave implementation
available with @code{mode_compute=7})
@item 9
Uses the CMA-ES (Covariance Matrix Adaptation Evolution Strategy) algorithm, an evolutionary algorithm for difficult non-linear non-convex optimization
@item @var{FUNCTION_NAME}
It is also possible to give a @var{FUNCTION_NAME} to this option,
instead of an @var{INTEGER}. In that case, Dynare takes the return
@ -3988,6 +3991,7 @@ and graphics that can be later directly included in LaTeX files (not
yet implemented)
@item kalman_algo = @var{INTEGER}
@anchor{kalman_algo}
@dots{}
@item kalman_tol = @var{DOUBLE}
@ -5263,6 +5267,9 @@ Critical value for correlation @math{\rho}: plot couples of parmaters with
@item mode_file = @var{FILENAME}
@xref{mode_file}.
@item kalman_algo = @var{INTEGER}
@xref{kalman_algo}.
@end table
@customhead{Identification Analysis Options}
@table @code
@ -5290,6 +5297,9 @@ for identification analysis. Default: @code{0}
@item ar = @var{INTEGER}
Maximum number of lags for moments in identification analysis. Default: @code{1}
@item lik_init = @var{INTEGER}
@xref{lik_init}.
@end table
@end deffn
@ -5357,6 +5367,9 @@ Specify the parameter set to use. Default: @code{prior_mean}
@item lik_init = @var{INTEGER}
@xref{lik_init}.
@item kalman_algo = @var{INTEGER}
@xref{kalman_algo}.
@end table
@end deffn
@ -5841,20 +5854,17 @@ model. Output @code{.eps} files are contained in
@xref{simulation_file_tag}.
@item horizon = @var{INTEGER}
The forecast horizon. Default: @code{12}
@anchor{horizon} The forecast horizon. Default: @code{12}
@item filtered_probabilities
@anchor{filtered_probabilities} Uses filtered probabilities at the end
of the sample as initial conditions for regime probabilities. Default:
@code{off}
@item no_error_bands
@anchor{no_error_bands} Do not output error bands. Default: @code{off}
(@i{i.e.} output error bands)
of the sample as initial conditions for regime probabilities. Only one
of @code{filtered_probabilities}, @code{regime} and @code{regimes} may
be passed. Default: @code{off}
@item error_band_percentiles = [@var{DOUBLE1} @dots{}]
@anchor{error_band_percentiles} The percentiles to compute. Default:
@code{[0.16 0.50 0.84]}. If @code{no_error_bands} is passed, the default
@code{[0.16 0.50 0.84]}. If @code{median} is passed, the default
is @code{[0.5]}
@item shock_draws = @var{INTEGER}
@ -5873,11 +5883,25 @@ draws in posterior draws file are used. Default: @code{1}
@anchor{free_parameters} A vector of free parameters to initialize theta
of the model. Default: use estimated parameters
@item median
@anchor{median}
@item parameter_uncertainty
@anchor{parameter_uncertainty} Calculate IRFs under parameter
uncertainty. Requires that @command{ms_simulation} has been
run. Default: @code{off}
A shortcut to setting @code{error_band_percentiles=[0.5]}. Default:
@code{off}
@item regime = @var{INTEGER}
@anchor{regime} Given the data and model parameters, what is the ergodic
probability of being in the specified regime. Only one of
@code{filtered_probabilities}, @code{regime} and @code{regimes} may be
passed. Default: @code{off}
@item regimes
@anchor{regimes} Describes the evolution of regimes. Only one of
@code{filtered_probabilities}, @code{regime} and @code{regimes} may be
passed. Default: @code{off}
@item median
@anchor{median} A shortcut to setting
@code{error_band_percentiles=[0.5]}. Default: @code{off}
@end table
@ -5909,9 +5933,6 @@ while data files are contained in @code{<output_file_tag/Forecast>}.
@item data_obs_nbr = @var{INTEGER}
The number of data points included in the output. Default: @code{0}
@item no_error_bands
@xref{no_error_bands}.
@item error_band_percentiles = [@var{DOUBLE1} @dots{}]
@xref{error_band_percentiles}.
@ -5927,6 +5948,16 @@ The number of data points included in the output. Default: @code{0}
@item free_parameters = @var{NUMERICAL_VECTOR}
@xref{free_parameters}.
@item parameter_uncertainty
@xref{parameter_uncertainty}.
@item regime = @var{INTEGER}
@xref{regime}.
@item regimes
@xref{regimes}.
@item median
@xref{median}.
@ -5959,11 +5990,15 @@ are contained in @code{<output_file_tag/Variance_Decomposition>}.
@item simulation_file_tag = @var{FILENAME}
@xref{simulation_file_tag}.
@item horizon = @var{INTEGER}
@xref{horizon}.
@item filtered_probabilities
@xref{filtered_probabilities}.
@item no_error_bands
@xref{no_error_bands}.
Do not output percentile error bands (@i{i.e.} compute mean). Default:
@code{off} (@i{i.e.} output error bands)
@item error_band_percentiles = [@var{DOUBLE1} @dots{}]
@xref{error_band_percentiles}.
@ -5980,9 +6015,15 @@ are contained in @code{<output_file_tag/Variance_Decomposition>}.
@item free_parameters = @var{NUMERICAL_VECTOR}
@xref{free_parameters}.
@item median
@item parameter_uncertainty
@xref{parameter_uncertainty}.
@xref{median}.
@item regime = @var{INTEGER}
@xref{regime}.
@item regimes
@xref{regimes}.
@end table
@ -6075,7 +6116,7 @@ line. The main directives are:
@item
@code{@@#define}, for defining a macro-processor variable,
@item
@code{@@#if}, @code{@@#then}, @code{@@#else}, @code{@@#endif} for
@code{@@#if}, @code{@@#ifdef}, @code{@@#else}, @code{@@#endif} for
conditional statements,
@item
@code{@@#for}, @code{@@#endfor} for constructing loops.
@ -6218,13 +6259,15 @@ end;
@end deffn
@deffn {Macro directive} @@#if @var{MACRO_EXPRESSION}
@deffnx {Macro directive} @@#ifdef @var{MACRO_VARIABLE}
@deffnx {Macro directive} @@#else
@deffnx {Macro directive} @@#endif
Conditional inclusion of some part of the @file{.mod} file.
The lines between @code{@@#if} and the next @code{@@#else} or
@code{@@#end} is executed only if the condition evaluates to a
non-null integer. The @code{@@#else} branch is optional and, if
present, is only evaluated if the condition evaluates to @code{0}.
The lines between @code{@@#if} or @code{@@#ifdef} and the next
@code{@@#else} or @code{@@#endif} is executed only if the condition
evaluates to a non-null integer. The @code{@@#else} branch is optional
and, if present, is only evaluated if the condition evaluates to
@code{0}.
@examplehead
@ -6242,6 +6285,23 @@ model;
end;
@end example
@examplehead
Choose between two alternative monetary policy rules using a
macro-variable. As @code{linear_mon_pol} was not previously defined in
this example, the second equation will be chosen:
@example
model;
@@#ifdef linear_mon_pol
i = w*i(-1) + (1-w)*i_ss + w2*(pie-piestar);
@@#else
i = i(-1)^w * i_ss^(1-w) * (pie/piestar)^w2;
@@#endif
...
end;
@end example
@end deffn
@deffn {Macro directive} @@#for @var{MACRO_VARIABLE} in @var{MACRO_EXPRESSION}

10
doc/gsa/gsa.tex
View File

@ -8,6 +8,8 @@
\usepackage{psfrag}
\usepackage{setspace}
\usepackage{rotating}
\usepackage{hyperref}
\hypersetup{breaklinks=true,pagecolor=white,colorlinks=true,linkcolor=blue,citecolor=blue,urlcolor=blue}
%\singlespacing (interlinea singola)
%\onehalfspacing (interlinea 1,5)
%\doublespacing (interlinea doppia)
@ -20,7 +22,13 @@
\begin{document}
% ----------------------------------------------------------------
\title{Sensitivity Analysis Toolbox for DYNARE}%
\title{Sensitivity Analysis Toolbox for DYNARE\thanks{Copyright \copyright~2012 Dynare
Team. Permission is granted to copy, distribute and/or modify
this document under the terms of the GNU Free Documentation
License, Version 1.3 or any later version published by the Free
Software Foundation; with no Invariant Sections, no Front-Cover
Texts, and no Back-Cover Texts. A copy of the license can be found
at: \url{http://www.gnu.org/licenses/fdl.txt}}}
\author{Marco Ratto\\
European Commission, Joint Research Centre \\

4
doc/macroprocessor/macroprocessor.tex
View File

@ -65,7 +65,7 @@
\begin{itemize}
\item file inclusion
\item loops (\textit{for} structure)
\item conditional inclusion (\textit{if/then/else} structures)
\item conditional inclusion (\textit{if/else} structures)
\item expression substitution
\end{itemize}
\item Implemented in Dynare starting from version 4.0
@ -95,7 +95,7 @@
\begin{itemize}
\item file inclusion: \verb+@#include+
\item definition a variable of the macro-processor: \verb+@#define+
\item conditional statements (\verb+@#if/@#then/@#else/@#endif+)
\item conditional statements (\verb+@#if/@#else/@#endif+)
\item loop statements (\verb+@#for/@#endfor+)
\end{itemize}
\item In most cases, directives occupy exactly one line of text. In case of need, two anti-slashes (\verb+\\+) at the end of the line indicates that the directive is continued on the next line.

424
license.txt
View File

@ -1,31 +1,21 @@
Format: http://dep.debian.net/deps/dep5/
Format: http://www.debian.org/doc/packaging-manuals/copyright-format/1.0/
Upstream-Name: Dynare
Upstream-Contact: Dynare Team, whose members in 2011 are:
Stéphane Adjemian <stephane.adjemian@ens.fr>
Houtan Bastani <houtan.bastani@ens.fr>
Michel Juillard <michel.juillard@mjui.fr>
Junior Maih <junior.maih@gmail.com>
Ferhat Mihoubi <fmihoubi@univ-evry.fr>
George Perendia <george@perendia.orangehome.co.uk>
Marco Ratto <marco.ratto@jrc.ec.europa.eu>
Sébastien Villemot <sebastien.villemot@ens.fr>
Upstream-Contact: Dynare Team, whose members in 2012 are:
Stéphane Adjemian <stephane.adjemian@ens.fr>
Houtan Bastani <houtan.bastani@ens.fr>
Michel Juillard <michel.juillard@mjui.fr>
Frédéric Karamé <frederic.karame@univ-evry.fr>
Junior Maih <junior.maih@gmail.com>
Ferhat Mihoubi <fmihoubi@univ-evry.fr>
George Perendia <george@perendia.orangehome.co.uk>
Johannes Pfeifer <jpfeifer@gmx.de>
Marco Ratto <marco.ratto@jrc.ec.europa.eu>
Sébastien Villemot <sebastien.villemot@ens.fr>
Source: http://www.dynare.org
Files: *
Copyright: 1996-2011 Dynare Team
Copyright: 1996-2012 Dynare Team
License: GPL-3+
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/>.
Files: matlab/AIM/SP*
Copyright: public-domain
@ -49,41 +39,18 @@ License: public-domain
Journal of Economic Dynamics and Control, 2010, vol. 34, issue 3,
pages 472-489
Files: matlab/bfgsi.m
Copyright: 1993-2009 Christopher Sims
License: GPL-3+
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/>.
Files: matlab/csolve.m matlab/csminit1.m matlab/numgrad2.m matlab/numgrad3.m
matlab/numgrad5.m matlab/csminwel.m matlab/bvar_density.m
matlab/bvar_toolbox.m matlab/partial_information/PI_gensys.m
matlab/qzswitch.m matlab/qzdiv.m
Files: matlab/bfgsi.m matlab/csolve.m matlab/csminit1.m matlab/numgrad2.m
matlab/numgrad3.m matlab/numgrad5.m matlab/csminwel.m matlab/bvar_density.m
matlab/bvar_toolbox.m matlab/partial_information/PI_gensys.m matlab/qzswitch.m
matlab/qzdiv.m
Copyright: 1993-2009 Christopher Sims
2006-2011 Dynare Team
License: GPL-3+
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/>.
Files: matlab/cmaes.m
Copyright: 2001-2012 Nikolaus Hansen
2012 Dynare Team
License: GPL-3+
Files: matlab/missing/stats/normpdf.m matlab/missing/stats/gamcdf.m
matlab/missing/stats/common_size.m matlab/missing/stats/chi2inv.m
@ -95,173 +62,54 @@ Files: matlab/missing/stats/normpdf.m matlab/missing/stats/gamcdf.m
Copyright: 1995-2007 Kurt Hornik
2008-2009 Dynare Team
License: GPL-3+
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/>.
Files: matlab/missing/bicgstab/bicgstab.m
Copyright: 2008 Radek Salac
License: GPL-3+
This file is part of Octave.
.
Octave 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.
.
Octave 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 Octave; see the file COPYING. If not, see
<http://www.gnu.org/licenses/>.
Files: doc/dynare.texi doc/*.tex doc/*.svg doc/*.dia doc/*.pdf doc/*.bib
Copyright: 1996-2011 Dynare Team
Copyright: 1996-2012 Dynare Team
License: GFDL-NIV-1.3+
Permission is granted to copy, distribute and/or modify this document
under the terms of the GNU Free Documentation License, Version 1.3 or
any later version published by the Free Software Foundation; with no
Invariant Sections, no Front-Cover Texts, and no Back-Cover Texts.
.
A copy of the license can be found at <http://www.gnu.org/licenses/fdl.txt>
Files: doc/userguide/*.tex doc/userguide/*.bib doc/userguide/*.pdf
Copyright: 2007-2011 Tommaso Mancini Griffoli
License: GFDL-NIV-1.3+
Permission is granted to copy, distribute and/or modify this document
under the terms of the GNU Free Documentation License, Version 1.3 or
any later version published by the Free Software Foundation; with no
Invariant Sections, no Front-Cover Texts, and no Back-Cover Texts.
.
A copy of the license can be found at <http://www.gnu.org/licenses/fdl.txt>
Files: doc/dr.tex doc/bvar_a_la_sims.tex
Copyright: 2007-2011 Sébastien Villemot
License: GFDL-NIV-1.3+
Permission is granted to copy, distribute and/or modify this document
under the terms of the GNU Free Documentation License, Version 1.3 or
any later version published by the Free Software Foundation; with no
Invariant Sections, no Front-Cover Texts, and no Back-Cover Texts.
.
A copy of the license can be found at <http://www.gnu.org/licenses/fdl.txt>
Files: dynare++/*.cweb dynare++/*.hweb dynare++/*.cpp dynare++/*.h
dynare++/*.tex dynare++/*.mod dynare++/*.m dynare++/*.web dynare++/*.lex
dynare++/*.y
Copyright: 2004-2011 Ondra Kamenik
License: GPL-3+
This program 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.
.
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 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/>.
Files: dynare++/utils/*.cpp dynare++/utils/*.h dynare++/parser/*.cpp
dynare++/parser/*.h dynare++/parser/*.lex dynare++/parser/*.y
Copyright: 2004-2011 Ondra Kamenik
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: m4/ax_blas.m4 m4/ax_lapack.m4
Copyright: 2008 Steven G. Johnson <stevenj@alum.mit.edu>
License: GPL-3+ with Autoconf exception
This program 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.
.
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 General
Public License for more details.
.
You should have received a copy of the GNU General Public License along
with this program. If not, see <http://www.gnu.org/licenses/>.
.
As a special exception, the respective Autoconf Macro's copyright owner
gives unlimited permission to copy, distribute and modify the configure
scripts that are the output of Autoconf when processing the Macro. You
need not follow the terms of the GNU General Public License when using
or distributing such scripts, even though portions of the text of the
Macro appear in them. The GNU General Public License (GPL) does govern
all other use of the material that constitutes the Autoconf Macro.
.
This special exception to the GPL applies to versions of the Autoconf
Macro released by the Autoconf Archive. When you make and distribute a
modified version of the Autoconf Macro, you may extend this special
exception to the GPL to apply to your modified version as well.
Files: m4/ax_pthread.m4
Copyright: 2008 Steven G. Johnson <stevenj@alum.mit.edu>
2010 Dynare Team
License: GPL-3+ with Autoconf exception
This program 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.
.
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 General
Public License for more details.
.
You should have received a copy of the GNU General Public License along
with this program. If not, see <http://www.gnu.org/licenses/>.
.
As a special exception, the respective Autoconf Macro's copyright owner
gives unlimited permission to copy, distribute and modify the configure
scripts that are the output of Autoconf when processing the Macro. You
need not follow the terms of the GNU General Public License when using
or distributing such scripts, even though portions of the text of the
Macro appear in them. The GNU General Public License (GPL) does govern
all other use of the material that constitutes the Autoconf Macro.
.
This special exception to the GPL applies to versions of the Autoconf
Macro released by the Autoconf Archive. When you make and distribute a
modified version of the Autoconf Macro, you may extend this special
exception to the GPL to apply to your modified version as well.
Files: m4/ax_boost_base.m4
Copyright: 2008 Thomas Porschberg <thomas@randspringer.de>
2009 Dynare Team
License:
License: other
Copying and distribution of this file, with or without modification, are
permitted in any medium without royalty provided the copyright notice
and this notice are preserved.
Files: m4/ax_compare_version.m4
Copyright: 2008 Tim Toolan <toolan@ele.uri.edu>
License:
License: other
Copying and distribution of this file, with or without modification, are
permitted in any medium without royalty provided the copyright notice
and this notice are preserved.
@ -270,22 +118,72 @@ Files: m4/ax_latex_bibtex_test.m4 m4/ax_latex_class.m4 m4/ax_tex_test.m4
Copyright: 2008 Boretti Mathieu <boretti@eig.unige.ch>
2009 Dynare Team
License: LGPL-2.1+
This library 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 2.1 of the License, or (at
your option) any later version.
.
This library 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 library. If not, see <http://www.gnu.org/licenses/>.
Files: m4/ax_matlab_arch.m4 m4/ax_matlab.m4 m4/ax_mexext.m4
Copyright: 2002-2003 Ralph Schleicher
2009 Dynare Team
License: GPL-2+ with Autoconf exception
Files: dynare.el
Copyright: 2010 Yannick Kalantzis
License: GPL-3+
Files: mex/sources/libslicot/*
Copyright: 2002-2009 NICONET e.V.
License: GPL-2+
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+
Files: ms-sbvar/utilities_dw/*
Copyright: 1996-2011 Daniel Waggoner
License: GPL-3+
Files: ms-sbvar/switch_dw/*
Copyright: 1996-2011 Daniel Waggoner and Tao Zha
License: GPL-3+
Files: ms-sbvar/switch_dw/state_space/sbvar/dw_csminwel.c
state_space/sbvar/dw_csminwel.h
Copyright: 1996 Christopher Sims
2003 Karibzhanov, Waggoner and Zha
License: GPL-3+
Files: matlab/ms-sbvar/cstz/*
Copyright: 1993-2011 Tao Zha
License: GPL-3+
Files: matlab/ms-sbvar/cstz/bfgsi.m matlab/ms-sbvar/cstz/csminit.m
matlab/ms-sbvar/cstz/csminwel.m
Copyright: 1993-2011 Tao Zha and Christopher Sims
License: GPL-3+
License: GFDL-NIV-1.3+
Permission is granted to copy, distribute and/or modify this document
under the terms of the GNU Free Documentation License, Version 1.3 or
any later version published by the Free Software Foundation; with no
Invariant Sections, no Front-Cover Texts, and no Back-Cover Texts.
.
A copy of the license can be found at <http://www.gnu.org/licenses/fdl.txt>
License: GPL-2+
This program 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 2 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 General Public License for more details.
.
You should have received a copy of the GNU General Public License
along with this program. If not, see
<http://www.gnu.org/licenses/>.
License: GPL-2+ with Autoconf exception
This program is free software; you can redistribute it and/or
modify it under the terms of the GNU General Public License as
@ -306,44 +204,61 @@ License: GPL-2+ with Autoconf exception
it under the same distribution terms that you use for the rest
of that program.
Files: dynare.el
Copyright: 2010 Yannick Kalantzis
License: GPL-3+
This program is free software; you can redistribute it and/or modify
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.
.
This program is distributed in the hope that it will be useful,
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 this program. If not, see
<http://www.gnu.org/licenses/>.
along with Dynare. If not, see <http://www.gnu.org/licenses/>.
Files: mex/sources/libslicot/*
Copyright: 2002-2009 NICONET e.V.
License: GPL-2+
This program 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 2 of
the License, or (at your option) any later version.
License: GPL-3+ with Autoconf exception
This program 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.
.
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 General Public License for more details.
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 General
Public License for more details.
.
You should have received a copy of the GNU General Public License
along with this program. If not, see
<http://www.gnu.org/licenses/>.
You should have received a copy of the GNU General Public License along
with this program. If not, see <http://www.gnu.org/licenses/>.
.
As a special exception, the respective Autoconf Macro's copyright owner
gives unlimited permission to copy, distribute and modify the configure
scripts that are the output of Autoconf when processing the Macro. You
need not follow the terms of the GNU General Public License when using
or distributing such scripts, even though portions of the text of the
Macro appear in them. The GNU General Public License (GPL) does govern
all other use of the material that constitutes the Autoconf Macro.
.
This special exception to the GPL applies to versions of the Autoconf
Macro released by the Autoconf Archive. When you make and distribute a
modified version of the Autoconf Macro, you may extend this special
exception to the GPL to apply to your modified version as well.
License: LGPL-2.1+
This library 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 2.1 of the License, or (at
your option) any later version.
.
This library 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 library. 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
@ -357,86 +272,3 @@ License: LGPL-3+
.
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+
This 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.
.
It 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.
.
If you did not received a copy of the GNU General Public License
with this software, see <http://www.gnu.org/licenses/>.
Files: ms-sbvar/switch_dw/*
Copyright: 1996-2011 Daniel Waggoner and Tao Zha
License: GPL-3+
This 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.
.
It 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.
.
If you did not received a copy of the GNU General Public License
with this software, see <http://www.gnu.org/licenses/>.
Files: ms-sbvar/switch_dw/state_space/sbvar/dw_csminwel.c
state_space/sbvar/dw_csminwel.h
Copyright: 1996 Christopher Sims
2003 Karibzhanov, Waggoner and Zha
License: GPL-3+
This 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.
.
It 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.
.
If you did not received a copy of the GNU General Public License
with this software, see <http://www.gnu.org/licenses/>.
Files: matlab/ms-sbvar/cstz/*
Copyright: 1993-2011 Tao Zha
License: GPL-3+
This 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.
.
It 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.
.
If you did not received a copy of the GNU General Public License
with this software, see <http://www.gnu.org/licenses/>.
Files: matlab/ms-sbvar/cstz/bfgsi.m matlab/ms-sbvar/cstz/csminit.m
matlab/ms-sbvar/cstz/csminwel.m
Copyright: 1993-2011 Tao Zha and Christopher Sims
License: GPL-3+
This 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.
.
It 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.
.
If you did not received a copy of the GNU General Public License
with this software, see <http://www.gnu.org/licenses/>.

3055
matlab/cmaes.m Normal file
View File

File diff suppressed because it is too large Load Diff

8
matlab/distributions/inverse_gamma_specification.m
View File

@ -44,7 +44,7 @@ function [s,nu] = inverse_gamma_specification(mu,sigma,type,use_fzero_flag)
%! @end deftypefn
%@eod:
% Copyright (C) 2003-2011 Dynare Team
% Copyright (C) 2003-2012 Dynare Team
%
% This file is part of Dynare.
%
@ -99,7 +99,7 @@ elseif type == 1; % Inverse Gamma 1
end
end
% Solve for nu using the secant method.
while abs(nu2-nu1) > 1e-8
while abs(nu2/nu1-1) > 1e-14
if err > 0
nu1 = nu;
if nu < nu2
@ -117,10 +117,10 @@ elseif type == 1; % Inverse Gamma 1
end
s = (sigma2+mu2)*(nu-2);
if check_solution_flag
if abs(mu-sqrt(s/2)*gamma((nu-1)/2)/gamma(nu/2))>1e-9
if abs(log(mu)-log(sqrt(s/2))-gammaln((nu-1)/2)+gammaln(nu/2))>1e-7
error('inverse_gamma_specification:: Failed in solving for the hyperparameters!');
end
if abs(sigma-sqrt(s/(nu-2)-mu^2))>1e-9
if abs(sigma-sqrt(s/(nu-2)-mu^2))>1e-7
error('inverse_gamma_specification:: Failed in solving for the hyperparameters!');
end
end

14
matlab/dr1.m
View File

@ -1,4 +1,4 @@
function [dr,info,M_,options_,oo_] = dr1(dr,task,M_,options_,oo_)
function [dr,info] = dr1(dr,task,M_,options_,oo_)
% function [dr,info,M_,options_,oo_] = dr1(dr,task,M_,options_,oo_)
% computes the reduced form solution of a rational expectation model (first or second order
% approximation of the stochastic model around the deterministic steady state).
@ -21,9 +21,6 @@ function [dr,info,M_,options_,oo_] = dr1(dr,task,M_,options_,oo_)
% indeterminacy.
% info=5: BK rank condition not satisfied.
% info=6: The jacobian matrix evaluated at the steady state is complex.
% M_ [matlab structure]
% options_ [matlab structure]
% oo_ [matlab structure]
%
% ALGORITHM
% ...
@ -61,7 +58,6 @@ end
if options_.k_order_solver;
dr = set_state_space(dr,M_);
[dr,info] = k_order_pert(dr,M_,options_,oo_);
oo_.dr = dr;
return;
end
@ -72,20 +68,12 @@ iyv = iyv(:);
iyr0 = find(iyv) ;
it_ = M_.maximum_lag + 1 ;
if M_.exo_nbr == 0
oo_.exo_steady_state = [] ;
end
klen = M_.maximum_lag + M_.maximum_lead + 1;
iyv = lead_lag_incidence';
iyv = iyv(:);
iyr0 = find(iyv) ;
it_ = M_.maximum_lag + 1 ;
if M_.exo_nbr == 0
oo_.exo_steady_state = [] ;
end
it_ = M_.maximum_lag + 1;
z = repmat(dr.ys,1,klen);
if ~options_.bytecode

6
matlab/dr_block.m
View File

@ -567,7 +567,9 @@ for i = 1:Size;
if block_type == 5
vghx_other = - inv(kron(eye(size(D_,2)), A_) + kron(C_', B_)) * vec(D_);
ghx_other = reshape(vghx_other, size(D_,1), size(D_,2));
else
elseif options_.sylvester_fp == 1
ghx_other = gensylv_fp(A_, B_, C_, D_, i);
else
[err, ghx_other] = gensylv(1, A_, B_, C_, -D_);
end;
if options_.aim_solver ~= 1 && options_.use_qzdiv
@ -650,7 +652,7 @@ for i = 1:Size;
end
end;
if task ~=1
if (maximum_lag > 0 && n_pred > 0)
if (maximum_lag > 0 && (n_pred > 0 || n_both > 0))
sorted_col_dr_ghx = M_.block_structure.block(i).sorted_col_dr_ghx;
dr.ghx(endo, sorted_col_dr_ghx) = dr.ghx(endo, sorted_col_dr_ghx) + ghx;
data(i).ghx = ghx;

116
matlab/dsge_likelihood.m
View File

@ -127,17 +127,21 @@ function [fval,exit_flag,ys,trend_coeff,info,Model,DynareOptions,BayesInfo,Dynar
% AUTHOR(S) stephane DOT adjemian AT univ DASH lemans DOT FR
% Declaration of the penalty as a persistent variable.
persistent penalty
% Initialization of the persistent variable.
if ~nargin || isempty(penalty)
penalty = 1e8;
if ~nargin, return, end
end
if nargin==1
penalty = xparam1;
return
end
% 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
% Initialization of the returned variables and others...
fval = [];
@ -146,6 +150,21 @@ trend_coeff = [];
exit_flag = 1;
info = 0;
singularity_flag = 0;
DLIK = [];
AHess = [];
if DynareOptions.estimation_dll
[fval,exit_flag,ys,trend_coeff,info,params,H,Q] ...
= logposterior(xparam1,DynareDataset, DynareOptions,Model, ...
EstimatedParameters,BayesInfo,DynareResults);
Model.params = params;
if ~isequal(Model.H,0)
Model.H = H;
end
Model.Sigma_e = Q;
DynareResults.dr.ys = ys;
return
end
% Set flag related to analytical derivatives.
if nargout > 9
@ -209,7 +228,7 @@ if EstimatedParameters.ncx
a = diag(eig(Q));
k = find(a < 0);
if k > 0
fval = BayesInfo.penalty+sum(-a(k));
fval = penalty+sum(-a(k));
exit_flag = 0;
info = 43;
return
@ -233,7 +252,7 @@ if EstimatedParameters.ncn
a = diag(eig(H));
k = find(a < 0);
if k > 0
fval = BayesInfo.penalty+sum(-a(k));
fval = penalty+sum(-a(k));
exit_flag = 0;
info = 44;
return
@ -339,13 +358,19 @@ end
diffuse_periods = 0;
correlated_errors_have_been_checked = 0;
singular_diffuse_filter = 0;
switch DynareOptions.lik_init
case 1% Standard initialization with the steady state of the state equation.
if kalman_algo~=2
% Use standard kalman filter except if the univariate filter is explicitely choosen.
kalman_algo = 1;
end
Pstar = lyapunov_symm(T,R*Q*R',DynareOptions.qz_criterium,DynareOptions.lyapunov_complex_threshold);
if DynareOptions.lyapunov_fp == 1
Pstar = lyapunov_symm(T,Q,DynareOptions.qz_criterium,DynareOptions.lyapunov_complex_threshold, 4, R);
else
Pstar = lyapunov_symm(T,R*Q*R',DynareOptions.qz_criterium,DynareOptions.lyapunov_complex_threshold);
end;
Pinf = [];
a = zeros(mm,1);
Zflag = 0;
@ -359,10 +384,14 @@ switch DynareOptions.lik_init
a = zeros(mm,1);
Zflag = 0;
case 3% Diffuse Kalman filter (Durbin and Koopman)
if kalman_algo ~= 4
% Use standard kalman filter except if the univariate filter is explicitely choosen.
if kalman_algo == 0
kalman_algo = 3;
elseif ~((kalman_algo == 3) || (kalman_algo == 4))
error(['diffuse filter: options_.kalman_algo can only be equal ' ...
'to 0 (default), 3 or 4'])
end
[Z,T,R,QT,Pstar,Pinf] = schur_statespace_transformation(Z,T,R,Q,DynareOptions.qz_criterium);
Zflag = 1;
% Run diffuse kalman filter on first periods.
@ -383,39 +412,39 @@ switch DynareOptions.lik_init
diffuse_periods = length(tmp);
if isinf(dLIK)
% Go to univariate diffuse filter if singularity problem.
kalman_algo = 4;
singularity_flag = 1;
singular_diffuse_filter = 1;
end
end
if (kalman_algo==4)
if singular_diffuse_filter || (kalman_algo==4)
% Univariate Diffuse Kalman Filter
if singularity_flag
if isequal(H,0)
H = zeros(nobs,1);
mmm = mm;
if isequal(H,0)
H1 = zeros(nobs,1);
mmm = mm;
else
if all(all(abs(H-diag(diag(H)))<1e-14))% ie, the covariance matrix is diagonal...
H1 = diag(H);
mmm = mm;
else
if all(all(abs(H-diag(diag(H)))<1e-14))% ie, the covariance matrix is diagonal...
H = diag(H);
mmm = mm;
else
Z = [Z, eye(pp)];
T = blkdiag(T,zeros(pp));
Q = blkdiag(Q,H);
R = blkdiag(R,eye(pp));
Pstar = blkdiag(Pstar,H);
Pinf = blckdiag(Pinf,zeros(pp));
H = zeros(nobs,1);
mmm = mm+pp;
end
Z = [Z, eye(pp)];
T = blkdiag(T,zeros(pp));
Q = blkdiag(Q,H);
R = blkdiag(R,eye(pp));
Pstar = blkdiag(Pstar,H);
Pinf = blckdiag(Pinf,zeros(pp));
H1 = zeros(nobs,1);
mmm = mm+pp;
end
% no need to test again for correlation elements
singularity_flag = 0;
end
[dLIK,tmp,a,Pstar] = univariate_kalman_filter_d(DynareDataset.missing.aindex,DynareDataset.missing.number_of_observations,DynareDataset.missing.no_more_missing_observations, ...
Y, 1, size(Y,2), ...
zeros(mmm,1), Pinf, Pstar, ...
kalman_tol, riccati_tol, DynareOptions.presample, ...
T,R,Q,H,Z,mmm,pp,rr);
% no need to test again for correlation elements
correlated_errors_have_been_checked = 1;
[dLIK,tmp,a,Pstar] = univariate_kalman_filter_d(DynareDataset.missing.aindex,...
DynareDataset.missing.number_of_observations,...
DynareDataset.missing.no_more_missing_observations, ...
Y, 1, size(Y,2), ...
zeros(mmm,1), Pinf, Pstar, ...
kalman_tol, riccati_tol, DynareOptions.presample, ...
T,R,Q,H1,Z,mmm,pp,rr);
diffuse_periods = length(tmp);
end
case 4% Start from the solution of the Riccati equation.
@ -586,7 +615,6 @@ if ((kalman_algo==1) || (kalman_algo==3))% Multivariate Kalman Filter
else
kalman_algo = 4;
end
singularity_flag = 1;
else
if DynareOptions.lik_init==3
LIK = LIK + dLIK;
@ -594,10 +622,10 @@ if ((kalman_algo==1) || (kalman_algo==3))% Multivariate Kalman Filter
end
end
if ( singularity_flag || (kalman_algo==2) || (kalman_algo==4) )
if (kalman_algo==2) || (kalman_algo==4)
% Univariate Kalman Filter
% resetting measurement error covariance matrix when necessary %
if singularity_flag
if ~correlated_errors_have_been_checked
if isequal(H,0)
H = zeros(nobs,1);
mmm = mm;

101
matlab/DsgeLikelihood_hh.m → matlab/dsge_likelihood_hh.m
View File

@ -1,4 +1,4 @@
function [fval,llik,cost_flag,ys,trend_coeff,info] = DsgeLikelihood_hh(xparam1,DynareDataset,DynareOptions,Model,EstimatedParameters,BayesInfo,DynareResults)
function [fval,llik,cost_flag,ys,trend_coeff,info] = dsge_likelihood_hh(xparam1,DynareDataset,DynareOptions,Model,EstimatedParameters,BayesInfo,DynareResults)
% function [fval,llik,cost_flag,ys,trend_coeff,info] = DsgeLikelihood_hh(xparam1,DynareDataset,DynareOptions,Model,EstimatedParameters,BayesInfo,DynareResults)
% Evaluates the posterior kernel of a dsge model.
%
@ -41,17 +41,20 @@ function [fval,llik,cost_flag,ys,trend_coeff,info] = DsgeLikelihood_hh(xparam1,D
% Declaration of the penalty as a persistent variable.
persistent penalty
% Initialization of the persistent variable.
if ~nargin || isempty(penalty)
penalty = 1e8;
if ~nargin, return, end
end
if nargin==1
penalty = xparam1;
return
end
% 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_hh() is called by an newrat() 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_hh() 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
% Initialization of the returned variables and others...
fval = [];
@ -121,7 +124,7 @@ if EstimatedParameters.ncx
a = diag(eig(Q));
k = find(a < 0);
if k > 0
fval = BayesInfo.penalty+sum(-a(k));
fval = penalty+sum(-a(k));
exit_flag = 0;
info = 43;
return
@ -145,7 +148,7 @@ if EstimatedParameters.ncn
a = diag(eig(H));
k = find(a < 0);
if k > 0
fval = BayesInfo.penalty+sum(-a(k));
fval = penalty+sum(-a(k));
exit_flag = 0;
info = 44;
return
@ -251,6 +254,8 @@ end
diffuse_periods = 0;
correlated_errors_have_been_checked = 0;
singular_diffuse_filter = 0;
switch DynareOptions.lik_init
case 1% Standard initialization with the steady state of the state equation.
if kalman_algo~=2
@ -271,10 +276,14 @@ switch DynareOptions.lik_init
a = zeros(mm,1);
Zflag = 0;
case 3% Diffuse Kalman filter (Durbin and Koopman)
if kalman_algo ~= 4
% Use standard kalman filter except if the univariate filter is explicitely choosen.
if kalman_algo == 0
kalman_algo = 3;
elseif ~((kalman_algo == 3) || (kalman_algo == 4))
error(['diffuse filter: options_.kalman_algo can only be equal ' ...
'to 0 (default), 3 or 4'])
end
[Z,T,R,QT,Pstar,Pinf] = schur_statespace_transformation(Z,T,R,Q,DynareOptions.qz_criterium);
Zflag = 1;
% Run diffuse kalman filter on first periods.
@ -282,9 +291,9 @@ switch DynareOptions.lik_init
% Multivariate Diffuse Kalman Filter
if no_missing_data_flag
[dLIK,dlik,a,Pstar] = kalman_filter_d(Y, 1, size(Y,2), ...
zeros(mm,1), Pinf, Pstar, ...
kalman_tol, riccati_tol, DynareOptions.presample, ...
T,R,Q,H,Z,mm,pp,rr);
zeros(mm,1), Pinf, Pstar, ...
kalman_tol, riccati_tol, DynareOptions.presample, ...
T,R,Q,H,Z,mm,pp,rr);
else
[dLIK,dlik,a,Pstar] = missing_observations_kalman_filter_d(DynareDataset.missing.aindex,DynareDataset.missing.number_of_observations,DynareDataset.missing.no_more_missing_observations, ...
Y, 1, size(Y,2), ...
@ -295,39 +304,39 @@ switch DynareOptions.lik_init
diffuse_periods = length(dlik);
if isinf(dLIK)
% Go to univariate diffuse filter if singularity problem.
kalman_algo = 4;
singularity_flag = 1;
singular_diffuse_filter
end
end
if (kalman_algo==4)
if singular_diffuse_filter || (kalman_algo==4)
% Univariate Diffuse Kalman Filter
if singularity_flag
if isequal(H,0)
H = zeros(nobs,1);
mmm = mm;
if isequal(H,0)
H1 = zeros(nobs,1);
mmm = mm;
else
if all(all(abs(H-diag(diag(H)))<1e-14))% ie, the covariance matrix is diagonal...
H1 = diag(H);
mmm = mm;
else
if all(all(abs(H-diag(diag(H)))<1e-14))% ie, the covariance matrix is diagonal...
H = diag(H);
mmm = mm;
else
Z = [Z, eye(pp)];
T = blkdiag(T,zeros(pp));
Q = blkdiag(Q,H);
R = blkdiag(R,eye(pp));
Pstar = blkdiag(Pstar,H);
Pinf = blckdiag(Pinf,zeros(pp));
H = zeros(nobs,1);
mmm = mm+pp;
end
Z = [Z, eye(pp)];
T = blkdiag(T,zeros(pp));
Q = blkdiag(Q,H);
R = blkdiag(R,eye(pp));
Pstar = blkdiag(Pstar,H);
Pinf = blckdiag(Pinf,zeros(pp));
H1 = zeros(nobs,1);
mmm = mm+pp;
end
% no need to test again for correlation elements
singularity_flag = 0;
end
[dLIK,dlik,a,Pstar] = univariate_kalman_filter_d(DynareDataset.missing.aindex,DynareDataset.missing.number_of_observations,DynareDataset.missing.no_more_missing_observations, ...
Y, 1, size(Y,2), ...
zeros(mmm,1), Pinf, Pstar, ...
kalman_tol, riccati_tol, DynareOptions.presample, ...
T,R,Q,H,Z,mmm,pp,rr);
% no need to test again for correlation elements
correlated_errors_have_been_checked = 1;
[dLIK,dlik,a,Pstar] = univariate_kalman_filter_d(DynareDataset.missing.aindex,...
DynareDataset.missing.number_of_observations,...
DynareDataset.missing.no_more_missing_observations, ...
Y, 1, size(Y,2), ...
zeros(mmm,1), Pinf, Pstar, ...
kalman_tol, riccati_tol, DynareOptions.presample, ...
T,R,Q,H1,Z,mmm,pp,rr);
diffuse_periods = length(dlik);
end
case 4% Start from the solution of the Riccati equation.
@ -382,10 +391,10 @@ if ((kalman_algo==1) || (kalman_algo==3))% Multivariate Kalman Filter
end
end
if ( singularity_flag || (kalman_algo==2) || (kalman_algo==4) )
if (kalman_algo==2) || (kalman_algo==4)
% Univariate Kalman Filter
% resetting measurement error covariance matrix when necessary %
if singularity_flag
if ~correlated_errors_have_been_checked
if isequal(H,0)
H = zeros(nobs,1);
mmm = mm;

36
matlab/dyn_figure.m Normal file
View File

@ -0,0 +1,36 @@
function h=dyn_figure(DynareOptions,varargin)
%function h=dyn_figure(DynareOptions,varargin)
% initializes figures for DYNARE
%
% INPUTS
% DynareOptions: dynare options
% varargin: the same list of possible inputs of the MATLAB function figure
%
% OUTPUTS
% h : figure handle
%
% SPECIAL REQUIREMENTS
% none
% 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/>.
if DynareOptions.nodisplay,
h = figure(varargin{:},'visible','off');
else
h = figure(varargin{:});
end

45
matlab/dyn_ramsey_static.m
View File

@ -56,12 +56,12 @@ if options_.steadystate_flag
[xx,params,check] = evaluate_steady_state_file(ys,exo_ss,params,...
M.fname,options_.steadystate_flag);
else
n_var = M.orig_endo_nbr+min(find([M.aux_vars.type] == 6)) - 1;
n_var = M.orig_endo_nbr;
xx = oo.steady_state(1:n_var);
[xx,info1] = dynare_solve(nl_func,xx,0);
steady_state = nl_func(xx);
end
steady_state = nl_func(xx);
[junk,junk,steady_state] = nl_func(xx);
@ -71,13 +71,17 @@ rJ = [];
mult = [];
% recovering usefull fields
params = M.params;
endo_nbr = M.endo_nbr;
endo_names = M.endo_names;
exo_nbr = M.exo_nbr;
orig_endo_nbr = M.orig_endo_nbr;
aux_vars_type = [M.aux_vars.type];
orig_endo_aux_nbr = orig_endo_nbr + min(find(aux_vars_type == 6)) - 1;
orig_eq_nbr = M.orig_eq_nbr;
inst_nbr = orig_endo_nbr - orig_eq_nbr;
inst_nbr = orig_endo_aux_nbr - orig_eq_nbr;
% indices of Lagrange multipliers
i_mult = [orig_endo_nbr+(1:orig_eq_nbr)]';
i_mult = [orig_endo_aux_nbr+(1:orig_eq_nbr)]';
fname = M.fname;
max_lead = M.maximum_lead;
max_lag = M.maximum_lag;
@ -86,24 +90,26 @@ max_lag = M.maximum_lag;
i_endo = [1:endo_nbr]';
% indices of endogenous variable except instruments
% i_inst = M.instruments;
% lead_lag incidence matrix for endogenous variables
% lead_lag incidence matrix
i_lag = M.lead_lag_incidence;
if options_.steadystate_flag
k_inst = [];
instruments = options_.instruments;
for i = 1:size(instruments,1)
k_inst = [k_inst; strmatch(options_.instruments(i,:), ...
M.endo_names,'exact')];
k_inst = [k_inst; strmatch(instruments(i,:), ...
endo_names,'exact')];
end
oo.steady_state(k_inst) = x;
[x,params,check] = evaluate_steady_state_file(oo.steady_state,...
[oo.exo_steady_state; ...
oo.exo_det_steady_state] ...
,M.params,M.fname,...
,params,fname,...
options_.steadystate_flag);
end
xx = zeros(endo_nbr,1);
xx(1:length(x)) = x;
% setting steady state of auxiliary variables
% that depends on original endogenous variables
if any([M.aux_vars.type] ~= 6)
@ -112,44 +118,43 @@ if any([M.aux_vars.type] ~= 6)
s_a_v_func = @(z) fh(z,...
[oo.exo_steady_state,...
oo.exo_det_steady_state],...
M.params);
x = s_a_v_func(x);
params);
xx = s_a_v_func(xx);
else
needs_set_auxiliary_variables = 0;
end
% value and Jacobian of objective function
ex = zeros(1,M.exo_nbr);
[U,Uy,Uyy] = feval([fname '_objective_static'],x,ex, M.params);
[U,Uy,Uyy] = feval([fname '_objective_static'],x,ex, params);
Uy = Uy';
Uyy = reshape(Uyy,endo_nbr,endo_nbr);
% set multipliers and auxiliary variables that
% depends on multipliers to 0 to compute residuals
xx = [x; zeros(M.endo_nbr - M.orig_eq_nbr,1)];
[res,fJ] = feval([fname '_static'],xx,[oo.exo_simul oo.exo_det_simul], ...
M.params);
% index of multipliers and corresponding equations
% the auxiliary variables before the Lagrange multipliers are treated
% as ordinary endogenous variables
n_var = M.orig_endo_nbr + min(find([M.aux_vars.type] == 6)) - 1;
aux_eq = [1:n_var orig_endo_nbr+orig_eq_nbr+1:size(fJ,1)];
A = fJ(aux_eq,n_var+1:end);
aux_eq = [1:orig_endo_aux_nbr, orig_endo_aux_nbr+orig_eq_nbr+1:size(fJ,1)];
A = fJ(aux_eq,orig_endo_aux_nbr+1:end);
y = res(aux_eq);
mult = -A\y;
aux_vars = -A\y;
resids1 = y+A*mult;
resids1 = y+A*aux_vars;
if inst_nbr == 1
r1 = sqrt(resids1'*resids1);
else
[q,r,e] = qr([A y]');
r1 = r(end,(orig_endo_nbr-inst_nbr+1:end))';
k = size(A,1)+(1-inst_nbr:0);
r1 = r(end,k)';
end
if options_.steadystate_flag
resids = r1;
else
resids = [res; r1];
resids = [res(orig_endo_nbr+(1:orig_endo_nbr-inst_nbr)); r1];
end
rJ = [];
steady_state = [x(1:n_var); mult];
steady_state = [xx(1:orig_endo_aux_nbr); aux_vars];

49
matlab/dyn_saveas.m Normal file
View File

@ -0,0 +1,49 @@
function dyn_saveas(h,fname,DynareOptions)
%function dyn_saveas(h,fname,DynareOptions)
% save figures for DYNARE
%
% INPUTS
% h : figure handle
% fname : name of the saved figure
% DynareOptions: dynare options
%
% OUTPUTS
% none
%
% SPECIAL REQUIREMENTS
% none
% 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/>.
if strmatch('eps',DynareOptions.graph_format)
eval(['print -depsc2 ' fname '.eps']);
end
if ~exist('OCTAVE_VERSION')
if strmatch('pdf',DynareOptions.graph_format)
eval(['print -dpdf ' fname]);
end
if strmatch('fig',DynareOptions.graph_format)
if DynareOptions.nodisplay
set(h, 'Visible','on');
end
saveas(h,[fname '.fig']);
end
end
if DynareOptions.nodisplay
close(h);
end

20
matlab/dynare_config.m
View File

@ -15,7 +15,7 @@ function dynareroot = dynare_config(path_to_dynare,verbose)
% SPECIAL REQUIREMENTS
% none
% Copyright (C) 2001-2011 Dynare Team
% Copyright (C) 2001-2012 Dynare Team
%
% This file is part of Dynare.
%
@ -70,8 +70,7 @@ if ~exist('OCTAVE_VERSION')
addpath([dynareroot '/missing/rows_columns'])
% Replacement for vec() (inexistent under MATLAB)
addpath([dynareroot '/missing/vec'])
[has_statistics_toolbox junk] = license('checkout','statistics_toolbox');
if ~has_statistics_toolbox
if ~user_has_matlab_license('statistics_toolbox')
% Replacements for functions of the stats toolbox
addpath([dynareroot '/missing/stats/'])
end
@ -101,8 +100,7 @@ if exist('OCTAVE_VERSION')
addpath([dynareroot '/missing/nanmean'])
end
else
[has_statistics_toolbox junk] = license('checkout','statistics_toolbox');
if ~has_statistics_toolbox
if ~user_has_matlab_license('statistics_toolbox')
addpath([dynareroot '/missing/nanmean'])
end
end
@ -113,8 +111,13 @@ if exist('OCTAVE_VERSION')
else
% Add win32 specific paths for Dynare Windows package
if strcmp(computer, 'PCWIN')
if matlab_ver_less_than('7.5')
mexpath = [dynareroot '../mex/matlab/win32-7.0-7.4'];
if matlab_ver_less_than('7.1')
mexpath = [dynareroot '../mex/matlab/win32-7.0-7.0.4'];
if exist(mexpath, 'dir')
addpath(mexpath)
end
elseif matlab_ver_less_than('7.5')
mexpath = [dynareroot '../mex/matlab/win32-7.1-7.4'];
if exist(mexpath, 'dir')
addpath(mexpath)
end
@ -190,6 +193,9 @@ mex_status(3,3) = {'Kronecker products'};
mex_status(4,1) = {'sparse_hessian_times_B_kronecker_C'};
mex_status(4,2) = {'kronecker'};
mex_status(4,3) = {'Sparse kronecker products'};
mex_status(5,1) = {'local_state_space_iteration_2'};
mex_status(5,2) = {'particle/local_state_space_iteration'};
mex_status(5,3) = {'Local state space iteraton (second order)'};
number_of_mex_files = size(mex_status,1);
%% Remove some directories from matlab's path. This is necessary if the user has
%% added dynare_v4/matlab with the subfolders. Matlab has to ignore these

34
matlab/dynare_estimation_1.m
View File

@ -12,7 +12,7 @@ function dynare_estimation_1(var_list_,dname)
% SPECIAL REQUIREMENTS
% none
% Copyright (C) 2003-2011 Dynare Team
% Copyright (C) 2003-2012 Dynare Team
%
% This file is part of Dynare.
%
@ -31,8 +31,29 @@ function dynare_estimation_1(var_list_,dname)
global M_ options_ oo_ estim_params_ bayestopt_ dataset_
% Set particle filter flag.
if options_.order > 1
if options_.particle.status && options_.order==2
disp(' ')
disp('Estimation using a non linear filter!')
disp(' ')
elseif options_.particle.status && options_.order>2
error(['Non linear filter are not implemented with order ' int2str(options_.order) ' approximation of the model!'])
elseif ~options_.particle.status && options_.order==2
disp('If you want to estimate the model with a second order approximation using a non linear filter, set options_.particle.status=1;')
disp('I set order=1!')
options_.order=1;
else
error(['Cannot estimate a model with an order ' int2str(options_.order) ' approximation!'])
end
end
if ~options_.dsge_var
objective_function = str2func('dsge_likelihood');
if options_.particle.status
objective_function = str2func('non_linear_dsge_likelihood');
else
objective_function = str2func('dsge_likelihood');
end
else
objective_function = str2func('DsgeVarLikelihood');
end
@ -267,7 +288,14 @@ if ~isequal(options_.mode_compute,0) && ~options_.mh_posterior_mode_estimation
[xparam1,fval,exitflag] = fminsearch(objective_function,xparam1,optim_options,dataset_,options_,M_,estim_params_,bayestopt_,oo_);
case 8
% Dynare implementation of the simplex algorithm.
[xparam1,fval,exitflag] = simplex_optimization-routine(objective_function,xparam1,optim_options,dataset_,options_,M_,estim_params_,bayestopt_,oo_);
[xparam1,fval,exitflag] = simplex_optimization_routine(objective_function,xparam1,options_.simplex,dataset_,options_,M_,estim_params_,bayestopt_,oo_);
case 9
H0 = 1e-4*ones(nx,1);
warning('off','CMAES:NonfinitenessRange');
warning('off','CMAES:InitialSigma');
[x, fval, COUNTEVAL, STOPFLAG, OUT, BESTEVER] = cmaes(func2str(objective_function),xparam1,H0,options_.cmaes,dataset_,options_,M_,estim_params_,bayestopt_,oo_);
xparam1=BESTEVER.x;
disp(sprintf('\n Objective function at mode: %f',fval))
case 101
myoptions=soptions;
myoptions(2)=1e-6; %accuracy of argument

14
matlab/dynare_estimation_init.m
View File

@ -58,12 +58,8 @@ for i = 1:size(M_.endo_names,1)
end
end
% Set the order of approximation to one (if needed).
if options_.order > 1 && isempty(options_.nonlinear_filter)
disp('This version of Dynare cannot estimate non linearized models!')
disp('Set "order" equal to 1.')
disp(' ')
options_.order = 1;
if options_.order>2
error(['I cannot estimate a model with a ' int2str(options_.order) ' order approximation of the model!'])
end
% Set options_.lik_init equal to 3 if diffuse filter is used or
@ -214,6 +210,7 @@ end
k2 = union(var_obs_index,[dr.nstatic+1:dr.nstatic+dr.npred]', 'rows');
% Set restrict_state to postion of observed + state variables in expanded state vector.
oo_.dr.restrict_var_list = k2;
bayestopt_.restrict_var_list = k2;
% set mf0 to positions of state variables in restricted state vector for likelihood computation.
[junk,bayestopt_.mf0] = ismember([dr.nstatic+1:dr.nstatic+dr.npred]',k2);
% Set mf1 to positions of observed variables in restricted state vector for likelihood computation.
@ -327,8 +324,9 @@ nvx = estim_params_.nvx;
ncx = estim_params_.ncx;
nvn = estim_params_.nvn;
ncn = estim_params_.ncn;
M.params(estim_params_.param_vals(:,1)) = xparam1(nvx+ncx+nvn+ncn+1:end);
if isfield(estim_params_,'param_vals')
M.params(estim_params_.param_vals(:,1)) = xparam1(nvx+ncx+nvn+ncn+1:end);
end;
oo_.steady_state = evaluate_steady_state(oo_.steady_state,M,options_,oo_,steadystate_check_flag);
if all(abs(oo_.steady_state(bayestopt_.mfys))<1e-9)
options_.noconstant = 1;

7
matlab/dynare_solve.m
View File

@ -16,7 +16,7 @@ function [x,info] = dynare_solve(func,x,jacobian_flag,varargin)
% SPECIAL REQUIREMENTS
% none
% Copyright (C) 2001-2011 Dynare Team
% Copyright (C) 2001-2012 Dynare Team
%
% This file is part of Dynare.
%
@ -39,8 +39,7 @@ options_ = set_default_option(options_,'solve_algo',2);
info = 0;
if options_.solve_algo == 0
if ~exist('OCTAVE_VERSION')
[has_optimization_toolbox junk] = license('checkout','optimization_toolbox');
if ~has_optimization_toolbox
if ~user_has_matlab_license('optimization_toolbox')
error('You can''t use solve_algo=0 since you don''t have MATLAB''s Optimization Toolbox')
end
end
@ -106,7 +105,7 @@ elseif options_.solve_algo == 2 || options_.solve_algo == 4
if ~jacobian_flag
fjac = zeros(nn,nn) ;
dh = max(abs(x),options_.gstep*ones(nn,1))*eps^(1/3);
dh = max(abs(x),options_.gstep(1)*ones(nn,1))*eps^(1/3);
for j = 1:nn
xdh = x ;
xdh(j) = xdh(j)+dh(j) ;

500
matlab/ep/extended_path.m
View File

@ -1,17 +1,17 @@
function time_series = extended_path(initial_conditions,sample_size)
% Stochastic simulation of a non linear DSGE model using the Extended Path method (Fair and Taylor 1983). A time
% series of size T is obtained by solving T perfect foresight models.
%
% series of size T is obtained by solving T perfect foresight models.
%
% INPUTS
% o initial_conditions [double] m*nlags array, where m is the number of endogenous variables in the model and
% nlags is the maximum number of lags.
% o sample_size [integer] scalar, size of the sample to be simulated.
%
%
% OUTPUTS
% o time_series [double] m*sample_size array, the simulations.
%
%
% ALGORITHM
%
%
% SPECIAL REQUIREMENTS
% Copyright (C) 2009, 2010, 2011 Dynare Team
@ -31,15 +31,66 @@ function time_series = extended_path(initial_conditions,sample_size)
% 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_ 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)
if ~( options_.bytecode && options_.block )
error('extended_path:: Options bytecode and block are mandatory!')
options_.verbosity = options_.ep.verbosity;
verbosity = options_.ep.verbosity+options_.ep.debug;
% Prepare a structure needed by the matlab implementation of the perfect foresight model solver
pfm.lead_lag_incidence = M_.lead_lag_incidence;
pfm.ny = M_.endo_nbr;
pfm.Sigma_e = M_.Sigma_e;
max_lag = M_.maximum_endo_lag;
pfm.max_lag = max_lag;
if pfm.max_lag > 0
pfm.nyp = nnz(pfm.lead_lag_incidence(1,:));
pfm.iyp = find(pfm.lead_lag_incidence(1,:)>0);
else
pfm.nyp = 0;
pfm.iyp = [];
end
pfm.ny0 = nnz(pfm.lead_lag_incidence(max_lag+1,:));
pfm.iy0 = find(pfm.lead_lag_incidence(max_lag+1,:)>0);
if M_.maximum_endo_lead
pfm.nyf = nnz(pfm.lead_lag_incidence(max_lag+2,:));
pfm.iyf = find(pfm.lead_lag_incidence(max_lag+2,:)>0);
else
pfm.nyf = 0;
pfm.iyf = [];
end
pfm.nd = pfm.nyp+pfm.ny0+pfm.nyf;
pfm.nrc = pfm.nyf+1;
pfm.isp = [1:pfm.nyp];
pfm.is = [pfm.nyp+1:pfm.ny+pfm.nyp];
pfm.isf = pfm.iyf+pfm.nyp;
pfm.isf1 = [pfm.nyp+pfm.ny+1:pfm.nyf+pfm.nyp+pfm.ny+1];
pfm.iz = [1:pfm.ny+pfm.nyp+pfm.nyf];
pfm.periods = options_.ep.periods;
pfm.steady_state = oo_.steady_state;
pfm.params = M_.params;
if M_.maximum_endo_lead
pfm.i_cols_1 = nonzeros(pfm.lead_lag_incidence(max_lag+(1:2),:)');
pfm.i_cols_A1 = find(pfm.lead_lag_incidence(max_lag+(1:2),:)');
else
pfm.i_cols_1 = nonzeros(pfm.lead_lag_incidence(max_lag+1,:)');
pfm.i_cols_A1 = find(pfm.lead_lag_incidence(max_lag+1,:)');
end
if max_lag > 0
pfm.i_cols_T = nonzeros(pfm.lead_lag_incidence(1:2,:)');
else
pfm.i_cols_T = nonzeros(pfm.lead_lag_incidence(1,:)');
end
pfm.i_cols_j = 1:pfm.nd;
pfm.i_upd = pfm.ny+(1:pfm.periods*pfm.ny);
pfm.dynamic_model = str2func([M_.fname,'_dynamic']);
pfm.verbose = options_.ep.verbosity;
pfm.maxit_ = options_.maxit_;
pfm.tolerance = options_.dynatol.f;
exo_nbr = M_.exo_nbr;
periods = options_.periods;
ep = options_.ep;
steady_state = oo_.steady_state;
dynatol = options_.dynatol;
% Set default initial conditions.
if isempty(initial_conditions)
@ -50,28 +101,32 @@ end
options_.maxit_ = options_.ep.maxit;
% Set the number of periods for the perfect foresight model
options_.periods = options_.ep.periods;
periods = options_.ep.periods;
pfm.periods = options_.ep.periods;
pfm.i_upd = pfm.ny+(1:pfm.periods*pfm.ny);
% keep a copy of pfm.i_upd
i_upd = pfm.i_upd;
% Set the algorithm for the perfect foresight solver
options_.stack_solve_algo = options_.ep.stack_solve_algo;
% Set check_stability flag
do_not_check_stability_flag = ~options_.ep.check_stability;
% Compute the first order reduced form if needed.
%
% REMARK. It is assumed that the user did run the same mod file with stoch_simul(order=1) and save
% all the globals in a mat file called linear_reduced_form.mat;
dr = struct();
if options_.ep.init
lrf = load('linear_reduced_form','oo_');
oo_.dr = lrf.oo_.dr; clear('lrf');
if options_.ep.init==2
lambda = .8;
end
options_.order = 1;
[dr,Info,M_,options_,oo_] = resol(1,M_,options_,oo_);
end
% Do not use a minimal number of perdiods for the perfect foresight solver (with bytecode and blocks)
options_.minimal_solving_period = options_.ep.periods;
% Get indices of variables with non zero steady state
idx = find(abs(oo_.steady_state)>0);
options_.minimal_solving_period = 100;%options_.ep.periods;
% Initialize the exogenous variables.
make_ex_;
@ -83,247 +138,250 @@ make_y_;
time_series = zeros(M_.endo_nbr,sample_size);
% Set the covariance matrix of the structural innovations.
variances = diag(M_.Sigma_e);
positive_var_indx = find(variances>0);
covariance_matrix = M_.Sigma_e(positive_var_indx,positive_var_indx);
number_of_structural_innovations = length(covariance_matrix);
covariance_matrix_upper_cholesky = chol(covariance_matrix);
variances = diag(M_.Sigma_e);
positive_var_indx = find(variances>0);
effective_number_of_shocks = length(positive_var_indx);
stdd = sqrt(variances(positive_var_indx));
covariance_matrix = M_.Sigma_e(positive_var_indx,positive_var_indx);
covariance_matrix_upper_cholesky = chol(covariance_matrix);
% (re)Set exo_nbr
%exo_nbr = effective_number_of_shocks;
% Set seed.
if options_.ep.set_dynare_seed_to_default
set_dynare_seed('default');
end
% Set bytecode flag
bytecode_flag = options_.ep.use_bytecode;
% Simulate shocks.
switch options_.ep.innovation_distribution
case 'gaussian'
oo_.ep.shocks = randn(sample_size,number_of_structural_innovations)*covariance_matrix_upper_cholesky;
oo_.ep.shocks = randn(sample_size,effective_number_of_shocks)*covariance_matrix_upper_cholesky;
otherwise
error(['extended_path:: ' options_.ep.innovation_distribution ' distribution for the structural innovations is not (yet) implemented!'])
end
% Set future shocks (Stochastic Extended Path approach)
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{:});
else
rrr = r;
www = w;
end
www = prod(www,2);
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('extended_path:: Unknown stochastic_method option!')
end
else
rrr = zeros(1,number_of_structural_innovations);
www = 1;
nnn = 1;
end
% Initializes some variables.
t = 0;
% Set waitbar (graphic or text mode)
graphic_waitbar_flag = ~( options_.console_mode || exist('OCTAVE_VERSION') );
if graphic_waitbar_flag
hh = waitbar(0,['Please wait. Extended Path simulations...']);
set(hh,'Name','EP simulations');
else
for i=1:2, disp(' '), end
if ~exist('OCTAVE_VERSION')
back = [];
end
end
hh = dyn_waitbar(0,'Please wait. Extended Path simulations...');
set(hh,'Name','EP simulations.');
% Main loop.
while (t<sample_size)
if ~mod(t,10)
if graphic_waitbar_flag
waitbar(t/sample_size);
else
if exist('OCTAVE_VERSION')
printf('Please wait. Extended Path simulations... %3.f%%\r done', 100*t/sample_size);
else
str = sprintf('Please wait. Extended Path simulations... %3.f%% done.', 100*t/sample_size);
fprintf([back '%s'],str);
back=repmat('\b',1,length(str));
end
end
dyn_waitbar(t/sample_size,hh,'Please wait. Extended Path simulations...');
end
% Set period index.
t = t+1;
shocks = oo_.ep.shocks(t,:);
% Put it in oo_.exo_simul (second line).
oo_.exo_simul(2,positive_var_indx) = shocks;
for s = 1:nnn
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;
periods1 = periods;
exo_simul_1 = zeros(periods1+2,exo_nbr);
exo_simul_1(2,:) = oo_.exo_simul(2,:);
pfm1 = pfm;
info_convergence = 0;
if ep.init% Compute first order solution (Perturbation)...
ex = zeros(size(endo_simul_1,2),size(exo_simul_1,2));
ex(1:size(exo_simul_1,1),:) = exo_simul_1;
exo_simul_1 = ex;
initial_path = simult_(initial_conditions,dr,exo_simul_1(2:end,:),1);
endo_simul_1(:,1:end-1) = initial_path(:,1:end-1)*ep.init+endo_simul_1(:,1:end-1)*(1-ep.init);
else
if t==1
endo_simul_1 = repmat(steady_state,1,periods1+2);
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.
elseif options_.ep.init==2
oo_.endo_simul(:,1:end-1) = initial_path(:,1:end-1)*lambda+oo_.endo_simul(:,1:end-1)*(1-lambda);
end
end
% Solve a perfect foresight model (using bytecoded version).
increase_periods = 0;
endo_simul = oo_.endo_simul;
while 1
if ~increase_periods
t0 = tic;
[flag,tmp] = bytecode('dynamic');
ctime = toc(t0);
info.convergence = ~flag;
info.time = ctime;
end
if verbosity
if info.convergence
if t<10
disp(['Time: ' int2str(t) '. Convergence of the perfect foresight model solver!'])
elseif t<100
disp(['Time: ' int2str(t) '. Convergence of the perfect foresight model solver!'])
elseif t<1000
disp(['Time: ' int2str(t) '. Convergence of the perfect foresight model solver!'])
else
disp(['Time: ' int2str(t) '. Convergence of the perfect foresight model solver!'])
end
else
if t<10
disp(['Time: ' int2str(t) '. No convergence of the perfect foresight model solver!'])
elseif t<100
disp(['Time: ' int2str(t) '. No convergence of the perfect foresight model solver!'])
elseif t<1000
disp(['Time: ' int2str(t) '. No convergence of the perfect foresight model solver!'])
else
disp(['Time: ' int2str(t) '. No convergence of the perfect foresight model solver!'])
end
end
end
% Test if periods is big enough.
if ~increase_periods && max(max(abs(tmp(idx,end-options_.ep.lp:end)./tmp(idx,end-options_.ep.lp-1:end-1)-1)))<options_.dynatol.x
break
else
options_.periods = options_.periods + options_.ep.step;
options_.minimal_solving_period = options_.periods;
increase_periods = increase_periods + 1;
if verbosity
if t<10
disp(['Time: ' int2str(t) '. I increase the number of periods to ' int2str(options_.periods) '.'])
elseif t<100
disp(['Time: ' int2str(t) '. I increase the number of periods to ' int2str(options_.periods) '.'])
elseif t<1000
disp(['Time: ' int2str(t) '. I increase the number of periods to ' int2str(options_.periods) '.'])
else
disp(['Time: ' int2str(t) '. I increase the number of periods to ' int2str(options_.periods) '.'])
end
end
if info.convergence
oo_.endo_simul = [ tmp , repmat(oo_.steady_state,1,options_.ep.step) ];
oo_.exo_simul = [ oo_.exo_simul ; zeros(options_.ep.step,size(shocks,2)) ];
tmp_old = tmp;
else
oo_.endo_simul = [ oo_.endo_simul , repmat(oo_.steady_state,1,options_.ep.step) ];
oo_.exo_simul = [ oo_.exo_simul ; zeros(options_.ep.step,size(shocks,2)) ];
end
t0 = tic;
end
% Solve a perfect foresight model.
increase_periods = 0;
% Keep a copy of endo_simul_1
endo_simul = endo_simul_1;
while 1
if ~increase_periods
if bytecode_flag
oo_.endo_simul = endo_simul_1;
oo_.exo_simul = exo_simul_1;
[flag,tmp] = bytecode('dynamic');
ctime = toc(t0);
info.time = info.time+ctime;
if info.convergence
maxdiff = max(max(abs(tmp(:,2:options_.ep.fp)-tmp_old(:,2:options_.ep.fp))));
if maxdiff<options_.dynatol.x
options_.periods = options_.ep.periods;
options_.minimal_solving_period = options_.periods;
oo_.exo_simul = oo_.exo_simul(1:(options_.periods+2),:);
break
end
else
flag = 1;
end
if flag
if options_.ep.stochastic.order == 0
[flag,tmp,err] = solve_perfect_foresight_model(endo_simul_1,exo_simul_1,pfm1);
else
info.convergence = ~flag;
if info.convergence
continue
else
if increase_periods==10;
if verbosity
if t<10
disp(['Time: ' int2str(t) '. Even with ' int2str(options_.periods) ', I am not able to solve the perfect foresight model. Use homotopy instead...'])
elseif t<100
disp(['Time: ' int2str(t) '. Even with ' int2str(options_.periods) ', I am not able to solve the perfect foresight model. Use homotopy instead...'])
elseif t<1000
disp(['Time: ' int2str(t) '. Even with ' int2str(options_.periods) ', I am not able to solve the perfect foresight model. Use homotopy instead...'])
else
disp(['Time: ' int2str(t) '. Even with ' int2str(options_.periods) ', I am not able to solve the perfect foresight model. Use homotopy instead...'])
end
end
break
[flag,tmp] = solve_stochastic_perfect_foresight_model(endo_simul_1,exo_simul_1,pfm1,options_.ep.stochastic.nodes,options_.ep.stochastic.order);
end
end
info_convergence = ~flag;
end
if verbosity
if info_convergence
if t<10
disp(['Time: ' int2str(t) '. Convergence of the perfect foresight model solver!'])
elseif t<100
disp(['Time: ' int2str(t) '. Convergence of the perfect foresight model solver!'])
elseif t<1000
disp(['Time: ' int2str(t) '. Convergence of the perfect foresight model solver!'])
else
disp(['Time: ' int2str(t) '. Convergence of the perfect foresight model solver!'])
end
else
if t<10
disp(['Time: ' int2str(t) '. No convergence of the perfect foresight model solver!'])
elseif t<100
disp(['Time: ' int2str(t) '. No convergence of the perfect foresight model solver!'])
elseif t<1000
disp(['Time: ' int2str(t) '. No convergence of the perfect foresight model solver!'])
else
disp(['Time: ' int2str(t) '. No convergence of the perfect foresight model solver!'])
end
end
end
if do_not_check_stability_flag
% Exit from the while loop.
endo_simul_1 = tmp;
break
else
% Test if periods is big enough.
% Increase the number of periods.
periods1 = periods1 + ep.step;
pfm1.periods = periods1;
pfm1.i_upd = pfm1.ny+(1:pfm1.periods*pfm1.ny);
% Increment the counter.
increase_periods = increase_periods + 1;
if verbosity
if t<10
disp(['Time: ' int2str(t) '. I increase the number of periods to ' int2str(periods1) '.'])
elseif t<100
disp(['Time: ' int2str(t) '. I increase the number of periods to ' int2str(periods1) '.'])
elseif t<1000
disp(['Time: ' int2str(t) '. I increase the number of periods to ' int2str(periods1) '.'])
else
disp(['Time: ' int2str(t) '. I increase the number of periods to ' int2str(periods1) '.'])
end
end
if info_convergence
% If the previous call to the perfect foresight model solver exited
% announcing that the routine converged, adapt the size of endo_simul_1
% and exo_simul_1.
endo_simul_1 = [ tmp , repmat(steady_state,1,ep.step) ];
exo_simul_1 = [ exo_simul_1 ; zeros(ep.step,exo_nbr)];
tmp_old = tmp;
else
% If the previous call to the perfect foresight model solver exited
% announcing that the routine did not converge, then tmp=1... Maybe
% should change that, because in some circonstances it may usefull
% to know where the routine did stop, even if convergence was not
% achieved.
endo_simul_1 = [ endo_simul_1 , repmat(steady_state,1,ep.step) ];
exo_simul_1 = [ exo_simul_1 ; zeros(ep.step,exo_nbr)];
end
% Solve the perfect foresight model with an increased number of periods.
if bytecode_flag
oo_.endo_simul = endo_simul_1;
oo_.exo_simul = exo_simul_1;
[flag,tmp] = bytecode('dynamic');
else
flag = 1;
end
if flag
if options_.ep.stochastic.order == 0
[flag,tmp,err] = solve_perfect_foresight_model(endo_simul_1,exo_simul_1,pfm1);
else
[flag,tmp] = solve_stochastic_perfect_foresight_model(endo_simul_1,exo_simul_1,pfm1,options_.ep.stochastic.nodes,options_.ep.stochastic.order);
end
end
info_convergence = ~flag;
if info_convergence
% If the solver achieved convergence, check that simulated paths did not
% change during the first periods.
% Compute the maximum deviation between old path and new path over the
% first periods
delta = max(max(abs(tmp(:,2)-tmp_old(:,2))));
if delta < dynatol.x
% If the maximum deviation is close enough to zero, reset the number
% of periods to ep.periods
periods1 = ep.periods;
pfm1.periods = periods1;
pfm1.i_upd = pfm1.ny+(1:pfm1.periods*pfm1.ny);
% Cut exo_simul_1 and endo_simul_1 consistently with the resetted
% number of periods and exit from the while loop.
exo_simul_1 = exo_simul_1(1:(periods1+2),:);
endo_simul_1 = endo_simul_1(:,1:(periods1+2));
break
end
else
% The solver did not converge... Try to solve the model again with a bigger
% number of periods, except if the number of periods has been increased more
% than 10 times.
if increase_periods==10;
if verbosity
if t<10
disp(['Time: ' int2str(t) '. Even with ' int2str(periods1) ', I am not able to solve the perfect foresight model. Use homotopy instead...'])
elseif t<100
disp(['Time: ' int2str(t) '. Even with ' int2str(periods1) ', I am not able to solve the perfect foresight model. Use homotopy instead...'])
elseif t<1000
disp(['Time: ' int2str(t) '. Even with ' int2str(periods1) ', I am not able to solve the perfect foresight model. Use homotopy instead...'])
else
disp(['Time: ' int2str(t) '. Even with ' int2str(periods1) ', I am not able to solve the perfect foresight model. Use homotopy instead...'])
end
end
% Exit from the while loop.
break
end
end
end% if info_convergence
end
if ~info.convergence% If the previous step was unsuccesfull, use an homotopic approach
[INFO,tmp] = homotopic_steps(.5,.01,t);
% Cumulate time.
info.time = ctime+INFO.time;
if (~isstruct(INFO) && isnan(INFO)) || ~INFO.convergence
end% while
if ~info_convergence% If exited from the while loop without achieving convergence, use an homotopic approach
if ~do_not_check_stability_flag
periods1 = ep.periods;
pfm1.periods = periods1;
pfm1.i_upd = i_upd;
exo_simul_1 = exo_simul_1(1:(periods1+2),:);
endo_simul_1 = endo_simul_1(:,1:(periods1+2));
end
[INFO,tmp] = homotopic_steps(endo_simul,exo_simul_1,.5,.01,pfm1);
if isstruct(INFO)
info_convergence = INFO.convergence;
else
info_convergence = 0;
end
if ~info_convergence
[INFO,tmp] = homotopic_steps(endo_simul,exo_simul_1,0,.01,pfm1);
if isstruct(INFO)
info_convergence = INFO.convergence;
else
info_convergence = 0;
end
if ~info_convergence
disp('Homotopy:: No convergence of the perfect foresight model solver!')
error('I am not able to simulate this model!');
else
info.convergence = 1;
oo_.endo_simul = tmp;
if verbosity && info.convergence
endo_simul_1 = tmp;
if verbosity && info_convergence
disp('Homotopy:: Convergence of the perfect foresight model solver!')
end
end
else
oo_.endo_simul = tmp;
info_convergence = 1;
endo_simul_1 = tmp;
if verbosity && info_convergence
disp('Homotopy:: Convergence of the perfect foresight model solver!')
end
end
% Save results of the perfect foresight model solver.
time_series(:,t) = time_series(:,t)+ www(s)*oo_.endo_simul(:,2);
%save('simulated_paths.mat','time_series');
% Set initial condition for the nex round.
%initial_conditions = oo_.endo_simul(:,2);
end
%oo_.endo_simul = oo_.endo_simul(:,1:options_.periods+M_.maximum_endo_lag+M_.maximum_endo_lead);
oo_.endo_simul(:,1:end-1) = oo_.endo_simul(:,2:end);
oo_.endo_simul(:,1) = time_series(:,t);
oo_.endo_simul(:,end) = oo_.steady_state;
end
% Save results of the perfect foresight model solver.
time_series(:,t) = endo_simul_1(:,2);
endo_simul_1(:,1:end-1) = endo_simul_1(:,2:end);
endo_simul_1(:,1) = time_series(:,t);
endo_simul_1(:,end) = oo_.steady_state;
end% (while) loop over t
if graphic_waitbar_flag
close(hh);
else
if ~exist('OCTAVE_VERSION')
fprintf(back);
end
end
dyn_waitbar_close(hh);
oo_.endo_simul = oo_.steady_state;
oo_.endo_simul = oo_.steady_state;

405
matlab/ep/extended_path_parfor.m Normal file
View File

@ -0,0 +1,405 @@
function time_series = extended_path(initial_conditions,sample_size)
% Stochastic simulation of a non linear DSGE model using the Extended Path method (Fair and Taylor 1983). A time
% series of size T is obtained by solving T perfect foresight models.
%
% INPUTS
% o initial_conditions [double] m*nlags array, where m is the number of endogenous variables in the model and
% nlags is the maximum number of lags.
% o sample_size [integer] scalar, size of the sample to be simulated.
%
% OUTPUTS
% o time_series [double] m*sample_size array, the simulations.
%
% ALGORITHM
%
% SPECIAL REQUIREMENTS
% Copyright (C) 2009, 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/>.
global M_ options_ oo_
options_.verbosity = options_.ep.verbosity;
verbosity = options_.ep.verbosity+options_.ep.debug;
% Prepare a structure needed by the matlab implementation of the perfect foresight model solver
pfm.lead_lag_incidence = M_.lead_lag_incidence;
pfm.ny = M_.endo_nbr;
pfm.max_lag = M_.maximum_endo_lag;
pfm.nyp = nnz(pfm.lead_lag_incidence(1,:));
pfm.iyp = find(pfm.lead_lag_incidence(1,:)>0);
pfm.ny0 = nnz(pfm.lead_lag_incidence(2,:));
pfm.iy0 = find(pfm.lead_lag_incidence(2,:)>0);
pfm.nyf = nnz(pfm.lead_lag_incidence(3,:));
pfm.iyf = find(pfm.lead_lag_incidence(3,:)>0);
pfm.nd = pfm.nyp+pfm.ny0+pfm.nyf;
pfm.nrc = pfm.nyf+1;
pfm.isp = [1:pfm.nyp];
pfm.is = [pfm.nyp+1:pfm.ny+pfm.nyp];
pfm.isf = pfm.iyf+pfm.nyp;
pfm.isf1 = [pfm.nyp+pfm.ny+1:pfm.nyf+pfm.nyp+pfm.ny+1];
pfm.iz = [1:pfm.ny+pfm.nyp+pfm.nyf];
pfm.periods = options_.ep.periods;
pfm.steady_state = oo_.steady_state;
pfm.params = M_.params;
pfm.i_cols_1 = nonzeros(pfm.lead_lag_incidence(2:3,:)');
pfm.i_cols_A1 = find(pfm.lead_lag_incidence(2:3,:)');
pfm.i_cols_T = nonzeros(pfm.lead_lag_incidence(1:2,:)');
pfm.i_cols_j = 1:pfm.nd;
pfm.i_upd = pfm.ny+(1:pfm.periods*pfm.ny);
pfm.dynamic_model = str2func([M_.fname,'_dynamic']);
pfm.verbose = options_.ep.verbosity;
pfm.maxit_ = options_.maxit_;
pfm.tolerance = options_.dynatol.f;
exo_nbr = M_.exo_nbr;
periods = options_.periods;
ep = options_.ep;
steady_state = oo_.steady_state;
dynatol = options_.dynatol;
% Set default initial conditions.
if isempty(initial_conditions)
initial_conditions = oo_.steady_state;
end
% Set maximum number of iterations for the deterministic solver.
options_.maxit_ = options_.ep.maxit;
% Set the number of periods for the perfect foresight model
periods = options_.ep.periods;
pfm.periods = options_.ep.periods;
pfm.i_upd = pfm.ny+(1:pfm.periods*pfm.ny);
% Set the algorithm for the perfect foresight solver
options_.stack_solve_algo = options_.ep.stack_solve_algo;
% Set check_stability flag
do_not_check_stability_flag = ~options_.ep.check_stability;
% Compute the first order reduced form if needed.
%
% REMARK. It is assumed that the user did run the same mod file with stoch_simul(order=1) and save
% all the globals in a mat file called linear_reduced_form.mat;
dr = struct();
if options_.ep.init
options_.order = 1;
[dr,Info,M_,options_,oo_] = resol(1,M_,options_,oo_);
end
% Do not use a minimal number of perdiods for the perfect foresight solver (with bytecode and blocks)
options_.minimal_solving_period = 100;%options_.ep.periods;
% Initialize the exogenous variables.
make_ex_;
% Initialize the endogenous variables.
make_y_;
% Initialize the output array.
time_series = zeros(M_.endo_nbr,sample_size);
% Set the covariance matrix of the structural innovations.
variances = diag(M_.Sigma_e);
positive_var_indx = find(variances>0);
effective_number_of_shocks = length(positive_var_indx);
stdd = sqrt(variances(positive_var_indx));
covariance_matrix = M_.Sigma_e(positive_var_indx,positive_var_indx);
covariance_matrix_upper_cholesky = chol(covariance_matrix);
% Set seed.
if options_.ep.set_dynare_seed_to_default
set_dynare_seed('default');
end
% Set bytecode flag
bytecode_flag = options_.ep.use_bytecode;
% Simulate shocks.
switch options_.ep.innovation_distribution
case 'gaussian'
oo_.ep.shocks = randn(sample_size,effective_number_of_shocks)*covariance_matrix_upper_cholesky;
otherwise
error(['extended_path:: ' options_.ep.innovation_distribution ' distribution for the structural innovations is not (yet) implemented!'])
end
% Set future shocks (Stochastic Extended Path approach)
if options_.ep.stochastic.status
switch options_.ep.stochastic.method
case 'tensor'
switch options_.ep.stochastic.ortpol
case 'hermite'
[r,w] = gauss_hermite_weights_and_nodes(options_.ep.stochastic.nodes);
otherwise
error('extended_path:: Unknown orthogonal polynomial option!')
end
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{:});
else
rrr = r;
www = w;
end
www = prod(www,2);
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('extended_path:: Unknown stochastic_method option!')
end
else
rrr = zeros(1,effective_number_of_shocks);
www = 1;
nnn = 1;
end
% Initializes some variables.
t = 0;
% Set waitbar (graphic or text mode)
hh = dyn_waitbar(0,'Please wait. Extended Path simulations...');
set(hh,'Name','EP simulations.');
if options_.ep.memory
mArray1 = zeros(M_.endo_nbr,100,nnn,sample_size);
mArray2 = zeros(M_.exo_nbr,100,nnn,sample_size);
end
% Main loop.
while (t<sample_size)
if ~mod(t,10)
dyn_waitbar(t/sample_size,hh,'Please wait. Extended Path simulations...');
end
% Set period index.
t = t+1;
shocks = oo_.ep.shocks(t,:);
% Put it in oo_.exo_simul (second line).
oo_.exo_simul(2,positive_var_indx) = shocks;
parfor s = 1:nnn
periods1 = periods;
exo_simul_1 = zeros(periods1+2,exo_nbr);
pfm1 = pfm;
info_convergence = 0;
switch ep.stochastic.ortpol
case 'hermite'
for u=1:ep.stochastic.order
exo_simul_1(2+u,positive_var_indx) = rrr(s,(((u-1)*effective_number_of_shocks)+1):(u*effective_number_of_shocks))*covariance_matrix_upper_cholesky;
end
otherwise
error('extended_path:: Unknown orthogonal polynomial option!')
end
if ep.stochastic.order && ep.stochastic.scramble
exo_simul_1(2+ep.stochastic.order+1:2+ep.stochastic.order+ep.stochastic.scramble,positive_var_indx) = ...
randn(ep.stochastic.scramble,effective_number_of_shocks)*covariance_matrix_upper_cholesky;
end
if ep.init% Compute first order solution (Perturbation)...
ex = zeros(size(endo_simul_1,2),size(exo_simul_1,2));
ex(1:size(exo_simul_1,1),:) = exo_simul_1;
exo_simul_1 = ex;
initial_path = simult_(initial_conditions,dr,exo_simul_1(2:end,:),1);
endo_simul_1(:,1:end-1) = initial_path(:,1:end-1)*ep.init+endo_simul_1(:,1:end-1)*(1-ep.init);
else
endo_simul_1 = repmat(steady_state,1,periods1+2);
end
% Solve a perfect foresight model.
increase_periods = 0;
endo_simul = endo_simul_1;
while 1
if ~increase_periods
if bytecode_flag
[flag,tmp] = bytecode('dynamic');
else
flag = 1;
end
if flag
[flag,tmp] = solve_perfect_foresight_model(endo_simul_1,exo_simul_1,pfm1);
end
info_convergence = ~flag;
end
if verbosity
if info_convergence
if t<10
disp(['Time: ' int2str(t) '. Convergence of the perfect foresight model solver!'])
elseif t<100
disp(['Time: ' int2str(t) '. Convergence of the perfect foresight model solver!'])
elseif t<1000
disp(['Time: ' int2str(t) '. Convergence of the perfect foresight model solver!'])
else
disp(['Time: ' int2str(t) '. Convergence of the perfect foresight model solver!'])
end
else
if t<10
disp(['Time: ' int2str(t) '. No convergence of the perfect foresight model solver!'])
elseif t<100
disp(['Time: ' int2str(t) '. No convergence of the perfect foresight model solver!'])
elseif t<1000
disp(['Time: ' int2str(t) '. No convergence of the perfect foresight model solver!'])
else
disp(['Time: ' int2str(t) '. No convergence of the perfect foresight model solver!'])
end
end
end
if do_not_check_stability_flag
% Exit from the while loop.
endo_simul_1 = tmp;
break
else
% Test if periods is big enough.
% Increase the number of periods.
periods1 = periods1 + ep.step;
pfm1.periods = periods1;
pfm1.i_upd = pfm1.ny+(1:pfm1.periods*pfm1.ny);
% Increment the counter.
increase_periods = increase_periods + 1;
if verbosity
if t<10
disp(['Time: ' int2str(t) '. I increase the number of periods to ' int2str(periods1) '.'])
elseif t<100
disp(['Time: ' int2str(t) '. I increase the number of periods to ' int2str(periods1) '.'])
elseif t<1000
disp(['Time: ' int2str(t) '. I increase the number of periods to ' int2str(periods1) '.'])
else
disp(['Time: ' int2str(t) '. I increase the number of periods to ' int2str(periods1) '.'])
end
end
if info_convergence
% If the previous call to the perfect foresight model solver exited
% announcing that the routine converged, adapt the size of endo_simul_1
% and exo_simul_1.
endo_simul_1 = [ tmp , repmat(steady_state,1,ep.step) ];
exo_simul_1 = [ exo_simul_1 ; zeros(ep.step,size(shocks,2)) ];
tmp_old = tmp;
else
% If the previous call to the perfect foresight model solver exited
% announcing that the routine did not converge, then tmp=1... Maybe
% should change that, because in some circonstances it may usefull
% to know where the routine did stop, even if convergence was not
% achieved.
endo_simul_1 = [ endo_simul_1 , repmat(steady_state,1,ep.step) ];
exo_simul_1 = [ exo_simul_1 ; zeros(ep.step,size(shocks,2)) ];
end
% Solve the perfect foresight model with an increased number of periods.
if bytecode_flag
[flag,tmp] = bytecode('dynamic');
else
flag = 1;
end
if flag
[flag,tmp] = solve_perfect_foresight_model(endo_simul_1,exo_simul_1,pfm1);
end
info_convergence = ~flag;
if info_convergence
% If the solver achieved convergence, check that simulated paths did not
% change during the first periods.
% Compute the maximum deviation between old path and new path over the
% first periods
delta = max(max(abs(tmp(:,2:ep.fp)-tmp_old(:,2:ep.fp))));
if delta < dynatol.x
% If the maximum deviation is close enough to zero, reset the number
% of periods to ep.periods
periods1 = ep.periods;
pfm1.periods = periods1;
pfm1.i_upd = pfm1.ny+(1:pfm1.periods*pfm1.ny);
% Cut exo_simul_1 and endo_simul_1 consistently with the resetted
% number of periods and exit from the while loop.
exo_simul_1 = exo_simul_1(1:(periods1+2),:);
endo_simul_1 = endo_simul_1(:,1:(periods1+2));
break
end
else
% The solver did not converge... Try to solve the model again with a bigger
% number of periods, except if the number of periods has been increased more
% than 10 times.
if increase_periods==10;
if verbosity
if t<10
disp(['Time: ' int2str(t) '. Even with ' int2str(periods1) ', I am not able to solve the perfect foresight model. Use homotopy instead...'])
elseif t<100
disp(['Time: ' int2str(t) '. Even with ' int2str(periods1) ', I am not able to solve the perfect foresight model. Use homotopy instead...'])
elseif t<1000
disp(['Time: ' int2str(t) '. Even with ' int2str(periods1) ', I am not able to solve the perfect foresight model. Use homotopy instead...'])
else
disp(['Time: ' int2str(t) '. Even with ' int2str(periods1) ', I am not able to solve the perfect foresight model. Use homotopy instead...'])
end
end
% Exit from the while loop.
break
end
end% if info_convergence
end
end% while
if ~info_convergence% If exited from the while loop without achieving convergence, use an homotopic approach
[INFO,tmp] = homotopic_steps(.5,.01,pfm1);
if (~isstruct(INFO) && isnan(INFO)) || ~info_convergence
[INFO,tmp] = homotopic_steps(0,.01,pfm1);
if ~info_convergence
disp('Homotopy:: No convergence of the perfect foresight model solver!')
error('I am not able to simulate this model!');
else
info_convergence = 1;
endo_simul_1 = tmp;
if verbosity && info_convergence
disp('Homotopy:: Convergence of the perfect foresight model solver!')
end
end
else
info_convergence = 1;
endo_simul_1 = tmp;
if verbosity && info_convergence
disp('Homotopy:: Convergence of the perfect foresight model solver!')
end
end
end
% Save results of the perfect foresight model solver.
if ep.memory
mArray1(:,:,s,t) = endo_simul_1(:,1:100);
mArrat2(:,:,s,t) = transpose(exo_simul_1(1:100,:));
end
results(:,s) = www(s)*endo_simul_1(:,2);
end
time_series(:,t) = sum(results,2);
oo_.endo_simul(:,1:end-1) = oo_.endo_simul(:,2:end);
oo_.endo_simul(:,1) = time_series(:,t);
oo_.endo_simul(:,end) = oo_.steady_state;
end% (while) loop over t
dyn_waitbar_close(hh);
oo_.endo_simul = oo_.steady_state;
if ep.memory
save([M_.fname '_memory'],'mArray1','mArray2','www');
end

294
matlab/ep/homotopic_steps.m
View File

@ -1,20 +1,21 @@
function [info,tmp] = homotopic_steps(initial_weight,step_length,time)
global oo_ options_ M_
function [info,tmp] = homotopic_steps(endo_simul0,exo_simul0,initial_weight,step_length,pfm)
global options_ oo_
%Set bytecode flag
bytecode_flag = options_.ep.use_bytecode;
% Set increase and decrease factors.
increase_factor = 5.0;
decrease_factor = 0.2;
% Save current state of oo_.endo_simul and oo_.exo_simul.
endo_simul = oo_.endo_simul;
exxo_simul = oo_.exo_simul;
% Reset exo_simul to zero.
oo_.exo_simul = zeros(size(oo_.exo_simul));
endo_simul = endo_simul0;
exxo_simul = exo_simul0;
initial_step_length = step_length;
max_iter = 1000/step_length;
weight = initial_weight;
verbose = 1;
iter = 0;
ctime = 0;
verbose = options_.ep.debug;
reduce_step_flag = 0;
@ -22,183 +23,114 @@ if verbose
format long
end
homotopy_1 = 1; % Only innovations are rescaled. Starting from weight equal to initial_weight.
homotopy_2 = 0; % Only innovations are rescaled. Starting from weight equal to zero.
disp(' ')
if homotopy_1
while weight<1
iter = iter+1;
oo_.exo_simul(2,:) = weight*exxo_simul(2,:);
t0 = tic;
% (re)Set iter.
iter = 0;
% (re)Set iter.
jter = 0;
% (re)Set weight.
weight = initial_weight;
% (re)Set exo_simul to zero.
exo_simul0 = zeros(size(exo_simul0));
while weight<1
iter = iter+1;
exo_simul0(2,:) = weight*exxo_simul(2,:);
if bytecode_flag
oo_.endo_simul = endo_simul_1;
oo_.exo_simul = exo_simul_1;
[flag,tmp] = bytecode('dynamic');
TeaTime = toc(t0);
ctime = ctime+TeaTime;
%old_weight = weight;
info.convergence = ~flag;
if verbose
if ~info.convergence
disp(['Iteration n° ' int2str(iter) ', weight is ' num2str(weight,8) ', Convergence problem!' ])
else
disp(['Iteration n° ' int2str(iter) ', weight is ' num2str(weight,8) ', Ok!' ])
end
end
if ~info.convergence
if abs(weight-initial_weight)<1e-12% First iterations.
if verbose
disp('I am reducing the initial weight!')
end
initial_weight = initial_weight/1.1;
weight = initial_weight;
if weight<1/4
homotopy_1 = 0;
homotopy_2 = 1;
break
end
continue
else% A good initial weight has been obtained. In case of convergence problem we have to reduce the step length.
if verbose
disp('I am reducing the step length!')
end
weight = weight-step_length;
step_length=step_length/10;
weight = weight+step_length;
if 10*step_length<options_.dynatol.x
homotopy_1 = 0;
homotopy_2 = 0;
break
end
continue
end
else
flag = 1;
end
if flag
[flag,tmp] = solve_perfect_foresight_model(endo_simul0,exo_simul0,pfm);
end
info.convergence = ~flag;% Equal to one if the perfect foresight solver converged for the current value of weight.
if verbose
if info.convergence
disp(['Iteration n° ' int2str(iter) ', weight is ' num2str(weight,8) ', Ok!' ])
else
oo_.endo_simul = tmp;
info.time = ctime;
if abs(1-weight)<=1e-12;
homotopy_1 = 0;
homotopy_2 = 0;
break
end
weight = weight+step_length;
%step_length = initial_step_length;
end
if iter>max_iter
info = NaN;
return
disp(['Iteration n° ' int2str(iter) ', weight is ' num2str(weight,8) ', Convergence problem!' ])
end
end
if weight<1 && homotopy_1
oo_.exo_simul(2,:) = exxo_simul(2,:);
t0 = tic;
if info.convergence
%if d<stochastic_extended_path_depth
endo_simul0 = tmp;
%end
jter = jter + 1;
if jter>3
if verbose
disp('I am increasing the step length!')
end
step_length=step_length*increase_factor;
jter = 0;
end
if abs(1-weight)<options_.dynatol.x;
break
end
weight = weight+step_length;
else% Perfect foresight solver failed for the current value of weight.
if initial_weight>0 && abs(weight-initial_weight)<1e-12% First iteration, the initial weight is too high.
if verbose
disp('I am reducing the initial weight!')
end
initial_weight = initial_weight/2;
weight = initial_weight;
if weight<1e-12
endo_simul0 = endo_simul;
exo_simul0 = exxo_simul;
info.convergence = 0;
info.depth = d;
tmp = [];
return
end
continue
else% Initial weight is OK, but the perfect foresight solver failed on some subsequent iteration.
if verbose
disp('I am reducing the step length!')
end
jter = 0;
if weight>0
weight = weight-step_length;
end
step_length=step_length*decrease_factor;
weight = weight+step_length;
if step_length<options_.dynatol.x
break
end
continue
end
end
if iter>max_iter
info = NaN;
return
end
end
if weight<1
exo_simul0 = exxo_simul;
if bytecode_flag
oo_.endo_simul = endo_simul_1;
oo_.exo_simul = exo_simul_1;
[flag,tmp] = bytecode('dynamic');
TeaTime = toc(t0);
ctime = ctime+TeaTime;
info.convergence = ~flag;
info.time = ctime;
if info.convergence
oo_.endo_simul = tmp;
homotopy_1 = 0;
homotopy_2 = 0;
else
flag = 1;
end
if flag
[flag,tmp] = solve_perfect_foresight_model(endo_simul0,exo_simul0,pfm);
end
info.convergence = ~flag;
if info.convergence
endo_simul0 = tmp;
return
else
if step_length>options_.dynatol.x
endo_simul0 = endo_simul;
exo_simul0 = exxo_simul;
info.convergence = 0;
info.depth = d;
tmp = [];
return
else
if step_length>1e-12
if verbose
disp('I am reducing step length!')
end
step_length=step_length/2;
else
weight = initial_weight;
step_length = initial_step_length;
info = NaN;
homotopy_2 = 1;
homotopy_1 = 0;
end
error('extended_path::homotopy: Oups! I did my best, but I am not able to simulate this model...')
end
end
end
iter = 0;
weight = 0;
if homotopy_2
while weight<1
iter = iter+1;
oo_.exo_simul(2,:) = weight*exxo_simul(2,:);
if time==1
oo_.endo_simul = repmat(oo_.steady_state,1,size(oo_.endo_simul,2));
else
oo_.endo_simul = endo_simul;
end
t0 = tic;
[flag,tmp] = bytecode('dynamic');
TeaTime = toc(t0);
ctime = ctime+TeaTime;
old_weight = weight;
info.convergence = ~flag;
if verbose
if ~info.convergence
disp(['Iteration n° ' int2str(iter) ', weight is ' num2str(old_weight,8) ', Convergence problem!' ])
else
disp(['Iteration n° ' int2str(iter) ', weight is ' num2str(old_weight,8) ', Ok!' ])
end
end
if ~info.convergence
if iter==1
disp('I am not able to simulate this model!')
disp('There is something wrong with the initial condition of the homotopic')
disp('approach...')
error(' ')
else
if verbose
disp('I am reducing the step length!')
end
step_length=step_length/10;
if 10*step_length<options_.dynatol.x
homotopy_1 = 0;
homotopy_2 = 0;
break
end
weight = old_weight+step_length;
end
else
oo_.endo_simul = tmp;
info.time = ctime;
if abs(1-weight)<=1e-12;
homotopy_2 = 1;
break
end
weight = weight+step_length;
step_length = initial_step_length;
end
if iter>max_iter
info = NaN;
return
end
end
if weight<1 && homotopy_2
oo_.exo_simul(2,:) = exxo_simul(2,:);
t0 = tic;
[flag,tmp] = bytecode('dynamic');
TeaTime = toc(t0);
ctime = ctime+TeaTime;
info.convergence = ~flag;
info.time = ctime;
if info.convergence
oo_.endo_simul = tmp;
homotopy_1 = 0;
else
if step_length>1e-12
if verbose
disp('I am reducing step length!')
end
step_length=step_length/2;
else
weight = initial_weight;
step_length = initial_step_length;
info = NaN;
homotopy_2 = 1;
homotopy_1 = 0;
end
end
end
end

4
matlab/evaluate_static_model.m
View File

@ -53,6 +53,10 @@ function [residuals,check1,jacob] = evaluate_static_model(ys,exo_ss,params,M,opt
% have zero residuals by construction
if ~isempty(mfsb)
residuals(mfsb) = feval(fh_static,b,ys,exo_ss,params);
else
%need to evaluate the recursive blocks to compute the
%temporary terms
feval(fh_static,b,ys,exo_ss,params);
end
end
else

32
matlab/evaluate_steady_state.m
View File

@ -1,7 +1,7 @@
function [ys,params,info] = evaluate_steady_state(ys_init,M,options,oo,steadystate_check_flag)
% function [ys,info] = evaluate_steady_state(M,options,oo)
% Computes the steady state
%
% Computes the steady state
%
% INPUTS
% ys_init vector initial values used to compute the steady
% state
@ -10,8 +10,8 @@ function [ys,params,info] = evaluate_steady_state(ys_init,M,options,oo,steadysta
% oo struct output results
% steadystate_check_flag boolean if true, check that the
% steadystate verifies the
% static model
%
% static model
%
% OUTPUTS
% ys vector steady state
% params vector model parameters possibly
@ -41,23 +41,27 @@ function [ys,params,info] = evaluate_steady_state(ys_init,M,options,oo,steadysta
info = 0;
check = 0;
steadystate_flag = options.steadystate_flag;
params = M.params;
exo_ss = [oo.exo_steady_state; oo.exo_det_steady_state];
updated_params_flag = 0;
if length(M.aux_vars) > 0
h_set_auxiliary_variables = str2func([M.fname '_set_auxiliary_variables']);
h_set_auxiliary_variables = str2func([M.fname '_set_auxiliary_variables']);
ys_init = h_set_auxiliary_variables(ys_init,exo_ss,M.params);
end
if options.ramsey_policy
[ys,params] = dyn_ramsey_static(ys_init,M,options,oo);
elseif steadystate_flag
% explicit steady state file
[ys,params1,check] = evaluate_steady_state_file(ys_init,exo_ss,params,M.fname,steadystate_flag);
updated_params_flag = max(abs(params1-params)) > 1e-12;
if ~length(find(isnan(params))) && ~length(find(isnan(params1)))
updated_params_flag = max(abs(params1-params))>1e-12;
elseif length(find(isnan(params))) && ~length(find(isnan(params1)))
updated_params_flag = 1;
end
if updated_params_flag
params = params1;
end
@ -88,7 +92,7 @@ function [ys,params,info] = evaluate_steady_state(ys_init,M,options,oo,steadysta
end
[ys,check] = dynare_solve('restricted_steadystate',...
ys(indv),...
options.jacobian_flag, ...
options.jacobian_flag, ...
exo_ss,indv);
end
elseif (options.bytecode == 0 && options.block == 0)
@ -96,7 +100,7 @@ function [ys,params,info] = evaluate_steady_state(ys_init,M,options,oo,steadysta
% non linear model
[ys,check] = dynare_solve([M.fname '_static'],...
ys_init,...
options.jacobian_flag, ...
options.jacobian_flag, ...
exo_ss, params);
else
% linear model
@ -140,13 +144,13 @@ function [ys,params,info] = evaluate_steady_state(ys_init,M,options,oo,steadysta
return
end
if options.steadystate_flag && updated_params_flag && ~isreal(M.params)
if options.steadystate_flag && updated_params_flag && ~isreal(params)
info(1) = 23;
info(2) = sum(imag(M.params).^2);
info(2) = sum(imag(params).^2);
return
end
if options.steadystate_flag && updated_params_flag && ~isempty(find(isnan(M.params)))
if options.steadystate_flag && updated_params_flag && ~isempty(find(isnan(params)))
info(1) = 24;
info(2) = NaN;
return

132
matlab/gauss_legendre_weights_and_nodes.m Normal file
View File

@ -0,0 +1,132 @@
function [nodes,weights] = gauss_legendre_weights_and_nodes(n,a,b)
% Computes the weights and nodes for a Legendre Gaussian quadrature rule.
%@info:
%! @deftypefn {Function File} {@var{nodes}, @var{weights} =} gauss_hermite_weights_and_nodes (@var{n})
%! @anchor{gauss_legendre_weights_and_nodes}
%! @sp 1
%! Computes the weights and nodes for a Legendre Gaussian quadrature rule. designed to approximate integrals
%! on the finite interval (-1,1) of an unweighted smooth function.
%! @sp 2
%! @strong{Inputs}
%! @sp 1
%! @table @ @var
%! @item n
%! Positive integer scalar, number of nodes (order of approximation).
%! @item a
%! Double scalar, lower bound.
%! @item b
%! Double scalar, upper bound.
%! @end table
%! @sp 1
%! @strong{Outputs}
%! @sp 1
%! @table @ @var
%! @item nodes
%! n*1 vector of doubles, the nodes (roots of an order n Legendre polynomial)
%! @item weights
%! n*1 vector of doubles, the associated weights.
%! @end table
%! @sp 2
%! @strong{Remarks:}
%! Only the first input argument (the number of nodes) is mandatory. The second and third input arguments
%! are used if a change of variables is necessary (ie if we need nodes over the interval [a,b] instead of
%! of the default interval [-1,1]).
%! @sp 2
%! @strong{This function is called by:}
%! @sp 2
%! @strong{This function calls:}
%! @sp 2
%! @end deftypefn
%@eod:
% 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/>.
% AUTHOR(S) stephane DOT adjemian AT univ DASH lemans DOT fr
bb = sqrt(1./(4-(1./transpose(1:n-1)).^2));
aa = zeros(n,1);
JacobiMatrix = diag(bb,1)+diag(aa)+diag(bb,-1);
[JacobiEigenVectors,JacobiEigenValues] = eig(JacobiMatrix);
[nodes,idx] = sort(diag(JacobiEigenValues));
JacobiEigenVector = JacobiEigenVectors(1,:);
JacobiEigenVector = transpose(JacobiEigenVector(idx));
weights = 2*JacobiEigenVector.^2;
if nargin==3
weights = .5*(b-a)*weights;
nodes = .5*(nodes+1)*(b-a)+a;
end
%@test:1
%$ [n2,w2] = gauss_legendre_weights_and_nodes(2);
%$ [n3,w3] = gauss_legendre_weights_and_nodes(3);
%$ [n4,w4] = gauss_legendre_weights_and_nodes(4);
%$ [n5,w5] = gauss_legendre_weights_and_nodes(5);
%$ [n7,w7] = gauss_legendre_weights_and_nodes(7);
%$
%$
%$ % Expected nodes (taken from Judd (1998, table 7.2)).
%$ e2 = .5773502691; e2 = [-e2; e2];
%$ e3 = .7745966692; e3 = [-e3; 0 ; e3];
%$ e4 = [.8611363115; .3399810435]; e4 = [-e4; flipud(e4)];
%$ e5 = [.9061798459; .5384693101]; e5 = [-e5; 0; flipud(e5)];
%$ e7 = [.9491079123; .7415311855; .4058451513]; e7 = [-e7; 0; flipud(e7)];
%$
%$ % Expected weights (taken from Judd (1998, table 7.2) and http://en.wikipedia.org/wiki/Gaussian_quadrature).
%$ f2 = [1; 1];
%$ f3 = [5; 8; 5]/9;
%$ f4 = [18-sqrt(30); 18+sqrt(30)]; f4 = [f4; flipud(f4)]/36;
%$ f5 = [322-13*sqrt(70); 322+13*sqrt(70)]/900; f5 = [f5; 128/225; flipud(f5)];
%$ f7 = [.1294849661; .2797053914; .3818300505]; f7 = [f7; .4179591836; flipud(f7)];
%$
%$ % Check the results.
%$ t(1) = dyn_assert(e2,n2,1e-9);
%$ t(2) = dyn_assert(e3,n3,1e-9);
%$ t(3) = dyn_assert(e4,n4,1e-9);
%$ t(4) = dyn_assert(e5,n5,1e-9);
%$ t(5) = dyn_assert(e7,n7,1e-9);
%$ t(6) = dyn_assert(w2,f2,1e-9);
%$ t(7) = dyn_assert(w3,f3,1e-9);
%$ t(8) = dyn_assert(w4,f4,1e-9);
%$ t(9) = dyn_assert(w5,f5,1e-9);
%$ t(10) = dyn_assert(w7,f7,1e-9);
%$ T = all(t);
%@eof:1
%@test:2
%$ nmax = 50;
%$
%$ t = zeros(nmax,1);
%$
%$ for i=1:nmax
%$ [n,w] = gauss_legendre_weights_and_nodes(i);
%$ t(i) = dyn_assert(sum(w),2,1e-12);
%$ end
%$
%$ T = all(t);
%@eof:2
%@test:3
%$ [n,w] = gauss_legendre_weights_and_nodes(9,pi,2*pi);
%$ % Check that the
%$ t(1) = all(n>pi);
%$ t(2) = all(n<2*pi);
%$ t(3) = dyn_assert(sum(w),pi,1e-12);
%$ T = all(t);
%@eof:3

73
matlab/gensylv_fp.m Normal file
View File

@ -0,0 +1,73 @@
function X = gensylv_fp(A, B, C, D, block)
% function X = gensylv_fp(A, B, C, D)
% Solve the Sylvester equation:
% A * X + B * X * C + D = 0
% INPUTS
% A
% B
% C
% D
% block : block number (for storage purpose)
% OUTPUTS
% X solution
%
% ALGORITHM
% fixed point method
% MARLLINY MONSALVE (2008): "Block linear method for large scale
% Sylvester equations", Computational & Applied Mathematics, Vol 27, n°1,
% p47-59
% ||A^-1||.||B||.||C|| < 1 is a suffisant condition:
% - to get a unique solution for the Sylvester equation
% - to get a convergent fixed-point algorithm
%
% SPECIAL REQUIREMENTS
% none.
% Copyright (C) 1996-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/>.
%tol = 1e-07;
%tol = 1e-13;
tol = 1e-12;
evol = 100;
A1 = inv(A);
eval(['persistent hxo_' int2str(block) ';']);
hxo = eval(['hxo_' int2str(block) ';']);
if isempty(hxo)
X = zeros(size(B, 2), size(C, 1));
else
X = hxo;
end;
it_fp = 0;
maxit_fp = 1000;
Z = - (B * X * C + D);
while it_fp < maxit_fp && evol > tol;
%X_old = X;
%X = - A1 * ( B * X * C + D);
%evol = max(max(abs(X - X_old)));
X = A1 * Z;
Z_old = Z;
Z = - (B * X * C + D);
evol = max(sum(abs(Z - Z_old))); %norm_1
%evol = max(sum(abs(Z - Z_old)')); %norm_inf
it_fp = it_fp + 1;
end;
%fprintf('sylvester it_fp=%d evol=%g | ',it_fp,evol);
if evol < tol
eval(['hxo_' int2str(block) ' = X;']);
else
error(['convergence not achieved in fixed point solution of Sylvester equation after ' int2str(it_fp) ' iterations']);
end;

108
matlab/global_initialization.m
View File

@ -11,7 +11,7 @@ function global_initialization()
% SPECIAL REQUIREMENTS
% none
% Copyright (C) 2003-2011 Dynare Team
% Copyright (C) 2003-2012 Dynare Team
%
% This file is part of Dynare.
%
@ -41,11 +41,13 @@ options_.terminal_condition = 0;
options_.rplottype = 0;
options_.smpl = 0;
options_.dynatol.f = 1e-5;
options_.dynatol.x = 1e-5;
options_.dynatol.x = 1e-5;
options_.maxit_ = 10;
options_.slowc = 1;
options_.timing = 0;
options_.gstep = 1e-2;
options_.gstep = ones(2,1);
options_.gstep(1) = 1e-2;
options_.gstep(2) = 1.0;
options_.scalv = 1;
options_.debug = 0;
options_.initval_file = 0;
@ -57,9 +59,12 @@ options_.solve_tolf = eps^(1/3);
options_.solve_tolx = eps^(2/3);
options_.solve_maxit = 500;
% Default number of threads for parallelized mex files.
options_.mode_check_neighbourhood_size = 0.5;
% Default number of threads for parallelized mex files.
options_.threads.kronecker.A_times_B_kronecker_C = 1;
options_.threads.kronecker.sparse_hessian_times_B_kronecker_C = 1;
options_.threads.local_state_space_iteration_2 = 1;
% steady state
options_.jacobian_flag = 1;
@ -109,10 +114,16 @@ options_.SpectralDensity = 0;
% Extended path options
%
% Set debug flag
ep.debug = 0;
% Set memory flag
ep.memory = 0;
% Set verbose mode
ep.verbosity = 0;
% Set bytecode flag
ep.use_bytecode = 0;
% Initialization of the perfect foresight equilibrium paths
% * init=0, previous solution is used.
% * 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.
ep.init = 0;
@ -122,10 +133,12 @@ ep.maxit = 500;
ep.periods = 200;
% Default step for increasing the number of periods if needed
ep.step = 50;
% Set check_stability flag
ep.check_stability = 0;
% Define last periods used to test if the solution is stable with respect to an increase in the number of periods.
ep.lp = 5;
% Define first periods used to test if the solution is stable with respect to an increase in the number of periods.
ep.fp = 100;
ep.fp = 2;
% Define the distribution for the structural innovations.
ep.innovation_distribution = 'gaussian';
% Set flag for the seed
@ -133,9 +146,9 @@ ep.set_dynare_seed_to_default = 1;
% Set algorithm for the perfect foresight solver
ep.stack_solve_algo = 4;
% Stochastic extended path related options.
ep.stochastic.status = 0;
ep.stochastic.method = 'tensor';
ep.stochastic.order = 1;
ep.stochastic.ortpol = 'hermite';
ep.stochastic.order = 0;
ep.stochastic.nodes = 5;
ep.stochastic.pruned.status = 0;
ep.stochastic.pruned.relative = 1e-5;
@ -146,37 +159,39 @@ options_.ep = ep;
% Particle filter
%
% Default is that we do not use the non linear kalman 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;
particle.initial_state_prior_std = .0001;
% Set the default order of approximation of the model (perturbation).
particle.perturbation = 2;
% Set the default number of particles.
particle_filter.number_of_particles = 500;
particle.number_of_particles = 500;
% Set the default approximation order (Smolyak)
particle_filter.smolyak_accuracy = 3;
particle.smolyak_accuracy = 3;
% By default we don't use pruning
particle_filter.pruning = 0;
particle.pruning = 0;
% Set default algorithm
particle_filter.algorithm = 'sequential_importance_particle_filter';
% Set the Gaussian approximation method
particle_filter.approximation_method = 'unscented';
particle.algorithm = 'sequential_importance_particle_filter';
% Set the Gaussian approximation method
particle.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 ;
particle.unscented.alpha = 1;
particle.unscented.beta = 2;
particle.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' ;
particle.resampling.status = 'systematic'; % 'generic'
particle.resampling.neff_threshold = .5;
% Choice of the resampling method
particle.resampling.method1 = 'traditional' ;
particle.resampling.method2 = 'kitagawa';
% Configuration of the mixture filters
particle.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 ;
particle.mixture_state_variables = 5 ;
particle.mixture_structural_shocks = 1 ;
particle.mixture_measurement_shocks = 1 ;
% Copy ep structure in options_ global structure
options_.particle = particle;
@ -205,9 +220,9 @@ options_.solve_algo = 2;
options_.linear = 0;
options_.replic = 50;
options_.drop = 100;
% if mjdgges.dll (or .mexw32 or ....) doesn't exist, matlab/qz is added to the path.
% There exists now qz/mjdgges.m that contains the calls to the old Sims code
% Hence, if mjdgges.m is visible exist(...)==2,
% if mjdgges.dll (or .mexw32 or ....) doesn't exist, matlab/qz is added to the path.
% There exists now qz/mjdgges.m that contains the calls to the old Sims code
% Hence, if mjdgges.m is visible exist(...)==2,
% this means that the DLL isn't avaiable and use_qzdiv is set to 1
if exist('mjdgges','file')==2
options_.use_qzdiv = 1;
@ -225,9 +240,9 @@ options_.ramsey_policy = 0;
options_.timeless = 0;
% estimation
estimation_info.prior = struct('name', {}, 'shape', {}, 'mean', {}, ...
'mode', {}, 'stdev', {}, 'date1', {}, ...
'date2', {}, 'shift', {}, 'variance', {});
estimation_info.parameters.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', {});
@ -240,10 +255,12 @@ estimation_info.measurement_error.prior = struct('name', {}, 'shape', {}, 'mean'
estimation_info.measurement_error_corr.prior = struct('name', {}, 'shape', {}, 'mean', {}, ...
'mode', {}, 'stdev', {}, 'date1', {}, ...
'date2', {}, 'shift', {}, 'variance', {});
estimation_info.parameters.prior_index = {};
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.parameters.options_index = {};
estimation_info.measurement_error.options_index = {};
estimation_info.structural_innovation.options_index = {};
estimation_info.measurement_error_corr.options_index = {};
@ -322,6 +339,7 @@ options_.filter_covariance = 0;
options_.filter_decomposition = 0;
options_.selected_variables_only = 0;
options_.initialize_estimated_parameters_with_the_prior_mode = 0;
options_.estimation_dll = 0;
% Misc
options_.conf_sig = 0.6;
oo_.exo_simul = [];
@ -341,9 +359,19 @@ M_.bvar = [];
options_.homotopy_mode = 0;
options_.homotopy_steps = 1;
% Simplex routine (variation on Nelder Mead algorithm)
% Simplex optimization routine (variation on Nelder Mead algorithm).
options_.simplex = [];
% CMAES optimization routine.
cmaes.SaveVariables='off';
cmaes.DispFinal='on';
cmaes.WarnOnEqualFunctionValues='no';
cmaes.DispModulo='10';
cmaes.LogModulo='0';
cmaes.LogTime='0';
options_.cmaes = cmaes;
% prior analysis
options_.prior_mc = 20000;
options_.prior_analysis_endo_var_list = [];
@ -357,6 +385,14 @@ options_.use_dll = 0;
% model evaluated using bytecode.dll
options_.bytecode = 0;
% use a fixed point method to solve Sylvester equation (for large scale
% models)
options_.sylvester_fp = 0;
% use a fixed point method to solve Lyapunov equation (for large scale
% models)
options_.lyapunov_fp = 0;
% dates for historical time series
options_.initial_date.freq = 1;
options_.initial_date.period = 1;

25
matlab/gsa/Morris_Measure_Groups.m
View File

@ -1,6 +1,4 @@
function [SAmeas, OutMatrix] = Morris_Measure_Groups(NumFact, Sample, Output, p, Group)
% [SAmeas, OutMatrix] = Morris_Measure_Groups(NumFact, Sample, Output, p, Group)
%
% Given the Morris sample matrix, the output values and the group matrix compute the Morris measures
@ -11,12 +9,12 @@ function [SAmeas, OutMatrix] = Morris_Measure_Groups(NumFact, Sample, Output, p,
% Each column represents one group.
% The element of each column are zero if the factor is not in the
% group. Otherwise it is 1.
%
% Sample := Matrix of the Morris sampled trajectories
%
% Output := Matrix of the output(s) values in correspondence of each point
% of each trajectory
%
% k = Number of factors
% -------------------------------------------------------------------------
% OUTPUTS
@ -24,6 +22,23 @@ function [SAmeas, OutMatrix] = Morris_Measure_Groups(NumFact, Sample, Output, p,
% for each output it gives the three measures of each factor
% -------------------------------------------------------------------------
% 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/>.
if nargin==0,
disp(' ')
disp('[SAmeas, OutMatrix] = Morris_Measure_Groups(NumFact, Sample, Output, p, Group);')

17
matlab/gsa/Sampling_Function_2.m
View File

@ -50,8 +50,23 @@ function [Outmatrix, OutFact] = Sampling_Function_2(p, k, r, UB, LB, GroupMat)
%
% F. Campolongo, J. Cariboni, JRC - IPSC Ispra, Varese, IT
% Last Update: 15 November 2005 by J.Cariboni
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 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/>.
% Parameters and initialisation of the output matrix
sizea = k;

17
matlab/gsa/cumplot.m
View File

@ -1,7 +1,22 @@
function h = cumplot(x);
%function h =cumplot(x)
% Copyright (C) 2005 Marco Ratto
% 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/>.
n=length(x);
x=[-inf; sort(x); Inf];

26
matlab/gsa/dat_fil_.m
View File

@ -1,20 +1,28 @@
function c = dat_fil_(data_file);
%
% Part of the Sensitivity Analysis Toolbox for DYNARE
%
% Written by Marco Ratto, 2006
% Written by Marco Ratto
% Joint Research Centre, The European Commission,
% (http://eemc.jrc.ec.europa.eu/),
% marco.ratto@jrc.it
%
% Disclaimer: This software is not subject to copyright protection and is in the public domain.
% It is an experimental system. The Joint Research Centre of European Commission
% assumes no responsibility whatsoever for its use by other parties
% and makes no guarantees, expressed or implied, about its quality, reliability, or any other
% characteristic. We would appreciate acknowledgement if the software is used.
% Reference:
% M. Ratto, Global Sensitivity Analysis for Macroeconomic models, MIMEO, 2006.
% 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/>.
try
eval(data_file);

27
matlab/gsa/filt_mc_.m
View File

@ -1,6 +1,5 @@
function [rmse_MC, ixx] = filt_mc_(OutDir,data_info)
% function [rmse_MC, ixx] = filt_mc_(OutDir)
% copyright Marco Ratto 2006
% inputs (from opt_gsa structure)
% vvarvecm = options_gsa_.var_rmse;
% loadSA = options_gsa_.load_rmse;
@ -10,21 +9,29 @@ function [rmse_MC, ixx] = filt_mc_(OutDir,data_info)
% istart = options_gsa_.istart_rmse;
% alphaPC = 0.5;
%
% Part of the Sensitivity Analysis Toolbox for DYNARE
%
% Written by Marco Ratto, 2006
% Written by Marco Ratto
% Joint Research Centre, The European Commission,
% (http://eemc.jrc.ec.europa.eu/),
% marco.ratto@jrc.it
%
% Disclaimer: This software is not subject to copyright protection and is in the public domain.
% It is an experimental system. The Joint Research Centre of European Commission
% assumes no responsibility whatsoever for its use by other parties
% and makes no guarantees, expressed or implied, about its quality, reliability, or any other
% characteristic. We would appreciate acknowledgement if the software is used.
% Reference:
% M. Ratto, Global Sensitivity Analysis for Macroeconomic models, MIMEO, 2006.
% 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/>.
global bayestopt_ estim_params_ M_ options_ oo_

28
matlab/gsa/ghx2transition.m
View File

@ -1,24 +1,30 @@
function [A,B] = ghx2transition(mm,iv,ic,aux)
% [A,B] = ghx2transition(mm,iv,ic,aux)
%
% Adapted by M. Ratto from kalman_transition_matrix.m
% (kalman_transition_matrix.m is part of DYNARE, copyright M. Juillard)
%
% Part of the Sensitivity Analysis Toolbox for DYNARE
%
% Written by Marco Ratto, 2006
% Adapted by M. Ratto (from kalman_transition_matrix.m)
% Joint Research Centre, The European Commission,
% (http://eemc.jrc.ec.europa.eu/),
% marco.ratto@jrc.it
%
% Disclaimer: This software is not subject to copyright protection and is in the public domain.
% It is an experimental system. The Joint Research Centre of European Commission
% assumes no responsibility whatsoever for its use by other parties
% and makes no guarantees, expressed or implied, about its quality, reliability, or any other
% characteristic. We would appreciate acknowledgement if the software is used.
% Reference:
% M. Ratto, Global Sensitivity Analysis for Macroeconomic models, MIMEO, 2006.
% 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/>.
global oo_ M_

2
matlab/gsa/gsa_plotmatrix.m
View File

@ -2,7 +2,7 @@ function gsa_plotmatrix(type,varargin)
% function gsa_plotmatrix(type,varargin)
% extended version of the standard MATLAB plotmatrix
% Copyright (C) 2011-2011 Dynare Team
% Copyright (C) 2011-2012 Dynare Team
%
% This file is part of Dynare.
%

17
matlab/gsa/gsa_skewness.m
View File

@ -1,5 +1,22 @@
function s=gsa_skewness(y),
% 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/>.
% y=stand_(y);
% s=mean(y.^3);
m2=mean((y-mean(y)).^2);

25
matlab/gsa/gsa_speed.m
View File

@ -1,21 +1,30 @@
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
%
% Written by Marco Ratto, 2006
% Written by Marco Ratto
% Joint Research Centre, The European Commission,
% (http://eemc.jrc.ec.europa.eu/),
% marco.ratto@jrc.it
%
% Disclaimer: This software is not subject to copyright protection and is in the public domain.
% It is an experimental system. The Joint Research Centre of European Commission
% assumes no responsibility whatsoever for its use by other parties
% and makes no guarantees, expressed or implied, about its quality, reliability, or any other
% characteristic. We would appreciate acknowledgement if the software is used.
% Reference:
% M. Ratto, Global Sensitivity Analysis for Macroeconomic models, MIMEO, 2006.
% 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/>.
nvar=length(mf);
nstate= size(A,1);

17
matlab/gsa/log_trans_.m
View File

@ -1,5 +1,22 @@
function [yy, xdir, isig, lam]=log_trans_(y0,xdir0)
% 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/>.
if nargin==1,
xdir0='';
end

18
matlab/gsa/map_ident_.m
View File

@ -1,4 +1,22 @@
function map_ident_(OutputDirectoryName)
% 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/>.
global bayestopt_ M_ options_ estim_params_ oo_
opt_gsa = options_.opt_gsa;

18
matlab/gsa/mc_moments.m
View File

@ -1,4 +1,22 @@
function [vdec, cc, ac] = mc_moments(mm, ss, dr)
% 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/>.
global options_ M_
[nr1, nc1, nsam] = size(mm);

19
matlab/gsa/myboxplot.m
View File

@ -1,8 +1,23 @@
function sout = myboxplot (data,notched,symbol,vertical,maxwhisker)
% sout = myboxplot (data,notched,symbol,vertical,maxwhisker)
% 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/>.
% % % % endif
if nargin < 5 | isempty(maxwhisker), maxwhisker = 1.5; end
if nargin < 4 | isempty(vertical), vertical = 1; end

18
matlab/gsa/myprctilecol.m
View File

@ -1,4 +1,22 @@
function y = myprctilecol(x,p);
% 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/>.
xx = sort(x);
[m,n] = size(x);

26
matlab/gsa/prior_draw_gsa.m
View File

@ -17,22 +17,30 @@ function pdraw = prior_draw_gsa(init,rdraw)
% SPECIAL REQUIREMENTS
% MATLAB Statistics Toolbox
%
%
% Part of the Sensitivity Analysis Toolbox for DYNARE
%
% Written by Marco Ratto, 2006
% Written by Marco Ratto
% Joint Research Centre, The European Commission,
% (http://eemc.jrc.ec.europa.eu/),
% marco.ratto@jrc.it
%
% Disclaimer: This software is not subject to copyright protection and is in the public domain.
% It is an experimental system. The Joint Research Centre of European Commission
% assumes no responsibility whatsoever for its use by other parties
% and makes no guarantees, expressed or implied, about its quality, reliability, or any other
% characteristic. We would appreciate acknowledgement if the software is used.
% Reference:
% M. Ratto, Global Sensitivity Analysis for Macroeconomic models, MIMEO, 2006.
% 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/>.
% global M_ options_ estim_params_ bayestopt_
global bayestopt_

17
matlab/gsa/priorcdf.m
View File

@ -8,6 +8,23 @@ function [xcum] = priorcdf(para, pshape, p6, p7, p3, p4)
% 5 is UNIFORM [p1,p2]
% Adapted by M. Ratto from MJ priordens.m
% 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/>.
nprio = length(pshape);
i = 1;

25
matlab/gsa/read_data.m
View File

@ -1,19 +1,28 @@
function [gend, data] = read_data
% Part of the Sensitivity Analysis Toolbox for DYNARE
%
% Written by Marco Ratto, 2006
% Written by Marco Ratto
% Joint Research Centre, The European Commission,
% (http://eemc.jrc.ec.europa.eu/),
% marco.ratto@jrc.it
%
% Disclaimer: This software is not subject to copyright protection and is in the public domain.
% It is an experimental system. The Joint Research Centre of European Commission
% assumes no responsibility whatsoever for its use by other parties
% and makes no guarantees, expressed or implied, about its quality, reliability, or any other
% characteristic. We would appreciate acknowledgement if the software is used.
% Reference:
% M. Ratto, Global Sensitivity Analysis for Macroeconomic models, MIMEO, 2006.
% 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/>.
global options_ bayestopt_

26
matlab/gsa/redform_map.m
View File

@ -11,21 +11,31 @@ function redform_map(dirname)
% ksstat = options_gsa_.ksstat_redform;
% alpha2 = options_gsa_.alpha2_redform;
%
% Part of the Sensitivity Analysis Toolbox for DYNARE
%
% Written by Marco Ratto, 2006
% Written by Marco Ratto
% Joint Research Centre, The European Commission,
% (http://eemc.jrc.ec.europa.eu/),
% marco.ratto@jrc.it
%
% Disclaimer: This software is not subject to copyright protection and is in the public domain.
% It is an experimental system. The Joint Research Centre of European Commission
% assumes no responsibility whatsoever for its use by other parties
% and makes no guarantees, expressed or implied, about its quality, reliability, or any other
% characteristic. We would appreciate acknowledgement if the software is used.
% Reference:
% M. Ratto, Global Sensitivity Analysis for Macroeconomic models, MIMEO, 2006.
% 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/>.
global M_ oo_ estim_params_ options_ bayestopt_

25
matlab/gsa/redform_screen.m
View File

@ -6,21 +6,30 @@ function redform_screen(dirname)
% anamexo = options_gsa_.namexo;
% iload = options_gsa_.load_redform;
%
% Part of the Sensitivity Analysis Toolbox for DYNARE
%
% Written by Marco Ratto, 2006
% Written by Marco Ratto
% Joint Research Centre, The European Commission,
% (http://eemc.jrc.ec.europa.eu/),
% marco.ratto@jrc.it
%
% Disclaimer: This software is not subject to copyright protection and is in the public domain.
% It is an experimental system. The Joint Research Centre of European Commission
% assumes no responsibility whatsoever for its use by other parties
% and makes no guarantees, expressed or implied, about its quality, reliability, or any other
% characteristic. We would appreciate acknowledgement if the software is used.
% Reference:
% M. Ratto, Global Sensitivity Analysis for Macroeconomic models, MIMEO, 2006.
% 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/>.
global M_ oo_ estim_params_ options_ bayestopt_

20
matlab/gsa/set_shocks_param.m
View File

@ -1,5 +1,23 @@
function set_shocks_param(xparam1)
global estim_params_ M_
% 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/>.
global estim_params_ M_
nvx = estim_params_.nvx;
ncx = estim_params_.ncx;

27
matlab/gsa/smirnov.m
View File

@ -2,23 +2,30 @@ function [H,prob,d] = smirnov(x1 , x2 , alpha, iflag )
% Smirnov test for 2 distributions
% [H,prob,d] = smirnov(x1 , x2 , alpha, iflag )
%
% Part of the Sensitivity Analysis Toolbox for DYNARE
%
% Written by Marco Ratto, 2006
% Written by Marco Ratto
% Joint Research Centre, The European Commission,
% (http://eemc.jrc.ec.europa.eu/),
% marco.ratto@jrc.it
%
% Disclaimer: This software is not subject to copyright protection and is in the public domain.
% It is an experimental system. The Joint Research Centre of European Commission
% assumes no responsibility whatsoever for its use by other parties
% and makes no guarantees, expressed or implied, about its quality, reliability, or any other
% characteristic. We would appreciate acknowledgement if the software is used.
% Reference:
% M. Ratto, Global Sensitivity Analysis for Macroeconomic models, MIMEO, 2006.
% 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/>.
if nargin<3
alpha = 0.05;

39
matlab/gsa/stab_map_.m
View File

@ -30,21 +30,30 @@ function x0 = stab_map_(OutputDirectoryName)
%
% USES qmc_sequence, stab_map_1, stab_map_2
%
% Part of the Sensitivity Analysis Toolbox for DYNARE
%
% Written by Marco Ratto, 2006
% Written by Marco Ratto
% Joint Research Centre, The European Commission,
% (http://eemc.jrc.ec.europa.eu/),
% marco.ratto@jrc.it
%
% Disclaimer: This software is not subject to copyright protection and is in the public domain.
% It is an experimental system. The Joint Research Centre of European Commission
% assumes no responsibility whatsoever for its use by other parties
% and makes no guarantees, expressed or implied, about its quality, reliability, or any other
% characteristic. We would appreciate acknowledgement if the software is used.
% Reference:
% M. Ratto, Global Sensitivity Analysis for Macroeconomic models, MIMEO, 2006.
% 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/>.
%global bayestopt_ estim_params_ dr_ options_ ys_ fname_
global bayestopt_ estim_params_ options_ oo_ M_
@ -79,6 +88,7 @@ nshock = nshock + estim_params_.nvn;
nshock = nshock + estim_params_.ncx;
nshock = nshock + estim_params_.ncn;
lpmat0=[];
xparam1=[];
pshape = bayestopt_.pshape(nshock+1:end);
p1 = bayestopt_.p1(nshock+1:end);
@ -387,6 +397,9 @@ else
end
load(filetoload,'lpmat','lpmat0','iunstable','istable','iindeterm','iwrong','egg','yys','nspred','nboth','nfwrd')
Nsam = size(lpmat,1);
if pprior==0,
eval(['load ' options_.mode_file '.mat;']);
end
if prepSA & isempty(strmatch('T',who('-file', filetoload),'exact')),
@ -522,18 +535,18 @@ if length(iunstable)>0 & length(iunstable)<Nsam,
c0=corrcoef(lpmat(istable,:));
c00=tril(c0,-1);
stab_map_2(lpmat(istable,:),alpha2, pvalue_corr, asname, OutputDirectoryName);
stab_map_2(lpmat(istable,:),alpha2, pvalue_corr, asname, OutputDirectoryName,xparam1);
if length(iunstable)>10,
stab_map_2(lpmat(iunstable,:),alpha2, pvalue_corr, auname, OutputDirectoryName);
stab_map_2(lpmat(iunstable,:),alpha2, pvalue_corr, auname, OutputDirectoryName,xparam1);
end
if length(iindeterm)>10,
stab_map_2(lpmat(iindeterm,:),alpha2, pvalue_corr, aindname, OutputDirectoryName);
stab_map_2(lpmat(iindeterm,:),alpha2, pvalue_corr, aindname, OutputDirectoryName,xparam1);
end
if length(ixun)>10,
stab_map_2(lpmat(ixun,:),alpha2, pvalue_corr, aunstname, OutputDirectoryName);
stab_map_2(lpmat(ixun,:),alpha2, pvalue_corr, aunstname, OutputDirectoryName,xparam1);
end
if length(iwrong)>10,
stab_map_2(lpmat(iwrong,:),alpha2, pvalue_corr, awrongname, OutputDirectoryName);
stab_map_2(lpmat(iwrong,:),alpha2, pvalue_corr, awrongname, OutputDirectoryName,xparam1);
end
x0=0.5.*(bayestopt_.ub(1:nshock)-bayestopt_.lb(1:nshock))+bayestopt_.lb(1:nshock);

25
matlab/gsa/stab_map_1.m
View File

@ -16,21 +16,30 @@ function [proba, dproba] = stab_map_1(lpmat, ibehaviour, inonbehaviour, aname, i
% solid lines for NON BEHAVIOURAL
% USES smirnov
%
% Part of the Sensitivity Analysis Toolbox for DYNARE
%
% Written by Marco Ratto, 2006
% Written by Marco Ratto
% Joint Research Centre, The European Commission,
% (http://eemc.jrc.ec.europa.eu/),
% marco.ratto@jrc.it
%
% Disclaimer: This software is not subject to copyright protection and is in the public domain.
% It is an experimental system. The Joint Research Centre of European Commission
% assumes no responsibility whatsoever for its use by other parties
% and makes no guarantees, expressed or implied, about its quality, reliability, or any other
% characteristic. We would appreciate acknowledgement if the software is used.
% Reference:
% M. Ratto, Global Sensitivity Analysis for Macroeconomic models, MIMEO, 2006.
% 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/>.
global estim_params_ bayestopt_ M_ options_

36
matlab/gsa/stab_map_2.m
View File

@ -1,21 +1,30 @@
function stab_map_2(x,alpha2, pvalue, fnam, dirname)
function stab_map_2(x,alpha2, pvalue, fnam, dirname,xparam1)
% function stab_map_2(x, alpha2, pvalue, fnam, dirname)
%
% Part of the Sensitivity Analysis Toolbox for DYNARE
%
% Written by Marco Ratto, 2006
% Written by Marco Ratto
% Joint Research Centre, The European Commission,
% (http://eemc.jrc.ec.europa.eu/),
% marco.ratto@jrc.it
%
% Disclaimer: This software is not subject to copyright protection and is in the public domain.
% It is an experimental system. The Joint Research Centre of European Commission
% assumes no responsibility whatsoever for its use by other parties
% and makes no guarantees, expressed or implied, about its quality, reliability, or any other
% characteristic. We would appreciate acknowledgement if the software is used.
% Reference:
% M. Ratto, Global Sensitivity Analysis for Macroeconomic models, MIMEO, 2006.
% 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/>.
%global bayestopt_ estim_params_ dr_ options_ ys_ fname_
global bayestopt_ estim_params_ options_ oo_ M_
@ -29,6 +38,9 @@ end
if nargin<5,
dirname='';
end
if nargin<6,
xparam1=[];
end
ys_ = oo_.dr.ys;
dr_ = oo_.dr;
@ -46,6 +58,9 @@ ifig=0;
j2=0;
if ishock==0
npar=estim_params_.np;
if ~isempty(xparam1),
xparam1=xparam1(nshock+1:end);
end
else
npar=estim_params_.np+nshock;
end
@ -68,6 +83,9 @@ for j=1:npar,
%plot(stock_par(ixx(nfilt+1:end,i),j),stock_par(ixx(nfilt+1:end,i),i2(jx)),'.k')
%hold on,
plot(x(:,j),x(:,i2(jx)),'.')
if ~isempty(xparam1)
hold on, plot(xparam1(j),xparam1(i2(jx)),'ro')
end
% xlabel(deblank(estim_params_.param_names(j,:)),'interpreter','none'),
% ylabel(deblank(estim_params_.param_names(i2(jx),:)),'interpreter','none'),
if ishock,

17
matlab/gsa/stand_.m
View File

@ -9,9 +9,24 @@ function [y, meany, stdy] = stand_(x)
% my: Vector of mean values for each column of y
% sy: Vector of standard deviations for each column of y
%
% Copyright (c) 2006 by JRC, European Commission, United Kingdom
% Author : Marco Ratto
% 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/>.
if nargin==0,
return;

6
matlab/gsa/tcrit.m
View File

@ -1,19 +1,15 @@
function t_crit = tcrit(n,pval0)
% function t_crit = tcrit(n,pval0)
%
% given the p-value pval0, the function givese the
% critical value t_crit of the t-distribution with n degress of freedom
%
% Part of the Sensitivity Analysis Toolbox for DYNARE
%
% Written by Marco Ratto
% Joint Research Centre, The European Commission,
% (http://eemc.jrc.ec.europa.eu/),
% marco.ratto@jrc.it
%
% Copyright (C) 2011-2011 Dynare Team
% Copyright (C) 2011-2012 Dynare Team
%
% This file is part of Dynare.
%

26
matlab/gsa/teff.m
View File

@ -1,22 +1,30 @@
function [yt, j0, ir, ic]=teff(T,Nsam,istable)
%
% [yt, j0, ir, ic]=teff(T,Nsam,istable)
%
% Part of the Sensitivity Analysis Toolbox for DYNARE
%
% Written by Marco Ratto, 2006
% Written by Marco Ratto
% Joint Research Centre, The European Commission,
% (http://eemc.jrc.ec.europa.eu/),
% marco.ratto@jrc.it
%
% Disclaimer: This software is not subject to copyright protection and is in the public domain.
% It is an experimental system. The Joint Research Centre of European Commission
% assumes no responsibility whatsoever for its use by other parties
% and makes no guarantees, expressed or implied, about its quality, reliability, or any other
% characteristic. We would appreciate acknowledgement if the software is used.
% Reference:
% M. Ratto, Global Sensitivity Analysis for Macroeconomic models, MIMEO, 2006.
% 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/>.
ndim = (length(size(T)));
if ndim==3,

20
matlab/gsa/th_moments.m
View File

@ -1,6 +1,22 @@
% Copyright (C) 2001 Michel Juillard
%
function [vdec, corr, autocorr, z, zz] = th_moments(dr,var_list)
% 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/>.
global M_ oo_ options_
nvar = size(var_list,1);

18
matlab/gsa/thet2tau.m
View File

@ -1,4 +1,22 @@
function tau = thet2tau(params, indx, indexo, flagmoments,mf,nlags,useautocorr)
% 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/>.
global M_ oo_ options_
if nargin==1,

25
matlab/gsa/trank.m
View File

@ -2,21 +2,30 @@ function yr = trank(y);
% yr = trank(y);
% yr is the rank transformation of y
%
% Part of the Sensitivity Analysis Toolbox for DYNARE
%
% Written by Marco Ratto, 2006
% Written by Marco Ratto
% Joint Research Centre, The European Commission,
% (http://eemc.jrc.ec.europa.eu/),
% marco.ratto@jrc.it
%
% Disclaimer: This software is not subject to copyright protection and is in the public domain.
% It is an experimental system. The Joint Research Centre of European Commission
% assumes no responsibility whatsoever for its use by other parties
% and makes no guarantees, expressed or implied, about its quality, reliability, or any other
% characteristic. We would appreciate acknowledgement if the software is used.
% Reference:
% M. Ratto, Global Sensitivity Analysis for Macroeconomic models, MIMEO, 2006.
% 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/>.
[nr, nc] = size(y);
for j=1:nc,

2
matlab/hessian.m
View File

@ -39,7 +39,7 @@ if ~isa(func, 'function_handle')
func = str2func(func);
end
n=size(x,1);
h1=max(abs(x),sqrt(gstep)*ones(n,1))*eps^(1/6);
h1=max(abs(x),sqrt(gstep(1))*ones(n,1))*eps^(1/6)*gstep(2);
h_1=h1;
xh1=x+h1;
h1=xh1-x;

2
matlab/homotopy1.m
View File

@ -77,5 +77,5 @@ for i=1:step_nbr+1
oo_.exo_steady_state(values(ix,2)) = points(ix,i);
oo_.exo_det_steady_state(values(ixd,2)) = points(ixd,i);
steady_(M_,options_,oo_);
oo_.steady_state = steady_(M_,options_,oo_);
end

2
matlab/homotopy2.m
View File

@ -100,6 +100,6 @@ for i = 1:nv
disp([ 'HOMOTOPY mode 2: lauching solver with ' deblank(varname) ' = ' num2str(v) ' ...'])
steady_(M_,options_,oo_);
oo_.steady_state = steady_(M_,options_,oo_);
end
end

2
matlab/homotopy3.m
View File

@ -91,7 +91,7 @@ while iter < step_nbr
old_ss = oo_.steady_state;
try
steady_(M_,options_,oo_);
oo_.steady_state = steady_(M_,options_,oo_);
if length([kplus; kminus]) == nv
return

6
matlab/initial_estimation_checks.m
View File

@ -42,6 +42,12 @@ if DynareOptions.dsge_var
[fval,cost_flag,info] = DsgeVarLikelihood(xparam1,DynareDataset,DynareOptions,Model,EstimatedParameters,BayesInfo,DynareResults);
else
[fval,cost_flag,ys,trend_coeff,info] = dsge_likelihood(xparam1,DynareDataset,DynareOptions,Model,EstimatedParameters,BayesInfo,DynareResults);
if DynareOptions.mode_compute == 5
% this call is necessary to initialized persistent variable
% 'penalty' in dsge_likelihood_hh
[fval,llik,cost_flag,ys,trend_coeff,info] = ...
dsge_likelihood_hh(xparam1,DynareDataset,DynareOptions,Model,EstimatedParameters,BayesInfo,DynareResults);
end
end
if info(1) > 0

11
matlab/k_order_pert.m
View File

@ -1,7 +1,7 @@
function [dr,info] = k_order_pert(dr,M,options,oo)
% Compute decision rules using the k-order DLL from Dynare++
% Copyright (C) 2009-2010 Dynare Team
% Copyright (C) 2009-2012 Dynare Team
%
% This file is part of Dynare.
%
@ -26,20 +26,17 @@ order = options.order;
switch(order)
case 1
[err, g_1] = k_order_perturbation(dr,M,options, ['.' ...
mexext]);
[err, g_1] = k_order_perturbation(dr,M,options);
mexErrCheck('k_order_perturbation', err);
dr.g_1 = g_1;
case 2
[err, g_0, g_1, g_2] = k_order_perturbation(dr,M,options, ['.' ...
mexext]);
[err, g_0, g_1, g_2] = k_order_perturbation(dr,M,options);
mexErrCheck('k_order_perturbation', err);
dr.g_0 = g_0;
dr.g_1 = g_1;
dr.g_2 = g_2;
case 3
[err, g_0, g_1, g_2, g_3] = k_order_perturbation(dr,M,options, ['.' ...
mexext]);
[err, g_0, g_1, g_2, g_3] = k_order_perturbation(dr,M,options);
mexErrCheck('k_order_perturbation', err);
dr.g_0 = g_0;
dr.g_1 = g_1;

11
matlab/kalman/likelihood/kalman_filter.m
View File

@ -169,9 +169,8 @@ if t<last
end
% Compute minus the log-likelihood.
if presample
if presample>=diffuse_periods
likk = likk(1+(presample-diffuse_periods):end);
end
end
LIK = sum(likk);
if presample > diffuse_periods
LIK = sum(likk(1+presample-diffuse_periods:end));
else
LIK = sum(likk);
end

12
matlab/kalman/likelihood/kalman_filter_d.m
View File

@ -111,14 +111,4 @@ end
dlik = dlik(1:s);
dlik = .5*(dlik + pp*log(2*pi));
if presample
if presample>=length(dlik)
dLIK = 0;
dlik = [];
else
dlik = dlik(1+presample:end);
dLIK = sum(dlik);% Minus the log-likelihood (for the initial periods).
end
else
dLIK = sum(dlik);
end
dLIK = sum(dlik(1+presample:end));

11
matlab/kalman/likelihood/missing_observations_kalman_filter.m
View File

@ -139,9 +139,8 @@ if t<last
end
% Compute minus the log-likelihood.
if presample
if presample>=diffuse_periods
lik = lik(1+(presample-diffuse_periods):end);
end
end
LIK = sum(lik);
if presample>=diffuse_periods
LIK = sum(lik(1+presample-diffuse_periods:end));
else
LIK = sum(lik);
end

12
matlab/kalman/likelihood/missing_observations_kalman_filter_d.m
View File

@ -125,14 +125,4 @@ end
dlik = .5*dlik(1:s);
if presample
if presample>=length(dlik)
dLIK = 0;
dlik = [];
else
dlik = dlik(1+presample:end);
dLIK = sum(dlik);% Minus the log-likelihood (for the initial periods).
end
else
dLIK = sum(dlik);
end
dLIK = sum(dlik(1+presample:end));

11
matlab/kalman/likelihood/univariate_kalman_filter.m
View File

@ -167,9 +167,8 @@ if t<last
end
% Compute minus the log-likelihood.
if presample
if presample>=diffuse_periods
lik = lik(1+(presample-diffuse_periods):end);
end
end
LIK = sum(lik);
if presample > diffuse_periods
LIK = sum(lik(1+presample-diffuse_periods:end));
else
LIK = sum(lik);
end

14
matlab/kalman/likelihood/univariate_kalman_filter_d.m
View File

@ -164,16 +164,4 @@ end
dlikk = .5*dlikk(1:s);
llik = .5*llik(1:s,:);
if presample
if presample>=length(dlik)
dLIK = 0;
dlikk= [];
llik = [];
else
dlikk= dlikk(1+presample:end);
llik = llik(1+presample:end,:);
dLIK = sum(dlikk);% Minus the log-likelihood (for the initial periods).
end
else
dLIK = sum(dlikk);
end
dLIK = sum(dlikk(1+presample:end));

52
matlab/lyapunov_symm.m
View File

@ -1,4 +1,4 @@
function [x,u] = lyapunov_symm(a,b,qz_criterium,lyapunov_complex_threshold,method)
function [x,u] = lyapunov_symm(a,b,qz_criterium,lyapunov_complex_threshold,method, R)
% Solves the Lyapunov equation x-a*x*a' = b, for b and x symmetric matrices.
% If a has some unit roots, the function computes only the solution of the stable subsystem.
%
@ -12,6 +12,7 @@ function [x,u] = lyapunov_symm(a,b,qz_criterium,lyapunov_complex_threshold,metho
% variables and the schur decomposition is triggered.
% method=2 then U, T, n and k are declared as persistent
% variables and the schur decomposition is not performed.
% method=3 fixed point method
% OUTPUTS
% x: [double] m*m solution matrix of the lyapunov equation, where m is the dimension of the stable subsystem.
% u: [double] Schur vectors associated with unit roots
@ -38,10 +39,55 @@ function [x,u] = lyapunov_symm(a,b,qz_criterium,lyapunov_complex_threshold,metho
%
% You should have received a copy of the GNU General Public License
% along with Dynare. If not, see <http://www.gnu.org/licenses/>.
if nargin<5
method = 0;
end
if method == 3
persistent X method1;
if ~isempty(method1)
method = method1;
end;
fprintf(' [methode=%d] ',method);
if method == 3
%tol = 1e-8;
tol = 1e-10;
it_fp = 0;
evol = 100;
if isempty(X)
X = b;
max_it_fp = 2000;
else
max_it_fp = 300;
end;
at = a';
%fixed point iterations
while evol > tol && it_fp < max_it_fp;
X_old = X;
X = a * X * at + b;
evol = max(sum(abs(X - X_old))); %norm_1
%evol = max(sum(abs(X - X_old)')); %norm_inf
it_fp = it_fp + 1;
end;
fprintf('lyapunov it_fp=%d evol=%g\n',it_fp,evol);
if it_fp >= max_it_fp
disp(['convergence not achieved in solution of Lyapunov equation after ' int2str(it_fp) ' iterations, switching method from 3 to 0']);
method1 = 0;
method = 0;
else
method1 = 3;
x = X;
return;
end;
end;
elseif method == 4
% works only with Matlab System Control toolbox
chol_b = R*chol(b,'lower');
Rx = dlyapchol(a,chol_b);
x = Rx' * Rx;
return;
end;
if method
persistent U T k n
else
@ -113,4 +159,4 @@ if i == 1
x(1,1) = (B(1,1)+c)/(1-T(1,1)*T(1,1));
end
x = U(:,k+1:end)*x*U(:,k+1:end)';
u = U(:,1:k);
u = U(:,1:k);

75
matlab/mode_check.m
View File

@ -1,24 +1,46 @@
function mode_check(fun,x,hessian,DynareDataset,DynareOptions,Model,EstimatedParameters,BayesInfo,DynareResults)
% Checks the estimated ML mode or Posterior mode.
% function mode_check(x,fval,hessian,gend,data,lb,ub)
% Checks the maximum likelihood mode
%
% INPUTS
% x: mode
% fval: value at the maximum likelihood mode
% hessian: matrix of second order partial derivatives
% gend: scalar specifying the number of observations
% data: matrix of data
% lb: lower bound
% ub: upper bound
%
% OUTPUTS
% none
%
% SPECIAL REQUIREMENTS
% none
%@info:
%! @deftypefn {Function File} mode_check (@var{fun}, @var{x}, @var{hessian}, @var{DynareDataset}, @var{DynareOptions}, @var{Model}, @var{EstimatedParameters}, @var{BayesInfo}, @var{DynareResults})
%! @anchor{mode_check}
%! @sp 1
%! Checks the estimated ML mode or Posterior mode by plotting sections of the likelihood/posterior kernel.
%! Each plot shows the variation of the likelihood implied by the variations of a single parameter, ceteris paribus)
%! @sp 2
%! @strong{Inputs}
%! @sp 1
%! @table @ @var
%! @item fun
%! Objective function.
%! @item x
%! Estimated mode.
%! @item start
%! Hessian of the objective function at the estimated mode @var{x}.
%! @item DynareDataset
%! Structure specifying the dataset used for estimation (dataset_).
%! @item DynareOptions
%! Structure defining dynare's options (options_).
%! @item Model
%! Structure specifying the (estimated) model (M_).
%! @item EstimatedParameters
%! Structure specifying the estimated parameters (estimated_params_).
%! @item BayesInfo
%! Structure containing information about the priors used for estimation (bayestopt_).
%! @item DynareResults
%! Structure gathering the results (oo_).
%! @end table
%! @sp 2
%! @strong{Outputs}
%! @sp 2
%! @strong{This function is called by:}
%! @sp 2
%! @strong{This function calls:}
%! The objective function (@var{func}).
%! @end deftypefn
%@eod:
% Copyright (C) 2003-2010 Dynare Team
% Copyright (C) 2003-2010, 2012 Dynare Team
%
% This file is part of Dynare.
%
@ -62,6 +84,8 @@ if TeX
fprintf(fidTeX,' \n');
end
ll = DynareOptions.mode_check_neighbourhood_size;
for plt = 1:nbplt,
if TeX
NAMES = [];
@ -82,9 +106,18 @@ for plt = 1:nbplt,
end
end
xx = x;
l1 = max(BayesInfo.lb(kk),0.5*x(kk));
l2 = min(BayesInfo.ub(kk),1.5*x(kk));
z = [l1:(l2-l1)/20:l2];
l1 = max(BayesInfo.lb(kk),(1-ll)*x(kk)); m1 = 0;
l2 = min(BayesInfo.ub(kk),(1+ll)*x(kk));
if l2<(1+ll)*x(kk)
l1 = x(kk) - (l2-x(kk));
m1 = 1;
end
if ~m1 && (l1>(1-ll)*x(kk)) && (x(kk)+(x(kk)-l1)<BayesInfo.ub(kk))
l2 = x(kk) + (x(kk)-l1);
end
z1 = l1:((x(kk)-l1)/10):x(kk);
z2 = x(kk):((l2-x(kk))/10):l2;
z = union(z1,z2);
if DynareOptions.mode_check_nolik==0,
y = zeros(length(z),2);
dy = priordens(xx,BayesInfo.pshape,BayesInfo.p6,BayesInfo.p7,BayesInfo.p3,BayesInfo.p4);

2
matlab/mr_gstep.m
View File

@ -59,6 +59,7 @@ while i<n
if gg(i)*(hh(i)*gg(i))/2 > htol
[f0 x fc retcode] = csminit(func0,x,f0,gg,0,diag(hh),DynareDataset,DynareOptions,Model,EstimatedParameters,BayesInfo,DynareResults);
ig(i)=1;
fprintf(['Done for param %s = %8.4f\n'],BayesInfo.name{i},x(i))
end
xh1=x;
end
@ -67,4 +68,3 @@ end
save gstep.mat x h1 f0

3
matlab/ms-sbvar/cstz/fn_foregraph.m
View File

@ -39,8 +39,7 @@ hornum = cell(length(vyrs),1); % horizontal year (number)
count=0;
for k=vyrs'
count=count+1;
jnk=num2str(k);
hornum{count}=jnk(3:4); % e.g., with '1990', we have '90'
hornum{count}=num2str(k);
end
count=0;

2
matlab/ms-sbvar/cstz/fn_rnrprior_covres_dobs.m
View File

@ -141,6 +141,8 @@ for i = 1:lags
sgpbid((i-1)*nvar+j) = lagdecay^2/sgsh(j); % ith equation
elseif (q_m==4)
sgpbid((i-1)*nvar+j) = (1/i^mu(4))^2/sgsh(j); % ith equation
elseif (q_m==1)
sgpbid((i-1)*nvar+j) = (1/(i*4)^mu(4))^2/sgsh(j); % ith equation
else
error('Incompatibility with lags, check the possible errors!!!')
%warning('Incompatibility with lags, check the possible errors!!!')

2
matlab/ms-sbvar/cstz/fn_rnrprior_covres_dobs_tv2.m
View File

@ -163,6 +163,8 @@ for i = 1:lags
sgpbid((i-1)*nvar+j) = lagdecay^2/sgsh(j); % ith equation
elseif (q_m==4)
sgpbid((i-1)*nvar+j) = (1/i^mu(4))^2/sgsh(j); % ith equation
elseif (q_m==1)
sgpbid((i-1)*nvar+j) = (1/(i*4)^mu(4))^2/sgsh(j); % ith equation
else
error('Incompatibility with lags, check the possible errors!!!')
%warning('Incompatibility with lags, check the possible errors!!!')

7
matlab/ms-sbvar/initialize_ms_sbvar_options.m
View File

@ -12,7 +12,7 @@ function options_=initialize_ms_sbvar_options(M_, options_)
% SPECIAL REQUIREMENTS
% none
% Copyright (C) 2011 Dynare Team
% Copyright (C) 2011-2012 Dynare Team
%
% This file is part of Dynare.
%
@ -100,12 +100,15 @@ options_.ms.real_time_smoothed_probabilities = 0;
% irf
options_.ms.horizon = 12;
options_.ms.filtered_probabilities = 0;
options_.ms.error_bands = 1;
options_.ms.percentiles = [.16 .5 .84];
options_.ms.parameter_uncertainty = 0;
options_.ms.shock_draws = 10000;
options_.ms.shocks_per_parameter = 10;
options_.ms.median = 0;
options_.ms.regime = 0;
options_.ms.regimes = 0;
% forecast
options_.ms.forecast_data_obs = 0;
% variance decomposition
options_.ms.error_bands = 1;
end

104
matlab/ms-sbvar/ms_forecast.m
View File

@ -14,7 +14,7 @@ function [options_, oo_]=ms_forecast(M_, options_, oo_)
% SPECIAL REQUIREMENTS
% none
% Copyright (C) 2011 Dynare Team
% Copyright (C) 2011-2012 Dynare Team
%
% This file is part of Dynare.
%
@ -34,46 +34,86 @@ function [options_, oo_]=ms_forecast(M_, options_, oo_)
disp('MS-SBVAR Forecasts');
options_ = set_file_tags(options_);
[options_, oo_] = set_ms_estimation_file(options_.ms.file_tag, options_, oo_);
options_ = set_ms_simulation_file(options_);
clean_files_for_second_type_of_mex(M_, options_, 'forecast')
clean_ms_forecast_files(options_.ms.output_file_tag);
forecastdir = [options_.ms.output_file_tag filesep 'Forecast'];
create_dir(forecastdir);
opt = { ...
{'file_tag', options_.ms.file_tag}, ...
{'seed', options_.DynareRandomStreams.seed}, ...
{'horizon', options_.ms.horizon}, ...
{'number_observations', options_.ms.forecast_data_obs}, ...
{'error_bands', options_.ms.error_bands}, ...
{'percentiles', options_.ms.percentiles}, ...
{'thin', options_.ms.thinning_factor}
};
% setup command line options
opt = ['-forecast -nodate -seed ' num2str(options_.DynareRandomStreams.seed)];
opt = [opt ' -do ' forecastdir];
opt = [opt ' -ft ' options_.ms.file_tag];
opt = [opt ' -fto ' options_.ms.output_file_tag];
opt = [opt ' -horizon ' num2str(options_.ms.horizon)];
opt = [opt ' -thin ' num2str(options_.ms.thinning_factor)];
opt = [opt ' -data ' num2str(options_.ms.forecast_data_obs)];
if options_.ms.regimes
opt = [opt ' -regimes'];
elseif options_.ms.regime
% regime-1 since regime is 0-indexed in C but 1-indexed in Matlab
opt = [opt ' -regime ' num2str(options_.ms.regime-1)];
end
if options_.ms.parameter_uncertainty
options_ = set_ms_simulation_file(options_);
opt = [opt ' -parameter_uncertainty'];
opt = [opt ' -shocks_per_parameter ' num2str(options_.ms.shocks_per_parameter)];
else
opt = [opt ' -shocks_per_parameter ' num2str(options_.ms.shock_draws)];
end
percentiles_size = 0;
if options_.ms.median
opt = [opt(:)' {{'median'}}];
percentiles_size = 1;
opt = [opt ' -percentiles ' num2str(percentiles_size) ' 0.5'];
else
percentiles_size = size(options_.ms.percentiles,2);
opt = [opt ' -percentiles ' num2str(percentiles_size)];
for i=1:size(options_.ms.percentiles,2)
opt = [opt ' ' num2str(options_.ms.percentiles(i))];
end
end
[err, forecast] = mex_ms_forecast([opt(:)', {{'free_parameters',oo_.ms.maxparams}, ...
{'shocks_per_parameter', options_.ms.shock_draws}}]);
mexErrCheck('mex_ms_forecast ergodic ', err);
plot_ms_forecast(M_,options_,forecast,'Forecast',options_.graph_save_formats,options_.TeX);
% forecast
[err] = ms_sbvar_command_line(opt);
mexErrCheck('ms_forecast',err);
[err, regime_forecast] = mex_ms_forecast([opt(:)', {{'free_parameters',oo_.ms.maxparams}, ...
{'shocks_per_parameter', options_.ms.shock_draws}, {'regimes'}}]);
mexErrCheck('mex_ms_forecast ergodic regimes', err);
save([forecastdir filesep 'ergodic_forecast.mat'], 'forecast', 'regime_forecast');
% Plot Forecasts
if options_.ms.regimes
n_chains = length(options_.ms.ms_chain);
n_regimes=1;
for i_chain=1:n_chains
n_regimes = n_regimes*length(options_.ms.ms_chain(i_chain).regime);
end
if exist(options_.ms.mh_file,'file') > 0
[err, forecast] = mex_ms_forecast([opt(:)', {{'free_parameters',oo_.ms.maxparams}, ...
{'shocks_per_parameter', options_.ms.shocks_per_parameter}, ...
{'simulation_file', options_.ms.mh_file}, {'parameter_uncertainty'}}]);
mexErrCheck('mex_ms_forecast bayesian ', err);
plot_ms_forecast(M_,options_,forecast,'Forecast w/ Parameter Uncertainty',options_.graph_save_formats,options_.TeX);
for regime_i=1:n_regimes
forecast_title = ['Forecast, Regimes ' num2str(regime_i)];
forecast_data = load([forecastdir filesep 'forecasts_percentiles_regime_' ...
num2str(regime_i-1) '_' options_.ms.output_file_tag ...
'.out'], '-ascii');
forecast_data = reshape_ascii_forecast_data(M_.endo_nbr, ...
percentiles_size, options_.ms.horizon, forecast_data);
save([forecastdir filesep 'forecast_regime_' num2str(regime_i-1)], ...
'forecast_data');
plot_ms_forecast(M_, options_, forecast_data, forecast_title);
end
else
if options_.ms.regime
forecast_data = load([forecastdir filesep 'forecasts_percentiles_regime_' ...
num2str(options_.ms.regime-1) '_' options_.ms.output_file_tag ...
'.out'], '-ascii');
forecast_title = ['Forecast, Regime ' num2str(options_.ms.regime)];
save_filename = ['forecast_regime_' num2str(options_.ms.regime-1)];
else
forecast_data = load([forecastdir filesep 'forecasts_percentiles_' ...
options_.ms.output_file_tag '.out'], '-ascii');
forecast_title = 'Forecast';
save_filename = 'forecast';
end
[err, regime_forecast] = mex_ms_forecast([opt(:)', {{'free_parameters',oo_.ms.maxparams}, ...
{'shocks_per_parameter', options_.ms.shocks_per_parameter}, ...
{'simulation_file', options_.ms.mh_file}, {'parameter_uncertainty','regimes'}}]);
mexErrCheck('mex_ms_forecast bayesian regimes ', err);
save([forecastdir filesep 'bayesian_forecast.mat'], 'forecast', 'regime_forecast');
forecast_data = reshape_ascii_forecast_data(M_.endo_nbr, ...
percentiles_size, options_.ms.horizon, forecast_data);
save([forecastdir filesep save_filename], 'forecast_data');
plot_ms_forecast(M_, options_, forecast_data, forecast_title);
end
end

121
matlab/ms-sbvar/ms_irf.m
View File

@ -15,7 +15,7 @@ function [options_, oo_]=ms_irf(varlist, M_, options_, oo_)
% SPECIAL REQUIREMENTS
% none
% Copyright (C) 2011 Dynare Team
% Copyright (C) 2011-2012 Dynare Team
%
% This file is part of Dynare.
%
@ -35,50 +35,91 @@ function [options_, oo_]=ms_irf(varlist, M_, options_, oo_)
disp('MS-SBVAR Impulse Response Function');
options_ = set_file_tags(options_);
[options_, oo_] = set_ms_estimation_file(options_.ms.file_tag, options_, oo_);
options_ = set_ms_simulation_file(options_);
clean_files_for_second_type_of_mex(M_, options_, 'irf')
clean_ms_irf_files(options_.ms.output_file_tag);
irfdir = [options_.ms.output_file_tag filesep 'IRF'];
create_dir(irfdir);
opt = { ...
{'file_tag', options_.ms.file_tag}, ...
{'seed', options_.DynareRandomStreams.seed}, ...
{'horizon', options_.ms.horizon}, ...
{'filtered', options_.ms.filtered_probabilities}, ...
{'error_bands', options_.ms.error_bands}, ...
{'percentiles', options_.ms.percentiles}, ...
{'thin', options_.ms.thinning_factor}
};
% setup command line options
opt = ['-ir -seed ' num2str(options_.DynareRandomStreams.seed)];
opt = [opt ' -do ' irfdir];
opt = [opt ' -ft ' options_.ms.file_tag];
opt = [opt ' -fto ' options_.ms.output_file_tag];
opt = [opt ' -horizon ' num2str(options_.ms.horizon)];
opt = [opt ' -thin ' num2str(options_.ms.thinning_factor)];
if options_.ms.regimes
opt = [opt ' -regimes'];
elseif options_.ms.regime
% regime-1 since regime is 0-indexed in C but 1-indexed in Matlab
opt = [opt ' -regime ' num2str(options_.ms.regime-1)];
elseif options_.ms.filtered_probabilities
opt = [opt ' -filtered'];
end
if options_.ms.parameter_uncertainty
options_ = set_ms_simulation_file(options_);
opt = [opt ' -parameter_uncertainty'];
opt = [opt ' -shocks_per_parameter ' num2str(options_.ms.shocks_per_parameter)];
else
opt = [opt ' -shocks_per_parameter ' num2str(options_.ms.shock_draws)];
end
percentiles_size = 0;
if options_.ms.median
opt = [opt(:)' {{'median'}}];
end
[err, irf] = mex_ms_irf([opt(:)', {{'free_parameters', oo_.ms.maxparams}, {'shocks_per_parameter', options_.ms.shock_draws}}]);
mexErrCheck('mex_ms_irf ergodic ', err);
plot_ms_irf(M_,options_,irf,options_.varobs,'Ergodic Impulse Responses',varlist);
[err, regime_irfs] = mex_ms_irf([opt(:)', {{'free_parameters',oo_.ms.maxparams}, {'shocks_per_parameter', options_.ms.shock_draws}, {'regimes'}}]);
mexErrCheck('mex_ms_irf ergodic regimes ',err);
for i=1:size(regime_irfs,1)
plot_ms_irf(M_,options_,squeeze(regime_irfs(i,:,:,:)),options_.varobs,['Ergodic ' ...
'Impulse Responses State ' int2str(i)],varlist);
end
save([irfdir filesep 'ergodic_irf.mat'], 'irf', 'regime_irfs');
if exist(options_.ms.mh_file,'file') > 0
[err, irf] = mex_ms_irf([opt(:)', {{'shocks_per_parameter', options_.ms.shocks_per_parameter}, ...
{'parameter_uncertainty'},{'simulation_file',options_.ms.mh_file}}]);
mexErrCheck('mex_ms_irf bayesian ',err);
plot_ms_irf(M_,options_,irf,options_.varobs,'Impulse Responses with Parameter Uncertainty',varlist);
[err, regime_irfs] = mex_ms_irf([opt(:)', {{'shocks_per_parameter', options_.ms.shocks_per_parameter}, ...
{'simulation_file',options_.ms.mh_file},{'parameter_uncertainty'},{'regimes'}}]);
mexErrCheck('mex_ms_irf bayesian regimes ',err);
for i=1:size(regime_irfs,1)
plot_ms_irf(M_,options_,squeeze(regime_irfs(i,:,:,:)),options_.varobs,['Impulse ' ...
'Responses with Parameter Uncertainty State ' int2str(i)],varlist);
percentiles_size = 1;
opt = [opt ' -percentiles ' num2str(percentiles_size) ' 0.5'];
else
percentiles_size = size(options_.ms.percentiles,2);
opt = [opt ' -percentiles ' num2str(percentiles_size)];
for i=1:size(options_.ms.percentiles,2)
opt = [opt ' ' num2str(options_.ms.percentiles(i))];
end
save([irfdir filesep 'bayesian_irf.mat'], 'irf', 'regime_irfs');
end
% irf
[err] = ms_sbvar_command_line(opt);
mexErrCheck('ms_irf',err);
% Plot IRFs
if options_.ms.regimes
n_chains = length(options_.ms.ms_chain);
n_regimes=1;
for i_chain=1:n_chains
n_regimes = n_regimes*length(options_.ms.ms_chain(i_chain).regime);
end
for regime_i=1:n_regimes
irf_title = ['Impulse Responses, Regime ' num2str(regime_i)];
irf_data = load([irfdir filesep 'ir_percentiles_regime_' ...
num2str(regime_i-1) '_' options_.ms.output_file_tag ...
'.out'], '-ascii');
irf_data = reshape_ascii_irf_data(M_.endo_nbr, percentiles_size, ...
options_.ms.horizon, irf_data);
save([irfdir filesep 'irf_regime_' num2str(regime_i-1)], 'irf_data');
plot_ms_irf(M_, options_, irf_data, irf_title, varlist);
end
else
if options_.ms.regime
irf_data = load([irfdir filesep 'ir_percentiles_regime_' ...
num2str(options_.ms.regime-1) '_' options_.ms.output_file_tag ...
'.out'], '-ascii');
irf_title = ['Impulse Response, Regime ' num2str(options_.ms.regime)];
save_filename = ['irf_regime_' num2str(options_.ms.regime-1)];
elseif options_.ms.filtered_probabilities
irf_data = load([irfdir filesep 'ir_percentiles_filtered_' ...
options_.ms.output_file_tag '.out'], '-ascii');
irf_title = 'Impulse Response Filtered';
save_filename = 'irf';
else
irf_data = load([irfdir filesep 'ir_percentiles_ergodic_' ...
options_.ms.output_file_tag '.out'], '-ascii');
irf_title = 'Impulse Response Ergodic';
save_filename = 'irf';
end
irf_data = reshape_ascii_irf_data(M_.endo_nbr, percentiles_size, ...
options_.ms.horizon, irf_data);
save([irfdir filesep save_filename], 'irf_data');
plot_ms_irf(M_, options_, irf_data, irf_title, varlist);
end
end

394
matlab/ms-sbvar/ms_mardd.m
View File

@ -1,197 +1,197 @@
function ms_mardd(options_)
% Applies to both linear and exclusion restrictions.
% (1) Marginal likelihood function p(Y) for constant structural VAR models, using Chib (1995)'s ``Marginal Likelihood from the Gibbs Output'' in JASA.
% (2) Conditional likelihood function f(Y|A0, A+) on the ML estimate for constant exclusion-identified models.
% See Forecast (II) pp.67-80.
%
% Tao Zha, September 1999. Quick revisions, May 2003. Final revision, September 2004.
% Copyright (C) 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/>.
msstart2 % start the program in which everyhting is initialized through msstart2.m
if ~options_.ms.indxestima
warning(' ')
disp('You must set IxEstima=1 in msstart to run this program')
disp('Press ctrl-c to abort now')
pause
end
A0xhat = zeros(size(A0hat));
Apxhat = zeros(size(Aphat));
if (0)
%Robustness check to see if the same result is obtained with the purterbation of the parameters.
for k=1:nvar
bk = Uiconst{k}'*A0hat(:,k);
gk = Viconst{k}'*Aphat(:,k);
A0xhat(:,k) = Uiconst{k}*(bk + 5.2*randn(size(bk))); % Perturbing the posterior estimate.
Apxhat(:,k) = Viconst{k}*(gk + 5.2*randn(size(gk))); % Perturbing the posterior estimate.
end
else
%At the posterior estimate.
A0xhat = A0hat; % ML estimate of A0
Apxhat = Aphat; % ML estimate of A+
end
%--- Rename variables.
YatYa = yty;
XatYa = xty;
ytx = xty';
YatXa = ytx;
XatXa = xtx;
%--------- The log value of p(A0,A+) at some point such as the peak ----------
vlog_a0p = 0;
Yexpt=0; % exponential term for Y in p(Y|A0,A+) at some point such as the peak
Apexpt=0.0; % 0.0 because we have chosen posterior estimate of A+ as A+*. Exponential term for A+ conditional on A0 and Y
%======= Computing the log prior pdf of a0a+ and the exponential term for Y in p(Y|A0,A+).
for k=1:nvar
a0k = A0xhat(:,k); % meaningful parameters in the kth equation.
apk = Apxhat(:,k); % meaningful parameters in the kth equation.
%--- Prior settings.
S0bar = H0invtld{k}; %See Claim 2 on p.69b.
Spbar = Hpinvtld{k};
bk = Uiconst{k}'*a0k; % free parameters in the kth equation.
gk = Viconst{k}'*apk; % free parameters in the kth equation.
gbark = Ptld{k}*bk; % bar: prior
%--- The exponential term for Y in p(Y|A0,A+)
Yexpt = Yexpt - 0.5*(a0k'*YatYa*a0k - 2*apk'*XatYa*a0k + apk'*XatXa*apk);
%--- The log prior pdf.
vlog_a0p = vlog_a0p - 0.5*(size(Uiconst{k},2)+size(Viconst{k},2))*log(2*pi) + 0.5*log(abs(det(S0bar))) + ...
0.5*log(abs(det(Spbar))) - 0.5*(bk'*S0bar*bk+(gk-gbark)'*Spbar*(gk-gbark));
%--- For p(A+|Y,a0) only.
tmpd = gk - Pmat{k}*bk;
Apexpt = Apexpt - 0.5*tmpd'*(Hpinv{k}*tmpd);
end
vlog_a0p
%--------- The log value of p(Y|A0,A+) at some point such as the peak. ----------
%--------- Note that logMarLHres is the same as vlog_Y_a, just to double check. ----------
vlog_Y_a = -0.5*nvar*fss*log(2*pi) + fss*log(abs(det(A0xhat))) + Yexpt
% a: given a0 and a+
logMarLHres = 0; % Initialize log of the marginal likelihood (restricted or constant parameters).
for ki=1:fss %ndobs+1:fss % Forward recursion to get the marginal likelihood. See F on p.19 and pp. 48-49.
%---- Restricted log marginal likelihood function (constant parameters).
[A0l,A0u] = lu(A0xhat);
ada = sum(log(abs(diag(A0u)))); % log|A0|
termexp = y(ki,:)*A0xhat - phi(ki,:)*Apxhat; % 1-by-nvar
logMarLHres = logMarLHres - (0.5*nvar)*log(2*pi) + ada - 0.5*termexp*termexp'; % log value
end
logMarLHres
%--------- The log value of p(A+|Y,A0) at some point such as the peak ----------
totparsp = 0.0;
tmpd = 0.0;
for k=1:nvar
totparsp = totparsp + size(Viconst{k},2);
tmpd = tmpd + 0.5*log(abs(det(Hpinv{k})));
end
vlog_ap_Ya0 = -0.5*totparsp*log(2*pi) + tmpd + Apexpt;
%===================================
% Compute p(a0,k|Y,ao) at some point such as the peak (in this situation, we simply
% generate results from the original Gibbs sampler). See FORECAST (2) pp.70-71
%===================================
%--- Global set up for Gibbs.
[Tinv,UT] = fn_gibbsrvar_setup(H0inv, Uiconst, Hpinv, Pmat, Viconst, nvar, fss);
%
vlog_a0_Yao = zeros(nvar,1);
% the log value of p(a0k|Y,ao) where ao: other a's at some point such as the peak of ONLY some a0's
vlog=zeros(ndraws2,1);
for k=1:nvar
bk = Uiconst{k}'*A0xhat(:,k);
indx_ks=[k:nvar]; % the columns that exclude 1-(k-1)th columns
A0gbs0 = A0hat; % starting at some point such as the peak
nk = n0(k);
if k<nvar
%--------- The 1st set of draws to be tossed away. ------------------
for draws = 1:ndraws1
if ~mod(draws,nbuffer)
disp(' ')
disp(sprintf('The %dth column or equation in A0 with %d 1st tossed-away draws in Gibbs',k,draws))
end
A0gbs1 = fn_gibbsrvar(A0gbs0,UT,nvar,fss,n0,indx_ks);
A0gbs0=A0gbs1; % repeat the Gibbs sampling
end
%--------- The 2nd set of draws to be used. ------------------
for draws = 1:ndraws2
if ~mod(draws,nbuffer)
disp(' ')
disp(sprintf('The %dth column or equation in A0 with %d usable draws in Gibbs',k,draws))
end
[A0gbs1, Wcell] = fn_gibbsrvar(A0gbs0,UT,nvar,fss,n0,indx_ks);
%------ See p.71, Forecast (II).
%------ Computing p(a0_k|Y,a_others) at some point such as the peak along the dimensions of indx_ks.
Vk = Tinv{k}\Wcell{k}; %V_k on p.71 of Forecast (II).
gbeta = Vk\bk; % inv(V_k)*b_k on p.71 of Forecast (II) where alpha_k = b_k in our notation.
[Vtq,Vtr]=qr(Vk',0); %To get inv(V_k)'*inv(V_k) in (*) on p.71 of Forecast (II).
%
vlog(draws) = 0.5*(fss+nk)*log(fss)-log(abs(det(Vk)))-0.5*(nk-1)*log(2*pi)-...
0.5*(fss+1)*log(2)-gammaln(0.5*(fss+1))+fss*log(abs(gbeta(1)))-...
0.5*fss*bk'*(Vtr\(Vtr'\bk));
A0gbs0=A0gbs1; % repeat the Gibbs sampling
end
vlogm=max(vlog);
qlog=vlog-vlogm;
vlogxhat=vlogm-log(ndraws2)+log(sum(exp(qlog)));
vlog_a0_Yao(k) = vlogxhat;
% The log value of p(a0_k|Y,a_others) where a_others: other a's at some point such as the peak of ONLY some a0's
else
disp(' ')
disp(sprintf('The last(6th) column or equation in A0 with no Gibbs draws'))
[A0gbs1, Wcell] = fn_gibbsrvar(A0gbs0,UT,nvar,fss,n0,indx_ks)
%------ See p.71, Forecast (II).
%------ Computing p(a0_k|Y,a_others) at some point such as the peak along the dimensions of indx_ks.
Vk = Tinv{k}\Wcell{k}; %V_k on p.71 of Forecast (II).
gbeta = Vk\bk; % inv(V_k)*b_k on p.71 of Forecast (II) where alpha_k = b_k in our notation.
[Vtq,Vtr]=qr(Vk',0); %To get inv(V_k)'*inv(V_k) in (*) on p.71 of Forecast (II).
%
vloglast = 0.5*(fss+nk)*log(fss)-log(abs(det(Vk)))-0.5*(nk-1)*log(2*pi)-...
0.5*(fss+1)*log(2)-gammaln(0.5*(fss+1))+fss*log(abs(gbeta(1)))-...
0.5*fss*bk'*(Vtr\(Vtr'\bk));
vlog_a0_Yao(k) = vloglast;
end
end
ndraws2
disp('Prior pdf -- log(p(a0hat, a+hat)):');
vlog_a0p
disp('LH pdf -- log(p(Y|a0hat, a+hat)):');
vlog_Y_a
disp('Posterior Kernal -- logp(ahat) + logp(Y|ahat):');
vlog_Y_a + vlog_a0p
disp('Posterior pdf -- log(p(a0_i_hat|a0_other_hat, Y)):');
vlog_a0_Yao
disp('Posterior pdf -- log(p(aphat|a0hat, Y)):');
vlog_ap_Ya0
%--------- The value of marginal density p(Y) ----------
disp(' ');
disp(' ');
disp('************ Marginal Likelihood of Y or Marginal Data Density: ************');
vlogY = vlog_a0p+vlog_Y_a-sum(vlog_a0_Yao)-vlog_ap_Ya0
function ms_mardd(options_)
% Applies to both linear and exclusion restrictions.
% (1) Marginal likelihood function p(Y) for constant structural VAR models, using Chib (1995)'s ``Marginal Likelihood from the Gibbs Output'' in JASA.
% (2) Conditional likelihood function f(Y|A0, A+) on the ML estimate for constant exclusion-identified models.
% See Forecast (II) pp.67-80.
%
% Tao Zha, September 1999. Quick revisions, May 2003. Final revision, September 2004.
% Copyright (C) 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/>.
msstart2 % start the program in which everyhting is initialized through msstart2.m
if ~options_.ms.indxestima
warning(' ')
disp('You must set IxEstima=1 in msstart to run this program')
disp('Press ctrl-c to abort now')
pause
end
A0xhat = zeros(size(A0hat));
Apxhat = zeros(size(Aphat));
if (0)
%Robustness check to see if the same result is obtained with the purterbation of the parameters.
for k=1:nvar
bk = Uiconst{k}'*A0hat(:,k);
gk = Viconst{k}'*Aphat(:,k);
A0xhat(:,k) = Uiconst{k}*(bk + 5.2*randn(size(bk))); % Perturbing the posterior estimate.
Apxhat(:,k) = Viconst{k}*(gk + 5.2*randn(size(gk))); % Perturbing the posterior estimate.
end
else
%At the posterior estimate.
A0xhat = A0hat; % ML estimate of A0
Apxhat = Aphat; % ML estimate of A+
end
%--- Rename variables.
YatYa = yty;
XatYa = xty;
ytx = xty';
YatXa = ytx;
XatXa = xtx;
%--------- The log value of p(A0,A+) at some point such as the peak ----------
vlog_a0p = 0;
Yexpt=0; % exponential term for Y in p(Y|A0,A+) at some point such as the peak
Apexpt=0.0; % 0.0 because we have chosen posterior estimate of A+ as A+*. Exponential term for A+ conditional on A0 and Y
%======= Computing the log prior pdf of a0a+ and the exponential term for Y in p(Y|A0,A+).
for k=1:nvar
a0k = A0xhat(:,k); % meaningful parameters in the kth equation.
apk = Apxhat(:,k); % meaningful parameters in the kth equation.
%--- Prior settings.
S0bar = H0invtld{k}; %See Claim 2 on p.69b.
Spbar = Hpinvtld{k};
bk = Uiconst{k}'*a0k; % free parameters in the kth equation.
gk = Viconst{k}'*apk; % free parameters in the kth equation.
gbark = Ptld{k}*bk; % bar: prior
%--- The exponential term for Y in p(Y|A0,A+)
Yexpt = Yexpt - 0.5*(a0k'*YatYa*a0k - 2*apk'*XatYa*a0k + apk'*XatXa*apk);
%--- The log prior pdf.
vlog_a0p = vlog_a0p - 0.5*(size(Uiconst{k},2)+size(Viconst{k},2))*log(2*pi) + 0.5*log(abs(det(S0bar))) + ...
0.5*log(abs(det(Spbar))) - 0.5*(bk'*S0bar*bk+(gk-gbark)'*Spbar*(gk-gbark));
%--- For p(A+|Y,a0) only.
tmpd = gk - Pmat{k}*bk;
Apexpt = Apexpt - 0.5*tmpd'*(Hpinv{k}*tmpd);
end
vlog_a0p
%--------- The log value of p(Y|A0,A+) at some point such as the peak. ----------
%--------- Note that logMarLHres is the same as vlog_Y_a, just to double check. ----------
vlog_Y_a = -0.5*nvar*fss*log(2*pi) + fss*log(abs(det(A0xhat))) + Yexpt
% a: given a0 and a+
logMarLHres = 0; % Initialize log of the marginal likelihood (restricted or constant parameters).
for ki=1:fss %ndobs+1:fss % Forward recursion to get the marginal likelihood. See F on p.19 and pp. 48-49.
%---- Restricted log marginal likelihood function (constant parameters).
[A0l,A0u] = lu(A0xhat);
ada = sum(log(abs(diag(A0u)))); % log|A0|
termexp = y(ki,:)*A0xhat - phi(ki,:)*Apxhat; % 1-by-nvar
logMarLHres = logMarLHres - (0.5*nvar)*log(2*pi) + ada - 0.5*termexp*termexp'; % log value
end
logMarLHres
%--------- The log value of p(A+|Y,A0) at some point such as the peak ----------
totparsp = 0.0;
tmpd = 0.0;
for k=1:nvar
totparsp = totparsp + size(Viconst{k},2);
tmpd = tmpd + 0.5*log(abs(det(Hpinv{k})));
end
vlog_ap_Ya0 = -0.5*totparsp*log(2*pi) + tmpd + Apexpt;
%===================================
% Compute p(a0,k|Y,ao) at some point such as the peak (in this situation, we simply
% generate results from the original Gibbs sampler). See FORECAST (2) pp.70-71
%===================================
%--- Global set up for Gibbs.
[Tinv,UT] = fn_gibbsrvar_setup(H0inv, Uiconst, Hpinv, Pmat, Viconst, nvar, fss);
%
vlog_a0_Yao = zeros(nvar,1);
% the log value of p(a0k|Y,ao) where ao: other a's at some point such as the peak of ONLY some a0's
vlog=zeros(ndraws2,1);
for k=1:nvar
bk = Uiconst{k}'*A0xhat(:,k);
indx_ks=[k:nvar]; % the columns that exclude 1-(k-1)th columns
A0gbs0 = A0hat; % starting at some point such as the peak
nk = n0(k);
if k<nvar
%--------- The 1st set of draws to be tossed away. ------------------
for draws = 1:ndraws1
if ~mod(draws,nbuffer)
disp(' ')
disp(sprintf('The %dth column or equation in A0 with %d 1st tossed-away draws in Gibbs',k,draws))
end
A0gbs1 = fn_gibbsrvar(A0gbs0,UT,nvar,fss,n0,indx_ks);
A0gbs0=A0gbs1; % repeat the Gibbs sampling
end
%--------- The 2nd set of draws to be used. ------------------
for draws = 1:ndraws2
if ~mod(draws,nbuffer)
disp(' ')
disp(sprintf('The %dth column or equation in A0 with %d usable draws in Gibbs',k,draws))
end
[A0gbs1, Wcell] = fn_gibbsrvar(A0gbs0,UT,nvar,fss,n0,indx_ks);
%------ See p.71, Forecast (II).
%------ Computing p(a0_k|Y,a_others) at some point such as the peak along the dimensions of indx_ks.
Vk = Tinv{k}\Wcell{k}; %V_k on p.71 of Forecast (II).
gbeta = Vk\bk; % inv(V_k)*b_k on p.71 of Forecast (II) where alpha_k = b_k in our notation.
[Vtq,Vtr]=qr(Vk',0); %To get inv(V_k)'*inv(V_k) in (*) on p.71 of Forecast (II).
%
vlog(draws) = 0.5*(fss+nk)*log(fss)-log(abs(det(Vk)))-0.5*(nk-1)*log(2*pi)-...
0.5*(fss+1)*log(2)-gammaln(0.5*(fss+1))+fss*log(abs(gbeta(1)))-...
0.5*fss*bk'*(Vtr\(Vtr'\bk));
A0gbs0=A0gbs1; % repeat the Gibbs sampling
end
vlogm=max(vlog);
qlog=vlog-vlogm;
vlogxhat=vlogm-log(ndraws2)+log(sum(exp(qlog)));
vlog_a0_Yao(k) = vlogxhat;
% The log value of p(a0_k|Y,a_others) where a_others: other a's at some point such as the peak of ONLY some a0's
else
disp(' ')
disp(sprintf('The last(6th) column or equation in A0 with no Gibbs draws'))
[A0gbs1, Wcell] = fn_gibbsrvar(A0gbs0,UT,nvar,fss,n0,indx_ks)
%------ See p.71, Forecast (II).
%------ Computing p(a0_k|Y,a_others) at some point such as the peak along the dimensions of indx_ks.
Vk = Tinv{k}\Wcell{k}; %V_k on p.71 of Forecast (II).
gbeta = Vk\bk; % inv(V_k)*b_k on p.71 of Forecast (II) where alpha_k = b_k in our notation.
[Vtq,Vtr]=qr(Vk',0); %To get inv(V_k)'*inv(V_k) in (*) on p.71 of Forecast (II).
%
vloglast = 0.5*(fss+nk)*log(fss)-log(abs(det(Vk)))-0.5*(nk-1)*log(2*pi)-...
0.5*(fss+1)*log(2)-gammaln(0.5*(fss+1))+fss*log(abs(gbeta(1)))-...
0.5*fss*bk'*(Vtr\(Vtr'\bk));
vlog_a0_Yao(k) = vloglast;
end
end
ndraws2
disp('Prior pdf -- log(p(a0hat, a+hat)):');
vlog_a0p
disp('LH pdf -- log(p(Y|a0hat, a+hat)):');
vlog_Y_a
disp('Posterior Kernal -- logp(ahat) + logp(Y|ahat):');
vlog_Y_a + vlog_a0p
disp('Posterior pdf -- log(p(a0_i_hat|a0_other_hat, Y)):');
vlog_a0_Yao
disp('Posterior pdf -- log(p(aphat|a0hat, Y)):');
vlog_ap_Ya0
%--------- The value of marginal density p(Y) ----------
disp(' ');
disp(' ');
disp('************ Marginal Likelihood of Y or Marginal Data Density: ************');
vlogY = vlog_a0p+vlog_Y_a-sum(vlog_a0_Yao)-vlog_ap_Ya0

115
matlab/ms-sbvar/ms_variance_decomposition.m
View File

@ -34,46 +34,97 @@ function [options_, oo_]=ms_variance_decomposition(M_, options_, oo_)
disp('MS-SBVAR Variance Decomposition');
options_ = set_file_tags(options_);
[options_, oo_] = set_ms_estimation_file(options_.ms.file_tag, options_, oo_);
options_ = set_ms_simulation_file(options_);
clean_files_for_second_type_of_mex(M_, options_, 'variance_decomposition')
clean_ms_variance_decomposition_files(options_.ms.output_file_tag);
vddir = [options_.ms.output_file_tag filesep 'Variance_Decomposition'];
create_dir(vddir);
% NOTICE THAT VARIANCE DECOMPOSITION DEFAULTS TO USING THE MEAN, NOT MEDIAN OR BANDED
% setup command line options
opt = ['-variance_decomposition -seed ' num2str(options_.DynareRandomStreams.seed)];
opt = [opt ' -do ' vddir];
opt = [opt ' -ft ' options_.ms.file_tag];
opt = [opt ' -fto ' options_.ms.output_file_tag];
opt = [opt ' -horizon ' num2str(options_.ms.horizon)];
opt = [opt ' -thin ' num2str(options_.ms.thinning_factor)];
opt = {
{'file_tag', options_.ms.file_tag}, ...
{'seed', options_.DynareRandomStreams.seed}, ...
{'horizon', options_.ms.horizon}, ...
{'filtered', options_.ms.filtered_probabilities}, ...
{'error_bands', options_.ms.error_bands}, ...
{'percentiles', options_.ms.percentiles}, ...
{'thin', options_.ms.thinning_factor}, ...
{'mean'} ...
};
if options_.ms.median
opt = [opt(:)' {{'median'}}];
if options_.ms.regimes
opt = [opt ' -regimes'];
elseif options_.ms.regime
% regime-1 since regime is 0-indexed in C but 1-indexed in Matlab
opt = [opt ' -regime ' num2str(options_.ms.regime-1)];
elseif options_.ms.filtered_probabilities
opt = [opt ' -filtered'];
end
[err, vd] = mex_ms_variance_decomposition([opt(:)', {{'free_parameters',oo_.ms.maxparams}, ...
{'shocks_per_parameter', options_.ms.shock_draws}}]);
mexErrCheck('mex_ms_variance_decomposition ergodic ', err);
plot_ms_variance_decomposition(M_,options_,vd, 'Ergodic Variance Decomposition',options_.graph_save_formats,options_.TeX);
if options_.ms.parameter_uncertainty
options_ = set_ms_simulation_file(options_);
opt = [opt ' -parameter_uncertainty'];
opt = [opt ' -shocks_per_parameter ' num2str(options_.ms.shocks_per_parameter)];
else
opt = [opt ' -shocks_per_parameter ' num2str(options_.ms.shock_draws)];
end
[err, regime_vd] = mex_ms_variance_decomposition([opt(:)', {{'free_parameters',oo_.ms.maxparams}, ...
{'shocks_per_parameter', options_.ms.shock_draws}, {'regimes'}}]);
mexErrCheck('mex_ms_variance_decomposition ergodic regimes', err);
save([vddir filesep 'ergodic_vd.mat'], 'vd', 'regime_vd');
percentiles_size = 1;
outfile = [vddir filesep 'var_decomp_mean_'];
if options_.ms.error_bands
% error_bands / percentiles used differently by
% Dan's variance decomposition code
% no_error_bands => mean is computed
percentiles_size = size(options_.ms.percentiles,2);
opt = [opt ' -percentiles ' num2str(percentiles_size)];
for i=1:size(options_.ms.percentiles,2)
opt = [opt ' ' num2str(options_.ms.percentiles(i))];
end
outfile = [vddir filesep 'var_decomp_percentiles_'];
end
if exist(options_.ms.mh_file,'file') > 0
[err, vd] = mex_ms_variance_decomposition([opt(:)', {{'simulation_file',options_.ms.mh_file}, ...
{'shocks_per_parameter', options_.ms.shocks_per_parameter}, {'parameter_uncertainty'}}]);
mexErrCheck('mex_ms_variance_decomposition bayesian ', err);
% variance_decomposition
[err] = ms_sbvar_command_line(opt);
mexErrCheck('ms_variance_decomposition',err);
[err, regime_vd] = mex_ms_variance_decomposition([opt(:)', {{'simulation_file',options_.ms.mh_file}, ...
{'shocks_per_parameter', options_.ms.shocks_per_parameter}, {'parameter_uncertainty'}, {'regimes'}}]);
mexErrCheck('mex_ms_variance_decomposition bayesian regimes ', err);
save([vddir filesep 'bayesian_vd.mat'], 'vd', 'regime_vd');
if options_.ms.regime || options_.ms.regimes
outfile = [outfile 'regime_'];
if options_.ms.regime
outfile = [outfile num2str(options_.ms.regime-1) ...
'_' options_.ms.output_file_tag '.out'];
end
elseif options_.ms.filtered_probabilities
outfile = [outfile 'filtered_' options_.ms.output_file_tag '.out'];
else
outfile = [outfile 'ergodic_' options_.ms.output_file_tag '.out'];
end
% Create plots
if options_.ms.regimes
n_chains = length(options_.ms.ms_chain);
n_regimes=1;
for i_chain=1:n_chains
n_regimes = n_regimes*length(options_.ms.ms_chain(i_chain).regime);
end
for regime_i=1:n_regimes
vd_title = ['Variance Decomposition, Regime ' num2str(regime_i)];
vd_data = load([outfile num2str(regime_i-1) '_' ...
options_.ms.output_file_tag '.out'], '-ascii');
vd_data = reshape_ascii_variance_decomposition_data( ...
M_.endo_nbr, percentiles_size, options_.ms.horizon, vd_data);
save([vddir filesep 'variance_decomposition_regime_' num2str(regime_i-1)], 'vd_data');
plot_ms_variance_decomposition(M_, options_, vd_data, vd_title);
end
else
if options_.ms.regime
vd_title = ['Variance Decomposition, Regime ' num2str(options_.ms.regime)];
save_filename = ['variance_decomposition_regime_' num2str(options_.ms.regime-1)];
else
save_filename = 'variance_decomposition';
if options_.ms.filtered_probabilities
vd_title = 'Variance Decomposition Filtered';
else
vd_title = 'Variance Decomposition Ergodic';
end
end
vd_data = load(outfile, '-ascii');
vd_data = reshape_ascii_variance_decomposition_data( ...
M_.endo_nbr, percentiles_size, options_.ms.horizon, vd_data);
save([vddir filesep save_filename], 'vd_data');
plot_ms_variance_decomposition(M_, options_, vd_data, vd_title);
end
end

324
matlab/ms-sbvar/ms_write_markov_file.m
View File

@ -1,162 +1,162 @@
function ms_write_markov_file(fname, options)
% function ms_write_markov_file(fname, options)
%
% INPUTS
% fname: (string) name of markov file
% options: (struct) options
%
% OUTPUTS
%
% SPECIAL REQUIREMENTS
% none
% Copyright (C) 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/>.
n_chains = length(options.ms.ms_chain);
nvars = size(options.varobs,1);
fh = fopen(fname,'w');
%/******************************************************************************/
%/********************* Markov State Variable Information **********************/
%/******************************************************************************/
fprintf(fh,'//== Markov State Variables with Simple Restrictions ==//\n\n');
%//This number is NOT used but read in.
fprintf(fh,'//== Number Observations ==//\n');
fprintf(fh,'0\n\n');
fprintf(fh,'//== Number independent state variables ==//\n');
fprintf(fh,'%d\n\n',n_chains);
for i_chain = 1:n_chains
%//=====================================================//
%//== state_variable[i] (1 <= i <= n_state_variables) ==//
%//=====================================================//
fprintf(fh,'//== Number of states for state_variable[%d] ==//\n', ...
i_chain);
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.
%//== Each column contains the parameters for a Dirichlet prior on the corresponding
%//== column of the transition matrix. Each element must be positive. For each column,
%//== the relative size of the prior elements determine the relative size of the elements
%//== of the transition matrix and overall larger sizes implies a tighter prior.
fprintf(fh,['//== Transition matrix prior for state_variable[%d] ==//\n'], ...
i_chain);
Alpha = ones(n_states,n_states);
for i_state = 1:n_states
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');
end
fprintf(fh,['\n//== Dirichlet dimensions for state_variable[%d] ' ...
'==//\n'],i_chain);
% fprintf(fh,'%d ',repmat(n_states,1,n_states));
fprintf(fh,'%d ',repmat(2,1,n_states));
fprintf(fh,'\n\n');
%//== The jth restriction matrix is n_states-by-free[j]. Each row of the restriction
%//== matrix has exactly one non-zero entry and the sum of each column must be one.
fprintf(fh,['//== Column restrictions for state_variable[%d] ' ...
'==//\n'],i_chain);
for i_state = 1:n_states
if i_state == 1
M = eye(n_states,2);
elseif i_state == n_states
M = [zeros(n_states-2,2); eye(2)];
else
M = zeros(n_states,2);
M(i_state+[-1 1],1) = ones(2,1)/2;
M(i_state,2) = 1;
disp(M)
end
for j_state = 1:n_states
fprintf(fh,'%f ',M(j_state,:));
fprintf(fh,'\n');
end
fprintf(fh,'\n');
end
end
%/******************************************************************************/
%/******************************* VAR Parameters *******************************/
%/******************************************************************************/
%//NOT read
fprintf(fh,'//== Number Variables ==//\n');
fprintf(fh,'%d\n\n',nvars);
%//NOT read
fprintf(fh,'//== Number Lags ==//\n');
fprintf(fh,'%d\n\n',options.ms.nlags);
%//NOT read
fprintf(fh,'//== Exogenous Variables ==//\n');
fprintf(fh,'1\n\n');
%//== nvar x n_state_variables matrix. In the jth row, a non-zero value implies that
%this state variable controls the jth column of A0 and Aplus
fprintf(fh,['//== Controlling states variables for coefficients ==//\' ...
'n']);
for i_var = 1:nvars
for i_chain = 1:n_chains
if ~isfield(options.ms.ms_chain(i_chain),'svar_coefficients') ...
|| isempty(options.ms.ms_chain(i_chain).svar_coefficients)
i_equations = 0;
else
i_equations = ...
options.ms.ms_chain(i_chain).svar_coefficients.equations;
end
if strcmp(i_equations,'ALL') || any(i_equations == i_var)
fprintf(fh,'%d ',1);
else
fprintf(fh,'%d ',0);
end
end
fprintf(fh,'\n');
end
%//== nvar x n_state_variables matrix. In the jth row, a non-zero value implies that
%this state variable controls the jth diagonal element of Xi
fprintf(fh,'\n//== Controlling states variables for variance ==//\n');
for i_var = 1:nvars
for i_chain = 1:n_chains
if ~isfield(options.ms.ms_chain(i_chain),'svar_variances') ...
|| isempty(options.ms.ms_chain(i_chain).svar_variances)
i_equations = 0;
else
i_equations = ...
options.ms.ms_chain(i_chain).svar_variances.equations;
end
if strcmp(i_equations,'ALL') || any(i_equations == i_var)
fprintf(fh,'%d ',1);
else
fprintf(fh,'%d ',0);
end
end
fprintf(fh,'\n');
end
fclose(fh);
function ms_write_markov_file(fname, options)
% function ms_write_markov_file(fname, options)
%
% INPUTS
% fname: (string) name of markov file
% options: (struct) options
%
% OUTPUTS
%
% SPECIAL REQUIREMENTS
% none
% Copyright (C) 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/>.
n_chains = length(options.ms.ms_chain);
nvars = size(options.varobs,1);
fh = fopen(fname,'w');
%/******************************************************************************/
%/********************* Markov State Variable Information **********************/
%/******************************************************************************/
fprintf(fh,'//== Markov State Variables with Simple Restrictions ==//\n\n');
%//This number is NOT used but read in.
fprintf(fh,'//== Number Observations ==//\n');
fprintf(fh,'0\n\n');
fprintf(fh,'//== Number independent state variables ==//\n');
fprintf(fh,'%d\n\n',n_chains);
for i_chain = 1:n_chains
%//=====================================================//
%//== state_variable[i] (1 <= i <= n_state_variables) ==//
%//=====================================================//
fprintf(fh,'//== Number of states for state_variable[%d] ==//\n', ...
i_chain);
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.
%//== Each column contains the parameters for a Dirichlet prior on the corresponding
%//== column of the transition matrix. Each element must be positive. For each column,
%//== the relative size of the prior elements determine the relative size of the elements
%//== of the transition matrix and overall larger sizes implies a tighter prior.
fprintf(fh,['//== Transition matrix prior for state_variable[%d] ==//\n'], ...
i_chain);
Alpha = ones(n_states,n_states);
for i_state = 1:n_states
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');
end
fprintf(fh,['\n//== Dirichlet dimensions for state_variable[%d] ' ...
'==//\n'],i_chain);
% fprintf(fh,'%d ',repmat(n_states,1,n_states));
fprintf(fh,'%d ',repmat(2,1,n_states));
fprintf(fh,'\n\n');
%//== The jth restriction matrix is n_states-by-free[j]. Each row of the restriction
%//== matrix has exactly one non-zero entry and the sum of each column must be one.
fprintf(fh,['//== Column restrictions for state_variable[%d] ' ...
'==//\n'],i_chain);
for i_state = 1:n_states
if i_state == 1
M = eye(n_states,2);
elseif i_state == n_states
M = [zeros(n_states-2,2); eye(2)];
else
M = zeros(n_states,2);
M(i_state+[-1 1],1) = ones(2,1)/2;
M(i_state,2) = 1;
disp(M)
end
for j_state = 1:n_states
fprintf(fh,'%f ',M(j_state,:));
fprintf(fh,'\n');
end
fprintf(fh,'\n');
end
end
%/******************************************************************************/
%/******************************* VAR Parameters *******************************/
%/******************************************************************************/
%//NOT read
fprintf(fh,'//== Number Variables ==//\n');
fprintf(fh,'%d\n\n',nvars);
%//NOT read
fprintf(fh,'//== Number Lags ==//\n');
fprintf(fh,'%d\n\n',options.ms.nlags);
%//NOT read
fprintf(fh,'//== Exogenous Variables ==//\n');
fprintf(fh,'1\n\n');
%//== nvar x n_state_variables matrix. In the jth row, a non-zero value implies that
%this state variable controls the jth column of A0 and Aplus
fprintf(fh,['//== Controlling states variables for coefficients ==//\' ...
'n']);
for i_var = 1:nvars
for i_chain = 1:n_chains
if ~isfield(options.ms.ms_chain(i_chain),'svar_coefficients') ...
|| isempty(options.ms.ms_chain(i_chain).svar_coefficients)
i_equations = 0;
else
i_equations = ...
options.ms.ms_chain(i_chain).svar_coefficients.equations;
end
if strcmp(i_equations,'ALL') || any(i_equations == i_var)
fprintf(fh,'%d ',1);
else
fprintf(fh,'%d ',0);
end
end
fprintf(fh,'\n');
end
%//== nvar x n_state_variables matrix. In the jth row, a non-zero value implies that
%this state variable controls the jth diagonal element of Xi
fprintf(fh,'\n//== Controlling states variables for variance ==//\n');
for i_var = 1:nvars
for i_chain = 1:n_chains
if ~isfield(options.ms.ms_chain(i_chain),'svar_variances') ...
|| isempty(options.ms.ms_chain(i_chain).svar_variances)
i_equations = 0;
else
i_equations = ...
options.ms.ms_chain(i_chain).svar_variances.equations;
end
if strcmp(i_equations,'ALL') || any(i_equations == i_var)
fprintf(fh,'%d ',1);
else
fprintf(fh,'%d ',0);
end
end
fprintf(fh,'\n');
end
fclose(fh);

134
matlab/ms-sbvar/ms_write_mhm_input.m
View File

@ -1,68 +1,68 @@
function ms_write_mhm_input(fname, options_ms)
% function ms_write_mhm_input(fname, options_ms)
%
% INPUTS
% fname: (string) filename
% options_ms: (struct) options
%
% OUTPUTS
% none
%
% SPECIAL REQUIREMENTS
% none
% Copyright (C) 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/>.
fh = fopen(fname,'w');
fprintf(fh,'/**********************************************************\n');
fprintf(fh,' *** This input file is read by swzmsbvar_mhm_1 and swzmsbvar_mhm_1.exe only, NOT by swzmsbvar_printdraws.exe.\n');
fprintf(fh,' ***\n');
fprintf(fh,' **********************************************************/\n');
fprintf(fh,'\n\n//------------- 1st set of posterior draws to find optimal scales for Metropolis (30000). ---------------\n');
fprintf(fh,'//== number draws for first burn-in ==// //For determining the Metropolis scales only.\n');
fprintf(fh,'%d\n\n',options_ms.draws_nbr_burn_in_1);
fprintf(fh,'//------------- MCMC burn-in draws once the Metropolis scales (previous stage) are fixed. --------------\n');
fprintf(fh,'//------------- 2nd set of standard burn-in posterior draws to throw away the initial draws (10000). ---------------\n');
fprintf(fh,'//== number draws for second burn-in ==//\n');
fprintf(fh,'%d\n\n',options_ms.draws_nbr_burn_in_2);
fprintf(fh,'//--------------- 1st set of posterior draws to compute the mean and variance for the weighting function in the MHM (200000) ----------------\n');
fprintf(fh,'//== number draws to estimate mean and variance ==//\n');
fprintf(fh,'%d\n\n',options_ms.draws_nbr_mean_var_estimate);
fprintf(fh,'//--------------- Only applied to mhm_2 process: total number of MCMC draws = thinning factor * 2nd set of saved posterior draws ----------------\n');
fprintf(fh,'//== thinning factor for modified harmonic mean process ==//\n');
fprintf(fh,'%d\n\n',options_ms.thinning_factor);
fprintf(fh,'//--------------- 2nd set of saved posterior draws from MHM_2 (second stage): saved draws AFTER thinning (1000000) ----------------\n');
fprintf(fh,'//== number draws for modified harmonic mean process ==//\n');
fprintf(fh,'%d\n\n',options_ms.draws_nbr_modified_harmonic_mean);
fprintf(fh,'//------- 1st stage: computing all three tightness factors for Dirichlet. ---------\n');
fprintf(fh,'//------- 2nd stage: hard-code the second scale factor (in principle, we can do all three). ---------\n');
fprintf(fh,'//------- It seems that Dan''s code only use the first element of the following scales. The scale applies to the Dirichlet''s hyperparameter alpha for the diagonal of the transition matrix in the weighting function. Note that the weighting function for the transition matrix parameters is Dirichlet. ---------\n');
fprintf(fh,'//== scale values for Dirichlet distribution ==//\n');
fprintf(fh,'3\n\n');
fprintf(fh,'%f ',options_ms.dirichlet_scale);
fprintf(fh,'\n');
function ms_write_mhm_input(fname, options_ms)
% function ms_write_mhm_input(fname, options_ms)
%
% INPUTS
% fname: (string) filename
% options_ms: (struct) options
%
% OUTPUTS
% none
%
% SPECIAL REQUIREMENTS
% none
% Copyright (C) 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/>.
fh = fopen(fname,'w');
fprintf(fh,'/**********************************************************\n');
fprintf(fh,' *** This input file is read by swzmsbvar_mhm_1 and swzmsbvar_mhm_1.exe only, NOT by swzmsbvar_printdraws.exe.\n');
fprintf(fh,' ***\n');
fprintf(fh,' **********************************************************/\n');
fprintf(fh,'\n\n//------------- 1st set of posterior draws to find optimal scales for Metropolis (30000). ---------------\n');
fprintf(fh,'//== number draws for first burn-in ==// //For determining the Metropolis scales only.\n');
fprintf(fh,'%d\n\n',options_ms.draws_nbr_burn_in_1);
fprintf(fh,'//------------- MCMC burn-in draws once the Metropolis scales (previous stage) are fixed. --------------\n');
fprintf(fh,'//------------- 2nd set of standard burn-in posterior draws to throw away the initial draws (10000). ---------------\n');
fprintf(fh,'//== number draws for second burn-in ==//\n');
fprintf(fh,'%d\n\n',options_ms.draws_nbr_burn_in_2);
fprintf(fh,'//--------------- 1st set of posterior draws to compute the mean and variance for the weighting function in the MHM (200000) ----------------\n');
fprintf(fh,'//== number draws to estimate mean and variance ==//\n');
fprintf(fh,'%d\n\n',options_ms.draws_nbr_mean_var_estimate);
fprintf(fh,'//--------------- Only applied to mhm_2 process: total number of MCMC draws = thinning factor * 2nd set of saved posterior draws ----------------\n');
fprintf(fh,'//== thinning factor for modified harmonic mean process ==//\n');
fprintf(fh,'%d\n\n',options_ms.thinning_factor);
fprintf(fh,'//--------------- 2nd set of saved posterior draws from MHM_2 (second stage): saved draws AFTER thinning (1000000) ----------------\n');
fprintf(fh,'//== number draws for modified harmonic mean process ==//\n');
fprintf(fh,'%d\n\n',options_ms.draws_nbr_modified_harmonic_mean);
fprintf(fh,'//------- 1st stage: computing all three tightness factors for Dirichlet. ---------\n');
fprintf(fh,'//------- 2nd stage: hard-code the second scale factor (in principle, we can do all three). ---------\n');
fprintf(fh,'//------- It seems that Dan''s code only use the first element of the following scales. The scale applies to the Dirichlet''s hyperparameter alpha for the diagonal of the transition matrix in the weighting function. Note that the weighting function for the transition matrix parameters is Dirichlet. ---------\n');
fprintf(fh,'//== scale values for Dirichlet distribution ==//\n');
fprintf(fh,'3\n\n');
fprintf(fh,'%f ',options_ms.dirichlet_scale);
fprintf(fh,'\n');
fclose(fh);

1619
matlab/ms-sbvar/msstart2.m
View File

File diff suppressed because it is too large Load Diff

331
matlab/ms-sbvar/msstart_setup.m
View File

@ -1,166 +1,165 @@
%function []= msstart_setup(options_)
% Copyright (C) 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/>.
% ** ONLY UNDER UNIX SYSTEM
%path(path,'/usr2/f1taz14/mymatlab')
%===========================================
% Exordium I
%===========================================
format short g % format
%
%options_.ms.freq = 4; % quarters or months
%options_.ms.initial_year=1959; % beginning of the year
%options_.ms.initial_subperiod=1; % begining of the quarter or month
%options_.ms.final_year=2005; % final year
%options_.ms.final_subperiod=4; % final month or quarter
nData=(options_.ms.final_year-options_.ms.initial_year)*options_.ms.freq + (options_.ms.final_subperiod-options_.ms.initial_subperiod+1);
% total number of the available data -- this is all you have
%*** Load data and series
%load datainf_argen.prn % the default name for the variable is "options_.ms.data".
%load datacbogdpffr.prn
%options_.ms.data = datacbogdpffr;
%clear datacbogdpffr;
[nt,ndv]=size(options_.data);
if nt~=nData
disp(' ')
warning(sprintf('nt=%d, Caution: not equal to the length in the data',nt));
%disp(sprintf('nt=%d, Caution: not equal to the length in the data',nt));
disp('Press ctrl-c to abort')
return
end
%--------
%1 CBO output gap -- log(x_t)-log(x_t potential)
%2 GDP deflator -- (P_t/P_{t-1})^4-1.0
%2 FFR/100.
options_.ms.vlist = [1:size(options_.varobs,1)]; % 1: U; 4: PCE inflation.
options_.ms.varlist=cellstr(options_.varobs);
options_.ms.log_var = sort(varlist_indices(options_.ms.vlistlog,options_.varobs)); % subset of "options_.ms.vlist. Variables in log level so that differences are in **monthly** growth, unlike R and U which are in annual percent (divided by 100 already).
options_.ms.percent_var =setdiff(options_.ms.vlist,options_.ms.log_var);
%options_.ms.restriction_fname='ftd_upperchol3v'; %Only used by msstart2.m.
ylab = options_.ms.varlist;
xlab = options_.ms.varlist;
%----------------
nvar = size(options_.varobs,1); % number of endogenous variables
nlogeno = length(options_.ms.log_var); % number of endogenous variables in options_.ms.log_var
npereno = length(options_.ms.percent_var); % number of endogenous variables in options_.ms.percent_var
if (nvar~=(nlogeno+npereno))
disp(' ')
warning('Check xlab, nlogeno or npereno to make sure of endogenous variables in options_.ms.vlist')
disp('Press ctrl-c to abort')
return
elseif (nvar==length(options_.ms.vlist))
nexo=1; % only constants as an exogenous variable. The default setting.
elseif (nvar<length(options_.ms.vlist))
nexo=length(options_.ms.vlist)-nvar+1;
else
disp(' ')
warning('Make sure there are only nvar endogenous variables in options_.ms.vlist')
disp('Press ctrl-c to abort')
return
end
%------- A specific sample is considered for estimation -------
yrStart=options_.ms.initial_year;
qmStart=options_.ms.initial_subperiod;
yrEnd=options_.ms.final_year;
qmEnd=options_.ms.final_subperiod;
%options_.forecast = 4; % number of years for forecasting
if options_.forecast<1
error('To be safe, the number of forecast years should be at least 1')
end
ystr=num2str(yrEnd);
forelabel = [ ystr(3:4) ':' num2str(qmEnd) ' Forecast'];
nSample=(yrEnd-yrStart)*options_.ms.freq + (qmEnd-qmStart+1);
if qmEnd==options_.ms.freq
E1yrqm = [yrEnd+1 1]; % first year and quarter (month) after the sample
else
E1yrqm = [yrEnd qmEnd+1]; % first year and quarter (month) after the sample
end
E2yrqm = [yrEnd+options_.forecast qmEnd]; % end at the last month (quarter) of a calendar year after the sample
[fdates,nfqm]=fn_calyrqm(options_.ms.freq,E1yrqm,E2yrqm); % forecast dates and number of forecast dates
[sdates,nsqm] = fn_calyrqm(options_.ms.freq,[yrStart qmStart],[yrEnd qmEnd]);
% sdates: dates for the whole sample (including options_.ms.nlags)
if nSample~=nsqm
warning('Make sure that nSample is consistent with the size of sdates')
disp('Hit any key to continue, or ctrl-c to abort')
pause
end
imstp = 4*options_.ms.freq; % <<>> impulse responses (4 years)
nayr = 4; %options_.forecast; % number of years before forecasting for plotting.
%------- Prior, etc. -------
%options_.ms.nlags = 4; % number of options_.ms.nlags
%options_.ms.cross_restrictions = 0; % 1: cross-A0-and-A+ restrictions; 0: options_.ms.restriction_fname is all we have
% Example for indxOres==1: restrictions of the form P(t) = P(t-1).
%options_.ms.contemp_reduced_form = 0; % 1: contemporaneous recursive reduced form; 0: restricted (non-recursive) form
%options_.ms.real_pseudo_forecast = 0; % 1: options_.ms.real_pseudo_forecast forecasts; 0: real time forecasts
%options_.ms.bayesian_prior = 1; % 1: Bayesian prior; 0: no prior
indxDummy = options_.ms.bayesian_prior; % 1: add dummy observations to the data; 0: no dummy added.
%options_.ms.dummy_obs = 0; % No dummy observations for xtx, phi, fss, xdatae, etc. Dummy observations are used as an explicit prior in fn_rnrprior_covres_dobs.m.
%if indxDummy
% options_.ms.dummy_obs=nvar+1; % number of dummy observations
%else
% options_.ms.dummy_obs=0; % no dummy observations
%end
%=== The following mu is effective only if options_.ms.bayesian_prior==1.
mu = options_.ms.coefficients_prior_hyperparameters;
% mu(1): overall tightness and also for A0;
% mu(2): relative tightness for A+;
% mu(3): relative tightness for the constant term;
% mu(4): tightness on lag decay; (1)
% mu(5): weight on nvar sums of coeffs dummy observations (unit roots);
% mu(6): weight on single dummy initial observation including constant
% (cointegration, unit roots, and stationarity);
%
%
hpmsmd = [0.0; 0.0];
indxmsmdeqn = [0; 0; 0; 0]; %This option disenable using this in fn_rnrprior_covres_dobs.m
tdf = 3; % degrees of freedom for t-dist for initial draw of the MC loop
nbuffer = 1000; % a block or buffer of draws (buffer) that is saved to the disk (not memory)
ndraws1=1*nbuffer; % 1st part of Monte Carlo draws
ndraws2=10*ndraws1; % 2nd part of Monte Carlo draws
% seednumber = options_.DynareRandomStreams.seed; %7910; %472534; % if 0, random state at each clock time
% % good one 420 for [29 45], [29 54]
% if seednumber
% randn('state',seednumber);
% rand('state',seednumber);
% else
% randn('state',fix(100*sum(clock)));
% rand('state',fix(100*sum(clock)));
% end
% nstarts=1 % number of starting points
% imndraws = nstarts*ndraws2; % total draws for impulse responses or forecasts
%<<<<<<<<<<<<<<<<<<<
%function []= msstart_setup(options_)
% Copyright (C) 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/>.
% ** ONLY UNDER UNIX SYSTEM
%path(path,'/usr2/f1taz14/mymatlab')
%===========================================
% Exordium I
%===========================================
format short g % format
%
%options_.ms.freq = 4; % quarters or months
%options_.ms.initial_year=1959; % beginning of the year
%options_.ms.initial_subperiod=1; % begining of the quarter or month
%options_.ms.final_year=2005; % final year
%options_.ms.final_subperiod=4; % final month or quarter
nData=(options_.ms.final_year-options_.ms.initial_year)*options_.ms.freq + (options_.ms.final_subperiod-options_.ms.initial_subperiod+1);
% total number of the available data -- this is all you have
%*** Load data and series
%load datainf_argen.prn % the default name for the variable is "options_.ms.data".
%load datacbogdpffr.prn
%options_.ms.data = datacbogdpffr;
%clear datacbogdpffr;
[nt,ndv]=size(options_.data);
if nt~=nData
disp(' ')
warning(sprintf('nt=%d, Caution: not equal to the length in the data',nt));
%disp(sprintf('nt=%d, Caution: not equal to the length in the data',nt));
disp('Press ctrl-c to abort')
return
end
%--------
%1 CBO output gap -- log(x_t)-log(x_t potential)
%2 GDP deflator -- (P_t/P_{t-1})^4-1.0
%2 FFR/100.
options_.ms.vlist = [1:size(options_.varobs,1)]; % 1: U; 4: PCE inflation.
options_.ms.varlist=cellstr(options_.varobs);
options_.ms.log_var = sort(varlist_indices(options_.ms.vlistlog,options_.varobs)); % subset of "options_.ms.vlist. Variables in log level so that differences are in **monthly** growth, unlike R and U which are in annual percent (divided by 100 already).
options_.ms.percent_var =setdiff(options_.ms.vlist,options_.ms.log_var);
%options_.ms.restriction_fname='ftd_upperchol3v'; %Only used by msstart2.m.
ylab = options_.ms.varlist;
xlab = options_.ms.varlist;
%----------------
nvar = size(options_.varobs,1); % number of endogenous variables
nlogeno = length(options_.ms.log_var); % number of endogenous variables in options_.ms.log_var
npereno = length(options_.ms.percent_var); % number of endogenous variables in options_.ms.percent_var
if (nvar~=(nlogeno+npereno))
disp(' ')
warning('Check xlab, nlogeno or npereno to make sure of endogenous variables in options_.ms.vlist')
disp('Press ctrl-c to abort')
return
elseif (nvar==length(options_.ms.vlist))
nexo=1; % only constants as an exogenous variable. The default setting.
elseif (nvar<length(options_.ms.vlist))
nexo=length(options_.ms.vlist)-nvar+1;
else
disp(' ')
warning('Make sure there are only nvar endogenous variables in options_.ms.vlist')
disp('Press ctrl-c to abort')
return
end
%------- A specific sample is considered for estimation -------
yrStart=options_.ms.initial_year;
qmStart=options_.ms.initial_subperiod;
yrEnd=options_.ms.final_year;
qmEnd=options_.ms.final_subperiod;
%options_.forecast = 4; % number of years for forecasting
if options_.forecast<1
error('To be safe, the number of forecast years should be at least 1')
end
forelabel = [num2str(yrEnd) ':' num2str(qmEnd) ' Forecast'];
nSample=(yrEnd-yrStart)*options_.ms.freq + (qmEnd-qmStart+1);
if qmEnd==options_.ms.freq
E1yrqm = [yrEnd+1 1]; % first year and quarter (month) after the sample
else
E1yrqm = [yrEnd qmEnd+1]; % first year and quarter (month) after the sample
end
E2yrqm = [yrEnd+options_.forecast qmEnd]; % end at the last month (quarter) of a calendar year after the sample
[fdates,nfqm]=fn_calyrqm(options_.ms.freq,E1yrqm,E2yrqm); % forecast dates and number of forecast dates
[sdates,nsqm] = fn_calyrqm(options_.ms.freq,[yrStart qmStart],[yrEnd qmEnd]);
% sdates: dates for the whole sample (including options_.ms.nlags)
if nSample~=nsqm
warning('Make sure that nSample is consistent with the size of sdates')
disp('Hit any key to continue, or ctrl-c to abort')
pause
end
imstp = 4*options_.ms.freq; % <<>> impulse responses (4 years)
nayr = 4; %options_.forecast; % number of years before forecasting for plotting.
%------- Prior, etc. -------
%options_.ms.nlags = 4; % number of options_.ms.nlags
%options_.ms.cross_restrictions = 0; % 1: cross-A0-and-A+ restrictions; 0: options_.ms.restriction_fname is all we have
% Example for indxOres==1: restrictions of the form P(t) = P(t-1).
%options_.ms.contemp_reduced_form = 0; % 1: contemporaneous recursive reduced form; 0: restricted (non-recursive) form
%options_.ms.real_pseudo_forecast = 0; % 1: options_.ms.real_pseudo_forecast forecasts; 0: real time forecasts
%options_.ms.bayesian_prior = 1; % 1: Bayesian prior; 0: no prior
indxDummy = options_.ms.bayesian_prior; % 1: add dummy observations to the data; 0: no dummy added.
%options_.ms.dummy_obs = 0; % No dummy observations for xtx, phi, fss, xdatae, etc. Dummy observations are used as an explicit prior in fn_rnrprior_covres_dobs.m.
%if indxDummy
% options_.ms.dummy_obs=nvar+1; % number of dummy observations
%else
% options_.ms.dummy_obs=0; % no dummy observations
%end
%=== The following mu is effective only if options_.ms.bayesian_prior==1.
mu = options_.ms.coefficients_prior_hyperparameters;
% mu(1): overall tightness and also for A0;
% mu(2): relative tightness for A+;
% mu(3): relative tightness for the constant term;
% mu(4): tightness on lag decay; (1)
% mu(5): weight on nvar sums of coeffs dummy observations (unit roots);
% mu(6): weight on single dummy initial observation including constant
% (cointegration, unit roots, and stationarity);
%
%
hpmsmd = [0.0; 0.0];
indxmsmdeqn = [0; 0; 0; 0]; %This option disenable using this in fn_rnrprior_covres_dobs.m
tdf = 3; % degrees of freedom for t-dist for initial draw of the MC loop
nbuffer = 1000; % a block or buffer of draws (buffer) that is saved to the disk (not memory)
ndraws1=1*nbuffer; % 1st part of Monte Carlo draws
ndraws2=10*ndraws1; % 2nd part of Monte Carlo draws
% seednumber = options_.DynareRandomStreams.seed; %7910; %472534; % if 0, random state at each clock time
% % good one 420 for [29 45], [29 54]
% if seednumber
% randn('state',seednumber);
% rand('state',seednumber);
% else
% randn('state',fix(100*sum(clock)));
% rand('state',fix(100*sum(clock)));
% end
% nstarts=1 % number of starting points
% imndraws = nstarts*ndraws2; % total draws for impulse responses or forecasts
%<<<<<<<<<<<<<<<<<<<

71
matlab/ms-sbvar/plot_ms_forecast.m
View File

@ -1,16 +1,21 @@
function plot_ms_forecast(M_,options_,forecast,title_,save_graph_formats,TeX)
% function plot_ms_forecast(M_,options_,forecast,title_,save_graph_formats,TeX)
function plot_ms_forecast(M_, options_, forecast, figure_name)
% function plot_ms_forecast(M_, options_, forecast, figure_name)
% plots the forecast from the output from a ms-sbvar
%
% INPUTS
% M_
% forecast should be in the form (percentile x horizon x nvar ), if banded otherwise
% ( horizon x nvar )
% M_: (struct) model structure
% options_: (struct) options
% forecast: (matrix) in the form (percentile x horizon x nvar ), if banded otherwise
% ( horizon x nvar )
% figure_name: (string) title
%
% title: optional super title
% OUTPUTS
% none
%
% SPECIAL REQUIREMENTS
% none
% Copyright (C) 2011 Dynare Team
% Copyright (C) 2011-2012 Dynare Team
%
% This file is part of Dynare.
%
@ -26,15 +31,12 @@ function plot_ms_forecast(M_,options_,forecast,title_,save_graph_formats,TeX)
%
% You should have received a copy of the GNU General Public License
% along with Dynare. If not, see <http://www.gnu.org/licenses/>.
nc = 2;
nr = 2;
nvars = M_.endo_nbr;
endo_names = M_.endo_names;
var_list = endo_names(1:M_.orig_endo_nbr,:);
i_var = [];
names = {};
tex_names = {};
m = 1;
@ -49,16 +51,9 @@ function plot_ms_forecast(M_,options_,forecast,title_,save_graph_formats,TeX)
tex_names{m} = tex_name;
m = m + 1;
end
i_var = [i_var; tmp];
end
nvar = length(i_var);
dims = size(forecast);
if nargin < 3
title_ = '';
end
if (length(dims) == 2)
% Point Forecast (horizon x nvars )
horizon = dims(1);
@ -66,40 +61,40 @@ function plot_ms_forecast(M_,options_,forecast,title_,save_graph_formats,TeX)
elseif (length(dims) == 3)
% Banded Forecast
horizon = dims(2);
num_percentiles = dims(1);
num_percentiles = dims(1);
else
error('The impulse response matrix passed to be plotted does not appear to be the correct size');
end
if num_percentiles == 1
plot_point_forecast(forecast, nvars, nr, nc, var_list, title_, ...
save_graph_formats, TeX, names, tex_names, ...
plot_point_forecast(forecast, nvars, nr, nc, var_list, figure_name, ...
options_.graph_save_formats, options_.TeX, names, tex_names, ...
[options_.ms.output_file_tag filesep 'Output' filesep 'Forecast']);
else
plot_banded_forecast(forecast, nvars, nr, nc, var_list, num_percentiles, ...
title_, save_graph_formats, TeX, names, tex_names, ...
figure_name, options_.graph_save_formats, options_.TeX, names, tex_names, ...
[options_.ms.output_file_tag filesep 'Output' filesep 'Forecast']);
end
end
function plot_point_forecast(forecast,nvars,nr,nc,endo_names,title_,save_graph_formats,TeX,names,tex_names,dirname)
function plot_point_forecast(forecast,nvars,nr,nc,endo_names,figure_name,save_graph_formats,TeX,names,tex_names,dirname)
if nvars > nr*nc
graph_name = 'MS-Forecast (1)';
fig = figure('Name','Forecast (I)');
graph_name = ['MS (1) ' figure_name];
figure('Name', graph_name);
else
graph_name = 'MS-Forecast';
fig = figure('Name','Forecast');
end
graph_name = figure_name;
figure('Name', graph_name);
end
m = 1;
n_fig = 1;
for j=1:nvars
if m > nr*nc
graph_name = ['MS-Forecast (' int2str(n_fig) ')']
graph_name = ['MS (' int2str(n_fig) ') ' figure_name];
dyn_save_graph(dirname,['MS-forecast-' int2str(n_fig)],...
save_graph_formats,TeX,names,tex_names,graph_name);
n_fig =n_fig+1;
figure('Name',['MS-Forecast (' int2str(n_fig) ')']);
figure('Name', graph_name);
m = 1;
end
subplot(nr,nc,m);
@ -115,19 +110,19 @@ function plot_point_forecast(forecast,nvars,nr,nc,endo_names,title_,save_graph_f
end
end
function plot_banded_forecast(forecast,nvars,nr,nc,endo_names,num_percentiles,title_,save_graph_formats,TeX,names,tex_names,dirname)
function plot_banded_forecast(forecast,nvars,nr,nc,endo_names,num_percentiles,figure_name,save_graph_formats,TeX,names,tex_names,dirname)
if nvars > nr*nc
graph_name = 'MS-Forecast (1)';
fig = figure('Name','Forecast (I)');
graph_name = ['MS (1) ' figure_name];
figure('Name', graph_name);
else
graph_name = 'MS-Forecast';
fig = figure('Name','Forecast');
end
graph_name = figure_name;
figure('Name', graph_name);
end
m = 1;
n_fig = 1;
for j=1:nvars
if m > nr*nc
graph_name = ['MS-Forecast (' int2str(n_fig) ')'];
graph_name = ['MS (' int2str(n_fig) ') ' figure_name];
dyn_save_graph(dirname,['MS-forecast-' int2str(n_fig)],...
save_graph_formats,TeX,names,tex_names,graph_name);
n_fig =n_fig+1;

68
matlab/ms-sbvar/plot_ms_irf.m
View File

@ -1,21 +1,21 @@
function plot_ms_irf(M_, options_, irf, names, title_, varlist)
% function plot_ms_irf(M_, options_, irf, names, title_, varlist)
function plot_ms_irf(M_, options_, irf, figure_name, varlist)
% function plot_ms_irf(M_, options_, irf, figure_name, varlist)
% plots the impulse responses from the output from a ms-sbvar
%
% INPUTS
% M_
% irf should be in the form (percentile x horizon x (nvar x nvar)), if banded otherwise
% ( horizon x (nvar x nvar) )
% M_: (struct) model structure
% options_: (struct) options
% irf: (matrix) in the form (percentile x horizon x (nvar x nvar)), if banded otherwise
% ( horizon x (nvar x nvar) )
% figure_name: (string) title
%
% names: character list of the names of the variables
% OUTPUTS
% none
%
% title: optional super title
%
% The element in position (k,i+j*nvars) of ir is the response of the ith
% variable to the jth shock at horizon k. Horizon 0 is the contemporaneous
% response.
% SPECIAL REQUIREMENTS
% none
% Copyright (C) 2011 Dynare Team
% Copyright (C) 2011-2012 Dynare Team
%
% This file is part of Dynare.
%
@ -33,19 +33,16 @@ function plot_ms_irf(M_, options_, irf, names, title_, varlist)
% along with Dynare. If not, see <http://www.gnu.org/licenses/>.
if nargin < 4
title_ = '';
figure_name = '';
end
nc = 2;
nr = 2;
nvars = M_.endo_nbr;
endo_names = M_.endo_names;
if isempty(varlist)
var_list = endo_names(1:M_.orig_endo_nbr,:);
end
i_var = [];
names = {};
tex_names = {};
m = 1;
@ -54,7 +51,6 @@ function plot_ms_irf(M_, options_, irf, names, title_, varlist)
if isempty(tmp)
error([var_list(i,:) ' isn''t and endogenous variable'])
end
i_var = [i_var; tmp];
tex_name = deblank(M_.endo_names_tex(tmp,:));
if ~isempty(tex_name)
names{m} = deblank(var_list(i,:));
@ -71,11 +67,8 @@ function plot_ms_irf(M_, options_, irf, names, title_, varlist)
m = m + 1;
end
end
nvar = length(i_var);
dims = size(irf);
if (length(dims) == 2)
% Point IRF (horizon x (nvarsxnvars) )
horizon = dims(1);
@ -87,11 +80,11 @@ function plot_ms_irf(M_, options_, irf, names, title_, varlist)
else
error('The impulse response matrix passed to be plotted does not appear to be the correct size');
end
if size(endo_names,1) ~= nvars
error('The names passed are not the same length as the number of variables')
error('The names passed are not the same length as the number of variables');
end
if num_percentiles == 1
% loop through the shocks
for s=1:nvars
@ -100,7 +93,7 @@ function plot_ms_irf(M_, options_, irf, names, title_, varlist)
shock(:,i) = irf(:,((i-1) + ((s-1)*nvars)+1));
end
plot_point_irf_for_shock(shock, nvars,endo_names, deblank(endo_names(s,:)), ...
title_, [options_.ms.output_file_tag filesep 'Output' filesep 'IRF'], options_, names, tex_names);
figure_name, [options_.ms.output_file_tag filesep 'Output' filesep 'IRF'], options_, names, tex_names);
end
else
for s=1:nvars
@ -111,26 +104,25 @@ function plot_ms_irf(M_, options_, irf, names, title_, varlist)
end
end
plot_banded_irf_for_shock(shock, nvars,endo_names, deblank(endo_names(s,:)), ...
title_, [options_.ms.output_file_tag filesep 'Output' filesep 'IRF'], options_, names, tex_names);
figure_name, [options_.ms.output_file_tag filesep 'Output' filesep 'IRF'], options_, names, tex_names);
end
end
end
function [fig] = plot_point_irf_for_shock(irf,nvars,endo_names,shock_name,title_,dirname,options_,names,tex_names)
fig = figure('Name',title_);
function [fig] = plot_point_irf_for_shock(irf,nvars,endo_names,shock_name,figure_name,dirname,options_,names,tex_names)
fig = figure('Name',figure_name);
for k=1:nvars
subplot(ceil(sqrt(nvars)), ceil(sqrt(nvars)),k);
plot(irf(:,k))
disp([endo_names(k,:) ' shock from ' shock_name]);
title([endo_names(k,:) ' shock from ' shock_name]);
end
dyn_save_graph(dirname,[title_ ' ' shock_name],options_.graph_save_formats, ...
options_.TeX,names,tex_names,[title_ ' ' shock_name]);
dyn_save_graph(dirname,[figure_name ' ' shock_name],options_.graph_save_formats, ...
options_.TeX,names,tex_names,[figure_name ' ' shock_name]);
end
function [fig] = plot_banded_irf_for_shock(irf,nvars, endo_names, shock_name,title_,dirname,options_,names,tex_names)
fig = figure('Name',title_);
function [fig] = plot_banded_irf_for_shock(irf,nvars, endo_names, shock_name,figure_name,dirname,options_,names,tex_names)
fig = figure('Name',figure_name);
npercentiles = size(irf,3);
for k=1:nvars
subplot(ceil(sqrt(nvars)), ceil(sqrt(nvars)),k);
@ -142,8 +134,6 @@ function [fig] = plot_banded_irf_for_shock(irf,nvars, endo_names, shock_name,tit
disp([endo_names(k,:) ' shock from ' shock_name]);
title([endo_names(k,:) ' shock from ' shock_name]);
end
dyn_save_graph(dirname,[title_ ' ' shock_name],options_.graph_save_formats, ...
options_.TeX,names,tex_names,[title_ ' ' shock_name]);
dyn_save_graph(dirname,[figure_name ' ' shock_name],options_.graph_save_formats, ...
options_.TeX,names,tex_names,[figure_name ' ' shock_name]);
end

70
matlab/ms-sbvar/plot_ms_variance_decomposition.m
View File

@ -1,12 +1,14 @@
function plot_ms_variance_decomposition(M_, options_, vd, title_, graph_save_formats, TeX, varargin)
% function plot_ms_variance_decomposition(M_, options_, vd, title_, graph_save_formats, TeX, varargin)
function plot_ms_variance_decomposition(M_, options_, vd, figure_name, varargin)
% function plot_ms_variance_decomposition(M_, options_, vd, figure_name, varargin)
% plot the variance decomposition of shocks
%
% Inputs
% M_
% shocks: matrix of the individual shocks Tx(KxK)with J=number of shocks
% INPUTS
% M_: (struct) model structure
% options_: (struct) options
% vd: (matrix) variance decomposition
% figure_name: (string) graph name
%
% Optional Inputs
% OPTIONAL INPUTS
% 'data': the actual data, TxK with K=number of data series
% 'steady': the steady state value, TxK
% 'shock_names': to specify the names of the shocks
@ -14,10 +16,13 @@ function plot_ms_variance_decomposition(M_, options_, vd, title_, graph_save_for
% 'dates': pass a date vector to use, otherwise will just index on 1:T
% 'colors': Jx3 list of the rgb colors to use for each shock
%
% Example:
% plot_historic_decomposition(shocks,'VD','shock_names',shock_names,'series_names',series_names)
% OUTPUTS
% none
%
% SPECIAL REQUIREMENTS
% none
% Copyright (C) 2011 Dynare Team
% Copyright (C) 2011-2012 Dynare Team
%
% This file is part of Dynare.
%
@ -34,6 +39,11 @@ function plot_ms_variance_decomposition(M_, options_, vd, title_, graph_save_for
% You should have received a copy of the GNU General Public License
% along with Dynare. If not, see <http://www.gnu.org/licenses/>.
if length(size(vd)) == 3
plot_ms_variance_decomposition_error_bands(M_, options_, vd, figure_name);
return;
end
nvars = M_.endo_nbr;
endo_names = M_.endo_names;
@ -48,11 +58,13 @@ function plot_ms_variance_decomposition(M_, options_, vd, title_, graph_save_for
m = m + 1;
end
end
dims = size(vd);
if length(dims) == 3
[T,K,J] = dims;
shocks = vd;
T = dims(1);
K = dims(2);
J = dims(3);
shocks = vd;
else
T = dims(1);
K = nvars;
@ -65,30 +77,25 @@ function plot_ms_variance_decomposition(M_, options_, vd, title_, graph_save_for
end
shocks = temp_vd;
end
for i=1:nvars
shock_names{i} = endo_names(i,:);
series_names{i} = endo_names(i,:);
end
if nargin < 2
title_ = '';
end
x = [1:T]; plot_dates = 0;
x = [1:T];
plot_dates = 0;
data = 0;
steady = 0;
colors = [ .1 .1 .75
.8 0 0
1 .7 .25
1 1 0
.5 1 .5
1 1 0
.5 1 .5
.7 .7 .1
.5 .6 .2
.1 .5 .1];
% overide the defaults with optional inputs
for i=1:length(varargin)
if strcmpi(varargin{i},'data')
@ -105,12 +112,11 @@ function plot_ms_variance_decomposition(M_, options_, vd, title_, graph_save_for
colors = varargin{i+1};
end
end
% add an extra period to the time series
x(T+1) = x(T) + (x(T) - x(T-1));
figure('Name',title_)
figure('Name',figure_name)
for k=1:K
% Go through each series
subplot(K,1,k);
@ -144,13 +150,13 @@ function plot_ms_variance_decomposition(M_, options_, vd, title_, graph_save_for
plot(x(2:end)',steady(:,k), '--k','LineWidth',2.25);
end
if k==K
if exist('OCTAVE_VERSION')
if exist('OCTAVE_VERSION')
legend(shock_names,'Location','SouthOutside');
else
legend(shock_names,'Location','BestOutside','Orientation','horizontal');
else
legend(shock_names,'Location','BestOutside','Orientation','horizontal');
end
end
hold off
if plot_dates
datetick 'x';
@ -159,9 +165,9 @@ function plot_ms_variance_decomposition(M_, options_, vd, title_, graph_save_for
ylim([0 , 1])
grid on
title(series_names{k});
%suptitle(title_);
end
dyn_save_graph([options_.ms.output_file_tag filesep 'Output' ...
filesep 'Variance_Decomposition'], 'MS-Variance-Decomposition', ...
graph_save_formats, TeX,names,tex_names,'Variance decomposition');
options_.graph_save_formats, options_.TeX, names, tex_names, ...
'Variance decomposition');
end

105
matlab/ms-sbvar/plot_ms_variance_decomposition_error_bands.m Normal file
View File

@ -0,0 +1,105 @@
function plot_ms_variance_decomposition_error_bands(M_, options_, vddata, figure_name)
% function plot_ms_variance_decomposition_error_bands(M_, options_, vddata, figure_name)
% plots the variance decomposition with percentiles
%
% INPUTS
% M_: (struct) model structure
% options_: (struct) options
% vddata: (matrix) variance_decomposition (percentile, options_.ms.horizon, nvar
% x nvar)
% figure_name: (string) title
%
% OUTPUTS
% none
%
% SPECIAL REQUIREMENTS
% none
% Copyright (C) 2011-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/>.
nvars = M_.endo_nbr;
endo_names = M_.endo_names;
var_list = endo_names(1:M_.orig_endo_nbr,:);
names = {};
tex_names = {};
m = 1;
for i = 1:size(var_list)
tmp = strmatch(var_list(i,:), endo_names, 'exact');
if isempty(tmp)
error([var_list(i,:) ' isn''t an endogenous variable'])
end
tex_name = deblank(M_.endo_names_tex(tmp,:));
if ~isempty(tex_name)
names{m} = deblank(var_list(i,:));
tex_names{m} = tex_name;
m = m + 1;
end
end
for i=1:M_.exo_nbr
tex_name = deblank(M_.exo_names_tex(i,:));
if ~isempty(tex_name)
names{m} = deblank(M_.exo_names(i,:));
tex_names{m} = tex_name;
m = m + 1;
end
end
dims = size(vddata);
if length(dims) ~= 3
error('The variance decomposition matrix passed to be plotted does not appear to be the correct size');
end
num_percentiles = dims(1);
if size(endo_names, 1) ~= nvars
error('The names passed are not the same length as the number of variables')
end
for s=1:nvars
shock = zeros(options_.ms.horizon, nvars, num_percentiles);
for n=1:num_percentiles
for i=1:nvars
shock(:,i,n) = vddata(n, :, ((i-1) + ((s-1)*nvars)+1));
end
end
plot_banded_vddata_for_shock(shock, nvars, endo_names, ...
deblank(endo_names(s,:)), figure_name, ...
[options_.ms.output_file_tag filesep 'Output' filesep 'Variance_Decomposition'], ...
options_, names, tex_names);
end
end
function [fig] = plot_banded_vddata_for_shock(vddata, nvars, endo_names, ...
shock_name, figure_name, dirname, options_, names, tex_names)
fig = figure('Name', figure_name);
npercentiles = size(vddata,3);
for k=1:nvars
subplot(ceil(sqrt(nvars)), ceil(sqrt(nvars)),k);
for nn=1:npercentiles
plot(vddata(:,k,nn))
hold on
end
hold off
disp([endo_names(k,:) ' contribution to ' shock_name]);
title([endo_names(k,:) ' contribution to ' shock_name]);
end
dyn_save_graph(dirname, [figure_name ' ' shock_name], ...
options_.graph_save_formats, options_.TeX, names, tex_names, ...
[figure_name ' ' shock_name]);
end

44
matlab/ms-sbvar/reshape_ascii_forecast_data.m Normal file
View File

@ -0,0 +1,44 @@
function forecast_data=reshape_ascii_forecast_data(endo_nbr, psize, horizon, ascii_data)
% function forecast_data=reshape_ascii_forecast_data(endo_nbr, psize, horizon, ascii_data)
%
% INPUTS
% endo_nbr: number of endogenous
% psize: number of percentiles
% horizon: forecast horizon
% ascii_data: data from .out file created by Dan's C code
%
% OUTPUTS
% forecast_data: new 3-d array holding data with error bands
%
% SPECIAL REQUIREMENTS
% none
% Copyright (C) 2011-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/>.
if psize <= 1
forecast_data = ascii_data;
return
end
forecast_data = zeros(psize, horizon, endo_nbr);
for i=1:endo_nbr
for j=1:psize
forecast_data(j,:,i) = ascii_data(1+horizon*(j-1):horizon*j,i)';
end
end
end

44
matlab/ms-sbvar/reshape_ascii_irf_data.m Normal file
View File

@ -0,0 +1,44 @@
function irf_data=reshape_ascii_irf_data(endo_nbr, psize, horizon, ascii_data)
% function irf_data=reshape_ascii_irf_data(endo_nbr, psize, horizon, ascii_data)
%
% INPUTS
% endo_nbr: number of endogenous
% psize: number of percentiles
% horizon: forecast horizon
% ascii_data: data from .out file created by Dan's C code
%
% OUTPUTS
% irf_data: new 3-d array holding data with error bands
%
% SPECIAL REQUIREMENTS
% none
% Copyright (C) 2011-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/>.
if psize <= 1
irf_data = ascii_data;
return
end
irf_data = zeros(psize, horizon, endo_nbr*endo_nbr);
for i=1:endo_nbr*endo_nbr
for j=1:psize
irf_data(j,:,i) = ascii_data(1+horizon*(j-1):horizon*j,i)';
end
end
end

44
matlab/ms-sbvar/reshape_ascii_variance_decomposition_data.m Normal file
View File

@ -0,0 +1,44 @@
function vd_data=reshape_ascii_variance_decomposition_data(endo_nbr, psize, horizon, ascii_data)
% function vd_data=reshape_ascii_vd_data(endo_nbr, psize, horizon, ascii_data)
%
% INPUTS
% endo_nbr: number of endogenous
% psize: number of percentiles
% horizon: forecast horizon
% ascii_data: data from .out file created by Dan's C code
%
% OUTPUTS
% vd_data: new 3-d array holding data with error bands
%
% SPECIAL REQUIREMENTS
% none
% Copyright (C) 2011-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/>.
if psize <= 1
vd_data = ascii_data;
return
end
vd_data = zeros(psize, horizon, endo_nbr*endo_nbr);
for i=1:endo_nbr*endo_nbr
for j=1:psize
vd_data(j,:,i) = ascii_data(1+horizon*(j-1):horizon*j,i)';
end
end
end

80
matlab/non_linear_dsge_likelihood.m
View File

@ -101,7 +101,7 @@ function [fval,exit_flag,ys,trend_coeff,info,Model,DynareOptions,BayesInfo,Dynar
%! @end deftypefn
%@eod:
% Copyright (C) 2010-2011 Dynare Team
% Copyright (C) 2010, 2011, 2012 Dynare Team
%
% This file is part of Dynare.
%
@ -127,7 +127,6 @@ 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;
@ -142,10 +141,10 @@ end
fval = [];
ys = [];
trend_coeff = [];
cost_flag = 1;
exit_flag = 1;
% Set the number of observed variables
nvobs = DynareDataset.info.vobs;
nvobs = DynareDataset.info.nvobs;
%------------------------------------------------------------------------------
% 1. Get the structural parameters & define penalties
@ -155,7 +154,7 @@ nvobs = DynareDataset.info.vobs;
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;
exit_flag = 0;
info = 41;
return
end
@ -164,7 +163,7 @@ end
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;
exit_flag = 0;
info = 42;
return
end
@ -172,50 +171,50 @@ 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);
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);
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;
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);
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.
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;
exit_flag = 0;
info = 43;
return
end
end
offset = offset+EstimatedParameters_.ncx;
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));
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
@ -227,17 +226,17 @@ if EstimatedParameters_.ncn
k = find(a < 0);
if k > 0
fval = penalty+sum(-a(k));
cost_flag = 0;
exit_flag = 0;
info = 44;
return
end
end
offset = offset+EstimatedParameters_.ncn;
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);
if EstimatedParameters.np > 0
Model.params(EstimatedParameters.param_vals(:,1)) = xparam1(offset+1:end);
end
% Update Model.Sigma_e and Model.H.
@ -253,11 +252,11 @@ Model.H = H;
if info(1) == 1 || info(1) == 2 || info(1) == 5
fval = penalty+1;
cost_flag = 0;
exit_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;
exit_flag = 0;
return
end
@ -266,8 +265,8 @@ BayesInfo.mf = BayesInfo.mf1;
% Define the deterministic linear trend of the measurement equation.
if DynareOptions.noconstant
constant = zeros(nvobs,1);
else
constant = zeros(nvobs,1);
else
if DynareOptions.loglinear
constant = log(SteadyState(BayesInfo.mfys));
else
@ -275,7 +274,6 @@ else
end
end
% Define the deterministic linear trend of the measurement equation.
if BayesInfo.with_trend
trend_coeff = zeros(DynareDataset.info.nvobs,1);
@ -294,7 +292,7 @@ end
start = DynareOptions.presample+1;
np = size(T,1);
mf = BayesInfo.mf;
Y = transpose(dataset_.rawdata);
Y = transpose(DynareDataset.rawdata);
%------------------------------------------------------------------------------
% 3. Initial condition of the Kalman filter
@ -332,10 +330,10 @@ 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.
case 1% Initial state vector covariance 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).
case 2% Initial state vector covariance 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;
@ -344,7 +342,7 @@ switch DynareOptions.particle.initialization
StateVectorVariance = cov(y_');
DynareOptions.periods = old_DynareOptionsperiods;
clear('old_DynareOptionsperiods','y_');
case 3
case 3% Initial state vector covariance is a diagonal matrix.
StateVectorMean = ReducedForm.constant(mf0);
StateVectorVariance = DynareOptions.particle.initial_state_prior_std*eye(number_of_state_variables);
otherwise
@ -357,13 +355,13 @@ ReducedForm.StateVectorVariance = StateVectorVariance;
% 4. Likelihood evaluation
%------------------------------------------------------------------------------
DynareOptions.warning_for_steadystate = 0;
LIK = feval(DynareOptions.particle.algorithm,ReducedForm,Y,[]);
LIK = feval(DynareOptions.particle.algorithm,ReducedForm,Y,[],DynareOptions);
if imag(LIK)
likelihood = penalty;
cost_flag = 0;
exit_flag = 0;
elseif isnan(LIK)
likelihood = penalty;
cost_flag = 0;
exit_flag = 0;
else
likelihood = LIK;
end

3
matlab/numgrad2.m
View File

@ -32,6 +32,7 @@ g=zeros(n,1);
badg=0;
goog=1;
scale=1;
g0 = 0;
for i=1:n
if size(x,1)>size(x,2)
tvecv=tvec(i,:);
@ -57,4 +58,4 @@ for i=1:n
g(i) = 0;
badg = 1;
end
end
end

3
matlab/parallel/AnalyseComputationalEnvironment.m
View File

@ -559,6 +559,7 @@ for Node=1:length(DataInput) % To obtain a recoursive function remove the 'for'
si0=[];
de0=[];
disp('Checking Hardware please wait ...');
if (DataInput(Node).Local == 1)
if Environment,
[si0 de0]=system('grep processor /proc/cpuinfo');
@ -579,7 +580,7 @@ for Node=1:length(DataInput) % To obtain a recoursive function remove the 'for'
RealCPUnbr='';
keyboard;
% keyboard;
RealCPUnbr=GiveCPUnumber(de0,OStargetUnix);
% Questo controllo penso che si possa MIGLIORARE!!!!!

30
matlab/parallel/closeSlave.m
View File

@ -1,4 +1,4 @@
function closeSlave(Parallel,TmpFolder),
function closeSlave(Parallel,TmpFolder,partial),
% PARALLEL CONTEXT
% In parallel context, this utility closes all remote matlab instances
% called by masterParallel when strategy (1) is active i.e. always open (which leaves
@ -32,6 +32,32 @@ function closeSlave(Parallel,TmpFolder),
% You should have received a copy of the GNU General Public License
% along with Dynare. If not, see <http://www.gnu.org/licenses/>.
if nargin<3,
partial=0;
end
s=warning('off');
if partial==1
save('slaveParallel_break','partial')
for indPC=1:length(Parallel),
if (Parallel(indPC).Local==0),
dynareParallelSendFiles('slaveParallel_break.mat',TmpFolder,Parallel(indPC));
end
end
% delete('slaveParallel_break')
return
end
if partial==-1
delete('slaveParallel_break.mat')
for indPC=1:length(Parallel),
if (Parallel(indPC).Local==0),
dynareParallelDelete( 'slaveParallel_break.mat',TmpFolder,Parallel(indPC));
end
end
% delete('slaveParallel_break')
return
end
for indPC=1:length(Parallel),
if (Parallel(indPC).Local==0),
@ -58,3 +84,5 @@ while(1)
end
end
s=warning('on');

2
matlab/parallel/dynareParallelDelete.m
View File

@ -40,7 +40,7 @@ else
end
for indPC=1:length(Parallel),
if isunix
if ~ispc || strcmpi('unix',Parallel(indPC).OperatingSystem),
[NonServeS NonServeD]=system(['ssh ',Parallel(indPC).UserName,'@',Parallel(indPC).ComputerName,' rm -f ',Parallel(indPC).RemoteDirectory,'/',pname,fname]);
else
delete(['\\',Parallel(indPC).ComputerName,'\',Parallel(indPC).RemoteDrive,'$\',Parallel(indPC).RemoteDirectory,'\',pname,fname]);

6
matlab/parallel/dynareParallelDeleteNewFiles.m
View File

@ -57,9 +57,9 @@ for indPC=1:length(Parallel),
fileaddress={sT(1:SlashNumberAndPosition(end)),sT(SlashNumberAndPosition(end)+1:end)};
dynareParallelDelete(fileaddress{2},[PRCDir,fS,fileaddress{1}],Parallel(indPC));
display('New file deleted in remote -->');
display(fileaddress{2});
display('<--');
disp('New file deleted in remote -->');
disp(fileaddress{2});
disp('<--');
end
else

6
matlab/parallel/dynareParallelGetNewFiles.m
View File

@ -58,9 +58,9 @@ for indPC=1:length(Parallel),
fileaddress={sT(1:SlashNumberAndPosition(end)),sT(SlashNumberAndPosition(end)+1:end)};
dynareParallelGetFiles(fileaddress,PRCDir,Parallel(indPC));
display('New file copied in local -->');
display(fileaddress{2});
display('<--');
disp('New file copied in local -->');
disp(fileaddress{2});
disp('<--');
end
else

18
matlab/parallel/dynareParallelRmDir.m
View File

@ -35,7 +35,25 @@ if nargin ==0,
return
end
% security check of remote folder delete
ok(1)=isempty(strfind(Parallel_info.RemoteTmpFolder,'..'));
tmp1=strfind(Parallel_info.RemoteTmpFolder,'2');
ok(2)=tmp1(1)==1;
ok(3)=~isempty(strfind(Parallel_info.RemoteTmpFolder,'-'));
ok(4)=~isempty(strfind(Parallel_info.RemoteTmpFolder,'h'));
ok(5)=~isempty(strfind(Parallel_info.RemoteTmpFolder,'m'));
ok(6)=~isempty(strfind(Parallel_info.RemoteTmpFolder,'s'));
ok(7)=~isempty(PRCDir);
if sum(ok)<7,
error('The name of the remote tmp folder does not comply the security standards!'),
end
for indPC=1:length(Parallel),
ok(1)=isempty(strfind(Parallel(indPC).RemoteDirectory,'..'));
if sum(ok)<7,
error('The remote folder path structure does not comply the security standards!'),
end
while (1)
if ~ispc || strcmpi('unix',Parallel(indPC).OperatingSystem)
[stat NonServe] = system(['ssh ',Parallel(indPC).UserName,'@',Parallel(indPC).ComputerName,' rm -fr ',Parallel(indPC).RemoteDirectory,'/',PRCDir,]);

3
matlab/parallel/fMessageStatus.m
View File

@ -43,7 +43,8 @@ catch
end
fslave = dir( ['slaveParallel_input',int2str(njob),'.mat']);
if isempty(fslave),
fbreak = dir( ['slaveParallel_break.mat']);
if isempty(fslave) || ~isempty(fbreak),
error('Master asked to break the job');
end

31
matlab/parallel/fParallel.m
View File

@ -81,20 +81,33 @@ try,
% Save the output result.
save([ fname,'_output_',int2str(whoiam),'.mat'],'fOutputVar' )
end
if isfield(fOutputVar,'CloseAllSlaves'),
CloseAllSlaves = 1;
fOutputVar = rmfield(fOutputVar,'CloseAllSlaves');
save([ fname,'_output_',int2str(whoiam),'.mat'],'fOutputVar' )
save(['comp_status_',funcName,int2str(whoiam),'.mat'],'CloseAllSlaves');
end
disp(['fParallel ',int2str(whoiam),' completed.'])
catch,
disp(['fParallel ',int2str(whoiam),' crashed.'])
fOutputVar.error = lasterror;
save([ fname,'_output_',int2str(whoiam),'.mat'],'fOutputVar' )
waitbarString = fOutputVar.error.message;
% waitbarTitle=['Metropolis-Hastings ',options_.parallel(ThisMatlab).ComputerName];
if Parallel(ThisMatlab).Local,
waitbarTitle='Local ';
theerror = lasterror;
if strfind(theerror.message,'Master asked to break the job')
fOutputVar.message = theerror;
save([ fname,'_output_',int2str(whoiam),'.mat'],'fOutputVar' )
waitbarString = theerror.message;
else
waitbarTitle=[Parallel(ThisMatlab).ComputerName];
disp(['fParallel ',int2str(whoiam),' crashed.'])
fOutputVar.error = theerror;
save([ fname,'_output_',int2str(whoiam),'.mat'],'fOutputVar' )
waitbarString = theerror.message;
% waitbarTitle=['Metropolis-Hastings ',options_.parallel(ThisMatlab).ComputerName];
if Parallel(ThisMatlab).Local,
waitbarTitle='Local ';
else
waitbarTitle=[Parallel(ThisMatlab).ComputerName];
end
fMessageStatus(NaN,whoiam,waitbarString, waitbarTitle, Parallel(ThisMatlab));
end
fMessageStatus(NaN,whoiam,waitbarString, waitbarTitle, Parallel(ThisMatlab));
end
diary off;

Some files were not shown because too many files have changed in this diff Show More