Merge branch 'master' into ecb-master

time-shift
Stéphane Adjemian (Charybdis) 2017-06-03 14:51:28 +02:00
commit 7d373f0d7e
741 changed files with 8340 additions and 7707 deletions

4
.gitignore vendored
View File

@ -60,6 +60,7 @@ checksum
/doc/dynare.info
/doc/dynare.info-1
/doc/dynare.info-2
/doc/dynare.info-3
/doc/dynare.cp
/doc/dynare.fn
/doc/dynare.fns
@ -212,3 +213,6 @@ tests/julia/rbc/rbc*.jl
# Octave variables saved when Octave crashes
octave-workspace
# VERSION generated file
VERSION

805
NEWS
View File

@ -1,3 +1,800 @@
Announcement for Dynare 4.5.0 (on 2013-12-16)
=============================================
We are pleased to announce the release of Dynare 4.5.0.
This major release adds new features and fixes various bugs.
The Windows packages are already available for download at:
http://www.dynare.org/download/dynare-stable
The Mac and Debian/Ubuntu packages should follow soon.
All users are strongly encouraged to upgrade.
This release is compatible with MATLAB versions ranging from 7.3 (R2006b) to
9.2 (R2017a) and with GNU Octave version 4.2.
Here is the list of major user-visible changes:
Dynare 4.5
==========
- Ramsey policy
+ Added command `ramsey_model` that builds the expanded model with
FOC conditions for the planner's problem but doesn't perform any
computation. Usefull to compute Ramsey policy in a perfect
foresight model,
+ `ramsey_policy` accepts multipliers in its variable list and
displays results for them.
- Perfect foresight models
+ New commands `perfect_foresight_setup` (for preparing the
simulation) and `perfect_foresight_solver` (for computing it). The
old `simul` command still exist and is now an alias for
`perfect_foresight_setup` + `perfect_foresight_solver`. It is no
longer possible to manipulate by hand the contents of
`oo_.exo_simul` when using `simul`. People who want to do
it must first call `perfect_foresight_setup`, then do the
manipulations, then call `perfect_foresight_solver`,
+ By default, the perfect foresight solver will try a homotopy
method if it fails to converge at the first try. The old behavior
can be restored with the `no_homotopy` option,
+ New option `stack_solve_algo=7` that allows specifying a
`solve_algo` solver for solving the model,
+ New option `solve_algo` that allows specifying a solver for
solving the model when using `stack_solve_algo=7`,
+ New option `lmmcp` that solves the model via a Levenberg-Marquardt
mixed complementarity problem (LMMCP) solver,
+ New option `robust_lin_solve` that triggers the use of a robust
linear solver for the default `solve_algo=4`,
+ New options `tolf` and `tolx` to control termination criteria of
solvers,
+ New option `endogenous_terminal_period` to `simul`,
+ Added the possibility to set the initial condition of the
(stochastic) extended path simulations with the histval block.
- Optimal simple rules
+ Saves the optimal value of parameters to `oo_.osr.optim_params`,
+ New block `osr_params_bounds` allows specifying bounds for the
estimated parameters,
+ New option `opt_algo` allows selecting different optimizers while
the new option `optim` allows specifying the optimizer options,
+ The `osr` command now saves the names, bounds, and indices for the
estimated parameters as well as the indices and weights of the
variables entering the objective function into `M_.osr`.
- Forecasts and Smoothing
+ The smoother and forecasts take uncertainty about trends and means
into account,
+ Forecasts accounting for measurement error are now saved in fields
of the form `HPDinf_ME` and `HPDsup_ME`,
+ New fields `oo_.Smoother.Trend` and `oo_.Smoother.Constant` that
save the trend and constant parts of the smoothed variables,
+ new field `oo_.Smoother.TrendCoeffs` that stores the trend
coefficients.
+ Rolling window forecasts allowed in `estimation` command by
passing a vector to `first_obs`,
+ The `calib_smoother` command now accepts the `loglinear`,
`prefilter`, `first_obs` and `filter_decomposition` options.
- Estimation
+ New options: `logdata`, `consider_all_endogenous`,
`consider_only_observed`, `posterior_max_subsample_draws`,
`mh_conf_sig`, `diffuse_kalman_tol`, `dirname`, `nodecomposition`
+ `load_mh_file` and `mh_recover` now try to load chain's proposal density,
+ New option `load_results_after_load_mh` that allows loading some
posterior results from a previous run if no new MCMC draws are
added,
+ New option `posterior_nograph` that suppresses the generation of
graphs associated with Bayesian IRFs, posterior smoothed objects,
and posterior forecasts,
+ Saves the posterior density at the mode in
`oo_.posterior.optimization.log_density`,
+ The `filter_covariance` option now also works with posterior
sampling like Metropolis-Hastings,
+ New option `no_posterior_kernel_density` to suppress computation
of kernel density of posterior objects,
+ Recursive estimation and forecasting now provides the individual
`oo_` structures for each sample in `oo_recursive_`,
+ The `trace_plot` command can now plot the posterior density,
+ New command `generate_trace_plots` allows generating all trace
plots for one chain,
+ New commands `prior_function` and `posterior_function` that
execute a user-defined function on parameter draws from the
prior/posterior distribution,
+ New option `huge_number` for replacement of infinite bounds with
large number during `mode_compute`,
+ New option `posterior_sampling_method` allows selecting the new
posterior sampling options:
`tailored_random_block_metropolis_hastings` (Tailored randomized
block (TaRB) Metropolis-Hastings), `slice` (Slice sampler),
`independent_metropolis_hastings` (Independent
Metropolis-Hastings),
+ New option `posterior_sampler_options` that allow controlling the
options of the `posterior_sampling_method`, its `scale_file`-option
pair allows loading the `_mh_scale.mat`-file storing the tuned
scale factor from a previous run of `mode_compute=6`,
+ New option `raftery_lewis_diagnostics` that computes Raftery/Lewis
(1992) convergence diagnostics,
+ New option `fast_kalman_filter` that provides fast Kalman filter
using Chandrasekhar recursions as described in Ed Herbst (2015),
+ The `dsge_var` option now saves results at the posterior mode into
`oo_.dsge_var`,
+ New option `smoothed_state_uncertainty` to provide the uncertainty
estimate for the smoothed state estimate from the Kalman smoother,
+ New prior density: generalized Weibull distribution,
+ Option `mh_recover` now allows continuing a crashed chain at the
last save mh-file,
+ New option `nonlinear_filter_initialization` for the
`estimation` command. Controls the initial covariance matrix
of the state variables in nonlinear filters.
+ The `conditional_variance_decomposition` option now displays
output and stores it as a LaTeX-table when the `TeX` option is
invoked,
+ The `use_calibration` to `estimated_params_init` now also works
with ML,
+ Improved initial estimation checks.
- Steady state
+ The default solver for finding the steady state is now a
trust-region solver (can be triggered explicitly with option
`solve_algo=4`),
+ New options `tolf` and `tolx` to control termination criteria of
solver,
+ The debugging mode now provides the termination values in steady
state finding.
- Stochastic simulations
+ New options `nodecomposition`,
+ New option `bandpass_filter` to compute bandpass-filtered
theoretical and simulated moments,
+ New option `one_sided_hp_filter` to compute one-sided HP-filtered
simulated moments,
+ `stoch_simul` displays a simulated variance decomposition when
simulated moments are requested,
+ `stoch_simul` saves skewness and kurtosis into respective fields
of `oo_` when simulated moments have been requested,
+ `stoch_simul` saves the unconditional variance decomposition in
`oo_.variance_decomposition`,
+ New option `dr_display_tol` that governs omission of small terms
in display of decision rules,
+ The `stoch_simul` command now prints the displayed tables as LaTeX
code when the new `TeX` option is enabled,
+ The `loglinear` option now works with lagged and leaded exogenous
variables like news shocks,
+ New option `spectral_density` that allows displaying the spectral
density of (filtered) endogenous variables,
+ New option `contemporaneous_correlation` that allows saving
contemporaneous correlations in addition to the covariances.
- Identification
+ New options `diffuse_filter` and `prior_trunc`,
+ The `identification` command now supports correlations via
simulated moments,
- Sensitivity analysis
+ New blocks `irf_calibration` and `moment_calibration`,
+ Outputs LaTeX tables if the new `TeX` option is used,
+ New option `relative_irf` to `irf_calibration` block.
- Conditional forecast
+ Command `conditional_forecast` now takes into account `histval`
block if present.
- Shock decomposition
+ New option `colormap` to `shocks_decomposition` for controlling
the color map used in the shocks decomposition graphs,
+ `shocks_decomposition` now accepts the `nograph` option,
+ New command `realtime_shock_decomposition` that for each period `T= [presample,...,nobs]`
allows computing the:
* realtime historical shock decomposition `Y(t|T)`, i.e. without observing data in `[T+1,...,nobs]`
* forecast shock decomposition `Y(T+k|T)`
* realtime conditional shock decomposition `Y(T+k|T+k)-Y(T+k|T)`
+ New block `shock_groups` that allows grouping shocks for the
`shock_decomposition` and `realtime_shock_decomposition` commands,
+ New command `plot_shock_decomposition` that allows plotting the
results from `shock_decomposition` and
`realtime_shock_decomposition` for different vintages and shock
groupings.
- Macroprocessor
+ Can now pass a macro-variable to the `@#include` macro directive,
+ New preprocessor flag `-I`, macro directive `@#includepath`, and
dynare config file block `[paths]` to pass a search path to the
macroprocessor to be used for file inclusion via `@#include`.
- Command line
+ New option `onlyclearglobals` (do not clear JIT compiled functions
with recent versions of Matlab),
+ New option `minimal_workspace` to use fewer variables in the
current workspace,
+ New option `params_derivs_order` allows limiting the order of the
derivatives with respect to the parameters that are calculated by
the preprocessor,
+ New command line option `mingw` to support the MinGW-w64 C/C++
Compiler from TDM-GCC for `use_dll`.
- dates/dseries/reporting classes
+ New methods `abs`, `cumprod` and `chain`,
+ New option `tableRowIndent` to `addTable`,
+ Reporting system revamped and made more efficient, dependency on
matlab2tikz has been dropped.
- Optimization algorithms
+ `mode_compute=2` Uses the simulated annealing as described by
Corana et al. (1987),
+ `mode_compute=101` Uses SOLVEOPT as described by Kuntsevich and
Kappel (1997),
+ `mode_compute=102` Uses `simulannealbnd` from Matlab's Global
Optimization Toolbox (if available),
+ New option `silent_optimizer` to shut off output from mode
computing/optimization,
+ New options `verbosity` and `SaveFiles` to control output and
saving of files during mode computing/optimization.
- LaTeX output
+ New command `write_latex_original_model`,
+ New option `write_equation_tags` to `write_latex_dynamic_model`
that allows printing the specified equation tags to the generate
LaTeX code,
+ New command `write_latex_parameter_table` that writes the names and
values of model parameters to a LaTeX table,
+ New command `write_latex_prior_table` that writes the descriptive
statistics about the prior distribution to a LaTeX table,
+ New command `collect_latex_files` that creates one compilable LaTeX
file containing all TeX-output.
- Misc.
+ Provides 64bit preprocessor,
+ Introduces new path management to avoid conflicts with other
toolboxes,
+ Full compatibility with Matlab 2014b's new graphic interface,
+ When using `model(linear)`, Dynare automatically checks
whether the model is truly linear,
+ `usedll`, the `msvc` option now supports `normcdf`, `acosh`,
`asinh`, and `atanh`,
+ New parallel option `NumberOfThreadsPerJob` for Windows nodes that
sets the number of threads assigned to each remote MATLAB/Octave
run,
+ Improved numerical performance of
`schur_statespace_transformation` for very large models,
+ The `all_values_required` option now also works with `histval`,
+ Add missing `horizon` option to `ms_forecast`,
+ BVAR now saves the marginal data density in
`oo_.bvar.log_marginal_data_density` and stores prior and
posterior information in `oo_.bvar.prior` and
`oo_.bvar.posterior`.
* Bugs and problems identified in version 4.4.3 and that have been fixed in version 4.5.0:
- BVAR models
+ `bvar_irf` could display IRFs in an unreadable way when they moved from
negative to positive values,
+ In contrast to what is stated in the documentation, the confidence interval
size `conf_sig` was 0.6 by default instead of 0.9.
- Conditional forecasts
+ The `conditional_forecast` command produced wrong results in calibrated
models when used at initial values outside of the steady state (given with
`initval`),
+ The `plot_conditional_forecast` option could produce unreadable figures if
the areas overlap,
+ The `conditional_forecast` command after MLE crashed,
+ In contrast to what is stated in the manual, the confidence interval size
`conf_sig` was 0.6 by default instead of 0.8.
+ Conditional forecasts were wrong when the declaration of endogenous
variables was not preceeding the declaration of the exogenous
variables and parameters.
- Discretionary policy
+ Dynare allowed running models where the number of instruments did not match
the number of omitted equations,
+ Dynare could crash in some cases when trying to display the solution,
+ Parameter dependence embedded via a `steady_state` was not taken into
account, typically resulting in crashes.
- dseries class
+ When subtracting a dseries object from a number, the number was instead
subtracted from the dseries object.
- DSGE-VAR models
+ Dynare crashed when estimation encountered non-finite values in the Jacobian
at the steady state,
+ The presence of a constant was not considered for degrees of freedom
computation of the Gamma function used during the posterior computation; due
to only affecting the constant term, results should be be unaffected, except
for model_comparison when comparing models with and without.
- Estimation command
+ In contrast to what was stated in the manual, the confidence interval size
`conf_sig` for `forecast` without MCMC was 0.6 by default instead of 0.9,
+ Calling estimation after identification could lead to crashes,
+ When using recursive estimation/forecasting and setting some elements of
`nobs` to be larger than the number of observations T in the data,
`oo_recursive_` contained additional cell entries that simply repeated the
results obtained for `oo_recursive_T`,
+ Computation of Bayesian smoother could crash for larger models when
requesting `forecast` or `filtered_variables`,
+ Geweke convergence diagnostics were not computed on the full MCMC chain when
the `load_mh_file` option was used,
+ The Geweke convergence diagnostics always used the default `taper_steps` and
`geweke_interval`,
+ Bayesian IRFs (`bayesian_irfs` option) could be displayed in an unreadable
way when they move from negative to positive values,
+ If `bayesian_irfs` was requested when `mh_replic` was too low to compute
HPDIs, plotting was crashing,
+ The x-axis value in `oo_.prior_density` for the standard deviation and
correlation of measurement errors was written into a field
`mearsurement_errors_*` instead of `measurement_errors_*`,
+ Using a user-defined `mode_compute` crashed estimation,
+ Option `mode_compute=10` did not work with infinite prior bounds,
+ The posterior variances and covariances computed by `moments_varendo` were
wrong for very large models due to a matrix erroneously being filled up with
zeros,
+ Using the `forecast` option with `loglinear` erroneously added the unlogged
steady state,
+ When using the `loglinear` option the check for the presence of a constant
was erroneously based on the unlogged steady state,
+ Estimation of `observation_trends` was broken as the trends specified as a
function of deep parameters were not correctly updated during estimation,
+ When using `analytic_derivation`, the parameter values were not set before
testing whether the steady state file changes parameter values, leading to
subsequent crashes,
+ If the steady state of an initial parameterization did not solve, the
observation equation could erroneously feature no constant when the
`use_calibration` option was used,
+ When computing posterior moments, Dynare falsely displayed that moment
computations are skipped, although the computation was performed correctly,
+ If `conditional_variance_decomposition` was requested, although all
variables contain unit roots, Dynare crashed instead of providing an error
message,
+ Computation of the posterior parameter distribution was erroneously based
on more draws than specified (there was one additional draw for every Markov
chain),
+ The estimation option `lyapunov=fixed_point` was broken,
+ Computation of `filtered_vars` with only one requested step crashed Dynare,
+ Option `kalman_algo=3` was broken with non-diagonal measurement error,
+ When using the diffuse Kalman filter with missing observations, an additive
factor log(2*pi) was missing in the last iteration step,
+ Passing of the `MaxFunEvals` and `InitialSimplexSize` options to
`mode_compute=8` was broken,
+ Bayesian forecasts contained initial conditions and had the wrong length in
both plots and stored variables,
+ Filtered variables obtained with `mh_replic=0`, ML, or
`calibrated_smoother` were padded with zeros at the beginning and end and
had the wrong length in stored variables,
+ Computation of smoothed measurement errors in Bayesian estimation was broken,
+ The `selected_variables_only` option (`mh_replic=0`, ML, or
`calibrated_smoother`) returned wrong results for smoothed, updated, and
filtered variables,
+ Combining the `selected_variables_only` option with forecasts obtained
using `mh_replic=0`, ML, or `calibrated_smoother` leaded to crashes,
+ `oo_.UpdatedVariables` was only filled when the `filtered_vars` option was specified,
+ When using Bayesian estimation with `filtered_vars`, but without
`smoother`, then `oo_.FilteredVariables` erroneously also contained filtered
variables at the posterior mean as with `mh_replic=0`,
+ Running an MCMC a second time in the same folder with a different number of
iterations could result in crashes due to the loading of stale files,
+ Results displayed after Bayesian estimation when not specifying
the `smoother` option were based on the parameters at the mode
from mode finding instead of the mean parameters from the
posterior draws. This affected the smoother results displayed, but
also calls to subsequent command relying on the parameters stored
in `M_.params` like `stoch_simul`,
+ The content of `oo_.posterior_std` after Bayesian estimation was based on
the standard deviation at the posterior mode, not the one from the MCMC, this
was not consistent with the reference manual,
+ When the initialization of an MCMC run failed, the metropolis.log file was
locked, requiring a restart of Matlab to restart estimation,
+ If the posterior mode was right at the corner of the prior bounds, the
initialization of the MCMC erroneously crashed,
+ If the number of dropped draws via `mh_drop` coincided with the number of
draws in a `_mh'-file`, `oo_.posterior.metropolis.mean` and
`oo_.posterior.metropolis.Variance` were NaN.
- Estimation and calibrated smoother
+ When using `observation_trends` with the `prefilter` option, the mean shift
due to the trend was not accounted for,
+ When using `first_obs`>1, the higher trend starting point of
`observation_trends` was not taken into account, leading, among other things,
to problems in recursive forecasting,
+ The diffuse Kalman smoother was crashing if the forecast error variance
matrix becomes singular,
+ The multivariate Kalman smoother provided incorrect state estimates when
all data for one observation are missing,
+ The multivariate diffuse Kalman smoother provided incorrect state estimates
when the `Finf` matrix becomes singular,
+ The univariate diffuse Kalman filter was crashing if the initial covariance
matrix of the nonstationary state vector is singular,
- Forecats
+ In contrast to what is stated in the manual, the confidence interval size
`conf_sig` was 0.6 by default instead of 0.9.
+ Forecasting with exogenous deterministic variables provided wrong decision
rules, yielding wrong forecasts.
+ Forecasting with exogenous deterministic variables crashed when the
`periods` option was not explicitly specified,
+ Option `forecast` when used with `initval` was using the initial values in
the `initval` block and not the steady state computed from these initial
values as the starting point of forecasts.
- Global Sensitivity Analysis
+ Sensitivity with ML estimation could result in crashes,
+ Option `mc` must be forced if `neighborhood_width` is used,
+ Fixed dimension of `stock_logpo` and `stock_ys`,
+ Incomplete variable initialization could lead to crashes with `prior_range=1`.
- Indentification
+ Identification did not correctly pass the `lik_init` option,
requiring the manual setting of `options_.diffuse_filter=1` in
case of unit roots,
+ Testing identification of standard deviations as the only
parameters to be estimated with ML leaded to crashes,
+ Automatic increase of the lag number for autocovariances when the
number of parameters is bigger than the number of non-zero moments
was broken,
+ When using ML, the asymptotic Hessian was not computed,
+ Checking for singular values when the eigenvectors contained only
one column did not work correctly,
- Model comparison
+ Selection of the `modifiedharmonicmean` estimator was broken,
- Optimal Simple Rules
+ When covariances were specified, variables that only entered with
their variance and no covariance term obtained a wrong weight,
resulting in wrong results,
+ Results reported for stochastic simulations after `osr` were based
on the last parameter vector encountered during optimization,
which does not necessarily coincide with the optimal parameter
vector,
+ Using only one (co)variance in the objective function resulted in crashes,
+ For models with non-stationary variables the objective function was computed wrongly.
- Ramsey policy
+ If a Lagrange multiplier appeared in the model with a lead or a lag
of more than one period, the steady state could be wrong.
+ When using an external steady state file, incorrect steady states
could be accepted,
+ When using an external steady state file with more than one
instrument, Dynare crashed,
+ When using an external steady state file and running `stoch_simul`
after `ramsey_planner`, an incorrect steady state was used,
+ When the number of instruments was not equal to the number of
omitted equations, Dynare crashed with a cryptic message,
+ The `planner_objective` accepted `varexo`, but ignored them for computations,
- Shock decomposition
+ Did not work with the `parameter_set=calibration` option if an
`estimated_params` block is present,
+ Crashed after MLE.
- Perfect foresight models
+ The perfect foresight solver could accept a complex solution
instead of continuing to look for a real-valued one,
+ The `initval_file` command only accepted column and not row vectors,
+ The `initval_file` command did not work with Excel files,
+ Deterministic simulations with one boundary condition crashed in
`solve_one_boundary` due to a missing underscore when passing
`options_.simul.maxit`,
+ Deterministic simulation with exogenous variables lagged by more
than one period crashed,
+ Termination criterion `maxit` was hard-coded for `solve_algo=0`
and could no be changed,
+ When using `block`/`bytecode`, relational operators could not be enforced,
+ When using `block` some exceptions were not properly handled,
leading to code crashes,
+ Using `periods=1` crashed the solver (bug only partially fixed).
- Smoothing
+ The univariate Kalman smoother returned wrong results when used
with correlated measurement error,
+ The diffuse smoother sometimes returned linear combinations of the
smoothed stochastic trend estimates instead of the original trend
estimates.
- Perturbation reduced form
+ In contrast to what is stated in the manual, the results of the
unconditional variance decomposition were only stored in
`oo_.gamma_y(nar+2)`, not in `oo_.variance_decomposition`,
+ Dynare could crash when the steady state could not be computed
when using the `loglinear` option,
+ Using `bytcode` when declared exogenous variables were not
used in the model leaded to crashes in stochastic simulations,
+ Displaying decision rules involving lags of auxiliary variables of
type 0 (leads>1) crashed.
+ The `relative_irf` option resulted in wrong output at `order>1` as
it implicitly relies on linearity.
- Displaying of the MH-history with the `internals` command crashed
if parameter names did not have same length.
- Dynare crashed when the user-defined steady state file returned an
error code, but not an conformable-sized steady state vector.
- Due to a bug in `mjdgges.mex` unstable parameter draws with
eigenvalues up to 1+1e-6 could be accepted as stable for the
purpose of the Blanchard-Kahn conditions, even if `qz_criterium<1`.
- The `use_dll` option on Octave for Windows required to pass a
compiler flag at the command line, despite the manual stating this
was not necessary.
- Dynare crashed for models with `block` option if the Blanchard-Kahn
conditions were not satisfied instead of generating an error
message.
- The `verbose` option did not work with `model(block)`.
- When falsely specifying the `model(linear)` for nonlinear models,
incorrect steady states were accepted instead of aborting.
- The `STEADY_STATE` operator called on model local variables
(so-called pound variables) did not work as expected.
- The substring operator in macro-processor was broken. The
characters of the substring could be mixed with random characters
from the memory space.
- Block decomposition could sometimes cause the preprocessor to crash.
- A bug when external functions were used in model local variables
that were contained in equations that required auxiliary
variable/equations led to crashes of Matlab.
- Sampling from the prior distribution for an inverse gamma II
distribution when `prior_trunc>0` could result in incorrect
sampling.
- Sampling from the prior distribution for a uniform distribution
when `prior_trunc>0` was ignoring the prior truncation.
- Conditional forecasts were wrong when the declaration of endogenous
variables was not preceeding the declaration of the exogenous
variables and parameters.
Announcement for Dynare 4.4.3 (on 2014-07-31)
=============================================
@ -988,7 +1785,7 @@ Here is a list of the main bugfixes since version 4.2.0:
* Option `conditional_variance_decomposition' of `stoch_simul' and
`estimation' has been fixed
* Automatic detrending now works in conjunction with the `EXPECTATION'
operator
@ -1029,7 +1826,7 @@ This release is compatible with MATLAB versions ranging from 6.5 (R13) to 7.11
Here is the list of major user-visible changes:
* New solution algorithms:
* New solution algorithms:
- Pruning for second order simulations has been added, as described in Kim,
Kim, Schaumburg and Sims (2008) [1,2]
@ -1066,7 +1863,7 @@ Here is the list of major user-visible changes:
- Syntax of deterministic shocks has changed: after the values keyword,
arbitrary expressions must be enclosed within parentheses (but numeric
constants are still accepted as is)
constants are still accepted as is)
* Various improvements:
@ -1095,7 +1892,7 @@ Here is the list of major user-visible changes:
from the console, it will replace graphical waitbars by text waitbars for
long computations
- Steady option "solve_algo=0" (uses fsolve()) now works under Octave
- Steady option "solve_algo=0" (uses fsolve()) now works under Octave
* For Emacs users:

View File

@ -303,7 +303,7 @@ After this, prepare the source and configure the build tree as described for Lin
- **NB**: If not compiling Dynare mex files for Octave, add ```--without-octave``` to the installation command
- **NB**: To compile the latest stable version of dynare, follow the same instructions as above, omitting the ```--HEAD``` argument
- **NB**: To update a ```--HEAD``` install of dynare you need to uninstall it then install it again: ```brew uninstall dynare; brew install dynare --HEAD```.
- **NB**: If you want to maintain a separate git directory of dynare, you can do a ```--HEAD``` install of dynare, then uninstall it. This will have the effect of bringing in all the dependencies you will need to then compile dynare from your git directory. Then, change to the git directory and type:
- **NB**: If you want to maintain a separate git directory of dynare, you can do a ```--HEAD``` install of dynare, then uninstall it. This will have the effect of bringing in all the dependencies you will need to then compile dynare from your git directory. (For `flex` and `bison` it may be necessary to symlink them via `brew link bison --force` and `brew link flex --force` as they are keg-only). Then, change to the git directory and type:
- ```autoreconf -si; ./configure --with-matlab=/Applications/MATLAB_R2015a.app MATLAB_VERSION=R2015a```, adjusting the Matlab path and version to accord with your version
- Once compilation is done, open Matlab and type the last line shown when you type ```brew info dynare``` in the Terminal window. With the typical Homebrew setup, this is:
- ```addpath /usr/local/opt/dynare/lib/dynare/matlab```

1
VERSION.in Normal file
View File

@ -0,0 +1 @@
@PACKAGE_VERSION@

View File

