From 6b2cee20177d438b688f67c34014793ebb00ad34 Mon Sep 17 00:00:00 2001 From: Marco Ratto Date: Sat, 5 Nov 2011 10:32:37 +0100 Subject: [PATCH 01/13] Added analytic derivatives for prior distributions. --- matlab/lpdfgam.m | 20 +++++++++++++++++++- matlab/lpdfgbeta.m | 20 +++++++++++++++++++- matlab/lpdfig1.m | 21 +++++++++++++++++++-- matlab/lpdfig2.m | 18 +++++++++++++++++- matlab/lpdfnorm.m | 12 ++++++++++-- 5 files changed, 84 insertions(+), 7 deletions(-) diff --git a/matlab/lpdfgam.m b/matlab/lpdfgam.m index 0a8505960..1246f2a9e 100644 --- a/matlab/lpdfgam.m +++ b/matlab/lpdfgam.m @@ -1,4 +1,4 @@ -function ldens = lpdfgam(x,a,b); +function [ldens,Dldens,D2ldens] = lpdfgam(x,a,b); % Evaluates the logged GAMMA PDF at x. % % INPUTS @@ -37,4 +37,22 @@ if length(a)==1 ldens(idx) = -gammaln(a) - a*log(b) + (a-1)*log(x(idx)) - x(idx)/b ; else ldens(idx) = -gammaln(a(idx)) - a(idx).*log(b(idx)) + (a(idx)-1).*log(x(idx)) - x(idx)./b(idx) ; +end + + + +if nargout >1 + if length(a)==1 + Dldens(idx) = (a-1)./(x(idx)) - ones(length(idx),1)/b ; + else + Dldens(idx) = (a(idx)-1)./(x(idx)) - ones(length(idx),1)./b(idx) ; + end +end + +if nargout == 3 + if length(a)==1 + D2ldens(idx) = -(a-1)./(x(idx)).^2; + else + D2ldens(idx) = -(a(idx)-1)./(x(idx)).^2; + end end \ No newline at end of file diff --git a/matlab/lpdfgbeta.m b/matlab/lpdfgbeta.m index 33ff3df48..1906161e5 100644 --- a/matlab/lpdfgbeta.m +++ b/matlab/lpdfgbeta.m @@ -1,4 +1,4 @@ -function ldens = lpdfgbeta(x,a,b,aa,bb); +function [ldens,Dldens,D2ldens] = lpdfgbeta(x,a,b,aa,bb); % Evaluates the logged BETA PDF at x. % % INPUTS @@ -38,4 +38,22 @@ if length(a)==1 ldens(idx) = -betaln(a,b) + (a-1)*log(x(idx)-aa) + (b-1)*log(bb-x(idx)) - (a+b-1)*log(bb-aa) ; else ldens(idx) = -betaln(a(idx),b(idx)) + (a(idx)-1).*log(x(idx)-aa(idx)) + (b(idx)-1).*log(bb(idx)-x(idx)) - (a(idx)+b(idx)-1).*log(bb(idx)-aa(idx)); +end + + +if nargout >1 + if length(a)==1 + Dldens(idx) = (a-1)./(x(idx)-aa) - (b-1)./(bb-x(idx)) ; + else + Dldens(idx) = (a(idx)-1)./(x(idx)-aa(idx)) - (b(idx)-1)./(bb(idx)-x(idx)); + end +end + + +if nargout == 3 + if length(a)==1 + D2ldens(idx) = -(a-1)./(x(idx)-aa).^2 + (b-1)./(bb-x(idx)).^2 ; + else + D2ldens(idx) = -(a(idx)-1)./(x(idx)-aa(idx)).^2 + (b(idx)-1)./(bb(idx)-x(idx)).^2; + end end \ No newline at end of file diff --git a/matlab/lpdfig1.m b/matlab/lpdfig1.m index 9c17277a8..bb0b54078 100644 --- a/matlab/lpdfig1.m +++ b/matlab/lpdfig1.m @@ -1,4 +1,4 @@ -function ldens = lpdfig1(x,s,nu) +function [ldens,Dldens,D2ldens] = lpdfig1(x,s,nu) % Evaluates the logged INVERSE-GAMMA-1 PDF at x. % % X ~ IG1(s,nu) if X = sqrt(Y) where Y ~ IG2(s,nu) and Y = inv(Z) with Z ~ G(nu/2,2/s) (Gamma distribution) @@ -41,4 +41,21 @@ if length(s)==1 ldens(idx) = log(2) - gammaln(.5*nu) - .5*nu*(log(2)-log(s)) - (nu+1)*log(x(idx)) - .5*s./(x(idx).*x(idx)) ; else ldens(idx) = log(2) - gammaln(.5*nu(idx)) - .5*nu(idx).*(log(2)-log(s(idx))) - (nu(idx)+1).*log(x(idx)) - .5*s(idx)./(x(idx).*x(idx)) ; -end \ No newline at end of file +end + +if nargout >1 + if length(s)==1 + Dldens(idx) = - (nu+1)./(x(idx)) + s./(x(idx).^3) ; + else + Dldens(idx) = - (nu(idx)+1)./(x(idx)) + s(idx)./(x(idx).^3) ; + end +end + +if nargout == 3 + if length(s)==1 + D2ldens(idx) = (nu+1)./(x(idx).^2) - 3*s.*x(idx).^2./(x(idx).^4) ; + else + D2ldens(idx) = (nu(idx)+1)./(x(idx).^2) - 3*s.*x(idx).^2./(x(idx).^4) ; + end +end + diff --git a/matlab/lpdfig2.m b/matlab/lpdfig2.m index 7f9558181..63e458db0 100644 --- a/matlab/lpdfig2.m +++ b/matlab/lpdfig2.m @@ -1,4 +1,4 @@ -function ldens = lpdfig2(x,s,nu) +function [ldens,Dldens,D2ldens] = lpdfig2(x,s,nu) % Evaluates the logged INVERSE-GAMMA-2 PDF at x. % % X ~ IG2(s,nu) if X = inv(Z) where Z ~ G(nu/2,2/s) (Gamma distribution) @@ -41,4 +41,20 @@ if length(s)==1 ldens(idx) = -gammaln(.5*nu) - (.5*nu)*(log(2)-log(s)) - .5*(nu+2)*log(x(idx)) -.5*s./x(idx); else ldens(idx) = -gammaln(.5*nu(idx)) - (.5*nu(idx)).*(log(2)-log(s(idx))) - .5*(nu(idx)+2).*log(x(idx)) -.5*s(idx)./x(idx); +end + +if nargout >1 + if length(s)==1 + Dldens(idx) = - .5*(nu+2)./(x(idx)) + .5*s./x(idx).^2; + else + Dldens(idx) = - .5*(nu(idx)+2)./(x(idx)) + .5*s(idx)./x(idx).^2; + end +end + +if nargout == 3 + if length(s)==1 + D2ldens(idx) = .5*(nu+2)./(x(idx)).^2 - s./x(idx).^3; + else + D2ldens(idx) = .5*(nu(idx)+2)./(x(idx)).^2 - s(idx)./x(idx).^3; + end end \ No newline at end of file diff --git a/matlab/lpdfnorm.m b/matlab/lpdfnorm.m index 0e232884a..9afdd2b6d 100644 --- a/matlab/lpdfnorm.m +++ b/matlab/lpdfnorm.m @@ -1,4 +1,4 @@ -function ldens = lpdfnorm(x,a,b) +function [ldens,Dldens,D2ldens] = lpdfnorm(x,a,b) % Evaluates the logged UNIVARIATE GAUSSIAN PDF at x. % % INPUTS @@ -32,4 +32,12 @@ function ldens = lpdfnorm(x,a,b) if nargin<3, b=1; end if nargin<2, a=0; end -ldens = -log(b) -.5*log(2*pi) - .5*((x-a)./b).*((x-a)./b) ; \ No newline at end of file +ldens = -log(b) -.5*log(2*pi) - .5*((x-a)./b).*((x-a)./b) ; + +if nargout >1 + Dldens = - (1/b)*((x-a)/b) ; +end + +if nargout == 3 + D2ldens = - (1/b)^2 ; +end \ No newline at end of file From 8313b158b68d762231c582af2eff3cce2fa87e94 Mon Sep 17 00:00:00 2001 From: Marco Ratto Date: Sat, 5 Nov 2011 10:46:16 +0100 Subject: [PATCH 02/13] Aligned with DsgeLikelihood.m, manual cherry-picks from 30afa5f415c0d9ad998c3b09ea7c919dc6ea02b9 --- matlab/DsgeLikelihood_hh.m | 115 ++++++++++++++++++++++++++----------- 1 file changed, 82 insertions(+), 33 deletions(-) diff --git a/matlab/DsgeLikelihood_hh.m b/matlab/DsgeLikelihood_hh.m index 896b474f6..3f0d0749a 100644 --- a/matlab/DsgeLikelihood_hh.m +++ b/matlab/DsgeLikelihood_hh.m @@ -1,5 +1,5 @@ function [fval,llik,cost_flag,ys,trend_coeff,info] = DsgeLikelihood_hh(xparam1,DynareDataset,DynareOptions,Model,EstimatedParameters,BayesInfo,DynareResults) -% function [fval,cost_flag,ys,trend_coeff,info] = DsgeLikelihood(xparam1,gend,data,data_index,number_of_observations,no_more_missing_observations) +% function [fval,llik,cost_flag,ys,trend_coeff,info] = DsgeLikelihood_hh(xparam1,DynareDataset,DynareOptions,Model,EstimatedParameters,BayesInfo,DynareResults) % Evaluates the posterior kernel of a dsge model. % % INPUTS @@ -11,6 +11,7 @@ function [fval,llik,cost_flag,ys,trend_coeff,info] = DsgeLikelihood_hh(xparam1,D % no_more_missing_observations [integer] % OUTPUTS % fval : value of the posterior kernel at xparam1. +% llik : probabilities at each time point % cost_flag : zero if the function returns a penalty, one otherwise. % ys : steady state of original endogenous variables % trend_coeff : @@ -52,11 +53,14 @@ if nargin==1 return end -fval = []; -ys = []; -trend_coeff = []; -cost_flag = 1; +% Initialization of the returned variables and others... +fval = []; +ys = []; +trend_coeff = []; +cost_flag = 1; llik=NaN; +info = 0; +singularity_flag = 0; if DynareOptions.block == 1 error('DsgeLikelihood_hh:: This routine (called if mode_compute==5) is not compatible with the block option!') @@ -159,8 +163,6 @@ end Model.Sigma_e = Q; Model.H = H; - - %------------------------------------------------------------------------------ % 2. call model setup & reduction program %------------------------------------------------------------------------------ @@ -223,8 +225,31 @@ Y = DynareDataset.data-trend; %------------------------------------------------------------------------------ % 3. Initial condition of the Kalman filter %------------------------------------------------------------------------------ - kalman_algo = DynareOptions.kalman_algo; + +% resetting measurement errors covariance matrix for univariate filters +if (kalman_algo == 2) || (kalman_algo == 4) + if isequal(H,0) + H = zeros(nobs,1); + mmm = mm; + else + if all(all(abs(H-diag(diag(H)))<1e-14))% ie, the covariance matrix is diagonal... + H = diag(H); + mmm = mm; + else + Z = [Z, eye(pp)]; + T = blkdiag(T,zeros(pp)); + Q = blkdiag(Q,H); + R = blkdiag(R,eye(pp)); + Pstar = blkdiag(Pstar,H); + Pinf = blckdiag(Pinf,zeros(pp)); + H = zeros(nobs,1); + mmm = mm+pp; + end + end +end + + diffuse_periods = 0; switch DynareOptions.lik_init case 1% Standard initialization with the steady state of the state equation. @@ -271,20 +296,32 @@ switch DynareOptions.lik_init if isinf(dLIK) % Go to univariate diffuse filter if singularity problem. kalman_algo = 4; + singularity_flag = 1; end end if (kalman_algo==4) % Univariate Diffuse Kalman Filter - if no_correlation_flag - mmm = mm; - else - Z = [Z, eye(pp)]; - T = blkdiag(T,zeros(pp)); - Q = blkdiag(Q,H); - R = blkdiag(R,eye(pp)); - Pstar = blkdiag(Pstar,H); - Pinf = blckdiag(Pinf,zeros(pp)); - mmm = mm+pp; + if singularity_flag + if isequal(H,0) + H = zeros(nobs,1); + mmm = mm; + else + if all(all(abs(H-diag(diag(H)))<1e-14))% ie, the covariance matrix is diagonal... + H = diag(H); + mmm = mm; + else + Z = [Z, eye(pp)]; + T = blkdiag(T,zeros(pp)); + Q = blkdiag(Q,H); + R = blkdiag(R,eye(pp)); + Pstar = blkdiag(Pstar,H); + Pinf = blckdiag(Pinf,zeros(pp)); + H = zeros(nobs,1); + mmm = mm+pp; + end + end + % no need to test again for correlation elements + singularity_flag = 0; end [dLIK,dlik,a,Pstar] = univariate_kalman_filter_d(DynareDataset.missing.aindex,DynareDataset.missing.number_of_observations,DynareDataset.missing.no_more_missing_observations, ... Y, 1, size(Y,2), ... @@ -294,13 +331,13 @@ switch DynareOptions.lik_init diffuse_periods = length(dlik); end case 4% Start from the solution of the Riccati equation. - if kalman_algo ~= 2 + if kalman_algo ~= 2 kalman_algo = 1; end if isequal(H,0) - [err,Pstar] = kalman_steady_state(transpose(T),R*Q*transpose(R),transpose(build_selection_matrix(mf,np,length(mf)))); + [err,Pstar] = kalman_steady_state(transpose(T),R*Q*transpose(R),transpose(build_selection_matrix(Z,np,length(Z)))); else - [err,Pstar] = kalman_steady_state(transpose(T),R*Q*transpose(R),transpose(build_selection_matrix(mf,np,length(mf))),H); + [err,Pstar] = kalman_steady_state(transpose(T),R*Q*transpose(R),transpose(build_selection_matrix(Z,np,length(Z))),H); end if err disp(['DsgeLikelihood:: I am not able to solve the Riccati equation, so I switch to lik_init=1!']); @@ -316,8 +353,6 @@ end % 4. Likelihood evaluation %------------------------------------------------------------------------------ -singularity_flag = 0; - if ((kalman_algo==1) || (kalman_algo==3))% Multivariate Kalman Filter if no_missing_data_flag [LIK,lik] = kalman_filter(Y,diffuse_periods+1,size(Y,2), ... @@ -333,6 +368,11 @@ if ((kalman_algo==1) || (kalman_algo==3))% Multivariate Kalman Filter T,Q,R,H,Z,mm,pp,rr,Zflag,diffuse_periods); end if isinf(LIK) + if kalman_algo == 1 + kalman_algo = 2; + else + kalman_algo = 4; + end singularity_flag = 1; else if DynareOptions.lik_init==3 @@ -342,21 +382,30 @@ if ((kalman_algo==1) || (kalman_algo==3))% Multivariate Kalman Filter end end -if ( (singularity_flag) || (kalman_algo==2) || (kalman_algo==4) )% Univariate Kalman Filter +if ( singularity_flag || (kalman_algo==2) || (kalman_algo==4) ) + % Univariate Kalman Filter + % resetting measurement error covariance matrix when necessary % if singularity_flag - if no_correlation + if isequal(H,0) + H = zeros(nobs,1); mmm = mm; else - Z = [Z, eye(pp)]; - T = blkdiag(T,zeros(pp)); - Q = blkdiag(Q,H); - R = blkdiag(R,eye(pp)); - Pstar = blkdiag(Pstar,H); - Pinf = blckdiag(Pinf,zeros(pp)); - mmm = mm+pp; - a = [a; zeros(pp,1)]; + if all(all(abs(H-diag(diag(H)))<1e-14))% ie, the covariance matrix is diagonal... + H = diag(H); + mmm = mm; + else + Z = [Z, eye(pp)]; + T = blkdiag(T,zeros(pp)); + Q = blkdiag(Q,H); + R = blkdiag(R,eye(pp)); + Pstar = blkdiag(Pstar,H); + Pinf = blckdiag(Pinf,zeros(pp)); + H = zeros(nobs,1); + mmm = mm+pp; + end end end + [LIK,lik] = univariate_kalman_filter(DynareDataset.missing.aindex,DynareDataset.missing.number_of_observations,DynareDataset.missing.no_more_missing_observations,Y,diffuse_periods+1,size(Y,2), ... a,Pstar, ... DynareOptions.kalman_tol, ... From db61c7c144e062b4a901d56914ba725582682c55 Mon Sep 17 00:00:00 2001 From: Marco Ratto Date: Sat, 5 Nov 2011 11:04:08 +0100 Subject: [PATCH 03/13] Extended for second order derivatives and full Hessian. --- matlab/getH.m | 381 ++++++++++++++++++++++++++++++++++++-------------- 1 file changed, 274 insertions(+), 107 deletions(-) diff --git a/matlab/getH.m b/matlab/getH.m index df7b1cea3..a457969bb 100644 --- a/matlab/getH.m +++ b/matlab/getH.m @@ -1,4 +1,4 @@ -function [H, dA, dOm, Hss, gp, info] = getH(A, B, M_,oo_,kronflag,indx,indexo) +function [H, dA, dOm, Hss, gp, d2A, d2Om, H2ss] = getH(A, B, M_,oo_,kronflag,indx,indexo) % computes derivative of reduced form linear model w.r.t. deep params % @@ -31,10 +31,29 @@ yy0=oo_.dr.ys(I); % yy0 = [ yy0; oo_.dr.ys(find(M_.lead_lag_incidence(j,:)))]; % end dyssdtheta=zeros(length(oo_.dr.ys),M_.param_nbr); -df = feval([M_.fname,'_params_derivs'],yy0, oo_.exo_steady_state', ... - M_.params, oo_.dr.ys, 1, dyssdtheta); -[residual, gg1] = feval([M_.fname,'_static'],oo_.dr.ys, oo_.exo_steady_state', M_.params); +if nargout>5, + [residual, gg1, gg2] = feval([M_.fname,'_static'],oo_.dr.ys, oo_.exo_steady_state', M_.params); + [df, gp, d2f] = feval([M_.fname,'_params_derivs'],yy0, oo_.exo_steady_state', ... + M_.params, oo_.dr.ys, 1, dyssdtheta); + d2f = get_all_resid_2nd_derivs(d2f,length(oo_.dr.ys),M_.param_nbr); + d2f = d2f(:,indx,indx); + if isempty(find(gg2)), + for j=1:length(indx), + d2yssdtheta(:,:,j) = -gg1\d2f(:,:,j); + end + else + gam = d2f*0; + for j=1:length(indx), + d2yssdtheta(:,:,j) = -gg1\(d2f(:,:,j)+gam(:,:,j)); + end + end +else + [residual, gg1] = feval([M_.fname,'_static'],oo_.dr.ys, oo_.exo_steady_state', M_.params); + df = feval([M_.fname,'_params_derivs'],yy0, oo_.exo_steady_state', ... + M_.params, oo_.dr.ys, 1, dyssdtheta); +end dyssdtheta = -gg1\df; + if any(any(isnan(dyssdtheta))), [U,T] = schur(gg1); global options_ @@ -45,15 +64,29 @@ if any(any(isnan(dyssdtheta))), [U,T] = ordschur(U,T,e1); T = T(k+1:end,k+1:end); dyssdtheta = -U(:,k+1:end)*(T\U(:,k+1:end)')*df; + if nargout>5, + for j=1:length(indx), + d2yssdtheta(:,:,j) = -U(:,k+1:end)*(T\U(:,k+1:end)')*d2f(:,:,j); + end + end end +if nargout>5, +[df, gp, d2f, gpp] = feval([M_.fname,'_params_derivs'],yy0, oo_.exo_steady_state', ... + M_.params, oo_.dr.ys, 1, dyssdtheta); +H2ss = d2yssdtheta(oo_.dr.order_var,:,:); +[residual, g1, g2, g3] = feval([M_.fname,'_dynamic'],yy0, oo_.exo_steady_state', ... + M_.params, oo_.dr.ys, 1); +else [df, gp] = feval([M_.fname,'_params_derivs'],yy0, oo_.exo_steady_state', ... M_.params, oo_.dr.ys, 1, dyssdtheta); [residual, g1, g2 ] = feval([M_.fname,'_dynamic'],yy0, oo_.exo_steady_state', ... M_.params, oo_.dr.ys, 1); +end Hss = dyssdtheta(oo_.dr.order_var,indx); dyssdtheta = dyssdtheta(I,:); [nr, nc]=size(g2); +[m, nelem]=size(g1); nc = sqrt(nc); ns = max(max(M_.lead_lag_incidence)); gp2 = gp*0; @@ -73,19 +106,9 @@ end gp = gp+gp2; gp = gp(:,:,indx); param_nbr = length(indx); - -% order_var = [oo_.dr.order_var; ... -% [size(oo_dr.ghx,2)+1:size(oo_dr.ghx,2)+size(oo_.dr.transition_auxiliary_variables,1)]' ]; -% [A(order_var,order_var),B(order_var,:)]=dynare_resolve; -% [A,B,ys,info]=dynare_resolve; -% if info(1) > 0 -% H = []; -% A0 = []; -% B0 = []; -% dA = []; -% dOm = []; -% return -% end +if nargout>5, + param_nbr_2 = param_nbr*(param_nbr+1)/2; +end m = size(A,1); n = size(B,2); @@ -93,9 +116,15 @@ n = size(B,2); klen = M_.maximum_endo_lag + M_.maximum_endo_lead + 1; -k1 = M_.lead_lag_incidence(find([1:klen] ~= M_.maximum_endo_lag+1),:); -a = g1(:,nonzeros(k1')); -da = gp(:,nonzeros(k1'),:); +k11 = M_.lead_lag_incidence(find([1:klen] ~= M_.maximum_endo_lag+1),:); +a = g1(:,nonzeros(k11')); +da = gp(:,nonzeros(k11'),:); +if nargout > 5, + nelem = size(g1,2); + g22 = get_all_2nd_derivs(gpp,m,nelem,M_.param_nbr); + g22 = g22(:,:,indx,indx); + d2a = g22(:,nonzeros(k11'),:,:); +end kstate = oo_.dr.kstate; GAM1 = zeros(M_.endo_nbr,M_.endo_nbr); @@ -103,22 +132,26 @@ Dg1 = zeros(M_.endo_nbr,M_.endo_nbr,param_nbr); % nf = find(M_.lead_lag_incidence(M_.maximum_endo_lag+2,:)); % GAM1(:, nf) = -g1(:,M_.lead_lag_incidence(M_.maximum_endo_lag+2,nf)); -k = find(kstate(:,2) == M_.maximum_endo_lag+2 & kstate(:,3)); -GAM1(:, kstate(k,1)) = -a(:,kstate(k,3)); -Dg1(:, kstate(k,1), :) = -da(:,kstate(k,3),:); -k = find(kstate(:,2) > M_.maximum_endo_lag+2 & kstate(:,3)); -kk = find(kstate(:,2) > M_.maximum_endo_lag+2 & ~kstate(:,3)); -nauxe = 0; -if ~isempty(k), - ax(:,k)= -a(:,kstate(k,3)); - ax(:,kk)= 0; - dax(:,k,:)= -da(:,kstate(k,3),:); - dax(:,kk,:)= 0; - nauxe=size(ax,2); - GAM1 = [ax GAM1]; - Dg1 = cat(2,dax,Dg1); - clear ax +k1 = find(kstate(:,2) == M_.maximum_endo_lag+2 & kstate(:,3)); +GAM1(:, kstate(k1,1)) = -a(:,kstate(k1,3)); +Dg1(:, kstate(k1,1), :) = -da(:,kstate(k1,3),:); +if nargout > 5, + D2g1 = zeros(M_.endo_nbr,M_.endo_nbr,param_nbr,param_nbr); + D2g1(:, kstate(k1,1), :, :) = -d2a(:,kstate(k1,3),:,:); end +% k = find(kstate(:,2) > M_.maximum_endo_lag+2 & kstate(:,3)); +% kk = find(kstate(:,2) > M_.maximum_endo_lag+2 & ~kstate(:,3)); +% nauxe = 0; +% if ~isempty(k), +% ax(:,k)= -a(:,kstate(k,3)); +% ax(:,kk)= 0; +% dax(:,k,:)= -da(:,kstate(k,3),:); +% dax(:,kk,:)= 0; +% nauxe=size(ax,2); +% GAM1 = [ax GAM1]; +% Dg1 = cat(2,dax,Dg1); +% clear ax +% end [junk,cols_b,cols_j] = find(M_.lead_lag_incidence(M_.maximum_endo_lag+1, ... @@ -127,84 +160,95 @@ GAM0 = zeros(M_.endo_nbr,M_.endo_nbr); Dg0 = zeros(M_.endo_nbr,M_.endo_nbr,param_nbr); GAM0(:,cols_b) = g1(:,cols_j); Dg0(:,cols_b,:) = gp(:,cols_j,:); +if nargout > 5, + D2g0 = zeros(M_.endo_nbr,M_.endo_nbr,param_nbr,param_nbr); + D2g0(:, cols_b, :, :) = g22(:,cols_j,:,:); +end % GAM0 = g1(:,M_.lead_lag_incidence(M_.maximum_endo_lag+1,:)); -k = find(kstate(:,2) == M_.maximum_endo_lag+1 & kstate(:,4)); +k2 = find(kstate(:,2) == M_.maximum_endo_lag+1 & kstate(:,4)); GAM2 = zeros(M_.endo_nbr,M_.endo_nbr); Dg2 = zeros(M_.endo_nbr,M_.endo_nbr,param_nbr); % nb = find(M_.lead_lag_incidence(1,:)); % GAM2(:, nb) = -g1(:,M_.lead_lag_incidence(1,nb)); -GAM2(:, kstate(k,1)) = -a(:,kstate(k,4)); -Dg2(:, kstate(k,1), :) = -da(:,kstate(k,4),:); -naux = 0; -k = find(kstate(:,2) < M_.maximum_endo_lag+1 & kstate(:,4)); -kk = find(kstate(:,2) < M_.maximum_endo_lag+1 ); -if ~isempty(k), - ax(:,k)= -a(:,kstate(k,4)); - ax = ax(:,kk); - dax(:,k,:)= -da(:,kstate(k,4),:); - dax = dax(:,kk,:); - naux = size(ax,2); - GAM2 = [GAM2 ax]; - Dg2 = cat(2,Dg2,dax); -end - -GAM0 = blkdiag(GAM0,eye(naux)); -Dg0 = cat(1,Dg0,zeros(naux+nauxe,M_.endo_nbr,param_nbr)); -Dg0 = cat(2,Dg0,zeros(m+nauxe,naux,param_nbr)); -Dg0 = cat(2,zeros(m+nauxe,nauxe,param_nbr),Dg0); - -GAM2 = [GAM2 ; A(M_.endo_nbr+1:end,:)]; -Dg2 = cat(1,Dg2,zeros(naux+nauxe,m,param_nbr)); -Dg2 = cat(2,zeros(m+nauxe,nauxe,param_nbr),Dg2); -GAM2 = [zeros(m+nauxe,nauxe) [GAM2; zeros(nauxe,m)]]; -GAM0 = [[zeros(m,nauxe);(eye(nauxe))] [GAM0; zeros(nauxe,m)]]; +GAM2(:, kstate(k2,1)) = -a(:,kstate(k2,4)); +Dg2(:, kstate(k2,1), :) = -da(:,kstate(k2,4),:); +if nargout > 5, + D2g2 = zeros(M_.endo_nbr,M_.endo_nbr,param_nbr,param_nbr); + D2g2(:, kstate(k2,1), :, :) = -d2a(:,kstate(k2,4),:,:); +end% naux = 0; +% k = find(kstate(:,2) < M_.maximum_endo_lag+1 & kstate(:,4)); +% kk = find(kstate(:,2) < M_.maximum_endo_lag+1 ); +% if ~isempty(k), +% ax(:,k)= -a(:,kstate(k,4)); +% ax = ax(:,kk); +% dax(:,k,:)= -da(:,kstate(k,4),:); +% dax = dax(:,kk,:); +% naux = size(ax,2); +% GAM2 = [GAM2 ax]; +% Dg2 = cat(2,Dg2,dax); +% end +% +% GAM0 = blkdiag(GAM0,eye(naux)); +% Dg0 = cat(1,Dg0,zeros(naux+nauxe,M_.endo_nbr,param_nbr)); +% Dg0 = cat(2,Dg0,zeros(m+nauxe,naux,param_nbr)); +% Dg0 = cat(2,zeros(m+nauxe,nauxe,param_nbr),Dg0); +% +% GAM2 = [GAM2 ; A(M_.endo_nbr+1:end,:)]; +% Dg2 = cat(1,Dg2,zeros(naux+nauxe,m,param_nbr)); +% Dg2 = cat(2,zeros(m+nauxe,nauxe,param_nbr),Dg2); +% GAM2 = [zeros(m+nauxe,nauxe) [GAM2; zeros(nauxe,m)]]; +% GAM0 = [[zeros(m,nauxe);(eye(nauxe))] [GAM0; zeros(nauxe,m)]]; GAM3 = -g1(:,length(yy0)+1:end); % GAM3 = -g1(oo_.dr.order_var,length(yy0)+1:end); -GAM3 = [GAM3; zeros(naux+nauxe,size(GAM3,2))]; +% GAM3 = [GAM3; zeros(naux+nauxe,size(GAM3,2))]; % Dg3 = -gp(oo_.dr.order_var,length(yy0)+1:end,:); Dg3 = -gp(:,length(yy0)+1:end,:); -Dg3 = cat(1,Dg3,zeros(naux+nauxe,size(GAM3,2),param_nbr)); - -auxe = zeros(0,1); -k0 = kstate(find(kstate(:,2) >= M_.maximum_endo_lag+2),:);; -i0 = find(k0(:,2) == M_.maximum_endo_lag+2); -for i=M_.maximum_endo_lag+3:M_.maximum_endo_lag+2+M_.maximum_endo_lead, - i1 = find(k0(:,2) == i); - n1 = size(i1,1); - j = zeros(n1,1); - for j1 = 1:n1 - j(j1) = find(k0(i0,1)==k0(i1(j1),1)); - end - auxe = [auxe; i0(j)]; - i0 = i1; +if nargout>5, + D2g3 = -g22(:,length(yy0)+1:end,:,:); + clear g22 d2a end -auxe = [(1:size(auxe,1))' auxe(end:-1:1)]; -n_ir1 = size(auxe,1); -nr = m + n_ir1; -GAM1 = [[GAM1 zeros(size(GAM1,1),naux)]; zeros(naux+nauxe,m+nauxe)]; -GAM1(m+1:end,:) = sparse(auxe(:,1),auxe(:,2),ones(n_ir1,1),n_ir1,nr); -Dg1 = cat(2,Dg1,zeros(M_.endo_nbr,naux,param_nbr)); -Dg1 = cat(1,Dg1,zeros(naux+nauxe,m+nauxe,param_nbr)); +% Dg3 = cat(1,Dg3,zeros(naux+nauxe,size(GAM3,2),param_nbr)); -Ax = A; -A1 = A; -Bx = B; -B1 = B; -for j=1:M_.maximum_endo_lead-1, - A1 = A1*A; - B1 = A*B1; - k = find(kstate(:,2) == M_.maximum_endo_lag+2+j ); - Ax = [A1(kstate(k,1),:); Ax]; - Bx = [B1(kstate(k,1),:); Bx]; -end -Ax = [zeros(m+nauxe,nauxe) Ax]; -A0 = A; -A=Ax; clear Ax A1; -B0=B; -B = Bx; clear Bx B1; +% auxe = zeros(0,1); +% k0 = kstate(find(kstate(:,2) >= M_.maximum_endo_lag+2),:);; +% i0 = find(k0(:,2) == M_.maximum_endo_lag+2); +% for i=M_.maximum_endo_lag+3:M_.maximum_endo_lag+2+M_.maximum_endo_lead, +% i1 = find(k0(:,2) == i); +% n1 = size(i1,1); +% j = zeros(n1,1); +% for j1 = 1:n1 +% j(j1) = find(k0(i0,1)==k0(i1(j1),1)); +% end +% auxe = [auxe; i0(j)]; +% i0 = i1; +% end +% auxe = [(1:size(auxe,1))' auxe(end:-1:1)]; +% n_ir1 = size(auxe,1); +% nr = m + n_ir1; +% GAM1 = [[GAM1 zeros(size(GAM1,1),naux)]; zeros(naux+nauxe,m+nauxe)]; +% GAM1(m+1:end,:) = sparse(auxe(:,1),auxe(:,2),ones(n_ir1,1),n_ir1,nr); +% Dg1 = cat(2,Dg1,zeros(M_.endo_nbr,naux,param_nbr)); +% Dg1 = cat(1,Dg1,zeros(naux+nauxe,m+nauxe,param_nbr)); + +% Ax = A; +% A1 = A; +% Bx = B; +% B1 = B; +% for j=1:M_.maximum_endo_lead-1, +% A1 = A1*A; +% B1 = A*B1; +% k = find(kstate(:,2) == M_.maximum_endo_lag+2+j ); +% Ax = [A1(kstate(k,1),:); Ax]; +% Bx = [B1(kstate(k,1),:); Bx]; +% end +% Ax = [zeros(m+nauxe,nauxe) Ax]; +% A0 = A; +% A=Ax; clear Ax A1; +% B0=B; +% B = Bx; clear Bx B1; m = size(A,1); n = size(B,2); @@ -312,25 +356,35 @@ else % generalized sylvester equation d(:,:,j) = Dg2(:,:,j)-elem(:,:,j)*A; end xx=sylvester3mr(a,b,c,d); + H=zeros(m*m+m*(m+1)/2,param_nbr+length(indexo)); + if nargout>1, + dOm = zeros(m,m,param_nbr+length(indexo)); + dA=zeros(m,m,param_nbr+length(indexo)); + dB=zeros(m,n,param_nbr); + end if ~isempty(indexo), - dSig = zeros(M_.exo_nbr,M_.exo_nbr); + dSig = zeros(M_.exo_nbr,M_.exo_nbr,length(indexo)); for j=1:length(indexo) - dSig(indexo(j),indexo(j)) = 2*sqrt(M_.Sigma_e(indexo(j),indexo(j))); - y = B*dSig*B'; - y = y(nauxe+1:end,nauxe+1:end); - H(:,j) = [zeros((m-nauxe)^2,1); dyn_vech(y)]; + dSig(indexo(j),indexo(j),j) = 2*sqrt(M_.Sigma_e(indexo(j),indexo(j))); + y = B*dSig(:,:,j)*B'; +% y = y(nauxe+1:end,nauxe+1:end); +% H(:,j) = [zeros((m-nauxe)^2,1); dyn_vech(y)]; + H(:,j) = [zeros(m^2,1); dyn_vech(y)]; if nargout>1, dOm(:,:,j) = y; end - dSig(indexo(j),indexo(j)) = 0; +% dSig(indexo(j),indexo(j)) = 0; end end for j=1:param_nbr, x = xx(:,:,j); y = inva * (Dg3(:,:,j)-(elem(:,:,j)-GAM1*x)*B); + if nargout>1, + dB(:,:,j) = y; + end y = y*M_.Sigma_e*B'+B*M_.Sigma_e*y'; - x = x(nauxe+1:end,nauxe+1:end); - y = y(nauxe+1:end,nauxe+1:end); +% x = x(nauxe+1:end,nauxe+1:end); +% y = y(nauxe+1:end,nauxe+1:end); if nargout>1, dA(:,:,j+length(indexo)) = x; dOm(:,:,j+length(indexo)) = y; @@ -353,7 +407,120 @@ else % generalized sylvester equation end +if nargout > 5, + tot_param_nbr = param_nbr + length(indexo); + tot_param_nbr_2 = tot_param_nbr*(tot_param_nbr+1)/2; + d = zeros(m,m,param_nbr_2); + d2A = zeros(m,m,tot_param_nbr,tot_param_nbr); + d2Om = zeros(m,m,tot_param_nbr,tot_param_nbr); + d2B = zeros(m,n,tot_param_nbr,tot_param_nbr); + cc=triu(ones(param_nbr,param_nbr)); + [i2,j2]=find(cc); + cc = blkdiag( zeros(length(indexo),length(indexo)), cc); + [ipar2]=find(cc); + ctot=triu(ones(tot_param_nbr,tot_param_nbr)); + ctot(1:length(indexo),1:length(indexo))=eye(length(indexo)); + [itot2, jtot2]=find(ctot); + jcount=0; + for j=1:param_nbr, + for i=j:param_nbr, + elem1 = (D2g0(:,:,j,i)-D2g1(:,:,j,i)*A); + elem1 = D2g2(:,:,j,i)-elem1*A; + elemj0 = Dg0(:,:,j)-Dg1(:,:,j)*A; + elemi0 = Dg0(:,:,i)-Dg1(:,:,i)*A; + elem2 = -elemj0*xx(:,:,i)-elemi0*xx(:,:,j); + elem2 = elem2 + ( Dg1(:,:,j)*xx(:,:,i) + Dg1(:,:,i)*xx(:,:,j) )*A; + elem2 = elem2 + GAM1*( xx(:,:,i)*xx(:,:,j) + xx(:,:,j)*xx(:,:,i)); + jcount=jcount+1; + d(:,:,jcount) = elem1+elem2; + end + end + xx2=sylvester3mr(a,b,c,d); + jcount = 0; + for j=1:param_nbr, + for i=j:param_nbr, + jcount=jcount+1; + x = xx2(:,:,jcount); + elem1 = (D2g0(:,:,j,i)-D2g1(:,:,j,i)*A); + elem1 = elem1 -( Dg1(:,:,j)*xx(:,:,i) + Dg1(:,:,i)*xx(:,:,j) ); + elemj0 = Dg0(:,:,j)-Dg1(:,:,j)*A-GAM1*xx(:,:,j); + elemi0 = Dg0(:,:,i)-Dg1(:,:,i)*A-GAM1*xx(:,:,i); + elem0 = elemj0*dB(:,:,i)+elemi0*dB(:,:,j); + y = inva * (D2g3(:,:,j,i)-elem0-(elem1-GAM1*x)*B); + d2B(:,:,j+length(indexo),i+length(indexo)) = y; + d2B(:,:,i+length(indexo),j+length(indexo)) = y; + y = y*M_.Sigma_e*B'+B*M_.Sigma_e*y'+ ... + dB(:,:,j)*M_.Sigma_e*dB(:,:,i)'+dB(:,:,i)*M_.Sigma_e*dB(:,:,j)'; +% x = x(nauxe+1:end,nauxe+1:end); +% y = y(nauxe+1:end,nauxe+1:end); + d2A(:,:,j+length(indexo),i+length(indexo)) = x; + d2A(:,:,i+length(indexo),j+length(indexo)) = x; + d2Om(:,:,j+length(indexo),i+length(indexo)) = y; + d2Om(:,:,i+length(indexo),j+length(indexo)) = y; + end + end + if ~isempty(indexo), + d2Sig = zeros(M_.exo_nbr,M_.exo_nbr,length(indexo)); + for j=1:length(indexo) + d2Sig(indexo(j),indexo(j),j) = 2; + y = B*d2Sig(:,:,j)*B'; + d2Om(:,:,j,j) = y; +% y = y(nauxe+1:end,nauxe+1:end); +% H(:,j) = [zeros((m-nauxe)^2,1); dyn_vech(y)]; +% H(:,j) = [zeros(m^2,1); dyn_vech(y)]; +% dOm(:,:,j) = y; + for i=1:param_nbr, + y = dB(:,:,i)*dSig(:,:,j)*B'+B*dSig(:,:,j)*dB(:,:,i)'; + d2Om(:,:,j,i+length(indexo)) = y; + d2Om(:,:,i+length(indexo),j) = y; + end + end + end +end return +function g22 = get_2nd_deriv(gpp,m,n,i,j), +g22=zeros(m,n); +is=find(gpp(:,3)==i); +is=is(find(gpp(is,4)==j)); + +if ~isempty(is), + g22(gpp(is,1),gpp(is,2))=gpp(is,5); +end +return + +function g22 = get_all_2nd_derivs(gpp,m,n,npar), + +g22=zeros(m,n,npar,npar); +% 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)); + g22(gpp(is,1),gpp(is,2),gpp(is,3),gpp(is,4))=gpp(is,5); + g22(gpp(is,1),gpp(is,2),gpp(is,4),gpp(is,3))=gpp(is,5); +end + +return + +function r22 = get_all_resid_2nd_derivs(rpp,m,npar), + +r22=zeros(m,npar,npar); +% c=ones(npar,npar); +% c=triu(c); +% ic=find(c); + +for is=1:length(rpp), +% d=zeros(npar,npar); +% d(rpp(is,2),rpp(is,3))=1; +% indx = find(ic==find(d)); + 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 \ No newline at end of file From 7517b51630d66922f1645bf4759846977ae2f63b Mon Sep 17 00:00:00 2001 From: Marco Ratto Date: Sat, 5 Nov 2011 11:05:10 +0100 Subject: [PATCH 04/13] Added scores and likelihood as optional output arguments --- matlab/AHessian.m | 39 +++++++++++++++++++++++++++++++++------ 1 file changed, 33 insertions(+), 6 deletions(-) diff --git a/matlab/AHessian.m b/matlab/AHessian.m index d03a01889..cc37f5305 100644 --- a/matlab/AHessian.m +++ b/matlab/AHessian.m @@ -1,5 +1,5 @@ -function [AHess] = AHessian(T,R,Q,H,P,Y,DT,DYss,DOm,DH,DP,start,mf,kalman_tol,riccati_tol) -% function [AHess] = AHessian(T,R,Q,H,P,Y,DT,DYss,DOm,DH,DP,start,mf,kalman_tol,riccati_tol) +function [AHess, DLIK, LIK] = AHessian(T,R,Q,H,P,Y,DT,DYss,DOm,DH,DP,start,mf,kalman_tol,riccati_tol) +% function [AHess, DLIK, LIK] = AHessian(T,R,Q,H,P,Y,DT,DYss,DOm,DH,DP,start,mf,kalman_tol,riccati_tol) % % computes the asymptotic hessian matrix of the log-likelihood function of % a state space model (notation as in kalman_filter.m in DYNARE @@ -28,6 +28,7 @@ function [AHess] = AHessian(T,R,Q,H,P,Y,DT,DYss,DOm,DH,DP,start,mf,kalman_tol,ri k = size(DT,3); % number of structural parameters smpl = size(Y,2); % Sample size. + pp = size(Y,1); % Maximum number of observed variables. mm = size(T,2); % Number of state variables. a = zeros(mm,1); % State vector. Om = R*Q*transpose(R); % Variance of R times the vector of structural innovations. @@ -36,6 +37,11 @@ function [AHess] = AHessian(T,R,Q,H,P,Y,DT,DYss,DOm,DH,DP,start,mf,kalman_tol,ri notsteady = 1; % Steady state flag. F_singular = 1; +lik = zeros(smpl,1); % Initialization of the vector gathering the densities. +LIK = Inf; % Default value of the log likelihood. +if nargout > 1, + DLIK = zeros(k,1); % Initialization of the score. +end AHess = zeros(k,k); % Initialization of the Hessian Da = zeros(mm,k); % State vector. Dv = zeros(length(mf),k); @@ -59,17 +65,21 @@ function [AHess] = AHessian(T,R,Q,H,P,Y,DT,DYss,DOm,DH,DP,start,mf,kalman_tol,ri F_singular = 0; iF = inv(F); K = P(:,mf)*iF; + lik(t) = log(det(F))+transpose(v)*iF*v; [DK,DF,DP1] = computeDKalman(T,DT,DOm,P,DP,DH,mf,iF,K); for ii = 1:k Dv(:,ii) = -Da(mf,ii) - DYss(mf,ii); Da(:,ii) = DT(:,:,ii)*(a+K*v) + T*(Da(:,ii)+DK(:,:,ii)*v + K*Dv(:,ii)); + if t>=start && nargout > 1 + DLIK(ii,1) = DLIK(ii,1) + trace( iF*DF(:,:,ii) ) + 2*Dv(:,ii)'*iF*v - v'*(iF*DF(:,:,ii)*iF)*v; + end end vecDPmf = reshape(DP(mf,mf,:),[],k); - iPmf = inv(P(mf,mf)); +% iPmf = inv(P(mf,mf)); if t>=start - AHess = AHess + Dv'*iPmf*Dv + .5*(vecDPmf' * kron(iPmf,iPmf) * vecDPmf); + AHess = AHess + Dv'*iF*Dv + .5*(vecDPmf' * kron(iF,iF) * vecDPmf); end a = T*(a+K*v); P = T*(P-K*P(mf,:))*transpose(T)+Om; @@ -92,13 +102,23 @@ function [AHess] = AHessian(T,R,Q,H,P,Y,DT,DYss,DOm,DH,DP,start,mf,kalman_tol,ri for ii = 1:k Dv(:,ii) = -Da(mf,ii)-DYss(mf,ii); Da(:,ii) = DT(:,:,ii)*(a+K*v) + T*(Da(:,ii)+DK(:,:,ii)*v + K*Dv(:,ii)); + if t>=start && nargout >1 + DLIK(ii,1) = DLIK(ii,1) + trace( iF*DF(:,:,ii) ) + 2*Dv(:,ii)'*iF*v - v'*(iF*DF(:,:,ii)*iF)*v; + end end if t>=start - AHess = AHess + Dv'*iPmf*Dv; + AHess = AHess + Dv'*iF*Dv; end a = T*(a+K*v); + lik(t) = transpose(v)*iF*v; end - AHess = AHess + .5*(smpl+t0-1)*(vecDPmf' * kron(iPmf,iPmf) * vecDPmf); + AHess = AHess + .5*(smpl+t0-1)*(vecDPmf' * kron(iF,iF) * vecDPmf); + if nargout > 1 + for ii = 1:k + DLIK(ii,1) = DLIK(ii,1) + (smpl-t0+1)*trace( iF*DF(:,:,ii) ); + end + end + lik(t0:smpl) = lik(t0:smpl) + log(det(F)); % for ii = 1:k; % for jj = 1:ii % H(ii,jj) = trace(iPmf*(.5*DP(mf,mf,ii)*iPmf*DP(mf,mf,jj) + Dv(:,ii)*Dv(:,jj)')); @@ -107,6 +127,13 @@ function [AHess] = AHessian(T,R,Q,H,P,Y,DT,DYss,DOm,DH,DP,start,mf,kalman_tol,ri end AHess = -AHess; +if nargout > 1, + DLIK = DLIK/2; +end +% adding log-likelihhod constants +lik = (lik + pp*log(2*pi))/2; + +LIK = sum(lik(start:end)); % Minus the log-likelihood. % end of main function function [DK,DF,DP1] = computeDKalman(T,DT,DOm,P,DP,DH,mf,iF,K) From 9ebb86784de2b8572ff12c3802029c6fd601336b Mon Sep 17 00:00:00 2001 From: Marco Ratto Date: Sat, 5 Nov 2011 11:08:00 +0100 Subject: [PATCH 05/13] Extension to full Hessian (partial commit, to be debugged). --- matlab/DsgeLikelihood.m | 60 +++++++++- matlab/get_Hessian.m | 258 ++++++++++++++++++++++++++++++++++++++++ 2 files changed, 316 insertions(+), 2 deletions(-) create mode 100644 matlab/get_Hessian.m diff --git a/matlab/DsgeLikelihood.m b/matlab/DsgeLikelihood.m index bd09d245c..3a7a61485 100644 --- a/matlab/DsgeLikelihood.m +++ b/matlab/DsgeLikelihood.m @@ -439,6 +439,7 @@ end if analytic_derivation no_DLIK = 0; + full_Hess = 0; DLIK = []; AHess = []; if nargin<7 || isempty(derivatives_info) @@ -453,12 +454,25 @@ if analytic_derivation else indparam=[]; end - [dum, DT, DOm, DYss] = getH(A,B,Model,DynareResults,0,indparam,indexo); + + if full_Hess, + [dum, DT, DOm, DYss, dum2, D2T, D2Om, D2Yss] = getH(A, B, Model,DynareResults,0,indparam,indexo); + else + [dum, DT, DOm, DYss] = getH(A, B, Model,DynareResults,0,indparam,indexo); + end else DT = derivatives_info.DT; DOm = derivatives_info.DOm; DYss = derivatives_info.DYss; - if isfield(derivatives_info,'no_DLIK') + if isfield(derivatives_info,'full_Hess'), + full_Hess = derivatives_info.full_Hess; + end + if full_Hess, + D2T = derivatives_info.D2T; + D2Om = derivatives_info.D2Om; + D2Yss = derivatives_info.D2Yss; + end + if isfield(derivatives_info,'no_DLIK'), no_DLIK = derivatives_info.no_DLIK; end clear('derivatives_info'); @@ -471,6 +485,17 @@ if analytic_derivation DH=zeros([size(H),length(xparam1)]); DQ=zeros([size(Q),length(xparam1)]); DP=zeros([size(T),length(xparam1)]); + if full_Hess, + for j=1:size(D2Yss,1), + tmp(j,:,:) = blkdiag(zeros(offset,offset), squeeze(D2Yss(j,:,:))); + end + D2Yss = tmp; + D2T = D2T(iv,iv,:,:); + D2Om = D2Om(iv,iv,:,:); + D2Yss = D2Yss(iv,:,:); + D2H=zeros([size(H),length(xparam1),length(xparam1)]); + D2P=zeros([size(T),length(xparam1),length(xparam1)]); + end for i=1:EstimatedParameters.nvx k =EstimatedParameters.var_exo(i,1); DQ(k,k,i) = 2*sqrt(Q(k,k)); @@ -478,11 +503,23 @@ if analytic_derivation kk = find(abs(dum) < 1e-12); dum(kk) = 0; DP(:,:,i)=dum; + if full_Hess + for j=1:i, + dum = lyapunov_symm(T,D2Om(:,:,i,j),DynareOptions.qz_criterium,DynareOptions.lyapunov_complex_threshold); + kk = (abs(dum) < 1e-12); + dum(kk) = 0; + D2P(:,:,i,j)=dum; + D2P(:,:,j,i)=dum; + end + end end offset = EstimatedParameters.nvx; for i=1:EstimatedParameters.nvn k = EstimatedParameters.var_endo(i,1); DH(k,k,i+offset) = 2*sqrt(H(k,k)); + if full_Hess + D2H(k,k,i+offset,i+offset) = 2; + end end offset = offset + EstimatedParameters.nvn; for j=1:EstimatedParameters.np @@ -490,6 +527,21 @@ if analytic_derivation kk = find(abs(dum) < 1e-12); dum(kk) = 0; DP(:,:,j+offset)=dum; + if full_Hess + DTj = DT(:,:,j+offset); + DPj = dum; + for i=1:j, + DTi = DT(:,:,i+offset); + DPi = DP(:,:,i+offset); + D2Tij = D2T(:,:,i,j); + D2Omij = D2Om(:,:,i,j); + tmp = D2Tij*Pstar*T' + T*Pstar*D2Tij' + DTi*DPj*T' + DTj*DPi*T' + T*DPj*DTi' + T*DPi*DTj' + DTi*Pstar*DTj' + DTj*Pstar*DTi' + D2Omij; + dum = lyapunov_symm(T,tmp,DynareOptions.qz_criterium,DynareOptions.lyapunov_complex_threshold); + dum(abs(dum)<1.e-12) = 0; + D2P(:,:,i+offset,j+offset) = dum; + D2P(:,:,j+offset,i+offset) = dum; + end + end end end @@ -515,6 +567,10 @@ if ((kalman_algo==1) || (kalman_algo==3))% Multivariate Kalman Filter end if nargout==11 [AHess] = AHessian(T,R,Q,H,Pstar,Y,DT,DYss,DOm,DH,DP,start,Z,kalman_tol,riccati_tol); + if full_Hess, + Hess = get_Hessian(T,R,Q,H,Pstar,Y,DT,DYss,DOm,DH,DP,D2T,D2Yss,D2Om,D2H,D2P,start,Z,kalman_tol,riccati_tol); + Hess0 = getHessian(Y,T,DT,D2T, R*Q*transpose(R),DOm,D2Om,Z,DYss,D2Yss); + end end end else diff --git a/matlab/get_Hessian.m b/matlab/get_Hessian.m new file mode 100644 index 000000000..a42b3c812 --- /dev/null +++ b/matlab/get_Hessian.m @@ -0,0 +1,258 @@ +function [Hess] = get_Hessian(T,R,Q,H,P,Y,DT,DYss,DOm,DH,DP,D2T,D2Yss,D2Om,D2H,D2P,start,mf,kalman_tol,riccati_tol) +% function [Hess] = get_Hessian(T,R,Q,H,P,Y,DT,DYss,DOm,DH,DP,D2T,D2Yss,D2Om,D2H,D2P,start,mf,kalman_tol,riccati_tol) +% +% computes the hessian matrix of the log-likelihood function of +% a state space model (notation as in kalman_filter.m in DYNARE +% Thanks to Nikolai Iskrev +% +% NOTE: the derivative matrices (DT,DR ...) are 3-dim. arrays with last +% dimension equal to the number of structural parameters +% NOTE: the derivative matrices (D2T,D2Om ...) are 4-dim. arrays with last +% two dimensions equal to the number of structural parameters + +% Copyright (C) 2011 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 . + + + k = size(DT,3); % number of structural parameters + smpl = size(Y,2); % Sample size. + pp = size(Y,1); % Maximum number of observed variables. + mm = size(T,2); % Number of state variables. + a = zeros(mm,1); % State vector. + Om = R*Q*transpose(R); % Variance of R times the vector of structural innovations. + t = 0; % Initialization of the time index. + oldK = 0; + notsteady = 1; % Steady state flag. + F_singular = 1; + + Hess = zeros(k,k); % Initialization of the Hessian + Da = zeros(mm,k); % State vector. + Dv = zeros(length(mf),k); + D2a = zeros(mm,k,k); % State vector. + D2v = zeros(length(mf),k,k); + + C = zeros(length(mf),mm); + for ii=1:length(mf); C(ii,mf(ii))=1;end % SELECTION MATRIX IN MEASUREMENT EQ. (FOR WHEN IT IS NOT CONSTANT) + dC = zeros(length(mf),mm,k); + d2C = zeros(length(mf),mm,k,k); + + s = zeros(pp,1); % CONSTANT TERM IN MEASUREMENT EQ. (FOR WHEN IT IS NOT CONSTANT) + ds = zeros(pp,1,k); + d2s = zeros(pp,1,k,k); + +% for ii = 1:k +% DOm = DR(:,:,ii)*Q*transpose(R) + R*DQ(:,:,ii)*transpose(R) + R*Q*transpose(DR(:,:,ii)); +% end + + while notsteady & t=start + Hess = Hess + Hesst; + end + a = T*(a+K*v); + P = T*(P-K*P(mf,:))*transpose(T)+Om; + DP = DP1; + D2P = D2P1; + end + notsteady = max(max(abs(K-oldK))) > riccati_tol; + oldK = K; + end + + if F_singular + error('The variance of the forecast error remains singular until the end of the sample') + end + + + if t < smpl + t0 = t+1; + while t < smpl + t = t+1; + v = Y(:,t)-a(mf); + tmp = (a+K*v); + for ii = 1:k, + Dv(:,ii) = -Da(mf,ii)-DYss(mf,ii); + dKi = DK(:,:,ii); + diFi = -iF*DF(:,:,ii)*iF; + dtmpi = Da(:,ii)+dKi*v+K*Dv(:,ii); + + for jj = 1:ii, + dFj = DF(:,:,jj); + diFj = -iF*DF(:,:,jj)*iF; + dKj = DK(:,:,jj); + d2Kij = D2K(:,:,jj,ii); + d2Fij = D2F(:,:,jj,ii); + d2iFij = -diFi*dFj*iF -iF*d2Fij*iF -iF*dFj*diFi; + dtmpj = Da(:,jj)+dKj*v+K*Dv(:,jj); + + d2vij = -D2Yss(mf,jj,ii) - D2a(mf,jj,ii); + d2tmpij = D2a(:,jj,ii) + d2Kij*v + dKj*Dv(:,ii) + dKi*Dv(:,jj) + K*d2vij; + D2a(:,jj,ii) = D2T(:,:,jj,ii)*tmp + DT(:,:,jj)*dtmpi + DT(:,:,ii)*dtmpj + T*d2tmpij; + + Hesst(ii,jj) = getHesst_ij(v,Dv(:,ii),Dv(:,jj),d2vij,iF,diFi,diFj,d2iFij,dFj,d2Fij); + end + Da(:,ii) = DT(:,:,ii)*tmp + T*dtmpi; + end + if t>=start + Hess = Hess + Hesst; + end + a = T*(a+K*v); + end +% Hess = Hess + .5*(smpl+t0-1)*(vecDPmf' * kron(iPmf,iPmf) * vecDPmf); + % for ii = 1:k; + % for jj = 1:ii + % H(ii,jj) = trace(iPmf*(.5*DP(mf,mf,ii)*iPmf*DP(mf,mf,jj) + Dv(:,ii)*Dv(:,jj)')); + % end + % end + end + +Hess = Hess + tril(Hess,-1)'; + +Hess = -Hess/2; +% end of main function + +function Hesst_ij = getHesst_ij(e,dei,dej,d2eij,iS,diSi,diSj,d2iSij,dSj,d2Sij); +% computes (i,j) term in the Hessian + +Hesst_ij = trace(diSi*dSj + iS*d2Sij) + e'*d2iSij*e + 2*(dei'*diSj*e + dei'*iS*dej + e'*diSi*dej + e'*iS*d2eij); + +% end of getHesst_ij + +function [DK,DF,DP1] = computeDKalman(T,DT,DOm,P,DP,DH,mf,iF,K) + + k = size(DT,3); + tmp = P-K*P(mf,:); + +for ii = 1:k + DF(:,:,ii) = DP(mf,mf,ii) + DH(:,:,ii); + DiF(:,:,ii) = -iF*DF(:,:,ii)*iF; + DK(:,:,ii) = DP(:,mf,ii)*iF + P(:,mf)*DiF(:,:,ii); + Dtmp = DP(:,:,ii) - DK(:,:,ii)*P(mf,:) - K*DP(mf,:,ii); + DP1(:,:,ii) = DT(:,:,ii)*tmp*T' + T*Dtmp*T' + T*tmp*DT(:,:,ii)' + DOm(:,:,ii); +end + +% end of computeDKalman + +function [d2K,d2S,d2P1] = computeD2Kalman(A,dA,d2A,d2Om,P0,dP0,d2P0,DH,mf,iF,K0,dK0); +% computes the second derivatives of the Kalman matrices +% note: A=T in main func. + + k = size(dA,3); + tmp = P0-K0*P0(mf,:); +[ns,no] = size(K0); + +% CPC = C*P0*C'; CPC = .5*(CPC+CPC');iF = inv(CPC); +% APC = A*P0*C'; +% APA = A*P0*A'; + + +d2K = zeros(ns,no,k,k); +d2S = zeros(no,no,k,k); +d2P1 = zeros(ns,ns,k,k); + +for ii = 1:k + dAi = dA(:,:,ii); + dFi = dP0(mf,mf,ii); + d2Omi = d2Om(:,:,ii); + diFi = -iF*dFi*iF; + dKi = dK0(:,:,ii); + for jj = 1:k + dAj = dA(:,:,jj); + dFj = dP0(mf,mf,jj); + d2Omj = d2Om(:,:,jj); + dFj = dP0(mf,mf,jj); + diFj = -iF*dFj*iF; + dKj = dK0(:,:,jj); + + d2Aij = d2A(:,:,jj,ii); + d2Pij = d2P0(:,:,jj,ii); + d2Omij = d2Om(:,:,jj,ii); + + % second order + + d2Fij = d2Pij(mf,mf) ; + +% d2APC = d2Aij*P0*C' + A*d2Pij*C' + A*P0*d2Cij' + dAi*dPj*C' + dAj*dPi*C' + A*dPj*dCi' + A*dPi*dCj' + dAi*P0*dCj' + dAj*P0*dCi'; + d2APC = d2Pij(:,mf); + + d2iF = -diFi*dFj*iF -iF*d2Fij*iF -iF*dFj*diFi; + + d2Kij= d2Pij(:,mf)*iF + P0(:,mf)*d2iF + dP0(:,mf,jj)*diFi + dP0(:,mf,ii)*diFj; + + d2KCP = d2Kij*P0(mf,:) + K0*d2Pij(mf,:) + dKi*dP0(mf,:,jj) + dKj*dP0(mf,:,ii) ; + + dtmpi = dP0(:,:,ii) - dK0(:,:,ii)*P0(mf,:) - K0*dP0(mf,:,ii); + dtmpj = dP0(:,:,jj) - dK0(:,:,jj)*P0(mf,:) - K0*dP0(mf,:,jj); + d2tmp = d2Pij - d2KCP; + + d2AtmpA = d2Aij*tmp*A' + A*d2tmp*A' + A*tmp*d2Aij' + dAi*dtmpj*A' + dAj*dtmpi*A' + A*dtmpj*dAi' + A*dtmpi*dAj' + dAi*tmp*dAj' + dAj*tmp*dAi'; + + d2K(:,:,ii,jj) = d2Kij; %#ok + d2P1(:,:,ii,jj) = d2AtmpA + d2Omij; %#ok<*NASGU> + d2S(:,:,ii,jj) = d2Fij; +% d2iS(:,:,ii,jj) = d2iF; + end +end + +% end of computeD2Kalman + + + \ No newline at end of file From 9be0749e1c636d82cbb6599a91a79bd309c82732 Mon Sep 17 00:00:00 2001 From: Marco Ratto Date: Sat, 5 Nov 2011 11:09:10 +0100 Subject: [PATCH 06/13] small bug fix --- matlab/disp_identification.m | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/matlab/disp_identification.m b/matlab/disp_identification.m index 71c226f8a..057981de8 100644 --- a/matlab/disp_identification.m +++ b/matlab/disp_identification.m @@ -188,7 +188,7 @@ if any(idemoments.ino), iweak = length(find(idemoments.jweak_pair(:,j))); if iweak, [jx,jy]=find(jmap_pair==j); - jstore=[jstore jx(1) jy(1)]'; + jstore=[jstore' jx(1) jy(1)]'; if SampleSize > 1 disp([' [',name{jx(1)},',',name{jy(1)},'] are PAIRWISE collinear (with tol = 1.e-10) for ',num2str((iweak)/SampleSize*100),'% of MC runs!' ]) else From 5bb2b2faadbe11960f08756594af7ecb91c6e9e2 Mon Sep 17 00:00:00 2001 From: Marco Ratto Date: Mon, 7 Nov 2011 09:12:54 +0100 Subject: [PATCH 07/13] Fixed call to mode_check --- matlab/dynare_estimation_1.m | 2 +- matlab/mode_check.m | 4 ++-- 2 files changed, 3 insertions(+), 3 deletions(-) diff --git a/matlab/dynare_estimation_1.m b/matlab/dynare_estimation_1.m index 60edfb505..a37a89d40 100644 --- a/matlab/dynare_estimation_1.m +++ b/matlab/dynare_estimation_1.m @@ -353,7 +353,7 @@ if ~options_.mh_posterior_mode_estimation && options_.cova_compute end if options_.mode_check == 1 && ~options_.mh_posterior_mode_estimation - mode_check('objective_function',xparam1,hh,dataset_,options_,M_,estim_params_,bayestopt_,oo_); + mode_check(objective_function,xparam1,hh,dataset_,options_,M_,estim_params_,bayestopt_,oo_); end if ~options_.mh_posterior_mode_estimation diff --git a/matlab/mode_check.m b/matlab/mode_check.m index 8e14a3e20..36d516c89 100644 --- a/matlab/mode_check.m +++ b/matlab/mode_check.m @@ -1,4 +1,4 @@ -function mode_check(func,x,fval,hessian,DynareDataset,DynareOptions,Model,EstimatedParameters,BayesInfo,DynareResults) +function mode_check(fun,x,hessian,DynareDataset,DynareOptions,Model,EstimatedParameters,BayesInfo,DynareResults) % function mode_check(x,fval,hessian,gend,data,lb,ub) % Checks the maximum likelihood mode @@ -91,7 +91,7 @@ for plt = 1:nbplt, end for i=1:length(z) xx(kk) = z(i); - [fval, exit_flag] = feval(fun,x,DynareDataset,DynareOptions,Model,EstimatedParameters,BayesInfo,DynareResults); + [fval, exit_flag] = feval(fun,xx,DynareDataset,DynareOptions,Model,EstimatedParameters,BayesInfo,DynareResults); if exit_flag y(i,1) = fval; else From 4acd50bdaf9a588b7a55ce3f0df00979936e1b83 Mon Sep 17 00:00:00 2001 From: Marco Ratto Date: Mon, 7 Nov 2011 09:14:26 +0100 Subject: [PATCH 08/13] When using steady state in the old format, M_ should be updated beforehand. --- matlab/evaluate_steady_state_file.m | 2 ++ 1 file changed, 2 insertions(+) diff --git a/matlab/evaluate_steady_state_file.m b/matlab/evaluate_steady_state_file.m index 9671ff06d..5ee03c0b5 100644 --- a/matlab/evaluate_steady_state_file.m +++ b/matlab/evaluate_steady_state_file.m @@ -41,6 +41,8 @@ function [ys,params1,check] = evaluate_steady_state_file(ys_init,exo_ss,params,f if steadystate_flag == 1 % old format + assignin('base','tmp_00_',params); + evalin('base','M_.params=tmp_00_; clear(''tmp_00_'')'); h_steadystate = str2func([fname '_steadystate']); [ys,check] = h_steadystate(ys_init, exo_ss); params1 = evalin('base','M_.params'); From bd26eb2251d60b065ab0f852c0cf8b6c56bb4421 Mon Sep 17 00:00:00 2001 From: Marco Ratto Date: Mon, 7 Nov 2011 09:17:38 +0100 Subject: [PATCH 09/13] provisions for analytic 1st and 2nd derivatives --- matlab/priordens.m | 37 ++++++++++++++++++++++++++++++++++++- 1 file changed, 36 insertions(+), 1 deletion(-) diff --git a/matlab/priordens.m b/matlab/priordens.m index c1b64fe1e..62cbda5d7 100644 --- a/matlab/priordens.m +++ b/matlab/priordens.m @@ -1,4 +1,4 @@ -function logged_prior_density = priordens(x, pshape, p6, p7, p3, p4,initialization) +function [logged_prior_density, dlprior, d2lprior] = priordens(x, pshape, p6, p7, p3, p4,initialization) % Computes a prior density for the structural parameters of DSGE models % % INPUTS @@ -81,6 +81,11 @@ if tt1 if isinf(logged_prior_density) return end + if nargout == 2, + [tmp, dlprior(id1)]=lpdfgbeta(x(id1),p6(id1),p7(id1),p3(id1),p4(id1)); + elseif nargout == 3 + [tmp, dlprior(id1), d2lprior(id1)]=lpdfgbeta(x(id1),p6(id1),p7(id1),p3(id1),p4(id1)); + end end if tt2 @@ -88,10 +93,20 @@ if tt2 if isinf(logged_prior_density) return end + if nargout == 2, + [tmp, dlprior(id2)]=lpdfgam(x(id2)-p3(id2),p6(id2),p7(id2)); + elseif nargout == 3 + [tmp, dlprior(id2), d2lprior(id2)]=lpdfgam(x(id2)-p3(id2),p6(id2),p7(id2)); + end end if tt3 logged_prior_density = logged_prior_density + sum(lpdfnorm(x(id3),p6(id3),p7(id3))) ; + if nargout == 2, + [tmp, dlprior(id3)]=lpdfnorm(x(id3),p6(id3),p7(id3)); + elseif nargout == 3 + [tmp, dlprior(id3), d2lprior(id3)]=lpdfnorm(x(id3),p6(id3),p7(id3)); + end end if tt4 @@ -99,6 +114,11 @@ if tt4 if isinf(logged_prior_density) return end + if nargout == 2, + [tmp, dlprior(id4)]=lpdfig1(x(id4)-p3(id4),p6(id4),p7(id4)); + elseif nargout == 3 + [tmp, dlprior(id4), d2lprior(id4)]=lpdfig1(x(id4)-p3(id4),p6(id4),p7(id4)); + end end if tt5 @@ -107,6 +127,12 @@ if tt5 return end logged_prior_density = logged_prior_density + sum(log(1./(p4(id5)-p3(id5)))) ; + if nargout >1, + dlprior(id5)=zeros(length(id5),1); + end + if nargout == 3 + d2lprior(id5)=zeros(length(id5),1); + end end if tt6 @@ -114,4 +140,13 @@ if tt6 if isinf(logged_prior_density) return end + if nargout == 2, + [tmp, dlprior(id6)]=lpdfig2(x(id6)-p3(id6),p6(id6),p7(id6)); + elseif nargout == 3 + [tmp, dlprior(id6), d2lprior(id6)]=lpdfig2(x(id6)-p3(id6),p6(id6),p7(id6)); + end +end + +if nargout==3, + d2lprior = diag(d2lprior); end \ No newline at end of file From 434157f611125026c5b569d45fc4346c5293af01 Mon Sep 17 00:00:00 2001 From: Marco Ratto Date: Mon, 7 Nov 2011 09:17:56 +0100 Subject: [PATCH 10/13] errors fixed in 2nd order derivatives --- matlab/lpdfgbeta.m | 4 ++-- matlab/lpdfig1.m | 4 ++-- 2 files changed, 4 insertions(+), 4 deletions(-) diff --git a/matlab/lpdfgbeta.m b/matlab/lpdfgbeta.m index 1906161e5..27a3c7fde 100644 --- a/matlab/lpdfgbeta.m +++ b/matlab/lpdfgbeta.m @@ -52,8 +52,8 @@ end if nargout == 3 if length(a)==1 - D2ldens(idx) = -(a-1)./(x(idx)-aa).^2 + (b-1)./(bb-x(idx)).^2 ; + D2ldens(idx) = -(a-1)./(x(idx)-aa).^2 - (b-1)./(bb-x(idx)).^2 ; else - D2ldens(idx) = -(a(idx)-1)./(x(idx)-aa(idx)).^2 + (b(idx)-1)./(bb(idx)-x(idx)).^2; + D2ldens(idx) = -(a(idx)-1)./(x(idx)-aa(idx)).^2 - (b(idx)-1)./(bb(idx)-x(idx)).^2; end end \ No newline at end of file diff --git a/matlab/lpdfig1.m b/matlab/lpdfig1.m index bb0b54078..c0bb8ff77 100644 --- a/matlab/lpdfig1.m +++ b/matlab/lpdfig1.m @@ -53,9 +53,9 @@ end if nargout == 3 if length(s)==1 - D2ldens(idx) = (nu+1)./(x(idx).^2) - 3*s.*x(idx).^2./(x(idx).^4) ; + D2ldens(idx) = (nu+1)./(x(idx).^2) - 3*s(idx)./(x(idx).^4) ; else - D2ldens(idx) = (nu(idx)+1)./(x(idx).^2) - 3*s.*x(idx).^2./(x(idx).^4) ; + D2ldens(idx) = (nu(idx)+1)./(x(idx).^2) - 3*s(idx)./(x(idx).^4) ; end end From 9d2e153ef22c8fe61781d30212bda2708b974e21 Mon Sep 17 00:00:00 2001 From: Marco Ratto Date: Mon, 7 Nov 2011 09:18:26 +0100 Subject: [PATCH 11/13] bug fixes in analytic likelihood scores. --- matlab/AHessian.m | 2 +- matlab/score.m | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/matlab/AHessian.m b/matlab/AHessian.m index cc37f5305..f1046c8dc 100644 --- a/matlab/AHessian.m +++ b/matlab/AHessian.m @@ -115,7 +115,7 @@ end AHess = AHess + .5*(smpl+t0-1)*(vecDPmf' * kron(iF,iF) * vecDPmf); if nargout > 1 for ii = 1:k - DLIK(ii,1) = DLIK(ii,1) + (smpl-t0+1)*trace( iF*DF(:,:,ii) ); +% DLIK(ii,1) = DLIK(ii,1) + (smpl-t0+1)*trace( iF*DF(:,:,ii) ); end end lik(t0:smpl) = lik(t0:smpl) + log(det(F)); diff --git a/matlab/score.m b/matlab/score.m index 9a41214af..dba04b84c 100644 --- a/matlab/score.m +++ b/matlab/score.m @@ -98,7 +98,7 @@ function [DLIK] = score(T,R,Q,H,P,Y,DT,DYss,DOm,DH,DP,start,mf,kalman_tol,riccat a = T*(a+K*v); end for ii = 1:k - DLIK(ii,1) = DLIK(ii,1) + (smpl-t0+1)*trace( iF*DF(:,:,ii) ); +% DLIK(ii,1) = DLIK(ii,1) + (smpl-t0+1)*trace( iF*DF(:,:,ii) ); end end From 3d2e55274de2c05e7717188bd40c3fbebc816302 Mon Sep 17 00:00:00 2001 From: Marco Ratto Date: Mon, 7 Nov 2011 09:19:36 +0100 Subject: [PATCH 12/13] bug fix + inclusion of prior derivatives in analytic computations. --- matlab/DsgeLikelihood.m | 16 ++++++++++++++-- 1 file changed, 14 insertions(+), 2 deletions(-) diff --git a/matlab/DsgeLikelihood.m b/matlab/DsgeLikelihood.m index 3a7a61485..e2099fca1 100644 --- a/matlab/DsgeLikelihood.m +++ b/matlab/DsgeLikelihood.m @@ -442,7 +442,7 @@ if analytic_derivation full_Hess = 0; DLIK = []; AHess = []; - if nargin<7 || isempty(derivatives_info) + if nargin<8 || isempty(derivatives_info) [A,B,nou,nou,Model,DynareOptions,DynareResults] = dynare_resolve(Model,DynareOptions,DynareResults); if ~isempty(EstimatedParameters.var_exo) indexo=EstimatedParameters.var_exo(:,1); @@ -643,7 +643,19 @@ end % ------------------------------------------------------------------------------ % 5. Adds prior if necessary % ------------------------------------------------------------------------------ -lnprior = priordens(xparam1,BayesInfo.pshape,BayesInfo.p6,BayesInfo.p7,BayesInfo.p3,BayesInfo.p4); +if analytic_derivation + if full_Hess, + [lnprior, dlnprior, d2lnprior] = priordens(xparam1,BayesInfo.pshape,BayesInfo.p6,BayesInfo.p7,BayesInfo.p3,BayesInfo.p4); + AHess = Hess + d2lnprior; + else + [lnprior, dlnprior] = priordens(xparam1,BayesInfo.pshape,BayesInfo.p6,BayesInfo.p7,BayesInfo.p3,BayesInfo.p4); + end + if no_DLIK==0 + DLIK = DLIK - dlnprior'; + end +else + lnprior = priordens(xparam1,BayesInfo.pshape,BayesInfo.p6,BayesInfo.p7,BayesInfo.p3,BayesInfo.p4); +end fval = (likelihood-lnprior); % Update DynareOptions.kalman_algo. From 9d5953403ed675871675e2bc168e9b999cbb60a0 Mon Sep 17 00:00:00 2001 From: Marco Ratto Date: Mon, 7 Nov 2011 09:21:36 +0100 Subject: [PATCH 13/13] allow tighter tolerance in htol if requested by the user. --- matlab/newrat.m | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/matlab/newrat.m b/matlab/newrat.m index bd59de7b9..619dd90a0 100644 --- a/matlab/newrat.m +++ b/matlab/newrat.m @@ -46,7 +46,7 @@ icount=0; nx=length(x); xparam1=x; %ftol0=1.e-6; -htol_base = max(1.e-5, ftol0); +htol_base = max(1.e-7, ftol0); flagit=0; % mode of computation of hessian in each iteration ftol=ftol0; gtol=1.e-3;