Commit 666c9b80 authored by Willi Mutschler's avatar Willi Mutschler Committed by Sébastien Villemot

Improvement of Identification Toolbox

# Improvements
  * heavily commented (also auxiliary functions) and changed notation to make all the functions (hopefully) more readable and understandable, and hence, easier to debug
  * added identification criteria of Komunjer and Ng (2011, Econometrica) and Qu and Tkachenko (2012, Quantitative Economics)
  * tests can be turned of, i.e. nostrength disables identification strenght, noreducedform disables reduced form criteria, nomoments disables moment criteria, nospectrum disables spectrum criteria, nominimal disables minimal system criteria
  * all kronflags (analytic_derivation_mode) actually work in all functions
  * added functionality when there is correlation in Sigma_e and when one wants to consider corr parameters of exogenous shocks. Previously, (1) corr parameters were not allowed when calling identification and (2) when Sigma_e was not diagonal then the toolbox relied on numerical derviatives only (kronflag=-1). Now it is possible to handle both identification of corr parameters as well as correct analytical derivatives when Sigma_e is not diagonal with all possible kronflag values (-1|-2|0|1)
  * all plots and results are stored in the same folder named identification (previously there was another one with a capital I (Identification))

# Needed changes to preprocessor
  * add as field to options_ident:

    - tex (same as in options_)
    - nostrength (to turn off identification strength)
    - noreducedform (to turn off reduced form criteria)
    - nomoments (to turn off Iskrev's moment criteria)
    - nominimal (to turn off Komunjer and Ng's minimal system criteria)
    - nospectrum (to turn off Qu and Tkachenko's spectrum criteria)

  * add to options_ident:
    - normalize_jacobians (whether to normalize Jacobians or not)
    - grid_nbr (integer used to discretize the interval [-pi;pi]
    - tol_rank (tolerance level to compute ranks)
    - tol_deriv (tolerance level to select nonzero columns in derivatives)
    - tol_sv (tolerance level to select nonzero singular values)
    - ChecksViaSubsets (for debugging purposes, uses different function to find problematic parameter sets)
    - max_dim_subsets_groups (for debugging purposes, used for ChecksViaSubsets)

# Further Suggestions
  * Rename getH.m into getParamsDerivReducedForm.m to make the purpose of this function evident
  * Rename getJJ.m into getIdentificationJacobians.m to make the purpose of this function evident
  * Rename thet2tau.m into IdentificationNumericalObjectiveFunction.m to make the purpose of this function evident
  * dYss, d2Yss, dg1 should also include derivatives wrt to stderr and corr parameters (even though these are just 0), as in other functions (getJJ, dynare_estimation) we always add these manually
  * I am pretty sure the current handling in getH.m of dYss and d2Yss is not correct in the case of nonstationary variables (if g2static is nonempty), I added a warning message, as I am not sure whether this is ever used
  * It would be straigthforward to also include stderr and corr parameters of measurement errors (these is not possible right now). Should I do this?
  * Computations of d2A and d2Om need to be checked, as the differences between computing these with analytically (kronflag=0|1) or numerically kronflag=-1|-2 is really large for the example model of AnSchorfheide.
  * I am not sure how to best normalize Qu and Tkachenko's G matrix. It looks (and in the Gaussian case actually is) very similar to the Ahess matrix. So I used the same normalization rule as for the Ahess matrix. See comments in identification_checks.m. Anyone has a better idea? Please also check the models in test/identification/cgg for differences.
  * parts that are unclear to me are marked by a [@wmutschl] tag
  * the run time of tests/identification/as2007.mod increases from 0h01m27s to 0h03m46s (as Qu and Tkachenko's G matrix takes a little while to compute). One could decrease prior_mc=250 to prior_mc=150.

# New functions
  * commutation: Returns Magnus and Neudecker's commutation matrix that solves k*vec(X)=vec(X')
  * DerivABCD: Derivative of X(p)=A(p)*B(p)*C(p)*D(p) w.r.t to p as in Magnus and Neudecker (1999), p. 175
  * DeriveMinimalState: Derives minimal state space system by checking observability and controllability of all possible combinations of variables
  * duplication: Duplication Matrix (and its Moore Penrose Inverse) as defined by Magnus and Neudecker (2002), p.49, Dp*vec(X) = X
  * identification_checks_via_subsets: finds problematic parameters in a bruteforce fashion: It computes the rank of the Jacobians for all possible parameter combinations, if the rank condition is not fullfilled, these parameter sets are flagged as non-identifiable. For debugging purposes only, as the current identification_checks.m (based on nullspace and multicorrelation coefficients) is much faster

# Detailed changes in getH.m
  * functionality improvements

    - heavily commented (also auxiliary functions) and changed notation of several variables to make this function (hopefully) more readable and understandable, and hence, easier to debug
    - added functionality when Sigma_e is not diagonal and/or when one wants to consider corr parameters of exogenous shocks independent of the value of kronflag
    - fixed function for all values of kronflag, i.e. kronflag=-2|-1|0|1. Previosuly, only kronflag=-2|0 were working, all other kronflags ran into errors (-1 was actually never called , but was dealt with in getJJ.m). I assume kronflag=-1|1 was used only for debugging issues, but still was not working. I fixed this now, the function now works out-of-the-box for all kronflag values.
    - I also outlined and documented what each kronflag does and point to the corresponding equations in Ratto and Iskrev (2012) or Iskrev (2010,Appendix A)
    - the function additionally outputs the Jacobians of B and Sig, which are needed for Qu and Tkachenko (2012) and Komunjer and Ng (2011)'s criteria
    - Moved computation of Jacobian of tau=[ys;vec(A);vech(B * M_.Sigma_e * B')] into getJJ.m to have all Jacobians which are needed for identification in one place. That is, getH.m computes first and second parameter derivatives of (1) reduced-form solution, (2) steady state and (3) Jacobian of dynamic model, whereas getJJ computes and sets up all Jacobians which are used for identification purposes. Therefore, getH might be useful more generally for other purposes than identification. For instance, when doing a GMM estimation, we could use this function to compute analytically the gradient of the moments and provide this to the optimizer used in a GMM context.

  * output arguments

    - renamed `H` (Jacobian wrt parameters of tau=[ys;vec(A);vech(B * M_.Sigma_e * B')] into dTAU, (as H is very confusing, e.g. in other functions it is a Hessian, or Hss and H2ss is also just the steady state. Morevoer, tau is used in Iskrev(2010) for the steady state and reduced-form solution)
    - renamed `Hss` (Jacobian of steady state wrt model parameters only) into `dYss` (as H is very confusing here, see above)
    - renamed `H2ss` (Hessian wrt model parameters only of ys) into d2Yss (as H is very confusing, see above)
    - renamed `gp` into `dg1`, where g1 corresponds to the same variable as in dynamic model files. Note that in params_deriv files gp lacks the contribution of Jacobian wrt steady state and dg1 includes this using the implicit function theorem as outlined in Ratto and Iskrev (2012). Hence, dg1 denotes Jacobian wrt to parameters. It is useful and important to distinguish gp and dg1.
    - added `dB` (Jacobian wrt parameters of solution matrix B) needed for Qu and Tkachenko (2012) as well as Komunjer and Ng (2011)
    - added `dSig` (Jacobian wrt parameters of M_.Sigma_e) needed for Qu and Tkachenko (2012) as well as Komunjer and Ng (2011)

  * input arguments

    - renamed `indx` (index of model parameters to be checked) into `indpmodel`, the p makes it more clear that this is a parameter index
    - renamed `indexo` (index of stderr parameters) into `indpstderr`, the p makes it more clear that this is a parameter index
    - renamed `iv` (index of variables to consider) into `indvar`
    - Renamed `M_` to `M`, `estim_params_` to `estim_params`, `options_` to `options` , `oo_` to `oo` to visualize that these are local and not global variables
    - included `indpcorr` a matrix of indices for corr parameters to be checked

  * misc

    - distinguished clearly between variables in DR or in declaration order without overwriting this in between
    - added which functions call getH.m
    - updated copyright to 2010-2019

# Detailed changes in getJJ.m

  * functionality improvements

    - heavily commented and changed notation of several variables to make this function (hopefully) more readable and understandable, and hence, easier to debug
    - added functionality when Sigma_e is not diagonal and/or when one wants to consider corr parameters of exogenous shocks independent of the value of kronflag
    - tidied the function up, such that it sets up all Jacobians which are needed for identification, i.e. Iskrev's J matrix, Qu and Tkachenko (2012)'s G matrix, Komunjer and Ng (2011)'s D matrix, reduced-form solution (dTAU), linear rational expectation (i.e. Jacobian of steady state and dynamic model equations dLRE).
    - dTAU is now constructed in getJJ instead of in getH (see comment above in getH.m)
    - works for all kronflags, i.e. for numerical derivatives (-1 and -2) as well as for analytical derivatives based on kronecker products (1) or Sylvester Equations (0)
    - added functionality for stderr and corr parameters independent of the value of kronflag (previously this was only possible with numerical derivatives, now it works for all kronflags)
    - finds minimal state vector needed for Komunjer and Ng (2011)'s criteria (function `DeriveMinimalState.m`)
    - moved computations from kronflag=-1 (which were used in case of corr in shock block) into getH.m, so that getJJ now only sets up the Jacobians for LRE, Iskrev's J, Qu and Tkachenko's G and Komunjer and Ng's D, whereas getH computes the Jacobians (wrt parameters) of A, B, Sigma_e, Om, Yss and g1. This should simplify debugging as everything is now in one place and not in two

  * output arguments

    - renamed `JJ` into `J`
    - renamed `H` into `dTAU` (as H is very confusing, e.g. in other functions it is a Hessian, or Hss and H2ss is also just the steady state. Morevoer, tau is used in Iskrev(2010) for the steady state and reduced-form solution)
    - renamed `gp` into `dLRE`, as this corresponds to Jacobian of LRE=[Yss;vec(g1)] where g1 is the Jacobian of the dynamic model equations.
    - renamed `gam` into `MOMENTS`
    - added `G` for Qu and Tkachenko's Jacobian matrix G
    - added `D` for Komunjer and Ng's Jacobian matrix D
    - reordered output arguments

  * input arguments

    - added `options_ident` as input argument; hence, `kronflag`, `nlags` and `useautocorr` are removed from input arguments as these are available in options_ident
    - Renamed `M_` to `M`, `estim_params_` to `estim_params`, `options_` to `options` , `oo_` to `oo` to visualize that these are local and not global variables
    - renamed `indx` (index of model parameters to be checked) into `indpmodel`, the p makes it more clear that this is a parameter index
    - renamed `indexo` (index of stderr parameters) into `indpstderr`, the p makes it more clear that this is a parameter index
    - added `indpcorr` (index of corr parameters)
    - renamed `mf` (index of VAROBS variables) into `indvobs`

  * misc

    - updated copyright to 2010-2019
    - provided some comments on several ways to compute the spectral density matrix
    - added which functions call getJJ.m

# Detailed changes in thet2tau.m

  * functionality improvements

    - heavily commented and changed notation of several variables to make this function (hopefully) more readable and understandable, and hence, easier to debug
    - Added output option to compute spectral density matrix
    - Reorded and added some output option.
    - Instead of Om, `outputflag=0` computes B and Sigma_e, which are needed for Qu and Tkachenko as well as Komunjer and Ng. The Jacobian of Om is then computed in getJJ or getH from Jacobian of B and Sigma_e. Due to some testing with An and Schorfheide model this seems to be more accurate when I compare these with the analytical derivatives. The old behavior (computing Om directly) can be restored by setting `outputflag=-2`.
    - In total this function can now be used to compute numerically Jacobians of Yss, A, B, Sigma_e, Om, g1, autocovariogram and spectral density
    - Clearly distinguished (and commented) on the different outputs of this function.
    - Works for all types of parameters, ie. model, stderr and corr.
    - This function can now also be used when there is no estimated_params block. Previously, there was an error when there was no estimated_params block when calling `set_all_parameters` as this requires some information in `estim_params`. I fixed this by providing a temporary local estim_parms structure with the necessary information on model, stderr and corr parameters. In this way, this can be easily extended to also include stderr and corr parameters of measurement errors.

  * output arguments

    - renamed `tau` into `out`, as this function computes *very* different things (and not only tau) depending on an input flag

  * input arguments

    - renamed `flagmoments` into `outputflag` as this function does not only compute moments but many other things (see above)
    - renamed `indx` (index of model parameters to be checked) into `indpmodel`, the p makes it more clear that this is a parameter index
    - renamed `indexo` (index of stderr parameters) into `indpstderr`, the p makes it more clear that this is a parameter index
    - added `indpcorr` (index of corr parameters)
    - merged `mf` (index of observable variables) and `iv` (index of variables to consider) into a single index `indvar` as there is no need to distinguish between these two indices (they were never used in combination)
    - added `grid_nbr` (number of grid points to compute spectral density)
    - reordered input arguments

  * misc

    - added which functions call thet2tau
    - updated copyright to 2010-2019

# Detailed changes in identification_analysis.m

  * functionality improvements

    - heavily commented and changed notation of several variables to make this function (hopefully) more readable and understandable, and hence, easier to debug
    - renamed `dg1` to `dLRE`, renamed `vecg1` to `lre`, renamed `H` to `dTAU` (see comments above)
    - added option `numzerotolderiv` with default `1.e-8` used for non-zero derivatives
    - added option `numzerotolrank` with default `1.e-10` used for rank computations
    - added theoretical identification analysis based on Komunjer and Ng (2011)'s method, i.e. steady state and observational equivalent spectral densities within a minimal system
    - added theoretical identification analysis based on Qu and Tkachenko (2012)'s method, i.e. steady state and spectral density
    - restructured the code slightly to combined chunks of code that belong together on the one hand, and on the other hand to differentiate between the different criteria
    - added call to new function `identification_checks_via_subsets.m` (see above for the definition of the functionality) to perform identification checks differently as find it more intuitive and (most likely) more precise.

  * input arguments

    - removed `bounds` and `dataset_` as input argument, because these are not needed
    - moved `name_tex` and `tittxt` into `options_ident` as these two inputs are only used in `ident_bruteforce.m` and already set in `dynare_identification.m`

  * output arguments

    - added `ide_spectrum` structure for Qu and Tkachenko's criteria based on the spectral density
    - added `ide_minimal` structure for Komunjer and Ng's criteria based on the minimal state space system
    - reordered output arguments

  * misc

    - added which functions call identification_analysis
    - updated copyright to 2010-2019

# Detailed changes in dynare_identification.m

  * functionality improvements

    - heavily commented and changed notation of several variables to make this function (hopefully) more readable and understandable, and hence, easier to debug
    - included more options (and default values) which can be set by the user, i.e. nostrength, nomoments, nominimal, nospectrum, tex, tol_rank, tol_deriv, tol_sv, grid_nbr, ChecksViaSubsets, max_dim_subsets_group
    - instead of turning warnings globally off, I specified the relevant warnings for matlab and octave, respectively, off
    - improved the warning messages slightly
    - restructured chunks of code with respect to different criteria

  * output arguments

    - renamed arguments: TAU to STO_TAU, GAM to STO_MOMENTS, LRE to STO_LRE, gp to STO_si_dLRE, H to STO_si_dTAU, JJ to STO_si_J
    - added arguments: STO_G and STO_D for the two new criteria

  * misc

    - added which functions call dynare_identification
    - updated copyright to 2010-2019

# Detailed changes in identification_checks.m

  * functionality improvements

    - added checks for Komunjer and Ng's D matrix. Note that the Jacobian D=[D_par D_rest], where D_par depends on the parameters and D_rest does not. So this is taken into account.
    - added checks for Qu and Tkachenko's G matrix. Note that the Jacobian G is a Gram matrix with dimension nparam x nparam, similar to Ahess. So this is taken into account. I am, however, not sure whether this is correct regarding the multicorrelation and pairwise correlation coefficients. Please double check.
    - the rank is now actually computed at the prespecified tolerance level (and not Matlab's default level), so this is in accordance to the further analysis of problematic parameter sets

  * output arguments

    - added the rank to output arguments which is later also displayed
    - replaced the J or JJ part in the variable names with X as this function is used for all sorts of Jacobians, not only Iskrev's J

  * input arguments

    - renamed hess_flag to output_flag (and clearly outlined what each value does)
    - added tol_rank and tol_sv as input arguments, such that the tolerance levels can be changed by the user and not preimplemented in this function
    - added param_nbr which is needed for Komunjer and Ng's D matrix

  * misc

    - updated copyright to 2010-2019

# Detailed changes in ident_bruteforce.m

  * functionality improvements

    - the output directory was set with a capital I, i.e. Identification, whereas in all other functions we rely on lower case i, i.e. identification. I changed this to lower-cases, so everything is now saved in the same folder.
    - changed displayed strings to be more precise with the corresponding papers and notation

  * input arguments

    - renamed `n` to `max_dim_cova_group` to name options the same across functions
    - renamed `pnames_TeX` to `name_tex` to name options the same across functions
    - added `tol_deriv` as tolerance level which can be changed by the user

  * misc

    - Added some comments
    - updated copyright to 2010-2019

# Detailed changes in disp_identification.m

  * functionality improvements

    - this function displays the same output for different Jacobians, hence I put the common code into a for loop. This should simplify changing the output that is printed to the console. Previously the code was simply repeated for the different criteria and only the strings changed.
    - some settings relevant for the computation are now printed as a summary to the console
    - the tolerance level, rank and required rank are always displayed on the command line to see how many problematic sets there are and which tolerance level was used
    - the function is also able to display problematic parameters computed by the new function `identification_checks_via_subsets.m` which is only used for debugging.

  * input arguments

    - added `idespectrum` structure for analysis based on Qu and Tkachenko
    - added `ideminimal` structure for analysis based on Komunjer and Ng
    - added `options_ident` to have all necessary settings in a structure

  * misc

    - Added some comments
    - Removed uncommented code that was not used as this was redundant and probably an artifact of the original programming?!
    - updated copyright to 2010-2019

# Detailed changes in dsge_likelihood.m

  * misc

    - adjusted call of getH due to changes of input and output arguments
    - updated copyright to 2010-2019

# Detailed changes in cosn.m

  * misc

    - commented functionality, input and output arguments of this function
    - updated copyright to 2010-2019
parent 41205954
......@@ -8768,51 +8768,128 @@ Performing identification analysis
.. command:: identification ;
identification (OPTIONS...);
|br| This command triggers identification analysis.
|br| This command triggers:
*Options*
1. Theoretical identification analysis based on
.. option:: ar = INTEGER
* moments as in *Iskrev (2010)*
* spectral density as in *Qu and Tkachenko (2012)*
* minimal system as in *Komunjer and Ng (2011)*
* reduced-form solution and linear rational expectation model
as in *Ratto and Iskrev (2011)*
Number of lags of computed autocorrelations (theoretical
moments). Default: ``1``.
2. Identification strength analysis based on sample information matrix as in
*Ratto and Iskrev (2011)*
.. option:: useautocorr = INTEGER
3. Parameter checks based on nullspace and multicorrelation coefficients to
determine which (combinations of) parameters are involved
If equal to ``1``, compute derivatives of autocorrelation. If
equal to ``0``, compute derivatives of
autocovariances. Default: ``0``.
*General Options*
.. option:: load_ident_files = INTEGER
.. option:: parameter_set = OPTION
If equal to ``1``, allow Dynare to load previously computed
analyzes. Default: ``0``.
Specify the parameter set to use. Possible values for OPTION are:
* ``calibration``
* ``prior_mode``
* ``prior_mean``
* ``posterior_mode``
* ``posterior_mean``
* ``posterior_median``
Default: ``prior_mean``.
.. option:: prior_mc = INTEGER
Size of Monte-Carlo sample. Default: ``1``.
Size of Monte-Carlo sample.
Default: ``1``.
.. option:: prior_range = INTEGER
Triggers uniform sample within the range implied by the prior
specifications (when ``prior_mc>1``). Default: ``0``.
specifications (when ``prior_mc>1``).
Default: ``0``.
.. option:: advanced = INTEGER
Shows a more detailed analysis, comprised of an analysis for
the linearized rational expectation model as well as the
associated reduced form solution. Further performs a brute
force search of the groups of parameters best reproducing the
behavior of each single parameter. The maximum dimension of
the group searched is triggered by
``max_dim_cova_group``. Default: ``0``.
If set to ``1``, shows a more detailed analysis, comprised of
an analysis for the linearized rational expectation model as
well as the associated reduced form solution. Further performs
a bruteforce search of the groups of parameters best reproducing
the behavior of each single parameter. The maximum dimension of
the group searched is triggered by ``max_dim_cova_group``.
Default: ``0``.
.. option:: max_dim_cova_group = INTEGER
In the brute force search (performed when ``advanced=1``) this
option sets the maximum dimension of groups of parameters that
best reproduce the behavior of each single model
parameter. Default: ``2``.
best reproduce the behavior of each single model parameter.
Default: ``2``.
.. option:: gsa_sample_file = INTEGER|FILENAME
If equal to ``0``, do not use sample file. If equal to ``1``,
triggers gsa prior sample. If equal to ``2``, triggers gsa
Monte-Carlo sample (i.e. loads a sample corresponding to
``pprior=0`` and ``ppost=0`` in the ``dynare_sensitivity``
options). If equal to ``FILENAME`` uses the provided path to
a specific user defined sample file.
Default: ``0``.
.. option:: diffuse_filter
Deals with non-stationary cases. See :opt:`diffuse_filter`.
*Numerical Options*
.. option:: analytic_derivation_mode = INTEGER
Different ways to compute derivatives either analytically or numerically.
Possible values are:
* ``0``: efficient sylvester equation method to compute
analytical derivatives
* ``1``: kronecker products method to compute analytical
derivatives
* ``-1``: numerical two-sided finite difference method
to compute all identification Jacobians
* ``-2``: numerical two-sided finite difference method
to compute derivatives of steady state and dynamic
model numerically, the identification Jacobians are
then computed analytically
Default: ``0``.
.. option:: normalize_jacobians = INTEGER
If set to ``1``: Normalize Jacobian matrices by rescaling
each row by its largest element in absolute value.
Normalize Gram (or Hessian-type) matrices by transforming
into correlation-type matrices.
Default: ``1``
.. option:: tol_rank = DOUBLE
Tolerance level used for rank computations.
Default: ``1.e-10``.
.. option:: tol_deriv = DOUBLE
Tolerance level for selecting non-zero columns in Jacobians.
Default: ``1.e-8``.
.. option:: tol_sv = DOUBLE
Tolerance level for selecting non-zero singular values.
Default: ``1.e-3``.
*Identification Strength Options*
.. option:: no_identification_strength
Disables computations of identification strength analysis
based on sample information matrix.
.. option:: periods = INTEGER
......@@ -8827,39 +8904,48 @@ Performing identification analysis
number of replicas to compute Simulated Moments
Uncertainty. Default: ``100``.
.. option:: gsa_sample_file = INTEGER
*Moments Options*
If equal to ``0``, do not use sample file. If equal to ``1``,
triggers gsa prior sample. If equal to ``2``, triggers gsa
Monte-Carlo sample (i.e. loads a sample corresponding to
``pprior=0`` and ``ppost=0`` in the ``dynare_sensitivity``
options). Default: ``0``.
.. option:: no_identification_moments
.. option:: gsa_sample_file = FILENAME
Disables computations of identification check based on
Iskrev (2010)'s J, i.e. derivative of first two moments.
Uses the provided path to a specific user defined sample
file. Default: ``0``.
.. option:: ar = INTEGER
.. option:: parameter_set = OPTION
Number of lags of computed autocovariances/autocorrelations
(theoretical moments) in Iskrev (2010)'s J criteria.
Default: ``1``.
Specify the parameter set to use. Possible values for OPTION are:
.. option:: useautocorr = INTEGER
* ``calibration``
* ``prior_mode``
* ``prior_mean``
* ``posterior_mode``
* ``posterior_mean``
* ``posterior_median``
If equal to ``1``, compute derivatives of autocorrelation. If
equal to ``0``, compute derivatives of
autocovariances. Default: ``0``.
Default: ``prior_mean``.
*Spectrum Options*
.. option:: lik_init = INTEGER
.. option:: no_identification_spectrum
See :opt:`lik_init <lik_init = INTEGER>`.
Disables computations of identification check based on
Qu and Tkachenko (2012)'s G, i.e. Gram matrix of derivatives of
first moment plus outer product of derivatives of spectral density.
.. option:: kalman_algo = INTEGER
.. option:: grid_nbr = INTEGER
See :opt:`kalman_algo <kalman_algo = INTEGER>`.
Number of grid points in [-pi;pi] to approximate the integral
to compute Qu and Tkachenko (2012)'s G criteria.
Default: ``5000``.
*Minimal State Space System Options*
.. option:: no_identification_minimal
Disables computations of identification check based on
Komunjer and Ng (2011)'s D, i.e. minimal state space system
and observational equivalent spectral density transformations.
*Misc Options*
.. option:: nograph
......@@ -8874,6 +8960,47 @@ Performing identification analysis
See :opt:`graph_format <graph_format = FORMAT>`.
.. option:: tex
See :opt:`tex`.
*Debug Options*
.. option:: load_ident_files = INTEGER
If equal to ``1``, allow Dynare to load previously computed
analyzes. Default: ``0``.
.. option:: lik_init = INTEGER
See :opt:`lik_init <lik_init = INTEGER>`.
.. option:: kalman_algo = INTEGER
See :opt:`kalman_algo <kalman_algo = INTEGER>`.
.. option:: no_identification_reducedform
Disables computations of identification check based on
steady state and reduced-form solution.
.. option:: checks_via_subsets = INTEGER
If equal to ``1``: finds problematic parameters in a bruteforce
fashion: It computes the rank of the Jacobians for all possible
parameter combinations. If the rank condition is not fullfilled,
these parameter sets are flagged as non-identifiable.
The maximum dimension of the group searched is triggered by
``max_dim_subsets_groups``.
Default: ``0``.
.. option:: max_dim_subsets_groups = INTEGER
Sets the maximum dimension of groups of parameters for which
the above bruteforce search is performed.
Default: ``4``.
Types of analysis and output files
----------------------------------
......@@ -9238,14 +9365,20 @@ For example, the placing::
dynare_sensitivity(identification=1, morris=2);
in the Dynare model file triggers identification analysis using
analytic derivatives *Iskrev (2010)*, jointly with the mapping of the
acceptable region.
analytic derivatives as in *Iskrev (2010)*, jointly with the mapping
of the acceptable region.
The identification analysis with derivatives can also be triggered by
the commands ``identification;`` This does not do the mapping of
acceptable regions for the model and uses the standard random sampler
of Dynare. It completely offsets any use of the sensitivity analysis
toolbox.
the single command::
identification;
This does not do the mapping of acceptable regions for the model and
uses the standard random sampler of Dynare. Additionally, using only
``identification;`` adds two additional identification checks: namely,
of *Qu and Tkachenko (2012)* based on the spectral density and of
*Komunjer and Ng (2011)* based on the minimal state space system.
It completely offsets any use of the sensitivity analysis toolbox.
......
function k = commutation(n, m, sparseflag)
% k = commutation(n, m, sparseflag)
% -------------------------------------------------------------------------
% Returns Magnus and Neudecker's commutation matrix of dimensions n by m,
% that solves k*vec(X)=vec(X')
% =========================================================================
% INPUTS
% n: [integer] row number of original matrix
% m: [integer] column number of original matrix
% sparseflag: [integer] whether to use sparse matrices (=1) or not (else)
% -------------------------------------------------------------------------
% OUTPUTS
% k: [n by m] commutation matrix
% -------------------------------------------------------------------------
% This function is called by
% * get_first_order_solution_params_deriv.m (previously getH.m)
% * get_identification_jacobians.m (previously getJJ.m)
% -------------------------------------------------------------------------
% This function calls
% * vec (embedded)
% =========================================================================
% Copyright (C) 2019 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/>.
% =========================================================================
% Original author: Thomas P Minka (tpminka@media.mit.edu), April 22, 2013
if nargin < 2
m = n(2);
n = n(1);
end
if nargin < 3
sparseflag = 0;
end
if 0
% first method
i = 1:(n*m);
a = reshape(i, n, m);
j = vec(transpose(a));
k = zeros(n*m,n*m);
for r = i
k(r, j(r)) = 1;
end
else
% second method
k = reshape(kron(vec(eye(n)), eye(m)), n*m, n*m);
end
if sparseflag ~= 0
k = sparse(k);
end
function V = vec(A)
V = A(:);
end
end
\ No newline at end of file
function [co, b, yhat] = cosn(H)
% function co = cosn(H);
% computes the cosine of the angle between the H(:,1) and its
% projection onto the span of H(:,2:end)
%
% Not the same as multiple correlation coefficient since the means are not
% zero
%
% Copyright (C) 2008-2017 Dynare Team
% function [co, b, yhat] = cosn(H)
% -------------------------------------------------------------------------
% computes the cosine of the angle between the (endogenous variable) H(:,1)
% and its projection onto the span of (exogenous variables) H(:,2:end)
% Note: This is not the same as multiple correlation coefficient since the
% means are not zero
% =========================================================================
% INPUTS
% * H [n by k]
% Data matrix, endogenous variable y is in the first column,
% exogenous variables X are in the remaining (k-1) columns
% -------------------------------------------------------------------------
% OUTPUTS
% * co [double] (approximate) multiple correlation coefficient
% * b [k by 1] ols estimator
% * y [n by 1] predicted endogenous values given ols estimation
% -------------------------------------------------------------------------
% This function is called by
% * identification_checks.m
% * ident_bruteforce.m
% =========================================================================
% Copyright (C) 2008-2019 Dynare Team
%
% This file is part of Dynare.
%
......@@ -23,15 +36,16 @@ function [co, b, yhat] = cosn(H)
%
% You should have received a copy of the GNU General Public License
% along with Dynare. If not, see <http://www.gnu.org/licenses/>.
% =========================================================================
y = H(:,1);
X = H(:,2:end);
b=(X\y);
b=(X\y); %ols estimator
if any(isnan(b)) || any(isinf(b))
b=0;
end
yhat = X*b;
yhat = X*b; %predicted values
if rank(yhat)
co = abs(y'*yhat/sqrt((y'*y)*(yhat'*yhat)));
else
......
This diff is collapsed.
......@@ -111,11 +111,11 @@ function [fval,info,exit_flag,DLIK,Hess,SteadyState,trend_coeff,Model,DynareOpti
%! @sp 2
%! @strong{This function calls:}
%! @sp 1
%! @ref{dynare_resolve}, @ref{lyapunov_symm}, @ref{lyapunov_solver}, @ref{compute_Pinf_Pstar}, @ref{kalman_filter_d}, @ref{missing_observations_kalman_filter_d}, @ref{univariate_kalman_filter_d}, @ref{kalman_steady_state}, @ref{getH}, @ref{kalman_filter}, @ref{score}, @ref{AHessian}, @ref{missing_observations_kalman_filter}, @ref{univariate_kalman_filter}, @ref{priordens}
%! @ref{dynare_resolve}, @ref{lyapunov_symm}, @ref{lyapunov_solver}, @ref{compute_Pinf_Pstar}, @ref{kalman_filter_d}, @ref{missing_observations_kalman_filter_d}, @ref{univariate_kalman_filter_d}, @ref{kalman_steady_state}, @ref{get_first_order_solution_params_deriv}, @ref{kalman_filter}, @ref{score}, @ref{AHessian}, @ref{missing_observations_kalman_filter}, @ref{univariate_kalman_filter}, @ref{priordens}
%! @end deftypefn
%@eod:
% Copyright (C) 2004-2018 Dynare Team
% Copyright (C) 2004-2019 Dynare Team
%
% This file is part of Dynare.
%
......@@ -523,10 +523,9 @@ if analytic_derivation
indparam=[];
end
if full_Hess
[dum, DT, DOm, DYss, dum2, D2T, D2Om, D2Yss] = getH(A, B, EstimatedParameters, Model,DynareResults,DynareOptions,kron_flag,indparam,indexo,iv);
clear dum dum2;
[DT, ~, ~, DOm, DYss, ~, D2T, D2Om, D2Yss] = get_first_order_solution_params_deriv(A, B, EstimatedParameters, Model,DynareResults,DynareOptions,kron_flag,indparam,indexo,[],iv);
else
[dum, DT, DOm, DYss] = getH(A, B, EstimatedParameters, Model,DynareResults,DynareOptions,kron_flag,indparam,indexo,iv);
[DT, ~, ~, DOm, DYss] = get_first_order_solution_params_deriv(A, B, EstimatedParameters, Model,DynareResults,DynareOptions,kron_flag,indparam,indexo,[],iv);
end
else
DT = derivatives_info.DT(iv,iv,:);
......
function [Dp,DpMPinv] = duplication(p)
% [Dp,DpMPinv] = duplication(p)
% -------------------------------------------------------------------------
% Duplication Matrix as defined by Magnus and Neudecker (2002), p.49
% =========================================================================
% INPUTS
% p: [integer] length of vector
% -------------------------------------------------------------------------
% OUTPUTS
% Dp: Duplication matrix
% DpMPinv: Moore-Penroze inverse of Dp
% -------------------------------------------------------------------------
% This function is called by
% * get_identification_jacobians.m (previously getJJ.m)
% =========================================================================
% Copyright (C) 2019 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/>.
% =========================================================================
% Original author: Thomas P Minka (tpminka@media.mit.edu), April 22, 2013
a = tril(ones(p));
i = find(a);
a(i) = 1:length(i);
a = a + transpose(tril(a,-1));
j = a(:);
m = p*(p+1)/2;
Dp = spalloc(p*p,m,p^2);
for r = 1:size(Dp,1)
Dp(r, j(r)) = 1;
end
if nargout > 1
DpMPinv = (transpose(Dp)*Dp)\transpose(Dp);
end
\ No newline at end of file
This diff is collapsed.
This diff is collapsed.
function [JJ, H, gam, gp, dA, dOm, dYss] = getJJ(A, B, estim_params_, M_,oo_,options_,kronflag,indx,indexo,mf,nlags,useautocorr)
% function [JJ, H, gam, gp, dA, dOm, dYss] = getJJ(A, B, estim_params_, M_,oo_,options_,kronflag,indx,indexo,mf,nlags,useautocorr)
% computes derivatives of 1st and 2nd order moments of observables with
% respect to estimated parameters
%
% Inputs:
% A: Transition matrix of lagged states from Kalman filter
% B: Matrix in state transition equation mapping shocks today to
% states today
% M_: structure storing the model information
% oo_: structure storing the results
% options_: structure storing the options
% kronflag: Indicator whether to rely on Kronecker products (1) or
% not (-1 or -2)
% indx: Index of estimated parameters in M_.params
% indexo: Index of estimated standard deviations in M_.exo_names
% mf: Index of observed variables
% nlags: Number of lags to consider for covariances and
% correlations
% useautocorr: Indicator on whether to use correlations (1) instead of
% covariances (0)
%
% Outputs:
% JJ: Jacobian of 1st and 2nd order moments of observables, i.e. dgam/dTHETA
% (see below for definition of gam)
% H: dTAU/dTHETA: Jacobian of TAU, vectorized form of
% linearized reduced form state space model, given ys [steady state],
% A [transition matrix], B [matrix of shocks], Sigma [covariance of shocks]
% TAU = [ys; vec(A); dyn_vech(B*Sigma*B')].
% gam: vector of theoretical moments of observed variables mf [JJ is the Jacobian of gam].
% gam = [ys(mf); dyn_vech(GAM{1}); vec(GAM{j+1})]; for j=1:ar and where
% GAM is the first output of th_autocovariances
% gp: Jacobian of linear rational expectation matrices [i.e.
% Jacobian of dynamic model] with respect to estimated
% structural parameters only (indx)
% dA: [endo_nbr by endo_nbr by (indx+indexo)] Jacobian of transition matrix A
% dOm: Jacobian of Omega = (B*Sigma*B')
% dYss Jacobian of steady state with respect to estimated structural parameters only (indx)
% Copyright (C) 2010-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 nargin<8 || isempty(indx)
% indx = [1:M_.param_nbr];
end
if nargin<9 || isempty(indexo)
indexo = [];
end
if nargin<11 || isempty(nlags)
nlags=3;
end
if nargin<12 || isempty(useautocorr)
useautocorr=0;
end
% if useautocorr,
warning('off','MATLAB:divideByZero')
% end
if kronflag == -1
fun = 'thet2tau';
params0 = M_.params;
para0 = get_all_parameters(estim_params_, M_);
JJ = fjaco(fun,para0,estim_params_,M_, oo_, indx,indexo,1,mf,nlags,useautocorr);
M_.params = params0;
params0 = M_.params;
H = fjaco(fun,para0,estim_params_,M_, oo_, indx,indexo,0,mf,nlags,useautocorr);
M_.params = params0;
params0 = M_.params;
gp = fjaco(fun,para0,estim_params_,M_, oo_, indx,indexo,-1);
M_.params = params0;
offset = length(para0)-length(indx);