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

Added parameter derivatives of perturbation solution up to 3 order 

# Preliminary comments
I finished the identification toolbox at orders two and three using the pruned state space system, but before I merge request this, I decided to first merge the new functionality to compute parameter derivatives of perturbation solution matrices at higher orders. So after this is approved, I merge the identification toolbox.
I guess @rattoma, @sebastien, and @michel are best choices to review this.
I outline the main idea first and then provide some more detailed changes I made to the functions.

***

# Main idea
This merge request is concerned with the *analytical*computation of the parameter derivatives of first, second and third order perturbation solution matrices, i.e. using _closed-form_ expressions to efficiently compute the derivative of  $g_x$ , $g_u$, $g_{xx}$, $g_{xu}$, $g_{uu}$, $g_{\sigma\sigma}$, $g_{xxx}$, $g_{xxu}$, $g_{xuu}$, $g_{uuu}$, $g_{x\sigma\sigma}$, $g_{u\sigma\sigma}$ *with respect to model parameters*  $\theta$.  Note that $\theta$ contains model parameters, stderr and corr parameters of shocks. stderr and corr parameters of measurement errors are not yet supported, (they can easily be included as exogenous shocks). The availability of such derivatives is beneficial in terms of more reliable analysis of model sensitivity and parameter identifiability as well as more efficient estimation methods, in particular for models solved up to third order, as it is well-known that numerical derivatives are a tricky business, especially for large models.

References for my approach are:
* Iskrev (2008, 2010) and Schmitt-Grohé and Uribe (2012, Appendix)  who were the first to compute the parameter derivatives analytically at first order, however, using inefficient (sparse) Kronecker products.
* Mutschler (2015) who provides the expressions for a second-order, but again using inefficient (sparse) Kronecker products.
* Ratto and Iskrev (2012) who show how the first-order system can be solved accurately, fast and efficiently using existing numerical algorithms for generalized Sylvester equations by taking the parameter derivative with respect to each parameter separately.
* Julliard and Kamenik (2004) who provide the perturbation solution equation system in tensor notation at any order k.
* Levintal (2017) who introduces permutation matrices to express the perturbation solution equation system in matrix notation up to fifth order.
Note that @rattoma already implemented the parameter derivatives of $g_x$ and $g_u$ analytically (and numerically), and I rely heavily on his work in `get_first_order_solution_params_derivs.m` (previously `getH.m`). My additions are mainly to this function and thus it is renamed to `get_perturbation_params_derivs.m`.

The basic idea of this merge request is to take the second- and third-order perturbation solution systems in Julliard and Kamenik (2004), unfold these into an equivalent matrix representation using permutation matrices as in Levintal (2017). Then extending Ratto and Iskrev (2012) one takes the derivative with respect to each parameter separately and gets a computational problem that is linear, albeit large, as it involves either solving generalized Sylvester equations or taking inverses of highly sparse matrices. I will now briefly summarize the perturbation solution system at third order and the system that results when taking the derivative with respect to parameters.

## Perturbation Solution
The following systems arise at first, second, and third order:
$(ghx): f_{x} z_{x} = f_{y_{-}^*} + f_{y_0} g_{x} + f_{y_{+}^{**}} g^{**}_{x} g^{*}_{x}= A g_{x} + f_{y_{-}^*}=0$

$(ghu): f_{z} z_{u} = f_{y_0} g_{u} + f_{y_{+}^{**}} g^{**}_{x} g^{*}_{u} + f_{u}= A g_u + f_u = 0$

$(ghxx) : A g_{xx} + B g_{xx} \left(g^{*}_{x} \otimes g^{*}_{x}\right) + f_{zz} \left( z_{x} \otimes z_{x} \right) = 0$

$(ghxu) : A g_{xu} + B g_{xx} \left(g^{*}_{x} \otimes g^{*}_{u}\right) + f_{zz} \left( z_{x} \otimes z_{u} \right) = 0$

$(ghuu) : A g_{uu} + B g_{xx} \left(g^{*}_{u} \otimes g^{*}_{u}\right) + f_{zz} \left( z_{u} \otimes z_{u} \right) = 0$

$(ghs2) : (A+B) g_{\sigma\sigma} +  \left( f_{y^{**}_{+}y^{**}_{+}} \left(g^{**}_{u} \otimes g^{**}_{u}\right) + f_{y^{**}_{+}} g^{**}_{uu}\right)vec(\Sigma) = 0$

$(ghxxx) : A g_{xxx} + B g_{xxx} \left(g^{*}_{x} \otimes g^{*}_{x} \otimes g^{*}_{x}\right) + f_{y_{+}}g^{**}_{xx} \left(g^{*}_x \otimes g^{*}_{xx}\right)P_{x\_xx} + f_{zz} \left( z_{x} \otimes z_{xx} \right)P_{x\_xx} + f_{zzz} \left( z_{x} \otimes z_{x} \otimes z_{x} \right) = 0$

$(ghxxu) : A g_{xxu} + B g_{xxx} \left(g^{*}_{x} \otimes g^{*}_{x} \otimes g^{*}_{u}\right) + f_{zzz} \left( z_{x} \otimes z_{x} \otimes z_{u} \right) + f_{zz} \left( \left( z_{x} \otimes z_{xu} \right)P_{x\_xu} + \left(z_{xx} \otimes z_{u}\right) \right) + f_{y_{+}}g^{**}_{xx} \left( \left(g^{*}_{x} \otimes g^{*}_{xu}\right)P_{x\_xu} + \left(g^{*}_{xx} \otimes g^{*}_{u}\right) \right) = 0$

$(ghxuu) : A g_{xuu} + B g_{xxx} \left(g^{*}_{x} \otimes g^{*}_{u} \otimes g^{*}_{u}\right) + f_{zzz} \left( z_{x} \otimes z_{u} \otimes z_{u} \right)+ f_{zz} \left( \left( z_{xu} \otimes z_{u} \right)P_{xu\_u} + \left(z_{x} \otimes z_{uu}\right) \right) + f_{y_{+}}g^{**}_{xx} \left( \left(g^{*}_{xu} \otimes g^{*}_{u}\right)P_{xu\_u} + \left(g^{*}_{x} \otimes g^{*}_{uu}\right) \right) = 0$

$(ghuuu) : A g_{uuu} + B g_{xxx} \left(g^{*}_{u} \otimes g^{*}_{u} \otimes g^{*}_{u}\right) + f_{zzz} \left( z_{u} \otimes z_{u} \otimes z_{u} \right)+ f_{zz} \left( z_{u} \otimes z_{uu} \right)P_{u\_uu} + f_{y_{+}}g^{**}_{xx} \left(g^{*}_{u} \otimes g^{*}_{uu}\right)P_{u\_uu}  = 0$

$(ghx\sigma\sigma) : A g_{x\sigma\sigma} + B g_{x\sigma\sigma} g^{*}_x + f_{y_{+}} g^{**}_{xx}\left(g^{*}_{x} \otimes g^{*}_{\sigma\sigma}\right) + f_{zz} \left(z_{x} \otimes z_{\sigma\sigma}\right) + F_{xu_{+}u_{+}}\left(I_{n_x} \otimes vec(\Sigma)\right) = 0$
$F_{xu_{+}u_{+}} = f_{y_{+}^{\ast\ast}} g_{xuu}^{\ast\ast} (g_x^{\ast} \otimes I_{n_u^2}) + f_{zz} \left( \left( z_{xu_{+}} \otimes z_{u_{+}} \right)P_{xu\_u} + \left(z_{x} \otimes z_{u_{+}u_{+}}\right) \right) + f_{zzz}\left(z_{x} \otimes z_{u_{+}} \otimes z_{u_{+}}\right)$

$(ghu\sigma\sigma) : A g_{u\sigma\sigma} + B g_{x\sigma\sigma} g^{*}_{u} + f_{y_{+}} g^{**}_{xx}\left(g^{*}_{u} \otimes g^{*}_{\sigma\sigma}\right) + f_{zz} \left(z_{u} \otimes z_{\sigma\sigma}\right) + F_{uu_{+}u_{+}}\left(I_{n_u} \otimes vec(\Sigma_u)\right) = 0$
$F_{uu_{+}u_{+}} = f_{y_{+}^{\ast\ast}} g_{xuu}^{\ast\ast} (g_u^{\ast} \otimes I_{n_u^2})  + f_{zz} \left( \left( z_{uu_{+}} \otimes z_{u_{+}} \right)P_{uu\_u} + \left(z_{u} \otimes z_{u_{+}u_{+}}\right) \right) + f_{zzz}\left(z_{u} \otimes z_{u_{+}} \otimes z_{u_{+}}\right)$

A and B are the common perturbation matrices:

$A = f_{y_0} + \begin{pmatrix} \underbrace{0}_{n\times n_{static}} &\vdots& \underbrace{f_{y^{**}_{+}} \cdot g^{**}_{x}}_{n \times n_{spred}} &\vdots& \underbrace{0}_{n\times n_{frwd}}  \end{pmatrix}$and $B = \begin{pmatrix} \underbrace{0}_{n \times n_{static}}&\vdots & \underbrace{0}_{n \times n_{pred}} & \vdots & \underbrace{f_{y^{**}_{+}}}_{n \times n_{sfwrd}} \end{pmatrix}$

and $z=(y_{-}^{\ast}; y; y_{+}^{\ast\ast}; u)$ denotes the dynamic model variables as in `M_.lead_lag_incidence`, $y^\ast$ denote state variables, $y^{\ast\ast}$ denote forward looking variables, $y_+$ denote the variables with a lead, $y_{-}$ denote variables with a lag, $y_0$ denote variables at period t, $f$ the model equations, and $f_z$ the first-order dynamic model derivatives, $f_{zz}$ the second-order dynamic derivatives, and $f_{zzz}$ the third-order dynamic model derivatives. Then:
$z_{x} = \begin{pmatrix}I\\g_{x}\\g^{**}_{x} g^{*}_{x}\\0\end{pmatrix}$, $z_{u} =\begin{pmatrix}0\\g_{u}\\g^{**}_{x} \cdot g^{*}_{u}\\I\end{pmatrix}$, $z_{u_{+}} =\begin{pmatrix}0\\0\\g^{**}_{u}\\0\end{pmatrix}$
$z_{xx} = \begin{pmatrix} 0\\g_{xx}\\g^{**}_{x} \left( g^{*}_x \otimes g^{*}_{x} \right) + g^{**}_{x} g^{*}_{x}\\0\end{pmatrix}$, $z_{xu} =\begin{pmatrix}0\\g_{xu}\\g^{**}_{xx} \left( g^{*}_x \otimes g^{*}_{u} \right) + g^{**}_{x} g^{*}_{xu}\\0\end{pmatrix}$, $z_{uu} =\begin{pmatrix}0\\g_{uu}\\g^{**}_{xx} \left( g^{*}_u \otimes g^{*}_{u} \right) + g^{**}_{x} g^{*}_{uu}\\0\end{pmatrix}$,
$z_{xu_{+}} =\begin{pmatrix}0\\0\\g^{**}_{xu} \left( g^{*}_x \otimes I \right)\\0\end{pmatrix}$, $z_{uu_{+}} =\begin{pmatrix}0\\0\\g^{**}_{xu} \left( g^{*}_{u} \otimes I \right)\\0\end{pmatrix}$, $z_{u_{+}u_{+}} =\begin{pmatrix}0\\0\\g^{\ast\ast}_{uu}\\0\end{pmatrix}$, $z_{\sigma\sigma} = \begin{pmatrix}0\\ g_{\sigma\sigma}\\ g^{\ast\ast}_{x}g^{\ast}_{\sigma\sigma} + g^{\ast\ast}_{\sigma\sigma}\\0 \end{pmatrix}$

$P$ are permutation matrices that can be computed using Matlab's `ipermute` function.

## Parameter derivatives of perturbation solutions
First, we need the parameter derivatives of first, second, third, and fourth derivatives of the dynamic model (i.e. g1,g2,g3,g4 in dynamic files), I make use of the implicit function theorem: Let $f_{z^k}$ denote the kth derivative (wrt all dynamic variables) of the dynamic model, then let $df_{z^k}$ denote the first-derivative (wrt all model parameters) of $f_{z^k}$ evaluated at the steady state. Note that $f_{z^k}$  is a function of both the model parameters $\theta$  and of the steady state of all dynamic variables $\bar{z}$, which also depend on the parameters. Hence, implicitly $f_{z^k}=f_{z^k}(\theta,\bar{z}(\theta))$  and $df_{z^k}$ consists of two parts:
1. direct derivative wrt to all model parameters given by the preprocessor in the `_params_derivs.m` files
2. contribution of derivative of steady state of dynamic variables (wrt all model parameters): $f_{z^{k+1}} \cdot d\bar{z}$
Note that we already have functionality to compute $d\bar{z}$ analytically.

Having this, the above perturbation systems are basically equations of the following types
$AX +BXC = RHS$ or $AX = RHS$
Now when taking the derivative (wrt to one single parameter $\theta_j$), we get
$A\mathrm{d}\{X\} + B\mathrm{d}\{X\}C = \mathrm{d}\{RHS\} - \mathrm{d}\{A\}X -  \mathrm{d}\{B\}XC - BX\mathrm{d}\{C\}$
or
$A\mathrm{d}\{X\}  = \mathrm{d}\{RHS\} - \mathrm{d}\{A\}X$
The first one is a Sylvester type equation, the second one can be solved by taking the inverse of $A$. The only diffculty and tedious work arrises in computing (the highly sparse) derivatives of $RHS$.

***

# New functions: `
## get_perturbation_params_derivs.m`and `get_perturbation_params_derivs_numerical_objective.m`
* The parameter derivatives up to third order are computed in the new function`get_perturbation_params_derivs.m` both analytically and numerically. For numerical derivatives `get_perturbation_params_derivs_numerical_objective.m` is the objective for `fjaco.m` or `hessian_sparse.m` or `hessian.m`.
* `get_perturbation_params_derivs.m` is basically an extended version of the previous `get_first_order_solution_params_derivs.m` function.
* * `get_perturbation_params_derivs_numerical_objective.m`builds upon `identification_numerical_objective.m`. It is used for numerical derivatives, whenever `analytic_derivation_mode=-1|-2`. It takes from `identification_numerical_objective.m` the parts that compute numerical parameter Jacobians of steady state, dynamic model equations, and perturbation solution matrices. Hence, these parts are removed in `identification_numerical_objective.m` and it only computes numerical parameter Jacobian of moments and spectrum which are needed for identification analysis in `get_identification_jacobians.m`, when `analytic_derivation_mode=-1` only.
* Detailed changes:
      * Most important: notation of this function is now in accordance to the k_order_solver, i.e. we do not compute derivatives of Kalman transition matrices A and B, but rather the solution matrices ghx,ghu,ghxx,ghxu,ghuu,ghs2,ghxxx,ghxxu,ghxuu,ghuuu,ghxss,ghuss in accordance with notation used in `oo_.dr`. As a byproduct at first-order, focusing on ghx and ghu instead of Kalman transition matrices A and B makes the computations slightly faster for large models (e.g. for Quest the computations were faster by a couple of seconds, not much, but okay).
      * Removed use of `kstate`, see also Dynare/dynare#1653 and Dynare/dynare!1656
      * Output arguments are stored in a structure `DERIVS`, there is also a flag `d2flag` that computes parameter hessians needed only in `dsge_likelihood.m`.
      * Removed `kronflag` as input. `options_.analytic_derivation_mode` is now used instead of `kronflag`.
      * Removed `indvar`, an index that was used to selected specific variables in the derivatives. This is not needed, as we always compute the parameter derivatives for all variables first and then select a subset of variables. The selection now takes place in other functions, like `dsge_likelihood.m`.
      * Introduced some checks: (i) deterministic exogenous variables are not supported, (ii) Kronecker method only compatible with first-order approximation so reset to sylvester method, (iii) for purely backward or forward models we need to be careful with the rows in `M_.lead_la	g_incidence`, (iv) if `_params_derivs.m` files are missing an error is thrown.
      * For numerical derivatives, if mod file does not contain an `estimated_params_block`, a temporary one with the most important parameter information is created.
## `unfold_g4.m`
* When evaluating g3 and g4 one needs to take into account that these do not contain symmetric elements, so one needs to use `unfold_g3.m` and the new function `unfold_g4.m`. This returns an unfolded version of the same matrix (i.e. with symmetric elements).

***

# New test models
`.gitignore` and `Makefile.am` are changed accordingly. Also now it is possible to run test suite on analytic_derivatives, i.e. run `make check m/analytic_derivatives`

## `analytic_derivatives/BrockMirman_PertParamsDerivs.mod`
* This is the Brock Mirman model, where we know the exact policy function $g$ for capital and consumption. As this does not imply a nonzero $g_{\sigma\sigma}$, $g_{x\sigma\sigma}$, $g_{u\sigma\sigma}$ I added some artificial equations to get nonzero solution matrices with respect to $\sigma$. The true perturbation solution matrices  $g_x$ , $g_u$, $g_{xx}$, $g_{xu}$, $g_{uu}$, $g_{\sigma\sigma}$, $g_{xxx}$, $g_{xxu}$, $g_{xuu}$, $g_{uuu}$, $g_{x\sigma\sigma}$, $g_{u\sigma\sigma}$ are then computed analytically with Matlab's symbolic toolbox and saved in `nBrockMirmanSYM.mat`. There is a preprocessor flag that recreates these analytical computations if changes are needed (and to check whether I made some errors here ;-) )
* Then solution matrices up to third order and their parameter Jacobians are then compared to the ones computed by Dynare's `k_order_solver` and by `get_perturbation_params_derivs` for all `analytic_derivation_mode`'s. There will be an error if the maximum absolute deviation is too large, i.e. for numerical derivatives (`analytic_derivation_mode=-1|-2`) the tolerance is choosen lower (around 1e-5); for analytical methods we are stricter: around 1e-13 for first-order,  1e-12 for second order, and 1e-11 for third-order.
* As a side note, this mod file also checks Dynare's `k_order_solver` algorithm and throws an error if something is wrong.
* This test model shows that the new functionality works well. And analytical derivatives perform way better and accurate than numerical ones, even for this small model.
## `analytic_derivatives/burnside_3_order_PertParamsDerivs.mod`
* This builds upon `tests/k_order_perturbation/burnside_k_order.mod` and computes the true parameter derivatives analytically by hand.
      * This test model also shows that the new functionality works well.

## `analytic_derivatives/LindeTrabandt2019.mod`
* Shows that the new functionality also works for medium-sized models, i.e. a SW type model solved at third order with 35 variables (11 states). 2 shocks and 20 parameters.
* This mod file can be used to tweak the speed of the computations in the future.
* Compares numerical versus analytical parameter derivatives (for first, second and third order). Note that this model clearly shows that numerical ones are quite different than analytical ones even at first order!
## `identification/LindeTrabandt2019_xfail.mod`
* This model is a check for issue Dynare/dynare#1595, see fjaco.m below, and will fail.
* Removed `analytic_derivatives/ls2003.mod` as this mod file is neither in the testsuite nor does it work.

***

# Detailed changes in other functions
## `get_first_order_solution_params_derivs.m`
* Deleted, or actually, renamed to `get_perturbation_params_derivs.m`, as this function now is able to compute the derivatives up to third order

## `identification_numerical_objective.m`
* `get_perturbation_params_derivs_numerical_objective.m`builds upon `identification_numerical_objective.m`. It takes from `identification_numerical_objective.m` the parts that compute numerical parameter Jacobians of steady state, dynamic model equations, and perturbation solution matrices. Hence, these parts are removed in `identification_numerical_objective.m` and it only computes numerical parameter Jacobian of moments and spectrum which are needed for identification analysis in `get_identification_jacobians.m`, when `analytic_derivation_mode=-1` only.

## `dsge_likelihood.m`
* As `get_first_order_solution_params_derivs.m`is renamed to `get_perturbation_params_derivs.m`, the call is adapted. That is,`get_perturbation_params_derivs` does not compute the derivatives of the Kalman transition `T`matrix anymore, but instead of the dynare solution matrix `ghx`. So we recreate `T` here as this amounts to adding some zeros and focusing on selected variables only.
* Added some checks to make sure the first-order approximation is selected.
* Removed `kron_flag` as input, as `get_perturbation_params_derivs` looks into `options_.analytic_derivation_mode` for `kron_flag`.

## `dynare_identification.m`
* make sure that setting `analytic_derivation_mode` is set both in `options_ident` and `options_`. Note that at the end of the function we restore the `options_` structure, so all changes are local. In a next merge request, I will remove the global variables to make all variables local.

## `get_identification_jacobians.m`
* As `get_first_order_solution_params_derivs.m`is renamed to `get_perturbation_params_derivs.m`, the call is adapted. That is,`get_perturbation_params_derivs` does not compute the derivatives of the Kalman transition `A` and `B` matrix anymore, but instead of the dynare solution matrix `ghx` and `ghu`. So we recreate these matrices here instead of in `get_perturbation_params_derivs.m`.
* Added `str2func` for better function handles in `fjaco.m`.

## `fjaco.m`
* make `tol`an option, which can be adjusted by changing `options_.dynatol.x`for identification and parameter derivatives purposes.
* include a check and an informative error message, if numerical derivatives (two-sided finite difference method) yield errors in `resol.m` for identification and parameter derivatives purposes. This closes issue  Dynare/dynare#1595.
* Changed year of copyright to 2010-2017,2019

