diff --git a/matlab/dsge_likelihood.m b/matlab/dsge_likelihood.m
index 1bdc7a9f0285ac2686187ed1acb2c09b71871996..1003c9a714b5b171f4132094981e97b7b383d012 100644
--- a/matlab/dsge_likelihood.m
+++ b/matlab/dsge_likelihood.m
@@ -498,6 +498,7 @@ if analytic_derivation,
end
DLIK = [];
AHess = [];
+ iv = DynareResults.dr.restrict_var_list;
if nargin<8 || isempty(derivatives_info)
[A,B,nou,nou,Model,DynareOptions,DynareResults] = dynare_resolve(Model,DynareOptions,DynareResults);
if ~isempty(EstimatedParameters.var_exo)
@@ -512,14 +513,15 @@ if analytic_derivation,
end
if full_Hess,
- [dum, DT, DOm, DYss, dum2, D2T, D2Om, D2Yss] = getH(A, B, Model,DynareResults,DynareOptions,kron_flag,indparam,indexo);
+ [dum, DT, DOm, DYss, dum2, D2T, D2Om, D2Yss] = getH(A, B, Model,DynareResults,DynareOptions,kron_flag,indparam,indexo,iv);
+ clear dum dum2;
else
- [dum, DT, DOm, DYss] = getH(A, B, Model,DynareResults,DynareOptions,kron_flag,indparam,indexo);
+ [dum, DT, DOm, DYss] = getH(A, B, Model,DynareResults,DynareOptions,kron_flag,indparam,indexo,iv);
end
else
- DT = derivatives_info.DT;
- DOm = derivatives_info.DOm;
- DYss = derivatives_info.DYss;
+ DT = derivatives_info.DT(iv,iv,:);
+ DOm = derivatives_info.DOm(iv,iv,:);
+ DYss = derivatives_info.DYss(iv,:);
if isfield(derivatives_info,'full_Hess'),
full_Hess = derivatives_info.full_Hess;
end
@@ -533,11 +535,7 @@ if analytic_derivation,
end
clear('derivatives_info');
end
- iv = DynareResults.dr.restrict_var_list;
DYss = [zeros(size(DYss,1),offset) DYss];
- DT = DT(iv,iv,:);
- DOm = DOm(iv,iv,:);
- DYss = DYss(iv,:);
DH=zeros([length(H),length(H),length(xparam1)]);
DQ=zeros([size(Q),length(xparam1)]);
DP=zeros([size(T),length(xparam1)]);
@@ -546,26 +544,25 @@ if analytic_derivation,
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)]);
+ D2H=sparse(size(D2Om,1),size(D2Om,2)); %zeros([size(H),length(xparam1),length(xparam1)]);
+ D2P=sparse(size(D2Om,1),size(D2Om,2)); %zeros([size(T),length(xparam1),length(xparam1)]);
+ jcount=0;
end
for i=1:EstimatedParameters.nvx
k =EstimatedParameters.var_exo(i,1);
DQ(k,k,i) = 2*sqrt(Q(k,k));
dum = lyapunov_symm(T,DOm(:,:,i),DynareOptions.qz_criterium,DynareOptions.lyapunov_complex_threshold);
- kk = find(abs(dum) < 1e-12);
- dum(kk) = 0;
+% 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;
+ jcount=jcount+1;
+ dum = lyapunov_symm(T,dyn_unvech(D2Om(:,jcount)),DynareOptions.qz_criterium,DynareOptions.lyapunov_complex_threshold);
+% kk = (abs(dum) < 1e-12);
+% dum(kk) = 0;
+ D2P(:,jcount)=dyn_vech(dum);
+% D2P(:,:,j,i)=dum;
end
end
end
@@ -580,22 +577,23 @@ if analytic_derivation,
offset = offset + EstimatedParameters.nvn;
for j=1:EstimatedParameters.np
dum = lyapunov_symm(T,DT(:,:,j+offset)*Pstar*T'+T*Pstar*DT(:,:,j+offset)'+DOm(:,:,j+offset),DynareOptions.qz_criterium,DynareOptions.lyapunov_complex_threshold);
- kk = find(abs(dum) < 1e-12);
- dum(kk) = 0;
+% 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+offset,
+ jcount=jcount+1;
DTi = DT(:,:,i);
DPi = DP(:,:,i);
- D2Tij = D2T(:,:,i,j+offset);
- D2Omij = D2Om(:,:,i,j+offset);
+ D2Tij = reshape(D2T(:,jcount),size(T));
+ D2Omij = dyn_unvech(D2Om(:,jcount));
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,j+offset) = dum;
- D2P(:,:,j+offset,i) = dum;
+% dum(abs(dum)<1.e-12) = 0;
+ D2P(:,jcount) = dyn_vech(dum);
+% D2P(:,:,j+offset,i) = dum;
end
end
end
@@ -603,6 +601,7 @@ if analytic_derivation,
analytic_deriv_info={analytic_derivation,DT,DYss,DOm,DH,DP,asy_Hess};
else
analytic_deriv_info={analytic_derivation,DT,DYss,DOm,DH,DP,D2T,D2Yss,D2Om,D2H,D2P};
+ clear DT DYss DOm DH DP D2T D2Yss D2Om D2H D2P,
end
else
analytic_deriv_info={0};
diff --git a/matlab/getH.m b/matlab/getH.m
index a0f14bc48604441f2d75fe01aba49ebddf7fd138..0f9fcfeae873e8451ab76f7400b874c3b8866752 100644
--- a/matlab/getH.m
+++ b/matlab/getH.m
@@ -1,4 +1,4 @@
-function [H, dA, dOm, Hss, gp, d2A, d2Om, H2ss] = getH(A, B, M_,oo_,options_,kronflag,indx,indexo)
+function [H, dA, dOm, Hss, gp, d2A, d2Om, H2ss] = getH(A, B, M_,oo_,options_,kronflag,indx,indexo,iv)
% computes derivative of reduced form linear model w.r.t. deep params
%
@@ -22,26 +22,89 @@ function [H, dA, dOm, Hss, gp, d2A, d2Om, H2ss] = getH(A, B, M_,oo_,options_,kro
if nargin<6 || isempty(kronflag), kronflag = 0; end
if nargin<7 || isempty(indx), indx = [1:M_.param_nbr]; end,
if nargin<8 || isempty(indexo), indexo = []; end,
+if nargin<9 || isempty(iv), iv = (1:length(A))'; end,
[I,J]=find(M_.lead_lag_incidence');
yy0=oo_.dr.ys(I);
param_nbr = length(indx);
+tot_param_nbr = param_nbr + length(indexo);
if nargout>5,
param_nbr_2 = param_nbr*(param_nbr+1)/2;
+ tot_param_nbr_2 = tot_param_nbr*(tot_param_nbr+1)/2;
end
m = size(A,1);
+m1=length(iv);
n = size(B,2);
+if kronflag==-1, % perturbation
+ gp=0;
+ fun = 'thet2tau';
+ params0 = M_.params;
+ H = fjaco(fun,[sqrt(diag(M_.Sigma_e(indexo,indexo))); M_.params(indx)], M_, oo_, indx, indexo,0);
+ if nargout>1,
+ dOm = zeros(m1,m1,tot_param_nbr);
+ dA=zeros(m1,m1,tot_param_nbr);
+ Hss=H(iv,length(indexo)+1:end);
+ da = H(m+1:m+m*m,:);
+ dom = H(m+m*m+1:end,:);
+ for j=1:tot_param_nbr,
+ tmp = dyn_unvech(dom(:,j));
+ dOm(:,:,j) = tmp(iv,iv);
+ tmp = reshape(da(:,j),m,m);
+ dA(:,:,j) = tmp(iv,iv);
+ end
+ clear da dom tmp
+ end
+ if nargout>5,
+ H2 = hessian_sparse('thet2tau',[sqrt(diag(M_.Sigma_e(indexo,indexo))); M_.params(indx)], ...
+ options_.gstep,M_, oo_, indx,indexo,0,[],[],[],iv);
+ H2ss = zeros(m1,tot_param_nbr,tot_param_nbr);
+ iax=find(triu(rand(tot_param_nbr,tot_param_nbr)));
+ H2 = H2(:,iax);
+ for j=1:m1,
+ H2ss(j,:,:)=dyn_unvech(full(H2(j,:)));
+ end
+ H2ss=H2ss(:,length(indexo)+1:end,length(indexo)+1:end);
+ d2A = sparse(m1*m1,tot_param_nbr_2);
+ d2Om = sparse(m1*(m1+1)/2,tot_param_nbr_2);
+ d2A(:,:) = H2(m1+1:m1+m1*m1,:);
+ d2Om(:,:) = H2(m1+m1*m1+1:end,:);
+ clear H2
+% tmp0=zeros(m1,m1);
+% tmp0(iv,iv)=1;
+% iax=find(tmp0);
+% d2A=d2a(iax,:);
+% iax=find(dyn_vech(tmp0));
+% d2Om=d2om(iax,:);
+
+ end
+% assignin('base','M_', M_);
+% assignin('base','oo_', oo_);
+ return
+end
+
if kronflag==-2,
if nargout>5,
[residual, g1, g2 ] = feval([M_.fname,'_dynamic'],yy0, oo_.exo_steady_state', ...
M_.params, oo_.dr.ys, 1);
- g22 = hessian('thet2tau',[M_.params(indx)],options_.gstep,M_, oo_, indx,[],-1);
- H2ss=g22(1:M_.endo_nbr,:);
+ g22 = hessian_sparse('thet2tau',[M_.params(indx)],options_.gstep,M_, oo_, indx,[],-1);
+ H2ss=full(g22(1:M_.endo_nbr,:));
H2ss = reshape(H2ss,[M_.endo_nbr param_nbr param_nbr]);
+ for j=1:M_.endo_nbr,
+ H2ss(j,:,:)=dyn_unvech(dyn_vech(H2ss(j,:,:)));
+ end
g22=g22(M_.endo_nbr+1:end,:);
- g22 = reshape(g22,[size(g1) param_nbr param_nbr]);
+ inx=find(g22);
+ gx22=zeros(length(inx),5);
+ for j=1:length(inx),
+ [i1, i2] = ind2sub(size(g22),inx(j));
+ [ig1, ig2] = ind2sub(size(g1),i1);
+ [ip1, ip2] = ind2sub([param_nbr param_nbr],i2);
+ gx22(j,:) = [ig1 ig2 ip1 ip2 g22(inx(j))];
+ end
+ g22 = gx22;
+ clear gx22;
else
[residual, g1 ] = feval([M_.fname,'_dynamic'],yy0, oo_.exo_steady_state', ...
M_.params, oo_.dr.ys, 1);
@@ -117,8 +180,9 @@ if nargout>5,
[df, gp, d2f, gpp, hp] = feval([M_.fname,'_params_derivs'],yy0, oo_.exo_steady_state', ...
M_.params, oo_.dr.ys, 1, dyssdtheta, d2yssdtheta);
H2ss = d2yssdtheta(oo_.dr.order_var,indx,indx);
- nelem=size(g1,2);
- g22 = get_all_2nd_derivs(gpp,m,nelem,M_.param_nbr);
+% nelem=size(g1,2);
+% g22 = get_all_2nd_derivs(gpp,m,nelem,M_.param_nbr);
+% g22 = g22(:,:,indx,indx);
else
[df, gp] = feval([M_.fname,'_params_derivs'],yy0, oo_.exo_steady_state', ...
M_.params, oo_.dr.ys, 1, dyssdtheta,d2yssdtheta);
@@ -148,25 +212,55 @@ gp = gp+gp2;
gp = gp(:,:,indx);
if nargout>5,
- h22 = get_all_hess_derivs(hp,nr,nc,M_.param_nbr);
- gp22 = g22*0;
+% h22 = get_all_hess_derivs(hp,nr,nc,M_.param_nbr);
+ g22 = gpp;
+ gp22 = sparse(nr*nc,param_nbr*param_nbr);
tmp1 = reshape(g3,[nr*nc*nc nc]);
- tmp2=tmp1*[dyssdtheta; zeros(nc-ns,M_.param_nbr)];
- tmp1=reshape(tmp2,[nr nc nc M_.param_nbr]);
-
+ tmp2=sparse(size(tmp1,1),M_.param_nbr);
+% tmp2=tmp1*[dyssdtheta; zeros(nc-ns,M_.param_nbr)];
+ tmpa=[dyssdtheta; zeros(nc-ns,M_.param_nbr)];
+ tmpa=sparse(tmpa);
+ for j=1:M_.param_nbr,
+ tmp2(:,j)=tmp1*tmpa(:,j);
+ end
+% tmp2=sparse(tmp2);
+% [i1 i2]=ind2sub([nc M_.param_nbr],[1:nc*M_.param_nbr]');
+
for j=1:nr,
+ tmp0=reshape(g2(j,:),[nc nc]);
+ tmp0 = tmp0(:,1:ns)*reshape(d2yssdtheta(I,:,:),[ns,M_.param_nbr*M_.param_nbr]);
for i=1:nc,
- gp22(j,i,:,:)=squeeze(tmp1(j,i,:,:))'*[dyssdtheta; zeros(nc-ns,M_.param_nbr)];
- tmpu = squeeze(h22(j,i,:,:))'*[dyssdtheta; zeros(nc-ns,M_.param_nbr)];
- gp22(j,i,:,:)=gp22(j,i,:,:)+reshape(tmpu+tmpu',[1 1 M_.param_nbr M_.param_nbr]);
+ indo = sub2ind([nr nc nc], ones(nc,1)*j ,ones(nc,1)*i, (1:nc)');
+ tmpx = (tmp2(indo,:))'*[dyssdtheta; zeros(nc-ns,M_.param_nbr)];
+% gp22(j,i,:,:)=squeeze(tmp1(j,i,:,:))'*[dyssdtheta; zeros(nc-ns,M_.param_nbr)];
+ tmpu = (get_hess_deriv(hp,j,i,nc,M_.param_nbr))'*[dyssdtheta; zeros(nc-ns,M_.param_nbr)];
+ tmpy = tmpx+tmpu+tmpu'+reshape(tmp0(i,:,:),[M_.param_nbr M_.param_nbr]);
+ tmpy = tmpy + get_2nd_deriv_mat(gpp,j,i,M_.param_nbr);
+ if any(any(tmpy)),
+ tmpy = triu(tmpy(indx,indx));
+ ina = find(tmpy);
+ gp22(sub2ind([nr nc],j,i),ina)=tmpy(ina);
+% gp22(j,i,:,:)= reshape(tmpy,[1 1 M_.param_nbr M_.param_nbr]);
+
+ end
end
- tmp0=reshape(g2(j,:),[nc nc]);
- gp22(j,:,:,:)=gp22(j,:,:,:)+reshape(tmp0(:,1:ns)*d2yssdtheta(I,:),[1 nc M_.param_nbr M_.param_nbr]);
+% gp22(j,:,:,:)=gp22(j,:,:,:)+reshape(tmp0(:,1:ns)*d2yssdtheta(I,:,:),[1 nc M_.param_nbr M_.param_nbr]);
end
- g22 = g22+gp22;
- g22 = g22(:,:,indx,indx);
- clear tmp0 tmp1 tmp2 tmpu,
+% g22 = g22+gp22;
+% g22 = g22(:,:,indx,indx);
+ clear tmp0 tmp1 tmp2 tmpu tmpx tmpy,
+ inx=find(gp22);
+ gx22=zeros(length(inx),5);
+ for j=1:length(inx),
+ [i1, i2] = ind2sub(size(gp22),inx(j));
+ [ig1, ig2] = ind2sub(size(g1),i1);
+ [ip1, ip2] = ind2sub([param_nbr param_nbr],i2);
+ gx22(j,:) = [ig1 ig2 ip1 ip2 gp22(inx(j))];
+ end
+ g22 = gx22;
+ clear gx22 gp22;
+
end
end
@@ -177,7 +271,14 @@ k11 = M_.lead_lag_incidence(find([1:klen] ~= M_.maximum_endo_lag+1),:);
a = g1(:,nonzeros(k11'));
da = gp(:,nonzeros(k11'),:);
if nargout > 5,
- d2a = g22(:,nonzeros(k11'),:,:);
+ indind = ismember(g22(:,2),nonzeros(k11'));
+ tmp = g22(indind,:);
+ d2a=tmp;
+ for j=1:size(tmp,1),
+ inxinx = find(nonzeros(k11')==tmp(j,2));
+ d2a(j,2) = inxinx;
+ end
+% d2a = g22(:,nonzeros(k11'),:,:);
end
kstate = oo_.dr.kstate;
@@ -188,8 +289,14 @@ 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),:,:);
+ indind = ismember(d2a(:,2),kstate(k1,3));
+ tmp = d2a(indind,:);
+ tmp(:,end)=-tmp(:,end);
+ D2g1 = tmp;
+ for j=1:size(tmp,1),
+ inxinx = (kstate(k1,3)==tmp(j,2));
+ D2g1(j,2) = kstate(k1(inxinx),1);
+ end
end
[junk,cols_b,cols_j] = find(M_.lead_lag_incidence(M_.maximum_endo_lag+1, ...
@@ -199,8 +306,13 @@ 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,:,:);
+ indind = ismember(g22(:,2),cols_j);
+ tmp = g22(indind,:);
+ D2g0=tmp;
+ for j=1:size(tmp,1),
+ inxinx = (cols_j==tmp(j,2));
+ D2g0(j,2) = cols_b(inxinx);
+ end
end
@@ -210,17 +322,32 @@ Dg2 = zeros(M_.endo_nbr,M_.endo_nbr,param_nbr);
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),:,:);
+ indind = ismember(d2a(:,2),kstate(k2,4));
+ tmp = d2a(indind,:);
+ tmp(:,end)=-tmp(:,end);
+ D2g2 = tmp;
+ for j=1:size(tmp,1),
+ inxinx = (kstate(k2,4)==tmp(j,2));
+ D2g2(j,2) = kstate(k2(inxinx),1);
+ end
end
GAM3 = -g1(:,length(yy0)+1:end);
Dg3 = -gp(:,length(yy0)+1:end,:);
if nargout>5,
- D2g3 = -g22(:,length(yy0)+1:end,:,:);
- clear g22 d2a
+ cols_ex = [length(yy0)+1:size(g1,2)];
+ indind = ismember(g22(:,2),cols_ex);
+ tmp = g22(indind,:);
+ tmp(:,end)=-tmp(:,end);
+ D2g3=tmp;
+ for j=1:size(tmp,1),
+ inxinx = find(cols_ex==tmp(j,2));
+ D2g3(j,2) = inxinx;
+ end
+ clear g22 d2a tmp
end
+clear g1 g2 g3 df d2f gpp hp residual gg1 gg2 gp2 dyssdtheta d2yssdtheta
if kronflag==1, % kronecker products
Dg0=reshape(Dg0,m^2,param_nbr);
@@ -294,13 +421,6 @@ if kronflag==1, % kronecker products
end
H = [ [zeros(M_.endo_nbr,length(indexo)) Hss]; [Hx H]];
-elseif kronflag==-1, % perturbation
- fun = 'thet2tau';
- params0 = M_.params;
- H = fjaco(fun,[sqrt(diag(M_.Sigma_e(indexo,indexo))); M_.params(indx)], M_, oo_, indx, indexo);
- assignin('base','M_', M_);
- assignin('base','oo_', oo_);
-
else % generalized sylvester equation
% solves a*x+b*x*c=d
@@ -321,10 +441,10 @@ else % generalized sylvester equation
[xx, flag]=sylvester3a(xx,a,b,c,d);
icount=icount+1;
end
- H=zeros(m*m+m*(m+1)/2,param_nbr+length(indexo));
+ H=zeros(m1*m1+m1*(m1+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));
+ dOm = zeros(m1,m1,param_nbr+length(indexo));
+ dA=zeros(m1,m1,param_nbr+length(indexo));
dB=zeros(m,n,param_nbr);
end
if ~isempty(indexo),
@@ -334,9 +454,9 @@ else % generalized sylvester equation
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)];
+ H(:,j) = [zeros(m1^2,1); dyn_vech(y(iv,iv))];
if nargout>1,
- dOm(:,:,j) = y;
+ dOm(:,:,j) = y(iv,iv);
end
% dSig(indexo(j),indexo(j)) = 0;
end
@@ -351,10 +471,10 @@ else % generalized sylvester equation
% 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;
+ dA(:,:,j+length(indexo)) = x(iv,iv);
+ dOm(:,:,j+length(indexo)) = y(iv,iv);
end
- H(:,j+length(indexo)) = [x(:); dyn_vech(y)];
+ H(:,j+length(indexo)) = [vec(x(iv,iv)); dyn_vech(y(iv,iv))];
end
% for j=1:param_nbr,
% disp(['Derivatives w.r.t. ',M_.param_names(indx(j),:),', ',int2str(j),'/',int2str(param_nbr)])
@@ -368,82 +488,132 @@ else % generalized sylvester equation
% y = y(nauxe+1:end,nauxe+1:end);
% H(:,j) = [x(:); dyn_vech(y)];
% end
- H = [[zeros(M_.endo_nbr,length(indexo)) Hss]; H];
+ Hss = Hss(iv,:);
+ H = [[zeros(m1,length(indexo)) Hss]; H];
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);
+ H2ss = H2ss(iv,:,:);
+ d = zeros(m,m,floor(sqrt(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;
+ cumjcount=0;
+ jinx = [];
+ x2x=sparse(m*m,param_nbr_2);
+% x2x=[];
+ for i=1:param_nbr,
+ for j=1:i,
+ elem1 = (get_2nd_deriv(D2g0,m,m,j,i)-get_2nd_deriv(D2g1,m,m,j,i)*A);
+ elem1 = get_2nd_deriv(D2g2,m,m,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;
+ jinx = [jinx; [j i]];
+ d(:,:,jcount) = elem1+elem2;
+ if jcount==floor(sqrt(param_nbr_2)) || (j*i)==param_nbr^2,
+ if (j*i)==param_nbr^2,
+ d = d(:,:,1:jcount);
+ end
+% d(find(abs(d)<1.e-12))=0;
+ xx2=sylvester3(a,b,c,d);
+ flag=1;
+ icount=0;
+ while flag && icount<4,
+ [xx2, flag]=sylvester3a(xx2,a,b,c,d);
+ icount = icount + 1;
+ end
+% inx = find(abs(xx2)>1.e-12);
+% xx2(find(abs(xx2)<1.e-12))=0;
+ x2x(:,cumjcount+1:cumjcount+jcount)=reshape(xx2,[m*m jcount]);
+ cumjcount=cumjcount+jcount;
+% [i1 i2 i3]=ind2sub(size(xx2),inx);
+% x2x = [x2x; [i1 i2 jinx(i3,:) xx2(inx)]];
+ jcount = 0;
+ jinx = [];
+ end
end
end
- xx2=sylvester3(a,b,c,d);
- flag=1;
- icount=0;
- while flag && icount<4,
- [xx2, flag]=sylvester3a(xx2,a,b,c,d);
- icount = icount + 1;
- end
+ clear d xx2;
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;
+ icount = 0;
+ cumjcount = 0;
+ MAX_DIM_MAT = 100000000;
+ ncol = max(1,floor(MAX_DIM_MAT/(8*m1*(m1+1)/2)));
+ ncol = min(ncol, tot_param_nbr_2);
+ d2A = sparse(m1*m1,tot_param_nbr_2);
+ d2Om = sparse(m1*(m1+1)/2,tot_param_nbr_2);
+ d2A_tmp = zeros(m1*m1,ncol);
+ d2Om_tmp = zeros(m1*(m1+1)/2,ncol);
+ tmpDir = CheckPath('tmp_derivs',M_.dname);
+ offset = length(indexo);
+ % d2B = zeros(m,n,tot_param_nbr,tot_param_nbr);
+ d2Sig = zeros(M_.exo_nbr,M_.exo_nbr,length(indexo));
+ for j=1:tot_param_nbr,
+ for i=1:j,
+ jcount=jcount+1;
+ if j<=offset,
+ if i==j,
+ d2Sig(indexo(j),indexo(j),j) = 2;
+ y = B*d2Sig(:,:,j)*B';
+% y(abs(y)<1.e-8)=0;
+ d2Om_tmp(:,jcount) = dyn_vech(y(iv,iv));
+ end
+ else
+ jind = j-offset;
+ iind = i-offset;
+ if i<=offset,
+ y = dB(:,:,jind)*dSig(:,:,i)*B'+B*dSig(:,:,i)*dB(:,:,jind)';
+% y(abs(y)<1.e-8)=0;
+ d2Om_tmp(:,jcount) = dyn_vech(y(iv,iv));
+ else
+ icount=icount+1;
+ x = reshape(x2x(:,icount),[m m]);
+% x = get_2nd_deriv(x2x,m,m,iind,jind);%xx2(:,:,jcount);
+ elem1 = (get_2nd_deriv(D2g0,m,m,iind,jind)-get_2nd_deriv(D2g1,m,m,iind,jind)*A);
+ elem1 = elem1 -( Dg1(:,:,jind)*xx(:,:,iind) + Dg1(:,:,iind)*xx(:,:,jind) );
+ elemj0 = Dg0(:,:,jind)-Dg1(:,:,jind)*A-GAM1*xx(:,:,jind);
+ elemi0 = Dg0(:,:,iind)-Dg1(:,:,iind)*A-GAM1*xx(:,:,iind);
+ elem0 = elemj0*dB(:,:,iind)+elemi0*dB(:,:,jind);
+ y = inva * (get_2nd_deriv(D2g3,m,n,iind,jind)-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(:,:,jind)*M_.Sigma_e*dB(:,:,iind)'+dB(:,:,iind)*M_.Sigma_e*dB(:,:,jind)';
+% x(abs(x)<1.e-8)=0;
+ d2A_tmp(:,jcount) = vec(x(iv,iv));
+% y(abs(y)<1.e-8)=0;
+ d2Om_tmp(:,jcount) = dyn_vech(y(iv,iv));
+ end
+ end
+ if jcount==ncol || i*j==tot_param_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)],'d2A')
+ save([tmpDir filesep 'd2Om_' int2str(cumjcount+1) '_' int2str(cumjcount+jcount)],'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(m1*m1,ncol);
+ d2Om_tmp = zeros(m1*(m1+1)/2,ncol);
+
end
end
end
@@ -462,9 +632,28 @@ if ~isempty(is),
end
return
-function g22 = get_all_2nd_derivs(gpp,m,n,npar),
+function g22 = get_2nd_deriv_mat(gpp,i,j,n),
+
+g22=zeros(n,n);
+is=find(gpp(:,1)==i);
+is=is(find(gpp(is,2)==j));
+
+if ~isempty(is),
+ g22(sub2ind([n,n],gpp(is,3),gpp(is,4)))=gpp(is,5)';
+ g22(sub2ind([n,n],gpp(is,4),gpp(is,3)))=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);
@@ -473,9 +662,14 @@ 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);
+ g22(sub2ind([m,n],gpp(is,1),gpp(is,2)),sub2ind([npar,npar],gpp(is,4),gpp(is,3)))=gpp(is,5);
+ else
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
+end
return
@@ -506,3 +700,21 @@ for is=1:length(hp),
end
return
+
+function h2 = get_hess_deriv(hp,i,j,m,npar),
+
+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
+
+is=is1(find(hp(is1,3)==j));
+
+if ~isempty(is),
+ h2(sub2ind([m,npar],hp(is,2),hp(is,4)))=hp(is,5)';
+end
+
+return
diff --git a/matlab/hessian_sparse.m b/matlab/hessian_sparse.m
new file mode 100644
index 0000000000000000000000000000000000000000..eb447f8d34d9ac9c9237c021dba01762a708a911
--- /dev/null
+++ b/matlab/hessian_sparse.m
@@ -0,0 +1,81 @@
+function hessian_mat = hessian_sparse(func,x,gstep,varargin)
+% function hessian_mat = hessian_sparse(func,x,gstep,varargin)
+% Computes second order partial derivatives
+%
+% INPUTS
+% func [string] name of the function
+% x [double] vector, the Hessian of "func" is evaluated at x.
+% gstep [double] scalar, size of epsilon.
+% varargin [void] list of additional arguments for "func".
+%
+% OUTPUTS
+% hessian_mat [double, sparse] Hessian matrix
+%
+% ALGORITHM
+% Uses Abramowitz and Stegun (1965) formulas 25.3.24 and 25.3.27 p. 884
+%
+% SPECIAL REQUIREMENTS
+% none
+%
+
+% Copyright (C) 2001-2012 Dynare Team
+%
+% This file is part of Dynare.
+%
+% Dynare is free software: you can redistribute it and/or modify
+% it under the terms of the GNU General Public License as published by
+% the Free Software Foundation, either version 3 of the License, or
+% (at your option) any later version.
+%
+% Dynare is distributed in the hope that it will be useful,
+% but WITHOUT ANY WARRANTY; without even the implied warranty of
+% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+% GNU General Public License for more details.
+%
+% You should have received a copy of the GNU General Public License
+% along with Dynare. If not, see <http://www.gnu.org/licenses/>.
+
+if ~isa(func, 'function_handle')
+ func = str2func(func);
+end
+n=size(x,1);
+h1=max(abs(x),sqrt(gstep(1))*ones(n,1))*eps^(1/6)*gstep(2);
+h_1=h1;
+xh1=x+h1;
+h1=xh1-x;
+xh1=x-h_1;
+h_1=x-xh1;
+xh1=x;
+f0=feval(func,x,varargin{:});
+f1=zeros(size(f0,1),n);
+f_1=f1;
+for i=1:n
+ xh1(i)=x(i)+h1(i);
+ f1(:,i)=feval(func,xh1,varargin{:});
+ xh1(i)=x(i)-h_1(i);
+ f_1(:,i)=feval(func,xh1,varargin{:});
+ xh1(i)=x(i);
+end
+xh_1=xh1;
+hessian_mat = sparse(size(f0,1),n*n);
+
+for i=1:n
+% if i > 1
+% k=[i:n:n*(i-1)];
+% hessian_mat(:,(i-1)*n+1:(i-1)*n+i-1)=hessian_mat(:,k);
+% hessian_mat(:,k)=0;
+% end
+ hessian_mat(:,(i-1)*n+i)=(f1(:,i)+f_1(:,i)-2*f0)./(h1(i)*h_1(i));
+ temp=f1+f_1-f0*ones(1,n);
+ for j=1:i-1
+ xh1(i)=x(i)+h1(i);
+ xh1(j)=x(j)+h_1(j);
+ xh_1(i)=x(i)-h1(i);
+ xh_1(j)=x(j)-h_1(j);
+ hessian_mat(:,(i-1)*n+j)=-(-feval(func,xh1,varargin{:})-feval(func,xh_1,varargin{:})+temp(:,i)+temp(:,j))./(2*h1(i)*h_1(j));
+ xh1(i)=x(i);
+ xh1(j)=x(j);
+ xh_1(i)=x(i);
+ xh_1(j)=x(j);
+ end
+end
\ No newline at end of file
diff --git a/matlab/kalman/likelihood/computeDLIK.m b/matlab/kalman/likelihood/computeDLIK.m
index 4efc95b33e89e6b78f94bb94279950c5747df830..1acb0796ff167c6112be36ad521eb0051a391635 100644
--- a/matlab/kalman/likelihood/computeDLIK.m
+++ b/matlab/kalman/likelihood/computeDLIK.m
@@ -1,4 +1,4 @@
-function [Da,DP1,DLIK,D2a,D2P1,Hesst] = computeDLIK(k,tmp,Z,Zflag,v,T,K,P,iF,Da,DYss,DT,DOm,DP,DH,notsteady,D2a,D2Yss,D2T,D2Om,D2P)
+function [Da,DP,DLIK,D2a,D2P,Hesst] = computeDLIK(k,tmp,Z,Zflag,v,T,K,P,iF,Da,DYss,DT,DOm,DP,DH,notsteady,D2a,D2Yss,D2T,D2Om,D2P)
% Copyright (C) 2012 Dynare Team
%
@@ -25,18 +25,20 @@ if notsteady
if Zflag
[DK,DF,DP1] = computeDKalmanZ(T,DT,DOm,P,DP,DH,Z,iF,K);
if nargout>4,
- [D2K,D2F,D2P1] = computeD2KalmanZ(T,DT,D2T,D2Om,P,DP,D2P,DH,Z,iF,K,DK);
+ [D2K,D2F,D2P] = computeD2KalmanZ(T,DT,D2T,D2Om,P,DP,D2P,DH,Z,iF,K,DK);
end
else
[DK,DF,DP1] = computeDKalman(T,DT,DOm,P,DP,DH,Z,iF,K);
if nargout>4,
- [D2K,D2F,D2P1] = computeD2Kalman(T,DT,D2T,D2Om,P,DP,D2P,DH,Z,iF,K,DK);
+ [D2K,D2F,D2P] = computeD2Kalman(T,DT,D2T,D2Om,P,DP,D2P,DH,Z,iF,K,DK);
end
end
+DP=DP1;
+clear DP1,
else
- DP1=DP;
+ DP=DP;
if nargout>4,
- D2P1=D2P;
+ D2P=D2P;
end
end
@@ -64,6 +66,7 @@ end
Hesst = zeros(k,k);
DLIK=zeros(k,1);
+jcount=0;
for ii = 1:k
% dai = da(:,:,ii);
dKi = DK(:,:,ii);
@@ -72,6 +75,7 @@ for ii = 1:k
if nargout>4,
diFi = -iF*DF(:,:,ii)*iF;
for jj = 1:ii
+ jcount=jcount+1;
dFj = DF(:,:,jj);
diFj = -iF*DF(:,:,jj)*iF;
dKj = DK(:,:,jj);
@@ -87,7 +91,7 @@ for ii = 1:k
d2vij = -D2Yss(Z,jj,ii) - D2a(Z,jj,ii);
end
d2tmpij = D2a(:,jj,ii) + d2Kij*v + dKj*Dv(:,ii) + dKi*Dv(:,jj) + K*d2vij;
- D2a(:,jj,ii) = D2T(:,:,jj,ii)*tmp + DT(:,:,jj)*dtmp(:,ii) + DT(:,:,ii)*dtmp(:,jj) + T*d2tmpij;
+ D2a(:,jj,ii) = reshape(D2T(:,jcount),size(T))*tmp + DT(:,:,jj)*dtmp(:,ii) + DT(:,:,ii)*dtmp(:,jj) + T*d2tmpij;
D2a(:,ii,jj) = D2a(:,jj,ii);
if nargout==6,
@@ -125,19 +129,19 @@ Hesst_ij = trace(diSi*dSj + iS*d2Sij) + e'*d2iSij*e + 2*(dei'*diSj*e + dei'*iS*d
% end of getHesst_ij
-function [DK,DF,DP1] = computeDKalman(T,DT,DOm,P,DP,DH,Z,iF,K)
+function [DK,DF,DP1] = computeDKalman(T,DT,DOm,P,DP1,DH,Z,iF,K)
k = size(DT,3);
tmp = P-K*P(Z,:);
DF = zeros([size(iF),k]);
DK = zeros([size(K),k]);
-DP1 = zeros([size(P),k]);
+% DP1 = zeros([size(P),k]);
for ii = 1:k
- DF(:,:,ii) = DP(Z,Z,ii) + DH(:,:,ii);
+ DF(:,:,ii) = DP1(Z,Z,ii) + DH(:,:,ii);
DiF = -iF*DF(:,:,ii)*iF;
- DK(:,:,ii) = DP(:,Z,ii)*iF + P(:,Z)*DiF;
- Dtmp = DP(:,:,ii) - DK(:,:,ii)*P(Z,:) - K*DP(Z,:,ii);
+ DK(:,:,ii) = DP1(:,Z,ii)*iF + P(:,Z)*DiF;
+ Dtmp = DP1(:,:,ii) - DK(:,:,ii)*P(Z,:) - K*DP1(Z,:,ii);
DP1(:,:,ii) = DT(:,:,ii)*tmp*T' + T*Dtmp*T' + T*tmp*DT(:,:,ii)' + DOm(:,:,ii);
end
@@ -161,7 +165,7 @@ end
% end of computeDKalmanZ
-function [d2K,d2S,d2P1] = computeD2Kalman(A,dA,d2A,d2Om,P0,dP0,d2P0,DH,Z,iF,K0,dK0);
+function [d2K,d2S,d2P1] = computeD2Kalman(A,dA,d2A,d2Om,P0,dP0,d2P1,DH,Z,iF,K0,dK0);
% computes the second derivatives of the Kalman matrices
% note: A=T in main func.
@@ -176,8 +180,9 @@ function [d2K,d2S,d2P1] = computeD2Kalman(A,dA,d2A,d2Om,P0,dP0,d2P0,DH,Z,iF,K0,d
d2K = zeros(ns,no,k,k);
d2S = zeros(no,no,k,k);
-d2P1 = zeros(ns,ns,k,k);
+% d2P1 = zeros(ns,ns,k,k);
+jcount=0;
for ii = 1:k
dAi = dA(:,:,ii);
dFi = dP0(Z,Z,ii);
@@ -185,6 +190,7 @@ for ii = 1:k
diFi = -iF*dFi*iF;
dKi = dK0(:,:,ii);
for jj = 1:ii,
+ jcount=jcount+1;
dAj = dA(:,:,jj);
dFj = dP0(Z,Z,jj);
% d2Omj = d2Om(:,:,jj);
@@ -192,9 +198,9 @@ for ii = 1:k
diFj = -iF*dFj*iF;
dKj = dK0(:,:,jj);
- d2Aij = d2A(:,:,jj,ii);
- d2Pij = d2P0(:,:,jj,ii);
- d2Omij = d2Om(:,:,jj,ii);
+ d2Aij = reshape(d2A(:,jcount),[ns ns]);
+ d2Pij = dyn_unvech(d2P1(:,jcount));
+ d2Omij = dyn_unvech(d2Om(:,jcount));
% second order
@@ -216,10 +222,10 @@ for ii = 1:k
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<NASGU>
- d2P1(:,:,ii,jj) = d2AtmpA + d2Omij; %#ok<*NASGU>
+ d2P1(:,jcount) = dyn_vech(d2AtmpA + d2Omij); %#ok<*NASGU>
d2S(:,:,ii,jj) = d2Fij;
d2K(:,:,jj,ii) = d2Kij; %#ok<NASGU>
- d2P1(:,:,jj,ii) = d2AtmpA + d2Omij; %#ok<*NASGU>
+% d2P1(:,:,jj,ii) = d2AtmpA + d2Omij; %#ok<*NASGU>
d2S(:,:,jj,ii) = d2Fij;
% d2iS(:,:,ii,jj) = d2iF;
end
@@ -227,7 +233,7 @@ end
% end of computeD2Kalman
-function [d2K,d2S,d2P1] = computeD2KalmanZ(A,dA,d2A,d2Om,P0,dP0,d2P0,DH,Z,iF,K0,dK0);
+function [d2K,d2S,d2P1] = computeD2KalmanZ(A,dA,d2A,d2Om,P0,dP0,d2P1,DH,Z,iF,K0,dK0);
% computes the second derivatives of the Kalman matrices
% note: A=T in main func.
@@ -242,15 +248,17 @@ function [d2K,d2S,d2P1] = computeD2KalmanZ(A,dA,d2A,d2Om,P0,dP0,d2P0,DH,Z,iF,K0,
d2K = zeros(ns,no,k,k);
d2S = zeros(no,no,k,k);
-d2P1 = zeros(ns,ns,k,k);
+% d2P1 = zeros(ns,ns,k,k);
-for ii = 1:k
+jcount=0;
+for ii = 1:k,
dAi = dA(:,:,ii);
dFi = Z*dP0(:,:,ii)*Z;
% d2Omi = d2Om(:,:,ii);
diFi = -iF*dFi*iF;
dKi = dK0(:,:,ii);
for jj = 1:ii,
+ jcount=jcount+1;
dAj = dA(:,:,jj);
dFj = Z*dP0(:,:,jj)*Z;
% d2Omj = d2Om(:,:,jj);
@@ -258,9 +266,9 @@ for ii = 1:k
diFj = -iF*dFj*iF;
dKj = dK0(:,:,jj);
- d2Aij = d2A(:,:,jj,ii);
- d2Pij = d2P0(:,:,jj,ii);
- d2Omij = d2Om(:,:,jj,ii);
+ d2Aij = reshape(d2A(:,jcount),[ns ns]);
+ d2Pij = dyn_unvech(d2P1(:,jcount));
+ d2Omij = dyn_unvech(d2Om(:,jcount));
% second order
@@ -282,11 +290,11 @@ for ii = 1:k
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<NASGU>
- d2P1(:,:,ii,jj) = d2AtmpA + d2Omij; %#ok<*NASGU>
+ d2P1(:,jcount) = dyn_vech(d2AtmpA + d2Omij); %#ok<*NASGU>
d2S(:,:,ii,jj) = d2Fij;
% d2iS(:,:,ii,jj) = d2iF;
d2K(:,:,jj,ii) = d2Kij; %#ok<NASGU>
- d2P1(:,:,jj,ii) = d2AtmpA + d2Omij; %#ok<*NASGU>
+% d2P1(:,:,jj,ii) = d2AtmpA + d2Omij; %#ok<*NASGU>
d2S(:,:,jj,ii) = d2Fij;
end
end
diff --git a/matlab/kalman/likelihood/univariate_computeDLIK.m b/matlab/kalman/likelihood/univariate_computeDLIK.m
index a9eb34459152b0b65fce003f2b99636c24e827c6..1338a84a6284bb5a30a5b4ac958d3c7e500beb0d 100644
--- a/matlab/kalman/likelihood/univariate_computeDLIK.m
+++ b/matlab/kalman/likelihood/univariate_computeDLIK.m
@@ -34,12 +34,14 @@ if notsteady,
D2F = zeros(k,k);
D2v = zeros(k,k);
D2K = zeros(rows(K),k,k);
+ jcount=0;
for j=1:k,
D2v(:,j) = -Z*D2a(:,:,j) - Z*D2Yss(:,:,j);
- for i=j:k,
- D2F(j,i)=Z*D2P(:,:,j,i)*Z';
+ for i=1:j,
+ jcount=jcount+1;
+ D2F(j,i)=Z*dyn_unvech(D2P(:,jcount))*Z';
D2F(i,j)=D2F(j,i);
- D2K(:,j,i) = (D2P(:,:,j,i)*Z')/F-(DP(:,:,j)*Z')*DF(i)/F^2-(DP(:,:,i)*Z')*DF(j)/F^2- ...
+ D2K(:,j,i) = (dyn_unvech(D2P(:,jcount))*Z')/F-(DP(:,:,j)*Z')*DF(i)/F^2-(DP(:,:,i)*Z')*DF(j)/F^2- ...
PZ*D2F(j,i)/F^2 + 2*PZ/F^3*DF(i)*DF(j);
D2K(:,i,j) = D2K(:,j,i);
end
@@ -51,12 +53,22 @@ if notsteady,
DF = squeeze(DP(Z,Z,:))+DH';
DK = squeeze(DP(:,Z,:))/F-PZ*transpose(DF)/F^2;
if nargout>4
+ D2F = zeros(k,k);
+ D2K = zeros(rows(K),k,k);
D2v = squeeze(-D2a(Z,:,:) - D2Yss(Z,:,:));
- D2F = squeeze(D2P(Z,Z,:,:));
- D2K = squeeze(D2P(:,Z,:,:))/F;
+ jcount=0;
for j=1:k,
- D2K(:,:,j) = D2K(:,:,j) -PZ*D2F(j,:)/F^2 - squeeze(DP(:,Z,:))*DF(j)/F^2 - ...
- squeeze(DP(:,Z,j))*DF'/F^2 + 2/F^3*PZ*DF'*DF(j);
+ for i=1:j,
+ jcount=jcount+1;
+ tmp = dyn_unvech(D2P(:,jcount));
+ D2F(j,i) = tmp(Z,Z);
+ D2F(i,j)=D2F(j,i);
+ D2K(:,i,j) = tmp(:,Z)/F;
+ D2K(:,i,j) = D2K(:,i,j) -PZ*D2F(j,i)/F^2 - squeeze(DP(:,Z,i))*DF(j)/F^2 - ...
+ squeeze(DP(:,Z,j))*DF(i)'/F^2 + 2/F^3*PZ*DF(i)'*DF(j);
+ D2K(:,j,i) = D2K(:,i,j);
+% D2K = squeeze(D2P(:,Z,:,:))/F;
+ end
end
end
end
@@ -129,15 +141,26 @@ if notsteady,
if nargout>4,
if Zflag,
for j=1:k,
- D2P = D2P;
+ for i=1:j,
+ jcount = jcount+1;
+ tmp = dyn_unvech(D2P(:,jcount));
+ tmp = tmp - (tmp*Z')*K' - (DP(:,:,j)*Z')*DK(:,i)' ...
+ - (DP(:,:,i)*Z')*DK(:,j)' -PZ*D2K(:,j,i)';
+ D2P(:,jcount) = dyn_vech(tmp);
+% D2P(:,:,i,j) = D2P(:,:,j,i);
+ end
end
else
- D2PZ = squeeze(D2P(:,Z,:,:));
DPZ = squeeze(DP(:,Z,:));
+ jcount = 0;
for j=1:k,
- for i=j:k,
- D2P(:,:,j,i) = D2P(:,:,j,i) - D2PZ(:,j,i)*K' - DPZ(:,j)*DK(:,i)'- DPZ(:,i)*DK(:,j)' - PZ*squeeze(D2K(:,j,i))';
- D2P(:,:,i,j) = D2P(:,:,j,i);
+ for i=1:j,
+ jcount = jcount+1;
+ tmp = dyn_unvech(D2P(:,jcount));
+ D2PZ = tmp(:,Z);
+ tmp = tmp - D2PZ*K' - DPZ(:,j)*DK(:,i)'- DPZ(:,i)*DK(:,j)' - PZ*squeeze(D2K(:,j,i))';
+ D2P(:,jcount) = dyn_vech(tmp);
+% D2P(:,:,i,j) = D2P(:,:,j,i);
end
end
diff --git a/matlab/kalman/likelihood/univariate_computeDstate.m b/matlab/kalman/likelihood/univariate_computeDstate.m
index 30a44febbfff49063050f01a601fdd9182961bd1..a535a0a5189919f72000f688f6d11b5e0dac4703 100644
--- a/matlab/kalman/likelihood/univariate_computeDstate.m
+++ b/matlab/kalman/likelihood/univariate_computeDstate.m
@@ -34,21 +34,23 @@ if notsteady,
DP1 = DP1 + DOm;
end
if nargout>2,
+ jcount=0;
for j=1:k,
for i=1:j,
- D2a(:,j,i) = DT(:,:,i)*Da(:,j) + DT(:,:,j)*Da(:,i) + T*D2a(:,j,i)+ D2T(:,:,j,i)*a;
+ jcount=jcount+1;
+ D2a(:,j,i) = DT(:,:,i)*Da(:,j) + DT(:,:,j)*Da(:,i) + T*D2a(:,j,i)+ reshape(D2T(:,jcount),size(T))*a;
D2a(:,i,j) = D2a(:,j,i);
if notsteady,
- D2P(:,:,j,i) = T*D2P(:,:,j,i)*T' +DT(:,:,i)*DP(:,:,j)*T'+T*DP(:,:,j)*DT(:,:,i)' + ...
+ tmp = dyn_unvech(D2P(:,jcount));
+ tmp = T*tmp*T' +DT(:,:,i)*DP(:,:,j)*T'+T*DP(:,:,j)*DT(:,:,i)' + ...
DT(:,:,j)*DP(:,:,i)*T'+T*DP(:,:,i)*DT(:,:,j)' + ...
DT(:,:,j)*P*DT(:,:,i)'+DT(:,:,i)*P*DT(:,:,j)'+ ...
- D2T(:,:,j,i)*P*T'+T*P*D2T(:,:,j,i)';
- D2P(:,:,i,j) = D2P(:,:,j,i);
+ reshape(D2T(:,jcount),size(T))*P*T'+T*P*reshape(D2T(:,jcount),size(T))' + ...
+ dyn_unvech(D2Om(:,jcount));
+ D2P(:,jcount) = dyn_vech(tmp);
+% D2P(:,:,i,j) = D2P(:,:,j,i);
end
end
end
- if notsteady,
- D2P = D2P + D2Om;
- end
end
diff --git a/matlab/thet2tau.m b/matlab/thet2tau.m
index 79c593fad8b9e50bf8f6cc8da21da07e39cb57fd..2fedb0806a66b5823aef63628c8231706039ff0e 100644
--- a/matlab/thet2tau.m
+++ b/matlab/thet2tau.m
@@ -1,4 +1,4 @@
-function tau = thet2tau(params, M_, oo_, indx, indexo, flagmoments,mf,nlags,useautocorr)
+function tau = thet2tau(params, M_, oo_, indx, indexo, flagmoments,mf,nlags,useautocorr,iv)
%
% Copyright (C) 2011 Dynare Team
@@ -31,6 +31,9 @@ end
if nargin<9 || isempty(useautocorr),
useautocorr=0;
end
+if nargin<10 || isempty(iv),
+ iv=[1:M_.endo_nbr];
+end
M_.params(indx) = params(length(indexo)+1:end);
if ~isempty(indexo)
@@ -38,7 +41,8 @@ if ~isempty(indexo)
end
[A,B,tele,tubbies,M_,options_,oo_] = dynare_resolve(M_,options_,oo_);
if flagmoments==0,
- tau = [oo_.dr.ys(oo_.dr.order_var); A(:); dyn_vech(B*M_.Sigma_e*B')];
+ ys=oo_.dr.ys(oo_.dr.order_var);
+ tau = [ys(iv); vec(A(iv,iv)); dyn_vech(B(iv,:)*M_.Sigma_e*B(iv,:)')];
elseif flagmoments==-1
[I,J]=find(M_.lead_lag_incidence');
yy0=oo_.dr.ys(I);