From acb9518deca048e796efa45969a8c371227acd7a Mon Sep 17 00:00:00 2001 From: Johannes Pfeifer Date: Mon, 25 Sep 2023 22:11:21 +0200 Subject: [PATCH] get_perturbation_params_derivs and friends: do not pass full oo_ structure --- matlab/+mom/objective_function.m | 2 +- matlab/+osr/objective.m | 2 +- matlab/dsge_likelihood.m | 4 +- matlab/fjaco.m | 8 +- matlab/get_identification_jacobians.m | 123 +++++++++--------- matlab/get_perturbation_params_derivs.m | 76 ++++++----- ...bation_params_derivs_numerical_objective.m | 27 ++-- matlab/identification_analysis.m | 4 +- matlab/identification_numerical_objective.m | 32 ++--- matlab/set_dynare_seed.m | 7 +- .../BrockMirman_PertParamsDerivs.mod | 2 +- .../burnside_3_order_PertParamsDerivs.mod | 2 +- 12 files changed, 151 insertions(+), 138 deletions(-) diff --git a/matlab/+mom/objective_function.m b/matlab/+mom/objective_function.m index b9fe38ba6..446c112b3 100644 --- a/matlab/+mom/objective_function.m +++ b/matlab/+mom/objective_function.m @@ -151,7 +151,7 @@ if strcmp(options_mom_.mom.mom_method,'GMM') stderrparam_nbr = estim_params_.nvx; % number of stderr parameters corrparam_nbr = estim_params_.ncx; % number of corr parameters totparam_nbr = stderrparam_nbr+corrparam_nbr+modparam_nbr; - oo_.dr.derivs = get_perturbation_params_derivs(M_, options_mom_, estim_params_, oo_, indpmodel, indpstderr, indpcorr, 0); %analytic derivatives of perturbation matrices + oo_.dr.derivs = get_perturbation_params_derivs(M_, options_mom_, estim_params_, oo_.dr, oo_.steady_state, oo_.exo_steady_state, oo_.exo_det_steady_state, indpmodel, indpstderr, indpcorr, 0); %analytic derivatives of perturbation matrices oo_.mom.model_moments_params_derivs = NaN(options_mom_.mom.mom_nbr,totparam_nbr); pruned_state_space = pruned_state_space_system(M_, options_mom_, oo_.dr, oo_.mom.obs_var, options_mom_.ar, 0, 1); else diff --git a/matlab/+osr/objective.m b/matlab/+osr/objective.m index 590591241..1bb3148ad 100644 --- a/matlab/+osr/objective.m +++ b/matlab/+osr/objective.m @@ -64,7 +64,7 @@ if ~options_.analytic_derivation loss = full(weights(:)'*vx(:)); else totparam_nbr=length(i_params); - oo_.dr.derivs = get_perturbation_params_derivs(M_, options_, [], oo_, i_params, [], [], 0); %analytic derivatives of perturbation matrices + oo_.dr.derivs = get_perturbation_params_derivs(M_, options_, [], oo_.dr, oo_.steady_state, oo_.exo_steady_state, oo_.exo_det_steady_state, i_params, [], [], 0); %analytic derivatives of perturbation matrices pruned_state_space = pruned_state_space_system(M_, options_, oo_.dr, i_var, 0, 0, 1); vx = pruned_state_space.Var_y + pruned_state_space.E_y*pruned_state_space.E_y'; diff --git a/matlab/dsge_likelihood.m b/matlab/dsge_likelihood.m index 3110a02f9..f7b41c42f 100644 --- a/matlab/dsge_likelihood.m +++ b/matlab/dsge_likelihood.m @@ -560,14 +560,14 @@ if analytic_derivation 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); + DERIVS = get_perturbation_params_derivs(Model, DynareOptions, EstimatedParameters, DynareResults.dr, DynareResults.steady_state, DynareResults.exo_steady_state, DynareResults.exo_det_steady_state, 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); + DERIVS = get_perturbation_params_derivs(Model, DynareOptions, EstimatedParameters, DynareResults.dr, DynareResults.steady_state, DynareResults.exo_steady_state, DynareResults.exo_det_steady_state, indparam, indexo, [], false); end DT = zeros(Model.endo_nbr, Model.endo_nbr, size(DERIVS.dghx,3)); DT(:,Model.nstatic+(1:Model.nspred),:) = DERIVS.dghx; diff --git a/matlab/fjaco.m b/matlab/fjaco.m index 25e517912..4f582d3db 100644 --- a/matlab/fjaco.m +++ b/matlab/fjaco.m @@ -10,7 +10,7 @@ function fjac = fjaco(f,x,varargin) % OUTPUT % fjac : finite difference Jacobian % -% Copyright © 2010-2017,2019-2020 Dynare Team +% Copyright © 2010-2017,2019-2023 Dynare Team % % This file is part of Dynare. % @@ -31,7 +31,7 @@ ff=feval(f,x,varargin{:}); 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; + tol= varargin{4}.dynatol.x; end h = tol.*max(abs(x),1); xh1=x+h; xh0=x-h; @@ -57,9 +57,9 @@ 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}); + probl_par = get_the_name(j,varargin{4}.TeX,varargin{3},varargin{2},varargin{4}); skipline() - message = get_error_message(info,varargin{5}); + message = get_error_message(info,varargin{4}); 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') diff --git a/matlab/get_identification_jacobians.m b/matlab/get_identification_jacobians.m index 7b8d2bd81..d04d49266 100644 --- a/matlab/get_identification_jacobians.m +++ b/matlab/get_identification_jacobians.m @@ -1,12 +1,11 @@ -function [MEAN, dMEAN, REDUCEDFORM, dREDUCEDFORM, DYNAMIC, dDYNAMIC, MOMENTS, dMOMENTS, dSPECTRUM, dSPECTRUM_NO_MEAN, dMINIMAL, derivatives_info] = get_identification_jacobians(estim_params, M, oo, options, options_ident, indpmodel, indpstderr, indpcorr, indvobs) -% function [MEAN, dMEAN, REDUCEDFORM, dREDUCEDFORM, DYNAMIC, dDYNAMIC, MOMENTS, dMOMENTS, dSPECTRUM, dMINIMAL, derivatives_info] = get_identification_jacobians(estim_params, M, oo, options, options_ident, indpmodel, indpstderr, indpcorr, indvobs) +function [MEAN, dMEAN, REDUCEDFORM, dREDUCEDFORM, DYNAMIC, dDYNAMIC, MOMENTS, dMOMENTS, dSPECTRUM, dSPECTRUM_NO_MEAN, dMINIMAL, derivatives_info] = get_identification_jacobians(estim_params, M, options, options_ident, indpmodel, indpstderr, indpcorr, indvobs, dr, endo_steady_state, exo_steady_state, exo_det_steady_state) +% [MEAN, dMEAN, REDUCEDFORM, dREDUCEDFORM, DYNAMIC, dDYNAMIC, MOMENTS, dMOMENTS, dSPECTRUM, dSPECTRUM_NO_MEAN, dMINIMAL, derivatives_info] = get_identification_jacobians(estim_params, M, options, options_ident, indpmodel, indpstderr, indpcorr, indvobs, dr, endo_steady_state, exo_steady_state, exo_det_steady_state) % previously getJJ.m in Dynare 4.5 % Sets up the Jacobians needed for identification analysis % ========================================================================= % INPUTS % 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 % options_ident: [structure] storing the options for identification % indpmodel: [modparam_nbr by 1] index of estimated parameters in M_.params; @@ -22,6 +21,10 @@ function [MEAN, dMEAN, REDUCEDFORM, dREDUCEDFORM, DYNAMIC, dDYNAMIC, MOMENTS, dM % in the estimated_params block; if estimated_params block % is not available, then no corr parameters are selected % indvobs: [obs_nbr by 1] index of observed (VAROBS) variables +% dr [structure] Reduced form model. +% endo_steady_state [vector] steady state value for endogenous variables +% exo_steady_state [vector] steady state value for exogenous variables +% exo_det_steady_state [vector] steady state value for exogenous deterministic variables % ------------------------------------------------------------------------- % OUTPUTS % @@ -150,94 +153,94 @@ obs_nbr = length(indvobs); d2flag = 0; % do not compute second parameter derivatives % Get Jacobians (wrt selected params) of steady state, dynamic model derivatives and perturbation solution matrices for all endogenous variables -oo.dr.derivs = get_perturbation_params_derivs(M, options, estim_params, oo, indpmodel, indpstderr, indpcorr, d2flag); +dr.derivs = get_perturbation_params_derivs(M, options, estim_params, dr, endo_steady_state, exo_steady_state, exo_det_steady_state, indpmodel, indpstderr, indpcorr, d2flag); [I,~] = find(lead_lag_incidence'); %I is used to select nonzero columns of the Jacobian of endogenous variables in dynamic model files -yy0 = oo.dr.ys(I); %steady state of dynamic (endogenous and auxiliary variables) in lead_lag_incidence order -Yss = oo.dr.ys(oo.dr.order_var); % steady state in DR order +yy0 = dr.ys(I); %steady state of dynamic (endogenous and auxiliary variables) in lead_lag_incidence order +Yss = dr.ys(dr.order_var); % steady state in DR order if order == 1 - [~, g1 ] = feval([fname,'.dynamic'], yy0, oo.exo_steady_state', params, oo.dr.ys, 1); + [~, g1 ] = feval([fname,'.dynamic'], yy0, exo_steady_state', params, dr.ys, 1); %g1 is [endo_nbr by yy0ex0_nbr first derivative (wrt all dynamic variables) of dynamic model equations, i.e. df/dyy0ex0, rows are in declaration order, columns in lead_lag_incidence order DYNAMIC = [Yss; vec(g1)]; %add steady state and put rows of g1 in DR order - dDYNAMIC = [oo.dr.derivs.dYss; - reshape(oo.dr.derivs.dg1,size(oo.dr.derivs.dg1,1)*size(oo.dr.derivs.dg1,2),size(oo.dr.derivs.dg1,3)) ]; %reshape dg1 in DR order and add steady state + dDYNAMIC = [dr.derivs.dYss; + reshape(dr.derivs.dg1,size(dr.derivs.dg1,1)*size(dr.derivs.dg1,2),size(dr.derivs.dg1,3)) ]; %reshape dg1 in DR order and add steady state REDUCEDFORM = [Yss; - vec(oo.dr.ghx); - dyn_vech(oo.dr.ghu*Sigma_e*transpose(oo.dr.ghu))]; %in DR order + vec(dr.ghx); + dyn_vech(dr.ghu*Sigma_e*transpose(dr.ghu))]; %in DR order dREDUCEDFORM = zeros(endo_nbr*nspred+endo_nbr*(endo_nbr+1)/2, totparam_nbr); for j=1:totparam_nbr - dREDUCEDFORM(:,j) = [vec(oo.dr.derivs.dghx(:,:,j)); - dyn_vech(oo.dr.derivs.dOm(:,:,j))]; + dREDUCEDFORM(:,j) = [vec(dr.derivs.dghx(:,:,j)); + dyn_vech(dr.derivs.dOm(:,:,j))]; end - dREDUCEDFORM = [ [zeros(endo_nbr, stderrparam_nbr+corrparam_nbr) oo.dr.derivs.dYss]; dREDUCEDFORM ]; % add steady state + dREDUCEDFORM = [ [zeros(endo_nbr, stderrparam_nbr+corrparam_nbr) dr.derivs.dYss]; dREDUCEDFORM ]; % add steady state elseif order == 2 - [~, g1, g2 ] = feval([fname,'.dynamic'], yy0, oo.exo_steady_state', params, oo.dr.ys, 1); + [~, g1, g2 ] = feval([fname,'.dynamic'], yy0, exo_steady_state', params, dr.ys, 1); %g1 is [endo_nbr by yy0ex0_nbr first derivative (wrt all dynamic variables) of dynamic model equations, i.e. df/dyy0ex0, rows are in declaration order, columns in lead_lag_incidence order %g2 is [endo_nbr by yy0ex0_nbr^2] second derivative (wrt all dynamic variables) of dynamic model equations, i.e. d(df/dyy0ex0)/dyy0ex0, rows are in declaration order, columns in lead_lag_incidence order DYNAMIC = [Yss; vec(g1); vec(g2)]; %add steady state and put rows of g1 and g2 in DR order - dDYNAMIC = [oo.dr.derivs.dYss; - reshape(oo.dr.derivs.dg1,size(oo.dr.derivs.dg1,1)*size(oo.dr.derivs.dg1,2),size(oo.dr.derivs.dg1,3)); %reshape dg1 in DR order - reshape(oo.dr.derivs.dg2,size(oo.dr.derivs.dg1,1)*size(oo.dr.derivs.dg1,2)^2,size(oo.dr.derivs.dg1,3))]; %reshape dg2 in DR order + dDYNAMIC = [dr.derivs.dYss; + reshape(dr.derivs.dg1,size(dr.derivs.dg1,1)*size(dr.derivs.dg1,2),size(dr.derivs.dg1,3)); %reshape dg1 in DR order + reshape(dr.derivs.dg2,size(dr.derivs.dg1,1)*size(dr.derivs.dg1,2)^2,size(dr.derivs.dg1,3))]; %reshape dg2 in DR order REDUCEDFORM = [Yss; - vec(oo.dr.ghx); - dyn_vech(oo.dr.ghu*Sigma_e*transpose(oo.dr.ghu)); - vec(oo.dr.ghxx); - vec(oo.dr.ghxu); - vec(oo.dr.ghuu); - vec(oo.dr.ghs2)]; %in DR order + vec(dr.ghx); + dyn_vech(dr.ghu*Sigma_e*transpose(dr.ghu)); + vec(dr.ghxx); + vec(dr.ghxu); + vec(dr.ghuu); + vec(dr.ghs2)]; %in DR order dREDUCEDFORM = zeros(endo_nbr*nspred+endo_nbr*(endo_nbr+1)/2+endo_nbr*nspred^2+endo_nbr*nspred*exo_nbr+endo_nbr*exo_nbr^2+endo_nbr, totparam_nbr); for j=1:totparam_nbr - dREDUCEDFORM(:,j) = [vec(oo.dr.derivs.dghx(:,:,j)); - dyn_vech(oo.dr.derivs.dOm(:,:,j)); - vec(oo.dr.derivs.dghxx(:,:,j)); - vec(oo.dr.derivs.dghxu(:,:,j)); - vec(oo.dr.derivs.dghuu(:,:,j)); - vec(oo.dr.derivs.dghs2(:,j))]; + dREDUCEDFORM(:,j) = [vec(dr.derivs.dghx(:,:,j)); + dyn_vech(dr.derivs.dOm(:,:,j)); + vec(dr.derivs.dghxx(:,:,j)); + vec(dr.derivs.dghxu(:,:,j)); + vec(dr.derivs.dghuu(:,:,j)); + vec(dr.derivs.dghs2(:,j))]; end - dREDUCEDFORM = [ [zeros(endo_nbr, stderrparam_nbr+corrparam_nbr) oo.dr.derivs.dYss]; dREDUCEDFORM ]; % add steady state + dREDUCEDFORM = [ [zeros(endo_nbr, stderrparam_nbr+corrparam_nbr) dr.derivs.dYss]; dREDUCEDFORM ]; % add steady state elseif order == 3 - [~, g1, g2, g3 ] = feval([fname,'.dynamic'], yy0, oo.exo_steady_state', params, oo.dr.ys, 1); + [~, g1, g2, g3 ] = feval([fname,'.dynamic'], yy0, exo_steady_state', params, dr.ys, 1); %g1 is [endo_nbr by yy0ex0_nbr first derivative (wrt all dynamic variables) of dynamic model equations, i.e. df/dyy0ex0, rows are in declaration order, columns in lead_lag_incidence order %g2 is [endo_nbr by yy0ex0_nbr^2] second derivative (wrt all dynamic variables) of dynamic model equations, i.e. d(df/dyy0ex0)/dyy0ex0, rows are in declaration order, columns in lead_lag_incidence order DYNAMIC = [Yss; vec(g1); vec(g2); vec(g3)]; %add steady state and put rows of g1 and g2 in DR order - dDYNAMIC = [oo.dr.derivs.dYss; - reshape(oo.dr.derivs.dg1,size(oo.dr.derivs.dg1,1)*size(oo.dr.derivs.dg1,2),size(oo.dr.derivs.dg1,3)); %reshape dg1 in DR order - reshape(oo.dr.derivs.dg2,size(oo.dr.derivs.dg1,1)*size(oo.dr.derivs.dg1,2)^2,size(oo.dr.derivs.dg1,3)); - reshape(oo.dr.derivs.dg2,size(oo.dr.derivs.dg1,1)*size(oo.dr.derivs.dg1,2)^2,size(oo.dr.derivs.dg1,3))]; %reshape dg3 in DR order + dDYNAMIC = [dr.derivs.dYss; + reshape(dr.derivs.dg1,size(dr.derivs.dg1,1)*size(dr.derivs.dg1,2),size(dr.derivs.dg1,3)); %reshape dg1 in DR order + reshape(dr.derivs.dg2,size(dr.derivs.dg1,1)*size(dr.derivs.dg1,2)^2,size(dr.derivs.dg1,3)); + reshape(dr.derivs.dg2,size(dr.derivs.dg1,1)*size(dr.derivs.dg1,2)^2,size(dr.derivs.dg1,3))]; %reshape dg3 in DR order REDUCEDFORM = [Yss; - vec(oo.dr.ghx); - dyn_vech(oo.dr.ghu*Sigma_e*transpose(oo.dr.ghu)); - vec(oo.dr.ghxx); vec(oo.dr.ghxu); vec(oo.dr.ghuu); vec(oo.dr.ghs2); - vec(oo.dr.ghxxx); vec(oo.dr.ghxxu); vec(oo.dr.ghxuu); vec(oo.dr.ghuuu); vec(oo.dr.ghxss); vec(oo.dr.ghuss)]; %in DR order + vec(dr.ghx); + dyn_vech(dr.ghu*Sigma_e*transpose(dr.ghu)); + vec(dr.ghxx); vec(dr.ghxu); vec(dr.ghuu); vec(dr.ghs2); + vec(dr.ghxxx); vec(dr.ghxxu); vec(dr.ghxuu); vec(dr.ghuuu); vec(dr.ghxss); vec(dr.ghuss)]; %in DR order dREDUCEDFORM = zeros(size(REDUCEDFORM,1)-endo_nbr, totparam_nbr); for j=1:totparam_nbr - dREDUCEDFORM(:,j) = [vec(oo.dr.derivs.dghx(:,:,j)); - dyn_vech(oo.dr.derivs.dOm(:,:,j)); - vec(oo.dr.derivs.dghxx(:,:,j)); vec(oo.dr.derivs.dghxu(:,:,j)); vec(oo.dr.derivs.dghuu(:,:,j)); vec(oo.dr.derivs.dghs2(:,j)) - vec(oo.dr.derivs.dghxxx(:,:,j)); vec(oo.dr.derivs.dghxxu(:,:,j)); vec(oo.dr.derivs.dghxuu(:,:,j)); vec(oo.dr.derivs.dghuuu(:,:,j)); vec(oo.dr.derivs.dghxss(:,:,j)); vec(oo.dr.derivs.dghuss(:,:,j))]; + dREDUCEDFORM(:,j) = [vec(dr.derivs.dghx(:,:,j)); + dyn_vech(dr.derivs.dOm(:,:,j)); + vec(dr.derivs.dghxx(:,:,j)); vec(dr.derivs.dghxu(:,:,j)); vec(dr.derivs.dghuu(:,:,j)); vec(dr.derivs.dghs2(:,j)) + vec(dr.derivs.dghxxx(:,:,j)); vec(dr.derivs.dghxxu(:,:,j)); vec(dr.derivs.dghxuu(:,:,j)); vec(dr.derivs.dghuuu(:,:,j)); vec(dr.derivs.dghxss(:,:,j)); vec(dr.derivs.dghuss(:,:,j))]; end - dREDUCEDFORM = [ [zeros(endo_nbr, stderrparam_nbr+corrparam_nbr) oo.dr.derivs.dYss]; dREDUCEDFORM ]; % add steady state + dREDUCEDFORM = [ [zeros(endo_nbr, stderrparam_nbr+corrparam_nbr) dr.derivs.dYss]; dREDUCEDFORM ]; % add steady state end % Get (pruned) state space representation: -pruned = pruned_state_space_system(M, options, oo.dr, indvobs, nlags, useautocorr, 1); +pruned = pruned_state_space_system(M, options, dr, indvobs, nlags, useautocorr, 1); MEAN = pruned.E_y; dMEAN = pruned.dE_y; %storage for Jacobians used in dsge_likelihood.m for analytical Gradient and Hession of likelihood (only at order=1) derivatives_info = struct(); if order == 1 dT = zeros(endo_nbr,endo_nbr,totparam_nbr); - dT(:,(nstatic+1):(nstatic+nspred),:) = oo.dr.derivs.dghx; + dT(:,(nstatic+1):(nstatic+nspred),:) = dr.derivs.dghx; derivatives_info.DT = dT; - derivatives_info.DOm = oo.dr.derivs.dOm; - derivatives_info.DYss = oo.dr.derivs.dYss; + derivatives_info.DOm = dr.derivs.dOm; + derivatives_info.DYss = dr.derivs.dYss; end %% Compute dMOMENTS @@ -255,7 +258,7 @@ if ~no_identification_moments if kronflag == -1 %numerical derivative of autocovariogram - dMOMENTS = fjaco(str2func('identification_numerical_objective'), xparam1, 1, estim_params, M, oo, options, indpmodel, indpstderr, indpcorr, indvobs, useautocorr, nlags, grid_nbr); %[outputflag=1] + dMOMENTS = fjaco(str2func('identification_numerical_objective'), xparam1, 1, estim_params, M, options, indpmodel, indpstderr, indvobs, useautocorr, nlags, grid_nbr, dr, endo_steady_state, exo_steady_state, exo_det_steady_state); %[outputflag=1] dMOMENTS = [dMEAN; dMOMENTS]; %add Jacobian of steady state of VAROBS variables else dMOMENTS = zeros(obs_nbr + obs_nbr*(obs_nbr+1)/2 + nlags*obs_nbr^2 , totparam_nbr); @@ -312,7 +315,7 @@ if ~no_identification_spectrum IA = eye(size(pruned.A,1)); if kronflag == -1 %numerical derivative of spectral density - dOmega_tmp = fjaco(str2func('identification_numerical_objective'), xparam1, 2, estim_params, M, oo, options, indpmodel, indpstderr, indpcorr, indvobs, useautocorr, nlags, grid_nbr); %[outputflag=2] + dOmega_tmp = fjaco(str2func('identification_numerical_objective'), xparam1, 2, estim_params, M, options, indpmodel, indpstderr, indvobs, useautocorr, nlags, grid_nbr, dr, endo_steady_state, exo_steady_state, exo_det_steady_state); %[outputflag=2] kk = 0; for ig = 1:length(freqs) kk = kk+1; @@ -389,14 +392,14 @@ if ~no_identification_minimal dMINIMAL = []; else % Derive and check minimal state vector of first-order - SYS.A = oo.dr.ghx(pruned.indx,:); - SYS.dA = oo.dr.derivs.dghx(pruned.indx,:,:); - SYS.B = oo.dr.ghu(pruned.indx,:); - SYS.dB = oo.dr.derivs.dghu(pruned.indx,:,:); - SYS.C = oo.dr.ghx(pruned.indy,:); - SYS.dC = oo.dr.derivs.dghx(pruned.indy,:,:); - SYS.D = oo.dr.ghu(pruned.indy,:); - SYS.dD = oo.dr.derivs.dghu(pruned.indy,:,:); + SYS.A = dr.ghx(pruned.indx,:); + SYS.dA = dr.derivs.dghx(pruned.indx,:,:); + SYS.B = dr.ghu(pruned.indx,:); + SYS.dB = dr.derivs.dghu(pruned.indx,:,:); + SYS.C = dr.ghx(pruned.indy,:); + SYS.dC = dr.derivs.dghx(pruned.indy,:,:); + SYS.D = dr.ghu(pruned.indy,:); + SYS.dD = dr.derivs.dghu(pruned.indy,:,:); [CheckCO,minnx,SYS] = get_minimal_state_representation(SYS,1); if CheckCO == 0 @@ -415,7 +418,7 @@ if ~no_identification_minimal dminB = reshape(dminB,size(dminB,1)*size(dminB,2),size(dminB,3)); dminC = reshape(dminC,size(dminC,1)*size(dminC,2),size(dminC,3)); dminD = reshape(dminD,size(dminD,1)*size(dminD,2),size(dminD,3)); - dvechSig = reshape(oo.dr.derivs.dSigma_e,exo_nbr*exo_nbr,totparam_nbr); + dvechSig = reshape(dr.derivs.dSigma_e,exo_nbr*exo_nbr,totparam_nbr); indvechSig= find(tril(ones(exo_nbr,exo_nbr))); dvechSig = dvechSig(indvechSig,:); Inx = eye(minnx); diff --git a/matlab/get_perturbation_params_derivs.m b/matlab/get_perturbation_params_derivs.m index 6475ac204..70540b1a3 100644 --- a/matlab/get_perturbation_params_derivs.m +++ b/matlab/get_perturbation_params_derivs.m @@ -1,5 +1,5 @@ -function DERIVS = get_perturbation_params_derivs(M, options, estim_params, oo, indpmodel, indpstderr, indpcorr, d2flag) -% DERIVS = get_perturbation_params_derivs(M, options, estim_params, oo, indpmodel, indpstderr, indpcorr, d2flag) +function DERIVS = get_perturbation_params_derivs(M, options, estim_params, dr, endo_steady_state, exo_steady_state, exo_det_steady_state, indpmodel, indpstderr, indpcorr, d2flag) +% DERIVS = get_perturbation_params_derivs(M, options, estim_params, dr, endo_steady_state, exo_steady_state, exo_det_steady_state, indpmodel, indpstderr, indpcorr, d2flag) % previously getH.m in dynare 4.5 % ------------------------------------------------------------------------- % Computes derivatives (with respect to parameters) of @@ -16,17 +16,24 @@ function DERIVS = get_perturbation_params_derivs(M, options, estim_params, oo, i % % ========================================================================= % INPUTS -% M: [structure] storing the model information -% options: [structure] storing the options -% estim_params: [structure] storing the estimation information -% oo: [structure] storing the solution results -% indpmodel: [modparam_nbr by 1] index of selected (estimated) parameters in M.params; -% corresponds to model parameters (no stderr and no corr) in estimated_params block -% indpstderr: [stderrparam_nbr by 1] index of selected (estimated) standard errors, -% i.e. for all exogenous variables where 'stderr' is given in the estimated_params block -% indpcorr: [corrparam_nbr by 2] matrix of selected (estimated) correlations, -% i.e. for all exogenous variables where 'corr' is given in the estimated_params block -% d2flag: [boolean] flag to compute second-order parameter derivatives of steady state and first-order Kalman transition matrices +% M: [structure] storing the model information +% options: [structure] storing the options +% estim_params: [structure] storing the estimation information +% dr [structure] Reduced form model. +% endo_steady_state [vector] steady state value for endogenous variables +% exo_steady_state [vector] steady state value for exogenous variables +% exo_det_steady_state [vector] steady state value for exogenous deterministic variables +% indpmodel: [modparam_nbr by 1] index of selected (estimated) parameters in M.params; +% corresponds to model parameters (no stderr and no corr) +% in estimated_params block +% indpstderr: [stderrparam_nbr by 1] index of selected (estimated) standard errors, +% i.e. for all exogenous variables where 'stderr' is given +% in the estimated_params block +% indpcorr: [corrparam_nbr by 2] matrix of selected (estimated) correlations, +% i.e. for all exogenous variables where 'corr' is given in +% the estimated_params block +% d2flag: [boolean] flag to compute second-order parameter derivatives of steady state +% and first-order Kalman transition matrices % ------------------------------------------------------------------------- % OUTPUTS % DERIVS: Structure with the following fields: @@ -154,26 +161,25 @@ end qz_criterium = options.qz_criterium; threads_BC = options.threads.kronecker.sparse_hessian_times_B_kronecker_C; -% Get fields from oo -exo_steady_state = oo.exo_steady_state; -ghx = oo.dr.ghx; -ghu = oo.dr.ghu; +% Get fields from dr +ghx = dr.ghx; +ghu = dr.ghu; if order > 1 - ghxx = oo.dr.ghxx; - ghxu = oo.dr.ghxu; - ghuu = oo.dr.ghuu; - ghs2 = oo.dr.ghs2; + ghxx = dr.ghxx; + ghxu = dr.ghxu; + ghuu = dr.ghuu; + ghs2 = dr.ghs2; end if order > 2 - ghxxx = oo.dr.ghxxx; - ghxxu = oo.dr.ghxxu; - ghxuu = oo.dr.ghxuu; - ghuuu = oo.dr.ghuuu; - ghxss = oo.dr.ghxss; - ghuss = oo.dr.ghuss; + ghxxx = dr.ghxxx; + ghxxu = dr.ghxxu; + ghxuu = dr.ghxuu; + ghuuu = dr.ghuuu; + ghxss = dr.ghxss; + ghuss = dr.ghuss; end -order_var = oo.dr.order_var; -ys = oo.dr.ys; +order_var = dr.order_var; +ys = dr.ys; % Some checks if exo_det_nbr > 0 @@ -289,7 +295,7 @@ if analytic_derivation_mode == -1 % - perturbation solution matrices: dghx, dghu, dghxx, dghxu, dghuu, dghs2, dghxxx, dghxxu, dghxuu, dghuuu, dghxss, dghuss %Parameter Jacobian of covariance matrix and solution matrices (wrt selected stderr, corr and model paramters) - dSig_gh = fjaco(numerical_objective_fname, xparam1, 'perturbation_solution', estim_params, M, oo, options); + dSig_gh = fjaco(numerical_objective_fname, xparam1, 'perturbation_solution', estim_params, M, options, dr, endo_steady_state, exo_steady_state, exo_det_steady_state); ind_Sigma_e = (1:exo_nbr^2); ind_ghx = ind_Sigma_e(end) + (1:endo_nbr*nspred); ind_ghu = ind_ghx(end) + (1:endo_nbr*exo_nbr); @@ -342,7 +348,7 @@ if analytic_derivation_mode == -1 end %Parameter Jacobian of dynamic model derivatives (wrt selected model parameters only) - dYss_g = fjaco(numerical_objective_fname, modparam1, 'dynamic_model', estim_params_model, M, oo, options); + dYss_g = fjaco(numerical_objective_fname, modparam1, 'dynamic_model', estim_params_model, M, options, dr, endo_steady_state, exo_steady_state, exo_det_steady_state); ind_Yss = 1:endo_nbr; if options.discretionary_policy || options.ramsey_policy ind_g1 = ind_Yss(end) + (1:M.eq_nbr*yy0ex0_nbr); @@ -368,7 +374,7 @@ if analytic_derivation_mode == -1 % Hessian (wrt paramters) of steady state and first-order solution matrices ghx and Om % note that hessian_sparse.m (contrary to hessian.m) does not take symmetry into account, but focuses already on unique values options.order = 1; %make sure only first order - d2Yss_KalmanA_Om = hessian_sparse(numerical_objective_fname, xparam1, gstep, 'Kalman_Transition', estim_params, M, oo, options); + d2Yss_KalmanA_Om = hessian_sparse(numerical_objective_fname, xparam1, gstep, 'Kalman_Transition', estim_params, M, options, dr, endo_steady_state, exo_steady_state, exo_det_steady_state); options.order = order; %make sure to set back ind_KalmanA = ind_Yss(end) + (1:endo_nbr^2); DERIVS.d2KalmanA = d2Yss_KalmanA_Om(ind_KalmanA, indp2tottot2); %only unique elements @@ -399,7 +405,7 @@ if analytic_derivation_mode == -2 % computation of d2Yss and d2g1 % 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 options.order = 1; %d2flag requires only first order - d2Yss_g1 = hessian_sparse(numerical_objective_fname, modparam1, gstep, 'dynamic_model', estim_params_model, M, oo, options); % d2flag requires only first-order + d2Yss_g1 = hessian_sparse(numerical_objective_fname, modparam1, gstep, 'dynamic_model', estim_params_model, M, options, dr, endo_steady_state, exo_steady_state, exo_det_steady_state); % d2flag requires only first-order options.order = order; %make sure to set back the order d2Yss = reshape(full(d2Yss_g1(1:endo_nbr,:)), [endo_nbr modparam_nbr modparam_nbr]); %put into tensor notation for j=1:endo_nbr @@ -425,7 +431,7 @@ if analytic_derivation_mode == -2 end %Parameter Jacobian of dynamic model derivatives (wrt selected model parameters only) - dYss_g = fjaco(numerical_objective_fname, modparam1, 'dynamic_model', estim_params_model, M, oo, options); + dYss_g = fjaco(numerical_objective_fname, modparam1, 'dynamic_model', estim_params_model, M, options, dr, endo_steady_state, exo_steady_state, exo_det_steady_state); ind_Yss = 1:endo_nbr; ind_g1 = ind_Yss(end) + (1:endo_nbr*yy0ex0_nbr); dYss = dYss_g(ind_Yss,:); %in tensor notation, wrt selected model parameters only @@ -592,7 +598,7 @@ elseif (analytic_derivation_mode == 0 || analytic_derivation_mode == 1) [II, JJ] = ind2sub([yy0ex0_nbr yy0ex0_nbr], find(g2(j,:))); %g2 is [endo_nbr by yy0ex0_nbr^2] for i = 1:yy0ex0_nbr is = find(II==i); - is = is(find(JJ(is)<=yy0_nbr)); %focus only on oo.dr.ys(I) derivatives as exogenous variables are 0 in a stochastic context + is = is(find(JJ(is)<=yy0_nbr)); %focus only on dr.ys(I) derivatives as exogenous variables are 0 in a stochastic context if ~isempty(is) tmp_g2 = full(g2(j,find(g2(j,:)))); dg1(j,i,:) = tmp_g2(is)*dyy0(JJ(is),:); %put into tensor notation diff --git a/matlab/get_perturbation_params_derivs_numerical_objective.m b/matlab/get_perturbation_params_derivs_numerical_objective.m index d9a4bbcd3..a4540250e 100644 --- a/matlab/get_perturbation_params_derivs_numerical_objective.m +++ b/matlab/get_perturbation_params_derivs_numerical_objective.m @@ -1,5 +1,5 @@ -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) +function [out,info] = get_perturbation_params_derivs_numerical_objective(params, outputflag, estim_params, M, options, dr, endo_steady_state, exo_steady_state, exo_det_steady_state) +%function [out,info] = get_perturbation_params_derivs_numerical_objective(params, outputflag, estim_params, M, options, dr, steady_state, exo_steady_state, exo_det_steady_state) % ------------------------------------------------------------------------- % Objective function to compute numerically the Jacobians used for get_perturbation_params_derivs % ========================================================================= @@ -9,8 +9,11 @@ function [out,info] = get_perturbation_params_derivs_numerical_objective(params, % 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 +% dr [structure] Reduced form model. +% endo_steady_state [vector] steady state value for endogenous variables +% exo_steady_state [vector] steady state value for exogenous variables +% exo_det_steady_state [vector] steady state value for exogenous deterministic variables % ------------------------------------------------------------------------- % % OUTPUT @@ -51,7 +54,7 @@ function [out,info] = get_perturbation_params_derivs_numerical_objective(params, %% Update stderr, corr and model parameters and compute perturbation approximation and steady state with updated parameters M = set_all_parameters(params,estim_params,M); -[oo.dr,info,M.params] = compute_decision_rules(M,options,oo.dr, oo.steady_state, oo.exo_steady_state, oo.exo_det_steady_state); +[dr,info,M.params] = compute_decision_rules(M,options,dr,endo_steady_state,exo_steady_state,exo_det_steady_state); Sigma_e = M.Sigma_e; if info(1) > 0 @@ -59,16 +62,16 @@ if info(1) > 0 out = []; return else - ys = oo.dr.ys; %steady state of model variables in declaration order - ghx = oo.dr.ghx; ghu = oo.dr.ghu; + ys = dr.ys; %steady state of model variables in declaration order + ghx = dr.ghx; ghu = dr.ghu; if options.order > 1 - ghxx = oo.dr.ghxx; ghxu = oo.dr.ghxu; ghuu = oo.dr.ghuu; ghs2 = oo.dr.ghs2; + ghxx = dr.ghxx; ghxu = dr.ghxu; ghuu = dr.ghuu; ghs2 = dr.ghs2; 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; + ghxxx = dr.ghxxx; ghxxu = dr.ghxxu; ghxuu = dr.ghxuu; ghxss = dr.ghxss; ghuuu = dr.ghuuu; ghuss = dr.ghuss; end end -Yss = ys(oo.dr.order_var); %steady state of model variables in DR order +Yss = ys(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') @@ -85,13 +88,13 @@ end 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); + [~, g1] = feval([M.fname,'.dynamic'], ys(I), 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); + [~, g1, g2] = feval([M.fname,'.dynamic'], ys(I), 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); + [~, g1, g2, g3] = feval([M.fname,'.dynamic'], ys(I), exo_steady_state', M.params, ys, 1); g3 = unfold_g3(g3, length(ys(I))+M.exo_nbr); out = [Yss; g1(:); g2(:); g3(:)]; end diff --git a/matlab/identification_analysis.m b/matlab/identification_analysis.m index 1feeb6552..b01235957 100644 --- a/matlab/identification_analysis.m +++ b/matlab/identification_analysis.m @@ -138,7 +138,7 @@ error_indicator.identification_spectrum=0; if info(1) == 0 %no errors in solution % Compute parameter Jacobians for identification analysis - [MEAN, dMEAN, REDUCEDFORM, dREDUCEDFORM, DYNAMIC, dDYNAMIC, MOMENTS, dMOMENTS, dSPECTRUM, dSPECTRUM_NO_MEAN, dMINIMAL, derivatives_info] = get_identification_jacobians(estim_params_, M_, oo_, options_, options_ident, indpmodel, indpstderr, indpcorr, indvobs); + [MEAN, dMEAN, REDUCEDFORM, dREDUCEDFORM, DYNAMIC, dDYNAMIC, MOMENTS, dMOMENTS, dSPECTRUM, dSPECTRUM_NO_MEAN, dMINIMAL, derivatives_info] = get_identification_jacobians(estim_params_, M_, options_, options_ident, indpmodel, indpstderr, indpcorr, indvobs, oo_.dr, oo_.steady_state, oo_.exo_steady_state, oo_.exo_det_steady_state); if isempty(dMINIMAL) % Komunjer and Ng is not computed if (1) minimality conditions are not fullfilled or (2) there are more shocks and measurement errors than observables, so we need to reset options error_indicator.identification_minimal = 1; @@ -202,7 +202,7 @@ if info(1) == 0 %no errors in solution options_ident_local.no_identification_spectrum = 1; %do not recompute dSPECTRUM options_ident_local.ar = nlags; %store new lag number options_.ar = nlags; %store new lag number - [~, ~, ~, ~, ~, ~, MOMENTS, dMOMENTS, ~, ~, ~, ~] = get_identification_jacobians(estim_params_, M_, oo_, options_, options_ident_local, indpmodel, indpstderr, indpcorr, indvobs); + [~, ~, ~, ~, ~, ~, MOMENTS, dMOMENTS, ~, ~, ~, ~] = get_identification_jacobians(estim_params_, M_, options_, options_ident_local, indpmodel, indpstderr, indpcorr, indvobs, oo_.dr, oo_.steady_state, oo_.exo_steady_state, oo_.exo_det_steady_state); ind_dMOMENTS = (find(max(abs(dMOMENTS'),[],1) > tol_deriv)); %new index with non-zero rows end diff --git a/matlab/identification_numerical_objective.m b/matlab/identification_numerical_objective.m index 172f1f6f2..4a0683fde 100644 --- a/matlab/identification_numerical_objective.m +++ b/matlab/identification_numerical_objective.m @@ -1,21 +1,23 @@ -function out = identification_numerical_objective(params, outputflag, estim_params, M, oo, options, indpmodel, indpstderr, indpcorr, indvar, useautocorr, nlags, grid_nbr) -%function out = identification_numerical_objective(params, outputflag, estim_params, M, oo, options, indpmodel, indpstderr, indpcorr, indvar, useautocorr, nlags, grid_nbr) +function out = identification_numerical_objective(params, outputflag, estim_params, M, options, indpmodel, indpstderr, indvar, useautocorr, nlags, grid_nbr, dr, steady_state, exo_steady_state, exo_det_steady_state) +%function out = identification_numerical_objective(params, outputflag, %estim_params, M, options, indpmodel, indpstderr, indvar, useautocorr, nlags, grid_nbr, dr, steady_state, exo_steady_state, exo_det_steady_state) % ------------------------------------------------------------------------- % Objective function to compute numerically the Jacobians used for identification analysis % Previously this function was called thet2tau.m % ========================================================================= % INPUTS -% params: [vector] parameter values at which to evaluate objective function -% stderr parameters come first, corr parameters second, model parameters third -% outputflag: [integer] flag which objective to compute (see below) -% 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 -% indpmodel: [vector] Index of model parameters -% indpstderr: [vector] Index of stderr parameters -% indpcorr: [matrix] Index of corr parameters -% indvar: [vector] Index of selected or observed variables +% params: [vector] parameter values at which to evaluate objective function +% stderr parameters come first, corr parameters second, model parameters third +% outputflag: [integer] flag which objective to compute (see below) +% estim_params: [structure] storing the estimation information +% M: [structure] storing the model information +% options: [structure] storing the options +% indpmodel: [vector] Index of model parameters +% indpstderr: [vector] Index of stderr parameters +% indvar: [vector] Index of selected or observed variables +% dr [structure] Reduced form model. +% endo_steady_state [vector] steady state value for endogenous variables +% exo_steady_state [vector] steady state value for exogenous variables +% exo_det_steady_state [vector] steady state value for exogenous deterministic variables % ------------------------------------------------------------------------- % OUTPUTS % out: dependent on outputflag @@ -76,9 +78,9 @@ else end %% compute Kalman transition matrices and steady state with updated parameters -[oo.dr,info,M.params] = compute_decision_rules(M,options,oo.dr, oo.steady_state, oo.exo_steady_state, oo.exo_det_steady_state); +[dr,info,M.params] = compute_decision_rules(M,options,dr, steady_state, exo_steady_state, exo_det_steady_state); options = rmfield(options,'options_ident'); -pruned = pruned_state_space_system(M, options, oo.dr, indvar, nlags, useautocorr, 0); +pruned = pruned_state_space_system(M, options, dr, indvar, nlags, useautocorr, 0); %% 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 diff --git a/matlab/set_dynare_seed.m b/matlab/set_dynare_seed.m index 9791bded4..4050e6555 100644 --- a/matlab/set_dynare_seed.m +++ b/matlab/set_dynare_seed.m @@ -1,9 +1,8 @@ function set_dynare_seed(varargin) -% Set seeds depending on matlab (octave) version. This routine is called is +% set_dynare_seed(varargin) +% Set seeds depending on Matlab (Octave) version. This routine is % a wrapper for set_dynare_seed_local_options -% dynare_config and can be called by the -% user in the mod file. -% + % Copyright © 2010-2023 Dynare Team % % This file is part of Dynare. diff --git a/tests/analytic_derivatives/BrockMirman_PertParamsDerivs.mod b/tests/analytic_derivatives/BrockMirman_PertParamsDerivs.mod index e1acfb16a..bcfc7735f 100644 --- a/tests/analytic_derivatives/BrockMirman_PertParamsDerivs.mod +++ b/tests/analytic_derivatives/BrockMirman_PertParamsDerivs.mod @@ -460,7 +460,7 @@ for jj = 1:2 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); + DERIVS = get_perturbation_params_derivs(M_, options_, estim_params_, oo_.dr, oo_.steady_state, oo_.exo_steady_state, oo_.exo_det_steady_state, 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); diff --git a/tests/analytic_derivatives/burnside_3_order_PertParamsDerivs.mod b/tests/analytic_derivatives/burnside_3_order_PertParamsDerivs.mod index da9951e12..c63091bf3 100644 --- a/tests/analytic_derivatives/burnside_3_order_PertParamsDerivs.mod +++ b/tests/analytic_derivatives/burnside_3_order_PertParamsDerivs.mod @@ -170,7 +170,7 @@ 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); + DERIVS = get_perturbation_params_derivs(M_, options_, estim_params_, oo_.dr, oo_.steady_state, oo_.exo_steady_state, oo_.exo_det_steady_state, 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);