***

# Further suggestions and questions
* Ones this is merged, I will merge request an improvement of the identification toolbox, which will work up to third order using the pruned state space. This will also remove some issues and bugs, and also I will remove global variables in this request.
* The third-order derivatives can be further improved by taking sparsity into account and use mex versions for kronecker products etc. I leave this for further testing (and if anybody actually uses this ;-) )
time-shift
Willi Mutschler 2019-12-17 18:17:09 +00:00 committed by Sébastien Villemot
parent d5f90c8376
commit 5a8c206760
18 changed files with 2856 additions and 1083 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

28
matlab/dsge_likelihood.m
View File

@ -111,7 +111,7 @@ 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{get_first_order_solution_params_deriv}, @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_perturbation_params_deriv}, @ref{kalman_filter}, @ref{score}, @ref{AHessian}, @ref{missing_observations_kalman_filter}, @ref{univariate_kalman_filter}, @ref{priordens}
%! @end deftypefn
%@eod:
@ -522,11 +522,29 @@ if analytic_derivation
else
indparam=[];
end
if full_Hess
[DT, ~, ~, DOm, DYss, ~, D2T, D2Om, D2Yss] = get_first_order_solution_params_deriv(A, B, EstimatedParameters, Model,DynareResults,DynareOptions,kron_flag,indparam,indexo,[],iv);
else
[DT, ~, ~, DOm, DYss] = get_first_order_solution_params_deriv(A, B, EstimatedParameters, Model,DynareResults,DynareOptions,kron_flag,indparam,indexo,[],iv);
old_order = DynareOptions.order;
if DynareOptions.order > 1%not sure whether this check is necessary
DynareOptions.order = 1; fprintf('Reset order to 1 for analytical parameter derivatives.\n');
end
old_analytic_derivation_mode = DynareOptions.analytic_derivation_mode;
DynareOptions.analytic_derivation_mode = kron_flag;
if full_Hess
DERIVS = get_perturbation_params_derivs(Model, DynareOptions, EstimatedParameters, DynareResults, indparam, indexo, [], true);
indD2T = reshape(1:Model.endo_nbr^2, Model.endo_nbr, Model.endo_nbr);
indD2Om = dyn_unvech(1:Model.endo_nbr*(Model.endo_nbr+1)/2);
D2T = DERIVS.d2KalmanA(indD2T(iv,iv),:);
D2Om = DERIVS.d2Om(dyn_vech(indD2Om(iv,iv)),:);
D2Yss = DERIVS.d2Yss(iv,:,:);
else
DERIVS = get_perturbation_params_derivs(Model, DynareOptions, EstimatedParameters, DynareResults, indparam, indexo, [], false);
end
DT = zeros(Model.endo_nbr, Model.endo_nbr, size(DERIVS.dghx,3));
DT(:,Model.nstatic+(1:Model.nspred),:) = DERIVS.dghx;
DT = DT(iv,iv,:);
DOm = DERIVS.dOm(iv,iv,:);
DYss = DERIVS.dYss(iv,:);
DynareOptions.order = old_order; %make sure order is reset (not sure if necessary)
DynareOptions.analytic_derivation_mode = old_analytic_derivation_mode;%make sure analytic_derivation_mode is reset (not sure if necessary)
else
DT = derivatives_info.DT(iv,iv,:);
DOm = derivatives_info.DOm(iv,iv,:);

7
matlab/dynare_identification.m
View File

