Merge branch 'master' into json

time-shift
Stéphane Adjemian (Charybdis) 2017-06-16 20:03:36 +02:00
commit 119b5a62f2
880 changed files with 38591 additions and 29405 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

797
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)
=============================================

View File

@ -91,7 +91,7 @@ If you have downloaded the sources from an official source archive or the source
If you want to use Git, do the following from a terminal:
git clone --recursive http://github.com/DynareTeam/dynare.git
git clone --recursive https://github.com/DynareTeam/dynare.git
cd dynare
autoreconf -si
@ -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
@ -26,8 +26,8 @@
}}
\author{S\'ebastien Villemot\thanks{Paris School of Economics and
CEPREMAP.} \and Johannes Pfeifer\thanks{University of Mannheim. E-mail: \href{mailto:pfeifer@uni-mannheim.de}{\texttt{pfeifer@uni-mannheim.de}}.}}
\date{First version: September 2007 \hspace{1cm} This version: October 2016}
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: 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}

File diff suppressed because it is too large Load Diff

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

@ -12,7 +12,8 @@
#include <vector>
#include <map>
namespace ogp {
namespace ogp
{
class AtomAsgnEvaluator;
@ -21,7 +22,8 @@ namespace ogp {
* all expressions on the right hand side, the parsed formulas of
* the right hand sides, and the information about the left hand
* sides. See documentation to the order member below. */
class AtomAssignings {
class AtomAssignings
{
friend class AtomAsgnEvaluator;
protected:
typedef std::map<const char *, int, ltstr> Tvarintmap;
@ -44,12 +46,14 @@ namespace ogp {
public:
/** Construct the object using the provided static atoms. */
AtomAssignings(StaticAtoms &a) : atoms(a), expr(atoms)
{}
{
}
/** Make a copy with provided reference to (posibly different)
* static atoms. */
AtomAssignings(const AtomAssignings &aa, StaticAtoms &a);
virtual ~AtomAssignings()
{}
{
}
/** Parse the assignments from the given string. */
void parse(int length, const char *stream);
/** Process a syntax error from bison. */
@ -82,7 +86,8 @@ namespace ogp {
class AtomAsgnEvaluator : public FormulaEvalLoader,
public AtomValues,
protected FormulaEvaluator,
public std::vector<double> {
public std::vector<double>
{
protected:
typedef std::map<int, double> Tusrvalmap;
Tusrvalmap user_values;
@ -90,8 +95,12 @@ namespace ogp {
public:
AtomAsgnEvaluator(const AtomAssignings &a)
: FormulaEvaluator(a.expr),
std::vector<double>(a.expr.nformulas()), aa(a) {}
virtual ~AtomAsgnEvaluator() {}
std::vector<double>(a.expr.nformulas()), aa(a)
{
}
virtual ~AtomAsgnEvaluator()
{
}
/** This sets all initial values to NaNs, all constants and
* all values set by user by call set_value. This is called by
* FormulaEvaluator::eval() method, which is called by eval()
@ -113,8 +122,11 @@ namespace ogp {
* result is that this object as std::vector<double> will
* contain the values. It is ordered given by formulas in
* expr. */
void eval()
{FormulaEvaluator::eval(*this, *this);}
void
eval()
{
FormulaEvaluator::eval(*this, *this);
}
/** This returns a value for a given name. If the name is not
* found among atoms, or there is no assignment for the atom,
* NaN is returned. */

View File

@ -9,7 +9,8 @@
#include <string>
namespace ogp {
namespace ogp
{
using std::string;
using std::map;
@ -18,14 +19,18 @@ namespace ogp {
/** This class tracts an information about the performed
* substitutions. In fact, there is only one number to keep track
* about, this is a number of substitutions. */
struct SubstInfo {
struct SubstInfo
{
int num_substs;
SubstInfo() : num_substs(0) {}
SubstInfo() : num_substs(0)
{
}
};
/** This class tracks all atom substitutions during the job and
* then builds structures when all substitutions are finished. */
class AtomSubstitutions {
class AtomSubstitutions
{
public:
typedef pair<const char *, int> Tshiftname;
typedef map<const char *, Tshiftname, ltstr> Tshiftmap;
@ -61,12 +66,16 @@ namespace ogp {
* substitution job is done by a substitution method of the
* new atoms. */
AtomSubstitutions(const FineAtoms &oa, FineAtoms &na)
: old_atoms(oa), new_atoms(na) {}
: old_atoms(oa), new_atoms(na)
{
}
/** Construct a copy of the object using a different instances
* of old atoms and new atoms, which are supposed to be
* semantically same as the atoms from as. */
AtomSubstitutions(const AtomSubstitutions &as, const FineAtoms &oa, FineAtoms &na);
virtual ~AtomSubstitutions() {}
virtual ~AtomSubstitutions()
{
}
/** This is called during the substitution job from the
* substitution method of the new atoms. This says that the
* new name, say "a_m3" is a substitution of old name "a"
@ -81,31 +90,53 @@ namespace ogp {
* substitution, it returns NULL. */
const char *get_new4old(const char *oldname, int tshift) const;
/** Return new2old. */
const Tshiftmap& get_new2old() const
{return new2old;}
const Tshiftmap &
get_new2old() const
{
return new2old;
}
/** Return old2new. */
const Toldnamemap& get_old2new() const
{return old2new;}
const Toldnamemap &
get_old2new() const
{
return old2new;
}
/** Return substitution info. */
const SubstInfo& get_info() const
{return info;}
const SubstInfo &
get_info() const
{
return info;
}
/** Return old atoms. */
const FineAtoms& get_old_atoms() const
{return old_atoms;}
const FineAtoms &
get_old_atoms() const
{
return old_atoms;
}
/** Return new atoms. */
const FineAtoms& get_new_atoms() const
{return new_atoms;}
const FineAtoms &
get_new_atoms() const
{
return new_atoms;
}
/** Debug print. */
void print() const;
};
class SAtoms : public FineAtoms {
class SAtoms : public FineAtoms
{
public:
SAtoms()
: FineAtoms() {}
: FineAtoms()
{
}
SAtoms(const SAtoms &sa)
: FineAtoms(sa) {}
virtual ~SAtoms() {}
: FineAtoms(sa)
{
}
virtual ~SAtoms()
{
}
/** This substitutes all lags and leads for all exogenous and
* all lags and leads greater than 1 for all endogenous
* variables. This is useful for perfect foresight problems
@ -118,12 +149,18 @@ namespace ogp {
protected:
/** This finds an endogenous variable name which occurs between
* ll1 and ll2 included. */
const char* findEndoWithLeadInInterval(int ll1, int ll2) const
{return findNameWithLeadInInterval(get_endovars(), ll1, ll2);}
const char *
findEndoWithLeadInInterval(int ll1, int ll2) const
{
return findNameWithLeadInInterval(get_endovars(), ll1, ll2);
}
/** This finds an exogenous variable name which occurs between
* ll1 and ll2 included. */
const char* findExoWithLeadInInterval(int ll1, int ll2) const
{return findNameWithLeadInInterval(get_exovars(), ll1, ll2);}
const char *
findExoWithLeadInInterval(int ll1, int ll2) const
{
return findNameWithLeadInInterval(get_exovars(), ll1, ll2);
}
/** This attempts to find a non registered name of the form
* <str>_m<abs(ll)> or <str>_p<abs(ll)>. A letter 'p' is

View File

@ -5,15 +5,20 @@
#ifndef OGP_CSV_PARSER
#define OGP_CSV_PARSER
namespace ogp {
namespace ogp
{
class CSVParserPeer {
class CSVParserPeer
{
public:
virtual ~CSVParserPeer() {}
virtual ~CSVParserPeer()
{
}
virtual void item(int irow, int icol, const char *str, int length) = 0;
};
class CSVParser {
class CSVParser
{
private:
CSVParserPeer &peer;
int row;
@ -21,21 +26,36 @@ namespace ogp {
const char *parsed_string;
public:
CSVParser(CSVParserPeer &p)
: peer(p), row(0), col(0), parsed_string(0) {}
: peer(p), row(0), col(0), parsed_string(0)
{
}
CSVParser(const CSVParser &csvp)
: peer(csvp.peer), row(csvp.row),
col(csvp.col), parsed_string(csvp.parsed_string) {}
virtual ~CSVParser() {}
col(csvp.col), parsed_string(csvp.parsed_string)
{
}
virtual ~CSVParser()
{
}
void csv_error(const char *mes);
void csv_parse(int length, const char *str);
void nextrow()
{row++; col = 0;}
void nextcol()
{col++;}
void item(int off, int length)
{peer.item(row, col, parsed_string+off, length);}
void
nextrow()
{
row++; col = 0;
}
void
nextcol()
{
col++;
}
void
item(int off, int length)
{
peer.item(row, col, parsed_string+off, length);
}
};
};

View File

@ -13,15 +13,20 @@
#include <string>
#include <cstring>
namespace ogp {
namespace ogp
{
using std::vector;
using std::map;
using std::set;
using std::string;
struct ltstr {
bool operator()(const char* a1, const char* a2) const
{ return strcmp(a1, a2) < 0; }
struct ltstr
{
bool
operator()(const char *a1, const char *a2) const
{
return strcmp(a1, a2) < 0;
}
};
/** Class storing names. We will keep names of variables in
@ -30,7 +35,8 @@ namespace ogp {
* deallocation. The main function of the class is to allocate
* space for names, and return a pointer of the stored name if
* required. */
class NameStorage {
class NameStorage
{
protected:
/** Vector of names allocated, this is the storage. */
vector<char *> name_store;
@ -38,24 +44,34 @@ namespace ogp {
* allocated or not. */
set<const char *, ltstr> name_set;
public:
NameStorage() {}
NameStorage()
{
}
NameStorage(const NameStorage &stor);
virtual ~NameStorage();
virtual
~NameStorage();
/** Query for the name. If the name has been stored, it
* returns its address, otherwise 0. */
const char *query(const char *name) const;
/** Insert the name if it has not been inserted yet, and
* return its new or old allocation. */
const char *insert(const char *name);
int num() const
{return (int)name_store.size();}
const char* get_name(int i) const
{return name_store[i];}
int
num() const
{
return (int) name_store.size();
}
const char *
get_name(int i) const
{
return name_store[i];
}
/** Debug print. */
void print() const;
};
class Constants : public AtomValues {
class Constants : public AtomValues
{
public:
/** Type for a map mapping tree indices to double values. */
typedef map<int, double> Tconstantmap;
@ -64,15 +80,21 @@ namespace ogp {
/** Map mapping a tree index of a constant to its double value. */
Tconstantmap cmap;
public:
Constants() {}
Constants()
{
}
/** Copy constructor. */
Constants(const Constants &c)
: cmap(c.cmap), cinvmap(c.cinvmap) {}
: cmap(c.cmap), cinvmap(c.cinvmap)
{
}
/** Copy constructor registering the constants in the given
* tree. The mapping from old tree indices to new ones is
* traced in tmap. */
Constants(const Constants &c, OperationTree &otree, Tintintmap &tmap)
{import_constants(c, otree, tmap);}
{
import_constants(c, otree, tmap);
}
/** Import constants registering their tree indices in the
* given tree. The mapping form old tree indices to new ones
* is traced in tmap. */
@ -96,8 +118,11 @@ namespace ogp {
int check(const char *str) const;
/** Debug print. */
void print() const;
const Tconstantmap& get_constantmap() const
{return cmap;}
const Tconstantmap &
get_constantmap() const
{
return cmap;
}
private:
/** Inverse map to Tconstantmap. */
typedef map<double, int> Tconstantinvmap;
@ -112,7 +137,8 @@ namespace ogp {
* leads. This abstraction does not distinguish between a parameter
* and a variable without lag or lead. In this sense, everything is a
* variable.*/
class DynamicAtoms : public Atoms, public Constants {
class DynamicAtoms : public Atoms, public Constants
{
public:
/** Definition of a type mapping lags to the indices of the variables. */
typedef map<int, int> Tlagmap;
@ -144,7 +170,9 @@ namespace ogp {
/** Construct empty DynamicAtoms. */
DynamicAtoms();
DynamicAtoms(const DynamicAtoms &da);
virtual ~DynamicAtoms() {}
virtual ~DynamicAtoms()
{
}
/** Check the nulary term identified by its string
* representation. The nulary term can be either a constant or
* a variable. If constant, -1 is returned so that it could be
@ -158,8 +186,11 @@ namespace ogp {
* returns -1. */
void assign(const char *name, int t);
/** Return a number of all variables. */
int nvar() const
{ return nv; }
int
nvar() const
{
return nv;
}
/** Return the vector of variable indices. */
vector<int> variables() const;
/** Return max lead and min lag for a variable given by the
@ -194,15 +225,24 @@ namespace ogp {
* exception if the tree index t is not a named atom. */
int lead(int t) const;
/** Return maximum lead. */
int get_maxlead() const
{return maxlead;}
int
get_maxlead() const
{
return maxlead;
}
/** Return minimum lag. */
int get_minlag() const
{return minlag;}
int
get_minlag() const
{
return minlag;
}
/** Return the name storage to allow querying to other
* classes. */
const NameStorage& get_name_storage() const
{return varnames;}
const NameStorage &
get_name_storage() const
{
return varnames;
}
/** Assign the variable with a given lead. The varname must be
* from the varnames storage. The method checks if the
* variable iwht the given lead/lag is not assigned. If so, an
@ -239,7 +279,6 @@ namespace ogp {
static bool is_string_constant(const char *str);
};
/** This class is a parent of all orderings of the dynamic atoms
* of variables which can appear before t, at t, or after t. It
* encapsulates the ordering, and the information about the number
@ -261,7 +300,8 @@ namespace ogp {
* preimplemented method, or plugging your own implementation. The
* method do_ordering() is called by the user after the constructor.
*/
class VarOrdering {
class VarOrdering
{
protected:
/** Number of static variables. */
int n_stat;
@ -315,39 +355,63 @@ namespace ogp {
* reimplemented. */
VarOrdering(const vector<const char *> &vnames, const DynamicAtoms &a)
: n_stat(0), n_pred(0), n_both(0), n_forw(0), varnames(vnames), atoms(a)
{}
{
}
VarOrdering(const VarOrdering &vo, const vector<const char *> &vnames,
const DynamicAtoms &a);
virtual VarOrdering *clone(const vector<const char *> &vnames,
const DynamicAtoms &a) const = 0;
/** Destructor does nothing here. */
virtual ~VarOrdering() {}
virtual ~VarOrdering()
{
}
/** This is the method setting the ordering and the map. A
* subclass must reimplement it, possibly using a
* preimplemented ordering. This method must be called by the
* user after the class has been created. */
virtual void do_ordering() = 0;
/** Return number of static. */
int nstat() const
{return n_stat;}
int
nstat() const
{
return n_stat;
}
/** Return number of predetermined. */
int npred() const
{return n_pred;}
int
npred() const
{
return n_pred;
}
/** Return number of both. */
int nboth() const
{return n_both;}
int
nboth() const
{
return n_both;
}
/** Return number of forward looking. */
int nforw() const
{return n_forw;}
int
nforw() const
{
return n_forw;
}
/** Return the set of tree indices for derivatives. */
const vector<int>& get_der_atoms() const
{return der_atoms;}
const vector<int> &
get_der_atoms() const
{
return der_atoms;
}
/** Return the y2outer. */
const vector<int>& get_y2outer() const
{return y2outer;}
const vector<int> &
get_y2outer() const
{
return y2outer;
}
/** Return the outer2y. */
const vector<int>& get_outer2y() const
{return outer2y;}
const vector<int> &
get_outer2y() const
{
return outer2y;
}
/** Query the atom given by the tree index. True is returned
* if the atom is one of the variables in the object. */
bool check(int t) const;
@ -358,8 +422,11 @@ namespace ogp {
/** This returns a length of ordered row of atoms. In all
* cases so far, it does not depend on the ordering and it is
* as follows. */
int length() const
{return n_stat+2*n_pred+3*n_both+2*n_forw;}
int
length() const
{
return n_stat+2*n_pred+3*n_both+2*n_forw;
}
/** Debug print. */
void print() const;
protected:
@ -373,16 +440,22 @@ namespace ogp {
* stat(t), pred(t), both(t), forw(t), both(t+1),
* forw(t+1). It builds the der_atoms, the map of positions,
* as well as y2outer and outer2y. */
void do_pbspbfbf()
{do_general(pbspbfbf);}
void
do_pbspbfbf()
{
do_general(pbspbfbf);
}
/** This is a preimplemented ordering for do_ordering()
* method. It assumes that the variables appear only at time
* t-1, t, t+1. It orders the atoms as both(t+1), forw(t+1),
* stat(t), pred(t), both(t), forw(t), pred(t-1),
* both(t-1). It builds the der_atoms, the map of positions,
* as well as y2outer and outer2y. */
void do_bfspbfpb()
{do_general(bfspbfpb);}
void
do_bfspbfpb()
{
do_general(bfspbfpb);
}
/** This is a preimplemented ordering for do_ordering()
* method. It makes no assumptions about occurences of
* variables at different times. It orders the atoms with

View File

@ -10,7 +10,8 @@
#include <vector>
#include <string>
namespace ogp {
namespace ogp
{
using std::vector;
using std::string;
@ -19,50 +20,83 @@ namespace ogp {
* assumes that we have only time t-1, t, and t+1, orders them as
* pred(t-1), both(t-1), stat(t), pred(t), both(t), forw(t),
* both(t+1), forw(t+1). */
class EndoVarOrdering1 : public VarOrdering {
class EndoVarOrdering1 : public VarOrdering
{
public:
EndoVarOrdering1(const vector<const char *> &vnames, const DynamicAtoms &a)
: VarOrdering(vnames, a) {}
: VarOrdering(vnames, a)
{
}
EndoVarOrdering1(const EndoVarOrdering1 &vo, const vector<const char *> &vnames,
const DynamicAtoms &a)
: VarOrdering(vo, vnames, a) {}
VarOrdering* clone(const vector<const char*>& vnames, const DynamicAtoms& a) const
{return new EndoVarOrdering1(*this, vnames, a);}
void do_ordering()
{do_pbspbfbf();}
: VarOrdering(vo, vnames, a)
{
}
VarOrdering *
clone(const vector<const char *> &vnames, const DynamicAtoms &a) const
{
return new EndoVarOrdering1(*this, vnames, a);
}
void
do_ordering()
{
do_pbspbfbf();
}
};
/** This is just another ordering used for endogenous
* variables. It assumes that we have only time t-1, t, and t+1,
* orders them as both(t+1), forw(t+1), pred(t-1), both(t-1),
* stat(t), pred(t), both(t), forw(t). */
class EndoVarOrdering2 : public VarOrdering {
class EndoVarOrdering2 : public VarOrdering
{
public:
EndoVarOrdering2(const vector<const char *> &vnames, const DynamicAtoms &a)
: VarOrdering(vnames, a) {}
: VarOrdering(vnames, a)
{
}
EndoVarOrdering2(const EndoVarOrdering2 &vo, const vector<const char *> &vnames,
const DynamicAtoms &a)
: VarOrdering(vo, vnames, a) {}
VarOrdering* clone(const vector<const char*>& vnames, const DynamicAtoms& a) const
{return new EndoVarOrdering2(*this, vnames, a);}
void do_ordering()
{do_bfspbfpb();}
: VarOrdering(vo, vnames, a)
{
}
VarOrdering *
clone(const vector<const char *> &vnames, const DynamicAtoms &a) const
{
return new EndoVarOrdering2(*this, vnames, a);
}
void
do_ordering()
{
do_bfspbfpb();
}
};
/** This is just ordering used for exogenous variables. It makes
* no assumptions about their timing. It orders them from the
* least time to the latest time. */
class ExoVarOrdering : public VarOrdering {
class ExoVarOrdering : public VarOrdering
{
public:
ExoVarOrdering(const vector<const char *> &vnames, const DynamicAtoms &a)
: VarOrdering(vnames, a) {}
: VarOrdering(vnames, a)
{
}
ExoVarOrdering(const ExoVarOrdering &vo, const vector<const char *> &vnames,
const DynamicAtoms &a)
: VarOrdering(vo, vnames, a) {}
VarOrdering* clone(const vector<const char*>& vnames, const DynamicAtoms& a) const
{return new ExoVarOrdering(*this, vnames, a);}
void do_ordering()
{do_increasing_time();}
: VarOrdering(vo, vnames, a)
{
}
VarOrdering *
clone(const vector<const char *> &vnames, const DynamicAtoms &a) const
{
return new ExoVarOrdering(*this, vnames, a);
}
void
do_ordering()
{
do_increasing_time();
}
};
class FineAtoms;
@ -71,7 +105,8 @@ namespace ogp {
* and exo). It maps the ordering to the particular outer
* orderings of endo and exo. It works tightly with the FineAtoms
* class. */
class AllvarOuterOrdering {
class AllvarOuterOrdering
{
protected:
/** Type for a map mapping a variable name to an integer. */
typedef map<const char *, int, ltstr> Tvarintmap;
@ -98,14 +133,23 @@ namespace ogp {
/** Copy constructor using the storage of provided atoms. */
AllvarOuterOrdering(const AllvarOuterOrdering &allvar_outer, const FineAtoms &a);
/** Return endo2all mapping. */
const vector<int>& get_endo2all() const
{return endo2all;}
const vector<int> &
get_endo2all() const
{
return endo2all;
}
/** Return exo2all mapping. */
const vector<int>& get_exo2all() const
{return exo2all;}
const vector<int> &
get_exo2all() const
{
return exo2all;
}
/** Return the allvar ordering. */
const vector<const char*>& get_allvar() const
{return allvar;}
const vector<const char *> &
get_allvar() const
{
return allvar;
}
};
/** This class refines the DynamicAtoms by distinguishing among
@ -127,7 +171,8 @@ namespace ogp {
* concatenation of endo and exo variables in their internal
* orderings. This is the ordering with respect to which all
* derivatives are taken. */
class FineAtoms : public DynamicAtoms {
class FineAtoms : public DynamicAtoms
{
friend class AllvarOuterOrdering;
protected:
typedef map<const char *, int, ltstr> Tvarintmap;
@ -184,14 +229,19 @@ namespace ogp {
vector<int> exo_atoms_map;
public:
FineAtoms()
: endo_order(NULL), exo_order(NULL), allvar_order(NULL) {}
: endo_order(NULL), exo_order(NULL), allvar_order(NULL)
{
}
FineAtoms(const FineAtoms &fa);
/** Deletes endo_order and exo_order. */
virtual ~FineAtoms()
{
if (endo_order) delete endo_order;
if (exo_order) delete exo_order;
if (allvar_order) delete allvar_order;
if (endo_order)
delete endo_order;
if (exo_order)
delete exo_order;
if (allvar_order)
delete allvar_order;
}
/** Overrides DynamicAtoms::check_variable so that the error
* would be raised if the variable name is not declared. A
@ -200,23 +250,38 @@ namespace ogp {
* subclass. */
int check_variable(const char *name) const;
/** This calculates min lag and max lead of endogenous variables. */
void endovarspan(int& mlead, int& mlag) const
{varspan(endovars, mlead, mlag);}
void
endovarspan(int &mlead, int &mlag) const
{
varspan(endovars, mlead, mlag);
}
/** This calculates mim lag and max lead of exogenous variables. */
void exovarspan(int& mlead, int& mlag) const
{varspan(exovars, mlead, mlag);}
void
exovarspan(int &mlead, int &mlag) const
{
varspan(exovars, mlead, mlag);
}
/** This calculates the number of periods in which at least
* one exogenous variable occurs. */
int num_exo_periods() const;
/** Return an (external) ordering of parameters. */
const vector<const char*>& get_params() const
{return params;}
const vector<const char *> &
get_params() const
{
return params;
}
/** Return an external ordering of endogenous variables. */
const vector<const char*>& get_endovars() const
{return endovars;}
const vector<const char *> &
get_endovars() const
{
return endovars;
}
/** Return an external ordering of exogenous variables. */
const vector<const char*>& get_exovars() const
{return exovars;}
const vector<const char *> &
get_exovars() const
{
return exovars;
}
/** This constructs internal orderings and makes the indices
* returned by variables method available. Further it
* constructs outer ordering of all variables by a simple
@ -302,20 +367,35 @@ namespace ogp {
int name2outer_allvar(const char *name) const;
/** Return the number of endogenous variables at time t-1, these are state
* variables. */
int nys() const
{return npred()+nboth();}
int
nys() const
{
return npred()+nboth();
}
/** Return the number of endogenous variables at time t+1. */
int nyss() const
{return nboth()+nforw();}
int
nyss() const
{
return nboth()+nforw();
}
/** Return the number of endogenous variables. */
int ny() const
{return endovars.size();}
int
ny() const
{
return endovars.size();
}
/** Return the number of exogenous variables. */
int nexo() const
{return (int)exovars.size();}
int
nexo() const
{
return (int) exovars.size();
}
/** Return the number of parameters. */
int np() const
{return (int)(params.size());}
int
np() const
{
return (int) (params.size());
}
/** Register unique endogenous variable name. The order of
* calls defines the endo outer ordering. The method is
* virtual, since a superclass may want to do some additional

View File

@ -5,15 +5,21 @@
#include "tree.h"
namespace ogp {
namespace ogp
{
using std::vector;
/** Pure virtual class defining a minimal interface for
* representation of nulary terms within FormulaParser. */
class Atoms {
class Atoms
{
public:
Atoms() {}
virtual ~Atoms() {}
Atoms()
{
}
virtual ~Atoms()
{
}
/** This returns previously assigned internal index to the
* given atom, or returns -1 if the atom has not been assigned
* yet. The method can raise an exception, if the Atoms
@ -39,16 +45,20 @@ namespace ogp {
* implementations of this class will have to be connected with
* Atoms to have knowledge about the atoms and their indices in
* the tree, and will call EvalTree::set_nulary. */
class AtomValues {
class AtomValues
{
public:
virtual ~AtomValues() {}
virtual ~AtomValues()
{
}
virtual void setValues(EvalTree &et) const = 0;
};
class FormulaDerEvaluator;
class FoldMultiIndex;
/** For ordering FoldMultiIndex in the std::map. */
struct ltfmi {
struct ltfmi
{
bool operator()(const FoldMultiIndex &i1, const FoldMultiIndex &i2) const;
};
@ -64,7 +74,8 @@ namespace ogp {
* iterators of the map do not survive the insertions to the map,
* and implementation of the constructor has to be very difficult.
*/
class FormulaDerivatives {
class FormulaDerivatives
{
friend class FormulaDerEvaluator;
protected:
/** Vector of derivatives. This is a list of derivatives (tree
@ -100,12 +111,17 @@ namespace ogp {
FormulaDerivatives(OperationTree &otree, const vector<int> &vars, int f, int max_order);
/** Copy constructor. */
FormulaDerivatives(const FormulaDerivatives &fd);
virtual ~FormulaDerivatives(){}
virtual ~FormulaDerivatives()
{
}
/** Random access to the derivatives via multiindex. */
int derivative(const FoldMultiIndex &mi) const;
/** Return the order. */
int get_order() const
{return order;}
int
get_order() const
{
return order;
}
/** Debug print. */
void print(const OperationTree &otree) const;
};
@ -122,7 +138,8 @@ namespace ogp {
* be used in constructors of FormulaEvaluator and
* FormulaDerEvaluator in order to evaluate formulas (zero
* derivatives) and higher derivatives resp. */
class FormulaParser {
class FormulaParser
{
friend class FormulaCustomEvaluator;
friend class FormulaDerEvaluator;
protected:
@ -140,10 +157,13 @@ namespace ogp {
public:
/** Construct an empty formula parser. */
FormulaParser(Atoms &a)
: atoms(a) {}
: atoms(a)
{
}
/** Copy constructor using a different instance of Atoms. */
FormulaParser(const FormulaParser &fp, Atoms &a);
virtual ~FormulaParser();
virtual
~FormulaParser();
/** Requires an addition of the formula; called from the
* parser. */
@ -162,18 +182,27 @@ namespace ogp {
/** Adds a derivative to the tree. This just calls
* OperationTree::add_derivative. */
int add_derivative(int t, int v)
{return otree.add_derivative(t, v);}
int
add_derivative(int t, int v)
{
return otree.add_derivative(t, v);
}
/** Adds a substitution. This just calls
* OperationTree::add_substitution. */
int add_substitution(int t, const map<int,int>& subst)
{return otree.add_substitution(t, subst);}
int
add_substitution(int t, const map<int, int> &subst)
{
return otree.add_substitution(t, subst);
}
/** Add the substitution given by the map where left sides of
* substitutions come from another parser. The right sides are
* from this object. The given t is from the given parser fp. */
int add_substitution(int t, const map<int,int>& subst,
int
add_substitution(int t, const map<int, int> &subst,
const FormulaParser &fp)
{return otree.add_substitution(t, subst, fp.otree);}
{
return otree.add_substitution(t, subst, fp.otree);
}
/** This adds formulas from the given parser with (possibly)
* different atoms applying substitutions from the given map
* mapping atoms from fp to atoms of the object. */
@ -187,12 +216,18 @@ namespace ogp {
* anything do with atoms, but usually some action is also
* needed (at leat to assign the tree index t to some
* atom). */
void nularify(int t)
{otree.nularify(t);}
void
nularify(int t)
{
otree.nularify(t);
}
/** Returns a set of nulary terms of the given term. Just
* calls OperationTree::nulary_of_term. */
const unordered_set<int>& nulary_of_term(int t) const
{return otree.nulary_of_term(t);}
const unordered_set<int> &
nulary_of_term(int t) const
{
return otree.nulary_of_term(t);
}
/** Parse a given string containing one or more formulas. The
* formulas are parsed and added to the OperationTree and to
@ -216,9 +251,11 @@ namespace ogp {
int last_formula() const;
/** This returns a tree index of the i-th formula in the
* vector. */
int formula(int i) const
{return formulas[i];}
int
formula(int i) const
{
return formulas[i];
}
/** This returns a tree index of the last formula and pops its
* item from the formulas vector. The number of formulas is
@ -228,19 +265,34 @@ namespace ogp {
int pop_last_formula();
/** This returns a number of formulas. */
int nformulas() const
{return (int)(formulas.size());}
int
nformulas() const
{
return (int) (formulas.size());
}
/** This returns a reference to atoms. */
const Atoms& getAtoms() const
{return atoms;}
Atoms& getAtoms()
{return atoms;}
const Atoms &
getAtoms() const
{
return atoms;
}
Atoms &
getAtoms()
{
return atoms;
}
/** This returns the tree. */
const OperationTree& getTree() const
{return otree;}
OperationTree& getTree()
{return otree;}
const OperationTree &
getTree() const
{
return otree;
}
OperationTree &
getTree()
{
return otree;
}
/** Debug print. */
void print() const;
@ -256,9 +308,12 @@ namespace ogp {
* classes which will load the results of formula (zero
* derivative) evaluations. A primitive implementation of this
* class can be a vector of doubles. */
class FormulaEvalLoader {
class FormulaEvalLoader
{
public:
virtual ~FormulaEvalLoader() {}
virtual ~FormulaEvalLoader()
{
}
/** Set the value res for the given formula. The formula is
* identified by an index corresponding to the ordering in
* which the formulas have been parsed (starting from
@ -272,7 +327,8 @@ namespace ogp {
* terms in the beginning. Using this constructor, one has to make
* sure, that the terms in the beginning do not refer to terms
* behind the initial part. */
class FormulaCustomEvaluator {
class FormulaCustomEvaluator
{
protected:
/** The evaluation tree. */
EvalTree etree;
@ -282,11 +338,13 @@ namespace ogp {
/** Construct from FormulaParser and given list of terms. */
FormulaCustomEvaluator(const FormulaParser &fp, const vector<int> &ts)
: etree(fp.otree), terms(ts)
{}
{
}
/** Construct from OperationTree and given list of terms. */
FormulaCustomEvaluator(const OperationTree &ot, const vector<int> &ts)
: etree(ot), terms(ts)
{}
{
}
/** Evaluate the terms using the given AtomValues and load the
* results using the given loader. The loader is called for
* each term in the order of the terms. */
@ -294,23 +352,30 @@ namespace ogp {
protected:
FormulaCustomEvaluator(const FormulaParser &fp)
: etree(fp.otree, fp.last_formula()), terms(fp.formulas)
{}
{
}
};
/** This class evaluates zero derivatives of the FormulaParser. */
class FormulaEvaluator : public FormulaCustomEvaluator {
class FormulaEvaluator : public FormulaCustomEvaluator
{
public:
/** Construct from FormulaParser. */
FormulaEvaluator(const FormulaParser &fp)
: FormulaCustomEvaluator(fp) {}
: FormulaCustomEvaluator(fp)
{
}
};
/** This is a pure virtual class defining an interface for all
* classes which will load the results of formula derivative
* evaluations. */
class FormulaDerEvalLoader {
class FormulaDerEvalLoader
{
public:
virtual ~FormulaDerEvalLoader() {}
virtual ~FormulaDerEvalLoader()
{
}
/** This loads the result of the derivative of the given
* order. The semantics of i is the same as in
* FormulaEvalLoader::load. The indices of variables with
@ -323,7 +388,8 @@ namespace ogp {
/** This class is a utility class representing the tensor
* multindex. It can basically increment itself, and calculate
* its offset in the folded tensor. */
class FoldMultiIndex {
class FoldMultiIndex
{
/** Number of variables. */
int nvar;
/** Dimension. */
@ -348,7 +414,9 @@ namespace ogp {
FoldMultiIndex(const FoldMultiIndex &fmi);
/** Desctructor. */
virtual ~FoldMultiIndex()
{delete [] data;}
{
delete [] data;
}
/** Assignment operator. */
const FoldMultiIndex &operator=(const FoldMultiIndex &fmi);
/** Operator < implementing lexicographic ordering within one
@ -359,23 +427,38 @@ namespace ogp {
void increment();
/** Return offset of the multiindex in the folded tensor. */
int offset() const;
const int& operator[](int i) const
{return data[i];}
const int &
operator[](int i) const
{
return data[i];
}
/** Return order of the multiindex, i.e. dimension of the
* tensor. */
int order() const
{return ord;}
int
order() const
{
return ord;
}
/** Return the number of variables. */
int nv() const
{return nvar;}
int
nv() const
{
return nvar;
}
/** Return the data. */
const int* ind() const
{return data;}
const int *
ind() const
{
return data;
}
/** Return true if the end of the tensor is reached. The
* result of a subsequent increment should be considered
* unpredictable. */
bool past_the_end() const
{return (ord == 0) || (data[0] == nvar);}
bool
past_the_end() const
{
return (ord == 0) || (data[0] == nvar);
}
/** Prints the multiindex in the brackets. */
void print() const;
private:
@ -383,7 +466,8 @@ namespace ogp {
};
/** This class evaluates derivatives of the FormulaParser. */
class FormulaDerEvaluator {
class FormulaDerEvaluator
{
/** Its own instance of EvalTree. */
EvalTree etree;
/** The indices of derivatives for each formula. This is a

View File

@ -15,16 +15,19 @@
// in EVERY action consuming material (this can be done with #define
// YY_USER_ACTION) and in bison you must use option %locations.
#ifndef OG_LOCATION_H
#define OG_LOCATION_H
namespace ogp {
namespace ogp
{
struct location_type {
struct location_type
{
int off; // offset of the token
int ll; // length ot the token
location_type() : off(0), ll(0) {}
location_type() : off(0), ll(0)
{
}
};
};

View File

@ -8,7 +8,8 @@
#include <cstdlib> // For NULL
#include <vector>
namespace ogp {
namespace ogp
{
using std::vector;
/** This class reads the given string and parses it as a
@ -21,7 +22,8 @@ namespace ogp {
* read items, the iterator provides information on row number and
* column number of the item. */
class MPIterator;
class MatrixParser {
class MatrixParser
{
friend class MPIterator;
protected:
/** Raw data as they were read. */
@ -32,16 +34,28 @@ namespace ogp {
int nc;
public:
MatrixParser()
: nc(0) {}
: nc(0)
{
}
MatrixParser(const MatrixParser &mp)
: data(mp.data), row_lengths(mp.row_lengths), nc(mp.nc) {}
virtual ~MatrixParser() {}
: data(mp.data), row_lengths(mp.row_lengths), nc(mp.nc)
{
}
virtual ~MatrixParser()
{
}
/** Return a number of read rows. */
int nrows() const
{return (int) row_lengths.size();}
int
nrows() const
{
return (int) row_lengths.size();
}
/** Return a maximum number of items in the rows. */
int ncols() const
{return nc;}
int
ncols() const
{
return nc;
}
/** Parses a given data. This initializes the object data. */
void parse(int length, const char *stream);
/** Adds newly read item. This should be called from bison
@ -66,7 +80,8 @@ namespace ogp {
/** This is an iterator intended to iterate through a matrix parsed
* by MatrixParser. The iterator provides only read-only access. */
class MPIterator {
class MPIterator
{
friend class MatrixParser;
protected:
/** Reference to the matrix parser. */
@ -79,7 +94,9 @@ namespace ogp {
int r;
public:
MPIterator() : p(NULL), i(0), c(0), r(0) {}
MPIterator() : p(NULL), i(0), c(0), r(0)
{
}
/** Constructs an iterator pointing to the beginning of the
* parsed matrix. */
MPIterator(const MatrixParser &mp);
@ -87,31 +104,47 @@ namespace ogp {
* parsed matrix. */
MPIterator(const MatrixParser &mp, const char *dummy);
/** Return read-only reference to the pointed item. */
const double& operator*() const
{return p->data[i];}
const double &
operator*() const
{
return p->data[i];
}
/** Return a row index of the pointed item. */
int row() const
{return r;}
int
row() const
{
return r;
}
/** Return a column index of the pointed item. */
int col() const
{return c;}
int
col() const
{
return c;
}
/** Assignment operator. */
const MPIterator& operator=(const MPIterator& it)
{p = it.p; i = it.i; c = it.c; r = it.r; return *this;}
const MPIterator &
operator=(const MPIterator &it)
{
p = it.p; i = it.i; c = it.c; r = it.r; return *this;
}
/** Return true if the iterators are the same, this is if they
* have the same underlying object and the same item index. */
bool operator==(const MPIterator& it) const
{return it.p == p && it.i == i;}
bool
operator==(const MPIterator &it) const
{
return it.p == p && it.i == i;
}
/** Negative of the operator==. */
bool operator!=(const MPIterator& it) const
{return ! (it == *this);}
bool
operator!=(const MPIterator &it) const
{
return !(it == *this);
}
/** Increment operator. */
MPIterator &operator++();
};
};
#endif
// Local Variables:

View File

@ -5,7 +5,8 @@
#ifndef OGP_NAMELIST
#define OGP_NAMELIST
namespace ogp {
namespace ogp
{
/** Parent class of all parsers parsing a namelist. They must
* implement add_name() method and error() method, which is called
@ -16,9 +17,12 @@ namespace ogp {
* NameListParser::namelist_parse(int lengt, const char*
* text). When implementing error(), one may consult global
* location_type namelist_lloc. */
class NameListParser {
class NameListParser
{
public:
virtual ~NameListParser() {}
virtual ~NameListParser()
{
}
virtual void add_name(const char *name) = 0;
virtual void namelist_error(const char *mes) = 0;
void namelist_parse(int length, const char *text);

View File

@ -7,7 +7,8 @@
#include <string>
namespace ogp {
namespace ogp
{
using std::string;
/** This is an easy exception, which, besides the message, stores
@ -17,7 +18,8 @@ namespace ogp {
* semantics here. They should be documented in the function which
* throws an exception and sets them. Their default value is -1,
* which means they have not been set. */
class ParserException {
class ParserException
{
protected:
char *mes;
int off;
@ -41,24 +43,49 @@ namespace ogp {
* others and forgetting the last three. */
ParserException(const ParserException &e, const char *dum, int i1, int i2, int i3);
ParserException(const ParserException &e);
virtual ~ParserException();
virtual
~ParserException();
void print(FILE *fd) const;
const char* message() const
{return mes;}
int offset() const
{return off;}
const int& i1() const
{return aux_i1;}
int& i1()
{return aux_i1;}
const int& i2() const
{return aux_i2;}
int& i2()
{return aux_i2;}
const int& i3() const
{return aux_i3;}
int& i3()
{return aux_i3;}
const char *
message() const
{
return mes;
}
int
offset() const
{
return off;
}
const int &
i1() const
{
return aux_i1;
}
int &
i1()
{
return aux_i1;
}
const int &
i2() const
{
return aux_i2;
}
int &
i2()
{
return aux_i2;
}
const int &
i3() const
{
return aux_i3;
}
int &
i3()
{
return aux_i3;
}
protected:
void copy(const ParserException &e);
};

View File

@ -7,9 +7,11 @@
#include "dynamic_atoms.h"
namespace ogp {
namespace ogp
{
class StaticAtoms : public Atoms, public Constants {
class StaticAtoms : public Atoms, public Constants
{
protected:
typedef map<const char *, int, ltstr> Tvarmap;
typedef map<int, const char *> Tinvmap;
@ -25,7 +27,8 @@ namespace ogp {
Tinvmap indices;
public:
StaticAtoms() : Atoms(), Constants(), varnames(), varorder(), vars()
{}
{
}
/* Copy constructor. */
StaticAtoms(const StaticAtoms &a);
/** Conversion from DynamicAtoms. This takes all atoms from
@ -35,9 +38,13 @@ namespace ogp {
* to new tree indices. */
StaticAtoms(const DynamicAtoms &da, OperationTree &otree, Tintintmap &tmap)
: Atoms(), Constants(), varnames(), varorder(), vars()
{import_atoms(da, otree, tmap);}
{
import_atoms(da, otree, tmap);
}
/* Destructor. */
virtual ~StaticAtoms() {}
virtual ~StaticAtoms()
{
}
/** This imports atoms from dynamic atoms inserting the new
* tree indices to the given tree (including constants). The
* mapping from old atoms to new atoms is traced in tmap. */
@ -51,8 +58,11 @@ namespace ogp {
/** This assigns a given tree index to the variable name. The
* name should have been checked before the call. */
void assign(const char *name, int t);
int nvar() const
{return varnames.num();}
int
nvar() const
{
return varnames.num();
}
/** This returns a vector of all variables. */
vector<int> variables() const;
/** This returns a tree index of the given variable. */
@ -61,8 +71,11 @@ namespace ogp {
* returned if the tree index doesn't exist. */
const char *inv_index(int t) const;
/** This returns a name in a outer ordering. (There is no other ordering.) */
const char* name(int i) const
{return varorder[i];}
const char *
name(int i) const
{
return varorder[i];
}
/** Debug print. */
void print() const;
/** This registers a variable. A subclass can reimplement
@ -72,8 +85,11 @@ namespace ogp {
virtual void register_name(const char *name);
/** Return the name storage to allow querying to other
* classes. */
const NameStorage& get_name_storage() const
{return varnames;}
const NameStorage &
get_name_storage() const
{
return varnames;
}
protected:
/** This checks the variable. The implementing subclass might
* want to throw an exception if the variable has not been

View File

@ -8,7 +8,8 @@
#include "static_atoms.h"
#include "fine_atoms.h"
namespace ogp {
namespace ogp
{
/** This class represents static atoms distinguishing between
* parameters, endogenous and exogenous variables. The class
@ -18,7 +19,8 @@ namespace ogp {
* the latter case, one can decide if the ordering of this static
* atoms should be internal or external ordering of the original
* dynamic fine atoms. */
class StaticFineAtoms : public StaticAtoms {
class StaticFineAtoms : public StaticAtoms
{
public:
typedef map<int, int> Tintintmap;
protected:
@ -59,7 +61,9 @@ namespace ogp {
* atoms of exogenous variables. */
vector<int> exo_atoms_map;
public:
StaticFineAtoms() {}
StaticFineAtoms()
{
}
/** Copy constructor making a new storage for atom names. */
StaticFineAtoms(const StaticFineAtoms &sfa);
/** Conversion from dynamic FineAtoms taking its outer
@ -68,7 +72,9 @@ namespace ogp {
* tree indices of the dynamic atoms to tree indices of the
* static atoms. */
StaticFineAtoms(const FineAtoms &fa, OperationTree &otree, Tintintmap &tmap)
{StaticFineAtoms::import_atoms(fa, otree, tmap);}
{
StaticFineAtoms::import_atoms(fa, otree, tmap);
}
/** Conversion from dynamic FineAtoms taking its internal
* ordering as ordering of parameters, endogenous and
* exogenous. A biproduct is an integer to integer map mapping
@ -76,8 +82,12 @@ namespace ogp {
* static atoms. */
StaticFineAtoms(const FineAtoms &fa, OperationTree &otree, Tintintmap &tmap,
const char *dummy)
{StaticFineAtoms::import_atoms(fa, otree, tmap, dummy);}
virtual ~StaticFineAtoms() {}
{
StaticFineAtoms::import_atoms(fa, otree, tmap, dummy);
}
virtual ~StaticFineAtoms()
{
}
/** This adds atoms from dynamic atoms inserting new tree
* indices to the given tree and tracing the mapping from old
* atoms to new atoms in tmap. The ordering of the static
@ -96,28 +106,46 @@ namespace ogp {
* methods. This a responsibility of a subclass. */
int check_variable(const char *name) const;
/** Return an (external) ordering of parameters. */
const vector<const char*>& get_params() const
{return params;}
const vector<const char *> &
get_params() const
{
return params;
}
/** Return an external ordering of endogenous variables. */
const vector<const char*>& get_endovars() const
{return endovars;}
const vector<const char *> &
get_endovars() const
{
return endovars;
}
/** Return an external ordering of exogenous variables. */
const vector<const char*>& get_exovars() const
{return exovars;}
const vector<const char *> &
get_exovars() const
{
return exovars;
}
/** This constructs der_atoms, and the endo_endoms_map and
* exo_atoms_map, which can be created only after the parsing
* is finished. */
void parsing_finished();
/** Return the atoms with respect to which we are going to
* differentiate. */
vector<int> variables() const
{return der_atoms;}
vector<int>
variables() const
{
return der_atoms;
}
/** Return the endo_atoms_map. */
const vector<int>& get_endo_atoms_map() const
{return endo_atoms_map;}
const vector<int> &
get_endo_atoms_map() const
{
return endo_atoms_map;
}
/** Return the exo_atoms_map. */
const vector<int>& get_exo_atoms_map() const
{return endo_atoms_map;}
const vector<int> &
get_exo_atoms_map() const
{
return endo_atoms_map;
}
/** Return an index in the outer ordering of a given
* parameter. An exception is thrown if the name is not a
* parameter. */
@ -131,14 +159,23 @@ namespace ogp {
* and exogenous variable. */
int name2outer_exo(const char *name) const;
/** Return the number of endogenous variables. */
int ny() const
{return endovars.size();}
int
ny() const
{
return endovars.size();
}
/** Return the number of exogenous variables. */
int nexo() const
{return (int)exovars.size();}
int
nexo() const
{
return (int) exovars.size();
}
/** Return the number of parameters. */
int np() const
{return (int)(params.size());}
int
np() const
{
return (int) (params.size());
}
/** Register unique endogenous variable name. The order of
* calls defines the endo outer ordering. The method is
* virtual, since a superclass may want to do some additional

View File

@ -10,7 +10,8 @@
#include <boost/unordered_set.hpp>
#include <cstdio>
namespace ogp {
namespace ogp
{
using boost::unordered_set;
using boost::unordered_map;
@ -27,7 +28,8 @@ namespace ogp {
ERFC, PLUS, MINUS, TIMES, DIVIDE, POWER};
/** Class representing a nulary, unary, or binary operation. */
class Operation {
class Operation
{
protected:
/** Code of the operation. */
code_t code;
@ -39,19 +41,28 @@ namespace ogp {
public:
/** Constructs a binary operation. */
Operation(code_t cd, int oper1, int oper2)
: code(cd), op1(oper1), op2(oper2) {}
: code(cd), op1(oper1), op2(oper2)
{
}
/** Constructs a unary operation. */
Operation(code_t cd, int oper1)
: code(cd), op1(oper1), op2(-1) {}
: code(cd), op1(oper1), op2(-1)
{
}
/** Constructs a nulary operation. */
Operation()
: code(NONE), op1(-1), op2(-1) {}
: code(NONE), op1(-1), op2(-1)
{
}
/** A copy constructor. */
Operation(const Operation &op)
: code(op.code), op1(op.op1), op2(op.op2) {}
: code(op.code), op1(op.op1), op2(op.op2)
{
}
/** Operator =. */
const Operation& operator=(const Operation& op)
const Operation &
operator=(const Operation &op)
{
code = op.code;
op1 = op.op1;
@ -59,55 +70,79 @@ namespace ogp {
return *this;
}
/** Operator ==. */
bool operator==(const Operation& op) const
bool
operator==(const Operation &op) const
{
return code == op.code && op1 == op.op1 && op2 == op.op2;
}
/** Operator < implementing lexicographic ordering. */
bool operator<(const Operation& op) const
bool
operator<(const Operation &op) const
{
return (code < op.code ||
code == op.code &&
(op1 < op.op1 || op1 == op.op1 && op2 < op.op2));
return (code < op.code
|| code == op.code
&& (op1 < op.op1 || op1 == op.op1 && op2 < op.op2));
}
/** Returns a number of operands. */
int nary() const
int
nary() const
{
return (op2 == -1) ? ((op1 == -1) ? 0 : 1) : 2;
}
/** Returns a hash value of the operation. */
size_t hashval() const
size_t
hashval() const
{
return op2+1 + (op1+1)^15 + code^30;
}
code_t getCode() const
{ return code; }
int getOp1() const
{ return op1; }
int getOp2() const
{ return op2; }
code_t
getCode() const
{
return code;
}
int
getOp1() const
{
return op1;
}
int
getOp2() const
{
return op2;
}
};
/** This struct is a predicate for ordering of the operations in
* OperationTree class. now obsolete */
struct ltoper {
bool operator()(const Operation& oper1, const Operation& oper2) const
{return oper1 < oper2;}
struct ltoper
{
bool
operator()(const Operation &oper1, const Operation &oper2) const
{
return oper1 < oper2;
}
};
/** Hash function object for Operation. */
struct ophash {
size_t operator()(const Operation& op) const
{ return op.hashval(); }
struct ophash
{
size_t
operator()(const Operation &op) const
{
return op.hashval();
}
};
/** This struct is a function object selecting some
* operations. The operation is given by a tree index. */
struct opselector {
struct opselector
{
virtual bool operator()(int t) const = 0;
virtual ~opselector() {}
virtual ~opselector()
{
}
};
/** Forward declaration of OperationFormatter. */
@ -143,7 +178,8 @@ namespace ogp {
* differentiate wrt more and more variables, these maps will
* become richer and richer.
*/
class OperationTree {
class OperationTree
{
friend class EvalTree;
friend class DefaultOperationFormatter;
protected:
@ -196,7 +232,8 @@ namespace ogp {
: terms(ot.terms), opmap(ot.opmap), nul_incidence(ot.nul_incidence),
derivatives(ot.derivatives),
last_nulary(ot.last_nulary)
{}
{
}
/** Add a nulary operation. The caller is responsible for not
* inserting two semantically equivalent nulary operations.
@ -256,8 +293,11 @@ namespace ogp {
void nularify(int t);
/** Return the set of nulary terms of the given term. */
const unordered_set<int>& nulary_of_term(int t) const
{return nul_incidence[t];}
const unordered_set<int> &
nulary_of_term(int t) const
{
return nul_incidence[t];
}
/** Select subterms of the given term according a given
* operation selector and return the set of terms that
@ -285,8 +325,11 @@ namespace ogp {
void forget_derivative_maps();
/** This returns an operation of a given term. */
const Operation& operation(int t) const
{return terms[t];}
const Operation &
operation(int t) const
{
return terms[t];
}
/** This outputs the operation to the given file descriptor
* using the given OperationFormatter. */
@ -296,12 +339,18 @@ namespace ogp {
void print_operation(int t) const;
/** Return the last tree index of a nulary term. */
int get_last_nulary() const
{return last_nulary;}
int
get_last_nulary() const
{
return last_nulary;
}
/** Get the number of all operations. */
int get_num_op() const
{return (int)(terms.size());}
int
get_num_op() const
{
return (int) (terms.size());
}
private:
/** This registers a calculated derivative of the term in the
* #derivatives vector.
@ -345,7 +394,8 @@ namespace ogp {
* subclasses of OperationTree and EvalTree, since we need a
* support for this in OperationTree.
*/
class EvalTree {
class EvalTree
{
protected:
/** Reference to the OperationTree over which all evaluations
* are done. */
@ -364,7 +414,9 @@ namespace ogp {
* (included). */
EvalTree(const OperationTree &otree, int last = -1);
virtual ~EvalTree()
{ delete [] values; delete [] flags; }
{
delete [] values; delete [] flags;
}
/** Set evaluation flag to all terms (besides the first
* special terms) to false. */
void reset_all();
@ -375,18 +427,24 @@ namespace ogp {
/** Debug print. */
void print() const;
/* Return the operation tree. */
const OperationTree& getOperationTree() const
{return otree;}
const OperationTree &
getOperationTree() const
{
return otree;
}
private:
EvalTree(const EvalTree &);
};
/** This is an interface describing how a given operation is
* formatted for output. */
class OperationFormatter {
class OperationFormatter
{
public:
/** Empty virtual destructor. */
virtual ~OperationFormatter() {}
virtual ~OperationFormatter()
{
}
/** Print the formatted operation op with a given tree index t
* to a given descriptor. (See class OperationTree to know
* what is a tree index.) This prints all the tree. This
@ -403,13 +461,16 @@ namespace ogp {
* reimplemented by a subclass. In addition, during its life, the
* object maintains a set of tree indices which have been output
* and they are not output any more. */
class DefaultOperationFormatter : public OperationFormatter {
class DefaultOperationFormatter : public OperationFormatter
{
protected:
const OperationTree &otree;
set<int> stop_set;
public:
DefaultOperationFormatter(const OperationTree &ot)
: otree(ot) {}
: otree(ot)
{
}
/** Format the operation with the default syntax. */
void format(const Operation &op, int t, FILE *fd);
/** This prints a string represenation of the given term, for
@ -423,24 +484,32 @@ namespace ogp {
virtual void print_delim(FILE *fd) const;
};
class NularyStringConvertor {
class NularyStringConvertor
{
public:
virtual ~NularyStringConvertor() {}
virtual ~NularyStringConvertor()
{
}
/** Return the string representation of the atom with the tree
* index t. */
virtual std::string convert(int t) const = 0;
};
/** This class converts the given term to its mathematical string representation. */
class OperationStringConvertor {
class OperationStringConvertor
{
protected:
const NularyStringConvertor &nulsc;
const OperationTree &otree;
public:
OperationStringConvertor(const NularyStringConvertor &nsc, const OperationTree &ot)
: nulsc(nsc), otree(ot) {}
: nulsc(nsc), otree(ot)
{
}
/** Empty virtual destructor. */
virtual ~OperationStringConvertor() {}
virtual ~OperationStringConvertor()
{
}
/** Convert the operation to the string mathematical
* representation. This does not write any equation, just
* returns a string representation of the formula. */

View File

@ -18,45 +18,67 @@
class Dynare;
class DynareNameList : public NameList {
class DynareNameList : public NameList
{
vector<const char *> names;
public:
DynareNameList(const Dynare &dynare);
int getNum() const
{return (int)names.size();}
const char* getName(int i) const
{return names[i];}
int
getNum() const
{
return (int) names.size();
}
const char *
getName(int i) const
{
return names[i];
}
/** This for each string of the input vector calculates its index
* in the names. And returns the resulting vector of indices. If
* the name cannot be found, then an exception is raised. */
vector<int> selectIndices(const vector<const char *> &ns) const;
};
class DynareExogNameList : public NameList {
class DynareExogNameList : public NameList
{
vector<const char *> names;
public:
DynareExogNameList(const Dynare &dynare);
int getNum() const
{return (int)names.size();}
const char* getName(int i) const
{return names[i];}
int
getNum() const
{
return (int) names.size();
}
const char *
getName(int i) const
{
return names[i];
}
};
class DynareStateNameList : public NameList {
class DynareStateNameList : public NameList
{
vector<const char *> names;
public:
DynareStateNameList(const Dynare &dynare, const DynareNameList &dnl,
const DynareExogNameList &denl);
int getNum() const
{return (int)names.size();}
const char* getName(int i) const
{return names[i];}
int
getNum() const
{
return (int) names.size();
}
const char *
getName(int i) const
{
return names[i];
}
};
// The following only implements DynamicModel with help of ogdyn::DynareModel
class DynareJacobian;
class Dynare : public DynamicModel {
class Dynare : public DynamicModel
{
friend class DynareNameList;
friend class DynareExogNameList;
friend class DynareStateNameList;
@ -83,59 +105,129 @@ public:
double sstol, Journal &jr);
/** Makes a deep copy of the object. */
Dynare(const Dynare &dyn);
DynamicModel* clone() const
{return new Dynare(*this);}
virtual ~Dynare();
int nstat() const
{return model->getAtoms().nstat();}
int nboth() const
{return model->getAtoms().nboth();}
int npred() const
{return model->getAtoms().npred();}
int nforw() const
{return model->getAtoms().nforw();}
int nexog() const
{return model->getAtoms().nexo();}
int nys() const
{return model->getAtoms().nys();}
int nyss() const
{return model->getAtoms().nyss();}
int ny() const
{return model->getAtoms().ny();}
int order() const
{return model->getOrder();}
DynamicModel *
clone() const
{
return new Dynare(*this);
}
virtual
~Dynare();
int
nstat() const
{
return model->getAtoms().nstat();
}
int
nboth() const
{
return model->getAtoms().nboth();
}
int
npred() const
{
return model->getAtoms().npred();
}
int
nforw() const
{
return model->getAtoms().nforw();
}
int
nexog() const
{
return model->getAtoms().nexo();
}
int
nys() const
{
return model->getAtoms().nys();
}
int
nyss() const
{
return model->getAtoms().nyss();
}
int
ny() const
{
return model->getAtoms().ny();
}
int
order() const
{
return model->getOrder();
}
const NameList& getAllEndoNames() const
{return *dnl;}
const NameList& getStateNames() const
{return *dsnl;}
const NameList& getExogNames() const
{return *denl;}
const NameList &
getAllEndoNames() const
{
return *dnl;
}
const NameList &
getStateNames() const
{
return *dsnl;
}
const NameList &
getExogNames() const
{
return *denl;
}
TwoDMatrix& getVcov()
{return model->getVcov();}
const TwoDMatrix& getVcov() const
{return model->getVcov();}
Vector& getParams()
{return model->getParams();}
const Vector& getParams() const
{return model->getParams();}
void setInitOuter(const Vector& x)
{model->setInitOuter(x);}
TwoDMatrix &
getVcov()
{
return model->getVcov();
}
const TwoDMatrix &
getVcov() const
{
return model->getVcov();
}
Vector &
getParams()
{
return model->getParams();
}
const Vector &
getParams() const
{
return model->getParams();
}
void
setInitOuter(const Vector &x)
{
model->setInitOuter(x);
}
const TensorContainer<FSSparseTensor>& getModelDerivatives() const
{return md;}
const Vector& getSteady() const
{return *ysteady;}
Vector& getSteady()
{return *ysteady;}
const ogdyn::DynareModel& getModel() const
{return *model;}
const TensorContainer<FSSparseTensor> &
getModelDerivatives() const
{
return md;
}
const Vector &
getSteady() const
{
return *ysteady;
}
Vector &
getSteady()
{
return *ysteady;
}
const ogdyn::DynareModel &
getModel() const
{
return *model;
}
// here is true public interface
void solveDeterministicSteady(Vector &steady);
void solveDeterministicSteady()
{solveDeterministicSteady(*ysteady);}
void
solveDeterministicSteady()
{
solveDeterministicSteady(*ysteady);
}
void evaluateSystem(Vector &out, const Vector &yy, const Vector &xx);
void evaluateSystem(Vector &out, const Vector &yym, const Vector &yy,
const Vector &yyp, const Vector &xx);
@ -148,14 +240,19 @@ private:
void writeModelInfo(Journal &jr) const;
};
class DynareEvalLoader : public ogp::FormulaEvalLoader, public Vector {
class DynareEvalLoader : public ogp::FormulaEvalLoader, public Vector
{
public:
DynareEvalLoader(const ogp::FineAtoms &a, Vector &out);
void load(int i, double res)
{operator[](i) = res;}
void
load(int i, double res)
{
operator[](i) = res;
}
};
class DynareDerEvalLoader : public ogp::FormulaDerEvalLoader {
class DynareDerEvalLoader : public ogp::FormulaDerEvalLoader
{
protected:
const ogp::FineAtoms &atoms;
TensorContainer<FSSparseTensor> &md;
@ -165,27 +262,41 @@ public:
void load(int i, int iord, const int *vars, double res);
};
class DynareJacobian : public ogu::Jacobian, public ogp::FormulaDerEvalLoader {
class DynareJacobian : public ogu::Jacobian, public ogp::FormulaDerEvalLoader
{
protected:
Dynare &d;
public:
DynareJacobian(Dynare &dyn);
virtual ~DynareJacobian() {}
virtual ~DynareJacobian()
{
}
void load(int i, int iord, const int *vars, double res);
void eval(const Vector &in);
};
class DynareVectorFunction : public ogu::VectorFunction {
class DynareVectorFunction : public ogu::VectorFunction
{
protected:
Dynare &d;
public:
DynareVectorFunction(Dynare &dyn)
: d(dyn) {}
virtual ~DynareVectorFunction() {}
int inDim() const
{return d.ny();}
int outDim() const
{return d.ny();}
: d(dyn)
{
}
virtual ~DynareVectorFunction()
{
}
int
inDim() const
{
return d.ny();
}
int
outDim() const
{
return d.ny();
}
void eval(const ConstVector &in, Vector &out);
};

View File

@ -15,7 +15,8 @@
#include <map>
#include <vector>
namespace ogdyn {
namespace ogdyn
{
using std::map;
using std::vector;
@ -25,13 +26,20 @@ namespace ogdyn {
* expressions represented by tree indices. */
typedef map<const char *, int, ogp::ltstr> Tsubstmap;
class DynareStaticAtoms : public ogp::StaticAtoms {
class DynareStaticAtoms : public ogp::StaticAtoms
{
public:
DynareStaticAtoms()
: StaticAtoms() {}
: StaticAtoms()
{
}
DynareStaticAtoms(const DynareStaticAtoms &a)
: StaticAtoms(a) {}
virtual ~DynareStaticAtoms() {}
: StaticAtoms(a)
{
}
virtual ~DynareStaticAtoms()
{
}
/** This registers a unique varname identifier. It throws an
* exception if the variable name is duplicate. It checks the
* uniqueness and then it calls StaticAtoms::register_name. */
@ -43,8 +51,8 @@ namespace ogdyn {
int check_variable(const char *name) const;
};
class DynareDynamicAtoms : public ogp::SAtoms, public ogp::NularyStringConvertor {
class DynareDynamicAtoms : public ogp::SAtoms, public ogp::NularyStringConvertor
{
public:
enum atype {endovar, exovar, param};
protected:
@ -53,9 +61,13 @@ namespace ogdyn {
Tatypemap atom_type;
public:
DynareDynamicAtoms()
: ogp::SAtoms() {}
: ogp::SAtoms()
{
}
DynareDynamicAtoms(const DynareDynamicAtoms &dda);
virtual ~DynareDynamicAtoms() {}
virtual ~DynareDynamicAtoms()
{
}
/** This parses a variable of the forms: varname(+3),
* varname(3), varname, varname(-3), varname(0), varname(+0),
* varname(-0). */
@ -74,11 +86,11 @@ namespace ogdyn {
std::string convert(int t) const;
};
/** This class represents the atom values for dynare, where
* exogenous variables can occur only at time t, and endogenous at
* times t-1, t, and t+1. */
class DynareAtomValues : public ogp::AtomValues {
class DynareAtomValues : public ogp::AtomValues
{
protected:
/** Reference to the atoms (we suppose that they are only at
* t-1,t,t+1. */
@ -99,10 +111,14 @@ namespace ogdyn {
public:
DynareAtomValues(const ogp::FineAtoms &a, const Vector &pvals, const Vector &ym,
const Vector &y, const Vector &yp, const Vector &x)
: atoms(a), paramvals(pvals), yym(ym), yy(y), yyp(yp), xx(x) {}
: atoms(a), paramvals(pvals), yym(ym), yy(y), yyp(yp), xx(x)
{
}
DynareAtomValues(const ogp::FineAtoms &a, const Vector &pvals, const ConstVector &ym,
const Vector &y, const ConstVector &yp, const Vector &x)
: atoms(a), paramvals(pvals), yym(ym), yy(y), yyp(yp), xx(x) {}
: atoms(a), paramvals(pvals), yym(ym), yy(y), yyp(yp), xx(x)
{
}
void setValues(ogp::EvalTree &et) const;
};
@ -110,7 +126,8 @@ namespace ogdyn {
* makes only appropriate subvector yym and yyp of the y vector,
* makes a vector of zero exogenous variables and uses
* DynareAtomValues with more general interface. */
class DynareSteadyAtomValues : public ogp::AtomValues {
class DynareSteadyAtomValues : public ogp::AtomValues
{
protected:
/** Subvector of yy. */
const ConstVector yym;
@ -126,12 +143,18 @@ namespace ogdyn {
yyp(y, a.nstat()+a.npred(), a.nyss()),
xx(a.nexo()),
av(a, pvals, yym, y, yyp, xx)
{xx.zeros();}
void setValues(ogp::EvalTree& et) const
{av.setValues(et);}
{
xx.zeros();
}
void
setValues(ogp::EvalTree &et) const
{
av.setValues(et);
}
};
class DynareStaticSteadyAtomValues : public ogp::AtomValues {
class DynareStaticSteadyAtomValues : public ogp::AtomValues
{
protected:
/** Reference to static atoms over which the tree, where the
* values go, is defined. */
@ -153,7 +176,9 @@ namespace ogdyn {
: atoms_static(sa),
atoms(a),
yy(yyy),
paramvals(pvals) {}
paramvals(pvals)
{
}
/** Set the values to the tree defined over the static atoms. */
void setValues(ogp::EvalTree &et) const;
};
@ -166,7 +191,8 @@ namespace ogdyn {
* them to calculated values. If a variable is already set, it
* does not override its value. It has no methods, everything is
* done in the constructor. */
class DynareSteadySubstitutions : public ogp::FormulaEvalLoader {
class DynareSteadySubstitutions : public ogp::FormulaEvalLoader
{
protected:
const ogp::FineAtoms &atoms;
public:
@ -185,7 +211,8 @@ namespace ogdyn {
* over the static tree. It also needs dynamic version of the
* atoms, since it defines ordering of the vectors pvals, and
* yy. */
class DynareStaticSteadySubstitutions : public ogp::FormulaEvalLoader {
class DynareStaticSteadySubstitutions : public ogp::FormulaEvalLoader
{
protected:
const ogp::FineAtoms &atoms;
const ogp::StaticFineAtoms &atoms_static;
@ -204,7 +231,6 @@ namespace ogdyn {
};
#endif
// Local Variables:

View File

@ -7,7 +7,8 @@
#include <string>
class DynareException {
class DynareException
{
char *mes;
public:
DynareException(const char *m, const char *fname, int line, int col)
@ -27,11 +28,18 @@ public:
}
DynareException(const DynareException &e)
: mes(new char[strlen(e.mes)+1])
{strcpy(mes, e.mes);}
{
strcpy(mes, e.mes);
}
virtual ~DynareException()
{delete [] mes;}
const char* message() const
{return mes;}
{
delete [] mes;
}
const char *
message() const
{
return mes;
}
};
#endif

View File

@ -15,7 +15,8 @@
#include <map>
#include <boost/unordered_set.hpp>
namespace ogdyn {
namespace ogdyn
{
using boost::unordered_set;
using std::map;
@ -23,33 +24,44 @@ namespace ogdyn {
* positions (including the first, excluding the second). A
* position is given by the line and the column within the line
* (both starting from 1). */
struct PosInterval {
struct PosInterval
{
int fl;
int fc;
int ll;
int lc;
PosInterval() {}
PosInterval()
{
}
PosInterval(int ifl, int ifc, int ill, int ilc)
: fl(ifl), fc(ifc), ll(ill), lc(ilc) {}
const PosInterval& operator=(const PosInterval& pi)
{fl = pi.fl; fc = pi.fc; ll = pi.ll; lc = pi.lc; return *this;}
: fl(ifl), fc(ifc), ll(ill), lc(ilc)
{
}
const PosInterval &
operator=(const PosInterval &pi)
{
fl = pi.fl; fc = pi.fc; ll = pi.ll; lc = pi.lc; return *this;
}
/** This returns the interval beginning and interval length
* within the given string. */
void translate(const char *beg, int len, const char * &ibeg, int &ilen) const;
/** Debug print. */
void print() const
{printf("fl=%d fc=%d ll=%d lc=%d\n",fl,fc,ll,lc);}
void
print() const
{
printf("fl=%d fc=%d ll=%d lc=%d\n", fl, fc, ll, lc);
}
};
/** This class is basically a GeneralMatrix but is created from
* parsed matrix data. */
class ParsedMatrix : public TwoDMatrix {
class ParsedMatrix : public TwoDMatrix
{
public:
/** Construct the object from the parsed data of ogp::MatrixParser. */
ParsedMatrix(const ogp::MatrixParser &mp);
};
class PlannerBuilder;
class PlannerInfo;
class ForwSubstBuilder;
@ -59,7 +71,8 @@ namespace ogdyn {
/** A subclass is responsible for creating param_vals, init_vals,
* and vcov_mat. */
class DynareModel {
class DynareModel
{
friend class PlannerBuilder;
friend class ForwSubstBuilder;
friend class MultInitSS;
@ -106,29 +119,57 @@ namespace ogdyn {
DynareModel();
/** Construct a new deep copy. */
DynareModel(const DynareModel &dm);
virtual ~DynareModel();
virtual
~DynareModel();
virtual DynareModel *clone() const = 0;
const DynareDynamicAtoms& getAtoms() const
{return atoms;}
const ogp::FormulaParser& getParser() const
{return eqs;}
int getOrder() const
{return order;}
const DynareDynamicAtoms &
getAtoms() const
{
return atoms;
}
const ogp::FormulaParser &
getParser() const
{
return eqs;
}
int
getOrder() const
{
return order;
}
/** Return the vector of parameter values. */
const Vector& getParams() const
{return *param_vals;}
Vector& getParams()
{return *param_vals;}
const Vector &
getParams() const
{
return *param_vals;
}
Vector &
getParams()
{
return *param_vals;
}
/** Return the vector of initial values of endo variables. */
const Vector& getInit() const
{return *init_vals;}
Vector& getInit()
{return *init_vals;}
const Vector &
getInit() const
{
return *init_vals;
}
Vector &
getInit()
{
return *init_vals;
}
/** Return the vcov matrix. */
const TwoDMatrix& getVcov() const
{return *vcov_mat;}
TwoDMatrix& getVcov()
{return *vcov_mat;}
const TwoDMatrix &
getVcov() const
{
return *vcov_mat;
}
TwoDMatrix &
getVcov()
{
return *vcov_mat;
}
/** Return planner info. */
const PlannerInfo *get_planner_info() const;
/** Return forward substitutions info. */
@ -199,7 +240,8 @@ namespace ogdyn {
* parses variable declarations, model equations, parameter
* assignments, initval assignments, vcov matrix and order of
* approximation. */
class DynareParser : public DynareModel {
class DynareParser : public DynareModel
{
protected:
/** Static atoms for parameter assignments. */
DynareStaticAtoms pa_atoms;
@ -218,41 +260,66 @@ namespace ogdyn {
* in the model file. */
DynareParser(const char *str, int len, int ord);
DynareParser(const DynareParser &p);
virtual ~DynareParser();
DynareModel* clone() const
{return new DynareParser(*this);}
virtual
~DynareParser();
DynareModel *
clone() const
{
return new DynareParser(*this);
}
/** Adds a name of endogenous, exogenous or a parameter. This
* addss the name to the parent class DynareModel and also
* registers the name to either paramset, or initval. */
void add_name(const char *name, int flag);
/** Sets position of the model section. Called from
* dynglob.y. */
void set_model_pos(int off1, int off2)
{model_beg = off1; model_end = off2;}
void
set_model_pos(int off1, int off2)
{
model_beg = off1; model_end = off2;
}
/** Sets position of the section setting parameters. Called
* from dynglob.y. */
void set_paramset_pos(int off1, int off2)
{paramset_beg = off1; paramset_end = off2;}
void
set_paramset_pos(int off1, int off2)
{
paramset_beg = off1; paramset_end = off2;
}
/** Sets position of the initval section. Called from
* dynglob.y. */
void set_initval_pos(int off1, int off2)
{initval_beg = off1; initval_end = off2;}
void
set_initval_pos(int off1, int off2)
{
initval_beg = off1; initval_end = off2;
}
/** Sets position of the vcov section. Called from
* dynglob.y. */
void set_vcov_pos(int off1, int off2)
{vcov_beg = off1; vcov_end = off2;}
void
set_vcov_pos(int off1, int off2)
{
vcov_beg = off1; vcov_end = off2;
}
/** Parser the given string as integer and set to as the
* order. */
void set_order_pos(int off1, int off2)
{order_beg = off1; order_end = off2;}
void
set_order_pos(int off1, int off2)
{
order_beg = off1; order_end = off2;
}
/** Sets position of the planner_objective section. Called
* from dynglob.y. */
void set_pl_objective_pos(int off1, int off2)
{plobjective_beg = off1; plobjective_end = off2;}
void
set_pl_objective_pos(int off1, int off2)
{
plobjective_beg = off1; plobjective_end = off2;
}
/** Sets position of the planner_discount section. Called from
* dynglob.y. */
void set_pl_discount_pos(int off1, int off2)
{pldiscount_beg = off1; pldiscount_end = off2;}
void
set_pl_discount_pos(int off1, int off2)
{
pldiscount_beg = off1; pldiscount_end = off2;
}
/** Processes a syntax error from bison. */
void error(const char *mes);
/** Debug print. */
@ -289,28 +356,39 @@ namespace ogdyn {
* no automatic substitutions cannot be done here, which in turn
* implies that we cannot do here a social planner nor substitutions
* of multiple lags. */
class DynareSPModel : public DynareModel {
class DynareSPModel : public DynareModel
{
public:
DynareSPModel(const char **endo, int num_endo,
const char **exo, int num_exo,
const char **par, int num_par,
const char *equations, int len, int ord);
DynareSPModel(const DynareSPModel &dm)
: DynareModel(dm) {}
~DynareSPModel() {}
virtual DynareModel* clone() const
{return new DynareSPModel(*this);}
: DynareModel(dm)
{
}
~DynareSPModel()
{
}
virtual DynareModel *
clone() const
{
return new DynareSPModel(*this);
}
};
/** This class implements a selector of operations which correspond
* to non-linear functions. This inherits from ogp::opselector and
* is used to calculate non-linear subterms in
* DynareModel::get_nonlinear_subterms(). */
class NLSelector : public ogp::opselector {
class NLSelector : public ogp::opselector
{
private:
const DynareModel &model;
public:
NLSelector(const DynareModel& m) : model(m) {}
NLSelector(const DynareModel &m) : model(m)
{
}
bool operator()(int t) const;
};
@ -318,13 +396,16 @@ namespace ogdyn {
* equations and the first derivatives at zero shocks and at the
* given (static) state. Static means that lags and leads are
* ignored. */
class ModelSSWriter : public ogp::DefaultOperationFormatter {
class ModelSSWriter : public ogp::DefaultOperationFormatter
{
protected:
const DynareModel &model;
public:
ModelSSWriter(const DynareModel &m)
: DefaultOperationFormatter(m.eqs.getTree()),
model(m) {}
model(m)
{
}
/** This writes the evaluation of the system. It calls pure
* virtual methods for writing a preamble, then assignment of
* atoms, and then assignment for resulting object. These are
@ -343,15 +424,17 @@ namespace ogdyn {
virtual void write_der1_assignment(FILE *fd) const = 0;
};
class MatlabSSWriter : public ModelSSWriter {
class MatlabSSWriter : public ModelSSWriter
{
protected:
/** Identifier used in function names. */
char *id;
public:
MatlabSSWriter(const DynareModel &dm, const char *idd);
virtual ~MatlabSSWriter()
{delete [] id;}
{
delete [] id;
}
protected:
// from ModelSSWriter
void write_der0_preamble(FILE *fd) const;
@ -377,13 +460,16 @@ namespace ogdyn {
/** This class implements OperationFormatter for debugging
* purposes. It renders atoms in a more friendly way than the
* ogp::DefaulOperationFormatter. */
class DebugOperationFormatter : public ogp::DefaultOperationFormatter {
class DebugOperationFormatter : public ogp::DefaultOperationFormatter
{
protected:
const DynareModel &model;
public:
DebugOperationFormatter(const DynareModel &m)
: DefaultOperationFormatter(m.getParser().getTree()),
model(m) {}
model(m)
{
}
void format_nulary(int t, FILE *fd) const;
};
};

View File

@ -13,7 +13,8 @@ simul: m max_evals (10*m)
#include <vector>
#include <string>
struct DynareParams {
struct DynareParams
{
const char *modname;
std::string basename;
int num_per;
@ -46,16 +47,31 @@ struct DynareParams {
bool version;
DynareParams(int argc, char **argv);
void printHelp() const;
int getCheckShockPoints() const
{return check_num;}
double getCheckShockScale() const
{return check_scale;}
int getCheckEllipsePoints() const
{return 10*check_num;}
double getCheckEllipseScale() const
{return 0.5*check_scale;}
int getCheckPathPoints() const
{return 10*check_num;}
int
getCheckShockPoints() const
{
return check_num;
}
double
getCheckShockScale() const
{
return check_scale;
}
int
getCheckEllipsePoints() const
{
return 10*check_num;
}
double
getCheckEllipseScale() const
{
return 0.5*check_scale;
}
int
getCheckPathPoints() const
{
return 10*check_num;
}
private:
enum {opt_per, opt_burn, opt_sim, opt_rtper, opt_rtsim, opt_condper, opt_condsim,
opt_prefix, opt_threads,

View File

@ -5,14 +5,15 @@
#ifndef FORW_SUBST_BUILDER_H
#define FORW_SUBST_BUILDER_H
#include "dynare_model.h"
namespace ogdyn {
namespace ogdyn
{
/** This struct encapsulates information about the process of
* forward substitutions. */
struct ForwSubstInfo {
struct ForwSubstInfo
{
int num_affected_equations;
int num_subst_terms;
int num_aux_variables;
@ -21,10 +22,13 @@ namespace ogdyn {
: num_affected_equations(0),
num_subst_terms(0),
num_aux_variables(0),
num_new_terms(0) {}
num_new_terms(0)
{
}
};
class ForwSubstBuilder {
class ForwSubstBuilder
{
typedef map<int, const char *> Ttermauxmap;
protected:
/** Reference to the model, to which we will add equations and
@ -49,11 +53,17 @@ namespace ogdyn {
/** Copy constructor with a new instance of the model. */
ForwSubstBuilder(const ForwSubstBuilder &b, DynareModel &m);
/** Return the auxiliary variable mapping. */
const Tsubstmap& get_aux_map() const
{return aux_map;}
const Tsubstmap &
get_aux_map() const
{
return aux_map;
}
/** Return the information. */
const ForwSubstInfo& get_info() const
{return info;}
const ForwSubstInfo &
get_info() const
{
return info;
}
private:
ForwSubstBuilder(const ForwSubstBuilder &b);
/** This method takes a nonlinear term t, and if it has leads

View File

@ -8,15 +8,20 @@
#include "twod_matrix.h"
#include "journal.h"
namespace ogu {
namespace ogu
{
class OneDFunction {
class OneDFunction
{
public:
virtual ~OneDFunction() {}
virtual ~OneDFunction()
{
}
virtual double eval(double) = 0;
};
class GoldenSectionSearch {
class GoldenSectionSearch
{
protected:
static double tol;
static double golden;
@ -38,10 +43,15 @@ namespace ogu {
static bool search_for_finite(OneDFunction &f, double x1, double &x2, double &b);
};
class VectorFunction {
class VectorFunction
{
public:
VectorFunction() {}
virtual ~VectorFunction() {}
VectorFunction()
{
}
virtual ~VectorFunction()
{
}
virtual int inDim() const = 0;
virtual int outDim() const = 0;
/** Check dimensions of eval parameters. */
@ -50,15 +60,21 @@ namespace ogu {
virtual void eval(const ConstVector &in, Vector &out) = 0;
};
class Jacobian : public TwoDMatrix {
class Jacobian : public TwoDMatrix
{
public:
Jacobian(int n)
: TwoDMatrix(n,n) {}
virtual ~Jacobian() {}
: TwoDMatrix(n, n)
{
}
virtual ~Jacobian()
{
}
virtual void eval(const Vector &in) = 0;
};
class NLSolver : public OneDFunction {
class NLSolver : public OneDFunction
{
protected:
Journal &journal;
VectorFunction &func;
@ -73,8 +89,12 @@ namespace ogu {
NLSolver(VectorFunction &f, Jacobian &j, int maxit, double tl, Journal &jr)
: journal(jr), func(f), jacob(j), max_iter(maxit), tol(tl),
xnewton(f.inDim()), xcauchy(f.inDim()), x(f.inDim())
{xnewton.zeros(); xcauchy.zeros(); x.zeros();}
virtual ~NLSolver() {}
{
xnewton.zeros(); xcauchy.zeros(); x.zeros();
}
virtual ~NLSolver()
{
}
/** Returns true if the problem has converged. xx as input is the
* starting value, as output it is a solution. */
bool solve(Vector &xx, int &iter);

View File

@ -5,7 +5,8 @@
#include "dynare_model.h"
namespace ogdyn {
namespace ogdyn
{
using boost::unordered_set;
using std::map;
@ -13,7 +14,8 @@ namespace ogdyn {
/** This is a two dimensional array of integers. Nothing
* difficult. */
class IntegerMatrix {
class IntegerMatrix
{
protected:
/** Number of rows. */
int nr;
@ -24,28 +26,47 @@ namespace ogdyn {
public:
/** Construct uninitialized array. */
IntegerMatrix(int nrr, int ncc)
: nr(nrr), nc(ncc), data(new int[nr*nc]) {}
: nr(nrr), nc(ncc), data(new int[nr*nc])
{
}
/** Copy constructor. */
IntegerMatrix(const IntegerMatrix &im)
: nr(im.nr), nc(im.nc), data(new int[nr*nc])
{memcpy(data, im.data, nr*nc*sizeof(int));}
{
memcpy(data, im.data, nr*nc*sizeof(int));
}
virtual ~IntegerMatrix()
{delete [] data;}
{
delete [] data;
}
/** Assignment operator. It can only assing array with the
* same dimensions. */
const IntegerMatrix &operator=(const IntegerMatrix &im);
int& operator()(int i, int j)
{return data[i+j*nr];}
const int& operator()(int i, int j) const
{return data[i+j*nr];}
int nrows() const
{return nr;}
int ncols() const
{return nc;}
int &
operator()(int i, int j)
{
return data[i+j*nr];
}
const int &
operator()(int i, int j) const
{
return data[i+j*nr];
}
int
nrows() const
{
return nr;
}
int
ncols() const
{
return nc;
}
};
/** The three dimensional array of integers. Nothing difficult. */
class IntegerArray3 {
class IntegerArray3
{
protected:
/** First dimension. */
int n1;
@ -58,37 +79,61 @@ namespace ogdyn {
public:
/** Constrcut unitialized array. */
IntegerArray3(int nn1, int nn2, int nn3)
: n1(nn1), n2(nn2), n3(nn3), data(new int[n1*n2*n3]) {}
: n1(nn1), n2(nn2), n3(nn3), data(new int[n1*n2*n3])
{
}
/** Copy constructor. */
IntegerArray3(const IntegerArray3 &ia3)
: n1(ia3.n1), n2(ia3.n2), n3(ia3.n3), data(new int[n1*n2*n3])
{memcpy(data, ia3.data, n1*n2*n3*sizeof(int));}
{
memcpy(data, ia3.data, n1*n2*n3*sizeof(int));
}
virtual ~IntegerArray3()
{delete [] data;}
{
delete [] data;
}
/** Assignment operator assigning the arrays with the same dimensions. */
const IntegerArray3 &operator=(const IntegerArray3 &ia3);
int& operator()(int i, int j, int k)
{return data[i+j*n1+k*n1*n2];}
const int& operator()(int i, int j, int k) const
{return data[i+j*n1+k*n1*n2];}
int dim1() const
{return n1;}
int dim2() const
{return n2;}
int dim3() const
{return n3;}
int &
operator()(int i, int j, int k)
{
return data[i+j*n1+k*n1*n2];
}
const int &
operator()(int i, int j, int k) const
{
return data[i+j*n1+k*n1*n2];
}
int
dim1() const
{
return n1;
}
int
dim2() const
{
return n2;
}
int
dim3() const
{
return n3;
}
};
/** This struct encapsulates information about the building of a
* planner's problem. */
struct PlannerInfo {
struct PlannerInfo
{
int num_lagrange_mults;
int num_aux_variables;
int num_new_terms;
PlannerInfo()
: num_lagrange_mults(0),
num_aux_variables(0),
num_new_terms(0) {}
num_new_terms(0)
{
}
};
class MultInitSS;
@ -102,7 +147,8 @@ namespace ogdyn {
* need to create static atoms and static versions of all the tree
* index matrices. The algorithm and algebra are documented in
* dynare++-ramsey.pdf. */
class PlannerBuilder {
class PlannerBuilder
{
friend class MultInitSS;
public:
/** Type for a set of variable names. */
@ -191,8 +237,11 @@ namespace ogdyn {
* is supposed to be the copy of the model in the builder. */
PlannerBuilder(const PlannerBuilder &pb, ogdyn::DynareModel &m);
/** Return the information. */
const PlannerInfo& get_info() const
{return info;}
const PlannerInfo &
get_info() const
{
return info;
}
protected:
/** Differentiate the planner objective wrt endogenous
* variables with different lags. */
@ -240,7 +289,8 @@ namespace ogdyn {
*
* The code can be run only after the parsing has been finished in
* atoms. */
class MultInitSS : public ogp::FormulaEvalLoader {
class MultInitSS : public ogp::FormulaEvalLoader
{
protected:
/** The constant reference to the builder. */
const PlannerBuilder &builder;
@ -271,7 +321,6 @@ namespace ogdyn {
};
};
#endif
// Local Variables:

View File

@ -7,8 +7,8 @@
#include "QuasiTriangular.h"
class BlockDiagonal : public QuasiTriangular {
class BlockDiagonal : public QuasiTriangular
{
int *const row_len;
int *const col_len;
public:
@ -16,10 +16,16 @@ public:
BlockDiagonal(int p, const BlockDiagonal &b);
BlockDiagonal(const BlockDiagonal &b);
BlockDiagonal(const QuasiTriangular &t);
const BlockDiagonal& operator=(const QuasiTriangular& t)
{GeneralMatrix::operator=(t); return *this;}
const BlockDiagonal &
operator=(const QuasiTriangular &t)
{
GeneralMatrix::operator=(t); return *this;
}
const BlockDiagonal &operator=(const BlockDiagonal &b);
~BlockDiagonal() {delete [] row_len; delete [] col_len;}
~BlockDiagonal()
{
delete [] row_len; delete [] col_len;
}
void setZeroBlockEdge(diag_iter edge);
int getNumZeros() const;
int getNumBlocks() const;
@ -33,8 +39,11 @@ public:
col_iter col_begin(const DiagonalBlock &b);
const_row_iter row_end(const DiagonalBlock &b) const;
row_iter row_end(const DiagonalBlock &b);
QuasiTriangular* clone() const
{return new BlockDiagonal(*this);}
QuasiTriangular *
clone() const
{
return new BlockDiagonal(*this);
}
private:
void setZerosToRU(diag_iter edge);
const_diag_iter findBlockStart(const_diag_iter from) const;
@ -47,7 +56,6 @@ private:
#endif /* BLOCK_DIAGONAL_H */
// Local Variables:
// mode:C++
// End:

View File

@ -11,7 +11,8 @@
class GeneralMatrix;
class ConstGeneralMatrix {
class ConstGeneralMatrix
{
friend class GeneralMatrix;
protected:
ConstVector data;
@ -20,19 +21,46 @@ protected:
int ld;
public:
ConstGeneralMatrix(const double *d, int m, int n)
: data(d, m*n), rows(m), cols(n), ld(m) {}
: data(d, m*n), rows(m), cols(n), ld(m)
{
}
ConstGeneralMatrix(const GeneralMatrix &m);
ConstGeneralMatrix(const GeneralMatrix &m, int i, int j, int nrows, int ncols);
ConstGeneralMatrix(const ConstGeneralMatrix &m, int i, int j, int nrows, int ncols);
virtual ~ConstGeneralMatrix() {}
virtual ~ConstGeneralMatrix()
{
}
const double& get(int i, int j) const
{return data[j*ld+i];}
int numRows() const {return rows;}
int numCols() const {return cols;}
int getLD() const {return ld;}
const double* base() const {return data.base();}
const ConstVector& getData() const {return data;}
const double &
get(int i, int j) const
{
return data[j*ld+i];
}
int
numRows() const
{
return rows;
}
int
numCols() const
{
return cols;
}
int
getLD() const
{
return ld;
}
const double *
base() const
{
return data.base();
}
const ConstVector &
getData() const
{
return data;
}
double getNormInf() const;
double getNorm1() const;
@ -41,17 +69,29 @@ public:
/* x = scalar(a)*x + scalar(b)*this'*d */
void multVecTrans(double a, Vector &x, double b, const ConstVector &d) const;
/* x = x + this*d */
void multaVec(Vector& x, const ConstVector& d) const
{multVec(1.0, x, 1.0, d);}
void
multaVec(Vector &x, const ConstVector &d) const
{
multVec(1.0, x, 1.0, d);
}
/* x = x + this'*d */
void multaVecTrans(Vector& x, const ConstVector& d) const
{multVecTrans(1.0, x, 1.0, d);}
void
multaVecTrans(Vector &x, const ConstVector &d) const
{
multVecTrans(1.0, x, 1.0, d);
}
/* x = x - this*d */
void multsVec(Vector& x, const ConstVector& d) const
{multVec(1.0, x, -1.0, d);}
void
multsVec(Vector &x, const ConstVector &d) const
{
multVec(1.0, x, -1.0, d);
}
/* x = x - this'*d */
void multsVecTrans(Vector& x, const ConstVector& d) const
{multVecTrans(1.0, x, -1.0, d);}
void
multsVecTrans(Vector &x, const ConstVector &d) const
{
multVecTrans(1.0, x, -1.0, d);
}
/* m = inv(this)*m */
void multInvLeft(GeneralMatrix &m) const;
/* m = inv(this')*m */
@ -70,8 +110,8 @@ protected:
void multInvLeft(const char *trans, int mrows, int mcols, int mld, double *d) const;
};
class GeneralMatrix {
class GeneralMatrix
{
friend class ConstGeneralMatrix;
protected:
Vector data;
@ -80,11 +120,17 @@ protected:
int ld;
public:
GeneralMatrix(int m, int n)
: data(m*n), rows(m), cols(n), ld(m) {}
: data(m*n), rows(m), cols(n), ld(m)
{
}
GeneralMatrix(const double *d, int m, int n)
: data(d, m*n), rows(m), cols(n), ld(m) {}
: data(d, m*n), rows(m), cols(n), ld(m)
{
}
GeneralMatrix(double *d, int m, int n)
: data(d, m*n), rows(m), cols(n), ld(m) {}
: data(d, m*n), rows(m), cols(n), ld(m)
{
}
GeneralMatrix(const GeneralMatrix &m);
GeneralMatrix(const ConstGeneralMatrix &m);
GeneralMatrix(const GeneralMatrix &m, const char *dummy); // transpose
@ -101,113 +147,208 @@ public:
GeneralMatrix(const GeneralMatrix &a, const char *dum1,
const GeneralMatrix &b, const char *dum2);
virtual ~GeneralMatrix();
const GeneralMatrix& operator=(const GeneralMatrix& m)
{data=m.data; rows=m.rows; cols=m.cols; ld=m.ld; return *this;}
virtual
~GeneralMatrix();
const GeneralMatrix &
operator=(const GeneralMatrix &m)
{
data = m.data; rows = m.rows; cols = m.cols; ld = m.ld; return *this;
}
const double& get(int i, int j) const
{return data[j*ld+i];}
double& get(int i, int j)
{return data[j*ld+i];}
int numRows() const {return rows;}
int numCols() const {return cols;}
int getLD() const {return ld;}
double* base() {return data.base();}
const double* base() const {return data.base();}
Vector& getData() {return data;}
const Vector& getData() const {return data;}
const double &
get(int i, int j) const
{
return data[j*ld+i];
}
double &
get(int i, int j)
{
return data[j*ld+i];
}
int
numRows() const
{
return rows;
}
int
numCols() const
{
return cols;
}
int
getLD() const
{
return ld;
}
double *
base()
{
return data.base();
}
const double *
base() const
{
return data.base();
}
Vector &
getData()
{
return data;
}
const Vector &
getData() const
{
return data;
}
double getNormInf() const
{return ConstGeneralMatrix(*this).getNormInf();}
double getNorm1() const
{return ConstGeneralMatrix(*this).getNorm1();}
double
getNormInf() const
{
return ConstGeneralMatrix(*this).getNormInf();
}
double
getNorm1() const
{
return ConstGeneralMatrix(*this).getNorm1();
}
/* place matrix m to the position (i,j) */
void place(const ConstGeneralMatrix &m, int i, int j);
void place(const GeneralMatrix& m, int i, int j)
{place(ConstGeneralMatrix(m), i, j);}
void
place(const GeneralMatrix &m, int i, int j)
{
place(ConstGeneralMatrix(m), i, j);
}
/* this = a*b */
void mult(const ConstGeneralMatrix &a, const ConstGeneralMatrix &b);
void mult(const GeneralMatrix& a, const GeneralMatrix& b)
{mult(ConstGeneralMatrix(a), ConstGeneralMatrix(b));}
void
mult(const GeneralMatrix &a, const GeneralMatrix &b)
{
mult(ConstGeneralMatrix(a), ConstGeneralMatrix(b));
}
/* this = this + scalar*a*b */
void multAndAdd(const ConstGeneralMatrix &a, const ConstGeneralMatrix &b,
double mult = 1.0);
void multAndAdd(const GeneralMatrix& a, const GeneralMatrix& b,
void
multAndAdd(const GeneralMatrix &a, const GeneralMatrix &b,
double mult = 1.0)
{multAndAdd(ConstGeneralMatrix(a), ConstGeneralMatrix(b), mult);}
{
multAndAdd(ConstGeneralMatrix(a), ConstGeneralMatrix(b), mult);
}
/* this = this + scalar*a*b' */
void multAndAdd(const ConstGeneralMatrix &a, const ConstGeneralMatrix &b,
const char *dum, double mult = 1.0);
void multAndAdd(const GeneralMatrix& a, const GeneralMatrix& b,
void
multAndAdd(const GeneralMatrix &a, const GeneralMatrix &b,
const char *dum, double mult = 1.0)
{multAndAdd(ConstGeneralMatrix(a), ConstGeneralMatrix(b), dum, mult);}
{
multAndAdd(ConstGeneralMatrix(a), ConstGeneralMatrix(b), dum, mult);
}
/* this = this + scalar*a'*b */
void multAndAdd(const ConstGeneralMatrix &a, const char *dum, const ConstGeneralMatrix &b,
double mult = 1.0);
void multAndAdd(const GeneralMatrix& a, const char* dum, const GeneralMatrix& b,
void
multAndAdd(const GeneralMatrix &a, const char *dum, const GeneralMatrix &b,
double mult = 1.0)
{multAndAdd(ConstGeneralMatrix(a), dum, ConstGeneralMatrix(b), mult);}
{
multAndAdd(ConstGeneralMatrix(a), dum, ConstGeneralMatrix(b), mult);
}
/* this = this + scalar*a'*b' */
void multAndAdd(const ConstGeneralMatrix &a, const char *dum1,
const ConstGeneralMatrix &b, const char *dum2, double mult = 1.0);
void multAndAdd(const GeneralMatrix& a, const char* dum1,
void
multAndAdd(const GeneralMatrix &a, const char *dum1,
const GeneralMatrix &b, const char *dum2, double mult = 1.0)
{multAndAdd(ConstGeneralMatrix(a), dum1, ConstGeneralMatrix(b),dum2, mult);}
{
multAndAdd(ConstGeneralMatrix(a), dum1, ConstGeneralMatrix(b), dum2, mult);
}
/* this = this + scalar*a*a' */
void addOuter(const ConstVector &a, double mult = 1.0);
void addOuter(const Vector& a, double mult=1.0)
{addOuter(ConstVector(a), mult);}
void
addOuter(const Vector &a, double mult = 1.0)
{
addOuter(ConstVector(a), mult);
}
/* this = this * m */
void multRight(const ConstGeneralMatrix &m);
void multRight(const GeneralMatrix& m)
{multRight(ConstGeneralMatrix(m));}
void
multRight(const GeneralMatrix &m)
{
multRight(ConstGeneralMatrix(m));
}
/* this = m * this */
void multLeft(const ConstGeneralMatrix &m);
void multLeft(const GeneralMatrix& m)
{multLeft(ConstGeneralMatrix(m));}
void
multLeft(const GeneralMatrix &m)
{
multLeft(ConstGeneralMatrix(m));
}
/* this = this * m' */
void multRightTrans(const ConstGeneralMatrix &m);
void multRightTrans(const GeneralMatrix& m)
{multRightTrans(ConstGeneralMatrix(m));}
void
multRightTrans(const GeneralMatrix &m)
{
multRightTrans(ConstGeneralMatrix(m));
}
/* this = m' * this */
void multLeftTrans(const ConstGeneralMatrix &m);
void multLeftTrans(const GeneralMatrix& m)
{multLeftTrans(ConstGeneralMatrix(m));}
void
multLeftTrans(const GeneralMatrix &m)
{
multLeftTrans(ConstGeneralMatrix(m));
}
/* x = scalar(a)*x + scalar(b)*this*d */
void multVec(double a, Vector& x, double b, const ConstVector& d) const
{ConstGeneralMatrix(*this).multVec(a, x, b, d);}
void
multVec(double a, Vector &x, double b, const ConstVector &d) const
{
ConstGeneralMatrix(*this).multVec(a, x, b, d);
}
/* x = scalar(a)*x + scalar(b)*this'*d */
void multVecTrans(double a, Vector& x, double b, const ConstVector& d) const
{ConstGeneralMatrix(*this).multVecTrans(a, x, b, d);}
void
multVecTrans(double a, Vector &x, double b, const ConstVector &d) const
{
ConstGeneralMatrix(*this).multVecTrans(a, x, b, d);
}
/* x = x + this*d */
void multaVec(Vector& x, const ConstVector& d) const
{ConstGeneralMatrix(*this).multaVec(x, d);}
void
multaVec(Vector &x, const ConstVector &d) const
{
ConstGeneralMatrix(*this).multaVec(x, d);
}
/* x = x + this'*d */
void multaVecTrans(Vector& x, const ConstVector& d) const
{ConstGeneralMatrix(*this).multaVecTrans(x, d);}
void
multaVecTrans(Vector &x, const ConstVector &d) const
{
ConstGeneralMatrix(*this).multaVecTrans(x, d);
}
/* x = x - this*d */
void multsVec(Vector& x, const ConstVector& d) const
{ConstGeneralMatrix(*this).multsVec(x, d);}
void
multsVec(Vector &x, const ConstVector &d) const
{
ConstGeneralMatrix(*this).multsVec(x, d);
}
/* x = x - this'*d */
void multsVecTrans(Vector& x, const ConstVector& d) const
{ConstGeneralMatrix(*this).multsVecTrans(x, d);}
void
multsVecTrans(Vector &x, const ConstVector &d) const
{
ConstGeneralMatrix(*this).multsVecTrans(x, d);
}
/* this = zero */
void zeros();
@ -226,55 +367,83 @@ public:
/* this = this + scalar*m */
void add(double a, const ConstGeneralMatrix &m);
void add(double a, const GeneralMatrix& m)
{add(a, ConstGeneralMatrix(m));}
void
add(double a, const GeneralMatrix &m)
{
add(a, ConstGeneralMatrix(m));
}
/* this = this + scalar*m' */
void add(double a, const ConstGeneralMatrix &m, const char *dum);
void add(double a, const GeneralMatrix& m, const char* dum)
{add(a, ConstGeneralMatrix(m), dum);}
void
add(double a, const GeneralMatrix &m, const char *dum)
{
add(a, ConstGeneralMatrix(m), dum);
}
bool isFinite() const
{return (ConstGeneralMatrix(*this)).isFinite();}
bool
isFinite() const
{
return (ConstGeneralMatrix(*this)).isFinite();
}
bool isZero() const
{return (ConstGeneralMatrix(*this)).isZero();}
bool
isZero() const
{
return (ConstGeneralMatrix(*this)).isZero();
}
virtual void print() const
{ConstGeneralMatrix(*this).print();}
virtual void
print() const
{
ConstGeneralMatrix(*this).print();
}
private:
void copy(const ConstGeneralMatrix &m, int ioff = 0, int joff = 0);
void copy(const GeneralMatrix& m, int ioff = 0, int joff = 0)
{copy(ConstGeneralMatrix(m), ioff, joff);}
void
copy(const GeneralMatrix &m, int ioff = 0, int joff = 0)
{
copy(ConstGeneralMatrix(m), ioff, joff);
}
void gemm(const char *transa, const ConstGeneralMatrix &a,
const char *transb, const ConstGeneralMatrix &b,
double alpha, double beta);
void gemm(const char* transa, const GeneralMatrix& a,
void
gemm(const char *transa, const GeneralMatrix &a,
const char *transb, const GeneralMatrix &b,
double alpha, double beta)
{gemm(transa, ConstGeneralMatrix(a), transb, ConstGeneralMatrix(b),
alpha, beta);}
{
gemm(transa, ConstGeneralMatrix(a), transb, ConstGeneralMatrix(b),
alpha, beta);
}
/* this = this * op(m) (without whole copy of this) */
void gemm_partial_right(const char *trans, const ConstGeneralMatrix &m,
double alpha, double beta);
void gemm_partial_right(const char* trans, const GeneralMatrix& m,
void
gemm_partial_right(const char *trans, const GeneralMatrix &m,
double alpha, double beta)
{gemm_partial_right(trans, ConstGeneralMatrix(m), alpha, beta);}
{
gemm_partial_right(trans, ConstGeneralMatrix(m), alpha, beta);
}
/* this = op(m) *this (without whole copy of this) */
void gemm_partial_left(const char *trans, const ConstGeneralMatrix &m,
double alpha, double beta);
void gemm_partial_left(const char* trans, const GeneralMatrix& m,
void
gemm_partial_left(const char *trans, const GeneralMatrix &m,
double alpha, double beta)
{gemm_partial_left(trans, ConstGeneralMatrix(m), alpha, beta);}
{
gemm_partial_left(trans, ConstGeneralMatrix(m), alpha, beta);
}
/* number of rows/columns for copy used in gemm_partial_* */
static int md_length;
};
class SVDDecomp {
class SVDDecomp
{
protected:
/** Minimum of number of rows and columns of the decomposed
* matrix. */
@ -294,13 +463,22 @@ public:
U(A.numRows(), A.numRows()),
VT(A.numCols(), A.numCols()),
conv(false)
{construct(A);}
const GeneralMatrix& getU() const
{return U;}
const GeneralMatrix& getVT() const
{return VT;}
{
construct(A);
}
const GeneralMatrix &
getU() const
{
return U;
}
const GeneralMatrix &
getVT() const
{
return VT;
}
void solve(const GeneralMatrix &B, GeneralMatrix &X) const;
void solve(const Vector& b, Vector& x) const
void
solve(const Vector &b, Vector &x) const
{
GeneralMatrix xmat(x.base(), x.length(), 1);
solve(GeneralMatrix(b.base(), b.length(), 1), xmat);
@ -309,12 +487,8 @@ private:
void construct(const GeneralMatrix &A);
};
#endif /* GENERAL_MATRIX_H */
// Local Variables:
// mode:C++
// End:

View File

@ -10,7 +10,8 @@
#include "SimilarityDecomp.h"
#include "SylvesterSolver.h"
class GeneralSylvester {
class GeneralSylvester
{
SylvParams pars;
SylvMemoryDriver mem_driver;
int order;
@ -41,12 +42,33 @@ public:
const double *da, const double *db,
const double *dc, double *dd,
const SylvParams &ps);
virtual ~GeneralSylvester();
int getM() const {return c.numRows();}
int getN() const {return a.numRows();}
const double* getResult() const {return d.base();}
const SylvParams& getParams() const {return pars;}
SylvParams& getParams() {return pars;}
virtual
~GeneralSylvester();
int
getM() const
{
return c.numRows();
}
int
getN() const
{
return a.numRows();
}
const double *
getResult() const
{
return d.base();
}
const SylvParams &
getParams() const
{
return pars;
}
SylvParams &
getParams()
{
return pars;
}
void solve();
void check(const double *ds);
private:
@ -55,7 +77,6 @@ private:
#endif /* GENERAL_SYLVESTER_H */
// Local Variables:
// mode:C++
// End:

View File

@ -10,14 +10,21 @@
#include "QuasiTriangular.h"
#include "SimilarityDecomp.h"
class IterativeSylvester : public SylvesterSolver {
class IterativeSylvester : public SylvesterSolver
{
public:
IterativeSylvester(const QuasiTriangular &k, const QuasiTriangular &f)
: SylvesterSolver(k, f) {}
: SylvesterSolver(k, f)
{
}
IterativeSylvester(const SchurDecompZero &kdecomp, const SchurDecomp &fdecomp)
: SylvesterSolver(kdecomp, fdecomp) {}
: SylvesterSolver(kdecomp, fdecomp)
{
}
IterativeSylvester(const SchurDecompZero &kdecomp, const SimilarityDecomp &fdecomp)
: SylvesterSolver(kdecomp, fdecomp) {}
: SylvesterSolver(kdecomp, fdecomp)
{
}
void solve(SylvParams &pars, KronVector &x) const;
private:
double performFirstStep(KronVector &x) const;
@ -27,7 +34,6 @@ private:
#endif /* ITERATIVE_SYLVESTER_H */
// Local Variables:
// mode:C++
// End:

View File

@ -8,7 +8,8 @@
#include "KronVector.h"
#include "QuasiTriangular.h"
class KronUtils {
class KronUtils
{
public:
/* multiplies I_m\otimes..\I_m\otimes T\otimes I_m...I_m\otimes I_n
with given b and returns x. T must be (m,m), number of

View File

@ -9,26 +9,47 @@
class ConstKronVector;
class KronVector : public Vector {
class KronVector : public Vector
{
protected:
int m;
int n;
int depth;
public:
KronVector() : Vector((double*)0, 0), m(0), n(0), depth(0) {}
KronVector() : Vector((double *) 0, 0), m(0), n(0), depth(0)
{
}
KronVector(int mm, int nn, int dp); // new instance
KronVector(Vector &v, int mm, int nn, int dp); // conversion
KronVector(KronVector &, int i); // picks i-th subvector
KronVector(const ConstKronVector &v); // new instance and copy
const KronVector& operator=(KronVector& v)
{Vector::operator=(v); m=v.m; n=v.n; depth = v.depth; return *this;}
const KronVector& operator=(const KronVector& v)
{Vector::operator=(v); m=v.m; n=v.n; depth = v.depth; return *this;}
const KronVector &
operator=(KronVector &v)
{
Vector::operator=(v); m = v.m; n = v.n; depth = v.depth; return *this;
}
const KronVector &
operator=(const KronVector &v)
{
Vector::operator=(v); m = v.m; n = v.n; depth = v.depth; return *this;
}
const KronVector &operator=(const ConstKronVector &v);
const KronVector &operator=(const Vector &v);
int getM() const {return m;}
int getN() const {return n;}
int getDepth() const {return depth;}
int
getM() const
{
return m;
}
int
getN() const
{
return n;
}
int
getDepth() const
{
return depth;
}
};
class ConstKronVector : public ConstVector
@ -44,9 +65,21 @@ public:
ConstKronVector(const ConstVector &v, int mm, int nn, int dp);
ConstKronVector(const KronVector &v, int i);
ConstKronVector(const ConstKronVector &v, int i);
int getM() const {return m;}
int getN() const {return n;}
int getDepth() const {return depth;}
int
getM() const
{
return m;
}
int
getN() const
{
return n;
}
int
getDepth() const
{
return depth;
}
};
int power(int m, int depth);

View File

@ -15,24 +15,46 @@ using namespace std;
class DiagonalBlock;
class Diagonal;
class DiagPair {
class DiagPair
{
private:
double *a1;
double *a2;
public:
DiagPair() {}
DiagPair(double* aa1, double* aa2) {a1 = aa1; a2 = aa2;}
DiagPair(const DiagPair& p) {a1 = p.a1; a2 = p.a2;}
const DiagPair& operator=(const DiagPair& p) {a1 = p.a1; a2 = p.a2; return *this;}
const DiagPair& operator=(double v) {*a1 = v; *a2 = v; return *this;}
const double& operator*() const {return *a1;}
DiagPair()
{
}
DiagPair(double *aa1, double *aa2)
{
a1 = aa1; a2 = aa2;
}
DiagPair(const DiagPair &p)
{
a1 = p.a1; a2 = p.a2;
}
const DiagPair &
operator=(const DiagPair &p)
{
a1 = p.a1; a2 = p.a2; return *this;
}
const DiagPair &
operator=(double v)
{
*a1 = v; *a2 = v; return *this;
}
const double &
operator*() const
{
return *a1;
}
/** here we must not define double& operator*(), since it wouldn't
rewrite both values, we use operator= for this */
friend class Diagonal;
friend class DiagonalBlock;
};
class DiagonalBlock {
class DiagonalBlock
{
private:
int jbar;
bool real;
@ -40,7 +62,9 @@ private:
double *beta1;
double *beta2;
void copy(const DiagonalBlock& b) {
void
copy(const DiagonalBlock &b)
{
jbar = b.jbar;
real = b.real;
alpha = b.alpha;
@ -49,7 +73,9 @@ private:
}
public:
DiagonalBlock() {}
DiagonalBlock()
{
}
DiagonalBlock(int jb, bool r, double *a1, double *a2,
double *b1, double *b2)
: alpha(a1, a2)
@ -78,21 +104,44 @@ public:
beta2 = 0;
}
DiagonalBlock(const DiagonalBlock &b)
{copy(b);}
const DiagonalBlock& operator=(const DiagonalBlock& b)
{copy(b); return *this;}
int getIndex() const
{return jbar;}
bool isReal() const
{return real;}
const DiagPair& getAlpha() const
{return alpha;}
DiagPair& getAlpha()
{return alpha;}
double& getBeta1() const
{return *beta1;}
double& getBeta2() const
{return *beta2;}
{
copy(b);
}
const DiagonalBlock &
operator=(const DiagonalBlock &b)
{
copy(b); return *this;
}
int
getIndex() const
{
return jbar;
}
bool
isReal() const
{
return real;
}
const DiagPair &
getAlpha() const
{
return alpha;
}
DiagPair &
getAlpha()
{
return alpha;
}
double &
getBeta1() const
{
return *beta1;
}
double &
getBeta2() const
{
return *beta2;
}
double getDeterminant() const;
double getSBeta() const;
double getSize() const;
@ -103,22 +152,54 @@ public:
};
template <class _Tdiag, class _Tblock, class _Titer>
struct _diag_iter {
struct _diag_iter
{
typedef _diag_iter<_Tdiag, _Tblock, _Titer> _Self;
_Tdiag diag;
_Titer it;
public:
_diag_iter(_Tdiag d, _Titer iter) : diag(d), it(iter) {}
_Tblock operator*() const {return *it;}
_Self& operator++() {++it; return *this;}
_Self& operator--() {--it; return *this;}
bool operator==(const _Self& x) const {return x.it == it;}
bool operator!=(const _Self& x) const {return x.it != it;}
const _Self& operator=(const _Self& x) {it = x.it; return *this;}
_Titer iter() const {return it;}
_diag_iter(_Tdiag d, _Titer iter) : diag(d), it(iter)
{
}
_Tblock
operator*() const
{
return *it;
}
_Self &
operator++()
{
++it; return *this;
}
_Self &
operator--()
{
--it; return *this;
}
bool
operator==(const _Self &x) const
{
return x.it == it;
}
bool
operator!=(const _Self &x) const
{
return x.it != it;
}
const _Self &
operator=(const _Self &x)
{
it = x.it; return *this;
}
_Titer
iter() const
{
return it;
}
};
class Diagonal {
class Diagonal
{
public:
typedef _diag_iter<const Diagonal &, const DiagonalBlock &, list<DiagonalBlock>::const_iterator> const_diag_iter;
typedef _diag_iter<Diagonal &, DiagonalBlock &, list<DiagonalBlock>::iterator> diag_iter;
@ -128,17 +209,44 @@ private:
int num_real;
void copy(const Diagonal &);
public:
Diagonal() : num_all(0), num_real(0) {}
Diagonal() : num_all(0), num_real(0)
{
}
Diagonal(double *data, int d_size);
Diagonal(double *data, const Diagonal &d);
Diagonal(const Diagonal& d) {copy(d);}
const Diagonal& operator =(const Diagonal& d) {copy(d); return *this;}
virtual ~Diagonal() {}
Diagonal(const Diagonal &d)
{
copy(d);
}
const Diagonal &
operator=(const Diagonal &d)
{
copy(d); return *this;
}
virtual ~Diagonal()
{
}
int getNumComplex() const {return num_all - num_real;}
int getNumReal() const {return num_real;}
int getSize() const {return getNumReal() + 2*getNumComplex();}
int getNumBlocks() const {return num_all;}
int
getNumComplex() const
{
return num_all - num_real;
}
int
getNumReal() const
{
return num_real;
}
int
getSize() const
{
return getNumReal() + 2*getNumComplex();
}
int
getNumBlocks() const
{
return num_all;
}
void getEigenValues(Vector &eig) const;
void swapLogically(diag_iter it);
void checkConsistency(diag_iter it);
@ -147,14 +255,26 @@ public:
diag_iter findNextLargerBlock(diag_iter start, diag_iter end, double a);
void print() const;
diag_iter begin()
{return diag_iter(*this, blocks.begin());}
const_diag_iter begin() const
{return const_diag_iter(*this, blocks.begin());}
diag_iter end()
{return diag_iter(*this, blocks.end());}
const_diag_iter end() const
{return const_diag_iter(*this, blocks.end());}
diag_iter
begin()
{
return diag_iter(*this, blocks.begin());
}
const_diag_iter
begin() const
{
return const_diag_iter(*this, blocks.begin());
}
diag_iter
end()
{
return diag_iter(*this, blocks.end());
}
const_diag_iter
end() const
{
return const_diag_iter(*this, blocks.end());
}
/* redefine pointers as data start at p */
void changeBase(double *p);
@ -165,74 +285,123 @@ private:
};
template <class _TRef, class _TPtr>
struct _matrix_iter {
struct _matrix_iter
{
typedef _matrix_iter<_TRef, _TPtr> _Self;
int d_size;
bool real;
_TPtr ptr;
public:
_matrix_iter(_TPtr base, int ds, bool r)
{ptr = base; d_size = ds; real = r;}
virtual ~_matrix_iter() {}
const _Self& operator=(const _Self& it)
{ptr = it.ptr; d_size = it.d_size; real = it.real; return *this;}
bool operator==(const _Self& it) const
{return ptr == it.ptr;}
bool operator!=(const _Self& it) const
{return ptr != it.ptr;}
_TRef operator*() const
{return *ptr;}
_TRef a() const
{return *ptr;}
{
ptr = base; d_size = ds; real = r;
}
virtual ~_matrix_iter()
{
}
const _Self &
operator=(const _Self &it)
{
ptr = it.ptr; d_size = it.d_size; real = it.real; return *this;
}
bool
operator==(const _Self &it) const
{
return ptr == it.ptr;
}
bool
operator!=(const _Self &it) const
{
return ptr != it.ptr;
}
_TRef
operator*() const
{
return *ptr;
}
_TRef
a() const
{
return *ptr;
}
virtual _Self &operator++() = 0;
};
template <class _TRef, class _TPtr>
class _column_iter : public _matrix_iter<_TRef, _TPtr> {
class _column_iter : public _matrix_iter<_TRef, _TPtr>
{
typedef _matrix_iter<_TRef, _TPtr> _Tparent;
typedef _column_iter<_TRef, _TPtr> _Self;
int row;
public:
_column_iter(_TPtr base, int ds, bool r, int rw)
: _matrix_iter<_TRef, _TPtr>(base, ds, r), row(rw) {};
_Self& operator++()
{_Tparent::ptr++; row++; return *this;}
_TRef b() const
: _matrix_iter<_TRef, _TPtr>(base, ds, r), row(rw)
{
};
_Self &
operator++()
{
_Tparent::ptr++; row++; return *this;
}
_TRef
b() const
{
if (_Tparent::real)
{
if (_Tparent::real) {
return *(_Tparent::ptr);
} else {
}
else
{
return *(_Tparent::ptr+_Tparent::d_size);
}
}
int getRow() const {return row;}
int
getRow() const
{
return row;
}
};
template <class _TRef, class _TPtr>
class _row_iter : public _matrix_iter<_TRef, _TPtr> {
class _row_iter : public _matrix_iter<_TRef, _TPtr>
{
typedef _matrix_iter<_TRef, _TPtr> _Tparent;
typedef _row_iter<_TRef, _TPtr> _Self;
int col;
public:
_row_iter(_TPtr base, int ds, bool r, int cl)
: _matrix_iter<_TRef, _TPtr>(base, ds, r), col(cl) {};
_Self& operator++()
{_Tparent::ptr += _Tparent::d_size; col++; return *this;}
virtual _TRef b() const
: _matrix_iter<_TRef, _TPtr>(base, ds, r), col(cl)
{
};
_Self &
operator++()
{
_Tparent::ptr += _Tparent::d_size; col++; return *this;
}
virtual _TRef
b() const
{
if (_Tparent::real)
{
if (_Tparent::real) {
return *(_Tparent::ptr);
}else {
}
else
{
return *(_Tparent::ptr+1);
}
}
int getCol() const {return col;}
int
getCol() const
{
return col;
}
};
class SchurDecomp;
class SchurDecompZero;
class QuasiTriangular : public SqSylvMatrix {
class QuasiTriangular : public SqSylvMatrix
{
public:
typedef _column_iter<const double &, const double *> const_col_iter;
typedef _column_iter<double &, double *> col_iter;
@ -251,8 +420,13 @@ public:
QuasiTriangular(const SchurDecomp &decomp);
QuasiTriangular(const SchurDecompZero &decomp);
QuasiTriangular(const QuasiTriangular &t);
virtual ~QuasiTriangular();
const Diagonal& getDiagonal() const {return diagonal;}
virtual
~QuasiTriangular();
const Diagonal &
getDiagonal() const
{
return diagonal;
}
int getNumOffdiagonal() const;
void swapDiagLogically(diag_iter it);
void checkDiagConsistency(diag_iter it);
@ -260,7 +434,6 @@ public:
diag_iter findClosestDiagBlock(diag_iter start, diag_iter end, double a);
diag_iter findNextLargerBlock(diag_iter start, diag_iter end, double a);
/* (I+T)y = x, y-->x */
virtual void solvePre(Vector &x, double &eig_min);
/* (I+T')y = x, y-->x */
@ -286,14 +459,26 @@ public:
/* A = T'*A */
virtual void multLeftOtherTrans(GeneralMatrix &a) const;
const_diag_iter diag_begin() const
{return diagonal.begin();}
diag_iter diag_begin()
{return diagonal.begin();}
const_diag_iter diag_end() const
{return diagonal.end();}
diag_iter diag_end()
{return diagonal.end();}
const_diag_iter
diag_begin() const
{
return diagonal.begin();
}
diag_iter
diag_begin()
{
return diagonal.begin();
}
const_diag_iter
diag_end() const
{
return diagonal.end();
}
diag_iter
diag_end()
{
return diagonal.end();
}
/* iterators for off diagonal elements */
virtual const_col_iter col_begin(const DiagonalBlock &b) const;
@ -306,14 +491,26 @@ public:
virtual row_iter row_end(const DiagonalBlock &b);
/* clone */
virtual QuasiTriangular* clone() const
{return new QuasiTriangular(*this);}
virtual QuasiTriangular* clone(int p, const QuasiTriangular& t) const
{return new QuasiTriangular(p, t);}
virtual QuasiTriangular* clone(double r) const
{return new QuasiTriangular(r, *this);}
virtual QuasiTriangular* clone(double r, double rr, const QuasiTriangular& tt) const
{return new QuasiTriangular(r, *this, rr, tt);}
virtual QuasiTriangular *
clone() const
{
return new QuasiTriangular(*this);
}
virtual QuasiTriangular *
clone(int p, const QuasiTriangular &t) const
{
return new QuasiTriangular(p, t);
}
virtual QuasiTriangular *
clone(double r) const
{
return new QuasiTriangular(r, *this);
}
virtual QuasiTriangular *
clone(double r, double rr, const QuasiTriangular &tt) const
{
return new QuasiTriangular(r, *this, rr, tt);
}
protected:
void setMatrix(double r, const QuasiTriangular &t);
void addMatrix(double r, const QuasiTriangular &t);
@ -333,7 +530,6 @@ private:
#endif /* QUASI_TRIANGULAR_H */
// Local Variables:
// mode:C++
// End:

View File

@ -8,7 +8,8 @@
#include "QuasiTriangular.h"
#include "GeneralMatrix.h"
class QuasiTriangularZero : public QuasiTriangular {
class QuasiTriangularZero : public QuasiTriangular
{
int nz; // number of zero columns
GeneralMatrix ru; // data in right upper part (nz,d_size)
public:
@ -30,14 +31,26 @@ public:
void multKronTrans(KronVector &x) const;
void multLeftOther(GeneralMatrix &a) const;
/* clone */
virtual QuasiTriangular* clone() const
{return new QuasiTriangularZero(*this);}
virtual QuasiTriangular* clone(int p, const QuasiTriangular& t) const
{return new QuasiTriangularZero(p, (const QuasiTriangularZero&)t);}
virtual QuasiTriangular* clone(double r) const
{return new QuasiTriangularZero(r, *this);}
virtual QuasiTriangular* clone(double r, double rr, const QuasiTriangular& tt) const
{return new QuasiTriangularZero(r, *this, rr, (const QuasiTriangularZero&)tt);}
virtual QuasiTriangular *
clone() const
{
return new QuasiTriangularZero(*this);
}
virtual QuasiTriangular *
clone(int p, const QuasiTriangular &t) const
{
return new QuasiTriangularZero(p, (const QuasiTriangularZero &) t);
}
virtual QuasiTriangular *
clone(double r) const
{
return new QuasiTriangularZero(r, *this);
}
virtual QuasiTriangular *
clone(double r, double rr, const QuasiTriangular &tt) const
{
return new QuasiTriangularZero(r, *this, rr, (const QuasiTriangularZero &) tt);
}
void print() const;
};

View File

@ -9,7 +9,8 @@
#include "QuasiTriangular.h"
class QuasiTriangular;
class SchurDecomp {
class SchurDecomp
{
bool q_destroy;
SqSylvMatrix *q;
bool t_destroy;
@ -18,26 +19,51 @@ public:
SchurDecomp(const SqSylvMatrix &m);
SchurDecomp(const QuasiTriangular &tr);
SchurDecomp(QuasiTriangular &tr);
const SqSylvMatrix& getQ() const {return *q;}
const QuasiTriangular& getT() const {return *t;}
SqSylvMatrix& getQ() {return *q;}
QuasiTriangular& getT() {return *t;}
const SqSylvMatrix &
getQ() const
{
return *q;
}
const QuasiTriangular &
getT() const
{
return *t;
}
SqSylvMatrix &
getQ()
{
return *q;
}
QuasiTriangular &
getT()
{
return *t;
}
virtual int getDim() const;
virtual ~SchurDecomp();
virtual
~SchurDecomp();
};
class SchurDecompZero : public SchurDecomp {
class SchurDecompZero : public SchurDecomp
{
GeneralMatrix ru; /* right upper matrix */
public:
SchurDecompZero(const GeneralMatrix &m);
const GeneralMatrix& getRU() const {return ru;}
const GeneralMatrix &
getRU() const
{
return ru;
}
int getDim() const;
int getZeroCols() const {return ru.numRows();}
int
getZeroCols() const
{
return ru.numRows();
}
};
#endif /* SCHUR_DECOMP_H */
// Local Variables:
// mode:C++
// End:

View File

@ -10,12 +10,19 @@
#include "SchurDecomp.h"
#include "QuasiTriangular.h"
class SchurDecompEig : public SchurDecomp {
class SchurDecompEig : public SchurDecomp
{
public:
typedef QuasiTriangular::diag_iter diag_iter;
SchurDecompEig(const SqSylvMatrix& m) : SchurDecomp(m) {}
SchurDecompEig(const QuasiTriangular& tr) : SchurDecomp(tr) {};
SchurDecompEig(QuasiTriangular& tr) : SchurDecomp(tr) {}
SchurDecompEig(const SqSylvMatrix &m) : SchurDecomp(m)
{
}
SchurDecompEig(const QuasiTriangular &tr) : SchurDecomp(tr)
{
};
SchurDecompEig(QuasiTriangular &tr) : SchurDecomp(tr)
{
}
diag_iter bubbleEigen(diag_iter from, diag_iter to);
void orderEigen();
protected:
@ -24,8 +31,6 @@ protected:
#endif /* SCHUR_DECOMP_EIG_H */
// Local Variables:
// mode:C++
// End:

View File

@ -9,20 +9,31 @@
#include "BlockDiagonal.h"
#include "SylvParams.h"
class SimilarityDecomp {
class SimilarityDecomp
{
SqSylvMatrix *q;
BlockDiagonal *b;
SqSylvMatrix *invq;
typedef BlockDiagonal::diag_iter diag_iter;
public:
SimilarityDecomp(const double *d, int d_size, double log10norm = 3.0);
virtual ~SimilarityDecomp();
const SqSylvMatrix& getQ() const
{return *q;}
const SqSylvMatrix& getInvQ() const
{return *invq;}
const BlockDiagonal& getB() const
{return *b;}
virtual
~SimilarityDecomp();
const SqSylvMatrix &
getQ() const
{
return *q;
}
const SqSylvMatrix &
getInvQ() const
{
return *invq;
}
const BlockDiagonal &
getB() const
{
return *b;
}
void check(SylvParams &pars, const GeneralMatrix &m) const;
void infoToPars(SylvParams &pars) const;
protected:
@ -35,7 +46,6 @@ protected:
#endif /* SIMILARITY_DECOMP_H */
// Local Variables:
// mode:C++
// End:

View File

@ -7,20 +7,22 @@
#include "SylvMemory.h"
class SylvException : public MallocAllocator {
class SylvException : public MallocAllocator
{
protected:
char file[50];
int line;
const SylvException *source;
public:
SylvException(const char *f, int l, const SylvException *s);
virtual ~SylvException();
virtual
~SylvException();
virtual int printMessage(char *str, int maxlen) const;
void printMessage() const;
};
class SylvExceptionMessage : public SylvException {
class SylvExceptionMessage : public SylvException
{
char message[500];
public:
SylvExceptionMessage(const char *f, int l, const char *mes);
@ -33,7 +35,6 @@ public:
#endif /* SYLV_EXCEPTION_H */
// Local Variables:
// mode:C++
// End:

View File

@ -10,22 +10,37 @@
class SqSylvMatrix;
class SylvMatrix : public GeneralMatrix {
class SylvMatrix : public GeneralMatrix
{
public:
SylvMatrix(int m, int n)
: GeneralMatrix(m,n) {}
: GeneralMatrix(m, n)
{
}
SylvMatrix(const double *d, int m, int n)
: GeneralMatrix(d, m, n) {}
: GeneralMatrix(d, m, n)
{
}
SylvMatrix(double *d, int m, int n)
: GeneralMatrix(d, m, n) {}
: GeneralMatrix(d, m, n)
{
}
SylvMatrix(const GeneralMatrix &m)
: GeneralMatrix(m) {}
: GeneralMatrix(m)
{
}
SylvMatrix(const GeneralMatrix &m, int i, int j, int nrows, int ncols)
: GeneralMatrix(m, i, j, nrows, ncols) {}
: GeneralMatrix(m, i, j, nrows, ncols)
{
}
SylvMatrix(GeneralMatrix &m, int i, int j, int nrows, int ncols)
: GeneralMatrix(m, i, j, nrows, ncols) {}
: GeneralMatrix(m, i, j, nrows, ncols)
{
}
SylvMatrix(const GeneralMatrix &a, const GeneralMatrix &b)
: GeneralMatrix(a, b) {}
: GeneralMatrix(a, b)
{
}
/* this = |I 0|* this
|0 m| */
@ -47,20 +62,35 @@ public:
void eliminateRight(int row, int col, Vector &x);
};
class SqSylvMatrix : public SylvMatrix {
class SqSylvMatrix : public SylvMatrix
{
public:
SqSylvMatrix(int m) : SylvMatrix(m, m) {}
SqSylvMatrix(const double* d, int m) : SylvMatrix(d, m, m) {}
SqSylvMatrix(double* d, int m) : SylvMatrix(d, m, m) {}
SqSylvMatrix(const SqSylvMatrix& m) : SylvMatrix(m) {}
SqSylvMatrix(int m) : SylvMatrix(m, m)
{
}
SqSylvMatrix(const double *d, int m) : SylvMatrix(d, m, m)
{
}
SqSylvMatrix(double *d, int m) : SylvMatrix(d, m, m)
{
}
SqSylvMatrix(const SqSylvMatrix &m) : SylvMatrix(m)
{
}
SqSylvMatrix(const GeneralMatrix &m, int i, int j, int nrows)
: SylvMatrix(m, i, j, nrows, nrows) {}
: SylvMatrix(m, i, j, nrows, nrows)
{
}
SqSylvMatrix(GeneralMatrix &m, int i, int j, int nrows)
: SylvMatrix(m, i, j, nrows, nrows) {}
: SylvMatrix(m, i, j, nrows, nrows)
{
}
SqSylvMatrix(const GeneralMatrix &a, const GeneralMatrix &b);
const SqSylvMatrix& operator=(const SqSylvMatrix& m)
{GeneralMatrix::operator=(m); return *this;}
const SqSylvMatrix &
operator=(const SqSylvMatrix &m)
{
GeneralMatrix::operator=(m); return *this;
}
/* x = (this \otimes this..\otimes this)*d */
void multVecKron(KronVector &x, const KronVector &d) const;
/* x = (this' \otimes this'..\otimes this')*d */
@ -72,10 +102,8 @@ public:
void setUnit();
};
#endif /* SYLV_MATRIX_H */
// Local Variables:
// mode:C++
// End:

View File

@ -9,7 +9,8 @@
#include <new>
class MallocAllocator {
class MallocAllocator
{
#ifdef USE_MEMORY_POOL
public:
void *operator new(size_t size);
@ -26,7 +27,8 @@ void operator delete(void* p);
void operator delete[](void *p);
#endif
class SylvMemoryPool {
class SylvMemoryPool
{
char *base;
size_t length;
size_t allocated;
@ -43,7 +45,8 @@ public:
void setStackMode(bool);
};
class SylvMemoryDriver {
class SylvMemoryDriver
{
SylvMemoryDriver(const SylvMemoryDriver &);
const SylvMemoryDriver &operator=(const SylvMemoryDriver &);
public:
@ -57,7 +60,6 @@ protected:
#endif /* SYLV_MEMORY_H */
// Local Variables:
// mode:C++
// End:

View File

@ -15,27 +15,47 @@
typedef enum {def, changed, undef} status;
template <class _Type>
struct ParamItem {
struct ParamItem
{
protected:
typedef ParamItem<_Type> _Self;
status s;
_Type value;
public:
ParamItem()
{s = undef;}
{
s = undef;
}
ParamItem(_Type val)
{value = val; s = def;}
{
value = val; s = def;
}
ParamItem(const _Self &item)
{value = item.value; s = item.s;}
const _Self& operator=(const _Self& item)
{value = item.value; s = item.s; return *this;}
const _Self& operator=(const _Type& val)
{value = val; s = changed; return *this;}
_Type operator*() const
{return value;}
status getStatus() const
{return s;}
void print(FILE* f, const char* prefix, const char* str, const char* fmt) const
{
value = item.value; s = item.s;
}
const _Self &
operator=(const _Self &item)
{
value = item.value; s = item.s; return *this;
}
const _Self &
operator=(const _Type &val)
{
value = val; s = changed; return *this;
}
_Type
operator*() const
{
return value;
}
status
getStatus() const
{
return s;
}
void
print(FILE *f, const char *prefix, const char *str, const char *fmt) const
{
if (s == undef)
return;
@ -51,54 +71,95 @@ public:
}
};
class SylvParams {
class SylvParams
{
public:
typedef enum {iter, recurse} solve_method;
protected:
class DoubleParamItem : public ParamItem<double> {
class DoubleParamItem : public ParamItem<double>
{
public:
DoubleParamItem() : ParamItem<double>() {}
DoubleParamItem(double val) : ParamItem<double>(val) {}
DoubleParamItem(const DoubleParamItem& item) : ParamItem<double>(item) {}
const DoubleParamItem& operator=(const double& val)
{ParamItem<double>::operator=(val); return *this;}
DoubleParamItem() : ParamItem<double>()
{
}
DoubleParamItem(double val) : ParamItem<double>(val)
{
}
DoubleParamItem(const DoubleParamItem &item) : ParamItem<double>(item)
{
}
const DoubleParamItem &
operator=(const double &val)
{
ParamItem<double>::operator=(val); return *this;
}
#if defined(MATLAB_MEX_FILE) || defined(OCTAVE_MEX_FILE)
mxArray *createMatlabArray() const;
#endif
};
class IntParamItem : public ParamItem<int> {
class IntParamItem : public ParamItem<int>
{
public:
IntParamItem() : ParamItem<int>() {}
IntParamItem(int val) : ParamItem<int>(val) {}
IntParamItem(const IntParamItem& item) : ParamItem<int>(item) {}
const IntParamItem& operator=(const int& val)
{ParamItem<int>::operator=(val); return *this;}
IntParamItem() : ParamItem<int>()
{
}
IntParamItem(int val) : ParamItem<int>(val)
{
}
IntParamItem(const IntParamItem &item) : ParamItem<int>(item)
{
}
const IntParamItem &
operator=(const int &val)
{
ParamItem<int>::operator=(val); return *this;
}
#if defined(MATLAB_MEX_FILE) || defined(OCTAVE_MEX_FILE)
mxArray *createMatlabArray() const;
#endif
};
class BoolParamItem : public ParamItem<bool> {
class BoolParamItem : public ParamItem<bool>
{
public:
BoolParamItem() : ParamItem<bool>() {}
BoolParamItem(bool val) : ParamItem<bool>(val) {}
BoolParamItem(const BoolParamItem& item) : ParamItem<bool>(item) {}
const BoolParamItem& operator=(const bool& val)
{ParamItem<bool>::operator=(val); return *this;}
BoolParamItem() : ParamItem<bool>()
{
}
BoolParamItem(bool val) : ParamItem<bool>(val)
{
}
BoolParamItem(const BoolParamItem &item) : ParamItem<bool>(item)
{
}
const BoolParamItem &
operator=(const bool &val)
{
ParamItem<bool>::operator=(val); return *this;
}
#if defined(MATLAB_MEX_FILE) || defined(OCTAVE_MEX_FILE)
mxArray *createMatlabArray() const;
#endif
};
class MethodParamItem : public ParamItem<solve_method> {
class MethodParamItem : public ParamItem<solve_method>
{
public:
MethodParamItem() : ParamItem<solve_method>() {}
MethodParamItem(solve_method val) : ParamItem<solve_method>(val) {}
MethodParamItem(const MethodParamItem& item) : ParamItem<solve_method>(item) {}
const MethodParamItem operator=(const solve_method& val)
{ParamItem<solve_method>::operator=(val); return *this;}
MethodParamItem() : ParamItem<solve_method>()
{
}
MethodParamItem(solve_method val) : ParamItem<solve_method>(val)
{
}
MethodParamItem(const MethodParamItem &item) : ParamItem<solve_method>(item)
{
}
const MethodParamItem
operator=(const solve_method &val)
{
ParamItem<solve_method>::operator=(val); return *this;
}
#if defined(MATLAB_MEX_FILE) || defined(OCTAVE_MEX_FILE)
mxArray *createMatlabArray() const;
#endif
@ -138,12 +199,21 @@ public:
SylvParams(bool wc = false)
: method(recurse), convergence_tol(1.e-30), max_num_iter(15),
bs_norm(1.3), want_check(wc) {}
bs_norm(1.3), want_check(wc)
{
}
SylvParams(const SylvParams &p)
{copy(p);}
const SylvParams& operator=(const SylvParams& p)
{copy(p); return *this;}
~SylvParams() {}
{
copy(p);
}
const SylvParams &
operator=(const SylvParams &p)
{
copy(p); return *this;
}
~SylvParams()
{
}
void print(const char *prefix) const;
void print(FILE *fdesc, const char *prefix) const;
void setArrayNames(int &num, const char **names) const;
@ -156,7 +226,6 @@ private:
#endif /* SYLV_PARAMS_H */
// Local Variables:
// mode:C++
// End:

View File

@ -12,40 +12,47 @@
#include "SylvParams.h"
#include "SchurDecomp.h"
class SylvesterSolver {
class SylvesterSolver
{
protected:
const QuasiTriangular *const matrixK;
const QuasiTriangular *const matrixF;
private:
/* return true when it is more efficient to use QuasiTriangular
* than QuasiTriangularZero */
static bool zeroPad(const SchurDecompZero& kdecomp) {
return ((kdecomp.getZeroCols()*3 < kdecomp.getDim()*2) ||
(kdecomp.getZeroCols() < 10));
static bool
zeroPad(const SchurDecompZero &kdecomp)
{
return ((kdecomp.getZeroCols()*3 < kdecomp.getDim()*2)
|| (kdecomp.getZeroCols() < 10));
}
public:
SylvesterSolver(const QuasiTriangular &k, const QuasiTriangular &f)
: matrixK(new QuasiTriangular(k)),
matrixF(new QuasiTriangular(f))
{}
{
}
SylvesterSolver(const SchurDecompZero &kdecomp, const SchurDecomp &fdecomp)
: matrixK((zeroPad(kdecomp)) ?
new QuasiTriangular(kdecomp) : new QuasiTriangularZero(kdecomp)),
matrixF(new QuasiTriangular(fdecomp))
{}
{
}
SylvesterSolver(const SchurDecompZero &kdecomp, const SimilarityDecomp &fdecomp)
: matrixK((zeroPad(kdecomp)) ?
new QuasiTriangular(kdecomp) : new QuasiTriangularZero(kdecomp)),
matrixF(new BlockDiagonal(fdecomp.getB()))
{}
{
}
virtual ~SylvesterSolver()
{delete matrixK; delete matrixF;}
{
delete matrixK; delete matrixF;
}
virtual void solve(SylvParams &pars, KronVector &x) const = 0;
};
#endif /* SYLVESTER_SOLVER_H */
// Local Variables:
// mode:C++
// End:

View File

@ -7,7 +7,8 @@
#include "SylvMatrix.h"
class SymSchurDecomp {
class SymSchurDecomp
{
protected:
Vector lambda;
SqSylvMatrix q;
@ -16,12 +17,22 @@ public:
* symmetric and Lambda real diagonal, hence a vector. */
SymSchurDecomp(const GeneralMatrix &a);
SymSchurDecomp(const SymSchurDecomp &ssd)
: lambda(ssd.lambda), q(ssd.q) {}
virtual ~SymSchurDecomp() {}
const Vector& getLambda() const
{return lambda;}
const SqSylvMatrix& getQ() const
{return q;}
: lambda(ssd.lambda), q(ssd.q)
{
}
virtual ~SymSchurDecomp()
{
}
const Vector &
getLambda() const
{
return lambda;
}
const SqSylvMatrix &
getQ() const
{
return q;
}
/** Return factor F*F^T = A, raises and exception if A is not
* positive semidefinite, F must be square. */
void getFactor(GeneralMatrix &f) const;
@ -35,7 +46,6 @@ public:
#endif
// Local Variables:
// mode:C++
// End:

View File

@ -11,14 +11,16 @@
#include "QuasiTriangularZero.h"
#include "SimilarityDecomp.h"
class TriangularSylvester : public SylvesterSolver {
class TriangularSylvester : public SylvesterSolver
{
const QuasiTriangular *const matrixKK;
const QuasiTriangular *const matrixFF;
public:
TriangularSylvester(const QuasiTriangular &k, const QuasiTriangular &f);
TriangularSylvester(const SchurDecompZero &kdecomp, const SchurDecomp &fdecomp);
TriangularSylvester(const SchurDecompZero &kdecomp, const SimilarityDecomp &fdecomp);
virtual ~TriangularSylvester();
virtual
~TriangularSylvester();
void print() const;
void solve(SylvParams &pars, KronVector &d) const;
@ -36,11 +38,14 @@ public:
void linEval(double alpha, double beta1, double beta2,
KronVector &x1, KronVector &x2,
const ConstKronVector &d1, const ConstKronVector &d2) const;
void linEval(double alpha, double beta1, double beta2,
void
linEval(double alpha, double beta1, double beta2,
KronVector &x1, KronVector &x2,
const KronVector &d1, const KronVector &d2) const
{linEval(alpha, beta1, beta2, x1, x2,
ConstKronVector(d1), ConstKronVector(d2));}
{
linEval(alpha, beta1, beta2, x1, x2,
ConstKronVector(d1), ConstKronVector(d2));
}
/* evaluates:
|x1| |d1| |gamma -delta1| |d1|
@ -55,12 +60,15 @@ public:
double gamma, double delta1, double delta2,
KronVector &x1, KronVector &x2,
const ConstKronVector &d1, const ConstKronVector &d2) const;
void quaEval(double alpha, double betas,
void
quaEval(double alpha, double betas,
double gamma, double delta1, double delta2,
KronVector &x1, KronVector &x2,
const KronVector &d1, const KronVector &d2) const
{quaEval(alpha, betas, gamma, delta1, delta2, x1, x2,
ConstKronVector(d1), ConstKronVector(d2));}
{
quaEval(alpha, betas, gamma, delta1, delta2, x1, x2,
ConstKronVector(d1), ConstKronVector(d2));
}
private:
/* returns square of size of minimal eigenvalue of the system solved,
now obsolete */
@ -109,7 +117,6 @@ private:
#endif /* TRIANGULAR_SYLVESTER_H */
// Local Variables:
// mode:C++
// End:

View File

@ -14,43 +14,74 @@
class GeneralMatrix;
class ConstVector;
class Vector {
class Vector
{
protected:
int len;
int s;
double *data;
bool destroy;
public:
Vector() : len(0), s(1), data(0), destroy(false) {}
Vector(int l) : len(l), s(1), data(new double[l]), destroy(true) {}
Vector(Vector& v) : len(v.length()), s(v.skip()), data(v.base()), destroy(false) {}
Vector() : len(0), s(1), data(0), destroy(false)
{
}
Vector(int l) : len(l), s(1), data(new double[l]), destroy(true)
{
}
Vector(Vector &v) : len(v.length()), s(v.skip()), data(v.base()), destroy(false)
{
}
Vector(const Vector &v);
Vector(const ConstVector &v);
Vector(const double *d, int l)
: len(l), s(1), data(new double[len]), destroy(true)
{copy(d, 1);}
{
copy(d, 1);
}
Vector(double *d, int l)
: len(l), s(1), data(d), destroy(false) {}
: len(l), s(1), data(d), destroy(false)
{
}
Vector(double *d, int skip, int l)
: len(l), s(skip), data(d), destroy(false) {}
: len(l), s(skip), data(d), destroy(false)
{
}
Vector(Vector &v, int off, int l);
Vector(const Vector &v, int off, int l);
Vector(GeneralMatrix &m, int col);
Vector(int row, GeneralMatrix &m);
const Vector &operator=(const Vector &v);
const Vector &operator=(const ConstVector &v);
double& operator[](int i)
{return data[s*i];}
const double& operator[](int i) const
{return data[s*i];}
const double* base() const
{return data;}
double* base()
{return data;}
int length() const
{return len;}
int skip() const
{return s;}
double &
operator[](int i)
{
return data[s*i];
}
const double &
operator[](int i) const
{
return data[s*i];
}
const double *
base() const
{
return data;
}
double *
base()
{
return data;
}
int
length() const
{
return len;
}
int
skip() const
{
return s;
}
/** Exact equality. */
bool operator==(const Vector &y) const;
@ -61,11 +92,16 @@ public:
bool operator>(const Vector &y) const;
bool operator>=(const Vector &y) const;
virtual ~Vector();
virtual
~Vector();
void zeros();
void nans();
void infs();
bool toBeDestroyed() const {return destroy;}
bool
toBeDestroyed() const
{
return destroy;
}
void rotatePair(double alpha, double beta1, double beta2, int i);
void add(double r, const Vector &v);
void add(double r, const ConstVector &v);
@ -91,63 +127,106 @@ public:
Vector &x1, Vector &x2,
const Vector &b1, const Vector &b2);
/* same, but subtracts instead of add */
static void mult2s(double alpha, double beta1, double beta2,
static void
mult2s(double alpha, double beta1, double beta2,
Vector &x1, Vector &x2,
const Vector &b1, const Vector &b2)
{mult2a(-alpha, -beta1, -beta2, x1, x2, b1, b2);}
{
mult2a(-alpha, -beta1, -beta2, x1, x2, b1, b2);
}
private:
void copy(const double *d, int inc);
const Vector &operator=(int); // must not be used (not implemented)
const Vector &operator=(double); // must not be used (not implemented)
};
class BaseConstVector {
class BaseConstVector
{
protected:
int len;
int s;
const double *data;
public:
BaseConstVector(int l, int si, const double* d) : len(l), s(si), data(d) {}
BaseConstVector(const BaseConstVector& v) : len(v.len), s(v.s), data(v.data) {}
const BaseConstVector& operator=(const BaseConstVector& v)
{len = v.len; s = v.s; data = v.data; return *this;}
const double& operator[](int i) const
{return data[s*i];}
const double* base() const
{return data;}
int length() const
{return len;}
int skip() const
{return s;}
BaseConstVector(int l, int si, const double *d) : len(l), s(si), data(d)
{
}
BaseConstVector(const BaseConstVector &v) : len(v.len), s(v.s), data(v.data)
{
}
const BaseConstVector &
operator=(const BaseConstVector &v)
{
len = v.len; s = v.s; data = v.data; return *this;
}
const double &
operator[](int i) const
{
return data[s*i];
}
const double *
base() const
{
return data;
}
int
length() const
{
return len;
}
int
skip() const
{
return s;
}
};
class ConstGeneralMatrix;
class ConstVector : public BaseConstVector {
class ConstVector : public BaseConstVector
{
public:
ConstVector(const Vector& v) : BaseConstVector(v.length(), v.skip(), v.base()) {}
ConstVector(const ConstVector& v) : BaseConstVector(v) {}
ConstVector(const double* d, int l) : BaseConstVector(l, 1, d) {}
ConstVector(const Vector &v) : BaseConstVector(v.length(), v.skip(), v.base())
{
}
ConstVector(const ConstVector &v) : BaseConstVector(v)
{
}
ConstVector(const double *d, int l) : BaseConstVector(l, 1, d)
{
}
ConstVector(const Vector &v, int off, int l);
ConstVector(const ConstVector &v, int off, int l);
ConstVector(const double *d, int skip, int l);
ConstVector(const ConstGeneralMatrix &m, int col);
ConstVector(int row, const ConstGeneralMatrix &m);
virtual ~ConstVector() {}
virtual ~ConstVector()
{
}
/** Exact equality. */
bool operator==(const ConstVector &y) const;
bool operator!=(const ConstVector& y) const
{return ! operator==(y);}
bool
operator!=(const ConstVector &y) const
{
return !operator==(y);
}
/** Lexicographic ordering. */
bool operator<(const ConstVector &y) const;
bool operator<=(const ConstVector& y) const
{return operator<(y) || operator==(y);}
bool operator>(const ConstVector& y) const
{return ! operator<=(y);}
bool operator>=(const ConstVector& y) const
{return ! operator<(y);}
bool
operator<=(const ConstVector &y) const
{
return operator<(y) || operator==(y);
}
bool
operator>(const ConstVector &y) const
{
return !operator<=(y);
}
bool
operator>=(const ConstVector &y) const
{
return !operator<(y);
}
double getNorm() const;
double getMax() const;
@ -157,19 +236,23 @@ public:
void print() const;
};
class ZeroPad {
class ZeroPad
{
public:
static const int length = 16;
private:
double pad[16];
public:
ZeroPad();
const double* getBase() const {return pad;}
const double *
getBase() const
{
return pad;
}
};
#endif /* VECTOR_H */
// Local Variables:
// mode:C++
// End:

View File

@ -12,28 +12,55 @@
using namespace std;
class MMException : public MallocAllocator {
class MMException : public MallocAllocator
{
string message;
public:
MMException(string mes) : message(mes) {}
MMException(const char* mes) : message(mes) {}
const char* getMessage() const {return message.data();}
MMException(string mes) : message(mes)
{
}
MMException(const char *mes) : message(mes)
{
}
const char *
getMessage() const
{
return message.data();
}
};
class MMMatrixIn : public MallocAllocator {
class MMMatrixIn : public MallocAllocator
{
double *data;
int rows;
int cols;
public:
MMMatrixIn(const char *fname);
~MMMatrixIn();
const double* getData() const {return data;}
int size() const {return rows*cols;}
int row() const {return rows;}
int col() const {return cols;}
const double *
getData() const
{
return data;
}
int
size() const
{
return rows*cols;
}
int
row() const
{
return rows;
}
int
col() const
{
return cols;
}
};
class MMMatrixOut : public MallocAllocator {
class MMMatrixOut : public MallocAllocator
{
public:
static void write(const char *fname, int rows, int cols, const double *data);
static void write(const char *fname, const GeneralMatrix &m);

View File

@ -11,7 +11,8 @@
#include "rfs_tensor.h"
#include "t_container.h"
class Factory {
class Factory
{
void init(const Symmetry &s, const IntSequence &nvs);
void init(int dim, int nv);
void fillMatrix(TwoDMatrix &m) const;
@ -19,7 +20,8 @@ public:
double get() const;
// this can be used with UGSTensor, FGSTensor
template <class _Ttype>
_Ttype* make(int r, const Symmetry& s, const IntSequence& nvs)
_Ttype *
make(int r, const Symmetry &s, const IntSequence &nvs)
{
_Ttype *res = new _Ttype(r, TensorDimens(s, nvs));
init(s, nvs);
@ -29,7 +31,8 @@ public:
// this can be used with FFSTensor, UFSTensor, FRTensor, URTensor
template <class _Ttype>
_Ttype* make(int r, int nv, int dim)
_Ttype *
make(int r, int nv, int dim)
{
_Ttype *res = new _Ttype(r, nv, dim);
init(dim, nv);
@ -38,19 +41,25 @@ public:
}
template <class _Ttype, class _Ctype>
_Ctype* makeCont(int r, const IntSequence& nvs, int maxdim)
_Ctype *
makeCont(int r, const IntSequence &nvs, int maxdim)
{
int symnum = nvs.size();
_Ctype *res = new _Ctype(symnum);
for (int dim = 1; dim <= maxdim; dim++) {
if (symnum == 1) {
for (int dim = 1; dim <= maxdim; dim++)
{
if (symnum == 1)
{
// full symmetry
Symmetry sym(dim);
_Ttype *t = make<_Ttype>(r, sym, nvs);
res->insert(t);
} else {
}
else
{
// general symmetry
for (int i = 0; i <= dim; i++) {
for (int i = 0; i <= dim; i++)
{
Symmetry sym(i, dim-i);
_Ttype *t = make<_Ttype>(r, sym, nvs);
res->insert(t);
@ -61,10 +70,12 @@ public:
}
template <class _Ttype, class _Ptype>
_Ptype* makePoly(int r, int nv, int maxdim)
_Ptype *
makePoly(int r, int nv, int maxdim)
{
_Ptype *p = new _Ptype(r, nv);
for (int d = 1; d <= maxdim; d++) {
for (int d = 1; d <= maxdim; d++)
{
_Ttype *t = make<_Ttype>(r, nv, d);
p->insert(t);
}

View File

@ -10,20 +10,23 @@
#include "sparse_tensor.h"
#include "Vector.h"
class IntGenerator {
class IntGenerator
{
int maxim;
double probab;
public:
IntGenerator()
: maxim(5), probab(0.3) {}
: maxim(5), probab(0.3)
{
}
void init(int nf, int ny, int nv, int nw, int nu, int mx, double prob);
int get() const;
};
extern IntGenerator intgen;
class Monom : public IntSequence {
class Monom : public IntSequence
{
public:
Monom(int len); // generate a random monom
Monom(int len, int item); // generate monom whose items are the given item
@ -34,7 +37,8 @@ public:
};
class Monom2Vector;
class Monom1Vector {
class Monom1Vector
{
friend class Monom2Vector;
int nx;
int len;
@ -48,7 +52,8 @@ public:
};
//class Monom3Vector;
class Monom2Vector {
class Monom2Vector
{
int ny;
int nu;
int len;
@ -66,7 +71,8 @@ public:
void print() const;
};
class Monom4Vector {
class Monom4Vector
{
int len;
int nx1;
int nx2;
@ -95,8 +101,8 @@ protected:
void init_random();
};
struct SparseDerivGenerator {
struct SparseDerivGenerator
{
int maxdimen;
FGSContainer *bigg;
FGSContainer *g;
@ -107,8 +113,8 @@ struct SparseDerivGenerator {
~SparseDerivGenerator();
};
struct DenseDerivGenerator {
struct DenseDerivGenerator
{
int maxdimen;
FGSContainer *xcont;
FGSContainer *rcont;

View File

@ -11,10 +11,12 @@
#include <string>
#include <algorithm>
namespace ogu {
namespace ogu
{
/** A primitive exception. */
class Exception {
class Exception
{
static const int file_length = 100;
static const int mes_length = 500;
protected:
@ -38,13 +40,24 @@ namespace ogu {
strncpy(mes, m.c_str(), std::min(mes_length-1, (int) m.length()));
mes[mes_length-1] = '\0';
}
virtual ~Exception() {}
void print(FILE* fd) const
{ fprintf(fd, "%s:%d: %s\n", file, line, mes); }
void print() const
{ print(stdout); }
const char* message() const
{ return mes; }
virtual ~Exception()
{
}
void
print(FILE *fd) const
{
fprintf(fd, "%s:%d: %s\n", file, line, mes);
}
void
print() const
{
print(stdout);
}
const char *
message() const
{
return mes;
}
};
};

View File

@ -5,7 +5,8 @@
#ifndef OGU_MEMORY_FILE
#define OGU_MEMORY_FILE
namespace ogu {
namespace ogu
{
/** This function calculates an offset of a given position in a
* given string. The position is given by the line number and by
* the offset in the line (both starting from 1). */
@ -22,34 +23,51 @@ namespace ogu {
* could be opened for reading, data is NULL and length is -1. If
* the file is empty but exists, len is zero and data points to a
* newly allocated memory containing '\0' character at the end. */
class MemoryFile {
class MemoryFile
{
protected:
int len;
char *data;
public:
MemoryFile(const char *fname);
virtual ~MemoryFile()
{if (data) delete [] data;}
int length() const
{return len;}
const char* base() const
{return data;}
bool exists() const
{return len != -1;}
{
if (data)
delete [] data;
}
int
length() const
{
return len;
}
const char *
base() const
{
return data;
}
bool
exists() const
{
return len != -1;
}
/** Return the offset of a character in the given line
* (starting from 1) with the given offset in the line. */
int offset(int line, int lineoff) const
{return calc_pos_offset(len, data, line, lineoff);}
int
offset(int line, int lineoff) const
{
return calc_pos_offset(len, data, line, lineoff);
}
/** Return the line number and column number of the character
* defined by the offset. */
void line_and_col(int offset, int& line, int& col) const
{calc_pos_line_and_col(len, data, offset, line, col);}
void
line_and_col(int offset, int &line, int &col) const
{
calc_pos_line_and_col(len, data, offset, line, col);
}
};
};
#endif
// Local Variables:

View File

@ -7,31 +7,43 @@
#include <vector>
namespace ogu {
namespace ogu
{
using std::vector;
class PascalRow : public vector<int> {
class PascalRow : public vector<int>
{
int k;
public:
PascalRow()
: vector<int>(), k(1)
{ push_back(2); }
{
push_back(2);
}
void setFromPrevious(const PascalRow &prev);
void prolong(const PascalRow &prev);
void prolongFirst(int n);
void print() const;
};
class PascalTriangle {
class PascalTriangle
{
vector<PascalRow> tr;
public:
PascalTriangle()
{tr.push_back(PascalRow());}
{
tr.push_back(PascalRow());
}
PascalTriangle(const PascalTriangle &triang)
: tr(triang.tr) {}
const PascalTriangle& operator=(const PascalTriangle& triang)
{ tr = triang.tr; return *this;}
: tr(triang.tr)
{
}
const PascalTriangle &
operator=(const PascalTriangle &triang)
{
tr = triang.tr; return *this;
}
int noverk(int n, int k);
void print() const;
protected:
@ -43,7 +55,6 @@ namespace ogu {
extern ogu::PascalTriangle ptriang;
#endif
// Local Variables:

View File

@ -155,18 +155,18 @@ sigma_m =-5.85;
Lambdamu=3.4e-3;
LambdaA = 2.8e-3;
LambdaYd= (LambdaA+alppha*Lambdamu)/(1-alppha);
/*
The following parameters are set in the steady state file as they depend on other
deep parameters that were estimated in the original study. Setting them in the
deep parameters (some were estimated in the original study). Setting them in the
steady state file means they are updated for every parameter draw in the MCMC
algorithm, while the parameters initialized here are only set once for the initial
values of the parameters they depend on:
gammma1 as it depends on LambdaA, alppha, Lambdamu, betta, and delta
Rbar =0 as it depends on PI, LambdaA, alppha, Lambdamu, and betta
Lambdax
gammma1=mu_z*mu_I/betta-(1-delta);
R=1+(PIbar*mu_z/betta-1);
Lambdax=exp(LambdaYd);
LambdaYd= (LambdaA+alppha*Lambdamu)/(1-alppha);
*/

View File

@ -36,6 +36,7 @@ d=1;
phi=1;
m=0;
zeta=1;
LambdaYd= (LambdaA+alppha*Lambdamu)/(1-alppha);
mu_z=exp(LambdaYd);
mu_I=exp(Lambdamu);
mu_A=exp(LambdaA);

View File

@ -6,22 +6,26 @@
* The data are in file "fsdat_simul.m", and have been artificially generated.
* They are therefore different from the original dataset used by Schorfheide.
*
* The prior distribution follows the one originally specified in Schorfheide's paper.
* Note that the beta prior for rho implies an asymptote and corresponding prior mode
* for rho at 0. It is generally recommended to avoid this extreme type of prior.
* The prior distribution follows the one originally specified in Schorfheide's
* paper, except for parameter rho. In the paper, the elicited beta prior for rho
* implies an asymptote and corresponding prior mode at 0. It is generally
* recommended to avoid this extreme type of prior. Some optimizers, for instance
* mode_compute=12 (Mathworks' particleswarm algorithm) may find a posterior mode
* with rho equal to zero. We lowered the value of the prior standard deviation
* (changing .223 to .100) to remove the asymptote.
*
* The equations are taken from J. Nason and T. Cogley (1994): "Testing the
* implications of long-run neutrality for monetary business cycle models",
* Journal of Applied Econometrics, 9, S37-S70.
* Note that there is an initial minus sign missing in equation (A1), p. S63.
*
* This implementation was written by Michel Juillard. Please note that the
* This implementation was originally written by Michel Juillard. Please note that the
* following copyright notice only applies to this Dynare implementation of the
* model.
*/
/*
* Copyright (C) 2004-2015 Dynare Team
* Copyright (C) 2004-2017 Dynare Team
*
* This file is part of Dynare.
*
@ -109,7 +113,7 @@ alp, beta_pdf, 0.356, 0.02;
bet, beta_pdf, 0.993, 0.002;
gam, normal_pdf, 0.0085, 0.003;
mst, normal_pdf, 1.0002, 0.007;
rho, beta_pdf, 0.129, 0.223;
rho, beta_pdf, 0.129, 0.100;
psi, beta_pdf, 0.65, 0.05;
del, beta_pdf, 0.01, 0.005;
stderr e_a, inv_gamma_pdf, 0.035449, inf;

View File

@ -1,16 +1,15 @@
Format: http://www.debian.org/doc/packaging-manuals/copyright-format/1.0/
Format: https://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: *
@ -18,8 +17,8 @@ Copyright: 1996-2017 Dynare Team
License: GPL-3+
Files: matlab/AIM/SP*
Copyright: public-domain
License: public-domain
Copyright: none
License: public-domain-aim
This code is in the public domain and may be used freely.
However the authors would appreciate acknowledgement of the source by
citation of any of the following papers:
@ -86,6 +85,57 @@ Copyright: 2016 Benjamin Born and Johannes Pfeifer
2016 Dynare Team
License: GPL-3+
Files: matlab/gsa/Morris_Measure_Groups.m
matlab/gsa/Sampling_Function_2.m
Copyright: 2005 European Commission
2012 Dynare Team
License: GPL-3+
Comment: Written by Jessica Cariboni and Francesca Campolongo
Joint Research Centre, The European Commission,
Files: matlab/gsa/cumplot.m
matlab/gsa/filt_mc_.m
matlab/gsa/gsa_plotmatrix.m
matlab/gsa/gsa_skewness.m
matlab/gsa/gsa_speed.m
matlab/gsa/log_trans_.m
matlab/gsa/map_calibration.m
matlab/gsa/map_ident_.m
matlab/gsa/mcf_analysis.m
matlab/gsa/myboxplot.m
matlab/gsa/myprctilecol.m
matlab/gsa/prior_draw_gsa.m
matlab/gsa/read_data.m
matlab/gsa/redform_map.m
matlab/gsa/redform_screen.m
matlab/gsa/scatter_mcf.m
matlab/gsa/smirnov.m
matlab/gsa/stab_map_.m
matlab/gsa/stab_map_1.m
matlab/gsa/stab_map_2.m
matlab/gsa/stand_.m
matlab/gsa/tcrit.m
matlab/gsa/teff.m
matlab/gsa/trank.m
Copyright: 2011-2017 European Commission
2011-2017 Dynare Team
License: GPL-3+
Files: matlab/gsa/pick.m
Copyright: none
License: public-domain-jrc
This software has been developed at the Joint Research Centre of European Commission
by officers in the course of their official duties. This software is not subject to copyright
protection and is in the public domain. It is an experimental system. The Joint Research Centre
of European Commission assumes no responsibility whatsoever for its use by other parties
and makes no guarantees, expressed or implied, about its quality, reliability, or any other
characteristic. We would appreciate acknowledgement if the software is used.
Comment: This file is part of GLUEWIN.
The program has been developed by M. Ratto, European Commission, Joint Research Centre,
Institute for the Protection and Security of The Citizen, Technological and Economic Risk Management,
Applied Statistics, as a deliverable of the IMPACT project
(EC Fifth Framework Programme, SCA Project, IST-1999-11313, DG-INFSO).
Files: matlab/optimization/simpsa.m matlab/optimization/simpsaget.m matlab/optimization/simpsaset.m
Copyright: 2005 Henning Schmidt, FCC, henning@fcc.chalmers.se
2006 Brecht Donckels, BIOMATH, brecht.donckels@ugent.be
@ -127,37 +177,24 @@ Copyright: 2005 Jos van der Geest <jos@jasen.nl>
2013 Christophe Gouel
2016 Dynare Team
License: BSD-2-clause
Redistribution and use in source and binary forms, with or without
modification, are permitted provided that the following conditions are
met:
.
* Redistributions of source code must retain the above copyright
notice, this list of conditions and the following disclaimer.
* Redistributions in binary form must reproduce the above copyright
notice, this list of conditions and the following disclaimer in
the documentation and/or other materials provided with the distribution
.
THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE
LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
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/lmmcp/lmmcp.m
Copyright: 2005 Christian Kanzow and Stefania Petra
2013 Christophe Gouel
2014 Dynare Team
License: permissive
License: permissive-lmmcp
Unlimited permission is granted to everyone to use, copy, modify or
distribute this software.
Files: doc/dynare.texi doc/*.tex doc/*.svg doc/*.dia doc/*.pdf doc/*.bib
Files: matlab/utilities/graphics/distinguishable_colors.m
Copyright: 2010-2011 Timothy E. Holy
License: BSD-2-clause
Files: matlab/utilities/graphics/colorspace.m
Copyright: 2005-2010 Pascal Getreuer
License: BSD-2-clause
Files: doc/dynare.texi doc/*.tex doc/*.svg doc/*.pdf doc/*.bib
Copyright: 1996-2017 Dynare Team
License: GFDL-NIV-1.3+
@ -206,7 +243,7 @@ Files: m4/ax_compare_version.m4
Copyright: 2008 Tim Toolan <toolan@ele.uri.edu>
License: permissive-autoconf
Files: m4/ax_latex_bibtex_test.m4 m4/ax_latex_class.m4 m4/ax_tex_test.m4
Files: m4/ax_latex_class.m4 m4/ax_tex_test.m4
Copyright: 2008 Boretti Mathieu <boretti@eig.unige.ch>
2009 Dynare Team
License: LGPL-2.1+
@ -235,7 +272,7 @@ Copyright: 1996-2011 Daniel Waggoner and Tao Zha
License: GPL-3+
Files: contrib/ms-sbvar/switch_dw/state_space/sbvar/dw_csminwel.c
state_space/sbvar/dw_csminwel.h
contrib/ms-sbvar/switch_dw/state_space/sbvar/dw_csminwel.h
Copyright: 1996 Christopher Sims
2003 Karibzhanov, Waggoner and Zha
License: GPL-3+
@ -316,7 +353,7 @@ License: GPL-3+
Files: contrib/dmm/randlib/*
Copyright: none
License: public-domain
License: public-domain-dmm
We place the Randlib code that we have written in the public domain.
.
NO WARRANTY
@ -336,6 +373,29 @@ License: public-domain
ITS ANALYSIS BEING RENDERED INACCURATE OR LOSSES SUSTAINED BY THIRD
PARTIES) THE PROGRAM.
License: BSD-2-clause
Redistribution and use in source and binary forms, with or without
modification, are permitted provided that the following conditions are
met:
.
* Redistributions of source code must retain the above copyright
notice, this list of conditions and the following disclaimer.
* Redistributions in binary form must reproduce the above copyright
notice, this list of conditions and the following disclaimer in
the documentation and/or other materials provided with the distribution
.
THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE
LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
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.
License: GFDL-NIV-1.3+
Permission is granted to copy, distribute and/or modify this document
under the terms of the GNU Free Documentation License, Version 1.3 or
@ -344,21 +404,6 @@ License: GFDL-NIV-1.3+
.
A copy of the license can be found at <http://www.gnu.org/licenses/fdl.txt>
License: GPL-2+
This program is free software: you can redistribute it and/or
modify it under the terms of the GNU General Public License as
published by the Free Software Foundation, either version 2 of
the License, or (at your option) any later version.
.
This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.
.
You should have received a copy of the GNU General Public License
along with this program. If not, see
<http://www.gnu.org/licenses/>.
License: GPL-2+ with Autoconf exception
This program is free software; you can redistribute it and/or
modify it under the terms of the GNU General Public License as

View File

@ -1,66 +0,0 @@
# ===========================================================================
# http://www.nongnu.org/autoconf-archive/ax_latex_test.html
# ===========================================================================
#
# OBSOLETE MACRO
#
# Deprecated because of licensing issues. The Lesser GPL imposes licensing
# restrictions on the generated configure script unless it is augmented
# with an Autoconf Exception clause.
#
# SYNOPSIS
#
# AX_LATEX_BIBTEX_TEST(FILEDATA,BIBDATA,VARIABLETOSET,[NOCLEAN])
#
# DESCRIPTION
#
# This macros creates a bib file called contest.bib with BIBDATA,
# executes the latex application with FILEDATA as input, then runs
# bibtex on the resulting aux file, and finally sets VARIABLETOSET
# to yes or no depending on the result. If NOCLEAN is set, the folder
# used for the test is not deleted after testing.
#
# The macro assumes that the variables PDFLATEX and BIBTEX are set.
#
# Adapted from the macro AX_LATEX_TEST by Sébastien Villemot.
#
# LICENSE
#
# Copyright (c) 2008 Boretti Mathieu <boretti@eig.unige.ch>
# Copyright (c) 2009 Dynare Team
#
# This library is free software; you can redistribute it and/or modify it
# under the terms of the GNU Lesser General Public License as published by
# the Free Software Foundation; either version 2.1 of the License, or (at
# your option) any later version.
#
# This library is distributed in the hope that it will be useful, but
# WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser
# General Public License for more details.
#
# You should have received a copy of the GNU Lesser General Public License
# along with this library. If not, see <http://www.gnu.org/licenses/>.
AC_DEFUN([AX_LATEX_BIBTEX_TEST],[
rm -rf conftest.dir/.acltx
AS_MKDIR_P([conftest.dir/.acltx])
cd conftest.dir/.acltx
m4_ifval([$3],[$3="no"; export $3;])
cat > conftest.tex << ACLEOF
$1
ACLEOF
cat > conftest.bib << ACLEOF
$2
ACLEOF
$PDFLATEX conftest 2>&1 1>output
$BIBTEX conftest 2>&1 1>output2 m4_ifval([$3],[&& $3=yes])
cd ..
cd ..
sed 's/^/| /' conftest.dir/.acltx/conftest.tex >&5
echo "$as_me:$LINENO: executing $PDFLATEX conftest" >&5
sed 's/^/| /' conftest.dir/.acltx/output >&5
echo "$as_me:$LINENO: executing $BIBTEX conftest" >&5
sed 's/^/| /' conftest.dir/.acltx/output2 >&5
m4_ifval([$4],,[rm -rf conftest.dir/.acltx])
])

View File

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

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,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 @@ function [AHess, DLIK, LIK] = AHessian(T,R,Q,H,P,Y,DT,DYss,DOm,DH,DP,start,mf,ka
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
@ -127,7 +127,7 @@ end
end
AHess = -AHess;
if nargout > 1,
if nargout > 1
DLIK = DLIK/2;
end
% adding log-likelihhod constants
@ -150,6 +150,3 @@ for ii = 1:k
end
% 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

@ -66,12 +66,12 @@ 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)
@ -81,7 +81,7 @@ if (ia ~= 0)
return
end
[w,rts,lgroots,flag_trouble]=SPEigensystem(a,uprbnd,min(length(js),qrows-iq+1));
if flag_trouble==1;
if flag_trouble==1
disp('Problem in SPEIG');
aimcode=64;
return

View File

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

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,12 +3,12 @@ function [dr,aimcode,rts]=dynAIMsolver1(jacobia_,M_,dr)
% 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, . . . ,?
% 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
% 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'= £
% where [fy'-$... fy'i ... fy'+&]=[H-$... Hi ... H+&] and fu'= £
%
% INPUTS
% jacobia_ [matrix] 1st order derivative of the model
@ -46,7 +46,7 @@ function [dr,aimcode,rts]=dynAIMsolver1(jacobia_,M_,dr)
%
% GP July 2008
% Copyright (C) 2008-2012 Dynare Team
% Copyright (C) 2008-2017 Dynare Team
%
% This file is part of Dynare.
%

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 @@ function [dr,info]=AIM_first_order_solver(jacobia,M,dr,qz_criterium)
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

@ -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.
%

View File

@ -58,7 +58,7 @@ function [alphahat,etahat,epsilonhat,ahat,SteadyState,trend_coeff,aK,T,R,P,PK,de
% SPECIAL REQUIREMENTS
% None
% Copyright (C) 2006-2016 Dynare Team
% Copyright (C) 2006-2017 Dynare Team
%
% This file is part of Dynare.
%
@ -255,8 +255,10 @@ if kalman_algo == 2 || kalman_algo == 4
[Pstar,Pinf] = compute_Pinf_Pstar(mf,ST,R1,Q,options_.qz_criterium);
else
Pstar = blkdiag(Pstar,H);
if ~isempty(Pinf)
Pinf = blkdiag(Pinf,zeros(vobs));
end
end
%now reset H to 0
H = zeros(vobs,vobs);
else

View File

@ -16,7 +16,7 @@ function Draws = GetAllPosteriorDraws(column, FirstMhFile, FirstLine, TotalNumbe
% SPECIAL REQUIREMENTS
% none
% Copyright (C) 2005-2011 Dynare Team
% Copyright (C) 2005-2017 Dynare Team
%
% This file is part of Dynare.
%

View File

@ -13,7 +13,7 @@ function [xparams, logpost] = GetOneDraw(type)
% 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

@ -16,7 +16,7 @@ function oo_ = GetPosteriorParametersStatistics(estim_params_, M_, options_, bay
% SPECIAL REQUIREMENTS
% None.
% Copyright (C) 2006-2016 Dynare Team
% Copyright (C) 2006-2017 Dynare Team
%
% This file is part of Dynare.
%
@ -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);
@ -311,7 +311,7 @@ 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);

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.
%

View File

@ -16,7 +16,7 @@ function oo_ = PlotPosteriorDistributions(estim_params_, M_, options_, bayestopt
% SPECIAL REQUIREMENTS
% none
% Copyright (C) 2005-2016 Dynare Team
% Copyright (C) 2005-2017 Dynare Team
%
% This file is part of Dynare.
%
@ -63,7 +63,7 @@ for i=1:npar
subplotnum = subplotnum+1;
if subplotnum == 1
figunumber = figunumber+1;
hfig=dyn_figure(options_,'Name',figurename);
hfig=dyn_figure(options_.nodisplay,'Name',figurename);
end
[nam,texnam] = get_the_name(i,TeX,M_,estim_params_,options_);
if subplotnum == 1
@ -151,8 +151,8 @@ for i=1:npar
title(nam,'Interpreter','none');
hold off;
drawnow
if subplotnum == MaxNumberOfPlotPerFigure || i == npar;
dyn_saveas(hfig,[OutputDirectoryName '/' M_.fname '_PriorsAndPosteriors' int2str(figunumber)],options_);
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)

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;
@ -235,12 +235,17 @@ else
% which files have to be copied to run remotely
NamFileInput(1,:) = {'',[M_.fname '_static.m']};
NamFileInput(2,:) = {'',[M_.fname '_dynamic.m']};
if options_.steadystate_flag,
NamFileInput(3,:) = {'',[M_.fname '_set_auxiliary_variables.m']};
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']};
end
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
@ -438,11 +443,11 @@ end
% 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

@ -30,7 +30,7 @@ function myoutput=PosteriorIRF_core2(myinputs,fpar,npar,whoiam,ThisMatlab)
% SPECIAL REQUIREMENTS.
% None.
%
% Copyright (C) 2006-2016 Dynare Team
% Copyright (C) 2006-2017 Dynare Team
%
% This file is part of Dynare.
%
@ -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,16 +96,16 @@ end
OutputFileName={};
subplotnum = 0;
for i=fpar:npar,
for i=fpar:npar
figunumber = 0;
for j=1:nvar
if max(abs(MeanIRF(:,j,i))) >= options_.impulse_responses.plot_threshold
subplotnum = subplotnum+1;
if subplotnum == 1 && options_.relative_irf
hh = dyn_figure(options_,'Name',['Relative response to orthogonalized shock to ' tit(i,:)]);
hh = dyn_figure(options_.nodisplay,'Name',['Relative response to orthogonalized shock to ' tit(i,:)]);
elseif subplotnum == 1 && ~options_.relative_irf
hh = dyn_figure(options_,'Name',['Orthogonalized shock to ' tit(i,:)]);
hh = dyn_figure(options_.nodisplay,'Name',['Orthogonalized shock to ' tit(i,:)]);
end
set(0,'CurrentFigure',hh)
@ -152,14 +152,14 @@ 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_);
if RemoteFlag==1,
dyn_saveas(hh,[DirectoryName '/' M_.fname '_Bayesian_IRF_' deblank(tit(i,:)) '_' int2str(figunumber)],options_.nodisplay,options_.graph_format);
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

@ -25,7 +25,7 @@ function ReshapeMatFiles(type, type2)
% 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 ];

View File

@ -21,7 +21,8 @@ function [fval,info,exit_flag,DLIK,Hess,SteadyState,trend_coeff] = TaRB_optimiz
% 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

@ -14,7 +14,7 @@ function [] = Tracing()
% SPECIAL REQUIREMENTS
% none
% Copyright (C) 2010 Dynare Team
% Copyright (C) 2010-2017 Dynare Team
%
% This file is part of Dynare.
%

View File

@ -19,7 +19,7 @@ function [oo_] = UnivariateSpectralDensity(M_,oo_,options_,var_list)
% Adapted from th_autocovariances.m.
% Copyright (C) 2006-2015 Dynare Team
% Copyright (C) 2006-2017 Dynare Team
%
% This file is part of Dynare.
%
@ -134,7 +134,7 @@ 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
@ -159,12 +159,12 @@ if options_.nograph == 0
end
for i= 1:nvar
hh = dyn_figure(options_,'Name',['Spectral Density of ' deblank(M_.endo_names(ivar(i),:)) '.']);
hh = dyn_figure(options_.nodisplay,'Name',['Spectral Density of ' deblank(M_.endo_names(ivar(i),:)) '.']);
plot(freqs,f(i,:),'-k','linewidth',2)
xlabel('0 \leq \omega \leq \pi')
ylabel('f(\omega)')
box on
axis tight
dyn_saveas(hh,[M_.fname ,filesep,'graphs', filesep, 'SpectralDensity_' deblank(M_.endo_names(ivar(i),:))],options_)
dyn_saveas(hh,[M_.fname ,filesep,'graphs', filesep, 'SpectralDensity_' deblank(M_.endo_names(ivar(i),:))],options_.nodisplay,options_.graph_format)
end
end

View File

@ -0,0 +1,126 @@
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
% endo_names [exo_nbr*string length] variable names from M_.endo_names
% i_var [n_var*1] vector indices of requested variables in M_.endo_names and z
% initial_date [dseries object] first period of decomposition to plot
% DynareModel [structure] Dynare model structure
% DynareOptions [structure] Dynare options structure
% Copyright (C) 2016-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/>.
SteadyState=[];
fig_mode='';
fig_mode1='';
fig_name='';
screen_shocks=0;
use_shock_groups = DynareOptions.plot_shock_decomp.use_shock_groups;
if use_shock_groups
shock_groups = DynareModel.shock_groups.(use_shock_groups);
shock_ind = fieldnames(shock_groups);
end
% number of components equals number of shocks + 1 (initial conditions)
comp_nbr = size(z,2)-1;
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 '_'];
end
if isfield(opts_decomp,'screen_shocks')
if use_shock_groups
screen_shocks=0;
elseif comp_nbr>18
screen_shocks = opts_decomp.screen_shocks;
end
end
if isfield(opts_decomp,'fig_name')
fig_name = opts_decomp.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_'];
end
end
gend = size(z,3);
if isempty(initial_date)
x = 1:gend;
else
freq = initial_date.freq;
initial_period = initial_date.time(1) + (initial_date.time(2)-1)/freq;
x = initial_period:(1/freq):initial_period+(gend-1)/freq;
end
nvar = length(i_var);
labels = char(char(shock_names),'Initial values');
if ~(screen_shocks && comp_nbr>18)
screen_shocks=0;
end
comp_nbr0=comp_nbr;
%%plot decomposition
for j=1:nvar
d0={};
z1 = squeeze(z(i_var(j),:,:));
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);
if use_shock_groups
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
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),:)));
else
[STATUS] = xlwrite([DynareModel.fname,'_shock_decomposition',fig_mode,fig_name1],d0,deblank(endo_names(i_var(j),:)));
end
warning on
clear d0
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

View File

@ -1,6 +1,6 @@
function mexpath = add_path_to_mex_files(dynareroot, modifypath)
% Copyright (C) 2015-2016 Dynare Team
% Copyright (C) 2015-2017 Dynare Team
%
% This file is part of Dynare.
%
@ -22,7 +22,11 @@ if nargin<2
end
if exist('OCTAVE_VERSION')
if ispc() && strcmpi(computer(), 'i686-w64-mingw32')
mexpath = {[dynareroot '../mex/octave32/']};
else
mexpath = {[dynareroot '../mex/octave/']};
end
if modifypath
addpath(mexpath{1});
end
@ -48,7 +52,7 @@ else
end
end
else
tmp = [dynareroot '../mex/matlab/win64-7.8-9.1/'];
tmp = [dynareroot '../mex/matlab/win64-7.8-9.2/'];
if exist(tmp, 'dir')
mexpath = tmp;
if modifypath

View File

@ -0,0 +1,333 @@
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.
% 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
% 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
% i_var: [array] index of vars
% t0: [integer] first period
% t1: [integer] last period
% 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
% endo_names: [char] updated var names
% endo_names_tex: [char] updated TeX var names
% steady_state: [array] updated stady state of vars
% i_var: [integer array] updated var indices to plot
% oo_: [structure] Storage of results
%
% SPECIAL REQUIREMENTS
% none
% Copyright (C) 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/>.
opts = options_.plot_shock_decomp;
nvar = length(i_var);
GYTREND0 = q2a.GYTREND0;
var_type = q2a.type;
islog = q2a.islog;
aux = q2a.aux;
aux0 = aux;
cumfix = q2a.cumfix;
% usual shock decomp
if isstruct(oo_)
% z = oo_.shock_decomposition;
myopts=options_;
myopts.plot_shock_decomp.type='qoq';
myopts.plot_shock_decomp.realtime=0;
[z, junk] = plot_shock_decomposition(M_,oo_,myopts,[]);
else
z = oo_;
end
z = z(i_var,:,:);
mytype=var_type;
if isfield(q2a,'name')
mytxt = q2a.name;
mytex = q2a.name;
if isfield(q2a,'tex_name')
mytex = q2a.tex_name;
end
if mytype==2
gtxt = ['PHI' mytxt]; % inflation rate
gtex = ['{\pi(' mytex ')}'];
elseif mytype
gtxt = ['G' mytxt]; % inflation rate
gtex = ['{g(' mytex ')}'];
end
if isfield(q2a,'gname')
gtxt = q2a.gname;
end
if isfield(q2a,'tex_gname')
gtex = q2a.tex_gname;
end
mytype=0;
end
if isstruct(aux)
if ischar(aux.y)
myopts=options_;
myopts.plot_shock_decomp.type='qoq';
myopts.plot_shock_decomp.realtime=0;
[y_aux, steady_state_aux] = plot_shock_decomposition(M_,oo_,myopts,aux.y);
aux.y=y_aux;
aux.yss=steady_state_aux;
end
yaux=aux.y;
end
if mytype==2
gtxt = 'PHI'; % inflation rate
gtex = '\pi';
elseif mytype
gtxt = 'G'; % growth rate
gtex = 'g';
end
steady_state=steady_state(i_var);
% endo_names = M_.endo_names(i_var,:);
% endo_names_tex = M_.endo_names_tex(i_var,:);
nterms = size(z,2);
nfrcst = opts.forecast/4;
for j=1:nvar
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,:)]);
gendo_names_tex = char(gendo_names_tex,[gtex '(' deblank(endo_names_tex(j,:)) ')']);
else
if nvar==1 && ~mytype
endo_names = mytxt;
endo_names_tex = mytex;
gendo_names = gtxt;
gendo_names_tex = gtex;
else
endo_names = [deblank(M_.endo_names(i_var(j),:)) '_A'];
endo_names_tex = ['{' deblank(M_.endo_names_tex(i_var(j),:)) '}^A'];
gendo_names = [gtxt endo_names(j,:)];
gendo_names_tex = [gtex '(' deblank(endo_names_tex(j,:)) ')'];
end
end
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
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
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
z=gza;
endo_names = gendo_names;
endo_names_tex = gendo_names_tex;
elseif q2a.plot == 2
z=za;
else
z=cat(1,za,gza);
endo_names = char(endo_names,gendo_names);
endo_names_tex = char(endo_names_tex,gendo_names_tex);
end
% if isstruct(oo_)
% oo_.annualized_shock_decomposition=z;
% end
% realtime
if realtime_ && isstruct(oo_) && isfield(oo_, 'realtime_shock_decomposition')
init=1;
for i=t0:4:t1
yr=floor(i/4);
za=[];
gza=[];
myopts=options_;
myopts.plot_shock_decomp.type='qoq';
myopts.plot_shock_decomp.realtime=1;
myopts.plot_shock_decomp.vintage=i;
[z, steady_state_aux] = plot_shock_decomposition(M_,oo_,myopts,[]);
z = z(i_var,:,:);
if isstruct(aux)
if ischar(aux0.y)
[y_aux, steady_state_aux] = plot_shock_decomposition(M_,oo_,myopts,aux0.y);
aux.y=y_aux;
aux.yss=steady_state_aux;
end
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
% 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)
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(ztmp,steady_state(j),GYTREND0,var_type,islog,aux);
% if k<nterms
% 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
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
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
z=gza;
elseif q2a.plot == 2
z=za;
else
z=cat(1,za,gza);
end
if init==1
oo_.annualized_realtime_shock_decomposition.pool = z;
else
oo_.annualized_realtime_shock_decomposition.pool(:,:,yr) = z(:,:,end-nfrcst);
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
oo_.annualized_realtime_conditional_shock_decomposition.(['yr_' int2str(yr-nfrcst)]) = ...
oo_.annualized_realtime_shock_decomposition.pool(:,:,yr-nfrcst:end) - ...
oo_.annualized_realtime_forecast_shock_decomposition.(['yr_' int2str(yr-nfrcst)]);
% fix others
oo_.annualized_realtime_conditional_shock_decomposition.(['yr_' int2str(yr-nfrcst)])(:,end-1,:) = ...
oo_.annualized_realtime_conditional_shock_decomposition.(['yr_' int2str(yr-nfrcst)])(:,end-1,:) + ...
oo_.annualized_realtime_forecast_shock_decomposition.(['yr_' int2str(yr-nfrcst)])(:,end,:);
% fix total
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
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);
oo_.annualized_realtime_conditional_shock_decomposition.(['yr_' int2str(yr-my_forecast_)])(:,end-1,:) = ...
oo_.annualized_realtime_forecast_shock_decomposition.(['yr_' int2str(yr-my_forecast_)])(:,end,1:my_forecast_+1);
oo_.annualized_realtime_conditional_shock_decomposition.(['yr_' int2str(yr-my_forecast_)])(:,end,:) = ...
oo_.annualized_realtime_shock_decomposition.pool(:,end,yr-my_forecast_:yr);
end
end
end
end
% ztmp=oo_.realtime_shock_decomposition.pool(:,:,21:29)-oo_.realtime_forecast_shock_decomposition.time_21;
init=init+1;
end
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))]);
else
error()
end
end
end
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
steady_state = steady_state_ga;
else
steady_state = steady_state_a;
end
end

View File

@ -43,7 +43,7 @@ function [InnovationVariance,AutoregressiveParameters] = autoregressive_process_
%
% \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

@ -0,0 +1,122 @@
function forecasts = backward_model_forecast(initialcondition, listofvariables, periods, withuncertainty)
% Returns unconditional forecasts.
%
% INPUTS
% - initialcondition [dseries] Initial conditions for the endogenous variables.
% - periods [integer] scalar, the number of (forecast) periods.
% - withuncertainty [logical] scalar, returns confidence bands if true.
%
% OUTPUTS
% - forecast [dseries]
% Copyright (C) 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/>.
global M_ options_ oo_
% Check that the model is actually backward
if M_.maximum_lead
error(['backward_model_irf:: The specified model is not backward looking!'])
end
% Initialize returned argument.
forecasts = struct();
% Set defaults.
if nargin<2
listofvariables = cellstr(M_.endo_names);
periods = 8;
withuncertainty = false;
end
if nargin<3
periods = 8;
withuncertainty = false;
end
if nargin<4
withuncertainty = false;
end
% Get full list of endogenous variables
endo_names = cellstr(M_.endo_names);
% Get vector of indices for the selected endogenous variables.
n = length(listofvariables);
idy = zeros(n,1);
for i=1:n
j = strmatch(listofvariables{i}, endo_names, 'exact');
if isempty(j)
error('backward_model_forecast:: Variable %s is unknown!', listofvariables{i})
else
idy(i) = j;
end
end
% Set the number of simulations (if required).
if withuncertainty
B = 1000;
end
% Get the covariance matrix of the shocks.
if withuncertainty
Sigma = M_.Sigma_e + 1e-14*eye(M_.exo_nbr);
sigma = transpose(chol(Sigma));
end
% Set initial condition.
if isdates(initialcondition)
if isempty(M_.endo_histval)
error('backward_model_irf: histval block for setting initial condition is missing!')
end
initialcondition = dseries(transpose(M_.endo_histval), initialcondition, endo_names, cellstr(M_.endo_names_tex));
end
% Put initial conditions in a vector of doubles
initialconditions = transpose(initialcondition{endo_names{:}}.data);
% Compute forecast without shock
innovations = zeros(periods+max(M_.maximum_exo_lag, 1), M_.exo_nbr);
if M_.maximum_exo_lag
if isempty(M_.exo_histval)
error('You need to set the past values for the exogenous variables!')
else
innovations(1:M_.maximum_exo_lag, :) = M_.exo_histval;
end
end
oo__0 = simul_backward_model(initialconditions, periods, options_, M_, oo_, innovations);
forecasts.pointforecast = dseries(transpose(oo__0.endo_simul(idy,:)), initialcondition.init, listofvariables);
if withuncertainty
% Preallocate an array gathering the simulations.
ArrayOfForectasts = zeros(n, periods+1, B);
for i=1:B
innovations(max(M_.maximum_exo_lag, 1)+1:end,:) = transpose(sigma*randn(M_.exo_nbr, periods));
oo__ = simul_backward_model(initialconditions, periods, options_, M_, oo_, innovations);
ArrayOfForecasts(:,:,i) = oo__.endo_simul(idy,:);
end
% Compute mean (over future uncertainty) forecast.
forecasts.meanforecast = dseries(transpose(mean(ArrayOfForecasts, 3)), initialcondition.init, listofvariables);
forecasts.medianforecast = dseries(transpose(median(ArrayOfForecasts, 3)), initialcondition.init, listofvariables);
forecasts.stdforecast = dseries(transpose(std(ArrayOfForecasts, 1,3)), initialcondition.init, listofvariables);
% Compute lower and upper 95% confidence bands
ArrayOfForecasts = sort(ArrayOfForecasts, 3);
forecasts.lb = dseries(transpose(ArrayOfForecasts(:,:,round(0.025*B))), initialcondition.init, listofvariables);
forecasts.ub = dseries(transpose(ArrayOfForecasts(:,:,round(0.975*B))), initialcondition.init, listofvariables);
end

View File

@ -0,0 +1,121 @@
function [endogenousvariables, exogenousvariables] = backward_model_inversion(constraints, exogenousvariables, initialconditions, endo_names, exo_names, freeinnovations, DynareModel, DynareOptions, DynareOutput)
% 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).
% - endo_names [cell] list of endogenous variable names.
% - exo_names [cell] list of exogenous variable names.
% - freeinstruments [cell] list of exogenous variable names used to control the constrained endogenous variables.
%
% OUTPUTS
% - endogenous [dseries]
% - exogenous [dseries]
%
% REMARKS
% Copyright (C) 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/>.
% Get indices for the calibrated and free innovations.
freeinnovations_id = zeros(length(freeinnovations), 1);
if length(freeinnovations)<DynareModel.exo_nbr
for i=1:length(freeinnovations)
freeinnovations_id(i) = strmatch(freeinnovations{i}, exo_names, 'exact');
end
calibratedinnovations_id = setdiff(transpose(1:length(exo_names)), freeinnovations_id);
else
freeinnovations_id = transpose(1:length(exo_names));
calibratedinnovations_id = [];
end
nxfree = length(freeinnovations_id);
nxcalb = length(calibratedinnovations_id);
% Get indices for the the controlled and free endogenous variables.
controlledendogenousvariables_id = zeros(length(freeinnovations), 1);
if length(freeinnovations)<DynareModel.endo_nbr
for i=1:length(freeinnovations)
controlledendogenousvariables_id(i) = strmatch(constraints.name{i}, endo_names, 'exact');
end
freeendogenousvariables_id = setdiff(transpose(1:length(endo_names)), controlledendogenousvariables_id);
else
controlledendogenousvariables_id = transpose(1:length(endo_names));
freeendogenousvariables_id = [];
end
nyfree = length(freeendogenousvariables_id);
nyctrl = length(controlledendogenousvariables_id);
% Get indices of variables appearing at time t-1.
iy1 = find(DynareModel.lead_lag_incidence(1,:)>0);
% Get indices of variables appearing at time t.
iy0 = find(DynareModel.lead_lag_incidence(2,:)>0);
% Set indices for trust_region algorithm.
idx = 1:DynareModel.endo_nbr;
jdx = 1:(nyfree+nxfree);
% Build structure to be passed to the objective function.
ModelInversion.nyfree = nyfree;
ModelInversion.nyctrl = nyctrl;
ModelInversion.nxfree = nxfree;
ModelInversion.nxcalb = nxcalb;
ModelInversion.y_constrained_id = vec(DynareModel.lead_lag_incidence(2,controlledendogenousvariables_id));
ModelInversion.y_free_id = vec(DynareModel.lead_lag_incidence(2,freeendogenousvariables_id));
ModelInversion.x_free_id = freeinnovations_id;
ModelInversion.J_id = [ModelInversion.y_free_id ; sum(DynareModel.lead_lag_incidence(:)>0)+ModelInversion.x_free_id];
% Get the name of the dynamic model routines.
model_dynamic = str2func([DynareModel.fname,'_dynamic']);
model_dtransf = str2func('dynamic_backward_model_for_inversion');
% Initialization of vector y (free endogenous variables and free innovations).
y = NaN(nyfree+nxfree);
% Initialization of the returned simulations (endogenous variables).
Y = NaN(DynareModel.endo_nbr, nobs(constraints)+1);
initialconditions
constraints.dates(1)
Y(:,1) = initialconditions(constraints.dates(1)-1).data(1:DynareModel.endo_nbr);
for i=1:nyctrl
Y(controlledendogenousvariables_id(i),2:end) = transpose(constraints.data(:,i));
end
% Initialization of the returned simulations (exogenous variables).
X = exogenousvariables{exo_names{:}}(constraints.dates(1)-max(1,DynareModel.maximum_exo_lag):constraints.dates(end)).data;
% Inversion of the model, solvers for the free endogenous and exogenous variables (call a Newton-like algorithm in each period).
for it = 2:nobs(constraints)+1
% Set the lagged values of the endogenous variables.
ylag = Y(iy1,it-1);
% Set the current values of the constrained endogenous variables.
ycur = Y(controlledendogenousvariables_id,it);
% Vector z gather the free endogenous variables (initialized with lagged
% values) and the free exogenous variables (initialized with 0).
z = [Y(freeendogenousvariables_id,it-1); zeros(nxfree, 1)];
% Solves for z.
z = dynare_solve(model_dtransf, z, DynareOptions, model_dynamic, ylag, ycur, X, DynareModel.params, DynareOutput.steady_state, it+DynareModel.maximum_exo_lag, ModelInversion);
% Update the matrix of exogenous variables.
X(it,freeinnovations_id) = z(nyfree+1:end);
% Update the matrix of endogenous variables.
Y(freeendogenousvariables_id,it) = z(1:nyfree);
end
endogenousvariables = dseries(Y', constraints.dates(1)-1, endo_names);
exogenousvariables = dseries(X(max(DynareModel.maximum_exo_lag,1)+1:end,:), constraints.dates(1), exo_names);

View File

@ -0,0 +1,119 @@
function irfs = backward_model_irf(initialcondition, listofshocks, listofvariables, varargin)
% Returns impulse response functions.
%
% INPUTS
% - initialcondition [dseries,dates] Initial conditions for the endogenous variables, or period 0.
% - listofshocks [cell of strings] The innovations for which the IRFs need to be computed.
% - listofvariables [cell of strings] The endogenous variables which will be returned.
% - periods [integer] scalar, the number of periods.
%
% OUTPUTS
% - irfs [struct of dseries]
%
% REMARKS
% The names of the fields in the returned structure are given by the name
% of the innovations listed in the second input argument. Each field gather
% the associated paths for endogenous variables listed in the third input
% argument.
% Copyright (C) 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/>.
global M_ options_ oo_
% Check that the model is actually backward
if M_.maximum_lead
error(['simul_model_irf:: The specified model is not backward looking!'])
end
% Set default value for the fourth input argument.
if nargin<4
periods = 40;
notransform = true;
else
periods = varargin{1};
end
% Set default value for the last input argument (no transformation).
if nargin<5
notransform = true;
else
notransform = false;
transform = varargin{2};
end
% Get list of all exogenous variables in a cell of strings
exo_names = cellstr(M_.exo_names);
% Get the list of all exogenous variables in a cell of strings
endo_names = cellstr(M_.endo_names);
% Set initial condition.
if isdates(initialcondition)
if isempty(M_.endo_histval)
error('backward_model_irf: histval block for setting initial condition is missing!')
end
initialcondition = dseries(transpose(M_.endo_histval), initialcondition, endo_names, cellstr(M_.endo_names_tex));
end
% Get the covariance matrix of the shocks.
Sigma = M_.Sigma_e + 1e-14*eye(M_.exo_nbr);
sigma = transpose(chol(Sigma));
% Put initial conditions in a vector of doubles
initialconditions = transpose(initialcondition{endo_names{:}}.data);
% Initialization of the returned argument. Each will be a dseries object containing the IRFS for the endogenous variables listed in the third input argument.
irfs = struct();
% Get the covariance matrix of the shocks.
Sigma = M_.Sigma_e + 1e-14*eye(M_.exo_nbr);
sigma = transpose(chol(Sigma));
% Put initial conditions in a vector of doubles
initialconditions = transpose(initialcondition{endo_names{:}}.data);
% Compute the IRFs (loop over innovations).
for i=1:length(listofshocks)
% Get transition paths induced by the initial condition.
innovations = zeros(periods+M_.maximum_exo_lag, M_.exo_nbr);
if ~isempty(M_.exo_histval)
innovations(1:M_.maximum_exo_lag,:) = M_.exo_histval;
end
oo__0 = simul_backward_model(initialconditions, periods, options_, M_, oo_, innovations);
% Add the shock.
j = strmatch(listofshocks{i}, exo_names);
if isempty(j)
error('backward_model_irf: Exogenous variable %s is unknown!', listofshocks{i})
end
innovations(1+M_.maximum_exo_lag,:) = transpose(sigma(:,j));
oo__1 = simul_backward_model(initialconditions, periods, options_, M_, oo_, innovations);
% Transform the endogenous variables
if notransform
endo_simul__0 = oo__0.endo_simul;
endo_simul__1 = oo__1.endo_simul;
else
endo_simul__0 = feval(transform, oo__0.endo_simul);
endo_simul__1 = feval(transform, oo__1.endo_simul);
end
% Instantiate a dseries object (with all the endogenous variables)
allirfs = dseries(transpose(endo_simul__1-endo_simul__0), initialcondition.init, cellstr(M_.endo_names), cellstr(M_.endo_names_tex));
% Extract a sub-dseries object
irfs.(listofshocks{i}) = allirfs{listofvariables{:}};
end

View File

@ -0,0 +1,39 @@
function [r, J] = dynamic_backward_model_for_inversion(z, dynamicmodel, ylag, ycur, x, params, steady_state, it_, ModelInversion)
% Copyright (C) 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/>.
% Set up y
y = zeros(length(ylag)+ModelInversion.nyfree+ModelInversion.nyctrl,1);
y(1:length(ylag)) = ylag;
y(ModelInversion.y_constrained_id) = ycur;
if ModelInversion.nyfree
y(ModelInversion.y_free_id) = z(1:ModelInversion.nyfree);
end
% Update x
x(it_, ModelInversion.x_free_id) = transpose(z(ModelInversion.nyfree+(1:ModelInversion.nxfree)));
if nargout>1
[r, Jacobian] = feval(dynamicmodel, y, x, params, steady_state, it_);
else
r = feval(dynamicmodel, y, x, params, steady_state, it_);
return
end
J = Jacobian(:,ModelInversion.J_id);

View File

@ -0,0 +1,40 @@
function [r, J] = dynamic_backward_model_for_simulation(z, dynamicmodel, ylag, x, params, steady_state, it_)
% Copyright (C) 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/>.
% Get indices of the variables appearing at time t.
% NOTE: It is assumed that all variables appear at time t in the model.
idy = length(ylag)+(1:length(z));
% Build y vector to be passed to the dynamic model.
y = zeros(length(ylag)+length(z), 1);
y(1:length(ylag)) = ylag;
y(idy) = z;
if nargout>1
% Compute residuals and jacobian of the full dynamic model.
[r, Jacobian] = feval(dynamicmodel, y, x, params, steady_state, it_);
else
% Compute residuals and return.
r = feval(dynamicmodel, y, x, params, steady_state, it_);
return
end
% If the jacobian is computed, remove the columns related to the innovations
% and the variables appearing at time t-1.
J = Jacobian(:,idy);

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.
%

View File

@ -1,4 +1,4 @@
function DynareOutput = simul_backward_linear_model(initial_conditions, sample_size, DynareOptions, DynareModel, DynareOutput, innovations)
function DynareOutput = simul_backward_model(initial_conditions, sample_size, DynareOptions, DynareModel, DynareOutput, innovations)
%@info:
%! @deftypefn {Function File} {@var{DynareOutput} =} simul_backward_nonlinear_model (@var{sample_size},@var{DynareOptions}, @var{DynareModel}, @var{DynareOutput})
@ -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.
%
@ -76,10 +76,14 @@ if nargin<6
otherwise
error(['simul_backward_nonlinear_model:: ' DynareOption.bnlms.innovation_distribution ' distribution for the structural innovations is not (yet) implemented!'])
end
% Put the simulated innovations in DynareOutput.exo_simul.
DynareOutput.exo_simul = zeros(sample_size+1,number_of_shocks);
DynareOutput.exo_simul = zeros(sample_size, number_of_shocks);
DynareOutput.exo_simul(2:end,positive_var_indx) = DynareOutput.bnlms.shocks;
if isfield(DynareModel,'exo_histval') && ~isempty(DynareModel.exo_histval)
DynareOutput.exo_simul = [M_.exo_histval; DynareOutput.exo_simul];
else
DynareOutput.exo_simul = [zeros(1,number_of_shocks); DynareOutput.exo_simul];
end
else
number_of_shocks = size(innovations,2);
DynareOutput.exo_simul = innovations;
@ -92,4 +96,3 @@ else
DynareOutput = simul_backward_nonlinear_model(initial_conditions, sample_size, DynareOptions, ...
DynareModel, DynareOutput, innovations);
end

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