@ -18,7 +18,7 @@ dnl You should have received a copy of the GNU General Public License
dnl along with Dynare. If not, see <http://www.gnu.org/licenses/>.
AC_PREREQ([2.62])
AC_INIT([dynare], [4.5-unstable])
AC_INIT([dynare], [4.6-unstable])
AC_CONFIG_SRCDIR([preprocessor/DynareMain.cc])
AM_INIT_AUTOMAKE([1.11 -Wall -Wno-portability foreign no-dist-gzip dist-xz tar-pax])
@ -172,6 +172,7 @@ esac
AX_PTHREAD
AC_CONFIG_FILES([Makefile
VERSION
preprocessor/macro/Makefile
preprocessor/Makefile
doc/Makefile

View File

@ -12,7 +12,7 @@
\begin{document}
\title{BVAR models ``\`a la Sims'' in Dynare\thanks{Copyright \copyright~2007--2015 S\'ebastien
Villemot; \copyright~2016 S\'ebastien
Villemot; \copyright~2016--2017 S\'ebastien
Villemot and Johannes Pfeifer. 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
@ -27,7 +27,7 @@
\author{S\'ebastien Villemot\thanks{Paris School of Economics and
CEPREMAP.} \and Johannes Pfeifer\thanks{University of Cologne. E-mail: \href{mailto:jpfeifer@uni-koeln.de}{\texttt{jpfeifer@uni-koeln.de}}.}}
\date{First version: September 2007 \hspace{1cm} This version: October 2016}
\date{First version: September 2007 \hspace{1cm} This version: May 2017}
\maketitle
@ -545,7 +545,7 @@ Most results are stored for future use:
The syntax for computing impulse response functions is:
\medskip
\texttt{bvar\_irf(}\textit{number\_of\_periods},\textit{identification\_scheme}\texttt{);}
\texttt{bvar\_irf(}\textit{number\_of\_lags},\textit{identification\_scheme}\texttt{);}
\medskip
The \textit{identification\_scheme} option has two potential values
@ -556,7 +556,25 @@ The \textit{identification\_scheme} option has two potential values
Keep in mind that the first factorization of the covariance matrix is sensible to the ordering of the variables (as declared in the mod file with \verb+var+). This is not the case of the second factorization, but its structural interpretation is, at best, unclear (the Matrix square root of a covariance matrix, $\Sigma$, is the unique symmetric matrix $A$ such that $\Sigma = AA$).\newline
The mean, median, variance and confidence intervals for IRFs are saved in \texttt{oo\_.bvar.irf}
If you want to change the length of the IRFs plotted by the command, you can put\\
\medskip
\texttt{options\_.irf=40;}\\
\medskip
before the \texttt{bvar\_irf}-command. Similarly, to change the coverage of the highest posterior density intervals to e.g. 60\% you can put the command\\
\medskip
\texttt{options\_.bvar.conf\_sig=0.6;}\\
\medskip
there.\newline
The mean, median, variance, and confidence intervals for IRFs are saved in \texttt{oo\_.bvar.irf}
\section{Examples}

View File

@ -6114,7 +6114,7 @@ singularity is encountered, Dynare by default automatically switches to the univ
recursions as described by @cite{Herbst, 2015}. This setting is only used with
@code{kalman_algo=1} or @code{kalman_algo=3}. In case of using the diffuse Kalman
filter (@code{kalman_algo=3/lik_init=3}), the observables must be stationary. This option
is not yet compatible with @code{analytical_derivation}.
is not yet compatible with @ref{analytic_derivation}.
@item kalman_tol = @var{DOUBLE}
@anchor{kalman_tol} Numerical tolerance for determining the singularity of the covariance matrix of the prediction errors during the Kalman filter (minimum allowed reciprocal of the matrix condition number). Default value is @code{1e-10}
@ -6260,9 +6260,10 @@ where the model is ill-behaved. By default the original objective function is
used.
@item analytic_derivation
@anchor{analytic_derivation}
Triggers estimation with analytic gradient. The final hessian is also
computed analytically. Only works for stationary models without
missing observations.
missing observations, i.e. for @code{kalman_algo<3}.
@item ar = @var{INTEGER}
@xref{ar}. Only useful in conjunction with option @code{moments_varendo}.

View File

@ -101,7 +101,11 @@ void SystemResources::getRUS(double& load_avg, long int& pg_avail,
majflt = -1;
#endif
#if !defined(__MINGW32__) && !defined(__CYGWIN32__) && !defined(__CYGWIN__) && !defined(__MINGW64__) && !defined(__CYGWIN64__)
#define MINGCYGTMP (!defined(__MINGW32__) && !defined(__CYGWIN32__) && !defined(__CYGWIN__))
#define MINGCYG (MINGCYGTMP && !defined(__MINGW64__) && !defined(__CYGWIN64__))
#if MINGCYG
getloadavg(&load_avg, 1);
#else
load_avg = -1.0;

View File

@ -1,16 +1,15 @@
Format: http://www.debian.org/doc/packaging-manuals/copyright-format/1.0/
Upstream-Name: Dynare
Upstream-Contact: Dynare Team, whose members in 2016 are:
Upstream-Contact: Dynare Team, whose members in 2017 are:
Stéphane Adjemian <stephane.adjemian@univ-lemans.fr>
Houtan Bastani <houtan@dynare.org>
Michel Juillard <michel.juillard@mjui.fr>
Frédéric Karamé <frederic.karame@univ-lemans.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@dynare.org>
Marco Ratto <marco.ratto@ec.europa.eu>
Sébastien Villemot <sebastien.villemot@sciencespo.fr>
Source: http://www.dynare.org
Files: *
@ -205,7 +204,7 @@ License: permissive
Unlimited permission is granted to everyone to use, copy, modify or
distribute this software.
Files: matlab/utilities\graphics\distinguishable_colors.m
Files: matlab/utilities/graphics/distinguishable_colors.m
Copyright 2010-2011 by Timothy E. Holy
All rights reserved.
Redistribution and use in source and binary forms, with or without
@ -228,7 +227,7 @@ CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
POSSIBILITY OF SUCH DAMAGE.
Files: matlab/utilities\graphics\colorspace.m
Files: matlab/utilities/graphics/colorspace.m
Pascal Getreuer 2005-2010
All rights reserved.
Redistribution and use in source and binary forms, with or without

View File

@ -1,20 +0,0 @@
function display(t)
% Copyright (C) 2013 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/>.
fprintf('%s = <dynTimeIndex: %s>\n', inputname(1), int2str(t.index));

View File

@ -1,61 +0,0 @@
function t = dynTimeIndex() % --*-- Unitary tests --*--
% t = dynTimeIndex()
%
% Constructor for the dynTimeIndex class.
%
% INPUTS:
% None.
%
% OUTPUTS:
% * t, dynTimeIndex object.
%
% DESCRIPTION:
% The dynTimeIndex object is used to shift backward or forward dseries objects. For instance, if ts
% is a dseries object and t is a dynTimeIndex object then the following expressions are equivalent:
%
% us = ts.lag()
% us = ts.lag(1)
% us = lag(ts,1)
% us = ts(t-1)
%
% This class has only one member: t = int8(0) when instantiated.
% Copyright (C) 2013 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/>.
t = struct();
t.index = int8(0);
t = class(t,'dynTimeIndex');
%@test:1
%$ % Instantiate a dynTimeIndex object
%$ try
%$ u = dynTimeIndex();
%$ t(1) = 1;
%$ catch
%$ t(1) = 0;
%$ end
%$
%$ if t(1)
%$ t(2) = isa(u,'dynTimeIndex');
%$ end
%$
%$ T = all(t);
%@eof:1

View File

@ -1,62 +0,0 @@
function C = minus(A,B) % --*-- Unitary tests --*--
% C = minus(A,B)
%
% Overloads binary minus operator.
%
% INPUTS:
% * A, dynTimeIndex object.
% * B, integer scalar.
%
% OUTPUTS:
% * C, dynTimeIndex object.
%
% EXAMPLE:
%
% >> t = dynTimeIndex();
% >> t.index
%
% ans =
%
% 0
%
% >> s = t-1;
% >> s.index
%
% ans =
%
% -1
%
% Copyright (C) 2013 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 ~(isa(A,'dynTimeIndex') || isint(B))
error(['dynTimeIndex::plus: Input arguments ''' inputname(1) ''' and ''' inputname(2) ''' must be a dynTimeIndex object and an integer!'])
end
C = struct();
C.index = A.index-B;
C = class(C,'dynTimeIndex');
%@test:1
%$ a = dynTimeIndex();
%$ b = a-1;
%$ t(1) = isa(b,'dynTimeIndex');
%$ t(2) = isequal(b.index,-int8(1));
%$ T = all(t);
%@eof:1

View File

@ -1,74 +0,0 @@
function C = mpower(A,B) % --*-- Unitary tests --*--
% C = mpower(A,B)
%
% Overloads binary mpower operator (^).
%
% INPUTS :
% * A, dynTimeIndex object.
% * B, integer scalar.
%
% OUTPUTS :
% * C, dynTimeIndex object.
%
% EXAMPLE 1 :
%
% >> B = dynTimeIndex()-1;
% >> B
% B = <dynTimeIndex: -1>
% >> B^4
% ans = <dynTimeIndex: -4>
% >>
%
% EXAMPLE 2 :
% This method can be used to apply the lead and lag methods an arbitrary number of times to a dseries object. For instance, if
% ts is a dseries object, and if we define
%
% >> B = dynTimeIndex()-1;
% >> F = dynTimeIndex()+1;
%
% B and F can be used as lag and lead operators and the following syntax:
%
% >> us = ts(F^2);
%
% is equivalent to
%
% >> us = ts.lead(2)
%
% or
%
% >> us = ts.lead.lead
%
% Copyright (C) 2013 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 ~(isa(A,'dynTimeIndex') || isint(B))
error(['dynTimeIndex::mpower: Input arguments ''' inputname(1) ''' and ''' inputname(2) ''' must be a dynTimeIndex object and an integer!'])
end
C = struct();
C.index = A.index*B;
C = class(C,'dynTimeIndex');
%@test:1
%$ a = dynTimeIndex()+1;
%$ b = a^2;
%$ t(1) = isa(b,'dynTimeIndex');
%$ t(2) = isequal(b.index,int8(2));
%$ T = all(t);
%@eof:1

View File

@ -1,62 +0,0 @@
function C = plus(A,B) % --*-- Unitary tests --*--
% C = plus(A,B)
%
% Overloads binary plus operator.
%
% INPUTS:
% * A, dynTimeIndex object.
% * B, integer scalar.
%
% OUTPUTS:
% * C, dynTimeIndex object.
%
% EXAMPLE:
%
% >> t = dynTimeIndex();
% >> t.index
%
% ans =
%
% 0
%
% >> s = t+1;
% >> s.index
%
% ans =
%
% 1
%
% Copyright (C) 2013 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 ~(isa(A,'dynTimeIndex') || isint(B))
error(['dynTimeIndex::plus: Input arguments ''' inputname(1) ''' and ''' inputname(2) ''' must be a dynTimeIndex object and an integer!'])
end
C = struct();
C.index = A.index+B;
C = class(C,'dynTimeIndex');
%@test:1
%$ a = dynTimeIndex();
%$ b = a+1;
%$ t(1) = isa(b,'dynTimeIndex');
%$ t(2) = isequal(b.index,int8(1));
%$ T = all(t);
%@eof:1

View File

@ -1,20 +0,0 @@
function val = subsasgn(val, idx, rhs)
% Copyright (C) 2013 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/>.
error('dynTimeIndex::subsasgn: Members of dynTimeIndex class are private!')

View File

@ -1,54 +0,0 @@
function B = subsref(A,S) % --*-- Unitary tests --*--
% B = subsref(A,S)
%
% Overloads the subsref method for dynTimeIndex class. This method only allows to get
% the value of the field `index`.
% Copyright (C) 2013 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 length(S)>1
error('dynTimeIndex::subsref: Something is wrong in your syntax!')
end
if isequal(S.type,'.')
if isequal(S.subs,'index')
B = builtin('subsref', A, S(1));
else
error(['dynTimeIndex::subsref: ' S.subs ' is not a known member!'])
end
else
error('dynTimeIndex::subsref: Something is wrong in your syntax!')
end
%@test:1
%$ % Instantiate a dynTimeIndex object
%$ u = dynTimeIndex();
%$ try
%$ v = u.index;
%$ t(1) = 1;
%$ catch
%$ t(1) = 0;
%$ end
%$
%$ if t(1)
%$ t(2) = isequal(v,int8(0));
%$ end
%$
%$ T = all(t);
%@eof:1

View File

@ -8,7 +8,7 @@ function [AHess, DLIK, LIK] = AHessian(T,R,Q,H,P,Y,DT,DYss,DOm,DH,DP,start,mf,ka
% NOTE: the derivative matrices (DT,DR ...) are 3-dim. arrays with last
% dimension equal to the number of structural parameters
% Copyright (C) 2011-2016 Dynare Team
% Copyright (C) 2011-2017 Dynare Team
%
% This file is part of Dynare.
%
@ -39,7 +39,7 @@ F_singular = 1;
lik = zeros(smpl,1); % Initialization of the vector gathering the densities.
LIK = Inf; % Default value of the log likelihood.
if nargout > 1,
if nargout > 1
DLIK = zeros(k,1); % Initialization of the score.
end
AHess = zeros(k,k); % Initialization of the Hessian
@ -47,7 +47,7 @@ Da = zeros(mm,k); % State vector.
Dv = zeros(length(mf),k);
% for ii = 1:k
% DOm = DR(:,:,ii)*Q*transpose(R) + R*DQ(:,:,ii)*transpose(R) + R*Q*transpose(DR(:,:,ii));
% DOm = DR(:,:,ii)*Q*transpose(R) + R*DQ(:,:,ii)*transpose(R) + R*Q*transpose(DR(:,:,ii));
% end
while notsteady && t<smpl
@ -66,9 +66,7 @@ while notsteady && t<smpl
iF = inv(F);
K = P(:,mf)*iF;
lik(t) = log(det(F))+transpose(v)*iF*v;
[DK,DF,DP1] = computeDKalman(T,DT,DOm,P,DP,DH,mf,iF,K);
for ii = 1:k
Dv(:,ii) = -Da(mf,ii) - DYss(mf,ii);
Da(:,ii) = DT(:,:,ii)*(a+K*v) + T*(Da(:,ii)+DK(:,:,ii)*v + K*Dv(:,ii));
@ -81,7 +79,7 @@ while notsteady && t<smpl
if t>=start
AHess = AHess + Dv'*iF*Dv + .5*(vecDPmf' * kron(iF,iF) * vecDPmf);
end
a = T*(a+K*v);
a = T*(a+K*v);
P = T*(P-K*P(mf,:))*transpose(T)+Om;
DP = DP1;
end
@ -107,8 +105,8 @@ if t < smpl
end
end
if t>=start
AHess = AHess + Dv'*iF*Dv;
end
AHess = AHess + Dv'*iF*Dv;
end
a = T*(a+K*v);
lik(t) = transpose(v)*iF*v;
end
@ -124,17 +122,17 @@ if t < smpl
% H(ii,jj) = trace(iPmf*(.5*DP(mf,mf,ii)*iPmf*DP(mf,mf,jj) + Dv(:,ii)*Dv(:,jj)'));
% end
% end
end
end
AHess = -AHess;
if nargout > 1,
AHess = -AHess;
if nargout > 1
DLIK = DLIK/2;
end
% adding log-likelihhod constants
lik = (lik + pp*log(2*pi))/2;
LIK = sum(lik(start:end)); % Minus the log-likelihood.
% end of main function
% end of main function
function [DK,DF,DP1] = computeDKalman(T,DT,DOm,P,DP,DH,mf,iF,K)
@ -142,13 +140,11 @@ k = size(DT,3);
tmp = P-K*P(mf,:);
for ii = 1:k
DF(:,:,ii) = DP(mf,mf,ii) + DH(:,:,ii);
DF(:,:,ii) = DP(mf,mf,ii) + DH(:,:,ii);
DiF(:,:,ii) = -iF*DF(:,:,ii)*iF;
DK(:,:,ii) = DP(:,mf,ii)*iF + P(:,mf)*DiF(:,:,ii);
Dtmp = DP(:,:,ii) - DK(:,:,ii)*P(mf,:) - K*DP(mf,:,ii);
DP1(:,:,ii) = DT(:,:,ii)*tmp*T' + T*Dtmp*T' + T*tmp*DT(:,:,ii)' + DOm(:,:,ii);
end
% end of computeDKalman
% end of computeDKalman

View File

@ -1,4 +1,4 @@
function e = SPAimerr(c);
function e = SPAimerr(c)
% e = aimerr(c);
%
% Interpret the return codes generated by the aim routines.

View File

@ -8,9 +8,9 @@ function [b,rts,ia,nexact,nnumeric,lgroots,aimcode] = ...
% roots. This procedure will fail if the companion matrix is
% defective and does not have a linearly independent set of
% eigenvectors associated with the big roots.
%
%
% Input arguments:
%
%
% h Structural coefficient matrix (neq,neq*(nlag+1+nlead)).
% neq Number of equations.
% nlag Number of lags.
@ -19,9 +19,9 @@ function [b,rts,ia,nexact,nnumeric,lgroots,aimcode] = ...
% by numeric_shift and reduced_form.
% uprbnd Inclusive upper bound for the modulus of roots
% allowed in the reduced form.
%
%
% Output arguments:
%
%
% b Reduced form coefficient matrix (neq,neq*nlag).
% rts Roots returned by eig.
% ia Dimension of companion matrix (number of non-trivial
@ -57,7 +57,7 @@ function [b,rts,ia,nexact,nnumeric,lgroots,aimcode] = ...
% pages 472-489
b=[];rts=[];ia=[];nexact=[];nnumeric=[];lgroots=[];aimcode=[];
if(nlag<1 || nlead<1)
if(nlag<1 || nlead<1)
error('Aim_eig: model must have at least one lag and one lead');
end
% Initialization.
@ -66,26 +66,26 @@ bcols=neq*nlag;q=zeros(qrows,qcols);rts=zeros(qcols,1);
[h,q,iq,nexact]=SPExact_shift(h,q,iq,qrows,qcols,neq);
if (iq>qrows)
aimcode = 61;
return;
return
end
[h,q,iq,nnumeric]=SPNumeric_shift(h,q,iq,qrows,qcols,neq,condn);
if (iq>qrows)
aimcode = 62;
return;
return
end
[a,ia,js] = SPBuild_a(h,qcols,neq);
if (ia ~= 0)
if any(any(isnan(a))) || any(any(isinf(a)))
if any(any(isnan(a))) || any(any(isinf(a)))
disp('A is NAN or INF')
aimcode=63;
return
end
aimcode=63;
return
end
[w,rts,lgroots,flag_trouble]=SPEigensystem(a,uprbnd,min(length(js),qrows-iq+1));
if flag_trouble==1;
disp('Problem in SPEIG');
if flag_trouble==1
disp('Problem in SPEIG');
aimcode=64;
return
end
end
q = SPCopy_w(q,w,js,iq,qrows);
end
test=nexact+nnumeric+lgroots;

View File

@ -51,7 +51,7 @@ a(hrows,:) = hs(:,left);
% Delete inessential lags and build index array js. js indexes the
% columns in the big transition matrix that correspond to the
% essential lags in the model. They are the columns of q that will
% get the unstable left eigenvectors.
% get the unstable left eigenvectors.
js = 1:qcols;
zerocols = sum(abs(a)) == 0;

View File

@ -2,7 +2,7 @@ function q = SPCopy_w(q,w,js,iq,qrows)
% q = SPCopy_w(q,w,js,iq,qrows)
%
% Copy the eigenvectors corresponding to the largest roots into the
% remaining empty rows and columns js of q
% remaining empty rows and columns js of q
% Author: Gary Anderson
% Original file downloaded from:

View File

@ -30,7 +30,7 @@ function [w,rts,lgroots,flag_trouble] = SPEigensystem(a,uprbnd,rowsLeft)
% Journal of Economic Dynamics and Control, 2010, vol. 34, issue 3,
% pages 472-489
opts.disp=0;
opts.disp=0;
% next block is commented out because eigs() intermitently returns different rts
%try
% [w,d] = eigs(a',rowsLeft,'LM',opts);
@ -56,7 +56,7 @@ mag = abs(rts);
[mag,k] = sort(-mag);
rts = rts(k);
%end
flag_trouble=0;
flag_trouble=0;
%ws=SPSparse(w);
ws=sparse(w);
@ -65,7 +65,7 @@ ws = ws(:,k);
% Given a complex conjugate pair of vectors W = [w1,w2], there is a
% nonsingular matrix D such that W*D = real(W) + imag(W). That is to
% say, W and real(W)+imag(W) span the same subspace, which is all
% that aim cares about.
% that aim cares about.
ws = real(ws) + imag(ws);

View File

@ -44,4 +44,3 @@ while( any(zerorows) && iq <= qrows )
zerorows = find( sum(abs( hs(:,right)' ))==0 );
end
h=full(hs);

View File

@ -2,7 +2,7 @@ function scof = SPObstruct(cof,cofb,neq,nlag,nlead)
% scof = SPObstruct(cof,cofb,neq,nlag,nlead)
%
% Construct the coefficients in the observable structure.
%
%
% Input arguments:
%
% cof structural coefficients
@ -51,7 +51,7 @@ qs=sparse(q);
qs(1:rc,1:cc) = sparse(cofb);
qcols = neq*(nlag+nlead);
if( nlead > 1 )
if( nlead > 1 )
for i = 1:nlead-1
rows = i*neq + (1:neq);
qs(rows,:) = SPShiftright( qs((rows-neq),:), neq );

View File

@ -1,4 +1,4 @@
function [nonsing,b] = SPReduced_form(q,qrows,qcols,bcols,neq,condn);
function [nonsing,b] = SPReduced_form(q,qrows,qcols,bcols,neq,condn)
% [nonsing,b] = SPReduced_form(q,qrows,qcols,bcols,neq,b,condn);
%
% Compute reduced-form coefficient matrix, b.
@ -48,4 +48,3 @@ else %rescale by dividing row by maximal qr element
b = full(b);
end
end

View File

@ -3,7 +3,7 @@ function [y] = SPShiftright(x,n)
% [y] = shiftright(x,n)
%
% Shift the rows of x to the right by n columns, leaving zeros in the
% first n columns.
% first n columns.
% Original author: Gary Anderson
% Original file downloaded from:

View File

@ -1,20 +1,20 @@
function [dr,aimcode,rts]=dynAIMsolver1(jacobia_,M_,dr)
% function [dr,aimcode]=dynAIMsolver1(jacobia_,M_,dr)
% Maps Dynare jacobia to AIM 1st order model solver designed and developed by Gary ANderson
% Maps Dynare jacobia to AIM 1st order model solver designed and developed by Gary ANderson
% and derives the solution for gy=dr.hgx and gu=dr.hgu from the AIM outputs
% AIM System is given as a sum:
% i.e. for i=-$...+& SUM(Hi*xt+i)= £*zt, t = 0, . . . ,?
% AIM System is given as a sum:
% i.e. for i=-$...+& SUM(Hi*xt+i)= £*zt, t = 0, . . . ,?
% and its input as single array of matrices: [H-$... Hi ... H+&]
% and its solution as xt=SUM( Bi*xt+i) + @*£*zt for i=-$...-1
% with the output in form bb=[B-$... Bi ... B-1] and @=inv(Ho+H1*B-1)
% Dynare jacobian = [fy'-$... fy'i ... fy'+& fu']
% where [fy'-$... fy'i ... fy'+&]=[H-$... Hi ... H+&] and fu'= £
% and its solution as xt=SUM( Bi*xt+i) + @*£*zt for i=-$...-1
% with the output in form bb=[B-$... Bi ... B-1] and @=inv(Ho+H1*B-1)
% Dynare jacobian = [fy'-$... fy'i ... fy'+& fu']
% where [fy'-$... fy'i ... fy'+&]=[H-$... Hi ... H+&] and fu'= £
%
% INPUTS
% jacobia_ [matrix] 1st order derivative of the model
% jacobia_ [matrix] 1st order derivative of the model
% dr [matlab structure] Decision rules for stochastic simulations.
% M_ [matlab structure] Definition of the model.
%
% M_ [matlab structure] Definition of the model.
%
% OUTPUTS
% dr [matlab structure] Decision rules for stochastic simulations.
% aimcode [integer] 1: the model defines variables uniquely
@ -31,22 +31,22 @@ function [dr,aimcode,rts]=dynAIMsolver1(jacobia_,M_,dr)
% (c==63) e='Aim: A is NAN or INF.';
% (c==64) e='Aim: Problem in SPEIG.';
% else e='Aimerr: return code not properly specified';
%
%
% SPECIAL REQUIREMENTS
% Dynare use:
% Dynare use:
% 1) the lognormal block in DR1 is being invoked for some models and changing
% values of ghx and ghy. We need to return the AIM output
% values before that block and run the block with the current returned values
% of gy (i.e. dr.ghx) and gu (dr.ghu) if it is needed even when the AIM is used
% of gy (i.e. dr.ghx) and gu (dr.ghu) if it is needed even when the AIM is used
% (it does not depend on mjdgges output).
%
% 2) passing in aa={Q'|1}*jacobia_ can produce ~ one order closer
% results to the Dynare solutiion then when if plain jacobia_ is passed,
% i.e. diff < e-14 for aa and diff < *e-13 for jacobia_ if Q' is used.
%
% GP July 2008
% 2) passing in aa={Q'|1}*jacobia_ can produce ~ one order closer
% results to the Dynare solutiion then when if plain jacobia_ is passed,
% i.e. diff < e-14 for aa and diff < *e-13 for jacobia_ if Q' is used.
%
% GP July 2008
% Copyright (C) 2008-2012 Dynare Team
% Copyright (C) 2008-2017 Dynare Team
%
% This file is part of Dynare.
%
@ -69,8 +69,8 @@ lags=M_.maximum_endo_lag; % no of lags and leads
leads=M_.maximum_endo_lead;
klen=(leads+lags+1); % total lenght
theAIM_H=zeros(neq, neq*klen); % alloc space
lli=M_.lead_lag_incidence';
% "sparse" the compact jacobia into AIM H aray of matrices
lli=M_.lead_lag_incidence';
% "sparse" the compact jacobia into AIM H aray of matrices
% without exogenous shoc
theAIM_H(:,find(lli(:)))=jacobia_(:,nonzeros(lli(:)));
condn = 1.e-10;%SPAmalg uses this in zero tests
@ -102,7 +102,7 @@ if aimcode==1 %if OK
col_order=[((i-1)*neq)+dr.order_var' col_order];
end
bb_ord= bb(dr.order_var,col_order); % derive ordered gy
% variables are present in the state space at the lag at which they
% appear and at all smaller lags. The are ordered from smaller to
% higher lag (reversed order of M_.lead_lag_incidence rows for lagged
@ -117,7 +117,7 @@ if aimcode==1 %if OK
%theH0= theAIM_H (:,lags*neq+1: (lags+1)*neq);
% theHP= theAIM_H (:,(M_.maximum_endo_lag+1)*neq+1: (M_.maximum_endo_lag+2)*neq);
%theHP= theAIM_H (:,(lags+1)*neq+1: (lags+2)*neq);
theAIM_Psi= - jacobia_(:, size(nonzeros(lli(:)))+1:end);%
theAIM_Psi= - jacobia_(:, size(nonzeros(lli(:)))+1:end);%
%? = inv(H0 + H1B1)
%phi= (theH0+theHP*sparse(bb(:,(lags-1)*neq+1:end)))\eye( neq);
%AIM_ghu=phi*theAIM_Psi;
@ -137,8 +137,8 @@ else
if aimcode < 1 || aimcode > 5 % too big exception, use mjdgges
error('Error in AIM: aimcode=%d ; %s', aimcode, err);
end
% if aimcode > 5
% if aimcode > 5
% disp(['Error in AIM: aimcode=' sprintf('%d : %s',aimcode, err)]);
% aimcode=5;
% end
% end
end

View File

@ -51,7 +51,7 @@ function [dr,info]=AIM_first_order_solver(jacobia,M,dr,qz_criterium)
%! @end deftypefn
%@eod:
% Copyright (C) 2001-2016 Dynare Team
% Copyright (C) 2001-2017 Dynare Team
%
% This file is part of Dynare.
%
@ -91,11 +91,10 @@ if nba ~= nsfwrd
if nba > nsfwrd
temp = temp(nd-nba+1:nd-nsfwrd)-1-qz_criterium;
info(1) = 3;
elseif nba < nsfwrd;
elseif nba < nsfwrd
temp = temp(nd-nsfwrd+1:nd-nba)-1-qz_criterium;
info(1) = 4;
end
info(2) = temp'*temp;
return
end

View File

@ -2,7 +2,7 @@ function [DirectoryName, info] = CheckPath(type,dname)
% Creates the subfolder "./M_.dname/type" if it does not exist yet.
%
% INPUTS
% type [string] Name of the subfolder.
% type [string] Name of the subfolder.
% dname [string] Name of the directory
%
% OUTPUTS
@ -12,7 +12,7 @@ function [DirectoryName, info] = CheckPath(type,dname)
% SPECIAL REQUIREMENTS
% none
% Copyright (C) 2005-2013 Dynare Team
% Copyright (C) 2005-2017 Dynare Team
%
% This file is part of Dynare.
%

View File

@ -16,7 +16,7 @@ function CutSample(M_, options_, estim_params_)
% SPECIAL REQUIREMENTS
% none
% Copyright (C) 2005-2015 Dynare Team
% Copyright (C) 2005-2017 Dynare Team
%
% This file is part of Dynare.
%
@ -66,7 +66,7 @@ if (TotalNumberOfMhFiles-1)-(FirstMhFile+1)+1 > 0
elseif TotalNumberOfMhFiles == 1
record.KeepedDraws.Distribution = [];
elseif TotalNumberOfMhFiles == 2 && FirstMhFile > 1
record.KeepedDraws.Distribution = [MAX_nruns-FirstLine+1 ; record.MhDraws(end,3)];
record.KeepedDraws.Distribution = [MAX_nruns-FirstLine+1 ; record.MhDraws(end,3)];
end
% Save updated mh-history file.

View File

@ -1,8 +1,8 @@
function [alphahat,etahat,epsilonhat,ahat,SteadyState,trend_coeff,aK,T,R,P,PK,decomp,trend_addition,state_uncertainty,M_,oo_,options_,bayestopt_] = DsgeSmoother(xparam1,gend,Y,data_index,missing_value,M_,oo_,options_,bayestopt_,estim_params_)
% Estimation of the smoothed variables and innovations.
%
% INPUTS
% o xparam1 [double] (p*1) vector of (estimated) parameters.
% Estimation of the smoothed variables and innovations.
%
% INPUTS
% o xparam1 [double] (p*1) vector of (estimated) parameters.
% o gend [integer] scalar specifying the number of observations ==> varargin{1}.
% o data [double] (n*T) matrix of data.
% o data_index [cell] 1*smpl cell of column vectors of indices.
@ -12,7 +12,7 @@ function [alphahat,etahat,epsilonhat,ahat,SteadyState,trend_coeff,aK,T,R,P,PK,de
% o options_ [structure] describing the options
% o bayestopt_ [structure] describing the priors
% o estim_params_ [structure] characterizing parameters to be estimated
%
%
% OUTPUTS
% o alphahat [double] (m*T) matrix, smoothed endogenous variables (a_{t|T}) (decision-rule order)
% o etahat [double] (r*T) matrix, smoothed structural shocks (r>=n is the number of shocks).
@ -29,36 +29,36 @@ function [alphahat,etahat,epsilonhat,ahat,SteadyState,trend_coeff,aK,T,R,P,PK,de
% matrices (meaningless for periods 1:d) (decision-rule order)
% o decomp (K*m*r*(T+K)) 4D array of shock decomposition of k-step ahead
% filtered variables (decision-rule order)
% o trend_addition [double] (n*T) pure trend component; stored in options_.varobs order
% o trend_addition [double] (n*T) pure trend component; stored in options_.varobs order
% o state_uncertainty [double] (K,K,T) array, storing the uncertainty
% about the smoothed state (decision-rule order)
% o M_ [structure] decribing the model
% o oo_ [structure] storing the results
% o options_ [structure] describing the options
% o bayestopt_ [structure] describing the priors
%
%
% Notes:
% m: number of endogenous variables (M_.endo_nbr)
% T: number of Time periods (options_.nobs)
% r: number of strucural shocks (M_.exo_nbr)
% n: number of observables (length(options_.varobs))
% K: maximum forecast horizon (max(options_.nk))
%
%
% To get variables that are stored in decision rule order in order of declaration
% as in M_.endo_names, ones needs code along the lines of:
% variables_declaration_order(dr.order_var,:) = alphahat
%
% Defines bayestopt_.mf = bayestopt_.smoother_mf (positions of observed variables
% and requested smoothed variables in decision rules (decision rule order)) and
%
% Defines bayestopt_.mf = bayestopt_.smoother_mf (positions of observed variables
% and requested smoothed variables in decision rules (decision rule order)) and
% passes it back via global variable
%
% ALGORITHM
% Diffuse Kalman filter (Durbin and Koopman)
%
% ALGORITHM
% Diffuse Kalman filter (Durbin and Koopman)
%
% SPECIAL REQUIREMENTS
% None
% Copyright (C) 2006-2016 Dynare Team
% Copyright (C) 2006-2017 Dynare Team
%
% This file is part of Dynare.
%
@ -97,7 +97,7 @@ end
%------------------------------------------------------------------------------
% 2. call model setup & reduction program
%------------------------------------------------------------------------------
oldoo.restrict_var_list = oo_.dr.restrict_var_list;
oldoo.restrict_var_list = oo_.dr.restrict_var_list;
oldoo.restrict_columns = oo_.dr.restrict_columns;
oo_.dr.restrict_var_list = bayestopt_.smoother_var_list;
oo_.dr.restrict_columns = bayestopt_.smoother_restrict_columns;
@ -133,8 +133,8 @@ mf = bayestopt_.mf;
% ------------------------------------------------------------------------------
% 3. Initial condition of the Kalman filter
% ------------------------------------------------------------------------------
%
% Here, Pinf and Pstar are determined. If the model is stationary, determine
%
% Here, Pinf and Pstar are determined. If the model is stationary, determine
% Pstar as the solution of the Lyapunov equation and set Pinf=[] (Notation follows
% Koopman/Durbin (2003), Journal of Time Series Analysis 24(1))
%
@ -256,7 +256,7 @@ if kalman_algo == 2 || kalman_algo == 4
else
Pstar = blkdiag(Pstar,H);
if ~isempty(Pinf)
Pinf = blkdiag(Pinf,zeros(vobs));
Pinf = blkdiag(Pinf,zeros(vobs));
end
end
%now reset H to 0
@ -265,7 +265,7 @@ if kalman_algo == 2 || kalman_algo == 4
%do nothing, state vector was already expanded
end
end
[alphahat,epsilonhat,etahat,ahat,P,aK,PK,decomp,state_uncertainty] = missing_DiffuseKalmanSmootherH3_Z(ST, ...
Z,R1,Q,diag(H), ...
Pinf,Pstar,data1,vobs,np,smpl,data_index, ...
@ -282,7 +282,7 @@ if expanded_state_vector_for_univariate_filter && (kalman_algo == 2 || kalman_al
ahat = ahat(k,:);
aK = aK(:,k,:,:);
epsilonhat=etahat(end-vobs+1:end,:);
etahat=etahat(1:end-vobs,:);
etahat=etahat(1:end-vobs,:);
if ~isempty(PK)
PK = PK(:,k,k,:);
end

View File

@ -5,18 +5,18 @@ function Draws = GetAllPosteriorDraws(column, FirstMhFile, FirstLine, TotalNumbe
%
% INPUTS
% column: column
% FirstMhFile: first mh file
% FirstMhFile: first mh file
% FirstLine: first line
% TotalNumberOfMhFile: total number of mh file
% TotalNumberOfMhFile: total number of mh file
% NumberOfDraws: number of draws
% OUTPUTS
% Draws: draws from posterior distribution
%
%
% SPECIAL REQUIREMENTS
% none
% Copyright (C) 2005-2011 Dynare Team
% Copyright (C) 2005-2017 Dynare Team
%
% This file is part of Dynare.
%
@ -55,7 +55,7 @@ if nblck>1 && nargin<6
iline = 1;
end
end
else
else
for blck = 1:nblck
iline=iline0;
for file = FirstMhFile:TotalNumberOfMhFile

View File

@ -5,15 +5,15 @@ function [xparams, logpost] = GetOneDraw(type)
% INPUTS
% type: [string] 'posterior': draw from MCMC draws
% 'prior': draw from prior
%
%
% OUTPUTS
% xparams: vector of estimated parameters (drawn from posterior or prior distribution)
% logpost: log of the posterior density of this parameter vector
%
%
% SPECIAL REQUIREMENTS
% none
% Copyright (C) 2005-2015 Dynare Team
% Copyright (C) 2005-2017 Dynare Team
%
% This file is part of Dynare.
%

View File

@ -1,6 +1,6 @@
function [mean,variance] = GetPosteriorMeanVariance(M,drop)
% Copyright (C) 2012-2016 Dynare Team
% Copyright (C) 2012-2017 Dynare Team
%
% This file is part of Dynare.
%

View File

@ -1,22 +1,22 @@
function oo_ = GetPosteriorParametersStatistics(estim_params_, M_, options_, bayestopt_, oo_, pnames)
% This function prints and saves posterior estimates after the mcmc
% (+updates of oo_ & TeX output).
%
% INPUTS
% estim_params_ [structure]
% (+updates of oo_ & TeX output).
%
% INPUTS
% estim_params_ [structure]
% M_ [structure]
% options_ [structure]
% bayestopt_ [structure]
% oo_ [structure]
% pnames [char] Array of char, names of the prior shapes available
%
% OUTPUTS
% oo_ [structure]
%
% OUTPUTS
% oo_ [structure]
%
% SPECIAL REQUIREMENTS
% None.
% Copyright (C) 2006-2016 Dynare Team
% Copyright (C) 2006-2017 Dynare Team
%
% This file is part of Dynare.
%
@ -34,7 +34,7 @@ function oo_ = GetPosteriorParametersStatistics(estim_params_, M_, options_, bay
% along with Dynare. If not, see <http://www.gnu.org/licenses/>.
%if ~options_.mh_replic && options_.load_mh_file
% load([M_.fname '_results.mat'],'oo_');
% load([M_.fname '_results.mat'],'oo_');
%end
TeX = options_.TeX;
@ -48,7 +48,7 @@ nx = nvx+nvn+ncx+ncn+np;
MetropolisFolder = CheckPath('metropolis',M_.dname);
OutputFolder = CheckPath('Output',M_.dname);
FileName = M_.fname;
FileName = M_.fname;
load_last_mh_history_file(MetropolisFolder,FileName);
@ -108,14 +108,14 @@ if np
[post_mean, post_median, post_var, hpd_interval, post_deciles, ...
density] = posterior_moments(Draws,1,options_.mh_conf_sig);
name = bayestopt_.name{ip};
oo_ = Filloo(oo_,name,type,post_mean,hpd_interval,post_median,post_var,post_deciles,density);
oo_ = Filloo(oo_,name,type,post_mean,hpd_interval,post_median,post_var,post_deciles,density);
end
end
disp(sprintf(pformat,header_width,name,bayestopt_.p1(ip),...
post_mean, ...
hpd_interval, ...
pnames(bayestopt_.pshape(ip)+1,:), ...
bayestopt_.p2(ip)));
bayestopt_.p2(ip)));
if TeX
k = estim_params_.param_vals(i,1);
name = deblank(M_.param_names_tex(k,:));
@ -163,7 +163,7 @@ if nvx
end
disp(sprintf(pformat,header_width,name,bayestopt_.p1(ip),post_mean,hpd_interval,...
pnames(bayestopt_.pshape(ip)+1,:),bayestopt_.p2(ip)));
if TeX,
if TeX
name = deblank(M_.exo_names_tex(k,:));
TeXCore(fid,name,deblank(pnames(bayestopt_.pshape(ip)+1,:)),bayestopt_.p1(ip),...
bayestopt_.p2(ip),post_mean,sqrt(post_var),hpd_interval);
@ -171,7 +171,7 @@ if nvx
ip = ip+1;
end
if TeX
TeXEnd(fid,2,'standard deviation of structural shocks');
TeXEnd(fid,2,'standard deviation of structural shocks');
end
end
if nvn
@ -213,7 +213,7 @@ if nvn
ip = ip+1;
end
if TeX
TeXEnd(fid,3,'standard deviation of measurement errors');
TeXEnd(fid,3,'standard deviation of measurement errors');
end
end
if ncx
@ -311,15 +311,15 @@ if ncn
end
disp(sprintf(pformat, header_width,name,bayestopt_.p1(ip),post_mean,hpd_interval, ...
pnames(bayestopt_.pshape(ip)+1,:),bayestopt_.p2(ip)));
if TeX,
if TeX
name = ['(',deblank(M_.endo_names_tex(k1,:)) ',' deblank(M_.endo_names_tex(k2,:)),')'];
TeXCore(fid,name,deblank(pnames(bayestopt_.pshape(ip)+1,:)),bayestopt_.p1(ip),...
bayestopt_.p2(ip),post_mean,sqrt(post_var),hpd_interval);
bayestopt_.p2(ip),post_mean,sqrt(post_var),hpd_interval);
end
ip = ip+1;
end
if TeX
TeXEnd(fid,5,'correlation of measurement errors');
TeXEnd(fid,5,'correlation of measurement errors');
end
end
@ -359,7 +359,7 @@ fid = fidTeX;
function TeXCore(fid,name,shape,priormean,priorstd,postmean,poststd,hpd)
fprintf(fid,['$%s$ & %s & %7.3f & %6.4f & %7.3f& %6.4f & %7.4f & %7.4f \\\\ \n'],...
fprintf(fid,['$%s$ & %s & %7.3f & %6.4f & %7.3f& %6.4f & %7.4f & %7.4f \\\\ \n'],...
name,...
shape,...
priormean,...
@ -371,7 +371,7 @@ fprintf(fid,['$%s$ & %s & %7.3f & %6.4f & %7.3f& %6.4f & %7.4f & %7.4f \\\\ \n']
function TeXEnd(fid,fnum,title)
fprintf(fid,'\\end{longtable}\n ');
fprintf(fid,'\\end{longtable}\n ');
fprintf(fid,'\\end{center}\n');
fprintf(fid,'%% End of TeX file.\n');
fclose(fid);

View File

@ -1,6 +1,6 @@
function MakeAllFigures(NumberOfPlots,Caption,FigureProperties,Info)
% Copyright (C) 2005-2009 Dynare Team
% Copyright (C) 2005-2017 Dynare Team
%
% This file is part of Dynare.
%
@ -19,11 +19,11 @@ function MakeAllFigures(NumberOfPlots,Caption,FigureProperties,Info)
global M_ options_
FigHandle = figure('Name',FigureProperties.Name);
FigHandle = figure('Name',FigureProperties.Name);
NAMES = cell(NumberOfPlots,1);
if options_.TeX
TeXNAMES = cell(NumberOfPlots,1);
TeXNAMES = cell(NumberOfPlots,1);
end
if NumberOfPlots == 9
@ -53,7 +53,7 @@ elseif NumberOfPlots == 2
elseif NumberOfPlots == 1
nr = 1;
nc = 1;
end
end
for plt = 1:NumberOfPlots
eval(['NumberOfCurves = Info.Box' int2str(plt) '.Number;'])
@ -138,7 +138,7 @@ for plt = 1:NumberOfPlots
set(hh,'Color','r','LineStyle','-','LineWidth',2)
%
%
end
end
end
axis([xmin xmax ymin ymax])
title(NAMES{plt})
@ -150,14 +150,14 @@ if Info.SaveFormat.Eps
if isempty(Info.SaveFormat.Name)
eval(['print -depsc2 ' M_.fname Info.SaveFormat.GenericName int2str(Info.SaveFormat.Number) '.eps']);
else
eval(['print -depsc2 ' M_.fname Info.SaveFormat.GenericName Info.SaveFormat.Name '.eps']);
eval(['print -depsc2 ' M_.fname Info.SaveFormat.GenericName Info.SaveFormat.Name '.eps']);
end
end
if Info.SaveFormat.Pdf && ~isoctave
if isempty(Info.SaveFormat.Name)
eval(['print -dpdf ' M_.fname Info.SaveFormat.GenericName int2str(Info.SaveFormat.Number)]);
else
eval(['print -dpdf ' M_.fname Info.SaveFormat.GenericName Info.SaveFormat.Name]);
eval(['print -dpdf ' M_.fname Info.SaveFormat.GenericName Info.SaveFormat.Name]);
end
end
if Info.SaveFormat.Fig && ~isoctave

View File

@ -4,15 +4,15 @@ function oo_ = PlotPosteriorDistributions(estim_params_, M_, options_, bayestopt
% plots posterior distributions
%
% INPUTS
% estim_params_ [structure]
% estim_params_ [structure]
% M_ [structure]
% options_ [structure]
% options_ [structure]
% bayestopt_ [structure]
% oo_ [structure]
%
%
% OUTPUTS
% oo_ [structure]
%
% oo_ [structure]
%
% SPECIAL REQUIREMENTS
% none
@ -49,7 +49,7 @@ nn = sqrt(MaxNumberOfPlotPerFigure);
figurename = 'Priors and posteriors';
if TeX && any(strcmp('eps',cellstr(options_.graph_format)))
if TeX && any(strcmp('eps',cellstr(options_.graph_format)))
fidTeX = fopen([OutputDirectoryName '/' M_.fname '_PriorsAndPosteriors.tex'],'w');
fprintf(fidTeX,'%% TeX eps-loader file generated by PlotPosteriorDistributions.m (Dynare).\n');
fprintf(fidTeX,['%% ' datestr(now,0) '\n']);
@ -78,9 +78,9 @@ for i=1:npar
end
end
[x2,f2,abscissa,dens,binf2,bsup2] = draw_prior_density(i,bayestopt_);
top2 = max(f2);
top2 = max(f2);
if i <= nvx
name = deblank(M_.exo_names(estim_params_.var_exo(i,1),:));
name = deblank(M_.exo_names(estim_params_.var_exo(i,1),:));
x1 = oo_.posterior_density.shocks_std.(name)(:,1);
f1 = oo_.posterior_density.shocks_std.(name)(:,2);
oo_.prior_density.shocks_std.(name)(:,1) = x2;
@ -96,18 +96,18 @@ for i=1:npar
oo_.prior_density.measurement_errors_std.(name)(:,2) = f2;
if ~options_.mh_posterior_mode_estimation
pmod = oo_.posterior_mode.measurement_errors_std.(name);
end
end
elseif i <= nvx+nvn+ncx
j = i - (nvx+nvn);
k1 = estim_params_.corrx(j,1);
k2 = estim_params_.corrx(j,2);
name = [deblank(M_.exo_names(k1,:)) '_' deblank(M_.exo_names(k2,:))];
name = [deblank(M_.exo_names(k1,:)) '_' deblank(M_.exo_names(k2,:))];
x1 = oo_.posterior_density.shocks_corr.(name)(:,1);
f1 = oo_.posterior_density.shocks_corr.(name)(:,2);
oo_.prior_density.shocks_corr.(name)(:,1) = x2;
oo_.prior_density.shocks_corr.(name)(:,2) = f2;
if ~options_.mh_posterior_mode_estimation
pmod = oo_.posterior_mode.shocks_corr.(name);
pmod = oo_.posterior_mode.shocks_corr.(name);
end
elseif i <= nvx+nvn+ncx+ncn
j = i - (nvx+nvn+ncx);
@ -151,13 +151,13 @@ for i=1:npar
title(nam,'Interpreter','none');
hold off;
drawnow
if subplotnum == MaxNumberOfPlotPerFigure || i == npar;
if subplotnum == MaxNumberOfPlotPerFigure || i == npar
dyn_saveas(hfig,[OutputDirectoryName '/' M_.fname '_PriorsAndPosteriors' int2str(figunumber)],options_.nodisplay,options_.graph_format);
if TeX && any(strcmp('eps',cellstr(options_.graph_format)))
fprintf(fidTeX,'\\begin{figure}[H]\n');
for j = 1:size(NAMES,1)
fprintf(fidTeX,'\\psfrag{%s}[1][][0.5][0]{%s}\n',deblank(NAMES(j,:)),deblank(TeXNAMES(j,:)));
end
end
fprintf(fidTeX,'\\centering\n');
fprintf(fidTeX,'\\includegraphics[width=%2.2f\\textwidth]{%s/%s_PriorsAndPosteriors%s}\n',options_.figures.textwidth*min(subplotnum/nn,1),OutputDirectoryName,M_.fname,int2str(figunumber));
fprintf(fidTeX,'\\caption{Priors and posteriors.}');

View File

@ -16,7 +16,7 @@ function PosteriorIRF(type)
% functions associated with it(the _core1 and _core2).
% See also the comments posterior_sampler.m funtion.
% Copyright (C) 2006-2016 Dynare Team
% Copyright (C) 2006-2017 Dynare Team
%
% This file is part of Dynare.
%
@ -178,7 +178,7 @@ if strcmpi(type,'posterior')
end
end
if ~strcmpi(type,'prior'),
if ~strcmpi(type,'prior')
localVars.x=x;
end
@ -202,16 +202,16 @@ localVars.ifil2=ifil2;
localVars.MhDirectoryName=MhDirectoryName;
% Like sequential execution!
if isnumeric(options_.parallel),
if isnumeric(options_.parallel)
[fout] = PosteriorIRF_core1(localVars,1,B,0);
nosaddle = fout.nosaddle;
else
% Parallel execution!
[nCPU, totCPU, nBlockPerCPU] = distributeJobs(options_.parallel, 1, B);
for j=1:totCPU-1,
for j=1:totCPU-1
nfiles = ceil(nBlockPerCPU(j)/MAX_nirfs_dsge);
NumberOfIRFfiles_dsge(j+1) =NumberOfIRFfiles_dsge(j)+nfiles;
if MAX_nirfs_dsgevar,
if MAX_nirfs_dsgevar
nfiles = ceil(nBlockPerCPU(j)/MAX_nirfs_dsgevar);
else
nfiles=0;
@ -236,8 +236,8 @@ else
NamFileInput(1,:) = {'',[M_.fname '_static.m']};
NamFileInput(2,:) = {'',[M_.fname '_dynamic.m']};
NamFileInput(3,:) = {'',[M_.fname '_set_auxiliary_variables.m']};
if options_.steadystate_flag,
if options_.steadystate_flag == 1,
if options_.steadystate_flag
if options_.steadystate_flag == 1
NamFileInput(length(NamFileInput)+1,:)={'',[M_.fname '_steadystate.m']};
else
NamFileInput(length(NamFileInput)+1,:)={'',[M_.fname '_steadystate2.m']};
@ -245,7 +245,7 @@ else
end
[fout] = masterParallel(options_.parallel, 1, B,NamFileInput,'PosteriorIRF_core1', localVars, globalVars, options_.parallel_info);
nosaddle=0;
for j=1:length(fout),
for j=1:length(fout)
nosaddle = nosaddle + fout(j).nosaddle;
end
@ -408,9 +408,9 @@ if ~options_.nograph && ~options_.no_graph.posterior
for jj=1:nvar
if max(abs(MeanIRF(:,jj,ii))) >= options_.impulse_responses.plot_threshold
subplotnum = subplotnum+1;
if subplotnum == 1
fprintf(fidTeX,'\\begin{figure}[H]\n');
if subplotnum == 1
fprintf(fidTeX,'\\begin{figure}[H]\n');
end
name = deblank(varlist(jj,:));
texname = deblank(varlist_TeX(jj,:));
@ -418,7 +418,7 @@ if ~options_.nograph && ~options_.no_graph.posterior
end
if subplotnum == MaxNumberOfPlotPerFigure || (jj == nvar && subplotnum> 0)
figunumber = figunumber+1;
fprintf(fidTeX,'\\centering \n');
fprintf(fidTeX,'\\includegraphics[width=%2.2f\\textwidth]{%s/%s_Bayesian_IRF_%s_%d}\n',options_.figures.textwidth*min(subplotnum/nn,1),DirectoryName,M_.fname,deblank(tit(ii,:)),figunumber);
if options_.relative_irf
@ -429,9 +429,9 @@ if ~options_.nograph && ~options_.no_graph.posterior
fprintf(fidTeX,'\\label{Fig:BayesianIRF:%s:%d}\n',deblank(tit(ii,:)),figunumber);
fprintf(fidTeX,'\\end{figure}\n');
fprintf(fidTeX,' \n');
subplotnum = 0;
end
end
end
end
fprintf(fidTeX,'%% End of TeX file.\n');
@ -443,11 +443,11 @@ if ~options_.nograph && ~options_.no_graph.posterior
% Comment for testing!
if ~isoctave
if isnumeric(options_.parallel) || (M_.exo_nbr*ceil(size(varlist,1)/MaxNumberOfPlotPerFigure))<8,
if isnumeric(options_.parallel) || (M_.exo_nbr*ceil(size(varlist,1)/MaxNumberOfPlotPerFigure))<8
[fout] = PosteriorIRF_core2(localVars,1,M_.exo_nbr,0);
else
isRemoteOctave = 0;
for indPC=1:length(options_.parallel),
for indPC=1:length(options_.parallel)
isRemoteOctave = isRemoteOctave + (findstr(options_.parallel(indPC).MatlabOctavePath, 'octave'));
end
if isRemoteOctave

View File

@ -23,7 +23,7 @@ function myoutput=PosteriorIRF_core1(myinputs,fpar,B,whoiam, ThisMatlab)
% SPECIAL REQUIREMENTS.
% None.
%
% Copyright (C) 2006-2016 Dynare Team
% Copyright (C) 2006-2017 Dynare Team
%
% This file is part of Dynare.
%
@ -43,7 +43,7 @@ function myoutput=PosteriorIRF_core1(myinputs,fpar,B,whoiam, ThisMatlab)
global options_ estim_params_ oo_ M_ bayestopt_ dataset_ dataset_info
if nargin<4,
if nargin<4
whoiam=0;
end
@ -55,7 +55,7 @@ irun =myinputs.irun;
irun2=myinputs.irun2;
npar=myinputs.npar;
type=myinputs.type;
if ~strcmpi(type,'prior'),
if ~strcmpi(type,'prior')
x=myinputs.x;
end
@ -102,7 +102,7 @@ end
RemoteFlag = 0;
if whoiam
if Parallel(ThisMatlab).Local==0,
if Parallel(ThisMatlab).Local==0
RemoteFlag =1;
end
prct0={0,whoiam,Parallel(ThisMatlab)};
@ -165,7 +165,7 @@ while fpar<B
elseif info(1) == 5
errordef = 'Rank condition is not satisfied';
end
if strcmpi(type,'prior'),
if strcmpi(type,'prior')
disp(['PosteriorIRF :: Dynare is unable to solve the model (' errordef ')'])
continue
else
@ -240,9 +240,9 @@ while fpar<B
else
stock_irf_bvardsge(:,:,:,IRUN) = reshape(tmp_dsgevar,options_.irf,dataset_.vobs,M_.exo_nbr);
instr = [MhDirectoryName '/' M_.fname '_irf_bvardsge' ...
int2str(NumberOfIRFfiles_dsgevar) '.mat stock_irf_bvardsge;'];,
int2str(NumberOfIRFfiles_dsgevar) '.mat stock_irf_bvardsge;'];
eval(['save ' instr]);
if RemoteFlag==1,
if RemoteFlag==1
OutputFileName_bvardsge = [OutputFileName_bvardsge; {[MhDirectoryName filesep], [M_.fname '_irf_bvardsge' int2str(NumberOfIRFfiles_dsgevar) '.mat']}];
end
NumberOfIRFfiles_dsgevar = NumberOfIRFfiles_dsgevar+1;
@ -259,14 +259,14 @@ while fpar<B
int2str(NumberOfIRFfiles_dsgevar) '.mat stock_irf_bvardsge;'];
eval(['save ' instr]);
NumberOfIRFfiles_dsgevar = NumberOfIRFfiles_dsgevar+1;
if RemoteFlag==1,
if RemoteFlag==1
OutputFileName_bvardsge = [OutputFileName_bvardsge; {[MhDirectoryName filesep], [M_.fname '_irf_bvardsge' int2str(NumberOfIRFfiles_dsgevar) '.mat']}];
end
irun = 0;
end
end
save([MhDirectoryName '/' M_.fname '_irf_dsge' int2str(NumberOfIRFfiles_dsge) '.mat'],'stock_irf_dsge');
if RemoteFlag==1,
if RemoteFlag==1
OutputFileName_dsge = [OutputFileName_dsge; {[MhDirectoryName filesep], [M_.fname '_irf_dsge' int2str(NumberOfIRFfiles_dsge) '.mat']}];
end
NumberOfIRFfiles_dsge = NumberOfIRFfiles_dsge+1;
@ -278,7 +278,7 @@ while fpar<B
end
stock = stock_param;
save([MhDirectoryName '/' M_.fname '_param_irf' int2str(ifil2) '.mat'],'stock');
if RemoteFlag==1,
if RemoteFlag==1
OutputFileName_param = [OutputFileName_param; {[MhDirectoryName filesep], [M_.fname '_param_irf' int2str(ifil2) '.mat']}];
end
ifil2 = ifil2 + 1;

View File

@ -2,10 +2,10 @@ function myoutput=PosteriorIRF_core2(myinputs,fpar,npar,whoiam,ThisMatlab)
% function myoutput=PosteriorIRF_core2(myinputs,fpar,npar,whoiam, ThisMatlab)
% Generates the Posterior IRFs plot from the IRFs generated in
% PosteriorIRF_core1
%
%
% PARALLEL CONTEXT
% Performs in parallel execution a portion of the PosteriorIRF.m code.
% For more information, see the comment in posterior_sampler_core.m
% For more information, see the comment in posterior_sampler_core.m
% function.
%
% INPUTS
@ -49,7 +49,7 @@ function myoutput=PosteriorIRF_core2(myinputs,fpar,npar,whoiam,ThisMatlab)
global options_ M_
if nargin<4,
if nargin<4
whoiam=0;
end
@ -85,8 +85,8 @@ end
DirectoryName = CheckPath('Output',M_.dname);
RemoteFlag = 0;
if whoiam,
if Parallel(ThisMatlab).Local==0,
if whoiam
if Parallel(ThisMatlab).Local==0
RemoteFlag =1;
end
prct0={0,whoiam,Parallel(ThisMatlab)};
@ -96,7 +96,7 @@ end
OutputFileName={};
subplotnum = 0;
for i=fpar:npar,
for i=fpar:npar
figunumber = 0;
for j=1:nvar
@ -153,13 +153,13 @@ for i=fpar:npar,
if subplotnum == MaxNumberOfPlotPerFigure || (j == nvar && subplotnum> 0)
figunumber = figunumber+1;
dyn_saveas(hh,[DirectoryName '/' M_.fname '_Bayesian_IRF_' deblank(tit(i,:)) '_' int2str(figunumber)],options_.nodisplay,options_.graph_format);
if RemoteFlag==1,
if RemoteFlag==1
OutputFileName = [OutputFileName; {[DirectoryName,filesep], [M_.fname '_Bayesian_IRF_' deblank(tit(i,:)) '_' int2str(figunumber) '.*']}];
end
subplotnum = 0;
end
end% loop over selected endo_var
if whoiam,
if whoiam
fprintf('Done! \n');
waitbarString = [ 'Exog. shocks ' int2str(i) '/' int2str(npar) ' done.'];
% fMessageStatus((i-fpar+1)/(npar-fpar+1),whoiam,waitbarString, waitbarTitle, Parallel(ThisMatlab));

View File

@ -18,14 +18,14 @@ function ReshapeMatFiles(type, type2)
% posterior
% gsa
% prior
%
%
% OUTPUTS:
% none
% none
%
% SPECIAL REQUIREMENTS
% none
% Copyright (C) 2003-2011 Dynare Team
% Copyright (C) 2003-2017 Dynare Team
%
% This file is part of Dynare.
%
@ -44,15 +44,15 @@ function ReshapeMatFiles(type, type2)
global M_ options_
if nargin==1,
if nargin==1
MhDirectoryName = [ CheckPath('metropolis',M_.dname) filesep ];
else
if strcmpi(type2,'posterior')
MhDirectoryName = [CheckPath('metropolis',M_.dname) filesep ];
elseif strcmpi(type2,'gsa')
if options_.opt_gsa.morris==1,
if options_.opt_gsa.morris==1
MhDirectoryName = [CheckPath('gsa/screen',M_.dname) filesep ];
elseif options_.opt_gsa.morris==2,
elseif options_.opt_gsa.morris==2
MhDirectoryName = [CheckPath('gsa/identif',M_.dname) filesep ];
elseif options_.opt_gsa.pprior
MhDirectoryName = [CheckPath(['gsa' filesep 'prior'],M_.dname) filesep ];
@ -61,17 +61,17 @@ else
end
else
MhDirectoryName = [CheckPath('prior',M_.dname) filesep ];
end
end
end
switch type
case 'irf_dsge'
CAPtype = 'IRF_DSGE';
TYPEsize = [ options_.irf , size(options_.varlist,1) , M_.exo_nbr ];
TYPEarray = 4;
TYPEarray = 4;
case 'irf_bvardsge'
CAPtype = 'IRF_BVARDSGE';
TYPEsize = [ options_.irf , length(options_.varobs) , M_.exo_nbr ];
TYPEarray = 4;
TYPEarray = 4;
case 'smooth'
CAPtype = 'SMOOTH';
TYPEsize = [ M_.endo_nbr , options_.nobs ];
@ -134,7 +134,7 @@ switch TYPEarray
eval(['idx = idx + size(stock_' type ',4);'])
end
%eval(['STOCK_' CAPtype ' = sort(STOCK_' CAPtype ',4);'])
save([MhDirectoryName M_.fname '_' CAPtype 's' int2str(NumberOfTYPEfiles-foffset+1) '.mat'],['STOCK_' CAPtype]);
save([MhDirectoryName M_.fname '_' CAPtype 's' int2str(NumberOfTYPEfiles-foffset+1) '.mat'],['STOCK_' CAPtype]);
end
else
load([MhDirectoryName M_.fname '_' type '1.mat']);
@ -176,7 +176,7 @@ switch TYPEarray
load([MhDirectoryName M_.fname '_' type '1.mat']);
%eval(['STOCK_' CAPtype ' = sort(stock_' type ',3);'])
eval(['STOCK_' CAPtype ' = stock_' type ';'])
save([MhDirectoryName M_.fname '_' CAPtype 's' int2str(1) '.mat'],['STOCK_' CAPtype ]);
save([MhDirectoryName M_.fname '_' CAPtype 's' int2str(1) '.mat'],['STOCK_' CAPtype ]);
end
% Original file format may be useful in some cases...
% for file = 1:NumberOfTYPEfiles

View File

@ -3,15 +3,15 @@ function [fval,info,exit_flag,DLIK,Hess,SteadyState,trend_coeff] = TaRB_optimiz
% Wrapper function for target function used in TaRB algorithm; reassembles
% full parameter vector before calling target function
%
% INPUTS
% o optpar [double] (p_opt*1) vector of subset of parameters to be considered
% o par_vector [double] (p*1) full vector of parameters
% INPUTS
% o optpar [double] (p_opt*1) vector of subset of parameters to be considered
% o par_vector [double] (p*1) full vector of parameters
% o parameterindices [double] (p_opt*1) index of optpar entries in
% par_vector
% o TargetFun [char] string specifying the name of the objective
% function (posterior kernel).
% o varargin [structure] other inputs of target function
%
%
% OUTPUTS
% o fval [scalar] value of (minus) the likelihood.
% o info [double] (p*2) error code vector
@ -20,8 +20,9 @@ function [fval,info,exit_flag,DLIK,Hess,SteadyState,trend_coeff] = TaRB_optimiz
% o Hess [double] (p*p) asymptotic Hessian matrix.
% o SteadyState [double] Vector of doubles, steady state level for the endogenous variables.
% o trend_coeff [double] Matrix of doubles, coefficients of the deterministic trend in the measurement equation
%
% Copyright (C) 2015-16 Dynare Team
%
% Copyright (C) 2015-2017 Dynare Team
%
% This file is part of Dynare.
%
@ -40,4 +41,3 @@ function [fval,info,exit_flag,DLIK,Hess,SteadyState,trend_coeff] = TaRB_optimiz
par_vector(parameterindices,:)=optpar; %reassemble parameter
[fval,info,exit_flag,DLIK,Hess,SteadyState,trend_coeff] = feval(TargetFun,par_vector,varargin{:}); %call target function

View File

@ -2,19 +2,19 @@ function [] = Tracing()
% DESCRIPTION
% This function is used to test the correct execution of a matlab section
% on remote machine.
%
%
% If no error happen the function simply create a file.
%
% INPUTS
% ...
%
%
% OUTPUTS
% ...
%
%
% SPECIAL REQUIREMENTS
% none
% Copyright (C) 2010 Dynare Team
% Copyright (C) 2010-2017 Dynare Team
%
% This file is part of Dynare.
%

View File

@ -1,15 +1,15 @@
function [oo_] = UnivariateSpectralDensity(M_,oo_,options_,var_list)
% This function computes the theoretical spectral density of each
% endogenous variable declared in var_list. Results are stored in
% oo_.SpectralDensity and may be plotted. Plots are saved into the
% graphs-folder.
%
% endogenous variable declared in var_list. Results are stored in
% oo_.SpectralDensity and may be plotted. Plots are saved into the
% graphs-folder.
%
% INPUTS
% M_ [structure] Dynare's model structure
% oo_ [structure] Dynare's results structure
% options_ [structure] Dynare's options structure
% var_list [integer] Vector of indices for a subset of variables.
%
%
% OUTPUTS
% oo_ [structure] Dynare's results structure,
% containing the subfield
@ -17,7 +17,7 @@ function [oo_] = UnivariateSpectralDensity(M_,oo_,options_,var_list)
% and density, which are of size nvar*ngrid.
%
% Adapted from th_autocovariances.m.
% Adapted from th_autocovariances.m.
% Copyright (C) 2006-2017 Dynare Team
%
@ -38,7 +38,7 @@ function [oo_] = UnivariateSpectralDensity(M_,oo_,options_,var_list)
if options_.order > 1
disp('UnivariateSpectralDensity :: I Cannot compute the theoretical spectral density')
disp('UnivariateSpectralDensity :: I Cannot compute the theoretical spectral density')
disp('with a second order approximation of the DSGE model!')
disp('Please set order = 1. I abort')
return
@ -102,7 +102,7 @@ if ~isempty(u)
ivar = oo_.dr.order_var(iky);
end
iky = iv(ivar);
iky = iv(ivar);
aa = ghx(iky,:);
bb = ghu(iky,:);
ngrid = options_.hp_ngrid; %number of grid points
@ -122,7 +122,7 @@ elseif ~(options_.hp_filter == 0 && ~options_.bandpass.indicator) && options_.ba
filter_gain(freqs<=-2*pi/lowest_periodicity+2*pi & freqs>=-2*pi/highest_periodicity+2*pi)=1;
elseif ~(options_.hp_filter == 0 && ~options_.bandpass.indicator) && ~options_.bandpass.indicator %filter with HP-filter
lambda = options_.hp_filter;
filter_gain = 4*lambda*(1 - cos(freqs)).^2 ./ (1 + 4*lambda*(1 - cos(freqs)).^2);
filter_gain = 4*lambda*(1 - cos(freqs)).^2 ./ (1 + 4*lambda*(1 - cos(freqs)).^2);
end
mathp_col = NaN(ngrid,length(ivar)^2);
@ -134,12 +134,12 @@ for ig = 1:ngrid
g_omega = [aa*tneg(ig) bb]*f_omega*[aa'*tpos(ig); bb']; % selected variables
f_hp = filter_gain(ig)^2*g_omega; % spectral density of selected filtered series
mathp_col(ig,:) = (f_hp(:))'; % store as matrix row
end;
end
f = zeros(nvar,ngrid);
for i=1:nvar
f(i,:) = real(mathp_col(:,(i-1)*nvar+i)); %read out spectral density
end
end
oo_.SpectralDensity.freqs=freqs;
oo_.SpectralDensity.density=f;
@ -164,7 +164,7 @@ if options_.nograph == 0
xlabel('0 \leq \omega \leq \pi')
ylabel('f(\omega)')
box on
axis tight
axis tight
dyn_saveas(hh,[M_.fname ,filesep,'graphs', filesep, 'SpectralDensity_' deblank(M_.endo_names(ivar(i),:))],options_.nodisplay,options_.graph_format)
end
end

View File

@ -1,7 +1,7 @@
function WriteShockDecomp2Excel(z,shock_names,endo_names,i_var,initial_date,DynareModel,DynareOptions,opts_decomp)
%function WriteShockDecomp2Excel(z,shock_names,endo_names,i_var,initial_date,DynareModel,DynareOptions)
% Saves the results from the shock_decomposition command to xls
%
%
% Inputs
% z [n_var*(nshock+2)*nperiods] shock decomposition array, see shock_decomposition.m for details
% shock_names [endo_nbr*string length] shock names from M_.exo_names
@ -42,14 +42,14 @@ end
% number of components equals number of shocks + 1 (initial conditions)
comp_nbr = size(z,2)-1;
if nargin==8 ,
if nargin==8
if isfield(opts_decomp,'steady_state')
SteadyState = opts_decomp.steady_state;
end
if isfield(opts_decomp,'fig_mode') && ~isempty(opts_decomp.fig_mode)
fig_mode = opts_decomp.fig_mode;
fig_mode1 = ['_' fig_mode];
fig_mode = [fig_mode '_'];
fig_mode1 = ['_' fig_mode];
fig_mode = [fig_mode '_'];
end
if isfield(opts_decomp,'screen_shocks')
if use_shock_groups
@ -61,14 +61,14 @@ if nargin==8 ,
if isfield(opts_decomp,'fig_name')
fig_name = opts_decomp.fig_name;
% fig_name = ['_' fig_name];
fig_name1 = [fig_name];
fig_name = [fig_name '_'];
fig_name1 = [fig_name];
fig_name = [fig_name '_'];
end
if screen_shocks
fig_name1 = [fig_name1 '_screen'];
fig_name = [fig_name 'screen_'];
fig_name1 = [fig_name1 '_screen'];
fig_name = [fig_name 'screen_'];
end
end
end
gend = size(z,3);
@ -84,7 +84,7 @@ end
nvar = length(i_var);
labels = char(char(shock_names),'Initial values');
if ~(screen_shocks && comp_nbr>18),
if ~(screen_shocks && comp_nbr>18)
screen_shocks=0;
end
comp_nbr0=comp_nbr;
@ -92,14 +92,14 @@ comp_nbr0=comp_nbr;
for j=1:nvar
d0={};
z1 = squeeze(z(i_var(j),:,:));
if screen_shocks,
if screen_shocks
[junk, isort] = sort(mean(abs(z1(1:end-2,:)')), 'descend');
labels = char(char(shock_names(isort(1:16),:)),'Others', 'Initial values');
zres = sum(z1(isort(17:end),:),1);
z1 = [z1(isort(1:16),:); zres; z1(comp_nbr0:end,:)];
comp_nbr=18;
end
d0(1,:)=[{'Decomposition'} cellstr(labels(1:comp_nbr,:))' {'Smoot Var'}];
d0=[d0; num2cell([x' z1'])];
LastRow=size(d0,1);
@ -107,12 +107,12 @@ for j=1:nvar
d0(LastRow+2,1)={'Legend.'};
d0(LastRow+2,2)={'Shocks include:'};
d0(LastRow+3:LastRow+3+comp_nbr-1,1)=cellstr(labels(1:comp_nbr,:));
for ic=1:comp_nbr,
for ic=1:comp_nbr
group_members = shock_groups.(shock_ind{ic}).shocks;
d0(LastRow+2+ic,2:1+length(group_members))=group_members;
end
end
warning off
if ~ismac
[STATUS,MESSAGE] = xlswrite([DynareModel.fname,'_shock_decomposition',fig_mode,fig_name1],d0,deblank(endo_names(i_var(j),:)));
@ -122,6 +122,5 @@ for j=1:nvar
warning on
clear d0
end
end

View File

@ -1,5 +1,22 @@
function title=add_filter_subtitle(title,options_)
% Copyright (C) 2015-2017 Dynare Team
%
% This file is part of Dynare.
%
% Dynare is free software: you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% Dynare is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with Dynare. If not, see <http://www.gnu.org/licenses/>.
if ~options_.hp_filter && ~options_.one_sided_hp_filter && ~options_.bandpass.indicator %do not filter
%nothing to add here
elseif ~options_.hp_filter && ~options_.one_sided_hp_filter && options_.bandpass.indicator
@ -10,6 +27,6 @@ elseif options_.hp_filter && ~options_.one_sided_hp_filter && ~options_.bandpas
num2str(options_.hp_filter) ')'];
elseif ~options_.hp_filter && options_.one_sided_hp_filter && ~options_.bandpass.indicator
title = [title ' (One-sided HP filter, lambda = ' ...
num2str(options_.one_sided_hp_filter) ')'];
num2str(options_.one_sided_hp_filter) ')'];
end
end

View File

@ -1,28 +1,28 @@
function [z, endo_names, endo_names_tex, steady_state, i_var, oo_] = annualized_shock_decomposition(oo_, M_, options_, i_var, t0, t1, realtime_, vintage_, steady_state, q2a, cumfix)
% function oo_ = annualized_shock_decomposition(oo_,t0,options_.nobs);
% Computes annualized shocks contribution to a simulated trajectory. The fields set are
% oo_.annualized_shock_decomposition, oo_.annualized_realtime_shock_decomposition,
% oo_.annualized_realtime_conditional_shock_decomposition and oo_.annualized_realtime_forecast_shock_decomposition.
% oo_.annualized_shock_decomposition, oo_.annualized_realtime_shock_decomposition,
% oo_.annualized_realtime_conditional_shock_decomposition and oo_.annualized_realtime_forecast_shock_decomposition.
% Subfields are arrays n_var by nshock+2 by nperiods. The
% first nshock columns store the respective shock contributions, column n+1
% stores the role of the initial conditions, while column n+2 stores the
% value of the smoothed variables. Both the variables and shocks are stored
% value of the smoothed variables. Both the variables and shocks are stored
% in the order of endo_names and M_.exo_names, respectively.
%
% INPUTS
% oo_: [structure] Storage of results
% M_: [structure] Storage of model
% opts: [structure] options for shock decomp
% opts: [structure] options for shock decomp
% i_var: [array] index of vars
% t0: [integer] first period
% t1: [integer] last period
% realtime_: [integer]
% realtime_: [integer]
% vintage_: [integer]
% steady_state: [array] steady state value of quarterly (log-) level vars
% q2a: [structure] info on q2a
%
% OUTPUTS
% z: [matrix] shock decomp to plot
% z: [matrix] shock decomp to plot
% endo_names: [char] updated var names
% endo_names_tex: [char] updated TeX var names
% steady_state: [array] updated stady state of vars
@ -57,7 +57,7 @@ islog = q2a.islog;
aux = q2a.aux;
aux0 = aux;
cumfix = q2a.cumfix;
% usual shock decomp
% usual shock decomp
if isstruct(oo_)
% z = oo_.shock_decomposition;
myopts=options_;
@ -75,7 +75,7 @@ if isfield(q2a,'name')
if isfield(q2a,'tex_name')
mytex = q2a.tex_name;
end
if mytype==2,
if mytype==2
gtxt = ['PHI' mytxt]; % inflation rate
gtex = ['{\pi(' mytex ')}'];
elseif mytype
@ -101,7 +101,7 @@ if isstruct(aux)
end
yaux=aux.y;
end
if mytype==2,
if mytype==2
gtxt = 'PHI'; % inflation rate
gtex = '\pi';
elseif mytype
@ -115,7 +115,7 @@ nterms = size(z,2);
nfrcst = opts.forecast/4;
for j=1:nvar
if j>1,
if j>1
endo_names = char(endo_names,[deblank(M_.endo_names(i_var(j),:)) '_A']);
endo_names_tex = char(endo_names_tex,['{' deblank(M_.endo_names_tex(i_var(j),:)) '}^A']);
gendo_names = char(gendo_names,[gtxt endo_names(j,:)]);
@ -133,15 +133,15 @@ for j=1:nvar
gendo_names_tex = [gtex '(' deblank(endo_names_tex(j,:)) ')'];
end
end
for k =1:nterms,
if isstruct(aux),
for k =1:nterms
if isstruct(aux)
aux.y = squeeze(yaux(j,k,min((t0-3):-4:1):end));
end
[za(j,k,:), steady_state_a(j,1), gza(j,k,:), steady_state_ga(j,1)] = ...
quarterly2annual(squeeze(z(j,k,min((t0-3):-4:1):end)),steady_state(j),GYTREND0,var_type,islog,aux);
end
ztmp=squeeze(za(j,:,:));
if cumfix==0,
if cumfix==0
zscale = sum(ztmp(1:end-1,:))./ztmp(end,:);
ztmp(1:end-1,:) = ztmp(1:end-1,:)./repmat(zscale,[nterms-1,1]);
else
@ -149,7 +149,7 @@ for j=1:nvar
ztmp(end-1,:) = ztmp(end-1,:) + zres;
end
gztmp=squeeze(gza(j,:,:));
if cumfix==0,
if cumfix==0
gscale = sum(gztmp(1:end-1,:))./ gztmp(end,:);
gztmp(1:end-1,:) = gztmp(1:end-1,:)./repmat(gscale,[nterms-1,1]);
else
@ -160,7 +160,7 @@ for j=1:nvar
gza(j,:,:) = gztmp;
end
if q2a.plot ==1,
if q2a.plot ==1
z=gza;
endo_names = gendo_names;
endo_names_tex = gendo_names_tex;
@ -176,9 +176,9 @@ end
% end
% realtime
if realtime_ && isstruct(oo_) && isfield(oo_, 'realtime_shock_decomposition'),
if realtime_ && isstruct(oo_) && isfield(oo_, 'realtime_shock_decomposition')
init=1;
for i=t0:4:t1,
for i=t0:4:t1
yr=floor(i/4);
za=[];
gza=[];
@ -197,18 +197,18 @@ if realtime_ && isstruct(oo_) && isfield(oo_, 'realtime_shock_decomposition'),
yaux=aux.y;
end
nterms = size(z,2);
% z = oo_.realtime_shock_decomposition.(['time_' int2str(i)]);
% z = z(i_var,:,:);
for j=1:nvar
for k =nterms:-1:1,
for k =nterms:-1:1
% if k<nterms
% ztmp = squeeze(sum(z(j,[1:k-1,k+1:end-1],t0-4:end)));
% else
ztmp = squeeze(z(j,k,min((t0-3):-4:1):end));
% end
if isstruct(aux),
if isstruct(aux)
aux.y = squeeze(yaux(j,k,min((t0-3):-4:1):end));
end
[za(j,k,:), steady_state_a(j,1), gza(j,k,:), steady_state_ga(j,1)] = ...
@ -217,47 +217,47 @@ if realtime_ && isstruct(oo_) && isfield(oo_, 'realtime_shock_decomposition'),
% za(j,k,:) = za(j,end,:) - za(j,k,:);
% gza(j,k,:) = gza(j,end,:) - gza(j,k,:);
% end
end
ztmp=squeeze(za(j,:,:));
if cumfix==0,
if cumfix==0
zscale = sum(ztmp(1:end-1,:))./ztmp(end,:);
ztmp(1:end-1,:) = ztmp(1:end-1,:)./repmat(zscale,[nterms-1,1]);
else
zres = ztmp(end,:)-sum(ztmp(1:end-1,:));
ztmp(end-1,:) = ztmp(end-1,:) + zres;
end
gztmp=squeeze(gza(j,:,:));
if cumfix==0,
if cumfix==0
gscale = sum(gztmp(1:end-1,:))./ gztmp(end,:);
gztmp(1:end-1,:) = gztmp(1:end-1,:)./repmat(gscale,[nterms-1,1]);
else
gres = gztmp(end,:) - sum(gztmp(1:end-1,:));
gztmp(end-1,:) = gztmp(end-1,:)+gres;
end
za(j,:,:) = ztmp;
gza(j,:,:) = gztmp;
end
if q2a.plot ==1,
if q2a.plot ==1
z=gza;
elseif q2a.plot == 2
z=za;
else
z=cat(1,za,gza);
end
if init==1,
if init==1
oo_.annualized_realtime_shock_decomposition.pool = z;
else
oo_.annualized_realtime_shock_decomposition.pool(:,:,yr) = z(:,:,end-nfrcst);
end
end
oo_.annualized_realtime_shock_decomposition.(['yr_' int2str(yr)]) = z;
if opts.forecast
oo_.annualized_realtime_forecast_shock_decomposition.(['yr_' int2str(yr)]) = z(:,:,end-nfrcst:end);
if init>nfrcst
@ -272,7 +272,7 @@ if realtime_ && isstruct(oo_) && isfield(oo_, 'realtime_shock_decomposition'),
oo_.annualized_realtime_conditional_shock_decomposition.(['yr_' int2str(yr-nfrcst)])(:,end,:) = ...
oo_.annualized_realtime_shock_decomposition.pool(:,end,yr-nfrcst:end);
if i==t1
for my_forecast_=(nfrcst-1):-1:1,
for my_forecast_=(nfrcst-1):-1:1
oo_.annualized_realtime_conditional_shock_decomposition.(['yr_' int2str(yr-my_forecast_)]) = ...
oo_.annualized_realtime_shock_decomposition.pool(:,:,yr-my_forecast_:yr) - ...
oo_.annualized_realtime_forecast_shock_decomposition.(['yr_' int2str(yr-my_forecast_)])(:,:,1:my_forecast_+1);
@ -293,24 +293,24 @@ if realtime_ && isstruct(oo_) && isfield(oo_, 'realtime_shock_decomposition'),
switch realtime_
case 0
z = oo_.annualized_shock_decomposition;
case 1 % realtime
if vintage_
z = oo_.annualized_realtime_shock_decomposition.(['yr_' int2str(floor(vintage_/4))]);
else
z = oo_.annualized_realtime_shock_decomposition.pool;
end
case 2 % conditional
if vintage_
z = oo_.annualized_realtime_conditional_shock_decomposition.(['yr_' int2str(floor(vintage_/4))]);
else
error();
end
case 3 % forecast
if vintage_
z = oo_.annualized_realtime_forecast_shock_decomposition.(['yr_' int2str(floor(vintage_/4))]);
@ -320,15 +320,14 @@ if realtime_ && isstruct(oo_) && isfield(oo_, 'realtime_shock_decomposition'),
end
end
if q2a.plot ==0,
if q2a.plot ==0
i_var=1:2*nvar;
steady_state = [steady_state_a;steady_state_ga];
else
i_var=1:nvar;
if q2a.plot ==1,
if q2a.plot ==1
steady_state = steady_state_ga;
else
steady_state = steady_state_a;
end
end

View File

@ -2,48 +2,48 @@ function [InnovationVariance,AutoregressiveParameters] = autoregressive_process_
% This function computes the parameters of an AR(p) process from the variance and the autocorrelation function
% (the first p terms) of this process.
%
% INPUTS
% INPUTS
% [1] Variance [double] scalar, variance of the variable.
% [2] Rho [double] p*1 vector, the autocorelation function: \rho(1), \rho(2), ..., \rho(p).
% [3] p [double] scalar, the number of lags in the AR process.
%
% OUTPUTS
% OUTPUTS
% [1] InnovationVariance [double] scalar, the variance of the innovation.
% [2] AutoregressiveParameters [double] p*1 vector of autoregressive parameters.
%
% NOTES
% NOTES
%
% The AR(p) model for {y_t} is:
%
% y_t = \phi_1 * y_{t-1} + \phi_2 * y_{t-2} + ... + \phi_p * y_{t-p} + e_t
%
% y_t = \phi_1 * y_{t-1} + \phi_2 * y_{t-2} + ... + \phi_p * y_{t-p} + e_t
%
% Let \gamma(0) and \rho(1), ..., \rho(2) be the variance and the autocorrelation function of {y_t}. This function
% compute the variance of {e_t} and the \phi_i (i=1,...,p) from the variance and the autocorrelation function of {y_t}.
% compute the variance of {e_t} and the \phi_i (i=1,...,p) from the variance and the autocorrelation function of {y_t}.
% We know that:
%
%
% \gamma(0) = \phi_1 \gamma(1) + ... + \phi_p \gamma(p) + \sigma^2
%
% where \sigma^2 is the variance of {e_t}. Equivalently we have:
%
% \sigma^2 = \gamma(0) (1-\rho(1)\phi_1 - ... - \rho(p)\phi_p)
% \sigma^2 = \gamma(0) (1-\rho(1)\phi_1 - ... - \rho(p)\phi_p)
%
% We also have for any integer h>0:
%
%
% \rho(h) = \phi_1 \rho(h-1) + ... + \phi_p \rho(h-p)
%
% We can write the equations for \rho(1), ..., \rho(p) using matrices. Let R be the p*p autocorelation
% matrix and v be the p*1 vector gathering the first p terms of the autocorrelation function. We have:
% matrix and v be the p*1 vector gathering the first p terms of the autocorrelation function. We have:
%
% v = R*PHI
%
%
% where PHI is a p*1 vector with the autoregressive parameters of the AR(p) process. We can recover the autoregressive
% parameters by inverting the autocorrelation matrix: PHI = inv(R)*v.
%
%
% This function first computes the vector PHI by inverting R and computes the variance of the innovation by evaluating
%
% \sigma^2 = \gamma(0)*(1-PHI'*v)
% Copyright (C) 2009 Dynare Team
% Copyright (C) 2009-2017 Dynare Team
%
% This file is part of Dynare.
%

View File

@ -1,6 +1,6 @@
function [endogenousvariables, exogenousvariables] = backward_model_inversion(constraints, exogenousvariables, initialconditions, endo_names, exo_names, freeinnovations, DynareModel, DynareOptions, DynareOutput)
% INPUTS
% INPUTS
% - constraints [dseries] with N constrained endogenous variables from t1 to t2.
% - exogenousvariables [dseries] with Q exogenous variables.
% - initialconditions [dseries] with M endogenous variables starting before t1 (M initialcond must contain at least the state variables).
@ -8,11 +8,11 @@ function [endogenousvariables, exogenousvariables] = backward_model_inversion(co
% - exo_names [cell] list of exogenous variable names.
% - freeinstruments [cell] list of exogenous variable names used to control the constrained endogenous variables.
%
% OUTPUTS
% OUTPUTS
% - endogenous [dseries]
% - exogenous [dseries]
%
% REMARKS
% REMARKS
% Copyright (C) 2017 Dynare Team
%

View File

@ -34,7 +34,7 @@ function DynareOutput = simul_backward_linear_model(initial_conditions, sample_s
%! @end deftypefn
%@eod:
% Copyright (C) 2012-2016 Dynare Team
% Copyright (C) 2012-2017 Dynare Team
%
% This file is part of Dynare.
%
@ -84,9 +84,9 @@ A0inv = inv(jacob(:,jdx));
A1 = jacob(:,nonzeros(DynareModel.lead_lag_incidence(1,:)));
B = jacob(:,end-number_of_shocks+1:end);
% Simulations
% Simulations
for it = 2:sample_size+1
Y(:,it) = -A0inv*(cst + A1*Y(iy1,it-1) + B*DynareOutput.exo_simul(it,:)');
Y(:,it) = -A0inv*(cst + A1*Y(iy1,it-1) + B*DynareOutput.exo_simul(it,:)');
end
DynareOutput.endo_simul = Y;

View File

@ -34,7 +34,7 @@ function DynareOutput = simul_backward_model(initial_conditions, sample_size, Dy
%! @end deftypefn
%@eod:
% Copyright (C) 2012-2016 Dynare Team
% Copyright (C) 2012-2017 Dynare Team
%
% This file is part of Dynare.
%

View File

@ -2,7 +2,7 @@ function DynareOutput = simul_backward_nonlinear_model(initial_conditions, sampl
% Simulates a stochastic non linear backward looking model with arbitrary precision (a deterministic solver is used).
%
% INPUTS
% INPUTS
% - initial_conditions [double] n*1 vector, initial conditions for the endogenous variables.
% - sample_size [integer] scalar, number of periods for the simulation.
% - DynareOptions [struct] Dynare's options_ global structure.
@ -10,10 +10,10 @@ function DynareOutput = simul_backward_nonlinear_model(initial_conditions, sampl
% - DynareOutput [struct] Dynare's oo_ global structure.
% - innovations [double] T*q matrix, innovations to be used for the simulation.
%
% OUTPUTS
% OUTPUTS
% - DynareOutput [struct] Dynare's oo_ global structure.
%
% REMARKS
% REMARKS
% [1] The innovations used for the simulation are saved in DynareOutput.exo_simul, and the resulting paths for the endogenous
% variables are saved in DynareOutput.endo_simul.
% [2] The last input argument is not mandatory. If absent we use random draws and rescale them with the informations provided

View File

@ -14,7 +14,7 @@ function plan = basic_plan(plan, exogenous, expectation_type, date, value)
% plan [structure] Returns a structure containing the updated forecast scenario.
%
%
% Copyright (C) 2013-2014 Dynare Team
% Copyright (C) 2013-2017 Dynare Team
%
% This file is part of Dynare.
%
@ -38,7 +38,7 @@ exogenous = strtrim(exogenous);
ix = find(strcmp(exogenous, plan.exo_names));
if isempty(ix)
error(['in basic_plan the second argument ' exogenous ' is not an exogenous variable']);
end;
end
sdate = length(date);
if sdate > 1
if date(1) < plan.date(1) || date(end) > plan.date(end)

View File

@ -4,14 +4,14 @@ function d = bksup0(c,ny,jcf,iyf,icf,periods)
% INPUTS
% ny: number of endogenous variables
% jcf: variables index forward
%
%
% OUTPUTS
% d: vector of backsubstitution results
%
% SPECIAL REQUIREMENTS
% none
% Copyright (C) 2003-2009 Dynare Team
% Copyright (C) 2003-2017 Dynare Team
%
% This file is part of Dynare.
%

View File

@ -5,14 +5,14 @@ function d = bksup1(c,ny,jcf,iyf,periods)
% INPUTS
% ny: number of endogenous variables
% jcf: variables index forward
%
%
% OUTPUTS
% d: vector of backsubstitution results
%
% SPECIAL REQUIREMENTS
% none
% Copyright (C) 2003-2010 Dynare Team
% Copyright (C) 2003-2017 Dynare Team
%
% This file is part of Dynare.
%
@ -40,4 +40,3 @@ for i = 2:periods
end
d = c(:,jcf) ;

View File

@ -15,7 +15,7 @@ function d1 = bksupk(ny,fid,jcf,icc1)
% SPECIAL REQUIREMENTS
% none
% Copyright (C) 2003-2011 Dynare Team
% Copyright (C) 2003-2017 Dynare Team
%
% This file is part of Dynare.
%
@ -68,7 +68,7 @@ while i <= options_.periods
c = fread(fid,[jcf,ny],'float64')' ;
d1(ir) = c(:,jcf)-c(:,icf)*d1(irf) ;
ir = ir-ny ;
ir = ir-ny ;
irf = irf-ny ;
i = i+1;
end

View File

@ -1,6 +1,6 @@
function x = bseastr(s1,s2)
% Copyright (C) 2001-2009 Dynare Team
% Copyright (C) 2001-2017 Dynare Team
%
% This file is part of Dynare.
%
@ -33,16 +33,15 @@ for im = 1:m
for i = 1:min(length(key),length(temp))
if temp(i) > key(i)
h = mid - 1 ;
break
break
else
l = mid + 1 ;
break
break
end
end
else
x(im) = mid ;
break
break
end
end
end

View File

@ -12,7 +12,7 @@ function bvar_density(maxnlags)
% none
% Copyright (C) 2003-2007 Christopher Sims
% Copyright (C) 2007-2016 Dynare Team
% Copyright (C) 2007-2017 Dynare Team
%
% This file is part of Dynare.
%
@ -37,16 +37,16 @@ for nlags = 1:maxnlags
[ny, nx, posterior, prior] = bvar_toolbox(nlags);
oo_.bvar.posterior{nlags}=posterior;
oo_.bvar.prior{nlags}=prior;
posterior_int = matrictint(posterior.S, posterior.df, posterior.XXi);
prior_int = matrictint(prior.S, prior.df, prior.XXi);
lik_nobs = posterior.df - prior.df;
log_dnsty = posterior_int - prior_int - 0.5*ny*lik_nobs*log(2*pi);
oo_.bvar.log_marginal_data_density(nlags)=log_dnsty;
skipline()
fprintf('The marginal log density of the BVAR(%g) model is equal to %10.4f\n', ...
nlags, log_dnsty);
@ -61,7 +61,7 @@ function w = matrictint(S, df, XXi)
% S: parameter of inverse-Wishart distribution
% df: number of degrees of freedom of inverse-Wishart distribution
% XXi: first component of VCV matrix of matrix-normal distribution
%
%
% Computes the integral over (Phi, Sigma) of:
%
% det(Sigma)^(-k/2)*exp(-0.5*Tr((Phi-PhiHat)'*(XXi)^(-1)*(Phi-PhiHat)*Sigma^(-1)))*

View File

@ -54,7 +54,7 @@ p = 0;
% Loop counter initialization
d = 0;
while d <= options_.bvar_replic
Sigma = rand_inverse_wishart(ny, posterior.df, S_inv_upper_chol);
% Option 'lower' of chol() not available in old versions of
@ -62,15 +62,15 @@ while d <= options_.bvar_replic
Sigma_lower_chol = chol(Sigma)';
Phi = rand_matrix_normal(k, ny, posterior.PhiHat, Sigma_lower_chol, XXi_lower_chol);
% All the eigenvalues of the companion matrix have to be on or inside the unit circle
Companion_matrix(1:ny,:) = Phi(1:ny*nlags,:)';
Companion_matrix(1:ny,:) = Phi(1:ny*nlags,:)';
test = (abs(eig(Companion_matrix)));
if any(test>1.0000000000001)
p = p+1;
end
d = d+1;
% Without shocks
lags_data = forecast_data.initval;
for t = 1:options_.forecast
@ -80,7 +80,7 @@ while d <= options_.bvar_replic
lags_data(end,:) = y;
sims_no_shock(t, :, d) = y;
end
% With shocks
lags_data = forecast_data.initval;
for t = 1:options_.forecast
@ -130,9 +130,9 @@ dyn_saveas(dyn_graph.fh,[OutputDirectoryName '/' M_.fname '_BVAR_forecast_',num2
% Compute RMSE
if ~isempty(forecast_data.realized_val)
sq_err_cumul = zeros(1, ny);
lags_data = forecast_data.initval;
for t = 1:size(forecast_data.realized_val, 1)
X = [ reshape(flipdim(lags_data, 1)', 1, ny*nlags) forecast_data.realized_xdata(t, :) ];
@ -141,14 +141,14 @@ if ~isempty(forecast_data.realized_val)
lags_data(end,:) = y;
sq_err_cumul = sq_err_cumul + (y - forecast_data.realized_val(t, :)) .^ 2;
end
rmse = sqrt(sq_err_cumul / size(forecast_data.realized_val, 1));
fprintf('RMSE of BVAR(%d):\n', nlags);
for i = 1:length(options_.varobs)
fprintf('%s: %10.4f\n', options_.varobs{i}, rmse(i));
end
end
end
% Store results

View File

@ -11,7 +11,7 @@ function bvar_irf(nlags,identification)
% SPECIAL REQUIREMENTS
% none
% Copyright (C) 2007-2012 Dynare Team
% Copyright (C) 2007-2017 Dynare Team
%
% This file is part of Dynare.
%
@ -54,7 +54,7 @@ p = 0;
sampled_irfs = NaN(ny, ny, options_.irf, options_.bvar_replic);
for draw=1:options_.bvar_replic
% Get a covariance matrix from an inverted Wishart distribution.
Sigma = rand_inverse_wishart(ny, posterior.df, S_inv_upper_chol);
Sigma_upper_chol = chol(Sigma);
@ -62,10 +62,10 @@ for draw=1:options_.bvar_replic
% Get the Autoregressive matrices from a matrix variate distribution.
Phi = rand_matrix_normal(k, ny, posterior.PhiHat, Sigma_lower_chol, XXi_lower_chol);
% Form the companion matrix.
Companion_matrix(1:ny,:) = transpose(Phi(1:ny*nlags,:));
Companion_matrix(1:ny,:) = transpose(Phi(1:ny*nlags,:));
% All the eigenvalues of the companion matrix have to be on or
% inside the unit circle to rule out explosive time series.
test = (abs(eig(Companion_matrix)));
@ -78,7 +78,7 @@ for draw=1:options_.bvar_replic
elseif strcmpi(identification,'SquareRoot')
StructuralMat = sqrtm(Sigma);
end
% Build the IRFs...
lags_data = zeros(ny,ny*nlags) ;
sampled_irfs(:,:,1,draw) = Sigma_lower_chol ;
@ -88,7 +88,7 @@ for draw=1:options_.bvar_replic
lags_data(:,ny+1:end) = lags_data(:,1:end-ny) ;
lags_data(:,1:ny) = sampled_irfs(:,:,t,draw) ;
end
end
if p > 0
@ -106,7 +106,7 @@ sort_idx = round((0.5 + [-options_.bvar.conf_sig, options_.bvar.conf_sig, .0]/2)
posterior_down_conf_irfs = sorted_irfs(:,:,:,sort_idx(1));
posterior_up_conf_irfs = sorted_irfs(:,:,:,sort_idx(2));
posterior_median_irfs = sorted_irfs(:,:,:,sort_idx(3));
posterior_median_irfs = sorted_irfs(:,:,:,sort_idx(3));
number_of_columns = fix(sqrt(ny));
number_of_rows = ceil(ny / number_of_columns) ;

View File

@ -29,7 +29,7 @@ function [ny, nx, posterior, prior, forecast_data] = bvar_toolbox(nlags)
% forecasting, of size options_.forecast*nx (actually only
% contains "1" values for the constant term if nx ~= 0)
% - realized_val: only non-empty if options_.nobs doesn't point
% to the end of sample
% to the end of sample
% In that case, contains values of endogenous variables after
% options_.nobs and up to the end of the sample
% - realized_xdata: contains values of exogenous variables after
@ -42,7 +42,7 @@ function [ny, nx, posterior, prior, forecast_data] = bvar_toolbox(nlags)
% - bvar_prior_{tau,decay,lambda,mu,omega,flat,train}
% Copyright (C) 2003-2007 Christopher Sims
% Copyright (C) 2007-2012 Dynare Team
% Copyright (C) 2007-2017 Dynare Team
%
% This file is part of Dynare.
%
@ -67,7 +67,7 @@ options_ = set_default_option(options_, 'nobs', size(dataset,1)-options_.first_o
if (options_.first_obs+options_.nobs-1)> size(dataset,1)
fprintf('Incorrect or missing specification of the number of observations. nobs can be at most %4u\n',size(dataset,1)-options_.first_obs+1);
error('Inconsistent number of observations.')
error('Inconsistent number of observations.')
end
% Parameters for prior
@ -160,7 +160,7 @@ function [ydum,xdum,breaks]=varprior(nv,nx,lags,mnprior,vprior)
%function [ydum,xdum,breaks]=varprior(nv,nx,lags,mnprior,vprior)
% ydum, xdum: dummy observation data that implement the prior
% breaks: vector of points in the dummy data after which new dummy obs's start
% Set breaks=T+[0;breaks], ydata=[ydata;ydum], xdum=[xdata;xdum], where
% Set breaks=T+[0;breaks], ydata=[ydata;ydum], xdum=[xdata;xdum], where
% actual data matrix has T rows, in preparing input for rfvar3
% nv,nx,lags: VAR dimensions
% mnprior.tight:Overall tightness of Minnesota prior
@ -175,8 +175,8 @@ function [ydum,xdum,breaks]=varprior(nv,nx,lags,mnprior,vprior)
% taken to include the sum-of-coefficients and co-persistence components
% that are implemented directly in rfvar3.m. The diagonal prior on v, combined
% with sum-of-coefficients and co-persistence components and with the unit own-first-lag
% prior mean generates larger prior variances for own than for cross-effects even in
% this formulation, but here there is no way to shrink toward a set of unconstrained
% prior mean generates larger prior variances for own than for cross-effects even in
% this formulation, but here there is no way to shrink toward a set of unconstrained
% univariate AR's.
% Original file downloaded from:
@ -216,3 +216,108 @@ dimy = size(ydum);
ydum = reshape(permute(ydum,[1 3 2]),dimy(1)*dimy(3),nv);
xdum = reshape(permute(xdum,[1 3 2]),dimy(1)*dimy(3),nx);
breaks = breaks(1:(end-1));
function var=rfvar3(ydata,lags,xdata,breaks,lambda,mu)
%function var=rfvar3(ydata,lags,xdata,breaks,lambda,mu)
% This algorithm goes for accuracy without worrying about memory requirements.
% ydata: dependent variable data matrix
% xdata: exogenous variable data matrix
% lags: number of lags
% breaks: rows in ydata and xdata after which there is a break. This allows for
% discontinuities in the data (e.g. war years) and for the possibility of
% adding dummy observations to implement a prior. This must be a column vector.
% Note that a single dummy observation becomes lags+1 rows of the data matrix,
% with a break separating it from the rest of the data. The function treats the
% first lags observations at the top and after each "break" in ydata and xdata as
% initial conditions.
% lambda: weight on "co-persistence" prior dummy observations. This expresses
% belief that when data on *all* y's are stable at their initial levels, they will
% tend to persist at that level. lambda=5 is a reasonable first try. With lambda<0,
% constant term is not included in the dummy observation, so that stationary models
% with means equal to initial ybar do not fit the prior mean. With lambda>0, the prior
% implies that large constants are unlikely if unit roots are present.
% mu: weight on "own persistence" prior dummy observation. Expresses belief
% that when y_i has been stable at its initial level, it will tend to persist
% at that level, regardless of the values of other variables. There is
% one of these for each variable. A reasonable first guess is mu=2.
% The program assumes that the first lags rows of ydata and xdata are real data, not dummies.
% Dummy observations should go at the end, if any. If pre-sample x's are not available,
% repeating the initial xdata(lags+1,:) row or copying xdata(lags+1:2*lags,:) into
% xdata(1:lags,:) are reasonable subsititutes. These values are used in forming the
% persistence priors.
% Original file downloaded from:
% http://sims.princeton.edu/yftp/VARtools/matlab/rfvar3.m
[T,nvar] = size(ydata);
nox = isempty(xdata);
if ~nox
[T2,nx] = size(xdata);
else
T2 = T;
nx = 0;
xdata = zeros(T2,0);
end
% note that x must be same length as y, even though first part of x will not be used.
% This is so that the lags parameter can be changed without reshaping the xdata matrix.
if T2 ~= T, error('Mismatch of x and y data lengths'),end
if nargin < 4
nbreaks = 0;
breaks = [];
else
nbreaks = length(breaks);
end
breaks = [0;breaks;T];
smpl = [];
for nb = 1:nbreaks+1
smpl = [smpl;[breaks(nb)+lags+1:breaks(nb+1)]'];
end
Tsmpl = size(smpl,1);
X = zeros(Tsmpl,nvar,lags);
for is = 1:length(smpl)
X(is,:,:) = ydata(smpl(is)-(1:lags),:)';
end
X = [X(:,:) xdata(smpl,:)];
y = ydata(smpl,:);
% Everything now set up with input data for y=Xb+e
% Add persistence dummies
if lambda ~= 0 || mu > 0
ybar = mean(ydata(1:lags,:),1);
if ~nox
xbar = mean(xdata(1:lags,:),1);
else
xbar = [];
end
if lambda ~= 0
if lambda>0
xdum = lambda*[repmat(ybar,1,lags) xbar];
else
lambda = -lambda;
xdum = lambda*[repmat(ybar,1,lags) zeros(size(xbar))];
end
ydum = zeros(1,nvar);
ydum(1,:) = lambda*ybar;
y = [y;ydum];
X = [X;xdum];
end
if mu>0
xdum = [repmat(diag(ybar),1,lags) zeros(nvar,nx)]*mu;
ydum = mu*diag(ybar);
X = [X;xdum];
y = [y;ydum];
end
end
% Compute OLS regression and residuals
[vl,d,vr] = svd(X,0);
di = 1./diag(d);
B = (vr.*repmat(di',nvar*lags+nx,1))*vl'*y;
u = y-X*B;
xxi = vr.*repmat(di',nvar*lags+nx,1);
xxi = xxi*xxi';
var.B = B;
var.u = u;
var.xxi = xxi;

View File

@ -5,7 +5,7 @@ function cprod = cartesian_product_of_sets(varargin)
%! @deftypefn {Function File} {@var{cprod} =} cartesian_product_of_sets (@var{a},@var{b}, ...)
%! @anchor{cartesian_product_of_sets}
%! @sp 1
%! Computes A_1 * A_2 * .... * A_n with a generic set A_i = {e_1,e_2,e_3,...} where e_i is a string
%! Computes A_1 * A_2 * .... * A_n with a generic set A_i = {e_1,e_2,e_3,...} where e_i is a string
%! or a number. It is assumed that each element e_i is unique in set A_i.
%! @sp 2
%! @strong{Inputs}
@ -31,7 +31,7 @@ function cprod = cartesian_product_of_sets(varargin)
%! @end deftypefn
%@eod:
% Copyright (C) 2011-2012 Dynare Team
% Copyright (C) 2011-2017 Dynare Team
%
% This file is part of Dynare.
%

View File

@ -2,14 +2,14 @@ function cellofchar2mfile(fname, c, cname)
% Write a cell of char in a matlab script.
%
% INPUTS
% INPUTS
% - fname [string] name of the file where c is to be saved.
% - c [cell] a two dimensional cell of char.
%
% OUTPUTS
% OUTPUTS
% None.
% Copyright (C) 2015 Dynare Team
% Copyright (C) 2015-2017 Dynare Team
%
% This file is part of Dynare.
%

View File

@ -40,7 +40,7 @@ function [eigenvalues_,result,info] = check(M, options, oo)
%! @end deftypefn
%@eod:
% Copyright (C) 2001-2013 Dynare Team
% Copyright (C) 2001-2014 Dynare Team
%
% This file is part of Dynare.
%

View File

@ -2,7 +2,7 @@ function check_dsge_var_model(Model, EstimatedParameters, BayesInfo)
% Check if the dsge model can be estimated with the DSGE-VAR approach.
% Copyright (C) 2013 Dynare Team
% Copyright (C) 2013-2014 Dynare Team
%
% This file is part of Dynare.
%

View File

@ -12,7 +12,7 @@ function estim_params=check_for_calibrated_covariances(xparam1,estim_params,M)
% Notes: M is local to this function and not updated when calling
% set_all_parameters
% Copyright (C) 2013 Dynare Team
% Copyright (C) 2013-2017 Dynare Team
%
% This file is part of Dynare.
%
@ -42,7 +42,7 @@ calibrated_covariance_pos=covariance_pos(~ismember(covariance_pos,correlation_po
if any(calibrated_covariance_pos)
[rows, columns]=ind2sub(size(M.Sigma_e),calibrated_covariance_pos); %find linear indices of lower triangular covariance entries
estim_params.calibrated_covariances.position=[calibrated_covariance_pos;sub2ind(size(M.Sigma_e),columns,rows)]; %get linear entries of upper triangular parts
estim_params.calibrated_covariances.cov_value=Sigma_e_calibrated(estim_params.calibrated_covariances.position);
estim_params.calibrated_covariances.cov_value=Sigma_e_calibrated(estim_params.calibrated_covariances.position);
end
correlation_pos_ME=find(tril(M.Correlation_matrix_ME,-1)); %off-diagonal elements set by correlations after accounting for estimation
@ -50,6 +50,5 @@ calibrated_covariance_pos_ME=covariance_pos_ME(~ismember(covariance_pos_ME,corre
if any(calibrated_covariance_pos_ME)
[rows, columns]=ind2sub(size(M.H),calibrated_covariance_pos_ME); %find linear indices of lower triangular covariance entries
estim_params.calibrated_covariances_ME.position=[calibrated_covariance_pos_ME;sub2ind(size(M.H),columns,rows)]; %get linear entries of upper triangular parts
estim_params.calibrated_covariances_ME.cov_value=H_calibrated(estim_params.calibrated_covariances_ME.position);
estim_params.calibrated_covariances_ME.cov_value=H_calibrated(estim_params.calibrated_covariances_ME.position);
end

View File

@ -1,20 +1,20 @@
function varlist = check_list_of_variables(options_, M_, varlist)
% This function defines, if necessary, the list of endogenous variables
% for which the posterior statistics have to be computed.
%
% for which the posterior statistics have to be computed.
%
% INPUTS
%
% INPUTS
%
% options_ [structure] Dynare structure.
% M_ [structure] Dynare structure (related to model definition).
% varlist [string] Array of strings with name of the endogenous variables.
%
% OUTPUTS
% varlist [string]
%
%
% OUTPUTS
% varlist [string]
%
% SPECIAL REQUIREMENTS
% Copyright (C) 2003-2014 Dynare Team
% Copyright (C) 2003-2017 Dynare Team
%
% This file is part of Dynare.
%
@ -33,8 +33,8 @@ function varlist = check_list_of_variables(options_, M_, varlist)
%get uniques
[junk1,junk2,index_uniqes] = varlist_indices(varlist,M_.endo_names);
varlist=varlist(index_uniqes,:);
[junk1,junk2,index_uniques] = varlist_indices(varlist,M_.endo_names);
varlist=varlist(index_uniques,:);
msg = 0;
if options_.dsge_var && options_.bayesian_irf
@ -63,7 +63,7 @@ if ~isempty(varlist) && ~isempty(options_.endo_vars_for_moment_computations_in_e
error('You cannot use the consider_all_endogenous or consider_all_observed options when listing variables after the estimation command')
elseif isempty(varlist) && ~isempty(options_.endo_vars_for_moment_computations_in_estimation)
if strcmp(options_.endo_vars_for_moment_computations_in_estimation,'all_endogenous_variables')
varlist = M_.endo_names(1:M_.orig_endo_nbr, :);
varlist = M_.endo_names(1:M_.orig_endo_nbr, :);
elseif strcmp(options_.endo_vars_for_moment_computations_in_estimation,'only_observed_variables')
varlist = char(options_.varobs');
else
@ -106,7 +106,7 @@ elseif isempty(varlist) && isempty(options_.endo_vars_for_moment_computations_in
end
if ~isempty(cas)
string = [ cas , ' will be computed for the ' num2str(M_.endo_nbr) ' endogenous variables'];
string = [ string ' of your model, this can be very long....'];
string = [ string ' of your model, this can be very long....'];
format_text(string, 10)
if options_.nointeractive
% Default behaviour is to consider all the endogenous variables.

View File

@ -1,6 +1,6 @@
function check_matlab_path(change_path_flag)
% Copyright (C) 2015-2016 Dynare Team
% Copyright (C) 2015-2017 Dynare Team
%
% This file is part of Dynare.
%
@ -82,7 +82,7 @@ else
warning(msg);
skipline()
rmpath(DYNARE_PATH)
addpath(DYNARE_PATH)
addpath(DYNARE_PATH)
end
warning on backtrace
end
@ -90,7 +90,7 @@ else
else
% Check that the user did not put all the subfolders in the path.
% => If DYNARE_PATH/qz is in the path while mjdgges dll is available
% it most likely means that user wrongly put all subfolders in the
% it most likely means that user wrongly put all subfolders in the
% matlab's path!
mexpath = add_path_to_mex_files([DYNARE_PATH filesep], false);
MATLAB_PATH = path2cell(MATLAB_PATH);
@ -107,14 +107,14 @@ else
end
function q = path2cell(p)
% Converts the output of path() to a cell
% Converts the output of path() to a cell
s = strfind(p,pathsep);
n = length(s)+1;
q = cell(n,1);
q(1) = {p(1:s(1)-1)};
q(n) = {p(s(end)+1:end)};
for i=2:n-1
q(i) = {p(s(i-1)+1:s(i)-1)};
q(i) = {p(s(i-1)+1:s(i)-1)};
end
function flist = getallroutinenames(p, excludedsubfolders)

View File

@ -1,6 +1,6 @@
function check_model(DynareModel)
% Copyright (C) 2005-2012 Dynare Team
% Copyright (C) 2005-2013 Dynare Team
%
% This file is part of Dynare.
%

View File

@ -1,13 +1,13 @@
function [info,description] = check_posterior_analysis_data(type,M_)
% function [info,description] = check_posterior_analysis_data(type,M_)
% Checks the status of posterior analysis and in particular if files need to be
% Checks the status of posterior analysis and in particular if files need to be
% created or updated; called by posterior_analysis.m
%
%
% Inputs:
% type [string] name of the posterior moment considered
% M_ [structure] Dynare model structure
%
% Outputs:
%
% Outputs:
% info [scalar] return code
% info = 1; % select_posterior_draws has to be called first.
% info = 2; % _posterior_draws files have to be updated.
@ -17,7 +17,7 @@ function [info,description] = check_posterior_analysis_data(type,M_)
% info = 6; % Ok (nothing to do ;-)
% description [string] Message corresponding to info
% Copyright (C) 2008-2015 Dynare Team
% Copyright (C) 2008-2017 Dynare Team
%
% This file is part of Dynare.
%
@ -99,7 +99,7 @@ else
number_of_the_last_post_data_file = length(pdfinfo);
name_of_the_last_post_data_file = ...
[ pwd filesep MetropolisFolder filesep ...
M_.fname '_' ...
M_.fname '_' ...
generic_post_data_file_name ...
int2str(number_of_the_last_post_data_file) ...
'.mat' ];
@ -108,7 +108,7 @@ else
info = 5; % posterior data files have to be updated.
if nargout>1
description = 'posterior data files have to be updated.';
end
end
else
info = 6; % Ok (nothing to do ;-)
if nargout>1

View File

@ -13,7 +13,7 @@ function [posterior_sampler_options, options_] = check_posterior_sampler_options
% SPECIAL REQUIREMENTS
% none
% Copyright (C) 2015 Dynare Team
% Copyright (C) 2015-2017 Dynare Team
%
% This file is part of Dynare.
%
@ -32,32 +32,32 @@ function [posterior_sampler_options, options_] = check_posterior_sampler_options
init=0;
if isempty(posterior_sampler_options),
if isempty(posterior_sampler_options)
init=1;
end
if init,
if init
% set default options and user defined options
posterior_sampler_options.posterior_sampling_method = options_.posterior_sampler_options.posterior_sampling_method;
posterior_sampler_options.bounds = bounds;
switch posterior_sampler_options.posterior_sampling_method
case 'random_walk_metropolis_hastings'
posterior_sampler_options.parallel_bar_refresh_rate=50;
posterior_sampler_options.serial_bar_refresh_rate=3;
posterior_sampler_options.parallel_bar_title='RWMH';
posterior_sampler_options.serial_bar_title='RW Metropolis-Hastings';
% default options
posterior_sampler_options = add_fields_(posterior_sampler_options,options_.posterior_sampler_options.rwmh);
% user defined options
if ~isempty(options_.posterior_sampler_options.sampling_opt)
options_list = read_key_value_string(options_.posterior_sampler_options.sampling_opt);
for i=1:rows(options_list)
switch options_list{i,1}
case 'proposal_distribution'
if ~(strcmpi(options_list{i,2}, 'rand_multivariate_student') || ...
strcmpi(options_list{i,2}, 'rand_multivariate_normal'))
@ -66,15 +66,15 @@ if init,
else
posterior_sampler_options.proposal_distribution=options_list{i,2};
end
case 'student_degrees_of_freedom'
if options_list{i,2} <= 0
error('initial_estimation_checks:: the student_degrees_of_freedom takes a positive integer argument');
else
posterior_sampler_options.student_degrees_of_freedom=options_list{i,2};
end
case 'use_mh_covariance_matrix'
% indicates to use the covariance matrix from previous iterations to
% define the covariance of the proposal distribution
@ -101,23 +101,23 @@ if init,
end
end
end
case 'tailored_random_block_metropolis_hastings'
posterior_sampler_options.parallel_bar_refresh_rate=5;
posterior_sampler_options.serial_bar_refresh_rate=1;
posterior_sampler_options.parallel_bar_title='TaRB-MH';
posterior_sampler_options.serial_bar_title='TaRB Metropolis-Hastings';
% default options
posterior_sampler_options = add_fields_(posterior_sampler_options,options_.posterior_sampler_options.tarb);
% user defined options
if ~isempty(options_.posterior_sampler_options.sampling_opt)
options_list = read_key_value_string(options_.posterior_sampler_options.sampling_opt);
for i=1:rows(options_list)
switch options_list{i,1}
case 'proposal_distribution'
if ~(strcmpi(options_list{i,2}, 'rand_multivariate_student') || ...
strcmpi(options_list{i,2}, 'rand_multivariate_normal'))
@ -126,18 +126,18 @@ if init,
else
posterior_sampler_options.proposal_distribution=options_list{i,2};
end
case 'student_degrees_of_freedom'
if options_list{i,2} <= 0
error('initial_estimation_checks:: the student_degrees_of_freedom takes a positive integer argument');
else
posterior_sampler_options.student_degrees_of_freedom=options_list{i,2};
end
case 'mode_compute'
posterior_sampler_options.mode_compute=options_list{i,2};
case 'optim'
posterior_sampler_options.optim_opt=options_list{i,2};
@ -162,31 +162,31 @@ if init,
end
case 'save_tmp_file'
posterior_sampler_options.save_tmp_file = options_list{i,2};
otherwise
warning(['tarb_sampler: Unknown option (' options_list{i,1} ')!'])
end
end
end
case 'independent_metropolis_hastings'
posterior_sampler_options.parallel_bar_refresh_rate=50;
posterior_sampler_options.serial_bar_refresh_rate=3;
posterior_sampler_options.parallel_bar_title='IMH';
posterior_sampler_options.serial_bar_title='Ind. Metropolis-Hastings';
% default options
posterior_sampler_options = add_fields_(posterior_sampler_options,options_.posterior_sampler_options.imh);
% user defined options
if ~isempty(options_.posterior_sampler_options.sampling_opt)
options_list = read_key_value_string(options_.posterior_sampler_options.sampling_opt);
for i=1:rows(options_list)
switch options_list{i,1}
case 'proposal_distribution'
if ~(strcmpi(options_list{i,2}, 'rand_multivariate_student') || ...
strcmpi(options_list{i,2}, 'rand_multivariate_normal'))
@ -195,22 +195,22 @@ if init,
else
posterior_sampler_options.proposal_distribution=options_list{i,2};
end
case 'student_degrees_of_freedom'
if options_list{i,2} <= 0
error('initial_estimation_checks:: the student_degrees_of_freedom takes a positive integer argument');
else
posterior_sampler_options.student_degrees_of_freedom=options_list{i,2};
end
case 'use_mh_covariance_matrix'
% indicates to use the covariance matrix from previous iterations to
% define the covariance of the proposal distribution
% default = 0
posterior_sampler_options.use_mh_covariance_matrix = options_list{i,2};
options_.use_mh_covariance_matrix = options_list{i,2};
case 'save_tmp_file'
posterior_sampler_options.save_tmp_file = options_list{i,2};
@ -219,17 +219,17 @@ if init,
end
end
end
case 'slice'
posterior_sampler_options.parallel_bar_refresh_rate=1;
posterior_sampler_options.serial_bar_refresh_rate=1;
posterior_sampler_options.parallel_bar_title='SLICE';
posterior_sampler_options.serial_bar_title='SLICE';
% default options
posterior_sampler_options = add_fields_(posterior_sampler_options,options_.posterior_sampler_options.slice);
% user defined options
if ~isempty(options_.posterior_sampler_options.sampling_opt)
options_list = read_key_value_string(options_.posterior_sampler_options.sampling_opt);
@ -242,7 +242,7 @@ if init,
% <use_mh_covariance_matrix> or <slice_initialize_with_mode>
% default = 0
posterior_sampler_options.rotated = options_list{i,2};
case 'mode'
% for multimodal posteriors, provide the list of modes as a
% matrix, ordered by column, i.e. [x1 x2 x3] for three
@ -253,19 +253,19 @@ if init,
% This will automatically trigger <rotated>
% default = []
tmp_mode = options_list{i,2};
for j=1:size(tmp_mode,2),
for j=1:size(tmp_mode,2)
posterior_sampler_options.mode(j).m = tmp_mode(:,j);
end
case 'mode_files'
% for multimodal posteriors provide the name of
% a file containing a variable array xparams = [nparam * nmodes]
% one column per mode. With this info, the code will automatically
% set the <mode> option.
% set the <mode> option.
% This will automatically trigger <rotated>
% default = []
posterior_sampler_options.mode_files = options_list{i,2};
case 'slice_initialize_with_mode'
% the default for slice is to set mode_compute = 0 in the
% preprocessor and start the chain(s) from a random location in the prior.
@ -275,7 +275,7 @@ if init,
% iterations.
% default = 0
posterior_sampler_options.slice_initialize_with_mode = options_list{i,2};
case 'initial_step_size'
% sets the initial size of the interval in the STEPPING-OUT PROCEDURE
% the initial_step_size must be a real number in [0, 1],
@ -298,18 +298,18 @@ if init,
case 'save_tmp_file'
posterior_sampler_options.save_tmp_file = options_list{i,2};
otherwise
warning(['slice_sampler: Unknown option (' options_list{i,1} ')!'])
end
end
end
% slice posterior sampler does not require mode or hessian to run
% needs to be set to 1 to skip parts in dynare_estimation_1.m
% requiring posterior maximization/calibrated smoother before MCMC
options_.mh_posterior_mode_estimation=1;
if ~ posterior_sampler_options.slice_initialize_with_mode
% by default, slice sampler should trigger
% mode_compute=0 and
@ -327,52 +327,52 @@ if init,
options_.mh_posterior_mode_estimation=0;
end
end
if any(isinf(bounds.lb)) || any(isinf(bounds.ub)),
if any(isinf(bounds.lb)) || any(isinf(bounds.ub))
skipline()
disp('some priors are unbounded and prior_trunc is set to zero')
error('The option "slice" is inconsistent with prior_trunc=0.')
end
% moreover slice must be associated to:
% options_.mh_posterior_mode_estimation = 0;
% this is done below, but perhaps preprocessing should do this?
if ~isempty(posterior_sampler_options.mode)
% multimodal case
posterior_sampler_options.rotated = 1;
posterior_sampler_options.WR=[];
end
% posterior_sampler_options = set_default_option(posterior_sampler_options,'mode_files',[]);
posterior_sampler_options.W1=posterior_sampler_options.initial_step_size*(bounds.ub-bounds.lb);
if options_.load_mh_file,
if options_.load_mh_file
posterior_sampler_options.slice_initialize_with_mode = 0;
else
if ~posterior_sampler_options.slice_initialize_with_mode,
if ~posterior_sampler_options.slice_initialize_with_mode
posterior_sampler_options.invhess=[];
end
end
if ~isempty(posterior_sampler_options.mode_files), % multimodal case
if ~isempty(posterior_sampler_options.mode_files) % multimodal case
modes = posterior_sampler_options.mode_files; % these can be also mean files from previous parallel slice chains
load(modes, 'xparams')
if size(xparams,2)<2,
if size(xparams,2)<2
error(['check_posterior_sampler_options:: Variable xparams loaded in file <' modes '> has size [' int2str(size(xparams,1)) 'x' int2str(size(xparams,2)) ']: it must contain at least two columns, to allow multi-modal sampling.'])
end
for j=1:size(xparams,2),
for j=1:size(xparams,2)
mode(j).m=xparams(:,j);
end
posterior_sampler_options.mode = mode;
posterior_sampler_options.rotated = 1;
posterior_sampler_options.WR=[];
end
otherwise
error('check_posterior_sampler_options:: Unknown posterior_sampling_method option %s ',posterior_sampler_options.posterior_sampling_method);
end
return
end
@ -386,7 +386,7 @@ if ~strcmp(posterior_sampler_options.posterior_sampling_method,'slice')
end
end
if options_.load_mh_file && posterior_sampler_options.use_mh_covariance_matrix,
if options_.load_mh_file && posterior_sampler_options.use_mh_covariance_matrix
[junk, invhess] = compute_mh_covariance_matrix;
posterior_sampler_options.invhess = invhess;
end
@ -396,10 +396,8 @@ end
% check specific options for slice sampler
if strcmp(posterior_sampler_options.posterior_sampling_method,'slice')
invhess = posterior_sampler_options.invhess;
if posterior_sampler_options.rotated,
if isempty(posterior_sampler_options.mode_files) && isempty(posterior_sampler_options.mode), % rotated unimodal
if posterior_sampler_options.rotated
if isempty(posterior_sampler_options.mode_files) && isempty(posterior_sampler_options.mode) % rotated unimodal
if ~options_.cova_compute && ~(options_.load_mh_file && posterior_sampler_options.use_mh_covariance_matrix)
skipline()
disp('check_posterior_sampler_options:: I cannot start rotated slice sampler because')
@ -419,19 +417,13 @@ if strcmp(posterior_sampler_options.posterior_sampling_method,'slice')
posterior_sampler_options.WR=sqrt(diag(D))*3;
end
else
if ~options_.load_mh_file && ~posterior_sampler_options.slice_initialize_with_mode,
if ~options_.load_mh_file && ~posterior_sampler_options.slice_initialize_with_mode
posterior_sampler_options.invhess=[];
end
end
% needs to be re-set to zero otherwise posterior analysis is filtered
% out in dynare_estimation_1.m
options_.mh_posterior_mode_estimation = 0;
else
end
return
@ -442,6 +434,3 @@ fnam = fieldnames(sampler_options);
for j=1:length(fnam)
posterior_sampler_options.(fnam{j}) = sampler_options.(fnam{j});
end

View File

@ -1,13 +1,13 @@
function [info,description] = check_prior_analysis_data(type,M_)
% function [info,description] = check_prior_analysis_data(type,M_)
% Checks the status of prior analysis and in particular if files need to be
% Checks the status of prior analysis and in particular if files need to be
% created or updated; called by prior_analysis.m
%
%
% Inputs:
% type [string] name of the posterior moment considered
% M_ [structure] Dynare model structure
%
% Outputs:
%
% Outputs:
% info [scalar] return code
% info = 1; % prior_sampler has to be called first.
% info = 2; % _prior_draws files have to be updated.
@ -18,7 +18,7 @@ function [info,description] = check_prior_analysis_data(type,M_)
% description [string] Message corresponding to info
% Copyright (C) 2009-2011 Dynare Team
% Copyright (C) 2009-2017 Dynare Team
%
% This file is part of Dynare.
%
@ -40,7 +40,7 @@ if nargout>1
description = '';
end
%% Get informations about prior draws files.
%% Get informations about prior draws files.
if ~exist([ M_.dname '/prior/draws'],'dir')
disp('check_prior_analysis_data:: Can''t find any prior draws file!')
return
@ -65,7 +65,7 @@ else
end
return
else
info = 3; % Nothing to do!
info = 3; % Nothing to do!
if nargout>1
description = 'prior draws files are up to date.';
end
@ -98,9 +98,9 @@ else
name_of_the_last_prior_data_file = pdinfo(end).name;
pdfdate = pdinfo(end).datenum;
% /!\ REMARK /!\
% The user can change the model or the value of a calibrated
% parameter without changing the prior. In this case the (prior)
% moments should be computed. But this case cannot be detected!!!
% The user can change the model or the value of a calibrated
% parameter without changing the prior. In this case the (prior)
% moments should be computed. But this case cannot be detected!!!
if pdfdate<date_of_the_last_prior_draw_file
info = 5; % prior data files have to be updated.
if nargout>1

View File

@ -5,12 +5,12 @@ function check_prior_bounds(xparam1,bounds,M_,estim_params_,options_,bayestopt_)
% -xparam1 [double] vector of parameters to be estimated (initial values)
% -bounds [vector] vector containing the lower and upper
% bounds
% -M_ [structure] characterizing the model.
% -M_ [structure] characterizing the model.
% -estim_params_ [structure] characterizing parameters to be estimated
% -options_ [structure] characterizing the options
% -bayestopt_ [structure] characterizing priors
% Copyright (C) 2013 Dynare Team
% Copyright (C) 2013-2017 Dynare Team
%
% This file is part of Dynare.
%

View File

@ -1,17 +1,17 @@
function [R,indef, E, P]=chol_SE(A,pivoting)
% [R,indef, E, P]=chol_SE(A,pivoting)
% Performs a (modified) Cholesky factorization of the form
%
% Performs a (modified) Cholesky factorization of the form
%
% P'*A*P + E = R'*R
%
%
% As detailed in Schnabel/Eskow (1990), the factorization has 2 phases:
% Phase 1: While A can still be positive definite, pivot on the maximum diagonal element.
% Check that the standard Cholesky update would result in a positive diagonal
% at the current iteration. If so, continue with the normal Cholesky update.
% Otherwise switch to phase 2.
% If A is safely positive definite, stage 1 is never left, resulting in
% the standard Cholesky decomposition.
%
% at the current iteration. If so, continue with the normal Cholesky update.
% Otherwise switch to phase 2.
% If A is safely positive definite, stage 1 is never left, resulting in
% the standard Cholesky decomposition.
%
% Phase 2: Pivot on the minimum of the negatives of the lower Gershgorin bound
% estimates. To prevent negative diagonals, compute the amount to add to the
% pivot element and add this. Then, do the Cholesky update and update the estimates of the
@ -21,39 +21,40 @@ function [R,indef, E, P]=chol_SE(A,pivoting)
% - During factorization, L=R' is stored in the lower triangle of the original matrix A,
% miminizing the memory requirements
% - Conforming with the original Schnabel/Eskow (1990) algorithm:
% - at each iteration the updated Gershgorin bounds are estimated instead of recomputed,
% - at each iteration the updated Gershgorin bounds are estimated instead of recomputed,
% reducing the computational requirements from o(n^3) to o (n^2)
% - For the last 2 by 2 submatrix, an eigenvalue-based decomposition is used
% - While pivoting is not necessary, it improves the size of E, the add-on on to the diagonal. But this comes at
% - While pivoting is not necessary, it improves the size of E, the add-on on to the diagonal. But this comes at
% the cost of introduding a permutation.
%
%
% Inputs
% A [n*n] Matrix to be factorized
% pivoting [scalar] dummy whether pivoting is used
%
% Outputs
% R [n*n] originally stored in lower triangular portion of A, including the main diagonal
%
% E [n*1] Elements added to the diagonal of A
% P [n*1] record of how the rows and columns of the matrix were permuted while
% performing the decomposition
% INPUTS
% - A [n*n] Matrix to be factorized
% - pivoting [scalar] dummy whether pivoting is used
%
% OUTPUTS
% - R [n*n] originally stored in lower triangular portion of A, including the main diagonal
%
% - E [n*1] Elements added to the diagonal of A
% - P [n*1] record of how the rows and columns of the matrix were permuted while
% performing the decomposition
%
% REFERENCES:
% This implementation is based on
% This implementation is based on
%
% Robert B. Schnabel and Elizabeth Eskow. 1990. "A New Modified Cholesky
% Factorization," SIAM Journal of Scientific Statistical Computating,
% 11, 6: 1136-58.
%
%
% Elizabeth Eskow and Robert B. Schnabel 1991. "Algorithm 695 - Software for a New Modified Cholesky
% Factorization," ACM Transactions on Mathematical Software, Vol 17, No 3: 306-312
%
%
%
%
% Author: Johannes Pfeifer based on Eskow/Schnabel (1991)
%
% Copyright (C) 2015 Johannes Pfeifer
% Copyright (C) 2015 Dynare Team
% Copyright (C) 2015-2017 Dynare Team
%
% This file is part of Dynare.
%
% Dynare is free software: you can redistribute it and/or modify
@ -148,7 +149,7 @@ for j = 1:n-1
% Swap elements of the permutation vector
itemp = P(j);
P(j) = P(imaxd);
P(imaxd) = itemp;
P(imaxd) = itemp;
end
end
% check to see whether the normal Cholesky update for this
@ -183,13 +184,13 @@ for j = 1:n-1
g=gersh_nested(A,j,n);
end
end
% PHASE 2
if ~phase1
if j ~= n-1
if pivoting
% Find the minimum negative Gershgorin bound
[tmp,iming] = min(g(j:n));
[tmp,iming] = min(g(j:n));
iming=iming+j-1;
% Pivot to the top the row and column with the
% minimum negative Gershgorin bound
@ -229,9 +230,9 @@ for j = 1:n-1
end
% Calculate delta and add to the diagonal. delta=max{0,-A(j,j) + max{normj,taugam},delta_previous}
% where normj=sum of |A(i,j)|,for i=1,n, delta_previous is the delta computed at the previous iter and taugam is tau1*gamma.
normj=sum(abs(A(j+1:n,j)));
delta = max([0;delta;-A(j,j)+normj;-A(j,j)+taugam]); % get adjustment based on formula on bottom of p. 309 of Eskow/Schnabel (1991)
E(j) = delta;
@ -259,9 +260,9 @@ for j = 1:n-1
% Find eigenvalues of final 2 by 2 submatrix
% Find delta such that:
% 1. the l2 condition number of the final 2X2 submatrix + delta*I <= tau2
% 2. delta >= previous delta,
% 2. delta >= previous delta,
% 3. min(eigvals) + delta >= tau2 * gamma, where min(eigvals) is the smallest eigenvalue of the final 2X2 submatrix
A(n-1,n)=A(n,n-1); %set value above diagonal for computation of eigenvalues
eigvals = eig(A(n-1:n,n-1:n));
delta = max([ 0 ; delta ; -min(eigvals)+tau2*max((max(eigvals)-min(eigvals))/(1-tau1),gammma) ]); %Formula 5.3.2 of Schnabel/Eskow (1990)
@ -272,7 +273,7 @@ for j = 1:n-1
E(n-1) = delta;
E(n) = delta;
end
% Final update
A(n-1,n-1) = sqrt(A(n-1,n-1));
A(n,n-1) = A(n,n-1)/A(n-1,n-1);
@ -293,17 +294,16 @@ function g=gersh_nested(A,j,n)
g=zeros(n,1);
for ii = j:n
if ii == 1;
if ii == 1
sum_up_to_i = 0;
else
sum_up_to_i = sum(abs(A(ii,j:(ii-1))));
end;
if ii == n;
end
if ii == n
sum_after_i = 0;
else
sum_after_i = sum(abs(A((ii+1):n,ii)));
end;
end
g(ii) = sum_up_to_i + sum_after_i- A(ii,ii);
end
end

View File

@ -2,7 +2,7 @@ function clear_persistent_variables(folder, writelistofroutinestobecleared)
% Clear all the functions with persistent variables in directory folder (and subdirectories).
% Copyright (C) 2015 Dynare Team
% Copyright (C) 2015-2017 Dynare Team
%
% This file is part of Dynare.
%
@ -58,5 +58,5 @@ if writelistofroutinestobecleared
return
end
list_of_functions_to_be_cleared;
list_of_functions_to_be_cleared;
clear(list_of_functions{:});

View File

@ -11,7 +11,7 @@ function varargout = prior(varargin)
% SPECIAL REQUIREMENTS
% none
% Copyright (C) 2015-2016 Dynare Team
% Copyright (C) 2015-2017 Dynare Team
%
% This file is part of Dynare.
%
@ -42,7 +42,7 @@ if isempty(varargin) || ( isequal(length(varargin), 1) && isequal(varargin{1},'h
return
end
global options_ M_ estim_params_ bayestopt_ oo_
global options_ M_ estim_params_ bayestopt_ oo_
donesomething = false;

View File

@ -1,14 +1,14 @@
function collect_latex_files
% function collect_LaTeX_Files;
% Creates TeX-File embedding all eps-loaders created for current mod-file
%
%
% Inputs: none
%
% Notes:
%
% Notes:
% - The packages loaded enable pdflatex to run
% - The _dynamic and _static TeX-model files are not included as they are standalone TeX-files
% Copyright (C) 2015-16 Dynare Team
% Copyright (C) 2015-2017 Dynare Team
%
% This file is part of Dynare.
%
@ -49,7 +49,7 @@ for ii=1:length(TeX_Files)
~strcmp(TeX_Files(ii).name,[M_.fname,'_static.tex']) && ...
~strcmp(TeX_Files(ii).name,[M_.fname,'_original.tex']) && ...
~strcmp(TeX_Files(ii).name,[M_.fname,'_TeX_binder.tex'])
fprintf(fid,'%s \n',['\include{',f_name,'}']);
fprintf(fid,'%s \n',['\include{',f_name,'}']);
end
end
@ -58,7 +58,7 @@ TeX_Files=dir([M_.dname filesep 'Output' filesep M_.fname '*.tex']);
for ii=1:length(TeX_Files)
[pathstr,f_name,ext] = fileparts(TeX_Files(ii).name);
if ~strcmp(TeX_Files(ii).name,f_name_binder)
fprintf(fid,'%s \n',['\include{', M_.dname '/Output' '/',f_name,'}']);
fprintf(fid,'%s \n',['\include{', M_.dname '/Output' '/',f_name,'}']);
end
end
@ -67,7 +67,7 @@ TeX_Files=dir([M_.dname filesep 'graphs' filesep M_.fname '*.tex']);
for ii=1:length(TeX_Files)
[pathstr,f_name,ext] = fileparts(TeX_Files(ii).name);
if ~strcmp(TeX_Files(ii).name,f_name_binder)
fprintf(fid,'%s \n',['\include{', M_.dname '/graphs' '/',f_name,'}']);
fprintf(fid,'%s \n',['\include{', M_.dname '/graphs' '/',f_name,'}']);
end
end
@ -76,7 +76,7 @@ TeX_Files=dir([M_.dname filesep 'identification' filesep M_.fname '*.tex']);
for ii=1:length(TeX_Files)
[pathstr,f_name,ext] = fileparts(TeX_Files(ii).name);
if ~strcmp(TeX_Files(ii).name,f_name_binder)
fprintf(fid,'%s \n',['\include{', M_.dname '/identification' '/',f_name,'}']);
fprintf(fid,'%s \n',['\include{', M_.dname '/identification' '/',f_name,'}']);
end
end
@ -86,7 +86,7 @@ TeX_Files=dir([M_.dname filesep 'identification' filesep 'Output' filesep M_.fna
for ii=1:length(TeX_Files)
[pathstr,f_name,ext] = fileparts(TeX_Files(ii).name);
if ~strcmp(TeX_Files(ii).name,f_name_binder)
fprintf(fid,'%s \n',['\include{', M_.dname '/identification/Output' '/',f_name,'}']);
fprintf(fid,'%s \n',['\include{', M_.dname '/identification/Output' '/',f_name,'}']);
end
end
@ -95,7 +95,7 @@ TeX_Files=dir([M_.dname filesep 'gsa' filesep M_.fname '*.tex']);
for ii=1:length(TeX_Files)
[pathstr,f_name,ext] = fileparts(TeX_Files(ii).name);
if ~strcmp(TeX_Files(ii).name,f_name_binder)
fprintf(fid,'%s \n',['\include{', M_.dname '/gsa' '/',f_name,'}']);
fprintf(fid,'%s \n',['\include{', M_.dname '/gsa' '/',f_name,'}']);
end
end
@ -104,7 +104,7 @@ TeX_Files=dir([M_.dname filesep 'gsa' filesep 'Output' filesep M_.fname '*.tex'
for ii=1:length(TeX_Files)
[pathstr,f_name,ext] = fileparts(TeX_Files(ii).name);
if ~strcmp(TeX_Files(ii).name,f_name_binder)
fprintf(fid,'%s \n',['\include{', M_.dname '/gsa/Output' '/',f_name,'}']);
fprintf(fid,'%s \n',['\include{', M_.dname '/gsa/Output' '/',f_name,'}']);
end
end

View File

@ -1,36 +1,36 @@
function [Pstar,Pinf] = compute_Pinf_Pstar(mf,T,R,Q,qz_criterium, restrict_columns)
% function [Z,ST,QT,R1,Pstar,Pinf] = schur_statespace_transformation(mf,T,R,Q,qz_criterium, restrict_columns)
% Kitagawa transformation of state space system with a quasi-triangular
% transition matrix with unit roots at the top, but excluding zero columns of the transition matrix.
% transition matrix with unit roots at the top, but excluding zero columns of the transition matrix.
% Computation of Pstar and Pinf for Durbin and Koopman Diffuse filter
%
% The transformed state space is
% The transformed state space is
% y = [ss; z; x];
% s = static variables (zero columns of T)
% z = unit roots
% x = stable roots
% ss = s - z = stationarized static variables
%
% INPUTS
%
% INPUTS
% mf [integer] vector of indices of observed variables in
% state vector
% T [double] matrix of transition
% R [double] matrix of structural shock effects
% Q [double] matrix of covariance of structural shocks
% qz_criterium [double] numerical criterium for unit roots
%
% OUTPUTS
% qz_criterium [double] numerical criterium for unit roots
%
% OUTPUTS
% Pstar [double] matrix of covariance of stationary part
% Pinf [double] matrix of covariance initialization for
% nonstationary part
%
% ALGORITHM
% nonstationary part
%
% ALGORITHM
% Real Schur transformation of transition equation
%
% SPECIAL REQUIREMENTS
% None
% Copyright (C) 2006-2011 Dynare Team
% Copyright (C) 2006-2017 Dynare Team
%
% This file is part of Dynare.
%
@ -49,7 +49,7 @@ function [Pstar,Pinf] = compute_Pinf_Pstar(mf,T,R,Q,qz_criterium, restrict_colum
np = size(T,1);
% perform Kitagawa transformation
% perform Kitagawa transformation
[QT,ST] = schur(T);
e1 = abs(ordeig(ST)) > 2-qz_criterium;
[QT,ST] = ordschur(QT,ST,e1);
@ -112,4 +112,3 @@ end
Pinf = QT*Pinf*QT';
Pstar = QT*Pstar*QT';

View File

@ -1,21 +1,21 @@
function [posterior_mean,posterior_covariance,posterior_mode,posterior_kernel_at_the_mode] = compute_mh_covariance_matrix()
% Estimation of the posterior covariance matrix, posterior mean, posterior mode and evaluation of the posterior kernel at the
% estimated mode, using the draws from a metropolis-hastings. The estimated posterior mode and covariance matrix are saved in
% a file <M_.fname>_mh_mode.mat.
%
% INPUTS
% a file <M_.fname>_mh_mode.mat.
%
% INPUTS
% None.
%
%
% OUTPUTS
% o posterior_mean [double] (n*1) vector, posterior expectation of the parameters.
% o posterior_covariance [double] (n*n) matrix, posterior covariance of the parameters (computed from previous metropolis hastings).
% o posterior_mode [double] (n*1) vector, posterior mode of the parameters.
% o posterior_mode [double] (n*1) vector, posterior mode of the parameters.
% o posterior_kernel_at_the_mode [double] scalar.
%
% SPECIAL REQUIREMENTS
% None.
% Copyright (C) 2006-2011 Dynare Team
% Copyright (C) 2006-2017 Dynare Team
%
% This file is part of Dynare.
%
@ -62,7 +62,7 @@ offset = 0;
for b=1:nblck
first_line = FirstLine;
for n = FirstMhFile:TotalNumberOfMhFiles
load([ BaseName '_mh' int2str(n) '_blck' int2str(b) '.mat'],'x2','logpo2');
load([ BaseName '_mh' int2str(n) '_blck' int2str(b) '.mat'],'x2','logpo2');
[tmp,idx] = max(logpo2);
if tmp>posterior_kernel_at_the_mode
posterior_kernel_at_the_mode = tmp;

View File

@ -4,14 +4,14 @@ function moments=compute_model_moments(dr,M_,options_)
% dr: structure describing model solution
% M_: structure of Dynare options
% options_
%
%
% OUTPUTS
% moments: a cell array containing requested moments
%
%
% SPECIAL REQUIREMENTS
% none
% Copyright (C) 2008-2009 Dynare Team
% Copyright (C) 2008-2017 Dynare Team
%
% This file is part of Dynare.
%

View File

@ -2,21 +2,21 @@ function oo_ = compute_moments_varendo(type,options_,M_,oo_,var_list_)
% Computes the second order moments (autocorrelation function, covariance
% matrix and variance decomposition) distributions for all the endogenous variables selected in
% var_list_. The results are saved in oo_
%
%
% INPUTS:
% type [string] 'posterior' or 'prior'
% options_ [structure] Dynare structure.
% M_ [structure] Dynare structure (related to model definition).
% oo_ [structure] Dynare structure (results).
% var_list_ [string] Array of string with endogenous variable names.
%
%
% OUTPUTS
% oo_ [structure] Dynare structure (results).
%
% SPECIAL REQUIREMENTS
% none
% Copyright (C) 2008-2010 Dynare Team
% Copyright (C) 2008-2017 Dynare Team
%
% This file is part of Dynare.
%
@ -31,7 +31,7 @@ function oo_ = compute_moments_varendo(type,options_,M_,oo_,var_list_)
% 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/>.
% along with Dynare. If not, see <http://www.gnu.org/licenses/>.
fprintf('Estimation::compute_moments_varendo: I''m computing endogenous moments (this may take a while)... ');

View File

@ -1,6 +1,6 @@
function overallacceptanceratio = compute_overall_acceptance_ratio(MetropolisFolder, ModelName)
% Copyright (C) 2013 Dynare Team
% Copyright (C) 2013-2017 Dynare Team
%
% This file is part of Dynare.
%
@ -24,12 +24,12 @@ n = length(mh_history_files);
load([BaseName '_mh_history_' num2str(0)]);
TotalNumberOfDraws = record.MhDraws(end,1);
TotalNumberOfAcceptedProposals = record.AcceptanceRatio*record.MhDraws(end,1);
TotalNumberOfAcceptedProposals = record.AcceptanceRatio*record.MhDraws(end,1);
for i=2:n
load([BaseName '_mh_history_' num2str(i-1)]);
TotalNumberOfDraws = TotalNumberOfDraws + record.MhDraws(end,1);
TotalNumberOfAcceptedProposals = TotalNumberOfAcceptedProposals + record.AcceptanceRatio*record.MhDraws(end,1);
TotalNumberOfAcceptedProposals = TotalNumberOfAcceptedProposals + record.AcceptanceRatio*record.MhDraws(end,1);
end
overallacceptanceratio = TotalNumberOfAcceptedProposals/TotalNumberOfDraws;

View File

@ -18,7 +18,7 @@ function [trend_addition, trend_coeff]=compute_trend_coefficients(M_,DynareOptio
% SPECIAL REQUIREMENTS
% none
% Copyright (C) 2014 Dynare Team
% Copyright (C) 2014-2016 Dynare Team
%
% This file is part of Dynare.
%

View File

@ -1,14 +1,14 @@
function ConditionalVarianceDecomposition = conditional_variance_decomposition(StateSpaceModel, Steps, SubsetOfVariables,sigma_e_is_diagonal)
% This function computes the conditional variance decomposition of a given state space model
% for a subset of endogenous variables.
%
% INPUTS
%
% INPUTS
% StateSpaceModel [structure] Specification of the state space model.
% Steps [integer] 1*h vector of dates.
% SubsetOfVariables [integer] 1*q vector of indices.
%
% OUTPUTS
% ConditionalVarianceDecomposition [double] [n h p] array, where
%
% OUTPUTS
% ConditionalVarianceDecomposition [double] [n h p] array, where
% n is equal to length(SubsetOfVariables)
% h is the number of Steps
% p is the number of state innovations and
@ -16,7 +16,7 @@ function ConditionalVarianceDecomposition = conditional_variance_decomposition(S
%
% [1] In this version, absence of measurement errors is assumed...
% Copyright (C) 2010-2011 Dynare Team
% Copyright (C) 2010-2017 Dynare Team
%
% This file is part of Dynare.
%
@ -77,7 +77,7 @@ for h = 1:length(Steps)
SumOfVariances(:,h) = sum(ConditionalVariance(:,h,:),3);
end
ConditionalVarianceDecomposition = zeros(NumberOfVariables,length(Steps),number_of_state_innovations);
ConditionalVarianceDecomposition = zeros(NumberOfVariables,length(Steps),number_of_state_innovations);
for i=1:number_of_state_innovations
for h = 1:length(Steps)
ConditionalVarianceDecomposition(:,h,i) = squeeze(ConditionalVariance(:,h,i))./SumOfVariances(:,h);

View File

@ -9,7 +9,7 @@ function oo_ = ...
% dname [string] directory name where to save
% fname [string] name of the mod-file
% Steps [integers] horizons at which to conduct decomposition
% exonames [string] (n_exo*char_length) character array with names of exogenous variables
% exonames [string] (n_exo*char_length) character array with names of exogenous variables
% exo [string] name of current exogenous
% variable
% var_list [string] (n_endo*char_length) character array with name
@ -23,7 +23,7 @@ function oo_ = ...
% OUTPUTS
% oo_ [structure] Dynare structure where the results are saved.
% Copyright (C) 2009-2015 Dynare Team
% Copyright (C) 2009-2017 Dynare Team
%
% This file is part of Dynare.
%
@ -106,7 +106,7 @@ for i=1:length(Steps)
p_density(:,:,i) = pp_density;
else
[pp_mean, pp_median, pp_var, hpd_interval, pp_deciles] = ...
posterior_moments(tmp(:,i),0,mh_conf_sig);
posterior_moments(tmp(:,i),0,mh_conf_sig);
end
p_mean(i) = pp_mean;
p_median(i) = pp_median;

View File

@ -1,14 +1,14 @@
function oo_ = McMCDiagnostics(options_, estim_params_, M_, oo_)
% function McMCDiagnostics
% Computes convergence tests
%
% INPUTS
% Computes convergence tests
%
% INPUTS
% options_ [structure]
% estim_params_ [structure]
% M_ [structure]
%
% OUTPUTS
% oo_ [structure]
% OUTPUTS
% oo_ [structure]
%
% SPECIAL REQUIREMENTS
% none
@ -82,12 +82,12 @@ param_name_tex=[];
for jj=1:npar
if options_.TeX
[par_name_temp,par_name_tex_temp]=get_the_name(jj,options_.TeX,M_,estim_params_,options_);
param_name = strvcat(param_name,par_name_temp);
param_name = strvcat(param_name,par_name_temp);
par_name_tex_temp = strrep(par_name_tex_temp,'$','');
param_name_tex = strvcat(param_name_tex,par_name_tex_temp);
else
[par_name_temp]=get_the_name(jj,options_.TeX,M_,estim_params_,options_);
param_name = strvcat(param_name,par_name_temp);
param_name = strvcat(param_name,par_name_temp);
end
Draws = GetAllPosteriorDraws(jj,FirstMhFile,FirstLine,TotalNumberOfMhFiles,NumberOfDraws);
Draws = reshape(Draws,[NumberOfDraws nblck]);
@ -101,7 +101,7 @@ end
my_title='MCMC Inefficiency factors per block';
IFAC_header='Parameter';
IFAC_header_tex='Parameter';
for j=1:nblck,
for j=1:nblck
IFAC_header = char(IFAC_header,['Block ' int2str(j)]);
IFAC_header_tex = char(IFAC_header_tex,['Block~' int2str(j)]);
end
@ -139,9 +139,9 @@ if nblck == 1 % Brooks and Gelman tests need more than one block
first_obs_end_sample = first_obs_begin_sample+round(options_.convergence.geweke.geweke_interval(2)*NumberOfDraws*(1-options_.mh_drop));
param_name=[];
if options_.TeX
param_name_tex=[];
param_name_tex=[];
end
for jj=1:npar
for jj=1:npar
if options_.TeX
[param_name_temp, param_name_tex_temp]= get_the_name(jj,options_.TeX,M_,estim_params_,options_);
param_name_tex = strvcat(param_name_tex,strrep(param_name_tex_temp,'$',''));
@ -152,7 +152,7 @@ if nblck == 1 % Brooks and Gelman tests need more than one block
end
end
fprintf('\nGeweke (1992) Convergence Tests, based on means of draws %d to %d vs %d to %d.\n',first_obs_begin_sample,last_obs_begin_sample,first_obs_end_sample,NumberOfDraws);
fprintf('p-values are for Chi2-test for equality of means.\n');
fprintf('p-values are for Chi2-test for equality of means.\n');
Geweke_header=char('Parameter', 'Post. Mean', 'Post. Std', 'p-val No Taper');
for ii=1:length(options_.convergence.geweke.taper_steps)
Geweke_header=char(Geweke_header,['p-val ' num2str(options_.convergence.geweke.taper_steps(ii)),'% Taper']);
@ -168,12 +168,12 @@ if nblck == 1 % Brooks and Gelman tests need more than one block
end
[results_vec, results_struct] = geweke_moments(param_draws,options_);
convergence_diagnostics_geweke(jj,:)=results_vec;
param_draws1 = param_draws(first_obs_begin_sample:last_obs_begin_sample,:);
param_draws2 = param_draws(first_obs_end_sample:end,:);
[results_vec1] = geweke_moments(param_draws1,options_);
[results_vec2] = geweke_moments(param_draws2,options_);
results_struct = geweke_chi2_test(results_vec1,results_vec2,results_struct,options_);
eval(['oo_.convergence.geweke.',param_name(jj,:),'=results_struct;'])
datamat(jj,:)=[results_struct.posteriormean,results_struct.posteriorstd,results_struct.prob_chi2_test];
@ -190,22 +190,22 @@ if nblck == 1 % Brooks and Gelman tests need more than one block
headers = char(Geweke_tex_header);
lh = size(param_name_tex,2)+2;
my_title=sprintf('Geweke (1992) Convergence Tests, based on means of draws %d to %d vs %d to %d. p-values are for $\\\\chi^2$-test for equality of means.',first_obs_begin_sample,last_obs_begin_sample,first_obs_end_sample,NumberOfDraws);
dyn_latex_table(M_,options_,my_title,'geweke',headers,param_name_tex,datamat,lh,12,4,additional_header);
dyn_latex_table(M_,options_,my_title,'geweke',headers,param_name_tex,datamat,lh,12,4,additional_header);
end
skipline(2);
if options_.convergence.rafterylewis.indicator
if any(options_.convergence.rafterylewis.qrs<0) || any(options_.convergence.rafterylewis.qrs>1) || length(options_.convergence.rafterylewis.qrs)~=3 ...
|| (options_.convergence.rafterylewis.qrs(1)-options_.convergence.rafterylewis.qrs(2)<=0)
fprintf('\nCONVERGENCE DIAGNOSTICS: Invalid option for raftery_lewis_qrs. Using the default of [0.025 0.005 0.95].\n')
options_.convergence.rafterylewis.qrs=[0.025 0.005 0.95];
end
end
Raftery_Lewis_q=options_.convergence.rafterylewis.qrs(1);
Raftery_Lewis_r=options_.convergence.rafterylewis.qrs(2);
Raftery_Lewis_s=options_.convergence.rafterylewis.qrs(3);
oo_.Raftery_Lewis = raftery_lewis(x2,Raftery_Lewis_q,Raftery_Lewis_r,Raftery_Lewis_s);
oo_.Raftery_Lewis.parameter_names=param_name;
my_title=sprintf('Raftery/Lewis (1992) Convergence Diagnostics, based on quantile q=%4.3f with precision r=%4.3f with probability s=%4.3f.',Raftery_Lewis_q,Raftery_Lewis_r,Raftery_Lewis_s);
my_title=sprintf('Raftery/Lewis (1992) Convergence Diagnostics, based on quantile q=%4.3f with precision r=%4.3f with probability s=%4.3f.',Raftery_Lewis_q,Raftery_Lewis_r,Raftery_Lewis_s);
headers = char('Variables','M (burn-in)','N (req. draws)','N+M (total draws)','k (thinning)');
raftery_data_mat=[oo_.Raftery_Lewis.M_burn,oo_.Raftery_Lewis.N_prec,oo_.Raftery_Lewis.N_total,oo_.Raftery_Lewis.k_thin];
@ -217,15 +217,15 @@ if nblck == 1 % Brooks and Gelman tests need more than one block
labels_Raftery_Lewis_tex=char(param_name_tex,'Maximum');
lh = size(labels_Raftery_Lewis_tex,2)+2;
dyn_latex_table(M_,options_,my_title,'raftery_lewis',headers,labels_Raftery_Lewis_tex,raftery_data_mat,lh,10,0);
end
end
end
return;
return
end
Origin = 1000;
StepSize = ceil((NumberOfDraws-Origin)/100);% So that the computational time does not
ALPHA = 0.2; % increase too much with the number of simulations.
StepSize = ceil((NumberOfDraws-Origin)/100);% So that the computational time does not
ALPHA = 0.2; % increase too much with the number of simulations.
time = 1:NumberOfDraws;
xx = Origin:StepSize:NumberOfDraws;
NumberOfLines = length(xx);
@ -262,7 +262,7 @@ localVars.M_ = M_;
% Like sequential execution!
if isnumeric(options_.parallel),
if isnumeric(options_.parallel)
fout = McMCDiagnostics_core(localVars,1,npar,0);
UDIAG = fout.UDIAG;
clear fout
@ -273,19 +273,19 @@ else
ModelName = [ModelName '_bvar'];
end
NamFileInput={[M_.dname '/metropolis/'],[ModelName '_mh*_blck*.mat']};
[fout, nBlockPerCPU, totCPU] = masterParallel(options_.parallel, 1, npar,NamFileInput,'McMCDiagnostics_core', localVars, [], options_.parallel_info);
UDIAG = fout(1).UDIAG;
for j=2:totCPU,
for j=2:totCPU
UDIAG = cat(3,UDIAG ,fout(j).UDIAG);
end
end
UDIAG(:,[2 4 6],:) = UDIAG(:,[2 4 6],:)/nblck;
skipline()
clear pmet temp moyenne CSUP CINF csup cinf n linea iter tmp;
clear pmet temp moyenne CSUP CINF csup cinf n linea iter tmp;
pages = floor(npar/3);
k = 0;
k = 0;
for i = 1:pages
h=dyn_figure(options_.nodisplay,'Name','MCMC univariate convergence diagnostic (Brooks and Gelman,1998)');
boxplot = 1;
@ -296,12 +296,12 @@ for i = 1:pages
if crit == 1
plt1 = UDIAG(:,1,k);
plt2 = UDIAG(:,2,k);
namnam = [nam , ' (Interval)'];
namnam = [nam , ' (Interval)'];
elseif crit == 2
plt1 = UDIAG(:,3,k);
plt2 = UDIAG(:,4,k);
namnam = [nam , ' (m2)'];
elseif crit == 3
elseif crit == 3
plt1 = UDIAG(:,5,k);
plt2 = UDIAG(:,6,k);
namnam = [nam , ' (m3)'];
@ -330,7 +330,7 @@ for i = 1:pages
fprintf(fidTeX,'\\begin{figure}[H]\n');
for jj = 1:size(NAMES,1)
fprintf(fidTeX,'\\psfrag{%s}[1][][0.5][0]{%s}\n',deblank(NAMES(jj,:)),deblank(TEXNAMES(jj,:)));
end
end
fprintf(fidTeX,'\\centering \n');
fprintf(fidTeX,'\\includegraphics[width=%2.2f\\textwidth]{%s_udiag%s}\n',options_.figures.textwidth*min((boxplot-1)/3,1),[OutputFolder '/' ModelName],int2str(i));
fprintf(fidTeX,'\\caption{Univariate convergence diagnostics for the Metropolis-Hastings.\n');
@ -346,7 +346,7 @@ if reste
if reste == 1
nr = 3;
nc = 1;
elseif reste == 2;
elseif reste == 2
nr = 2;
nc = 3;
end
@ -359,12 +359,12 @@ if reste
if crit == 1
plt1 = UDIAG(:,1,k);
plt2 = UDIAG(:,2,k);
namnam = [nam , ' (Interval)'];
namnam = [nam , ' (Interval)'];
elseif crit == 2
plt1 = UDIAG(:,3,k);
plt2 = UDIAG(:,4,k);
namnam = [nam , ' (m2)'];
elseif crit == 3
elseif crit == 3
plt1 = UDIAG(:,5,k);
plt2 = UDIAG(:,6,k);
namnam = [nam , ' (m3)'];
@ -391,9 +391,9 @@ if reste
dyn_saveas(h,[ OutputFolder '/' ModelName '_udiag' int2str(pages+1)],options_.nodisplay,options_.graph_format);
if TeX && any(strcmp('eps',cellstr(options_.graph_format)))
fprintf(fidTeX,'\\begin{figure}[H]\n');
for jj = 1:size(NAMES,1);
for jj = 1:size(NAMES,1)
fprintf(fidTeX,'\\psfrag{%s}[1][][0.5][0]{%s}\n',deblank(NAMES(jj,:)),deblank(TEXNAMES(jj,:)));
end
end
fprintf(fidTeX,'\\centering \n');
fprintf(fidTeX,'\\includegraphics[width=%2.2f\\textwidth]{%s_udiag%s}\n',options_.figures.textwidth*min((boxplot-1)/nc,1),[OutputFolder '/' ModelName],int2str(pages+1));
if reste == 2
@ -435,7 +435,7 @@ for b = 1:nblck
end
clear logpo2;
tmp(:,2) = kron(transpose(1:nblck),ones(NumberOfDraws,1));
tmp(:,3) = kron(ones(nblck,1),time');
tmp(:,3) = kron(ones(nblck,1),time');
tmp = sortrows(tmp,1);
ligne = 0;
for iter = Origin:StepSize:NumberOfDraws
@ -454,12 +454,12 @@ for iter = Origin:StepSize:NumberOfDraws
for i=1:nblck
pmet = temp(find(temp(:,2)==i));
MDIAG(ligne,2) = MDIAG(ligne,2) + pmet(csup,1)-pmet(cinf,1);
moyenne = mean(pmet,1); %% Within mean.
moyenne = mean(pmet,1); %% Within mean.
MDIAG(ligne,4) = MDIAG(ligne,4) + sum((pmet(:,1)-moyenne).^2)/(n-1);
MDIAG(ligne,6) = MDIAG(ligne,6) + sum(abs(pmet(:,1)-moyenne).^3)/(n-1);
end
end
MDIAG(:,[2 4 6],:) = MDIAG(:,[2 4 6],:)/nblck;
MDIAG(:,[2 4 6],:) = MDIAG(:,[2 4 6],:)/nblck;
h = dyn_figure(options_.nodisplay,'Name','Multivariate convergence diagnostic');
boxplot = 1;
@ -467,12 +467,12 @@ for crit = 1:3
if crit == 1
plt1 = MDIAG(:,1);
plt2 = MDIAG(:,2);
namnam = 'Interval';
namnam = 'Interval';
elseif crit == 2
plt1 = MDIAG(:,3);
plt2 = MDIAG(:,4);
namnam = 'm2';
elseif crit == 3
elseif crit == 3
plt1 = MDIAG(:,5);
plt2 = MDIAG(:,6);
namnam = 'm3';
@ -499,7 +499,7 @@ if TeX && any(strcmp('eps',cellstr(options_.graph_format)))
fprintf(fidTeX,'\\begin{figure}[H]\n');
for jj = 1:3
fprintf(fidTeX,'\\psfrag{%s}[1][][0.5][0]{%s}\n',deblank(NAMES(jj,:)),' ');
end
end
fprintf(fidTeX,'\\centering \n');
fprintf(fidTeX,'\\includegraphics[width=0.8\\textwidth]{%s_mdiag}\n',[OutputFolder '/' ModelName]);
fprintf(fidTeX,'\\caption{Multivariate convergence diagnostics for the Metropolis-Hastings.\n');
@ -512,4 +512,3 @@ if TeX && any(strcmp('eps',cellstr(options_.graph_format)))
fprintf(fidTeX,'%% End Of TeX file.');
fclose(fidTeX);
end

View File

@ -2,12 +2,12 @@ function myoutput = McMCDiagnostics_core(myinputs,fpar,npar,whoiam, ThisMatlab)
% function myoutput = McMCDiagnostics_core(myinputs,fpar,npar,whoiam, ThisMatlab)
% Computes the Brooks/Gelman (1998) convergence diagnostics, both the
% parameteric and the non-parameteric versions
%
%
% PARALLEL CONTEXT
% Core functionality for MCMC Diagnostics, which can be parallelized.
% See also the comment in posterior_sampler_core.m funtion.
%
%
%
%
% INPUTS
% See See the comment in posterior_sampler_core.m funtion.
@ -25,15 +25,15 @@ function myoutput = McMCDiagnostics_core(myinputs,fpar,npar,whoiam, ThisMatlab)
%
% ALGORITHM
% Computes part of the convergence diagnostics, the rest is computed in McMCDiagnostics.m .
% The methodology and terminology is based on: Brooks/Gelman (1998): General
% Methods for Monitoring Convergence of Iterative Simulations, Journal of Computational
% The methodology and terminology is based on: Brooks/Gelman (1998): General
% Methods for Monitoring Convergence of Iterative Simulations, Journal of Computational
% and Graphical Statistics, Volume 7, Number 4, Pages 434-455
%
%
%
% SPECIAL REQUIREMENTS.
% None.
% Copyright (C) 2006-2016 Dynare Team
% Copyright (C) 2006-2017 Dynare Team
%
% This file is part of Dynare.
%
@ -50,7 +50,7 @@ function myoutput = McMCDiagnostics_core(myinputs,fpar,npar,whoiam, ThisMatlab)
% 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<4,
if nargin<4
whoiam=0;
end
@ -71,7 +71,7 @@ M_=myinputs.M_;
if whoiam
Parallel=myinputs.Parallel;
end
if ~exist('MetropolisFolder'),
if ~exist('MetropolisFolder')
MetropolisFolder = CheckPath('metropolis',M_.dname);
end
@ -81,16 +81,16 @@ UDIAG = zeros(NumberOfLines,6,npar-fpar+1);
if whoiam
waitbarString = ['Please wait... McMCDiagnostics (' int2str(fpar) 'of' int2str(npar) ')...'];
if Parallel(ThisMatlab).Local,
if Parallel(ThisMatlab).Local
waitbarTitle=['Local '];
else
waitbarTitle=[Parallel(ThisMatlab).ComputerName];
end
fMessageStatus(0,whoiam,waitbarString, waitbarTitle, Parallel(ThisMatlab));
end
for j=fpar:npar,
for j=fpar:npar
if isoctave
if (whoiam==0),
if (whoiam==0)
printf(' Parameter %d... ',j);
end
else
@ -131,13 +131,13 @@ for j=fpar:npar,
end
end
if isoctave
if (whoiam==0),
if (whoiam==0)
printf('Done! \n');
end
else
fprintf('Done! \n');
end
if whoiam,
if whoiam
waitbarString = [ 'Parameter ' int2str(j) '/' int2str(npar) ' done.'];
fMessageStatus((j-fpar+1)/(npar-fpar+1),whoiam,waitbarString, waitbarTitle, Parallel(ThisMatlab))
end

View File

@ -2,7 +2,7 @@ function results_struct = geweke_chi2_test(results1,results2,results_struct,opti
% results_struct = geweke_chi2_test(results1,results2,results_struct,options)
% PURPOSE: computes Geweke's chi-squared test for two sets of MCMC sample draws
%
% INPUTS
% INPUTS
% results1 [1 by (4+n_taper*2) vector] vector with post. mean,
% std, NSE_iid, RNE_iid, and tapered NSE and RNE
% for chain part 1
@ -12,7 +12,7 @@ function results_struct = geweke_chi2_test(results1,results2,results_struct,opti
% results_struct [structure] results structure generated by geweke_moments
% Dynareoptions [structure]
%
% OUTPUTS
% OUTPUTS
% results_struct [structure] containing the following fields:
% pooled_mean Pooled mean of the chain parts, weighted
% with precision
@ -26,7 +26,7 @@ function results_struct = geweke_chi2_test(results1,results2,results_struct,opti
% SPECIAL REQUIREMENTS
% None.
% Copyright (C) 2013 Dynare Team
% Copyright (C) 2013-2017 Dynare Team
%
% This file is part of Dynare.
%
@ -49,26 +49,24 @@ function results_struct = geweke_chi2_test(results1,results2,results_struct,opti
% J.M. Bernardo, A.P. Dawid, and A.F.M. Smith (eds.) Proceedings of
% the Fourth Valencia International Meeting on Bayesian Statistics,
% pp. 169-194, Oxford University Press
% Geweke (1999): `Using simulation methods for Bayesian econometric models:
% Geweke (1999): `Using simulation methods for Bayesian econometric models:
% Inference, development and communication', Econometric Reviews, 18(1),
% 1-73
% written by: Johannes Pfeifer,
% based on code by James P. LeSage, who in turn
% drew on MATLAB programs written by Siddartha Chib
% written by: Johannes Pfeifer,
% based on code by James P. LeSage, who in turn
% drew on MATLAB programs written by Siddartha Chib
for k=1:length(options.convergence.geweke.taper_steps)+1;
for k=1:length(options.convergence.geweke.taper_steps)+1
NSE=[results1(:,3+(k-1)*2) results2(:,3+(k-1)*2)];
means=[results1(:,1) results2(:,1)];
diff_Means=means(:,1)-means(:,2);
sum_of_weights=sum(1./(NSE.^2),2);
pooled_mean=sum(means./(NSE.^2),2)./sum_of_weights;
pooled_NSE=1./sqrt(sum_of_weights);
test_stat=diff_Means.^2./sum(NSE.^2,2);
test_stat=diff_Means.^2./sum(NSE.^2,2);
p = 1-chi2cdf(test_stat,1);
results_struct.pooled_mean(:,k) = pooled_mean;
results_struct.pooled_nse(:,k) = pooled_NSE;
results_struct.prob_chi2_test(:,k) = p;
end;
end

View File

@ -1,13 +1,13 @@
function [results_vec, results_struct] = geweke_moments(draws,Dynareoptions)
%[results_vec, results_struct] = geweke_moments(draws,Dynareoptions)
% PURPOSE: computes Gewke's convergence diagnostics NSE and RNE
% PURPOSE: computes Gewke's convergence diagnostics NSE and RNE
% (numerical std error and relative numerical efficiencies)
% INPUTS
% draws [ndraws by 1 vector]
% INPUTS
% draws [ndraws by 1 vector]
% Dynareoptions [structure]
%
% OUTPUTS
%
% OUTPUTS
% results_vec
% results_struct [structure] containing the following fields:
% posteriormean= posterior parameter mean
@ -22,7 +22,7 @@ function [results_vec, results_struct] = geweke_moments(draws,Dynareoptions)
% SPECIAL REQUIREMENTS
% None.
% Copyright (C) 2013 Dynare Team
% Copyright (C) 2013-2017 Dynare Team
%
% This file is part of Dynare.
%
@ -44,14 +44,14 @@ function [results_vec, results_struct] = geweke_moments(draws,Dynareoptions)
% J.M. Bernardo, A.P. Dawid, and A.F.M. Smith (eds.) Proceedings of
% the Fourth Valencia International Meeting on Bayesian Statistics,
% pp. 169-194, Oxford University Press
% Geweke (1999): `Using simulation methods for Bayesian econometric models:
% Geweke (1999): `Using simulation methods for Bayesian econometric models:
% Inference, development and communication', Econometric Reviews, 18(1),
% 1-73
% -----------------------------------------------------------------
% written by: Johannes Pfeifer,
% based on code by James P. LeSage, who in turn
% drew on MATLAB programs written by Siddartha Chib
% written by: Johannes Pfeifer,
% based on code by James P. LeSage, who in turn
% drew on MATLAB programs written by Siddartha Chib
ndraw = size(draws,1);
@ -64,10 +64,10 @@ n_draws_used = ns*n_groups; %effective number of draws used after rounding down
window_means= zeros(n_groups,1);
window_uncentered_variances= zeros(n_groups,1);
for ig=1:n_groups;
for ig=1:n_groups
window_means(ig,1)=sum(draws((ig-1)*ns+1:ig*ns,1))/ns;
window_uncentered_variances(ig,1)=sum(draws((ig-1)*ns+1:ig*ns,1).^2)/ns;
end; %for ig
window_uncentered_variances(ig,1)=sum(draws((ig-1)*ns+1:ig*ns,1).^2)/ns;
end %for ig
total_mean=mean(window_means);
total_variance=mean(window_uncentered_variances)-total_mean^2;
@ -88,9 +88,9 @@ results_struct.rne_iid = results_vec(1,4);
%get autocovariance of grouped means
centered_window_means=window_means-total_mean;
autocov_grouped_means=zeros(n_groups,1);
for lag=0:n_groups-1;
for lag=0:n_groups-1
autocov_grouped_means(lag+1)=centered_window_means(lag+1:n_groups,1)'*centered_window_means(1:n_groups-lag,1)/100;
end;
end
% numerical standard error with tapered autocovariance functions
for taper_index=1:length(taper_steps)
@ -105,5 +105,4 @@ for taper_index=1:length(taper_steps)
eval(['results_struct.nse_taper_',num2str(taper),'= NSE_taper;']);
eval(['results_struct.rne_taper_',num2str(taper),'= total_variance/(n_draws_used*NSE_taper^2);']);
end; % end of for mm loop
end % end of for mm loop

View File

@ -14,12 +14,12 @@ function Ifac = mcmc_ifac(X, Nc)
% ALGORITHM:
% Inefficiency factors are computed as
% \[
% Ifac = 1 + 2\sum\limits_{i=1}^{Nc} {\hat \rho(i)}
% Ifac = 1 + 2\sum\limits_{i=1}^{Nc} {\hat \rho(i)}
% \]
% where $\hat \rho(i)$ denotes the autocorrelation at lag i and the terms
% of the sum are truncated using a Parzen window.
%
% For inefficiency factors, see Section 6.1 of Paolo Giordani, Michael Pitt, and Robert Kohn (2011):
%
% For inefficiency factors, see Section 6.1 of Paolo Giordani, Michael Pitt, and Robert Kohn (2011):
% "Bayesian Inference for Time Series State Space Models" in : John Geweke, Gary Koop,
% Herman van Dijk (editors): "The Oxford Handbook of Bayesian
% Econometrics", Oxford University Press
@ -36,7 +36,7 @@ function Ifac = mcmc_ifac(X, Nc)
% consistent covariance matrix estimation", Econometrica, 59(3), p. 817-858
% Copyright (C) 2015-16 Dynare Team
% Copyright (C) 2015-2017 Dynare Team
%
% This file is part of Dynare.
%
@ -54,7 +54,7 @@ function Ifac = mcmc_ifac(X, Nc)
% along with Dynare. If not, see <http://www.gnu.org/licenses/>.
Nc = floor(min(Nc, length(X)/2));
if mod(Nc,2),
if mod(Nc,2)
Nc=Nc-1;
end
AcorrXSIM = dyn_autocorr(X(:), Nc);

View File

@ -1,15 +1,15 @@
function [raftery_lewis] = raftery_lewis(runs,q,r,s)
% function raftery_lewis = raftery_lewis(runs,q,r,s)
% Computes the convergence diagnostics of Raftery and Lewis (1992), i.e. the
% number of draws needed in MCMC to estimate the posterior cdf of the q-quantile
% Computes the convergence diagnostics of Raftery and Lewis (1992), i.e. the
% number of draws needed in MCMC to estimate the posterior cdf of the q-quantile
% within an accuracy r with probability s
%
%
% Inputs:
% - draws [n_draws by n_var] double matrix of draws from the sampler
% - draws [n_draws by n_var] double matrix of draws from the sampler
% - q [scalar] quantile of the quantity of interest
% - r [scalar] level of desired precision
% - r [scalar] level of desired precision
% - s [scalar] probability associated with r
%
%
% Output:
% raftery_lewis [structure] containing the fields:
% - M_burn [n_draws by 1] number of draws required for burn-in
@ -21,7 +21,7 @@ function [raftery_lewis] = raftery_lewis(runs,q,r,s)
% iterations due to dependence in chain
% - N_min [scalar] # draws if the chain is white noise
% - N_total [n_draws by 1] nburn + nprec
%
%
% ---------------------------------------------------------------------
% NOTES: Example values of q, r, s:
@ -30,23 +30,23 @@ function [raftery_lewis] = raftery_lewis(runs,q,r,s)
%
% - The result is quite sensitive to r, being proportional to the
% inverse of r^2.
% - For epsilon (closeness of probabilities to equilibrium values),
% - For epsilon (closeness of probabilities to equilibrium values),
% Raftery/Lewis use 0.001 and argue that the results
% are quite robust to changes in this value
%
% ---------------------------------------------------------------------
% REFERENCES:
% REFERENCES:
% Raftery, Adrien E./Lewis, Steven (1992a): "How many iterations in the Gibbs sampler?"
% in: Bernardo/Berger/Dawid/Smith (eds.): Bayesian Statistics, Vol. 4, Clarendon Press: Oxford,
% pp. 763-773.
% Raftery, Adrien E./Lewis, Steven (1992b): "Comment: One long run with diagnostics:
% Implementation strategies for Markov chain Monte Carlo." Statistical Science,
% 7(4), pp. 493-497.
% in: Bernardo/Berger/Dawid/Smith (eds.): Bayesian Statistics, Vol. 4, Clarendon Press: Oxford,
% pp. 763-773.
% Raftery, Adrien E./Lewis, Steven (1992b): "Comment: One long run with diagnostics:
% Implementation strategies for Markov chain Monte Carlo." Statistical Science,
% 7(4), pp. 493-497.
%
% ----------------------------------------------------
% Copyright (C) 2016 Benjamin Born and Johannes Pfeifer
% Copyright (C) 2016 Dynare Team
% Copyright (C) 2016-2017 Dynare Team
%
% This file is part of Dynare.
%
@ -80,42 +80,42 @@ thinned_chain = zeros(n_runs,1);
Phi = norminv((s+1)/2); %note the missing ^{-1} at the Phi in equation top page 5, see RL (1995)
raftery_lewis.N_min = fix(Phi^2*(1-q)*q/r^2+1);
for nv = 1:n_vars % big loop over variables
for nv = 1:n_vars % big loop over variables
if q > 0 && q < 1
work = (runs(:,nv) <= quantile(runs(:,nv),q));
else
error('Quantile must be between 0 and 1');
end;
k_thin_current_var = 1;
bic = 1;
error('Quantile must be between 0 and 1');
end
k_thin_current_var = 1;
bic = 1;
epss = 0.001;
% Find thinning factor for which first-order Markov Chain is preferred to second-order one
while(bic > 0)
thinned_chain=work(1:k_thin_current_var:n_runs,1);
thinned_chain=work(1:k_thin_current_var:n_runs,1);
[g2, bic] = first_vs_second_order_MC_test(thinned_chain);
k_thin_current_var = k_thin_current_var+1;
end;
end
k_thin_current_var = k_thin_current_var-1; %undo last step
%compute transition probabilities
%compute transition probabilities
transition_matrix = zeros(2,2);
for i1 = 2:size(thinned_chain,1)
transition_matrix(thinned_chain(i1-1)+1,thinned_chain(i1)+1) = transition_matrix(thinned_chain(i1-1)+1,thinned_chain(i1)+1)+1;
end;
end
alpha = transition_matrix(1,2)/(transition_matrix(1,1)+transition_matrix(1,2)); %prob of going from 1 to 2
beta = transition_matrix(2,1)/(transition_matrix(2,1)+transition_matrix(2,2)); %prob of going from 2 to 1
kmind=k_thin_current_var;
[g2, bic]=independence_chain_test(thinned_chain);
while(bic > 0)
thinned_chain=work(1:kmind:n_runs,1);
thinned_chain=work(1:kmind:n_runs,1);
[g2, bic] = independence_chain_test(thinned_chain);
kmind = kmind+1;
end;
end
m_star = log((alpha + beta)*epss/max(alpha,beta))/log(abs(1 - alpha - beta)); %equation bottom page 4
raftery_lewis.M_burn(nv) = fix((m_star+1)*k_thin_current_var);
n_star = (2 - (alpha + beta))*alpha*beta*(Phi^2)/((alpha + beta)^3 * r^2); %equation top page 5
@ -124,18 +124,18 @@ for nv = 1:n_vars % big loop over variables
raftery_lewis.k_ind(nv) = max(fix(raftery_lewis.I_stat(nv)+1),kmind);
raftery_lewis.k_thin(nv) = k_thin_current_var;
raftery_lewis.N_total(nv)= raftery_lewis.M_burn(nv)+raftery_lewis.N_prec(nv);
end;
end
end
function [g2, bic] = first_vs_second_order_MC_test(d)
%conducts a test of first vs. second order Markov Chain via BIC criterion
n_obs=size(d,1);
g2 = 0;
g2 = 0;
tran=zeros(2,2,2);
for t_iter=3:n_obs % count state transitions
tran(d(t_iter-2,1)+1,d(t_iter-1,1)+1,d(t_iter,1)+1)=tran(d(t_iter-2,1)+1,d(t_iter-1,1)+1,d(t_iter,1)+1)+1;
end;
end
% Compute the log likelihood ratio statistic for second-order MC vs first-order MC. G2 statistic of Bishop, Fienberg and Holland (1975)
for ind_1 = 1:2
for ind_2 = 1:2
@ -146,9 +146,9 @@ for ind_1 = 1:2
focus = tran(ind_1,ind_2,ind_3);
g2 = g2 + log(focus/fitted)*focus;
end
end; % end of for i3
end; % end of for i2
end; % end of for i1
end % end of for i3
end % end of for i2
end % end of for i1
g2 = g2*2;
bic = g2 - log(n_obs-2)*2;
@ -161,19 +161,19 @@ n_obs=size(d,1);
trans = zeros(2,2);
for ind_1 = 2:n_obs
trans(d(ind_1-1)+1,d(ind_1)+1)=trans(d(ind_1-1)+1,d(ind_1)+1)+1;
end;
dcm1 = n_obs - 1;
end
dcm1 = n_obs - 1;
g2 = 0;
% Compute the log likelihood ratio statistic for second-order MC vs first-order MC. G2 statistic of Bishop, Fienberg and Holland (1975)
for ind_1 = 1:2
for ind_2 = 1:2
if trans(ind_1,ind_2) ~= 0
fitted = ((trans(ind_1,1) + trans(ind_1,2))*(trans(1,ind_2) + trans(2,ind_2)))/dcm1;
fitted = ((trans(ind_1,1) + trans(ind_1,2))*(trans(1,ind_2) + trans(2,ind_2)))/dcm1;
focus = trans(ind_1,ind_2);
g2 = g2 + log(focus/fitted)*focus;
end;
end;
end;
g2 = g2*2;
end
end
end
g2 = g2*2;
bic = g2 - log(dcm1);
end

View File

@ -19,7 +19,7 @@ function [info] = convertAimCodeToInfo(aimCode)
% OUTPUTS
% info [integer] Code to be used to print error in print_info.m
% Copyright (C) 2011-2012 Dynare Team
% Copyright (C) 2011-2017 Dynare Team
%
% This file is part of Dynare.
%
@ -55,7 +55,7 @@ switch aimCode
info = 161;
case 62
info = 162;
case 63;
case 63
info = 163;
case 64
info = 164;

View File

@ -13,7 +13,7 @@ function oo_ = convert_dyn_45_to_44(M_, options_, oo_,bayestopt_)
% SPECIAL REQUIREMENTS
% none
% Copyright (C) 2015-2016 Dynare Team
% Copyright (C) 2015-2017 Dynare Team
%
% This file is part of Dynare.
%
@ -40,7 +40,7 @@ if isfield(oo_,'PointForecast')
oo_.MeanForecast.(moment_names{moment_iter}).(var_names{var_iter})=...
[oo_.SmoothedVariables.(moment_names{moment_iter}).(var_names{var_iter})(:,end)*ones(M_.maximum_endo_lag,1) oo_.MeanForecast.(moment_names{moment_iter}).(var_names{var_iter})];
oo_.PointForecast.(moment_names{moment_iter}).(var_names{var_iter})=...
[oo_.SmoothedVariables.(moment_names{moment_iter}).(var_names{var_iter})(:,end)*ones(M_.maximum_endo_lag,1) oo_.PointForecast.(moment_names{moment_iter}).(var_names{var_iter})];
[oo_.SmoothedVariables.(moment_names{moment_iter}).(var_names{var_iter})(:,end)*ones(M_.maximum_endo_lag,1) oo_.PointForecast.(moment_names{moment_iter}).(var_names{var_iter})];
else
oo_.MeanForecast.(moment_names{moment_iter}).(var_names{var_iter})=...
[oo_.SmoothedVariables.(moment_names{moment_iter}).(var_names{var_iter})(end)*ones(M_.maximum_endo_lag,1); oo_.MeanForecast.(moment_names{moment_iter}).(var_names{var_iter})];
@ -106,7 +106,7 @@ if isfield(oo_,'UpdatedVariables')
for ii=1:length(names)
%make sure Bayesian fields are not affected
if ~strcmp(names{ii},'Mean') && ~strcmp(names{ii},'Median') && ~strcmp(names{ii},'deciles') ...
&& ~strcmp(names{ii},'Var') && ~strcmp(names{ii},'HPDinf') && ~strcmp(names{ii},'HPDsup')
&& ~strcmp(names{ii},'Var') && ~strcmp(names{ii},'HPDinf') && ~strcmp(names{ii},'HPDsup')
current_var_index=find(strmatch(names{ii},deblank(M_.endo_names),'exact'));
if options_.loglinear == 1 %logged steady state must be used
constant_current_variable=log(oo_.dr.ys(current_var_index));
@ -127,7 +127,7 @@ if isfield(oo_,'FilteredVariables')
for ii=1:length(names)
%make sure Bayesian fields are not affected
if ~strcmp(names{ii},'Mean') && ~strcmp(names{ii},'Median') && ~strcmp(names{ii},'deciles') ...
&& ~strcmp(names{ii},'Var') && ~strcmp(names{ii},'HPDinf') && ~strcmp(names{ii},'HPDsup')
&& ~strcmp(names{ii},'Var') && ~strcmp(names{ii},'HPDinf') && ~strcmp(names{ii},'HPDsup')
current_var_index=find(strmatch(names{ii},deblank(M_.endo_names),'exact'));
if options_.loglinear == 1 %logged steady state must be used
constant_current_variable=log(oo_.dr.ys(current_var_index));

View File

@ -15,7 +15,7 @@ function oo_ = convert_oo_(M_, options_, oo_, from_ver, to_ver)
% SPECIAL REQUIREMENTS
% none
% Copyright (C) 2015 Dynare Team
% Copyright (C) 2015-2017 Dynare Team
%
% This file is part of Dynare.
%
@ -63,7 +63,7 @@ else
end
if strcmp(from_ver, to_ver)
return;
return
end
if ver_greater_than(to_ver, from_ver)

View File

@ -2,7 +2,7 @@ function oo_ = correlation_mc_analysis(SampleSize,type,dname,fname,vartan,nvar,v
% This function analyses the (posterior or prior) distribution of the
% endogenous variables correlation function.
% Copyright (C) 2008-2013 Dynare Team
% Copyright (C) 2008-2017 Dynare Team
%
% This file is part of Dynare.
%
@ -137,7 +137,7 @@ end
function oo_ = fill_output_structure(var1,var2,type,oo_,moment,lag,result)
switch moment
case {'Mean','Median','Variance','HPDinf','HPDsup'}
case {'Mean','Median','Variance','HPDinf','HPDsup'}
oo_.([type, 'TheoreticalMoments']).dsge.correlation.(moment).(var1).(var2)(lag,1) = result;
case {'deciles','density'}
oo_.([type, 'TheoreticalMoments']).dsge.correlation.(moment).(var1).(var2)(lag,1) = {result};

View File

@ -1,4 +1,4 @@
function [co, b, yhat] = cosn(H);
function [co, b, yhat] = cosn(H)
% function co = cosn(H);
% computes the cosine of the angle between the H(:,1) and its
@ -7,7 +7,7 @@ function [co, b, yhat] = cosn(H);
% Not the same as multiple correlation coefficient since the means are not
% zero
%
% Copyright (C) 2008-2012 Dynare Team
% Copyright (C) 2008-2017 Dynare Team
%
% This file is part of Dynare.
%
@ -28,15 +28,12 @@ y = H(:,1);
X = H(:,2:end);
b=(X\y);
if any(isnan(b)) || any(isinf(b)),
if any(isnan(b)) || any(isinf(b))
b=0;
end
yhat = X*b;
if rank(yhat),
if rank(yhat)
co = abs(y'*yhat/sqrt((y'*y)*(yhat'*yhat)));
else
co=0;
end

View File

@ -1,7 +1,7 @@
function oo_ = covariance_mc_analysis(NumberOfSimulations,type,dname,fname,vartan,nvar,var1,var2,mh_conf_sig,oo_,options_)
% This function analyses the (posterior or prior) distribution of the
% endogenous variables' covariance matrix.
%
%
% INPUTS
% NumberOfSimulations [integer] scalar, number of simulations.
% type [string] 'prior' or 'posterior'
@ -14,12 +14,12 @@ function oo_ = covariance_mc_analysis(NumberOfSimulations,type,dname,fname,varta
% mh_conf_sig [double] 2 by 1 vector with upper
% and lower bound of HPD intervals
% oo_ [structure] Dynare structure where the results are saved.
% options_ [structure] Dynare options structure
% options_ [structure] Dynare options structure
%
% OUTPUTS
% oo_ [structure] Dynare structure where the results are saved.
% Copyright (C) 2008-2015 Dynare Team
% Copyright (C) 2008-2017 Dynare Team
%
% This file is part of Dynare.
%
@ -64,7 +64,7 @@ var1=deblank(var1);
var2=deblank(var2);
if isfield(oo_,[ TYPE 'TheoreticalMoments'])
temporary_structure = oo_.([TYPE, 'TheoreticalMoments']);
temporary_structure = oo_.([TYPE, 'TheoreticalMoments']);
if isfield(temporary_structure,'dsge')
temporary_structure = oo_.([TYPE, 'TheoreticalMoments']).dsge;
if isfield(temporary_structure,'covariance')
@ -104,7 +104,7 @@ for file = 1:length(ListOfFiles)
temp=Covariance_matrix(:,cov_pos)./(sqrt(Covariance_matrix(:,var_pos_1)).*sqrt(Covariance_matrix(:,var_pos_2)));
temp(Covariance_matrix(:,cov_pos)==0)=0; %filter out 0 correlations that would result in 0/0
tmp_corr_mat(i1:i2)=temp;
end
end
i1 = i2+1;
end

View File

@ -12,7 +12,7 @@ function [x,rc] = csolve(FUN,x,gradfun,crit,itmax,varargin)
% the equation is solved.
% itmax: the solver stops when this number of iterations is reached, with rc=4
% varargin: in this position the user can place any number of additional arguments, all
% of which are passed on to FUN and gradfun (when it is non-empty) as a list of
% of which are passed on to FUN and gradfun (when it is non-empty) as a list of
% arguments following x.
% rc: 0 means normal solution, 1 and 3 mean no solution despite extremely fine adjustments
% in step length (very likely a numerical problem, or a discontinuity). 4 means itmax
@ -22,7 +22,7 @@ function [x,rc] = csolve(FUN,x,gradfun,crit,itmax,varargin)
% http://sims.princeton.edu/yftp/optimize/mfiles/csolve.m
% Copyright (C) 1993-2007 Christopher Sims
% Copyright (C) 2007-2011 Dynare Team
% Copyright (C) 2007-2017 Dynare Team
%
% This file is part of Dynare.
%
@ -60,7 +60,7 @@ if isempty(varargin)
f0=feval(FUN,x);
else
f0=feval(FUN,x,varargin{:});
end
end
af0=sum(abs(f0));
af00=af0;
itct=0;
@ -122,7 +122,7 @@ while ~done
factor=factor^.6;
shrink=1;
end
if abs(lambda*(1-factor))*dxSize > .1*delta;
if abs(lambda*(1-factor))*dxSize > .1*delta
lambda = factor*lambda;
elseif (lambda > 0) && (factor==.6) %i.e., we've only been shrinking
lambda=-.3;
@ -162,7 +162,7 @@ while ~done
if itct >= itmax
done=1;
rc=4;
elseif af0<crit;
elseif af0<crit
done=1;
rc=0;
end

View File

@ -29,12 +29,12 @@ function [nodes, weights] = cubature_with_gaussian_weight(d,n,method)
%! The routine returns nodes and associated weights to compute a multivariate integral of the form:
%!
%! \int_D f(x)*\exp(-<x,x>) dx
%!
%!
%!
%!
%! @end deftypefn
%@eod:
% Copyright (C) 2012-2013 Dynare Team
% Copyright (C) 2012-2017 Dynare Team
%
% This file is part of Dynare.
%
@ -90,7 +90,7 @@ if strcmp(method,'Stroud') && isequal(n,5)
nodes = zeros(d,m);
weights = zeros(m,1);
% Set the weight for the first node (0)
weights(1) = A;
weights(1) = A;
skip = 1;
% Set the remaining nodes and associated weights.
nodes(:,skip+(1:d)) = r*eye(d);
@ -130,7 +130,7 @@ m(:,4) = -m(:,1);
%@test:1
%$ % Set problem
%$ d = 4;
%$
%$
%$ t = zeros(5,1);
%$
%$ % Call the tested routine
@ -289,7 +289,7 @@ m(:,4) = -m(:,1);
%@test:5
%$ % Set problem
%$ d = 5;
%$
%$
%$ t = zeros(6,1);
%$
%$ % Call the tested routine
@ -333,7 +333,7 @@ m(:,4) = -m(:,1);
%@test:6
%$ % Set problem
%$ d = 3;
%$
%$
%$ t = zeros(4,1);
%$
%$ % Call the tested routine

View File

@ -43,7 +43,7 @@ function [X, info] = cycle_reduction(A0, A1, A2, cvg_tol, ch) % --*-- Unitary te
%! @end deftypefn
%@eod:
% Copyright (C) 2013 Dynare Team
% Copyright (C) 2013-2017 Dynare Team
%
% This file is part of Dynare.
%
@ -76,7 +76,7 @@ id0 = 1:n;
id2 = id0+n;
cont = 1;
while cont
while cont
tmp = ([A0; A2]/A1)*[A0 A2];
A1 = A1 - tmp(id0,id2) - tmp(id2,id0);
A0 = -tmp(id0,id0);
@ -97,7 +97,7 @@ while cont
info(2) = log(norm(A1,1));
end
return
end
end
it = it + 1;
end

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