@ -284,7 +284,6 @@ options_ = set_default_option(options_,'datafile','');
options_.mode_compute = 0;
options_.plot_priors = 0;
options_.smoother = 1;
options_.options_ident = options_ident; %store identification options into global options
[dataset_, dataset_info, xparam1, hh, M_, options_, oo_, estim_params_, bayestopt_, bounds] = dynare_estimation_init(M_.endo_names, fname, 1, M_, options_, oo_, estim_params_, bayestopt_);
%outputting dataset_ is needed for Octave
@ -294,6 +293,7 @@ options_ident = set_default_option(options_ident,'analytic_derivation_mode', opt
% 1: kronecker products method to compute analytical derivatives as in Iskrev (2010)
% -1: numerical two-sided finite difference method to compute numerical derivatives of all Jacobians using function identification_numerical_objective.m (previously thet2tau.m)
% -2: numerical two-sided finite difference method to compute numerically dYss, dg1, d2Yss and d2g1, the Jacobians are then computed analytically as in options 0
options_.analytic_derivation_mode = options_ident.analytic_derivation_mode; %overwrite setting
% initialize persistent variables in prior_draw
if prior_exist
@ -385,12 +385,13 @@ options_ident = set_default_option(options_ident,'max_dim_subsets_groups',min([4
% In identification_checks_via_subsets.m, when checks_via_subsets=1, this
% option sets the maximum dimension of groups of parameters for which
% the corresponding rank criteria is checked
options_.options_ident = options_ident; %store identification options into global options
MaxNumberOfBytes = options_.MaxNumberOfBytes; %threshold when to save from memory to files
store_options_ident = options_ident;
iload = options_ident.load_ident_files;
SampleSize = options_ident.prior_mc;
if iload <=0
%% Perform new identification analysis, i.e. do not load previous analysis
if prior_exist

33
matlab/fjaco.m
View File

@ -10,7 +10,7 @@ function fjac = fjaco(f,x,varargin)
% OUTPUT
% fjac : finite difference Jacobian
%
% Copyright (C) 2010-2017 Dynare Team
% Copyright (C) 2010-2017,2019 Dynare Team
%
% This file is part of Dynare.
%
@ -29,7 +29,10 @@ function fjac = fjaco(f,x,varargin)
ff=feval(f,x,varargin{:});
tol = eps.^(1/3);
tol = eps.^(1/3); %some default value
if strcmp(func2str(f),'get_perturbation_params_derivs_numerical_objective') || strcmp(func2str(f),'identification_numerical_objective')
tol= varargin{5}.dynatol.x;
end
h = tol.*max(abs(x),1);
xh1=x+h; xh0=x-h;
h=xh1-xh0;
@ -37,8 +40,34 @@ fjac = NaN(length(ff),length(x));
for j=1:length(x)
xx = x;
xx(j) = xh1(j); f1=feval(f,xx,varargin{:});
if isempty(f1) && (strcmp(func2str(f),'get_perturbation_params_derivs_numerical_objective') || strcmp(func2str(f),'identification_numerical_objective') )
[~,info]=feval(f,xx,varargin{:});
disp_info_error_identification_perturbation(info,j);
end
xx(j) = xh0(j); f0=feval(f,xx,varargin{:});
if isempty(f0) && (strcmp(func2str(f),'get_perturbation_params_derivs_numerical_objective') || strcmp(func2str(f),'identification_numerical_objective') )
[~,info]=feval(f,xx,varargin{:});
disp_info_error_identification_perturbation(info,j)
end
fjac(:,j) = (f1-f0)/h(j);
end
feval(f,x,varargin{:});
%Auxiliary functions
function disp_info_error_identification_perturbation(info,j)
% there are errors in the solution algorithm
probl_par = get_the_name(j,varargin{5}.TeX,varargin{3},varargin{2},varargin{5});
skipline()
message = get_error_message(info,0,varargin{5});
fprintf('Parameter error in numerical two-sided difference method:\n')
fprintf('Cannot solve the model for %s (info = %d, %s)\n', probl_par, info(1), message);
fprintf('Possible solutions:\n')
fprintf(' -- check your mod file, calibration and steady state computations carefully\n');
fprintf(' -- use analytic derivatives, i.e. set analytic_derivation_mode=0\n');
fprintf(' -- use an estimated_params block without %s or change its value\n', probl_par);
fprintf(' -- change numerical tolerance level in fjaco.m (you can tune ''options_.dynatol.x'' or change fjaco.m function directly)\n');
error('fjaco.m: numerical two-sided difference method yields errors in solution algorithm');
end
end %main function end

959
matlab/get_first_order_solution_params_deriv.m
View File

@ -1,959 +0,0 @@
function [dA, dB, dSigma_e, dOm, dYss, dg1, d2A, d2Om, d2Yss] = get_first_order_solution_params_deriv(A, B, estim_params, M, oo, options, kronflag, indpmodel, indpstderr, indpcorr, indvar)
%[dA, dB, dSigma_e, dOm, dYss, dg1, d2A, d2Om, d2Yss] = get_first_order_solution_params_deriv(A, B, estim_params, M, oo, options, kronflag, indpmodel, indpstderr, indpcorr, indvar)
% previously getH.m
% -------------------------------------------------------------------------
% Computes first and second derivatives (with respect to parameters) of
% (1) reduced-form solution (dA, dB, dSigma_e, dOm, d2A, d2Om)
% (1) steady-state (dYss, d2Yss)
% (3) Jacobian (wrt to dynamic variables) of dynamic model (dg1)
% Note that the order in the parameter Jacobians is the following:
% first stderr parameters, second corr parameters, third model parameters
% =========================================================================
% INPUTS
% A: [endo_nbr by endo_nbr] Transition matrix from Kalman filter
% for all endogenous declared variables, in DR order
% B: [endo_nbr by exo_nbr] Transition matrix from Kalman filter
% mapping shocks today to endogenous variables today, in DR order
% estim_params: [structure] storing the estimation information
% M: [structure] storing the model information
% oo: [structure] storing the reduced-form solution results
% options: [structure] storing the options
% kronflag: [scalar] method to compute Jacobians (equal to analytic_derivation_mode in options_ident). Default:0
% * 0: efficient sylvester equation method to compute
% analytical derivatives as in Ratto & Iskrev (2011)
% * 1: kronecker products method to compute analytical
% derivatives as in Iskrev (2010)
% * -1: numerical two-sided finite difference method to
% compute numerical derivatives of all output arguments
% using function identification_numerical_objective.m
% (previously thet2tau.m)
% * -2: numerical two-sided finite difference method to
% compute numerically dYss, dg1, d2Yss and d2g1, the other
% output arguments are computed analytically as in kronflag=0
% indpmodel: [modparam_nbr by 1] index of estimated parameters in M_.params;
% corresponds to model parameters (no stderr and no corr)
% in estimated_params block; if estimated_params block is
% not available, then all model parameters are selected
% indpstderr: [stderrparam_nbr by 1] index of estimated standard errors,
% i.e. for all exogenous variables where "stderr" is given
% in the estimated_params block; if estimated_params block
% is not available, then all stderr parameters are selected
% indpcorr: [corrparam_nbr by 2] matrix of estimated correlations,
% i.e. for all exogenous variables where "corr" is given
% in the estimated_params block; if estimated_params block
% is not available, then no corr parameters are selected
% indvar: [var_nbr by 1] index of considered (or observed) variables
% -------------------------------------------------------------------------
% OUTPUTS
% dA: [var_nbr by var_nbr by totparam_nbr] in DR order
% Jacobian (wrt to all parameters) of transition matrix A
% dB: [var_nbr by exo_nbr by totparam_nbr] in DR order
% Jacobian (wrt to all parameters) of transition matrix B
% dSigma_e: [exo_nbr by exo_nbr by totparam_nbr] in declaration order
% Jacobian (wrt to all paramters) of M_.Sigma_e
% dOm: [var_nbr by var_nbr by totparam_nbr] in DR order
% Jacobian (wrt to all paramters) of Om = (B*M_.Sigma_e*B')
% dYss: [var_nbr by modparam_nbr] in DR order
% Jacobian (wrt model parameters only) of steady state
% dg1: [endo_nbr by (dynamicvar_nbr + exo_nbr) by modparam_nbr] in DR order
% Jacobian (wrt to model parameters only) of Jacobian of dynamic model
% d2A: [var_nbr*var_nbr by totparam_nbr*(totparam_nbr+1)/2] in DR order
% Unique entries of Hessian (wrt all parameters) of transition matrix A
% d2Om: [var_nbr*(var_nbr+1)/2 by totparam_nbr*(totparam_nbr+1)/2] in DR order
% Unique entries of Hessian (wrt all parameters) of Omega
% d2Yss: [var_nbr by modparam_nbr by modparam_nbr] in DR order
% Unique entries of Hessian (wrt model parameters only) of steady state
% -------------------------------------------------------------------------
% This function is called by
% * dsge_likelihood.m
% * get_identification_jacobians.m (previously getJJ.m)
% -------------------------------------------------------------------------
% This function calls
% * [fname,'.dynamic']
% * [fname,'.dynamic_params_derivs']
% * [fname,'.static']
% * [fname,'.static_params_derivs']
% * commutation
% * dyn_vech
% * dyn_unvech
% * fjaco
% * get_2nd_deriv (embedded)
% * get_2nd_deriv_mat(embedded)
% * get_all_parameters
% * get_all_resid_2nd_derivs (embedded)
% * get_hess_deriv (embedded)
% * hessian_sparse
% * sylvester3
% * sylvester3a
% * identification_numerical_objective.m (previously thet2tau.m)
% =========================================================================
% Copyright (C) 2010-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/>.
% =========================================================================
fname = M.fname;
dname = M.dname;
maximum_exo_lag = M.maximum_exo_lag;
maximum_exo_lead = M.maximum_exo_lead;
maximum_endo_lag = M.maximum_endo_lag;
maximum_endo_lead = M.maximum_endo_lead;
lead_lag_incidence = M.lead_lag_incidence;
[I,~] = find(lead_lag_incidence'); %I is used to select nonzero columns of the Jacobian of endogenous variables in dynamic model files
ys = oo.dr.ys; %steady state of endogenous variables in declaration order
yy0 = oo.dr.ys(I); %steady state of dynamic (endogenous and auxiliary variables) in DR order
ex0 = oo.exo_steady_state'; %steady state of exogenous variables in declaration order
params0 = M.params; %values at which to evaluate dynamic, static and param_derivs files
Sigma_e0 = M.Sigma_e; %covariance matrix of exogenous shocks
Corr_e0 = M.Correlation_matrix; %correlation matrix of exogenous shocks
stderr_e0 = sqrt(diag(Sigma_e0)); %standard errors of exogenous shocks
param_nbr = M.param_nbr; %number of all declared model parameters in mod file
modparam_nbr = length(indpmodel); %number of model parameters to be used
stderrparam_nbr = length(indpstderr); %number of stderr parameters to be used
corrparam_nbr = size(indpcorr,1); %number of stderr parameters to be used
totparam_nbr = modparam_nbr + stderrparam_nbr + corrparam_nbr; %total number of parameters to be used
if nargout > 6
modparam_nbr2 = modparam_nbr*(modparam_nbr+1)/2; %number of unique entries of model parameters only in second-order derivative matrix
totparam_nbr2 = totparam_nbr*(totparam_nbr+1)/2; %number of unique entries of all parameters in second-order derivative matrix
%get indices of elements in second derivatives of parameters
indp2tottot = reshape(1:totparam_nbr^2,totparam_nbr,totparam_nbr);
indp2stderrstderr = indp2tottot(1:stderrparam_nbr , 1:stderrparam_nbr);
indp2stderrcorr = indp2tottot(1:stderrparam_nbr , stderrparam_nbr+1:stderrparam_nbr+corrparam_nbr);
indp2modmod = indp2tottot(stderrparam_nbr+corrparam_nbr+1:stderrparam_nbr+corrparam_nbr+modparam_nbr , stderrparam_nbr+corrparam_nbr+1:stderrparam_nbr+corrparam_nbr+modparam_nbr);
if totparam_nbr ~=1
indp2tottot2 = dyn_vech(indp2tottot); %index of unique second-order derivatives
else
indp2tottot2 = indp2tottot;
end
if modparam_nbr ~= 1
indp2modmod2 = dyn_vech(indp2modmod); %get rid of cross derivatives
else
indp2modmod2 = indp2modmod;
end
end
endo_nbr = size(A,1); %number of all declared endogenous variables
var_nbr = length(indvar); %number of considered variables
exo_nbr = size(B,2); %number of exogenous shocks in model
if kronflag == -1
% numerical two-sided finite difference method using function identification_numerical_objective.m (previously thet2tau.m) for Jacobian (wrt parameters) of A, B, Sig, Om, Yss, and g1
para0 = get_all_parameters(estim_params, M); %get all selected parameters in estimated_params block, stderr and corr come first, then model parameters
if isempty(para0)
%if there is no estimated_params block, consider all stderr and all model parameters, but no corr parameters
para0 = [stderr_e0', params0'];
end
%Jacobians (wrt paramters) of steady state, solution matrices A and B, as well as Sigma_e for ALL variables [outputflag=0]
dYssABSige = fjaco('identification_numerical_objective', para0, 0, estim_params, M, oo, options, indpmodel, indpstderr, indpcorr, indvar);
M.params = params0; %make sure values are set back
M.Sigma_e = Sigma_e0; %make sure values are set back
M.Correlation_matrix = Corr_e0 ; %make sure values are set back
% get Jacobians for Yss, A, and B from dYssABSige
indYss = 1:var_nbr;
indA = (var_nbr+1):(var_nbr+var_nbr^2);
indB = (var_nbr+var_nbr^2+1):(var_nbr+var_nbr^2+var_nbr*exo_nbr);
indSigma_e = (var_nbr+var_nbr^2+var_nbr*exo_nbr+1):(var_nbr+var_nbr^2+var_nbr*exo_nbr+exo_nbr*(exo_nbr+1)/2);
dYss = dYssABSige(indYss , stderrparam_nbr+corrparam_nbr+1:end); %in tensor notation only wrt model parameters
dA = reshape(dYssABSige(indA , :) , [var_nbr var_nbr totparam_nbr]); %in tensor notation
dB = reshape(dYssABSige(indB , :) , [var_nbr exo_nbr totparam_nbr]); %in tensor notation
dOm = zeros(var_nbr,var_nbr,totparam_nbr); %initialize in tensor notation
dSigma_e = zeros(exo_nbr,exo_nbr,totparam_nbr); %initialize in tensor notation
% get Jacobians of Sigma_e and Om wrt stderr parameters
if ~isempty(indpstderr)
for jp=1:stderrparam_nbr
dSigma_e(:,:,jp) = dyn_unvech(dYssABSige(indSigma_e , jp));
dOm(:,:,jp) = B*dSigma_e(:,:,jp)*B'; %note that derivatives of B wrt stderr parameters are zero by construction
end
end
% get Jacobians of Sigma_e and Om wrt corr parameters
if ~isempty(indpcorr)
for jp=1:corrparam_nbr
dSigma_e(:,:,stderrparam_nbr+jp) = dyn_unvech(dYssABSige(indSigma_e , stderrparam_nbr+jp));
dOm(:,:,stderrparam_nbr+jp) = B*dSigma_e(:,:,stderrparam_nbr+jp)*B'; %note that derivatives of B wrt corr parameters are zero by construction
end
end
% get Jacobian of Om wrt model parameters
if ~isempty(indpmodel)
for jp=1:modparam_nbr
dOm(:,:,stderrparam_nbr+corrparam_nbr+jp) = dB(:,:,stderrparam_nbr+corrparam_nbr+jp)*Sigma_e0*B' + B*Sigma_e0*dB(:,:,stderrparam_nbr+corrparam_nbr+jp)'; %note that derivatives of Sigma_e wrt model parameters are zero by construction
end
end
%Jacobian (wrt model parameters ONLY) of steady state and of Jacobian of all dynamic model equations [outputflag=-1]
dYssg1 = fjaco('identification_numerical_objective', params0(indpmodel), -1, estim_params, M, oo, options, indpmodel, [], [], (1:endo_nbr)');
M.params = params0; %make sure values are set back
M.Sigma_e = Sigma_e0; %make sure values are set back
M.Correlation_matrix = Corr_e0 ; %make sure values are set back
dg1 = reshape(dYssg1(endo_nbr+1:end,:),[endo_nbr, length(yy0)+length(ex0), modparam_nbr]); %get rid of steady state and in tensor notation
if nargout > 6
%Hessian (wrt paramters) of steady state, solution matrices A and Om [outputflag=-2]
% note that hessian_sparse does not take symmetry into account, i.e. compare hessian_sparse.m to hessian.m, but focuses already on unique values, which are duplicated below
d2YssAOm = hessian_sparse('identification_numerical_objective', para0', options.gstep, -2, estim_params, M, oo, options, indpmodel, indpstderr, indpcorr, indvar);
M.params = params0; %make sure values are set back
M.Sigma_e = Sigma_e0; %make sure values are set back
M.Correlation_matrix = Corr_e0 ; %make sure values are set back
d2A = d2YssAOm(indA , indp2tottot2); %only unique elements
d2Om = d2YssAOm(indA(end)+1:end , indp2tottot2); %only unique elements
d2Yss = zeros(var_nbr,modparam_nbr,modparam_nbr); %initialize
for j = 1:var_nbr
d2Yss(j,:,:) = dyn_unvech(full(d2YssAOm(j,indp2modmod2))); %full Hessian for d2Yss, note that here we duplicate unique values for model parameters
end
clear d2YssAOm
end
return %[END OF MAIN FUNCTION]!!!!!
end
if kronflag == -2
% numerical two-sided finite difference method to compute numerically
% dYss, dg1, d2Yss and d2g1, the rest is computed analytically (kronflag=0) below
modpara0 = params0(indpmodel); %focus only on model parameters for dYss, d2Yss and dg1
[~, g1] = feval([fname,'.dynamic'], yy0, ex0, params0, ys, 1);
%g1 is [endo_nbr by (dynamicvar_nbr+exo_nbr)] first derivative (wrt all endogenous, exogenous and auxiliary variables) of dynamic model equations, i.e. df/d[yy0;ex0], in DR order
if nargout > 6
% computation of d2Yss and d2g1, i.e. second derivative (wrt. parameters) of Jacobian (wrt endogenous and auxilary variables) of dynamic model [outputflag = -1]
% note that hessian_sparse does not take symmetry into account, i.e. compare hessian_sparse.m to hessian.m, but focuses already on unique values, which are duplicated below
d2Yssg1 = hessian_sparse('identification_numerical_objective', modpara0, options.gstep, -1, estim_params, M, oo, options, indpmodel, [], [], (1:endo_nbr)');
M.params = params0; %make sure values are set back
M.Sigma_e = Sigma_e0; %make sure values are set back
M.Correlation_matrix = Corr_e0 ; %make sure values are set back
d2Yss = reshape(full(d2Yssg1(1:endo_nbr,:)), [endo_nbr modparam_nbr modparam_nbr]); %put into tensor notation
for j=1:endo_nbr
d2Yss(j,:,:) = dyn_unvech(dyn_vech(d2Yss(j,:,:))); %add duplicate values to full hessian
end
d2g1_full = d2Yssg1(endo_nbr+1:end,:);
%store only nonzero unique entries and the corresponding indices of d2g1:
% rows: respective derivative term
% 1st column: equation number of the term appearing
% 2nd column: column number of variable in Jacobian of the dynamic model
% 3rd column: number of the first parameter in derivative
% 4th column: number of the second parameter in derivative
% 5th column: value of the Hessian term
ind_d2g1 = find(d2g1_full);
d2g1 = zeros(length(ind_d2g1),5);
for j=1:length(ind_d2g1)
[i1, i2] = ind2sub(size(d2g1_full),ind_d2g1(j));
[ig1, ig2] = ind2sub(size(g1),i1);
[ip1, ip2] = ind2sub([modparam_nbr modparam_nbr],i2);
d2g1(j,:) = [ig1 ig2 ip1 ip2 d2g1_full(ind_d2g1(j))];
end
clear d2g1_full;
end
%Jacobian (wrt parameters) of steady state and Jacobian of dynamic model equations [outputflag=-1]
dg1 = fjaco('identification_numerical_objective', modpara0, -1, estim_params, M, oo, options, indpmodel, [], [], (1:endo_nbr)');
M.params = params0; %make sure values are set back
M.Sigma_e = Sigma_e0; %make sure values are set back
M.Correlation_matrix = Corr_e0 ; %make sure values are set back
dYss = dg1(1:endo_nbr , :);
dg1 = reshape(dg1(endo_nbr+1 : end , :),[endo_nbr, length(yy0)+length(ex0), modparam_nbr]); %get rid of steady state
elseif (kronflag == 0 || kronflag == 1)
% Analytical method to compute dYss, dg1, d2Yss and d2g1
[~, g1_static] = feval([fname,'.static'], ys, ex0, params0);
%g1_static is [endo_nbr by endo_nbr] first-derivative (wrt variables) of static model equations f, i.e. df/dys, in declaration order
rp_static = feval([fname,'.static_params_derivs'], ys, repmat(ex0, maximum_exo_lag+maximum_exo_lead+1), params0);
%rp_static is [endo_nbr by param_nbr] first-derivative (wrt parameters) of static model equations f, i.e. df/dparams, in declaration order
dys = -g1_static\rp_static;
%use implicit function theorem (equation 5 of Ratto and Iskrev (2011) to compute [endo_nbr by param_nbr] first-derivative (wrt parameters) of steady state analytically, note that dys is in declaration order
d2ys = zeros(length(ys), param_nbr, param_nbr); %initialize in tensor notation
if nargout > 6
[~, ~, g2_static] = feval([fname,'.static'], ys, ex0, params0);
%g2_static is [endo_nbr by endo_nbr^2] second derivative (wrt variables) of static model equations f, i.e. d(df/dys)/dys, in declaration order
[~, g1, g2, g3] = feval([fname,'.dynamic'], yy0, ex0, params0, ys, 1);
%g1 is [endo_nbr by (dynamicvar_nbr+exo_nbr)] first derivative (wrt all endogenous, exogenous and auxiliary variables) of dynamic model equations, i.e. df/d[yy0;ex0], in DR order
%g2 is [endo_nbr by (dynamicvar_nbr+exo_nbr)^2] second derivative (wrt all endogenous, exogenous and auxiliary variables) of dynamic model equations, i.e. d(df/d[yy0;ex0])/d[yy0;ex0], in DR order
%g3 is [endo_nbr by (dynamicvar_nbr+exo_nbr)^2] third-derivative (wrt all endogenous, exogenous and auxiliary variables) of dynamic model equations, i.e. d(df/d[yy0;ex0])/d[yy0;ex0], in DR order
[~, gp_static, rpp_static] = feval([fname,'.static_params_derivs'], ys, ex0, params0);
%gp_static is [endo_nbr by endo_nbr by param_nbr] first derivative (wrt parameters) of first-derivative (wrt variables) of static model equations f, i.e. (df/dys)/dparams, in declaration order
%rpp_static are nonzero values and corresponding indices of second derivative (wrt parameters) of static model equations f, i.e. d(df/dparams)/dparams, in declaration order
rpp_static = get_all_resid_2nd_derivs(rpp_static, length(ys), param_nbr); %make full matrix out of nonzero values and corresponding indices
%rpp_static is [endo_nbr by param_nbr by param_nbr] second derivative (wrt parameters) of static model equations, i.e. d(df/dparams)/dparams, in declaration order
if isempty(find(g2_static))
%auxiliary expression on page 8 of Ratto and Iskrev (2011) is zero, i.e. gam = 0
for j = 1:param_nbr
%using the implicit function theorem, equation 15 on page 7 of Ratto and Iskrev (2011)
d2ys(:,:,j) = -g1_static\rpp_static(:,:,j);
%d2ys is [endo_nbr by param_nbr by param_nbr] second-derivative (wrt parameters) of steady state, i.e. d(dys/dparams)/dparams, in declaration order
end
else
gam = rpp_static*0; %initialize auxiliary expression on page 8 of Ratto and Iskrev (2011)
for j = 1:endo_nbr
tmp_gp_static_dys = (squeeze(gp_static(j,:,:))'*dys);
gam(j,:,:) = transpose(reshape(g2_static(j,:),[endo_nbr endo_nbr])*dys)*dys + tmp_gp_static_dys + tmp_gp_static_dys';
end
for j = 1:param_nbr
%using the implicit function theorem, equation 15 on page 7 of Ratto and Iskrev (2011)
d2ys(:,:,j) = -g1_static\(rpp_static(:,:,j)+gam(:,:,j));
%d2ys is [endo_nbr by param_nbr by param_nbr] second-derivative (wrt parameters) of steady state, i.e. d(dys/dparams)/dparams, in declaration order
end
clear gp_static g2_static tmp_gp_static_dys gam
end
end
%handling of steady state for nonstationary variables
if any(any(isnan(dys)))
[U,T] = schur(g1_static);
qz_criterium = options.qz_criterium;
e1 = abs(ordeig(T)) < qz_criterium-1;
k = sum(e1); % Number of non stationary variables.
% Number of stationary variables: n = length(e1)-k
[U,T] = ordschur(U,T,e1);
T = T(k+1:end,k+1:end);
%using implicit function theorem, equation 5 of Ratto and Iskrev (2011), in declaration order
dys = -U(:,k+1:end)*(T\U(:,k+1:end)')*rp_static;
if nargout > 6
disp('Computation of d2ys for nonstationary variables is not yet correctly handled if g2_static is nonempty, but continue anyways...')
for j = 1:param_nbr
%using implicit function theorem, equation 15 of Ratto and Iskrev (2011), in declaration order
d2ys(:,:,j) = -U(:,k+1:end)*(T\U(:,k+1:end)')*rpp_static(:,:,j); %THIS IS NOT CORRECT, IF g2_static IS NONEMPTY. WE NEED TO ADD GAM [willi]
end
end
end
if nargout > 6
[~, gp, ~, gpp, hp] = feval([fname,'.dynamic_params_derivs'], yy0, ex0, params0, ys, 1, dys, d2ys);
%gp is [endo_nbr by (dynamicvar_nbr + exo_nbr) by param_nbr] first-derivative (wrt parameters) of first-derivative (wrt all endogenous, auxiliary and exogenous variables) of dynamic model equations, i.e. d(df/dvars)/dparam, in DR order
%gpp are nonzero values and corresponding indices of second-derivative (wrt parameters) of first-derivative (wrt all endogenous, auxiliary and exogenous variables) of dynamic model equations, i.e. d(d(df/dvars)/dparam)/dparam, in DR order
%hp are nonzero values and corresponding indices of first-derivative (wrt parameters) of second-derivative (wrt all endogenous, auxiliary and exogenous variables) of dynamic model equations, i.e. d(d(df/dvars)/dvars)/dparam, in DR order
d2Yss = d2ys(oo.dr.order_var,indpmodel,indpmodel);
%[endo_nbr by mod_param_nbr by mod_param_nbr], i.e. put into DR order and focus only on model parameters
else
[~, gp] = feval([fname,'.dynamic_params_derivs'], yy0, repmat(ex0, [maximum_exo_lag+maximum_exo_lead+1,1]), params0, ys, 1, dys, d2ys);
%gp is [endo_nbr by (dynamicvar_nbr + exo_nbr) by param_nbr] first-derivative (wrt parameters) of first-derivative (wrt all endogenous, auxiliary and exogenous variables) of dynamic model equations, i.e. d(df/dvars)/dparam, in DR order
[~, g1, g2 ] = feval([fname,'.dynamic'], yy0, repmat(ex0, [maximum_exo_lag+maximum_exo_lead+1,1]), params0, ys, 1);
%g1 is [endo_nbr by (dynamicvar_nbr+exo_nbr)] first derivative (wrt all endogenous, exogenous and auxiliary variables) of dynamic model equations, i.e. df/d[yy0;ex0], in DR order
%g2 is [endo_nbr by (dynamicvar_nbr+exo_nbr)^2] second derivative (wrt all endogenous, exogenous and auxiliary variables) of dynamic model equations, i.e. d(df/d[yy0;ex0])/d[yy0;ex0], in DR order
end
yy0ex0_nbr = sqrt(size(g2,2)); % number of dynamic variables + exogenous variables (length(yy0)+length(ex0))
dYss = dys(oo.dr.order_var, indpmodel); %focus only on model parameters, note dys is in declaration order, dYss is in DR-order
dyy0 = dys(I,:);
yy0_nbr = max(max(lead_lag_incidence)); % retrieve the number of states excluding columns for shocks
% Computation of dg1, i.e. first derivative (wrt. parameters) of Jacobian (wrt endogenous and auxilary variables) of dynamic model using the implicit function theorem
% Let g1 denote the Jacobian of dynamic model equations, i.e. g1 = df/d[yy0ex0], evaluated at the steady state
% Let dg1 denote the first-derivative (wrt parameters) of g1 evaluated at the steady state
% Note that g1 is a function of both the parameters and of the steady state, which also depends on the parameters.
% Hence, implicitly g1=g1(p,yy0ex0(p)) and dg1 consists of two parts (see Ratto and Iskrev (2011) formula 7):
% (1) direct derivative wrt to parameters given by the preprocessor, i.e. gp
% and
% (2) contribution of derivative of steady state (wrt parameters), i.e. g2*dyy0
% Note that in a stochastic context ex0 is always zero and hence can be skipped in the computations
dg1_part2 = gp*0; %initialize part 2, it has dimension [endo_nbr by (dynamicvar_nbr+exo_nbr) by param_nbr]
for j = 1:endo_nbr
[II, JJ] = ind2sub([yy0ex0_nbr yy0ex0_nbr], find(g2(j,:)));
%g2 is [endo_nbr by (dynamicvar_nbr+exo_nbr)^2]
for i = 1:yy0ex0_nbr
is = find(II==i);
is = is(find(JJ(is)<=yy0_nbr)); %focus only on yy0 derivatives as ex0 variables are 0 in a stochastic context
if ~isempty(is)
tmp_g2 = full(g2(j,find(g2(j,:))));
dg1_part2(j,i,:) = tmp_g2(is)*dyy0(JJ(is),:); %put into tensor notation
end
end
end
dg1 = gp + dg1_part2; %dg is sum of two parts due to implicit function theorem
dg1 = dg1(:,:,indpmodel); %focus only on model parameters
if nargout > 6
% Computation of d2g1, i.e. second derivative (wrt. parameters) of Jacobian (wrt endogenous and auxilary variables) of dynamic model using the implicit function theorem
% Let g1 denote the Jacobian of dynamic model equations, i.e. g1 = df/d[yy0ex0], evaluated at the steady state
% Let d2g1 denote the second-derivative (wrt parameters) of g1
% Note that g1 is a function of both the parameters and of the steady state, which also depends on the parameters.
% Hence, implicitly g1=g1(p,yy0ex0(p)) and the first derivative is given by dg1 = gp + g2*dyy0ex0 (see above)
% Accordingly, d2g1, the second-derivative (wrt parameters), consists of five parts (ignoring transposes, see Ratto and Iskrev (2011) formula 16)
% (1) d(gp)/dp = gpp
% (2) d(gp)/dyy0ex0*d(yy0ex0)/dp = hp * dyy0ex0
% (3) d(g2)/dp * dyy0ex0 = hp * dyy0ex0
% (4) d(g2)/dyy0ex0*d(dyy0ex0)/dp * dyy0ex0 = g3 * dyy0ex0 * dyy0ex0
% (5) g2 * d(dyy0ex0)/dp = g2 * d2yy0ex0
% Note that part 2 and 3 are equivalent besides the use of transpose (see Ratto and Iskrev (2011) formula 16)
d2g1_full = sparse(endo_nbr*yy0ex0_nbr, param_nbr*param_nbr); %initialize
dyy0ex0 = sparse([dyy0; zeros(yy0ex0_nbr-yy0_nbr,param_nbr)]); %Jacobian (wrt model parameters) of steady state of dynamic (endogenous and auxiliary) and exogenous variables
g3 = unfold_g3(g3, yy0ex0_nbr);
g3_tmp = reshape(g3,[endo_nbr*yy0ex0_nbr*yy0ex0_nbr yy0ex0_nbr]);
d2g1_part4_left = sparse(endo_nbr*yy0ex0_nbr*yy0ex0_nbr,param_nbr);
for j = 1:param_nbr
%compute first two terms of part 4
d2g1_part4_left(:,j) = g3_tmp*dyy0ex0(:,j);
end
for j=1:endo_nbr
%Note that in the following we focus only on dynamic variables as exogenous variables are 0 by construction in a stochastic setting
d2g1_part5 = reshape(g2(j,:), [yy0ex0_nbr yy0ex0_nbr]);
d2g1_part5 = d2g1_part5(:,1:yy0_nbr)*reshape(d2ys(I,:,:),[yy0_nbr,param_nbr*param_nbr]);
for i=1:yy0ex0_nbr
ind_part4 = sub2ind([endo_nbr yy0ex0_nbr yy0ex0_nbr], ones(yy0ex0_nbr,1)*j ,ones(yy0ex0_nbr,1)*i, (1:yy0ex0_nbr)');
d2g1_part4 = (d2g1_part4_left(ind_part4,:))'*dyy0ex0;
d2g1_part2_and_part3 = (get_hess_deriv(hp,j,i,yy0ex0_nbr,param_nbr))'*dyy0ex0;
d2g1_part1 = get_2nd_deriv_mat(gpp,j,i,param_nbr);
d2g1_tmp = d2g1_part1 + d2g1_part2_and_part3 + d2g1_part2_and_part3' + d2g1_part4 + reshape(d2g1_part5(i,:,:),[param_nbr param_nbr]);
d2g1_tmp = d2g1_tmp(indpmodel,indpmodel); %focus only on model parameters
if any(any(d2g1_tmp))
ind_d2g1_tmp = find(triu(d2g1_tmp));
d2g1_full(sub2ind([endo_nbr yy0ex0_nbr],j,i), ind_d2g1_tmp) = transpose(d2g1_tmp(ind_d2g1_tmp));
end
end
end
clear d2g1_tmp d2g1_part1 d2g1_part2_and_part3 d2g1_part4 d2g1_part4_left d2g1_part5
%store only nonzero entries and the corresponding indices of d2g1:
% rows: respective derivative term
% 1st column: equation number of the term appearing
% 2nd column: column number of variable in Jacobian of the dynamic model
% 3rd column: number of the first parameter in derivative
% 4th column: number of the second parameter in derivative
% 5th column: value of the Hessian term
ind_d2g1 = find(d2g1_full);
d2g1 = zeros(length(ind_d2g1),5);
for j=1:length(ind_d2g1)
[i1, i2] = ind2sub(size(d2g1_full),ind_d2g1(j));
[ig1, ig2] = ind2sub(size(g1),i1);
[ip1, ip2] = ind2sub([modparam_nbr modparam_nbr],i2);
d2g1(j,:) = [ig1 ig2 ip1 ip2 d2g1_full(ind_d2g1(j))];
end
clear d2g1_full;
end
end
% clear variables that are not used any more
clear rp_static g1_static
clear ys dys dyy0 dyy0ex0
clear dg1_part2 tmp_g2
clear g2 gp rpp_static g2_static gp_static d2ys
clear hp g3 g3_tmp gpp
clear ind_d2g1 ind_d2g1_tmp ind_part4 i j i1 i2 ig1 ig2 I II JJ ip1 ip2 is
% Construct nonzero derivatives wrt to t+1, t, and t-1 variables using kstate
klen = maximum_endo_lag + maximum_endo_lead + 1; %total length
k11 = lead_lag_incidence(find([1:klen] ~= maximum_endo_lag+1),:);
g1nonzero = g1(:,nonzeros(k11'));
dg1nonzero = dg1(:,nonzeros(k11'),:);
if nargout > 6
indind = ismember(d2g1(:,2),nonzeros(k11'));
tmp = d2g1(indind,:);
d2g1nonzero = tmp;
for j = 1:size(tmp,1)
inxinx = find(nonzeros(k11')==tmp(j,2));
d2g1nonzero(j,2) = inxinx;
end
end
kstate = oo.dr.kstate;
% Construct nonzero derivatives wrt to t+1, i.e. GAM1=-f_{y^+} in Villemot (2011)
GAM1 = zeros(endo_nbr,endo_nbr);
dGAM1 = zeros(endo_nbr,endo_nbr,modparam_nbr);
k1 = find(kstate(:,2) == maximum_endo_lag+2 & kstate(:,3));
GAM1(:, kstate(k1,1)) = -g1nonzero(:,kstate(k1,3));
dGAM1(:, kstate(k1,1), :) = -dg1nonzero(:,kstate(k1,3),:);
if nargout > 6
indind = ismember(d2g1nonzero(:,2),kstate(k1,3));
tmp = d2g1nonzero(indind,:);
tmp(:,end)=-tmp(:,end);
d2GAM1 = tmp;
for j = 1:size(tmp,1)
inxinx = (kstate(k1,3)==tmp(j,2));
d2GAM1(j,2) = kstate(k1(inxinx),1);
end
end
% Construct nonzero derivatives wrt to t, i.e. GAM0=f_{y^0} in Villemot (2011)
[~,cols_b,cols_j] = find(lead_lag_incidence(maximum_endo_lag+1, oo.dr.order_var));
GAM0 = zeros(endo_nbr,endo_nbr);
dGAM0 = zeros(endo_nbr,endo_nbr,modparam_nbr);
GAM0(:,cols_b) = g1(:,cols_j);
dGAM0(:,cols_b,:) = dg1(:,cols_j,:);
if nargout > 6
indind = ismember(d2g1(:,2),cols_j);
tmp = d2g1(indind,:);
d2GAM0 = tmp;
for j = 1:size(tmp,1)
inxinx = (cols_j==tmp(j,2));
d2GAM0(j,2) = cols_b(inxinx);
end
end
% Construct nonzero derivatives wrt to t-1, i.e. GAM2=-f_{y^-} in Villemot (2011)
k2 = find(kstate(:,2) == maximum_endo_lag+1 & kstate(:,4));
GAM2 = zeros(endo_nbr,endo_nbr);
dGAM2 = zeros(endo_nbr,endo_nbr,modparam_nbr);
GAM2(:, kstate(k2,1)) = -g1nonzero(:,kstate(k2,4));
dGAM2(:, kstate(k2,1), :) = -dg1nonzero(:,kstate(k2,4),:);
if nargout > 6
indind = ismember(d2g1nonzero(:,2),kstate(k2,4));
tmp = d2g1nonzero(indind,:);
tmp(:,end) = -tmp(:,end);
d2GAM2 = tmp;
for j = 1:size(tmp,1)
inxinx = (kstate(k2,4)==tmp(j,2));
d2GAM2(j,2) = kstate(k2(inxinx),1);
end
end
% Construct nonzero derivatives wrt to u_t, i.e. GAM3=-f_{u} in Villemot (2011)
GAM3 = -g1(:,length(yy0)+1:end);
dGAM3 = -dg1(:,length(yy0)+1:end,:);
if nargout > 6
cols_ex = [length(yy0)+1:size(g1,2)];
indind = ismember(d2g1(:,2),cols_ex);
tmp = d2g1(indind,:);
tmp(:,end) = -tmp(:,end);
d2GAM3 = tmp;
for j = 1:size(tmp,1)
inxinx = find(cols_ex==tmp(j,2));
d2GAM3(j,2) = inxinx;
end
clear d2g1 d2g1nonzero tmp
end
clear cols_b cols_ex cols_j k1 k11 k2 klen kstate
clear g1nonzero dg1nonzero g1 yy0
%% Construct first derivative of Sigma_e
dSigma_e = zeros(exo_nbr,exo_nbr,totparam_nbr); %initialize
% note that derivatives wrt model parameters are zero by construction
% Compute first derivative of Sigma_e wrt stderr parameters (these come first)
if ~isempty(indpstderr)
for jp = 1:stderrparam_nbr
dSigma_e(indpstderr(jp),indpstderr(jp),jp) = 2*stderr_e0(indpstderr(jp));
if isdiag(Sigma_e0) == 0 % if there are correlated errors add cross derivatives
indotherex0 = 1:exo_nbr;
indotherex0(indpstderr(jp)) = [];
for kk = indotherex0
dSigma_e(indpstderr(jp), kk, jp) = Corr_e0(indpstderr(jp),kk)*stderr_e0(kk);
dSigma_e(kk, indpstderr(jp), jp) = dSigma_e(indpstderr(jp), kk, jp); %symmetry
end
end
end
end
% Compute first derivative of Sigma_e wrt corr parameters (these come second)
if ~isempty(indpcorr)
for jp = 1:corrparam_nbr
dSigma_e(indpcorr(jp,1),indpcorr(jp,2),stderrparam_nbr+jp) = stderr_e0(indpcorr(jp,1))*stderr_e0(indpcorr(jp,2));
dSigma_e(indpcorr(jp,2),indpcorr(jp,1),stderrparam_nbr+jp) = dSigma_e(indpcorr(jp,1),indpcorr(jp,2),stderrparam_nbr+jp); %symmetry
end
end
%% Construct second derivative of Sigma_e
if nargout > 6
% note that derivatives wrt (mod x mod) and (corr x corr) parameters
% are zero by construction; hence we only need to focus on (stderr x stderr), and (stderr x corr)
d2Sigma_e = zeros(exo_nbr,exo_nbr,totparam_nbr^2); %initialize full matrix, even though we'll reduce it later on to unique upper triangular values
% Compute upper triangular values of Hessian of Sigma_e wrt (stderr x stderr) parameters
if ~isempty(indp2stderrstderr)
for jp = 1:stderrparam_nbr
for ip = 1:jp
if jp == ip %same stderr parameters
d2Sigma_e(indpstderr(jp),indpstderr(jp),indp2stderrstderr(ip,jp)) = 2;
else %different stderr parameters
if isdiag(Sigma_e0) == 0 % if there are correlated errors
d2Sigma_e(indpstderr(jp),indpstderr(ip),indp2stderrstderr(ip,jp)) = Corr_e0(indpstderr(jp),indpstderr(ip));
d2Sigma_e(indpstderr(ip),indpstderr(jp),indp2stderrstderr(ip,jp)) = Corr_e0(indpstderr(jp),indpstderr(ip)); %symmetry
end
end
end
end
end
% Compute upper triangular values of Hessian of Sigma_e wrt (stderr x corr) parameters
if ~isempty(indp2stderrcorr)
for jp = 1:stderrparam_nbr
for ip = 1:corrparam_nbr
if indpstderr(jp) == indpcorr(ip,1) %if stderr equal to first index of corr parameter, derivative is equal to stderr corresponding to second index
d2Sigma_e(indpstderr(jp),indpcorr(ip,2),indp2stderrcorr(jp,ip)) = stderr_e0(indpcorr(ip,2));
d2Sigma_e(indpcorr(ip,2),indpstderr(jp),indp2stderrcorr(jp,ip)) = stderr_e0(indpcorr(ip,2)); % symmetry
end
if indpstderr(jp) == indpcorr(ip,2) %if stderr equal to second index of corr parameter, derivative is equal to stderr corresponding to first index
d2Sigma_e(indpstderr(jp),indpcorr(ip,1),indp2stderrcorr(jp,ip)) = stderr_e0(indpcorr(ip,1));
d2Sigma_e(indpcorr(ip,1),indpstderr(jp),indp2stderrcorr(jp,ip)) = stderr_e0(indpcorr(ip,1)); % symmetry
end
end
end
end
d2Sigma_e = d2Sigma_e(:,:,indp2tottot2); %focus on upper triangular hessian values
end
if kronflag == 1
% The following derivations are based on Iskrev (2010) and its online appendix A.
% Basic idea is to make use of the implicit function theorem.
% Let F = GAM0*A - GAM1*A*A - GAM2 = 0
% Note that F is a function of parameters p and A, which is also a
% function of p,therefore, F = F(p,A(p)), and hence,
% dF = Fp + dF_dA*dA or dA = - Fp/dF_dA
% Some auxiliary matrices
I_endo = speye(endo_nbr);
I_exo = speye(exo_nbr);
% Reshape to write derivatives in the Magnus and Neudecker style, i.e. dvec(X)/dp
dGAM0 = reshape(dGAM0, endo_nbr^2, modparam_nbr);
dGAM1 = reshape(dGAM1, endo_nbr^2, modparam_nbr);
dGAM2 = reshape(dGAM2, endo_nbr^2, modparam_nbr);
dGAM3 = reshape(dGAM3, endo_nbr*exo_nbr, modparam_nbr);
dSigma_e = reshape(dSigma_e, exo_nbr^2, totparam_nbr);
% Compute dA via implicit function
dF_dA = kron(I_endo,GAM0) - kron(A',GAM1) - kron(I_endo,GAM1*A); %equation 31 in Appendix A of Iskrev (2010)
Fp = kron(A',I_endo)*dGAM0 - kron( (A')^2,I_endo)*dGAM1 - dGAM2; %equation 32 in Appendix A of Iskrev (2010)
dA = -dF_dA\Fp;
% Compute dB from expressions 33 in Iskrev (2010) Appendix A
MM = GAM0-GAM1*A; %this corresponds to matrix M in Ratto and Iskrev (2011, page 6) and will be used if nargout > 6 below
invMM = MM\eye(endo_nbr);
dB = - kron( (invMM*GAM3)' , invMM ) * ( dGAM0 - kron( A' , I_endo ) * dGAM1 - kron( I_endo , GAM1 ) * dA ) + kron( I_exo, invMM ) * dGAM3 ;
dBt = commutation(endo_nbr, exo_nbr)*dB; %transose of derivative using the commutation matrix
% Add derivatives for stderr and corr parameters, which are zero by construction
dA = [zeros(endo_nbr^2, stderrparam_nbr+corrparam_nbr) dA];
dB = [zeros(endo_nbr*exo_nbr, stderrparam_nbr+corrparam_nbr) dB];
dBt = [zeros(endo_nbr*exo_nbr, stderrparam_nbr+corrparam_nbr) dBt];
% Compute dOm = dvec(B*Sig*B') from expressions 34 in Iskrev (2010) Appendix A
dOm = kron(I_endo,B*Sigma_e0)*dBt + kron(B,B)*dSigma_e + kron(B*Sigma_e0,I_endo)*dB;
% Put into tensor notation
dA = reshape(dA, endo_nbr, endo_nbr, totparam_nbr);
dB = reshape(dB, endo_nbr, exo_nbr, totparam_nbr);
dOm = reshape(dOm, endo_nbr, endo_nbr, totparam_nbr);
dSigma_e = reshape(dSigma_e, exo_nbr, exo_nbr, totparam_nbr);
if nargout > 6
% Put back into tensor notation as these will be reused later
dGAM0 = reshape(dGAM0, endo_nbr, endo_nbr, modparam_nbr);
dGAM1 = reshape(dGAM1, endo_nbr, endo_nbr, modparam_nbr);
dGAM2 = reshape(dGAM2, endo_nbr, endo_nbr, modparam_nbr);
dGAM3 = reshape(dGAM3, endo_nbr, exo_nbr, modparam_nbr);
dAA = dA(:, :, stderrparam_nbr+corrparam_nbr+1:end); %this corresponds to matrix dA in Ratto and Iskrev (2011, page 6), i.e. derivative of A with respect to model parameters only in tensor notation
dBB = dB(:, :, stderrparam_nbr+corrparam_nbr+1:end); %dBB is for all endogenous variables, whereas dB is only for selected variables
N = -GAM1; %this corresponds to matrix N in Ratto and Iskrev (2011, page 6)
P = A; %this corresponds to matrix P in Ratto and Iskrev (2011, page 6)
end
% Focus only on selected variables
dYss = dYss(indvar,:);
dA = dA(indvar,indvar,:);
dB = dB(indvar,:,:);
dOm = dOm(indvar,indvar,:);
elseif (kronflag == 0 || kronflag == -2)
% generalized sylvester equation solves MM*dAA+N*dAA*P=Q from Ratto and Iskrev (2011) equation 11 where
% dAA is derivative of A with respect to model parameters only in tensor notation
MM = (GAM0-GAM1*A);
N = -GAM1;
P = A;
Q_rightpart = zeros(endo_nbr,endo_nbr,modparam_nbr); %initialize
Q = Q_rightpart; %initialize and compute matrix Q in Ratto and Iskrev (2011, page 6)
for j = 1:modparam_nbr
Q_rightpart(:,:,j) = (dGAM0(:,:,j)-dGAM1(:,:,j)*A);
Q(:,:,j) = dGAM2(:,:,j)-Q_rightpart(:,:,j)*A;
end
%use iterated generalized sylvester equation to compute dAA
dAA = sylvester3(MM,N,P,Q);
flag = 1; icount = 0;
while flag && icount < 4
[dAA, flag] = sylvester3a(dAA,MM,N,P,Q);
icount = icount+1;
end
%stderr parameters come first, then corr parameters, model parameters come last
%note that stderr and corr derivatives are:
% - zero by construction for A and B
% - depend only on dSig for Om
dOm = zeros(var_nbr, var_nbr, totparam_nbr);
dA = zeros(var_nbr, var_nbr, totparam_nbr);
dB = zeros(var_nbr, exo_nbr, totparam_nbr);
if nargout > 6
dBB = zeros(endo_nbr, exo_nbr, modparam_nbr); %dBB is always for all endogenous variables, whereas dB is only for selected variables
end
%compute derivative of Om=B*Sig*B' that depends on Sig (other part is added later)
if ~isempty(indpstderr)
for j = 1:stderrparam_nbr
BSigjBt = B*dSigma_e(:,:,j)*B';
dOm(:,:,j) = BSigjBt(indvar,indvar);
end
end
if ~isempty(indpcorr)
for j = 1:corrparam_nbr
BSigjBt = B*dSigma_e(:,:,stderrparam_nbr+j)*B';
dOm(:,:,stderrparam_nbr+j) = BSigjBt(indvar,indvar);
end
end
%compute derivative of B and the part of Om=B*Sig*B' that depends on B (other part is computed above)
invMM = inv(MM);
for j = 1:modparam_nbr
dAAj = dAA(:,:,j);
dBj = invMM * ( dGAM3(:,:,j) - (Q_rightpart(:,:,j) -GAM1*dAAj ) * B ); %equation 14 in Ratto and Iskrev (2011), except in the paper there is a typo as the last B is missing
dOmj = dBj*Sigma_e0*B'+B*Sigma_e0*dBj';
%store derivatives in tensor notation
dA(:, :, stderrparam_nbr+corrparam_nbr+j) = dAAj(indvar,indvar);
dB(:, :, stderrparam_nbr+corrparam_nbr+j) = dBj(indvar,:);
dOm(:, :, stderrparam_nbr+corrparam_nbr+j) = dOmj(indvar,indvar);
if nargout > 6
dBB(:, :, j) = dBj;
end
end
dYss = dYss(indvar,:); % Focus only on relevant variables
end
%% Compute second-order derivatives (wrt params) of solution matrices using generalized sylvester equations, see equations 17 and 18 in Ratto and Iskrev (2011)
if nargout > 6
% solves MM*d2AA+N*d2AA*P = QQ where d2AA are second order derivatives (wrt model parameters) of A
d2Yss = d2Yss(indvar,:,:);
QQ = zeros(endo_nbr,endo_nbr,floor(sqrt(modparam_nbr2)));
jcount=0;
cumjcount=0;
jinx = [];
x2x=sparse(endo_nbr*endo_nbr,modparam_nbr2);
for i=1:modparam_nbr
for j=1:i
elem1 = (get_2nd_deriv(d2GAM0,endo_nbr,endo_nbr,j,i)-get_2nd_deriv(d2GAM1,endo_nbr,endo_nbr,j,i)*A);
elem1 = get_2nd_deriv(d2GAM2,endo_nbr,endo_nbr,j,i)-elem1*A;
elemj0 = dGAM0(:,:,j)-dGAM1(:,:,j)*A;
elemi0 = dGAM0(:,:,i)-dGAM1(:,:,i)*A;
elem2 = -elemj0*dAA(:,:,i)-elemi0*dAA(:,:,j);
elem2 = elem2 + ( dGAM1(:,:,j)*dAA(:,:,i) + dGAM1(:,:,i)*dAA(:,:,j) )*A;
elem2 = elem2 + GAM1*( dAA(:,:,i)*dAA(:,:,j) + dAA(:,:,j)*dAA(:,:,i));
jcount=jcount+1;
jinx = [jinx; [j i]];
QQ(:,:,jcount) = elem1+elem2;
if jcount==floor(sqrt(modparam_nbr2)) || (j*i)==modparam_nbr^2
if (j*i)==modparam_nbr^2
QQ = QQ(:,:,1:jcount);
end
xx2=sylvester3(MM,N,P,QQ);
flag=1;
icount=0;
while flag && icount<4
[xx2, flag]=sylvester3a(xx2,MM,N,P,QQ);
icount = icount + 1;
end
x2x(:,cumjcount+1:cumjcount+jcount)=reshape(xx2,[endo_nbr*endo_nbr jcount]);
cumjcount=cumjcount+jcount;
jcount = 0;
jinx = [];
end
end
end
clear d xx2;
jcount = 0;
icount = 0;
cumjcount = 0;
MAX_DIM_MAT = 100000000;
ncol = max(1,floor(MAX_DIM_MAT/(8*var_nbr*(var_nbr+1)/2)));
ncol = min(ncol, totparam_nbr2);
d2A = sparse(var_nbr*var_nbr,totparam_nbr2);
d2Om = sparse(var_nbr*(var_nbr+1)/2,totparam_nbr2);
d2A_tmp = zeros(var_nbr*var_nbr,ncol);
d2Om_tmp = zeros(var_nbr*(var_nbr+1)/2,ncol);
tmpDir = CheckPath('tmp_derivs',dname);
offset = stderrparam_nbr+corrparam_nbr;
% d2B = zeros(m,n,tot_param_nbr,tot_param_nbr);
for j=1:totparam_nbr
for i=1:j
jcount=jcount+1;
if j<=offset %stderr and corr parameters
y = B*d2Sigma_e(:,:,jcount)*B';
d2Om_tmp(:,jcount) = dyn_vech(y(indvar,indvar));
else %model parameters
jind = j-offset;
iind = i-offset;
if i<=offset
y = dBB(:,:,jind)*dSigma_e(:,:,i)*B'+B*dSigma_e(:,:,i)*dBB(:,:,jind)';
% y(abs(y)<1.e-8)=0;
d2Om_tmp(:,jcount) = dyn_vech(y(indvar,indvar));
else
icount=icount+1;
dAAj = reshape(x2x(:,icount),[endo_nbr endo_nbr]);
% x = get_2nd_deriv(x2x,m,m,iind,jind);%xx2(:,:,jcount);
elem1 = (get_2nd_deriv(d2GAM0,endo_nbr,endo_nbr,iind,jind)-get_2nd_deriv(d2GAM1,endo_nbr,endo_nbr,iind,jind)*A);
elem1 = elem1 -( dGAM1(:,:,jind)*dAA(:,:,iind) + dGAM1(:,:,iind)*dAA(:,:,jind) );
elemj0 = dGAM0(:,:,jind)-dGAM1(:,:,jind)*A-GAM1*dAA(:,:,jind);
elemi0 = dGAM0(:,:,iind)-dGAM1(:,:,iind)*A-GAM1*dAA(:,:,iind);
elem0 = elemj0*dBB(:,:,iind)+elemi0*dBB(:,:,jind);
y = invMM * (get_2nd_deriv(d2GAM3,endo_nbr,exo_nbr,iind,jind)-elem0-(elem1-GAM1*dAAj)*B);
% d2B(:,:,j+length(indexo),i+length(indexo)) = y;
% d2B(:,:,i+length(indexo),j+length(indexo)) = y;
y = y*Sigma_e0*B'+B*Sigma_e0*y'+ ...
dBB(:,:,jind)*Sigma_e0*dBB(:,:,iind)'+dBB(:,:,iind)*Sigma_e0*dBB(:,:,jind)';
% x(abs(x)<1.e-8)=0;
d2A_tmp(:,jcount) = vec(dAAj(indvar,indvar));
% y(abs(y)<1.e-8)=0;
d2Om_tmp(:,jcount) = dyn_vech(y(indvar,indvar));
end
end
if jcount==ncol || i*j==totparam_nbr^2
d2A(:,cumjcount+1:cumjcount+jcount) = d2A_tmp(:,1:jcount);
% d2A(:,:,j+length(indexo),i+length(indexo)) = x;
% d2A(:,:,i+length(indexo),j+length(indexo)) = x;
d2Om(:,cumjcount+1:cumjcount+jcount) = d2Om_tmp(:,1:jcount);
% d2Om(:,:,j+length(indexo),i+length(indexo)) = y;
% d2Om(:,:,i+length(indexo),j+length(indexo)) = y;
save([tmpDir filesep 'd2A_' int2str(cumjcount+1) '_' int2str(cumjcount+jcount) '.mat'],'d2A')
save([tmpDir filesep 'd2Om_' int2str(cumjcount+1) '_' int2str(cumjcount+jcount) '.mat'],'d2Om')
cumjcount = cumjcount+jcount;
jcount=0;
% d2A = sparse(m1*m1,tot_param_nbr*(tot_param_nbr+1)/2);
% d2Om = sparse(m1*(m1+1)/2,tot_param_nbr*(tot_param_nbr+1)/2);
d2A_tmp = zeros(var_nbr*var_nbr,ncol);
d2Om_tmp = zeros(var_nbr*(var_nbr+1)/2,ncol);
end
end
end
end
return
function g22 = get_2nd_deriv(gpp,m,n,i,j)
% inputs:
% - gpp: [#second_order_Jacobian_terms by 5] double Hessian matrix (wrt parameters) of a matrix
% rows: respective derivative term
% 1st column: equation number of the term appearing
% 2nd column: column number of variable in Jacobian
% 3rd column: number of the first parameter in derivative
% 4th column: number of the second parameter in derivative
% 5th column: value of the Hessian term
% - m: scalar number of equations
% - n: scalar number of variables
% - i: scalar number for which first parameter
% - j: scalar number for which second parameter
g22=zeros(m,n);
is=find(gpp(:,3)==i);
is=is(find(gpp(is,4)==j));
if ~isempty(is)
g22(sub2ind([m,n],gpp(is,1),gpp(is,2)))=gpp(is,5)';
end
return
function g22 = get_2nd_deriv_mat(gpp,i,j,npar)
% inputs:
% - gpp: [#second_order_Jacobian_terms by 5] double Hessian matrix of (wrt parameters) of dynamic Jacobian
% rows: respective derivative term
% 1st column: equation number of the term appearing
% 2nd column: column number of variable in Jacobian of the dynamic model
% 3rd column: number of the first parameter in derivative
% 4th column: number of the second parameter in derivative
% 5th column: value of the Hessian term
% - i: scalar number for which model equation
% - j: scalar number for which variable in Jacobian of dynamic model
% - npar: scalar Number of model parameters, i.e. equals M_.param_nbr
%
% output:
% g22: [npar by npar] Hessian matrix (wrt parameters) of Jacobian of dynamic model for equation i
% rows: first parameter in Hessian
% columns: second paramater in Hessian
g22=zeros(npar,npar);
is=find(gpp(:,1)==i);
is=is(find(gpp(is,2)==j));
if ~isempty(is)
g22(sub2ind([npar,npar],gpp(is,3),gpp(is,4)))=gpp(is,5)';
end
return
function g22 = get_all_2nd_derivs(gpp,m,n,npar,fsparse)
if nargin==4 || isempty(fsparse)
fsparse=0;
end
if fsparse
g22=sparse(m*n,npar*npar);
else
g22=zeros(m,n,npar,npar);
end
% c=ones(npar,npar);
% c=triu(c);
% ic=find(c);
for is=1:length(gpp)
% d=zeros(npar,npar);
% d(gpp(is,3),gpp(is,4))=1;
% indx = find(ic==find(d));
if fsparse
g22(sub2ind([m,n],gpp(is,1),gpp(is,2)),sub2ind([npar,npar],gpp(is,3),gpp(is,4)))=gpp(is,5);
else
g22(gpp(is,1),gpp(is,2),gpp(is,3),gpp(is,4))=gpp(is,5);
end
end
return
function r22 = get_all_resid_2nd_derivs(rpp,m,npar)
% inputs:
% - rpp: [#second_order_residual_terms by 4] double Hessian matrix (wrt paramters) of model equations
% rows: respective derivative term
% 1st column: equation number of the term appearing
% 2nd column: number of the first parameter in derivative
% 3rd column: number of the second parameter in derivative
% 4th column: value of the Hessian term
% - m: scalar Number of residuals (or model equations), i.e. equals endo_nbr
% - npar: scalar Number of model parameters, i.e. equals param_nbr
%
% output:
% r22: [endo_nbr by param_nbr by param_nbr] Hessian matrix of model equations with respect to parameters
% rows: equations in order of declaration
% 1st columns: first parameter number in derivative
% 2nd columns: second parameter in derivative
r22=zeros(m,npar,npar);
for is=1:length(rpp)
% Keep symmetry in hessian, hence 2 and 3 as well as 3 and 2, i.e. d2f/(dp1 dp2) = d2f/(dp2 dp1)
r22(rpp(is,1),rpp(is,2),rpp(is,3))=rpp(is,4);
r22(rpp(is,1),rpp(is,3),rpp(is,2))=rpp(is,4);
end
return
function h2 = get_all_hess_derivs(hp,r,m,npar)
h2=zeros(r,m,m,npar);
for is=1:length(hp)
h2(hp(is,1),hp(is,2),hp(is,3),hp(is,4))=hp(is,5);
end
return
function h2 = get_hess_deriv(hp,i,j,m,npar)
% inputs:
% - hp: [#first_order_Hessian_terms by 5] double Jacobian matrix (wrt paramters) of dynamic Hessian
% rows: respective derivative term
% 1st column: equation number of the term appearing
% 2nd column: column number of first variable in Hessian of the dynamic model
% 3rd column: column number of second variable in Hessian of the dynamic model
% 4th column: number of the parameter in derivative
% 5th column: value of the Hessian term
% - i: scalar number for which model equation
% - j: scalar number for which first variable in Hessian of dynamic model variable
% - m: scalar Number of dynamic model variables + exogenous vars, i.e. dynamicvar_nbr + exo_nbr
% - npar: scalar Number of model parameters, i.e. equals M_.param_nbr
%
% output:
% h2: [(dynamicvar_nbr + exo_nbr) by M_.param_nbr] Jacobian matrix (wrt parameters) of dynamic Hessian
% rows: second dynamic or exogenous variables in Hessian of specific model equation of the dynamic model
% columns: parameters
h2=zeros(m,npar);
is1=find(hp(:,1)==i);
is=is1(find(hp(is1,2)==j));
if ~isempty(is)
h2(sub2ind([m,npar],hp(is,3),hp(is,4)))=hp(is,5)';
end
return

16
matlab/get_identification_jacobians.m
View File

@ -71,7 +71,7 @@ function [J, G, D, dTAU, dLRE, dA, dOm, dYss, MOMENTS] = get_identification_jaco
% * get_minimal_state_representation
% * duplication
% * dyn_vech
% * get_first_order_solution_params_deriv (previously getH)
% * get_perturbation_params_deriv (previously getH)
% * get_all_parameters
% * fjaco
% * lyapunov_symm
@ -124,7 +124,13 @@ exo_nbr = M.exo_nbr;
endo_nbr = M.endo_nbr;
%% Construct dTAU, dLRE, dA, dOm, dYss
[dA, dB, dSig, dOm, dYss, dg1] = get_first_order_solution_params_deriv(A, B, estim_params, M, oo, options, kronflag, indpmodel, indpstderr, indpcorr, (1:endo_nbr)');
DERIVS = get_perturbation_params_derivs(M, options, estim_params, oo, indpmodel, indpstderr, indpcorr, 0);%d2flag=0
dA = zeros(M.endo_nbr, M.endo_nbr, size(DERIVS.dghx,3));
dA(:,M.nstatic+(1:M.nspred),:) = DERIVS.dghx;
dB = DERIVS.dghu;
dSig = DERIVS.dSigma_e;
dOm = DERIVS.dOm;
dYss = DERIVS.dYss;
% Collect terms for derivative of tau=[Yss; vec(A); vech(Om)]
dTAU = zeros(endo_nbr*endo_nbr+endo_nbr*(endo_nbr+1)/2, totparam_nbr);
@ -132,7 +138,7 @@ for j=1:totparam_nbr
dTAU(:,j) = [vec(dA(:,:,j)); dyn_vech(dOm(:,:,j))];
end
dTAU = [ [zeros(endo_nbr, stderrparam_nbr+corrparam_nbr) dYss]; dTAU ]; % add steady state
dLRE = [dYss; reshape(dg1,size(dg1,1)*size(dg1,2),size(dg1,3)) ]; %reshape dg1 and add steady state
dLRE = [dYss; reshape(DERIVS.dg1,size(DERIVS.dg1,1)*size(DERIVS.dg1,2),size(DERIVS.dg1,3)) ]; %reshape dg1 and add steady state
%State Space Matrices for VAROBS variables
C = A(indvobs,:);
@ -144,7 +150,7 @@ dD = dB(indvobs,:,:);
if ~no_identification_moments
if kronflag == -1
%numerical derivative of autocovariogram [outputflag=1]
J = fjaco('identification_numerical_objective', para0, 1, estim_params, M, oo, options, indpmodel, indpstderr, indpcorr, indvobs, useautocorr, nlags, grid_nbr);
J = fjaco(str2func('identification_numerical_objective'), para0, 1, estim_params, M, oo, options, indpmodel, indpstderr, indpcorr, indvobs, useautocorr, nlags, grid_nbr);
M.params = params0; %make sure values are set back
M.Sigma_e = Sigma_e0; %make sure values are set back
M.Correlation_matrix = Corr_e0 ; %make sure values are set back
@ -254,7 +260,7 @@ if ~no_identification_spectrum
IA = eye(size(A,1));
if kronflag == -1
%numerical derivative of spectral density [outputflag=2]
dOmega_tmp = fjaco('identification_numerical_objective', para0, 2, estim_params, M, oo, options, indpmodel, indpstderr, indpcorr, indvobs, useautocorr, nlags, grid_nbr);
dOmega_tmp = fjaco(str2func('identification_numerical_objective'), para0, 2, estim_params, M, oo, options, indpmodel, indpstderr, indpcorr, indvobs, useautocorr, nlags, grid_nbr);
M.params = params0; %make sure values are set back
M.Sigma_e = Sigma_e0; %make sure values are set back
M.Correlation_matrix = Corr_e0 ; %make sure values are set back

1533
matlab/get_perturbation_params_derivs.m Normal file
View File

File diff suppressed because it is too large Load Diff

109
matlab/get_perturbation_params_derivs_numerical_objective.m Normal file
View File

@ -0,0 +1,109 @@
function [out,info] = get_perturbation_params_derivs_numerical_objective(params, outputflag, estim_params, M, oo, options)
%function [out,info] = get_perturbation_params_derivs_numerical_objective(params, outputflag, estim_params, M, oo, options)
% -------------------------------------------------------------------------
% Objective function to compute numerically the Jacobians used for get_perturbation_params_derivs
% =========================================================================
% INPUTS
% params: [vector] parameter values at which to evaluate objective function
% stderr parameters come first, corr parameters second, model parameters third
% outputflag: [string] flag which objective to compute (see below)
% estim_params: [structure] storing the estimation information
% M: [structure] storing the model information
% oo: [structure] storing the solution results
% options: [structure] storing the options
% -------------------------------------------------------------------------
%
% OUTPUT (dependent on outputflag and order of approximation):
% - 'perturbation_solution': out = [vec(Sigma_e);vec(ghx);vec(ghu);vec(ghxx);vec(ghxu);vec(ghuu);vec(ghs2);vec(ghxxx);vec(ghxxu);vec(ghxuu);vec(ghuuu);vec(ghxss);vec(ghuss)]
% - 'dynamic_model': out = [Yss; vec(g1); vec(g2); vec(g3)]
% - 'Kalman_Transition': out = [Yss; vec(KalmanA); dyn_vech(KalmanB*Sigma_e*KalmanB')];
% all in DR-order
% -------------------------------------------------------------------------
% This function is called by
% * get_solution_params_deriv.m (previously getH.m)
% * get_identification_jacobians.m (previously getJJ.m)
% -------------------------------------------------------------------------
% This function calls
% * [M.fname,'.dynamic']
% * resol
% * dyn_vech
% =========================================================================
% 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/>.
% =========================================================================
%% Update stderr, corr and model parameters and compute perturbation approximation and steady state with updated parameters
M = set_all_parameters(params,estim_params,M);
Sigma_e = M.Sigma_e;
[~,info,M,options,oo] = resol(0,M,options,oo);
if info(1) > 0
% there are errors in the solution algorithm
out = [];
return
else
ys = oo.dr.ys; %steady state of model variables in declaration order
ghx = oo.dr.ghx; ghu = oo.dr.ghu;
if options.order > 1
ghxx = oo.dr.ghxx; ghxu = oo.dr.ghxu; ghuu = oo.dr.ghuu; ghs2 = oo.dr.ghs2;
%ghxs = oo.dr.ghxs; ghus = oo.dr.ghus; %these are zero due to certainty equivalence and Gaussian shocks
end
if options.order > 2
ghxxx = oo.dr.ghxxx; ghxxu = oo.dr.ghxxu; ghxuu = oo.dr.ghxuu; ghxss = oo.dr.ghxss; ghuuu = oo.dr.ghuuu; ghuss = oo.dr.ghuss;
%ghxxs = oo.dr.ghxxs; ghxus = oo.dr.ghxus; ghuus = oo.dr.ghuus; ghsss = oo.dr.ghsss; %these are zero due to certainty equivalence and Gaussian shocks
end
end
Yss = ys(oo.dr.order_var); %steady state of model variables in DR order
%% out = [vec(Sigma_e);vec(ghx);vec(ghu);vec(ghxx);vec(ghxu);vec(ghuu);vec(ghs2);vec(ghxxx);vec(ghxxu);vec(ghxuu);vec(ghuuu);vec(ghxss);vec(ghuss)]
if strcmp(outputflag,'perturbation_solution')
out = [Sigma_e(:); ghx(:); ghu(:)];
if options.order > 1
out = [out; ghxx(:); ghxu(:); ghuu(:); ghs2(:);];
end
if options.order > 2
out = [out; ghxxx(:); ghxxu(:); ghxuu(:); ghuuu(:); ghxss(:); ghuss(:)];
end
end
%% out = [Yss; vec(g1); vec(g2); vec(g3)]; of all endogenous variables, in DR order
if strcmp(outputflag,'dynamic_model')
[I,~] = find(M.lead_lag_incidence'); %I is used to evaluate dynamic model files
if options.order == 1
[~, g1] = feval([M.fname,'.dynamic'], ys(I), oo.exo_steady_state', M.params, ys, 1);
out = [Yss; g1(:)];
elseif options.order == 2
[~, g1, g2] = feval([M.fname,'.dynamic'], ys(I), oo.exo_steady_state', M.params, ys, 1);
out = [Yss; g1(:); g2(:)];
elseif options.order == 3
[~, g1, g2, g3] = feval([M.fname,'.dynamic'], ys(I), oo.exo_steady_state', M.params, ys, 1);
g3 = unfold_g3(g3, length(ys(I))+M.exo_nbr);
out = [Yss; g1(:); g2(:); g3(:)];
end
end
%% out = [Yss; vec(KalmanA); dyn_vech(KalmanB*Sigma_e*KalmanB')]; in DR order, where A and B are Kalman transition matrices
if strcmp(outputflag,'Kalman_Transition')
if options.order == 1
KalmanA = zeros(M.endo_nbr,M.endo_nbr);
KalmanA(:,M.nstatic+(1:M.nspred)) = ghx;
Om = ghu*Sigma_e*transpose(ghu);
out = [Yss; KalmanA(:); dyn_vech(Om)];
else
error('''get_perturbation_params_derivs_numerical_objective.m'': Kalman_Transition works only at order=1');
end
end

23
matlab/identification_numerical_objective.m
View File

@ -92,17 +92,6 @@ end
ys = oo.dr.ys; %steady state of model variables in declaration order
y0 = ys(oo.dr.order_var); %steady state of model variables in DR order
%% out = [Yss; vec(A); vec(B); dyn_vech(Sig_e)]; of indvar variables only, in DR order
if outputflag == 0
out = [y0(indvar); vec(A(indvar,indvar)); vec(B(indvar,:)); dyn_vech(M.Sigma_e)];
end
%% out = [Yss; vec(A); dyn_vech(Om)]; of indvar variables only, in DR order
if outputflag == -2
Om = B*M.Sigma_e*transpose(B);
out = [y0(indvar); vec(A(indvar,indvar)); dyn_vech(Om(indvar,indvar))];
end
%% out = [vech(cov(Y_t,Y_t)); vec(cov(Y_t,Y_{t-1}); ...; vec(cov(Y_t,Y_{t-nlags})] of indvar variables, in DR order. This is Iskrev (2010)'s J matrix.
if outputflag == 1
% Denote Ezz0 = E_t(z_t * z_t'), then the following Lyapunov equation defines the autocovariagrom: Ezz0 -A*Ezz*A' = B*Sig_e*B'
@ -145,14 +134,4 @@ if outputflag == 2
kk = kk+1;
out(1 + (kk-1)*var_nbr^2 : kk*var_nbr^2) = g_omega(:);
end
end
%% out = [Yss; vec(g1)]; of all endogenous variables, in DR order
if outputflag == -1
[I,~] = find(M.lead_lag_incidence'); %I is used to evaluate dynamic model files
yy0 = oo.dr.ys(I); %steady state of dynamic model variables in DR order
ex0 = oo.exo_steady_state';
[~, g1] = feval([M.fname,'.dynamic'], yy0, ex0, M.params, ys, 1);
out = [y0; g1(:)];
end
end

46
matlab/unfold_g4.m Normal file
View File

@ -0,0 +1,46 @@
function g4_unfolded = unfold_g4(g4, ny)
% Given the 4th order derivatives stored in a sparse matrix and without
% symmetric elements (as returned by the static/dynamic files) and the number
% of (static or dynamic) variables in the jacobian, returns
% an unfolded version of the same matrix (i.e. with symmetric elements).
% 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/>.
[i, j, v] = find(g4);
i_unfolded = [];
j_unfolded = [];
v_unfolded = [];
for k = 1:length(v)
l1 = rem(j(k)-1, ny);
j2 = floor((j(k)-1)/ny);
l2 = rem(j2, ny);
j3 = floor((j(k)-1)/ny^2);
l3 = rem(j3,ny);
l4 = floor(j3/ny);
p = unique(perms([l1 l2 l3 l4]), 'rows');
np = rows(p);
i_unfolded = [i_unfolded; repmat(i(k), np, 1)];
j_unfolded = [j_unfolded; 1 + p(:,1) + ny*(p(:,2) + ny*(p(:,3) + ny*p(:,4)))];
v_unfolded = [v_unfolded; repmat(v(k), np, 1)];
end
g4_unfolded = sparse(i_unfolded, j_unfolded, v_unfolded, size(g4, 1), size(g4, 2));

2
tests/.gitignore vendored
View File

@ -33,6 +33,7 @@ wsOct
!/AIM/fs2000_b1L1L_steadystate.m
!/AIM/fsdat.m
!/analytic_derivatives/fsdat_simul.m
!/analytic_derivatives/nBrockMirmanSYM.mat
!/block_bytecode/run_ls2003.m
!/bvar_a_la_sims/bvar_sample.m
!/conditional_forecasts/2/fsdat_simul.m
@ -67,6 +68,7 @@ wsOct
!/gsa/ls2003scr_results.mat
!/gsa/morris/nk_est_data.m
!/identification/as2007/as2007_steadystate.m
!/identification/as2007/G_QT.mat
!/identification/kim/kim2_steadystate.m
!/initval_file/ramst_initval_file_data.m
!/internals/tests.m

9
tests/Makefile.am
View File

@ -19,6 +19,8 @@ MODFILES = \
optimizers/fs2000_9.mod \
optimizers/fs2000_10.mod \
analytic_derivatives/fs2000_analytic_derivation.mod \
analytic_derivatives/BrockMirman_PertParamsDerivs.mod \
analytic_derivatives/burnside_3_order_PertParamsDerivs.mod \
measurement_errors/fs2000_corr_me_ml_mcmc/fs2000_corr_ME.mod \
TeX/fs2000_corr_ME.mod \
estimation/MH_recover/fs2000_recover_tarb.mod \
@ -383,6 +385,7 @@ XFAIL_MODFILES = ramst_xfail.mod \
estimation/fs2000_mixed_ML_xfail.mod \
estimation/fs2000_stochastic_singularity_xfail.mod \
identification/ident_unit_root/ident_unit_root_xfail.mod \
identification/LindeTrabandt/LindeTrabandt2019_xfail.mod \
steady_state/Linear_steady_state_xfail.mod \
optimal_policy/Ramsey/ramsey_histval_xfail.mod \
example1_extra_exo_xfail.mod
@ -696,6 +699,10 @@ identification: m/identification o/identification
m/identification: $(patsubst %.mod, %.m.trs, $(filter identification/%.mod, $(MODFILES)))
o/identification: $(patsubst %.mod, %.o.trs, $(filter identification/%.mod, $(MODFILES)))
analytic_derivatives: m/analytic_derivatives o/analytic_derivatives
m/analytic_derivatives: $(patsubst %.mod, %.m.trs, $(filter analytic_derivatives/%.mod, $(MODFILES)))
o/analytic_derivatives: $(patsubst %.mod, %.o.trs, $(filter analytic_derivatives/%.mod, $(MODFILES)))
fs2000: m/fs2000 o/fs2000
m/fs2000: $(patsubst %.mod, %.m.trs, $(filter fs2000/%.mod, $(MODFILES)))
o/fs2000: $(patsubst %.mod, %.o.trs, $(filter fs2000/%.mod, $(MODFILES)))
@ -813,6 +820,7 @@ EXTRA_DIST = \
data/test.xlsx \
gsa/morris/nk_est_data.m \
analytic_derivatives/fsdat_simul.m \
analytic_derivatives/nBrockMirmanSYM.mat \
fs2000/fsdat_simul.m \
ls2003/data_ca1.m \
measurement_errors/data_ca1.m \
@ -848,6 +856,7 @@ EXTRA_DIST = \
kalman/likelihood_from_dynare/fs2000_estimation_check.inc \
kalman/likelihood_from_dynare/fs2000ns_model.inc \
kalman/likelihood_from_dynare/fs2000ns_estimation_check.inc \
identification/as2007/G_QT.mat \
estimation/fsdat_simul.m \
ep/mean_preserving_spread.m \
decision_rules/example1_results_dyn_432.mat \

474
tests/analytic_derivatives/BrockMirman_PertParamsDerivs.mod Normal file
View File

@ -0,0 +1,474 @@
% "Augmented" Brock Mirman model
% True policy functions and their exact derivatives (ghx,ghu,ghxx,ghxu,ghuu,ghs2,ghxxx,ghxxu,ghxuu,ghuuu,ghxss,ghuss) are computed using Matlab's symbolic toolbox and saved to a mat file
% Created by @wmutschl (Willi Mutschler, willi@mutschler.eu)
% =========================================================================
% 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/>.
% =========================================================================
@#define CREATE_SYMBOLIC = 0
@#define ORDER = 3
%define parameter values which are used for calibration and estimated_params block
@#define value_SE_A = 0.3
@#define value_SE_X = 0.1
@#define value_SE_Y = 0.9
@#define value_rho_AX = 0.4
@#define value_alph = 0.35
@#define value_betta = 0.99
@#define value_rhoA = 0.9
@#define value_sigA = 0.6
@#define value_sigX = 0.25
@#define value_Xss = 1.2
@#define value_dumA = 2
@#define value_dumK = 1
@#define value_dumepsA = 4
@#define value_dumepsX = 5
%define parameter values which are used for estimated_params block, note that all are different
@#define prior_value_SE_A = 0.8
@#define prior_value_SE_X = 0.6
@#define prior_value_SE_Y = 0.1
@#define prior_value_rho_AX = 0.1
@#define prior_value_alph = 0.15
@#define prior_value_betta = 0.92
@#define prior_value_rhoA = 0.3
@#define prior_value_sigA = 0.3
@#define prior_value_sigX = 0.15
@#define prior_value_Xss = 1.1
@#define prior_value_dumA = 2
@#define prior_value_dumK = 1
@#define prior_value_dumepsA = 5
@#define prior_value_dumepsX = 6
var Y X K C A W Z; %note that declaration order is different from DR order on purpose
varobs C K A;
varexo epsA epsX epsY;
parameters alph betta rhoA sigA sigX Xss dumA dumK dumepsA dumepsX;
%Calibration
alph = @{value_alph};
betta = @{value_betta};
rhoA = @{value_rhoA};
sigA = @{value_sigA};
sigX = @{value_sigX};
Xss = @{value_Xss};
dumA = @{value_dumA};
dumK = @{value_dumK};
dumepsA = @{value_dumepsA};
dumepsX = @{value_dumepsX};
model;
%this is original Brock-Mirman model equations with AR(1) shock
C^(-1) = alph*betta*C(+1)^(-1)*A(+1)*K^(alph-1);
K = A*K(-1)^alph - C;
log(A) = rhoA*log(A(-1)) + sigA*epsA;
Y = A*K(-1)^alph + epsY;
%augmented auxiliary equations to get nonzero ghs2, ghxss and ghuss, they have no economic meaning
log(X) = log(Xss) + sigX*epsX;
W = X(+1)*exp(dumA*A(-1))*exp(dumK*K(-1)); %this will get a nonzero ghs2 and ghxss
Z = X(+1)*exp(dumepsA*epsA)*exp(dumepsX*epsX); %this will get a nonzero ghs2 and ghuss
end;
shocks;
var epsA = (@{value_SE_A})^2;
var epsX = (@{value_SE_X})^2;
var epsA, epsX = @{value_rho_AX}*@{value_SE_A}*@{value_SE_X};
var epsY = (@{value_SE_Y})^2;
end;
steady_state_model;
X = Xss;
A = 1;
K = (alph*betta*A)^(1/(1-alph));
C = A*K^alph-K;
Y = A*K^alph;
W = Xss*exp(dumA*A)*exp(dumK*K);
Z = Xss*exp(dumepsA*0)*exp(dumepsX*0);
end;
estimated_params;%note that the parameter ordering is different than declaration order
rhoA, normal_pdf, @{prior_value_rhoA}, 0.1;
betta, normal_pdf, @{prior_value_betta}, 0.1;
alph, normal_pdf, @{prior_value_alph}, 0.1;
corr epsA, epsX, normal_pdf, @{prior_value_rho_AX}, 0.1;
sigA, normal_pdf, @{prior_value_sigA}, 0.1;
stderr epsA, normal_pdf, @{prior_value_SE_A}, 0.1;
stderr epsX, normal_pdf, @{prior_value_SE_X}, 0.1;
stderr epsY, normal_pdf, @{prior_value_SE_Y}, 0.1;
sigX, normal_pdf, @{prior_value_sigX}, 0.1;
Xss, normal_pdf, @{prior_value_Xss}, 0.1;
dumepsX, normal_pdf, @{prior_value_dumepsX}, 0.1;
dumK, normal_pdf, @{prior_value_dumK}, 0.1;
dumepsA, normal_pdf, @{prior_value_dumepsA}, 0.1;
dumA, normal_pdf, @{prior_value_dumA}, 0.1;
end;
%save calibration parameters as these get overwritten through stoch_simul and identification command
calib_params = M_.params;
calib_Sigma_e = M_.Sigma_e;
stoch_simul(order=@{ORDER},nograph,irf=0,periods=0);
identification(order=@{ORDER},nograph,no_identification_strength);
indpmodel = estim_params_.param_vals(:,1);
indpstderr = estim_params_.var_exo(:,1);
indpcorr = estim_params_.corrx(:,1:2);
[I,~] = find(M_.lead_lag_incidence');
%% Parameter derivatives of perturbation
@#if CREATE_SYMBOLIC == 1
syms Y_ Y0 Yp C_ C0 Cp K_ K0 Kp A_ A0 Ap X_ X0 Xp W_ W0 Wp Z_ Z0 Zp;
syms epsA0 epsX0 epsY0;
syms RHO_AX SE_A SE_X SE_Y ALPH BETTA RHOA SIGA SIGX XSS DUMA DUMK DUMEPSA DUMEPSX;
syms sig;
SYM.corr_params = [RHO_AX];
SYM.stderr_params_decl = [SE_A SE_X SE_Y];
SYM.modparams_decl = [ALPH BETTA RHOA SIGA SIGX XSS DUMA DUMK DUMEPSA DUMEPSX];
SYM.endo_decl_ = [Y_ X_ K_ C_ A_ W_ Z_];
SYM.endo_decl0 = [Y0 X0 K0 C0 A0 W0 Z0];
SYM.endo_declp = [Yp Xp Kp Cp Ap Wp Zp];
SYM.ex0 = [epsA0 epsX0 epsY0];
SYM.model_eqs = [C0^(-1)-ALPH*BETTA*Cp^(-1)*Ap*K0^(ALPH-1);
K0-A0*K_^ALPH+C0;
log(A0)-RHOA*log(A_)-SIGA*epsA0;
Y0 - A0*K_^ALPH - epsY0;
log(X0)-log(XSS)-SIGX*epsX0;
W0 - Xp*exp(DUMA*A_)*exp(DUMK*K_);
Z0 - Xp*exp(DUMEPSA*epsA0)*exp(DUMEPSX*epsX0);
];
SYM.Sigma_e = [SE_A^2 RHO_AX*SE_A*SE_X 0; RHO_AX*SE_A*SE_X SE_X^2 0; 0 0 SE_Y^2];
SYM.Correlation_matrix = [1 RHO_AX 0; RHO_AX 1 0; 0 0 1];
%steady states
SYM.epsAbar = sym(0);
SYM.epsXbar = sym(0);
SYM.epsYbar = sym(0);
SYM.Abar = sym(1);
SYM.Xbar = XSS;
SYM.Kbar = (ALPH*BETTA*SYM.Abar)^(1/(1-ALPH));
SYM.Cbar = SYM.Abar*SYM.Kbar^ALPH-SYM.Kbar;
SYM.Ybar = SYM.Abar*SYM.Kbar^ALPH;
SYM.Wbar = XSS*exp(DUMA*SYM.Abar)*exp(DUMK*SYM.Kbar);
SYM.Zbar = XSS*exp(DUMEPSA*0)*exp(DUMEPSX*0);
SYM.endo_bar_decl = [SYM.Ybar SYM.Xbar SYM.Kbar SYM.Cbar SYM.Abar SYM.Wbar SYM.Zbar];
SYM.ex0bar = [SYM.epsAbar SYM.epsXbar SYM.epsYbar];
%True policy functions (in declaration order): gs is used for derivatives wrt to perturbation parameter
SYM.g = [A_^RHOA*K_^ALPH*exp(SIGA*epsA0) + epsY0; %Y
XSS*exp(SIGX*epsX0); %X
ALPH*BETTA*A_^RHOA*K_^ALPH*exp(SIGA*epsA0); %K
(1-ALPH*BETTA)*A_^RHOA*K_^ALPH*exp(SIGA*epsA0); %C
A_^RHOA*exp(SIGA*epsA0); %A
XSS*exp(SIGX*0)*exp(DUMA*A_)*exp(DUMK*K_); %W
XSS*exp(SIGX*0)*exp(DUMEPSA*epsA0)*exp(DUMEPSX*epsX0); %Z
];
SYM.gs = [A_^RHOA*K_^ALPH*exp(SIGA*sig*epsA0) + sig*epsY0; %Y
XSS*exp(SIGX*sig*epsX0); %X
ALPH*BETTA*A_^RHOA*K_^ALPH*exp(SIGA*sig*epsA0); %K
(1-ALPH*BETTA)*A_^RHOA*K_^ALPH*exp(SIGA*sig*epsA0); %C
A_^RHOA*exp(SIGA*sig*epsA0); %A
XSS*(1+1/2*SIGX^2*sig^2*SE_X^2)*exp(DUMA*A_)*exp(DUMK*K_); %W
XSS*(1+1/2*SIGX^2*sig^2*SE_X^2)*exp(DUMEPSA*epsA0)*exp(DUMEPSX*epsX0); %Z
];
%put into in DR-order
SYM.g = SYM.g(oo_.dr.order_var);
SYM.gs = SYM.gs(oo_.dr.order_var);
SYM.endo_DR_ = SYM.endo_decl_(oo_.dr.order_var);
SYM.endo_DR0 = SYM.endo_decl0(oo_.dr.order_var);
SYM.endo_DRp = SYM.endo_declp(oo_.dr.order_var);
SYM.Yss = SYM.endo_bar_decl(oo_.dr.order_var);
SYM.yy0 = [SYM.endo_decl_(M_.lead_lag_incidence(1,:)~=0) SYM.endo_decl0(M_.lead_lag_incidence(2,:)~=0) SYM.endo_declp(M_.lead_lag_incidence(3,:)~=0)];
SYM.yy0ex0 = [SYM.yy0 SYM.ex0];
SYM.yy0ex0bar = [SYM.endo_bar_decl(I) SYM.ex0bar];
SYM.x = SYM.endo_DR_(M_.nstatic + (1:M_.nspred));
SYM.xbar = SYM.Yss(M_.nstatic + (1:M_.nspred));
SYM.stderr_params = SYM.stderr_params_decl(indpstderr);
SYM.modparams = SYM.modparams_decl(indpmodel);
SYM.params = [SYM.stderr_params SYM.corr_params SYM.modparams];
SYM.yy0ex0_nbr = length(SYM.yy0ex0);
%% Parameter derivatives
SYM.dYss = jacobian(SYM.Yss,SYM.modparams);
SYM.d2Yss = jacobian(SYM.dYss(:),SYM.modparams);
SYM.g1 = jacobian(SYM.model_eqs,SYM.yy0ex0);
SYM.g2 = reshape(jacobian(vec(SYM.g1),SYM.yy0ex0),M_.endo_nbr,SYM.yy0ex0_nbr^2);
for j = 1:M_.endo_nbr
SYM.g3(j,:) = vec(transpose(jacobian(vec(SYM.g2(j,:)),SYM.yy0ex0))); %dynare ordering
end
SYM.dg1 = jacobian(vec(subs(SYM.g1,SYM.yy0ex0,SYM.yy0ex0bar)),SYM.modparams);
SYM.dg2 = jacobian(vec(subs(SYM.g2,SYM.yy0ex0,SYM.yy0ex0bar)),SYM.modparams);
SYM.dg3 = jacobian(vec(subs(SYM.g3,SYM.yy0ex0,SYM.yy0ex0bar)),SYM.modparams);
SYM.dSigma_e = jacobian(SYM.Sigma_e(:),SYM.params);
SYM.dCorrelation_matrix = jacobian(SYM.Correlation_matrix(:),SYM.params);
SYM.ghx = jacobian(SYM.g,SYM.x);
SYM.KalmanA = sym(zeros(M_.endo_nbr,M_.endo_nbr));
SYM.KalmanA(:,M_.nstatic + (1:M_.nspred)) = SYM.ghx;
SYM.dghx = jacobian(vec(subs(SYM.ghx, [SYM.x SYM.ex0], [SYM.xbar SYM.ex0bar])),SYM.params);
SYM.dKalmanA = jacobian(vec(subs(SYM.KalmanA, [SYM.x SYM.ex0], [SYM.xbar SYM.ex0bar])),SYM.params);
SYM.d2KalmanA = jacobian(SYM.dKalmanA(:),SYM.params);
SYM.ghu = jacobian(SYM.g,SYM.ex0);
SYM.dghu = jacobian(vec(subs(SYM.ghu, [SYM.x SYM.ex0], [SYM.xbar SYM.ex0bar])),SYM.params);
SYM.Om = SYM.ghu*SYM.Sigma_e*transpose(SYM.ghu);
SYM.dOm = jacobian(vec(subs(SYM.Om,[SYM.x SYM.ex0], [SYM.xbar SYM.ex0bar])),SYM.params);
SYM.d2Om = jacobian(SYM.dOm(:),SYM.params);
SYM.ghxx = jacobian(SYM.ghx(:),SYM.x);
SYM.ghxx = reshape(SYM.ghxx,M_.endo_nbr,M_.nspred^2);
SYM.dghxx = jacobian(vec(subs(SYM.ghxx, [SYM.x SYM.ex0], [SYM.xbar SYM.ex0bar])),SYM.params);
SYM.ghxu = jacobian(SYM.ghu(:),SYM.x);
SYM.ghxu = reshape(SYM.ghxu,M_.endo_nbr,M_.nspred*M_.exo_nbr);
SYM.dghxu = jacobian(vec(subs(SYM.ghxu, [SYM.x SYM.ex0], [SYM.xbar SYM.ex0bar])),SYM.params);
SYM.ghuu = jacobian(SYM.ghu(:),SYM.ex0);
SYM.ghuu = reshape(SYM.ghuu,M_.endo_nbr,M_.exo_nbr*M_.exo_nbr);
SYM.dghuu = jacobian(vec(subs(SYM.ghuu, [SYM.x SYM.ex0], [SYM.xbar SYM.ex0bar])),SYM.params);
SYM.ghs2 = jacobian(jacobian(SYM.gs,sig),sig);
SYM.dghs2 = jacobian(vec(subs(SYM.ghs2, [SYM.x SYM.ex0], [SYM.xbar SYM.ex0bar])),SYM.params);
for j = 1:M_.endo_nbr
SYM.ghxxx(j,:) = vec(transpose(jacobian(vec(SYM.ghxx(j,:)),SYM.x))); %dynare ordering
SYM.ghxxu(j,:) = vec(transpose(jacobian(vec(SYM.ghxx(j,:)),SYM.ex0))); %dynare ordering
SYM.ghxuu(j,:) = vec(transpose(jacobian(vec(SYM.ghxu(j,:)),SYM.ex0))); %dynare ordering
SYM.ghuuu(j,:) = vec(transpose(jacobian(vec(SYM.ghuu(j,:)),SYM.ex0))); %dynare ordering
SYM.ghxss(j,:) = vec(transpose(jacobian(vec(SYM.ghs2(j,:)),SYM.x))); %dynare ordering
SYM.ghuss(j,:) = vec(transpose(jacobian(vec(SYM.ghs2(j,:)),SYM.ex0))); %dynare ordering
end
SYM.dghxxx = jacobian(vec(subs(SYM.ghxxx, [SYM.x SYM.ex0], [SYM.xbar SYM.ex0bar])),SYM.params);
SYM.dghxxu = jacobian(vec(subs(SYM.ghxxu, [SYM.x SYM.ex0], [SYM.xbar SYM.ex0bar])),SYM.params);
SYM.dghxuu = jacobian(vec(subs(SYM.ghxuu, [SYM.x SYM.ex0], [SYM.xbar SYM.ex0bar])),SYM.params);
SYM.dghuuu = jacobian(vec(subs(SYM.ghuuu, [SYM.x SYM.ex0], [SYM.xbar SYM.ex0bar])),SYM.params);
SYM.dghxss = jacobian(vec(subs(SYM.ghxss, [SYM.x SYM.ex0], [SYM.xbar SYM.ex0bar])),SYM.params);
SYM.dghuss = jacobian(vec(subs(SYM.ghuss, [SYM.x SYM.ex0], [SYM.xbar SYM.ex0bar])),SYM.params);
% Evaluate numerically
for jj = 1:2
if jj == 1
RHO_AX = @{prior_value_rho_AX};
SE_A = @{prior_value_SE_A};
SE_X = @{prior_value_SE_X};
SE_Y = @{prior_value_SE_Y};
ALPH = @{prior_value_alph};
BETTA = @{prior_value_betta};
RHOA = @{prior_value_rhoA};
SIGA = @{prior_value_sigA};
SIGX = @{prior_value_sigX};
XSS = @{prior_value_Xss};
DUMA = @{prior_value_dumA};
DUMK = @{prior_value_dumK};
DUMEPSA = @{prior_value_dumepsA};
DUMEPSX = @{prior_value_dumepsX};
elseif jj == 2
RHO_AX = @{value_rho_AX};
SE_A = @{value_SE_A};
SE_X = @{value_SE_X};
SE_Y = @{value_SE_Y};
ALPH = @{value_alph};
BETTA = @{value_betta};
RHOA = @{value_rhoA};
SIGA = @{value_sigA};
SIGX = @{value_sigX};
XSS = @{value_Xss};
DUMA = @{value_dumA};
DUMK = @{value_dumK};
DUMEPSA = @{value_dumepsA};
DUMEPSX = @{value_dumepsX};
end
sig = 1;
nSYM.Yss = transpose(double(subs(SYM.Yss)));
nSYM.dYss = reshape(double(subs(SYM.dYss)), M_.endo_nbr, length(SYM.modparams));
nSYM.d2Yss = double(reshape(subs(SYM.d2Yss), [M_.endo_nbr length(SYM.modparams) length(SYM.modparams)]));
nSYM.g1 = double(subs(subs(SYM.g1, SYM.yy0ex0, SYM.yy0ex0bar)));
nSYM.dg1 = reshape(double(subs(SYM.dg1)),M_.endo_nbr, SYM.yy0ex0_nbr, length(SYM.modparams));
nSYM.g2 = sparse(double(subs(subs(SYM.g2, SYM.yy0ex0, SYM.yy0ex0bar))));
nSYM.dg2 = reshape(sparse(double(subs(SYM.dg2))), M_.endo_nbr, SYM.yy0ex0_nbr^2*length(SYM.modparams));
nSYM.g3 = sparse(double(subs(subs(SYM.g3, SYM.yy0ex0, SYM.yy0ex0bar))));
nSYM.dg3 = reshape(sparse(double(subs(SYM.dg3))), M_.endo_nbr, SYM.yy0ex0_nbr^3*length(SYM.modparams));
nSYM.Sigma_e = double(subs(SYM.Sigma_e));
nSYM.dSigma_e = double(reshape(subs(SYM.dSigma_e), M_.exo_nbr, M_.exo_nbr,length(SYM.params)));
nSYM.Correlation_matrix = double(subs(SYM.Correlation_matrix));
nSYM.dCorrelation_matrix = double(reshape(subs(SYM.dCorrelation_matrix), M_.exo_nbr, M_.exo_nbr,length(SYM.params)));
nSYM.ghx = double(subs(subs(SYM.ghx, SYM.yy0ex0, SYM.yy0ex0bar)));
nSYM.dghx = reshape(double(subs(SYM.dghx)), M_.endo_nbr, M_.nspred, length(SYM.params));
nSYM.ghu = double(subs(subs(SYM.ghu, SYM.yy0ex0, SYM.yy0ex0bar)));
nSYM.dghu = reshape(double(subs(SYM.dghu)), M_.endo_nbr, M_.exo_nbr, length(SYM.params));
nSYM.Om = double(subs(subs(SYM.Om, SYM.yy0ex0, SYM.yy0ex0bar)));
nSYM.dOm = reshape(double(subs(SYM.dOm)), M_.endo_nbr, M_.endo_nbr, length(SYM.params));
nSYM.d2Om = reshape(sparse(double(subs(SYM.d2Om))), M_.endo_nbr^2, length(SYM.params)^2);
SYM.indp2 = reshape(1:length(SYM.params)^2,length(SYM.params),length(SYM.params));
SYM.indx2 = reshape(1:M_.endo_nbr^2, M_.endo_nbr, M_.endo_nbr);
nSYM.d2Om = nSYM.d2Om(dyn_vech(SYM.indx2),dyn_vech(SYM.indp2));
nSYM.KalmanA = double(subs(subs(SYM.KalmanA, SYM.yy0ex0, SYM.yy0ex0bar)));
nSYM.dKalmanA = reshape(double(subs(SYM.dKalmanA)), M_.endo_nbr, M_.endo_nbr, length(SYM.params));
nSYM.d2KalmanA = reshape(double(subs(SYM.d2KalmanA)), M_.endo_nbr^2, length(SYM.params)^2);
nSYM.d2KalmanA = nSYM.d2KalmanA(:,dyn_vech(SYM.indp2));
nSYM.ghxx = double(subs(subs(SYM.ghxx, SYM.yy0ex0, SYM.yy0ex0bar)));
nSYM.dghxx = reshape(double(subs(SYM.dghxx)), M_.endo_nbr, M_.nspred^2, length(SYM.params));
nSYM.ghxu = double(subs(subs(SYM.ghxu, SYM.yy0ex0, SYM.yy0ex0bar)));
nSYM.dghxu = reshape(double(subs(SYM.dghxu)), M_.endo_nbr, M_.nspred*M_.exo_nbr, length(SYM.params));
nSYM.ghuu = double(subs(subs(SYM.ghuu, SYM.yy0ex0, SYM.yy0ex0bar)));
nSYM.dghuu = reshape(double(subs(SYM.dghuu)), M_.endo_nbr, M_.exo_nbr^2, length(SYM.params));
nSYM.ghs2 = double(subs(subs(SYM.ghs2, SYM.yy0ex0, SYM.yy0ex0bar)));
nSYM.dghs2 = reshape(double(subs(SYM.dghs2)), M_.endo_nbr, length(SYM.params));
nSYM.ghxxx = double(subs(subs(SYM.ghxxx, SYM.yy0ex0, SYM.yy0ex0bar)));
nSYM.dghxxx = reshape(double(subs(SYM.dghxxx)), M_.endo_nbr, M_.nspred^3, length(SYM.params));
nSYM.ghxxu = double(subs(subs(SYM.ghxxu, SYM.yy0ex0, SYM.yy0ex0bar)));
nSYM.dghxxu = reshape(double(subs(SYM.dghxxu)), M_.endo_nbr, M_.nspred^2*M_.exo_nbr, length(SYM.params));
nSYM.ghxuu = double(subs(subs(SYM.ghxuu, SYM.yy0ex0, SYM.yy0ex0bar)));
nSYM.dghxuu = reshape(double(subs(SYM.dghxuu)), M_.endo_nbr, M_.nspred*M_.exo_nbr^2, length(SYM.params));
nSYM.ghuuu = double(subs(subs(SYM.ghuuu, SYM.yy0ex0, SYM.yy0ex0bar)));
nSYM.dghuuu = reshape(double(subs(SYM.dghuuu)), M_.endo_nbr, M_.exo_nbr^3, length(SYM.params));
nSYM.ghxss = double(subs(subs(SYM.ghxss, SYM.yy0ex0, SYM.yy0ex0bar)));
nSYM.dghxss = reshape(double(subs(SYM.dghxss)), M_.endo_nbr, M_.nspred, length(SYM.params));
nSYM.ghuss = double(subs(subs(SYM.ghuss, SYM.yy0ex0, SYM.yy0ex0bar)));
nSYM.dghuss = reshape(double(subs(SYM.dghuss)), M_.endo_nbr, M_.exo_nbr, length(SYM.params));
if jj == 1
nSYMprior = nSYM;
save('nBrockMirmanSYM.mat','nSYMprior')
elseif jj==2
nSYMcalib = nSYM;
save('nBrockMirmanSYM.mat','nSYMcalib','-append')
end
end
@#endif
clc;
tol_vars.Yss = 1e-14;
tol_vars.Sigma_e = 1e-15;
tol_vars.Correlation_matrix = 1e-16;
tol_vars.g1 = 1e-13;
tol_vars.ghx = 1e-13;
tol_vars.ghu = 1e-13;
@#if ORDER > 1
tol_vars.g2 = 1e-12;
tol_vars.ghxx = 1e-12;
tol_vars.ghxu = 1e-12;
tol_vars.ghuu = 1e-12;
tol_vars.ghs2 = 1e-12;
@#endif
@#if ORDER > 2
tol_vars.g3 = 1e-11;
tol_vars.ghxxx = 1e-11;
tol_vars.ghxxu = 1e-11;
tol_vars.ghxuu = 1e-11;
tol_vars.ghuuu = 1e-11;
tol_vars.ghxss = 1e-11;
tol_vars.ghuss = 1e-11;
@#endif
tol_dvars.dYss = [1e-9 1e-9 1e-14 1e-14];
tol_dvars.dSigma_e = [1e-9 1e-15 1e-14 1e-14];
tol_dvars.dCorrelation_matrix = [1e-9 1e-15 1e-14 1e-14];
tol_dvars.dg1 = [1e-6 1e-6 1e-13 1e-13];
tol_dvars.dghx = [1e-8 1e-8 1e-13 1e-13];
tol_dvars.dghu = [1e-8 1e-8 1e-13 1e-13];
@#if ORDER > 1
tol_dvars.dg2 = [1e-5 1e-5 1e-11 1e-11];
tol_dvars.dghxx = [1e-6 1e-6 1e-12 1e-12];
tol_dvars.dghxu = [1e-6 1e-6 1e-12 1e-12];
tol_dvars.dghuu = [1e-6 1e-6 1e-12 1e-12];
tol_dvars.dghs2 = [1e-6 1e-6 1e-12 1e-12];
@#endif
@#if ORDER > 2
tol_dvars.dg3 = [1e-3 1e-3 1e-9 1e-9];
tol_dvars.dghxxx = [1e-5 1e-5 1e-11 1e-11];
tol_dvars.dghxxu = [1e-5 1e-5 1e-11 1e-11];
tol_dvars.dghxuu = [1e-5 1e-5 1e-11 1e-11];
tol_dvars.dghuuu = [1e-5 1e-5 1e-11 1e-11];
tol_dvars.dghxss = [1e-5 1e-5 1e-11 1e-11];
tol_dvars.dghuss = [1e-5 1e-5 1e-11 1e-11];
@#endif
tol_dvars.d2KalmanA = [1e-3 1e-3 1e-13 1e-13];
tol_dvars.d2Om = [1e-3 1e-3 1e-13 1e-13];
tol_dvars.d2Yss = [1e-3 1e-3 1e-13 1e-13];
options_.dynatol.x = eps.^(1/3); %set numerical differentiation step in fjaco.m
for jj = 1:2
lst_vars = {'Yss', 'Sigma_e', 'Correlation_matrix','g1','ghx','ghu'};
@#if ORDER > 1
lst_vars = [lst_vars, 'g2','ghxx','ghxu','ghuu','ghs2'];
@#endif
@#if ORDER > 2
lst_vars = [lst_vars, 'g3','ghxxx','ghxxu','ghxuu','ghuuu','ghxss','ghuss'];
@#endif
lst_dvars = {'dYss','dSigma_e','dg1','dghx','dghu'};
@#if ORDER > 1
lst_dvars = [lst_dvars, 'dg2','dghxx','dghxu','dghuu','dghs2'];
@#endif
@#if ORDER > 2
lst_dvars = [lst_dvars, 'dg3','dghxxx','dghxxu','dghxuu','dghuuu','dghxss','dghuss'];
@#endif
load('nBrockMirmanSYM.mat');
if jj==1
strparamset = 'PRIOR';
nSYM = nSYMprior;
xparam_prior = set_prior(estim_params_,M_,options_);
M_ = set_all_parameters(xparam_prior,estim_params_,M_);
elseif jj==2
strparamset = 'CALIBRATION';
nSYM = nSYMcalib;
xparam1_calib = [];
for j = 1:length(indpstderr)
xparam1_calib = [xparam1_calib; sqrt(calib_Sigma_e(j,j))];
end
for j = 1:size(indpcorr,1)
xparam1_calib = [xparam1_calib; calib_Sigma_e(indpcorr(j,1),indpcorr(j,2))/( sqrt(calib_Sigma_e(indpcorr(j,1),indpcorr(j,1))) * sqrt(calib_Sigma_e(indpcorr(j,2),indpcorr(j,2))) )];
end
xparam1_calib = [xparam1_calib; calib_params(indpmodel)];
M_ = set_all_parameters(xparam1_calib,estim_params_,M_);
end
[~,info,M_,options_,oo_] = resol(0,M_, options_, oo_);
%For convenience we save the objects to compare into oo_.dr.
oo_.dr.Yss = oo_.dr.ys(oo_.dr.order_var);
oo_.dr.Sigma_e = M_.Sigma_e;
oo_.dr.Correlation_matrix = M_.Correlation_matrix;
ex0 = oo_.exo_steady_state';
[~, oo_.dr.g1, oo_.dr.g2, oo_.dr.g3] = feval([M_.fname,'.dynamic'], oo_.dr.ys(I), oo_.exo_steady_state', M_.params, oo_.dr.ys, 1);
oo_.dr.g3 = unfold_g3(oo_.dr.g3, length(oo_.dr.ys(I))+length(oo_.exo_steady_state')); %add symmetric elements to g3
fprintf('***** %s: SOME COMMON OBJECTS *****\n', strparamset)
for id_var = 1:size(lst_vars,2)
dx = norm( nSYM.(sprintf('%s',lst_vars{id_var})) - oo_.dr.(sprintf('%s',lst_vars{id_var})), Inf);
fprintf('Max absolute deviation for %s: %e\n', lst_vars{id_var}, dx);
if dx > tol_vars.(sprintf('%s',lst_vars{id_var}))
error('Something wrong in steady state computation, solution algorithm or preprocessor')
end
end
for d2flag = [0 1];
if d2flag
lst_dvars = [lst_dvars {'d2KalmanA', 'd2Om', 'd2Yss'}];
end
KRONFLAG = [-1 -2 0 1];
for id_kronflag = 1:length(KRONFLAG)
fprintf('***** %s: d2flag=%d and kronflag=%d *****\n',strparamset, d2flag,KRONFLAG(id_kronflag))
options_.analytic_derivation_mode = KRONFLAG(id_kronflag);
DERIVS = get_perturbation_params_derivs(M_, options_, estim_params_, oo_, indpmodel, indpstderr, indpcorr, d2flag);
for id_var = 1:size(lst_dvars,2)
dx = norm( vec(nSYM.(sprintf('%s',lst_dvars{id_var}))) - vec(DERIVS.(sprintf('%s',lst_dvars{id_var}))), Inf);
fprintf('Max absolute deviation for %s: %e\n', lst_dvars{id_var}, dx);
if dx > tol_dvars.(sprintf('%s',lst_dvars{id_var}))(id_kronflag)
error('Something wrong in get_perturbation_params_derivs.m')
end
end
end
end
end

230
tests/analytic_derivatives/burnside_3_order_PertParamsDerivs.mod Normal file
View File

@ -0,0 +1,230 @@
% =========================================================================
% 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/>.
% =========================================================================
/*
Check the policy functions obtained by perturbation at a high approximation
order, using the Burnside (1998, JEDC) model (for which the analytical form of
the policy function is known).
As shown by Burnside, the policy function for yₜ is:
yₜ = βⁱ exp[aᵢ+bᵢ(xₜ−xₛₛ)]
where:
θ² ⎛ 2ρ 1−ρ²ⁱ⎞
— aᵢ = iθxₛₛ + σ² ─────── ⎢i − ────(1−ρⁱ) + ρ² ─────⎥
2(1−ρ)² ⎝ 1−ρ 1−ρ² ⎠
θρ
— bᵢ = ───(1−ρⁱ)
1−ρ
— xₛₛ is the steady state of x
— σ is the standard deviation of e.
With some algebra, it can be shown that the derivative of yₜ at the deterministic
steady state is equal to:
∂ᵐ⁺ⁿ⁺²ᵖ yₜ ∞ (2p)!
──────────────── = ∑ βⁱ bᵢᵐ⁺ⁿ ρᵐ ───── cᵢᵖ exp(iθxₛₛ)
∂ᵐxₜ₋₁ ∂ⁿeₜ ∂²ᵖs ⁱ⁼¹ p!
where:
— s is the stochastic scale factor
θ² ⎛ 2ρ 1−ρ²ⁱ⎞
— cᵢ = ─────── ⎢i − ────(1−ρⁱ) + ρ² ─────⎥
2(1−ρ)² ⎝ 1−ρ 1−ρ² ⎠
Note that derivatives with respect to an odd order for s (i.e. ∂²ᵖ⁺¹s) are always
equal to zero.
The policy function as returned in the oo_.dr.g_* matrices has the following properties:
— its elements are pre-multiplied by the Taylor coefficients;
— derivatives w.r.t. the stochastic scale factor have already been summed up;
— symmetric elements are folded (and they are not pre-multiplied by the number of repetitions).
As a consequence, the element gₘₙ corresponding to the m-th derivative w.r.t.
to xₜ₋₁ and the n-th derivative w.r.t. to eₜ is given by:
1 ∞ cᵢᵖ
gₘₙ = ────── ∑ ∑ βⁱ bᵢᵐ⁺ⁿ ρᵐ ──── exp(iθxₛₛ)
(m+n)! 0≤2p≤k-m-n ⁱ⁼¹ p!
where k is the order of approximation.
*/
@#define ORDER = 3
var y x;
varobs y;
varexo e;
parameters beta theta rho xbar;
xbar = 0.0179;
rho = -0.139;
theta = -1.5;
theta = -10;
beta = 0.95;
model;
y = beta*exp(theta*x(+1))*(1+y(+1));
x = (1-rho)*xbar + rho*x(-1)+e;
end;
shocks;
var e; stderr 0.0348;
end;
steady_state_model;
x = xbar;
y = beta*exp(theta*xbar)/(1-beta*exp(theta*xbar));
end;
estimated_params;
stderr e, normal_pdf, 0.0348,0.01;
beta, normal_pdf, 0.95, 0.01;
theta, normal_pdf, -10, 0.01;
rho, normal_pdf, -0.139, 0.01;
xbar, normal_pdf, 0.0179, 0.01;
end;
steady;check;model_diagnostics;
stoch_simul(order=@{ORDER},k_order_solver,irf=0,drop=0,periods=0,nograph);
identification(order=@{ORDER},nograph,no_identification_strength);
%make sure everything is computed at prior mean
xparam_prior = set_prior(estim_params_,M_,options_);
M_ = set_all_parameters(xparam_prior,estim_params_,M_);
[~,info,M_,options_,oo_] = resol(0,M_, options_, oo_);
indpmodel = estim_params_.param_vals(:,1);
indpstderr = estim_params_.var_exo(:,1);
indpcorr = estim_params_.corrx(:,1:2);
totparam_nbr = length(indpmodel) + length(indpstderr) + size(indpcorr,1);
%% Verify that the policy function coefficients are correct
i = 1:800;
SE_e=sqrt(M_.Sigma_e);
aux1 = rho*(1-rho.^i)/(1-rho);
aux2 = (1-rho.^(2*i))/(1-rho^2);
aux3 = 1/((1-rho)^2);
aux4 = (i-2*aux1+rho^2*aux2);
aux5=aux3*aux4;
b = theta*aux1;
c = 1/2*theta^2*SE_e^2*aux5;
%derivatives wrt to rho only
daux1_drho = zeros(1,length(i));
daux2_drho = zeros(1,length(i));
daux3_drho = 2/((1-rho)^3);
daux4_drho = zeros(1,length(i));
daux5_drho = zeros(1,length(i));
for ii = 1:length(i)
if ii == 1
daux1_drho(ii) = 1;
daux2_drho(ii) = 0;
else
daux1_drho(ii) = rho/(rho^2 - 2*rho + 1) - 1/(rho - 1) + ((ii+1)*rho^ii)/(rho - 1) - rho^(ii+1)/(rho^2 - 2*rho + 1);
daux2_drho(ii) = (2*rho)/(rho^4 - 2*rho^2 + 1) + (2*ii*rho^(2*ii-1))/(rho^2 - 1) - (2*rho^(2*ii+1))/(rho^4 - 2*rho^2 + 1);
end
daux4_drho(ii) = -2*daux1_drho(ii) + 2*rho*aux2(ii) + rho^2*daux2_drho(ii);
daux5_drho(ii) = daux3_drho*aux4(ii) + aux3*daux4_drho(ii);
end
%derivatives of b and c wrt to all parameters
db = zeros(size(b,1),size(b,2),M_.exo_nbr+M_.param_nbr);
db(:,:,3) = aux1;%wrt theta
db(:,:,4) = theta*daux1_drho;%wrt rho
dc = zeros(size(c,1),size(c,2),M_.exo_nbr+M_.param_nbr);
dc(:,:,1) = theta^2*SE_e*aux3*aux4;%wrt SE_e
dc(:,:,3) = theta*SE_e^2*aux3*aux4;%wrt theta
dc(:,:,4) = 1/2*theta^2*SE_e^2*daux5_drho; %wrt rho
d2flag=0;
g_0 = 1/2*oo_.dr.ghs2; if ~isequal(g_0,oo_.dr.g_0); error('something wrong'); end
g_1 = [oo_.dr.ghx oo_.dr.ghu] +3/6*[oo_.dr.ghxss oo_.dr.ghuss]; if ~isequal(g_1,oo_.dr.g_1); error('something wrong'); end
g_2 = 1/2*[oo_.dr.ghxx oo_.dr.ghxu oo_.dr.ghuu]; if ~isequal(g_2,oo_.dr.g_2); error('something wrong'); end
g_3 = 1/6*[oo_.dr.ghxxx oo_.dr.ghxxu oo_.dr.ghxuu oo_.dr.ghuuu]; if ~isequal(g_3,oo_.dr.g_3); error('something wrong'); end
tols = [1e-4 1e-4 1e-12 1e-12];
KRONFLAGS = [-1 -2 0 1];
for k = 1:length(KRONFLAGS)
fprintf('KRONFLAG=%d\n',KRONFLAGS(k));
options_.analytic_derivation_mode = KRONFLAGS(k);
DERIVS = get_perturbation_params_derivs(M_, options_, estim_params_, oo_, indpmodel, indpstderr, indpcorr, d2flag);
oo_.dr.dg_0 = permute(1/2*DERIVS.dghs2,[1 3 2]);
oo_.dr.dg_1 = cat(2,DERIVS.dghx,DERIVS.dghu) + 3/6*cat(2,DERIVS.dghxss,DERIVS.dghuss);
oo_.dr.dg_2 = 1/2*cat(2,DERIVS.dghxx,DERIVS.dghxu,DERIVS.dghuu);
oo_.dr.dg_3 = 1/6*[DERIVS.dghxxx DERIVS.dghxxu DERIVS.dghxuu DERIVS.dghuuu];
for ord = 0:@{ORDER}
g = oo_.dr.(['g_' num2str(ord)])(2,:); % Retrieve computed policy function for variable y
dg = oo_.dr.(['dg_' num2str(ord)])(2,:,:);
for m = 0:ord % m is the derivation order with respect to x(-1)
v = 0;
dv = zeros(1,M_.exo_nbr + M_.param_nbr);
for p = 0:floor((@{ORDER}-ord)/2) % 2p is the derivation order with respect to s
if ord+2*p > 0 % Skip the deterministic steady state constant
v = v + sum(beta.^i.*exp(theta*xbar*i).*b.^ord.*rho^m.*c.^p)/factorial(ord)/factorial(p);
%derivatives
dv(:,1) = dv(:,1) + sum( beta.^i.*exp(theta*xbar*i).*ord.*b.^(ord-1).*db(:,:,1).*rho^m.*c.^p...
+beta.^i.*exp(theta*xbar*i).*b.^ord.*rho^m.*p.*c.^(p-1).*dc(:,:,1)...
)/factorial(ord)/factorial(p);%wrt SE_E
dv(:,2) = dv(:,2) + sum( i.*beta.^(i-1).*exp(theta*xbar*i).*b.^ord.*rho^m.*c.^p...
+beta.^i.*exp(theta*xbar*i).*ord.*b.^(ord-1).*db(:,:,2).*rho^m.*c.^p...
+beta.^i.*exp(theta*xbar*i).*b.^ord.*rho^m.*p.*c.^(p-1).*dc(:,:,2)...
)/factorial(ord)/factorial(p);%wrt beta
dv(:,3) = dv(:,3) + sum( beta.^i.*exp(theta*xbar*i).*xbar.*i.*b.^ord.*rho^m.*c.^p...
+beta.^i.*exp(theta*xbar*i).*ord.*b.^(ord-1).*db(:,:,3).*rho^m.*c.^p...
+beta.^i.*exp(theta*xbar*i).*b.^ord.*rho^m.*p.*c.^(p-1).*dc(:,:,3)...
)/factorial(ord)/factorial(p);%wrt theta
dv(:,4) = dv(:,4) + sum( beta.^i.*exp(theta*xbar*i).*b.^ord.*m.*rho^(m-1).*c.^p...
+beta.^i.*exp(theta*xbar*i).*ord.*b.^(ord-1).*db(:,:,4).*rho^m.*c.^p...
+beta.^i.*exp(theta*xbar*i).*b.^ord.*rho^m.*p.*c.^(p-1).*dc(:,:,4)...
)/factorial(ord)/factorial(p);%wrt rho
dv(:,5) = dv(:,5) + sum( beta.^i.*exp(theta*xbar*i).*theta.*i.*b.^ord.*rho^m.*c.^p...
+beta.^i.*exp(theta*xbar*i).*ord.*b.^(ord-1).*db(:,:,5).*rho^m.*c.^p...
+beta.^i.*exp(theta*xbar*i).*b.^ord.*rho^m.*p.*c.^(p-1).*dc(:,:,5)...
)/factorial(ord)/factorial(p);%wrt xbar
end
end
if abs(v-g(ord+1-m)) > 1e-14
error(['Error in matrix oo_.dr.g_' num2str(ord)])
end
chk_dg = squeeze(dg(:,ord+1-m,:))';
if isempty(indpstderr)
chk_dv = dv(:,M_.exo_nbr+indpmodel);
elseif isempty(indpmodel)
chk_dv = dv(:,1:M_.exo_nbr);
else
chk_dv = dv;
end
fprintf('Max absolute deviation for dg_%d(2,%d,:): %e\n',ord,ord+1-m,norm( chk_dv - chk_dg, Inf));
if norm( chk_dv - chk_dg, Inf) > tols(k)
error(['Error in matrix dg_' num2str(ord)])
chk_dv
chk_dg
end
end
end
fprintf('\n');
end

66
tests/analytic_derivatives/ls2003.mod
View File

@ -1,66 +0,0 @@
var y y_s R pie dq pie_s de A y_obs pie_obs R_obs;
varexo e_R e_q e_ys e_pies e_A;
parameters psi1 psi2 psi3 rho_R tau alpha rr k rho_q rho_A rho_ys rho_pies;
psi1 = 1.54;
psi2 = 0.25;
psi3 = 0.25;
rho_R = 0.5;
alpha = 0.3;
rr = 2.51;
k = 0.5;
tau = 0.5;
rho_q = 0.4;
rho_A = 0.2;
rho_ys = 0.9;
rho_pies = 0.7;
model(linear);
y = y(+1) - (tau +alpha*(2-alpha)*(1-tau))*(R-pie(+1))-alpha*(tau +alpha*(2-alpha)*(1-tau))*dq(+1) + alpha*(2-alpha)*((1-tau)/tau)*(y_s-y_s(+1))-A(+1);
pie = exp(-rr/400)*pie(+1)+alpha*exp(-rr/400)*dq(+1)-alpha*dq+(k/(tau+alpha*(2-alpha)*(1-tau)))*y+k*alpha*(2-alpha)*(1-tau)/(tau*(tau+alpha*(2-alpha)*(1-tau)))*y_s;
pie = de+(1-alpha)*dq+pie_s;
R = rho_R*R(-1)+(1-rho_R)*(psi1*pie+psi2*(y+alpha*(2-alpha)*((1-tau)/tau)*y_s)+psi3*de)+e_R;
dq = rho_q*dq(-1)+e_q;
y_s = rho_ys*y_s(-1)+e_ys;
pie_s = rho_pies*pie_s(-1)+e_pies;
A = rho_A*A(-1)+e_A;
y_obs = y-y(-1)+A;
pie_obs = 4*pie;
R_obs = 4*R;
end;
shocks;
var e_R = 1.25^2;
var e_q = 2.5^2;
var e_A = 1.89;
var e_ys = 1.89;
var e_pies = 1.89;
end;
varobs y_obs R_obs pie_obs dq de;
estimated_params;
psi1 , gamma_pdf,1.5,0.5;
psi2 , gamma_pdf,0.25,0.125;
psi3 , gamma_pdf,0.25,0.125;
rho_R ,beta_pdf,0.5,0.2;
alpha ,beta_pdf,0.3,0.1;
rr ,gamma_pdf,2.5,1;
k , gamma_pdf,0.5,0.25;
tau ,gamma_pdf,0.5,0.2;
rho_q ,beta_pdf,0.4,0.2;
rho_A ,beta_pdf,0.5,0.2;
rho_ys ,beta_pdf,0.8,0.1;
rho_pies,beta_pdf,0.7,0.15;
stderr e_R,inv_gamma_pdf,1.2533,0.6551;
stderr e_q,inv_gamma_pdf,2.5066,1.3103;
stderr e_A,inv_gamma_pdf,1.2533,0.6551;
stderr e_ys,inv_gamma_pdf,1.2533,0.6551;
stderr e_pies,inv_gamma_pdf,1.88,0.9827;
end;
estimation(datafile='../ls2003/data_ca1.m',first_obs=8,nobs=79,mh_nblocks=2,
prefilter=1,mh_jscale=0.5,mh_replic=0, mode_compute=1, analytic_derivation);

BIN
tests/analytic_derivatives/nBrockMirmanSYM.mat Normal file
View File

Binary file not shown.

363
tests/identification/LindeTrabandt/LindeTrabandt2019_xfail.mod Normal file
View File

@ -0,0 +1,363 @@
% Original model by: J. Linde, M. Trabandt (2019: Should We Use Linearized Models to Calculate Fiscal Multipliers?
% Journal of Applied Econometrics, 2018, 33: 937-965. http://dx.doi.org/10.1002/jae.2641
% This version has some additional dynamics for capital and investment
% Created by @wmutschl (Willi Mutschler, willi@mutschler.eu)
% =========================================================================
% Declare endogenous variables
% =========================================================================
var
% Staggered-price economy
c ${c}$ (long_name='consumption')
lam ${\lambda}$ (long_name='lagrange multiplier budget')
n ${n}$ (long_name='labor')
w ${w}$ (long_name='real wage')
r ${r}$ (long_name='interest rate')
pie ${\pi}$ (long_name='inflation')
y ${y}$ (long_name='output')
pstar ${p^\ast}$ (long_name='price dispersion')
vartheta ${\vartheta}$ (long_name='aggregate price index')
mc ${mc}$ (long_name='marginal costs')
s ${s}$ (long_name='auxiliary variable for nonlinear pricing 1')
f ${f}$ (long_name='auxiliary variable for nonlinear pricing 2')
a ${a}$ (long_name='auxiliary variable for nonlinear pricing 3')
ptilde ${\tilde{p}}$ (long_name='reoptimized price')
delta1 ${\Delta_1}$ (long_name='price dispersion 1')
delta2 ${\Delta_2}$ (long_name='price dispersion 2')
delta3 ${\Delta_3}$ (long_name='price dispersion 3')
b ${b}$ (long_name='bonds')
tau ${\tau}$ (long_name='lump-sum tax')
g
nu %consumption preference shock
% Flex-price economy
cpot ${c^{pot}}$ (long_name='flex-price consumption')
npot ${n^{pot}}$ (long_name='flex-price labor')
wpot ${w^{pot}}$ (long_name='flex-price real wage')
rrpot ${r^{pot}}$ (long_name='flex-price interest rate')
ypot ${y^{pot}}$ (long_name='flex-price output')
bpot ${b^{pot}}$ (long_name='flex-price bonds')
taupot ${\tau^{pot}}$ (long_name='flex-price lump-sum tax')
% Added variables for capital and investment
k ${k}$ (long_name='capital')
rk ${r^{K}}$ (long_name='rental rate on capital')
iv ${i}$ (long_name='investment')
q ${q}$ (long_name='Tobins Q')
kpot ${k^{pot}}$ (long_name='flex-price capital')
rkpot ${r^{K,pot}}$ (long_name='flex-price rental rate on capital')
ivpot ${i^{pot}}$ (long_name='flex-price investment')
;
% =========================================================================
% Declare observable variables
% =========================================================================
varobs c y;
% =========================================================================
% Declare exogenous variables
% =========================================================================
varexo
epsg
epsnu
;
% =========================================================================
% Declare parameters
% =========================================================================
parameters
ALPHA ${\alpha}$ (long_name='capital share')
BETA ${\beta}$ (long_name='discount rate')
PIEBAR ${\bar{\pi}}$ (long_name='target inflation rate')
IOTA ${\iota}$ (long_name='indexation')
CHI ${\chi}$ (long_name='inverse frisch elasticity')
THETAP ${\theta_p}$ (long_name='net markup')
XIP ${\xi_p}$ (long_name='Calvo probability')
PSI ${\psi}$ (long_name='Kimbal curvature')
GAMMAPIE ${\gamma_\pi}$ (long_name='Taylor rule parameter inflation')
GAMMAX ${\gamma_x}$ (long_name='Taylor rule parameter output')
DELTA ${\delta}$ (long_name='depreciation rate')
QBAR ${\bar{Q}}$ (long_name='steady state Q')
NUBAR ${\bar{\nu}}$ (long_name='steady state consumption preference shock')
BBAR %government debt to annualized output ratio
GBAR_o_YBAR % government consumption ratio
TAUBAR
TAUNBAR
PHI
RHONU
RHOG
SIGNU
SIGG
;
% =========================================================================
% Calibration
% =========================================================================
BETA = 0.995;
PIEBAR = 1 + 0.005;
ALPHA = 0.3;
NUBAR = 0;%0.01;
BBAR = 0;%0.6;
GBAR_o_YBAR = 0;%0.2;
TAUBAR = 0;
PHI = 0.0101;
IOTA = 0.5;
CHI = 1/0.4;
XIP = 2/3;
PSI = -12.2;
THETAP = 0.1;
GAMMAPIE = 1.5;
GAMMAX = 0.125;
QBAR=1;
DELTA = 0;%0.025;
RHONU = 0.95;
RHOG = 0.95;
SIGG = 0.6/100;
SIGNU = 0.6/100;
TAUNBAR = 0;
% ========================================================================
% Model equations
% =========================================================================
model;
% -------------------------------------------------------------------------
% Auxiliary expressions
% -------------------------------------------------------------------------
#OMEGAP = (1+THETAP)*(1+PSI)/(1+PSI+THETAP*PSI); %gross markup
#gammap = pie^IOTA*PIEBAR^(1-IOTA); %indexation rule
%#gammap = PIEBAR;
% -------------------------------------------------------------------------
% Added equations for capital and investment
% -------------------------------------------------------------------------
[name='foc household wrt i']
lam = lam*q;
[name='foc household wrt k']
lam*q = BETA*lam(+1)*(rk(+1)+ (1-DELTA)*q(+1) );
[name='capital accumulation']
k = (1-DELTA)*k(-1) + iv;
[name='capital demand']
rk=ALPHA*mc*y/k(-1);
[name='flex price foc household wrt k']
(cpot-steady_state(cpot)*nu)^(-1) = BETA*(cpot(+1)-steady_state(cpot)*nu(+1))^(-1)*(rkpot(+1) + (1-DELTA) );
[name='flex price capital accumulation']
kpot = (1-DELTA)*kpot(-1) + ivpot;
[name='flex price capital demand']
rkpot=ALPHA/(1+THETAP)*(npot/kpot(-1))^(1-ALPHA);
% -------------------------------------------------------------------------
% Actual model equations
% -------------------------------------------------------------------------
[name='foc household wrt c (marginal utility)']
(c-steady_state(c)*nu)^(-1) = lam;
[name='foc household wrt l (leisure vs. labor)']
n^CHI = (1-TAUNBAR)*lam*w;
[name='foc household wrt b (Euler equation)']
lam = BETA*r/pie(+1)*lam(+1);
[name='aggregate resource constraint']
y = c + g + iv;
[name='production function']
y=pstar^(-1)*k(-1)^ALPHA*n^(1-ALPHA);
[name='Nonlinear pricing 1']
s = OMEGAP*lam*y*vartheta^(OMEGAP/(OMEGAP-1))*mc + BETA*XIP*(gammap/pie(+1))^(OMEGAP/(1-OMEGAP))*s(+1);
[name='Nonlinear pricing 2']
f = lam*y*vartheta^(OMEGAP/(OMEGAP-1)) + BETA*XIP*(gammap/pie(+1))^(1/(1-OMEGAP))*f(+1);
[name='Nonlinear pricing 3']
a = PSI*(OMEGAP-1)*lam*y + BETA*XIP*(gammap/pie(+1))*a(+1);
[name='Nonlinear pricing 4']
s = f*ptilde -a*ptilde^(1+OMEGAP/(OMEGAP-1));
[name='zero profit condition']
vartheta = 1 + PSI - PSI*delta2;
[name='aggregate price index']
vartheta = delta3;
[name='overall price dispersion']
pstar = 1/(1+PSI)*vartheta^(OMEGAP/(OMEGAP-1))*delta1^(OMEGAP/(1-OMEGAP)) + PSI/(1+PSI);
[name='price dispersion 1']
delta1^(OMEGAP/(1-OMEGAP)) = (1-XIP)*ptilde^(OMEGAP/(1-OMEGAP)) + XIP*(gammap/pie*delta1(-1))^(OMEGAP/(1-OMEGAP));
[name='price dispersion 2']
delta2 = (1-XIP)*ptilde+XIP*gammap/pie*delta2(-1);
[name='price dispersion 3']
delta3^(1/(1-OMEGAP)) = (1-XIP)*ptilde^(1/(1-OMEGAP)) + XIP*(gammap/pie*delta3(-1))^(1/(1-OMEGAP));
[name='marginal costs']
(1-ALPHA)*mc = w*k(-1)^(-ALPHA)*n^ALPHA;
[name='Taylor Rule']
r = steady_state(r)*(pie/PIEBAR)^GAMMAPIE*( (y/steady_state(y))/(ypot/steady_state(ypot)) )^GAMMAX;
[name='government budget']
b = r(-1)/pie*b(-1) + g/steady_state(y) -TAUNBAR*w*n/steady_state(y) - tau;
[name='fiscal rule']
tau = TAUBAR + PHI*(b(-1)-BBAR);
[name='flex price Euler equation']
(cpot-steady_state(cpot)*nu)^(-1) = BETA*rrpot*(cpot(+1)-steady_state(cpot)*nu(+1))^(-1);
[name='flex price foc household wrt l (leisure vs. labor)']
npot^CHI = (1-TAUNBAR)*(cpot-steady_state(cpot)*nu)^(-1)*wpot;
[name='flex price wage']
(1-ALPHA)/(1+THETAP)*kpot(-1)^ALPHA = wpot*npot^ALPHA;
[name='flex price aggregate resource constraint']
ypot = cpot + g + ivpot;
[name='flex price production function']
ypot=kpot(-1)^ALPHA*npot^(1-ALPHA);
[name='flex price government budget']
bpot = rrpot(-1)*bpot(-1) + g/steady_state(y) -TAUNBAR*wpot*npot/steady_state(y) - taupot;
[name='fiscal rule ']
taupot = TAUBAR + PHI*(bpot(-1)-BBAR);
g/steady_state(y) - GBAR_o_YBAR = RHOG*(g(-1)/steady_state(y) - GBAR_o_YBAR) + epsg;
nu - NUBAR = RHONU*(nu(-1) - NUBAR) + epsnu;
end;
% =========================================================================
% Steady state using a steady_state_model block
% =========================================================================
steady_state_model;
OMEGP = (1+THETAP)*(1+PSI)/(1+PSI+THETAP*PSI);
q = QBAR;
pie = PIEBAR;
GAMMAP = PIEBAR^IOTA*PIEBAR^(1-IOTA);
b = BBAR;
tau = TAUBAR;
taun = TAUNBAR;
nu = NUBAR;
r = pie/BETA;
RK = (1/BETA+DELTA-1)*q; % foc k
aux1 = ( (1-XIP)/(1-XIP*(GAMMAP/pie)^(1/(1-OMEGP)) ) )^(1-OMEGP);
aux2 = (1-XIP)/(1-XIP*(GAMMAP/pie));
ptilde = (1+PSI)/ ( aux1 + PSI*aux2 );
delta2 = aux2*ptilde;
delta3 = aux1*ptilde;
vartheta = delta3;
delta1 = ( (1-XIP)/(1-XIP*(GAMMAP/pie)^(OMEGP/(1-OMEGP) )) )^((1-OMEGP)/OMEGP)*ptilde;
pstar = 1/(1+PSI)*vartheta^(OMEGP/(OMEGP-1))*delta1^(OMEGP/(1-OMEGP)) + PSI/(1+PSI);
f_o_a = ( vartheta^(OMEGP/(OMEGP-1)) / (1-BETA*XIP*(GAMMAP/pie)^(1/(1-OMEGP))) ) / ( PSI*(OMEGP-1) /(1-BETA*XIP*GAMMAP/pie) );
s_o_a = ( OMEGP*vartheta^(OMEGP/(OMEGP-1)) / (1-BETA*XIP*(GAMMAP/pie)^(OMEGP/(1-OMEGP))) ) / ( PSI*(OMEGP-1) /(1-BETA*XIP*GAMMAP/pie) );
mc = (f_o_a*ptilde-ptilde^(1+OMEGP/(OMEGP-1)))*s_o_a^(-1);
k_o_n = (RK/(mc*ALPHA))^(1/(ALPHA-1));
w = (1-ALPHA)*mc*k_o_n^ALPHA; % labor demand
y_o_n = pstar^(-1)*k_o_n^ALPHA; % production function
iv_o_n = DELTA*k_o_n;
iv_o_y = iv_o_n / y_o_n;
g_o_y = GBAR_o_YBAR;
c_o_n = y_o_n - iv_o_n - g_o_y; % market clearing
% The labor level, closed-form solution for labor
n = (c_o_n^(-1)*w)^(1/(CHI+1)) ;
% Value of remaining variables
c = c_o_n*n;
rk = RK;
iv = iv_o_n *n;
k = k_o_n*n;
y = y_o_n*n;
lam = n^CHI/w;
g = g_o_y*y;
a = PSI*(OMEGP-1)*y*lam/(1-BETA*XIP*GAMMAP/pie);
f = f_o_a*a;
s = f*ptilde - a*ptilde^(1+OMEGP/(OMEGP-1));
%flex price economy
bpot = BBAR;
taupot = TAUBAR;
taun = TAUNBAR;
rrpot = 1/BETA;
RKpot = (1/BETA+DELTA-1)*QBAR; % foc k
rkpot = RKpot;
kpot_o_npot = (RKpot/(ALPHA/(1+THETAP)))^(1/(ALPHA-1));
wpot = (1-ALPHA)/(1+THETAP)*kpot_o_npot^ALPHA; % labor demand
ypot_o_npot = kpot_o_npot^ALPHA; % production function
ivpot_o_npot = DELTA*kpot_o_npot;
ivpot_o_ypot = ivpot_o_npot / ypot_o_npot;
cpot_o_npot = ypot_o_npot - ivpot_o_npot - g_o_y; % market clearing
% The labor level, closed-form solution for labor
npot = (cpot_o_npot^(-1)*wpot)^(1/(CHI+1)) ;
% Value of remaining variables
cpot = cpot_o_npot*npot;
ivpot = ivpot_o_npot *npot;
kpot = kpot_o_npot*npot;
ypot = ypot_o_npot*npot;
end;
% =========================================================================
% Declare settings for shocks
% =========================================================================
shocks;
var epsg = 1;
var epsnu = 1;
end;
estimated_params;
ALPHA, 0.3, normal_pdf, 0.3 , 0.1;
BETA, 0.995, normal_pdf, 0.995, 0.1;
PIEBAR, 1.005, normal_pdf, 1.005, 0.1;
IOTA, 0.5, normal_pdf, 0.5, 0.1;
CHI, 1/0.4, normal_pdf, 1/0.4, 0.1;
THETAP, 0.1, normal_pdf, 0.1, 0.1;
XIP, 2/3, normal_pdf, 2/3, 0.1;
PSI, -12.2, normal_pdf, -12.2, 0.1;
GAMMAPIE, 1.5, normal_pdf, 1.5, 0.1;
GAMMAX, 0.125, normal_pdf, 0.125, 0.1;
DELTA, 0, normal_pdf, 0, 0.1;
QBAR, 1, normal_pdf, 1, 0.1;
NUBAR, 0, normal_pdf, 0, 0.1;
BBAR, 0, normal_pdf, 0, 0.1;
GBAR_o_YBAR, 0, normal_pdf, 0, 0.1; %commenting this solves the analytic_derivation_mode=-2|-1 problem
TAUBAR, 0, normal_pdf, 0, 0.1;
TAUNBAR, 0, normal_pdf, 0, 0.1;
PHI, 0.0101, normal_pdf, 0.0101, 0.1;
RHONU, 0.95, normal_pdf, 0.95, 0.1;
RHOG, 0.95, normal_pdf, 0.95, 0.1;
SIGNU, 0.6/100, normal_pdf, 0.6/100, 0.1;
SIGG, 0.6/100, normal_pdf, 0.6/100, 0.1;
end;
% =========================================================================
% Computations
% =========================================================================
model_diagnostics;
steady;
resid;
check;
identification(order=1,no_identification_strength,analytic_derivation_mode= 0,ar=5); %works
%identification(no_identification_strength,analytic_derivation_mode= 1,ar=5); %works, this takes the longest due to Kronecker products
options_.dynatol.x = 1e-9;
identification(order=1,no_identification_strength,analytic_derivation_mode=-1,ar=5); %works, but tolerance is way too small
identification(order=1,no_identification_strength,analytic_derivation_mode=-2,ar=5); %works, but tolerance is way to small
options_.dynatol.x = 1e-5; %this is the default value
identification(order=1,no_identification_strength,analytic_derivation_mode=-1,ar=5); %works only if GBAR_o_YBAR is uncommented
identification(order=1,no_identification_strength,analytic_derivation_mode=-2,ar=5); %works only if GBAR_o_YBAR is uncommented

BIN
tests/identification/as2007/G_QT.mat Normal file
View File

Binary file not shown.

41
tests/identification/as2007/as2007_QT.mod
View File

@ -74,24 +74,23 @@ identification(parameter_set=calibration,
checks_via_subsets=1,
max_dim_subsets_groups=4);
% load('G_QT'); %note that this is computed using replication files of Qu and Tkachenko (2012)
% temp = load([M_.dname filesep 'identification' filesep M_.fname '_identif']);
% G_dynare = temp.ide_spectrum_point.G;
%
% % Compare signs
% if ~isequal(sign(G_dynare),sign(G_QT))
% error('signs of normalized G are note equal');
% end
%
% % Compare normalized versions
% tilda_G_dynare = temp.ide_spectrum_point.tilda_G;
% ind_G_QT = (find(max(abs(G_QT'),[],1) > temp.store_options_ident.tol_deriv));
% tilda_G_QT = zeros(size(G_QT));
% delta_G_QT = sqrt(diag(G_QT(ind_G_QT,ind_G_QT)));
% tilda_G_QT(ind_G_QT,ind_G_QT) = G_QT(ind_G_QT,ind_G_QT)./((delta_G_QT)*(delta_G_QT'));
% if ~isequal(rank(tilda_G_QT,temp.store_options_ident.tol_rank),rank(tilda_G_dynare,temp.store_options_ident.tol_rank))
% error('ranks are not the same for normalized version')
% end
%
% max(max(abs(abs(tilda_G_dynare)-abs(tilda_G_QT))))
% norm(tilda_G_dynare - tilda_G_QT)
load('G_QT'); %note that this is computed using replication files of Qu and Tkachenko (2012)
temp = load([M_.dname filesep 'identification' filesep M_.fname '_identif']);
G_dynare = temp.ide_spectrum_point.G;
% Compare signs
if ~isequal(sign(G_dynare),sign(G_QT))
error('signs of normalized G are note equal');
end
% Compare normalized versions
tilda_G_dynare = temp.ide_spectrum_point.tilda_G;
ind_G_QT = (find(max(abs(G_QT'),[],1) > temp.store_options_ident.tol_deriv));
tilda_G_QT = zeros(size(G_QT));
delta_G_QT = sqrt(diag(G_QT(ind_G_QT,ind_G_QT)));
tilda_G_QT(ind_G_QT,ind_G_QT) = G_QT(ind_G_QT,ind_G_QT)./((delta_G_QT)*(delta_G_QT'));
if ~isequal(rank(tilda_G_QT,temp.store_options_ident.tol_rank),rank(tilda_G_dynare,temp.store_options_ident.tol_rank))
error('ranks are not the same for normalized version')
end
max(max(abs(abs(tilda_G_dynare)-abs(tilda_G_QT))))
norm(tilda_G_dynare - tilda_G_QT)