Description of getH

0001 function [H, dA, dOm, Hss, gp, d2A, d2Om, H2ss] = getH(A, B, M_,oo_,kronflag,indx,indexo)
0002 
0003 % computes derivative of reduced form linear model w.r.t. deep params
0004 %
0005 % Copyright (C) 2010-2011 Dynare Team
0006 %
0007 % This file is part of Dynare.
0008 %
0009 % Dynare is free software: you can redistribute it and/or modify
0010 % it under the terms of the GNU General Public License as published by
0011 % the Free Software Foundation, either version 3 of the License, or
0012 % (at your option) any later version.
0013 %
0014 % Dynare is distributed in the hope that it will be useful,
0015 % but WITHOUT ANY WARRANTY; without even the implied warranty of
0016 % MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
0017 % GNU General Public License for more details.
0018 %
0019 % You should have received a copy of the GNU General Public License
0020 % along with Dynare.  If not, see <http://www.gnu.org/licenses/>.
0021 
0022 if nargin<3 || isempty(kronflag), kronflag = 0; end
0023 if nargin<4 || isempty(indx), indx = [1:M_.param_nbr];, end,
0024 if nargin<5 || isempty(indexo), indexo = [];, end,
0025 
0026 
0027 [I,J]=find(M_.lead_lag_incidence');
0028 yy0=oo_.dr.ys(I);
0029 % yy0=[];
0030 % for j=1:size(M_.lead_lag_incidence,1);
0031 %     yy0 = [ yy0; oo_.dr.ys(find(M_.lead_lag_incidence(j,:)))];
0032 % end
0033 dyssdtheta=zeros(length(oo_.dr.ys),M_.param_nbr);
0034 if nargout>5,
0035     [residual, gg1, gg2] = feval([M_.fname,'_static'],oo_.dr.ys, oo_.exo_steady_state', M_.params);
0036     [df, gp, d2f] = feval([M_.fname,'_params_derivs'],yy0, oo_.exo_steady_state', ...
0037         M_.params, oo_.dr.ys, 1, dyssdtheta);
0038     d2f = get_all_resid_2nd_derivs(d2f,length(oo_.dr.ys),M_.param_nbr);
0039     d2f = d2f(:,indx,indx);
0040     if isempty(find(gg2)),
0041         for j=1:length(indx),
0042         d2yssdtheta(:,:,j) = -gg1\d2f(:,:,j);
0043         end
0044     else
0045         gam = d2f*0;
0046         for j=1:length(indx),
0047         d2yssdtheta(:,:,j) = -gg1\(d2f(:,:,j)+gam(:,:,j));
0048         end
0049     end
0050 else
0051     [residual, gg1] = feval([M_.fname,'_static'],oo_.dr.ys, oo_.exo_steady_state', M_.params);
0052     df = feval([M_.fname,'_params_derivs'],yy0, oo_.exo_steady_state', ...
0053         M_.params, oo_.dr.ys, 1, dyssdtheta);
0054 end
0055 dyssdtheta = -gg1\df;
0056 
0057 if any(any(isnan(dyssdtheta))),    
0058     [U,T] = schur(gg1);
0059     global options_
0060     qz_criterium=options_.qz_criterium;
0061     e1 = abs(ordeig(T)) < qz_criterium-1;
0062     k = sum(e1);       % Number non stationary variables.
0063     n = length(e1)-k;  % Number of stationary variables.
0064     [U,T] = ordschur(U,T,e1);
0065     T = T(k+1:end,k+1:end);
0066     dyssdtheta = -U(:,k+1:end)*(T\U(:,k+1:end)')*df;
0067     if nargout>5,
0068         for j=1:length(indx),
0069             d2yssdtheta(:,:,j) = -U(:,k+1:end)*(T\U(:,k+1:end)')*d2f(:,:,j);
0070         end
0071     end
0072 end
0073 if nargout>5,
0074 [df, gp, d2f, gpp] = feval([M_.fname,'_params_derivs'],yy0, oo_.exo_steady_state', ...
0075            M_.params, oo_.dr.ys, 1, dyssdtheta);
0076 H2ss = d2yssdtheta(oo_.dr.order_var,:,:);
0077 [residual, g1, g2, g3] = feval([M_.fname,'_dynamic'],yy0, oo_.exo_steady_state', ...
0078                             M_.params, oo_.dr.ys, 1);
0079 else
0080 [df, gp] = feval([M_.fname,'_params_derivs'],yy0, oo_.exo_steady_state', ...
0081            M_.params, oo_.dr.ys, 1, dyssdtheta);
0082 [residual, g1, g2 ] = feval([M_.fname,'_dynamic'],yy0, oo_.exo_steady_state', ...
0083                             M_.params, oo_.dr.ys, 1);
0084 end
0085 
0086 Hss = dyssdtheta(oo_.dr.order_var,indx);
0087 dyssdtheta = dyssdtheta(I,:);
0088 [nr, nc]=size(g2);
0089 [m, nelem]=size(g1);
0090 nc = sqrt(nc);
0091 ns = max(max(M_.lead_lag_incidence));
0092 gp2 = gp*0;
0093 for j=1:nr,
0094     [I J]=ind2sub([nc nc],find(g2(j,:)));
0095     for i=1:nc,
0096         is = find(I==i);
0097         is = is(find(J(is)<=ns));
0098         if ~isempty(is),
0099             g20=full(g2(j,find(g2(j,:))));
0100             gp2(j,i,:)=g20(is)*dyssdtheta(J(is),:);
0101         end
0102     end
0103 end
0104 
0105 
0106 gp = gp+gp2;
0107 gp = gp(:,:,indx);
0108 param_nbr = length(indx);
0109 if nargout>5,
0110     param_nbr_2 = param_nbr*(param_nbr+1)/2;
0111 end
0112 
0113 m = size(A,1);
0114 n = size(B,2);
0115 
0116 
0117 
0118 klen = M_.maximum_endo_lag + M_.maximum_endo_lead + 1;
0119 k11 = M_.lead_lag_incidence(find([1:klen] ~= M_.maximum_endo_lag+1),:);
0120 a = g1(:,nonzeros(k11'));
0121 da = gp(:,nonzeros(k11'),:);
0122 if nargout > 5,
0123     nelem = size(g1,2);
0124     g22 = get_all_2nd_derivs(gpp,m,nelem,M_.param_nbr);
0125     g22 = g22(:,:,indx,indx);
0126     d2a = g22(:,nonzeros(k11'),:,:);
0127 end
0128 kstate = oo_.dr.kstate;
0129 
0130 GAM1 = zeros(M_.endo_nbr,M_.endo_nbr);
0131 Dg1 = zeros(M_.endo_nbr,M_.endo_nbr,param_nbr);
0132 % nf = find(M_.lead_lag_incidence(M_.maximum_endo_lag+2,:));
0133 % GAM1(:, nf) = -g1(:,M_.lead_lag_incidence(M_.maximum_endo_lag+2,nf));
0134 
0135 k1 = find(kstate(:,2) == M_.maximum_endo_lag+2 & kstate(:,3));
0136 GAM1(:, kstate(k1,1)) = -a(:,kstate(k1,3));
0137 Dg1(:, kstate(k1,1), :) = -da(:,kstate(k1,3),:);
0138 if nargout > 5,
0139     D2g1 = zeros(M_.endo_nbr,M_.endo_nbr,param_nbr,param_nbr);
0140     D2g1(:, kstate(k1,1), :, :) = -d2a(:,kstate(k1,3),:,:);
0141 end
0142 % k = find(kstate(:,2) > M_.maximum_endo_lag+2 & kstate(:,3));
0143 % kk = find(kstate(:,2) > M_.maximum_endo_lag+2 & ~kstate(:,3));
0144 % nauxe = 0;
0145 % if ~isempty(k),
0146 %     ax(:,k)= -a(:,kstate(k,3));
0147 %     ax(:,kk)= 0;
0148 %     dax(:,k,:)= -da(:,kstate(k,3),:);
0149 %     dax(:,kk,:)= 0;
0150 %     nauxe=size(ax,2);
0151 %     GAM1 = [ax GAM1];
0152 %     Dg1 = cat(2,dax,Dg1);
0153 %     clear ax
0154 % end
0155 
0156 
0157 [junk,cols_b,cols_j] = find(M_.lead_lag_incidence(M_.maximum_endo_lag+1, ...
0158                                                   oo_.dr.order_var));
0159 GAM0 = zeros(M_.endo_nbr,M_.endo_nbr);
0160 Dg0 = zeros(M_.endo_nbr,M_.endo_nbr,param_nbr);
0161 GAM0(:,cols_b) = g1(:,cols_j);
0162 Dg0(:,cols_b,:) = gp(:,cols_j,:);
0163 if nargout > 5,
0164     D2g0 = zeros(M_.endo_nbr,M_.endo_nbr,param_nbr,param_nbr);
0165     D2g0(:, cols_b, :, :) = g22(:,cols_j,:,:);
0166 end
0167 % GAM0 = g1(:,M_.lead_lag_incidence(M_.maximum_endo_lag+1,:));
0168 
0169 
0170 k2 = find(kstate(:,2) == M_.maximum_endo_lag+1 & kstate(:,4));
0171 GAM2 = zeros(M_.endo_nbr,M_.endo_nbr);
0172 Dg2 = zeros(M_.endo_nbr,M_.endo_nbr,param_nbr);
0173 % nb = find(M_.lead_lag_incidence(1,:));
0174 % GAM2(:, nb) = -g1(:,M_.lead_lag_incidence(1,nb));
0175 GAM2(:, kstate(k2,1)) = -a(:,kstate(k2,4));
0176 Dg2(:, kstate(k2,1), :) = -da(:,kstate(k2,4),:);
0177 if nargout > 5,
0178     D2g2 = zeros(M_.endo_nbr,M_.endo_nbr,param_nbr,param_nbr);
0179     D2g2(:, kstate(k2,1), :, :) = -d2a(:,kstate(k2,4),:,:);
0180 end% naux = 0;
0181 % k = find(kstate(:,2) < M_.maximum_endo_lag+1 & kstate(:,4));
0182 % kk = find(kstate(:,2) < M_.maximum_endo_lag+1 );
0183 % if ~isempty(k),
0184 %     ax(:,k)= -a(:,kstate(k,4));
0185 %     ax = ax(:,kk);
0186 %     dax(:,k,:)= -da(:,kstate(k,4),:);
0187 %     dax = dax(:,kk,:);
0188 %     naux = size(ax,2);
0189 %     GAM2 = [GAM2 ax];
0190 %     Dg2 = cat(2,Dg2,dax);
0191 % end
0192 %
0193 % GAM0 = blkdiag(GAM0,eye(naux));
0194 % Dg0 = cat(1,Dg0,zeros(naux+nauxe,M_.endo_nbr,param_nbr));
0195 % Dg0 = cat(2,Dg0,zeros(m+nauxe,naux,param_nbr));
0196 % Dg0 = cat(2,zeros(m+nauxe,nauxe,param_nbr),Dg0);
0197 %
0198 % GAM2 = [GAM2 ; A(M_.endo_nbr+1:end,:)];
0199 % Dg2 = cat(1,Dg2,zeros(naux+nauxe,m,param_nbr));
0200 % Dg2 = cat(2,zeros(m+nauxe,nauxe,param_nbr),Dg2);
0201 % GAM2 = [zeros(m+nauxe,nauxe) [GAM2; zeros(nauxe,m)]];
0202 % GAM0 = [[zeros(m,nauxe);(eye(nauxe))] [GAM0; zeros(nauxe,m)]];
0203 
0204 GAM3 = -g1(:,length(yy0)+1:end);
0205 % GAM3 = -g1(oo_.dr.order_var,length(yy0)+1:end);
0206 % GAM3 = [GAM3; zeros(naux+nauxe,size(GAM3,2))];
0207 % Dg3 = -gp(oo_.dr.order_var,length(yy0)+1:end,:);
0208 Dg3 = -gp(:,length(yy0)+1:end,:);
0209 if nargout>5,
0210     D2g3 = -g22(:,length(yy0)+1:end,:,:);
0211     clear g22 d2a
0212 end
0213 % Dg3 = cat(1,Dg3,zeros(naux+nauxe,size(GAM3,2),param_nbr));
0214 
0215 % auxe = zeros(0,1);
0216 % k0 = kstate(find(kstate(:,2) >= M_.maximum_endo_lag+2),:);;
0217 % i0 = find(k0(:,2) == M_.maximum_endo_lag+2);
0218 % for i=M_.maximum_endo_lag+3:M_.maximum_endo_lag+2+M_.maximum_endo_lead,
0219 %     i1 = find(k0(:,2) == i);
0220 %     n1 = size(i1,1);
0221 %     j = zeros(n1,1);
0222 %     for j1 = 1:n1
0223 %         j(j1) = find(k0(i0,1)==k0(i1(j1),1));
0224 %     end
0225 %     auxe = [auxe; i0(j)];
0226 %     i0 = i1;
0227 % end
0228 % auxe = [(1:size(auxe,1))' auxe(end:-1:1)];
0229 % n_ir1 = size(auxe,1);
0230 % nr = m + n_ir1;
0231 % GAM1 = [[GAM1 zeros(size(GAM1,1),naux)]; zeros(naux+nauxe,m+nauxe)];
0232 % GAM1(m+1:end,:) = sparse(auxe(:,1),auxe(:,2),ones(n_ir1,1),n_ir1,nr);
0233 % Dg1 = cat(2,Dg1,zeros(M_.endo_nbr,naux,param_nbr));
0234 % Dg1 = cat(1,Dg1,zeros(naux+nauxe,m+nauxe,param_nbr));
0235 
0236 % Ax = A;
0237 % A1 = A;
0238 % Bx = B;
0239 % B1 = B;
0240 % for j=1:M_.maximum_endo_lead-1,
0241 %     A1 = A1*A;
0242 %     B1 = A*B1;
0243 %     k = find(kstate(:,2) == M_.maximum_endo_lag+2+j );
0244 %     Ax = [A1(kstate(k,1),:); Ax];
0245 %     Bx = [B1(kstate(k,1),:); Bx];
0246 % end
0247 % Ax = [zeros(m+nauxe,nauxe) Ax];
0248 % A0 = A;
0249 % A=Ax; clear Ax A1;
0250 % B0=B;
0251 % B = Bx; clear Bx B1;
0252 
0253 m = size(A,1);
0254 n = size(B,2);
0255 
0256 % Dg1 = zeros(m,m,param_nbr);
0257 % Dg1(:, nf, :) = -gp(:,M_.lead_lag_incidence(3,nf),:);
0258 
0259 % Dg0 = gp(:,M_.lead_lag_incidence(2,:),:);
0260 
0261 % Dg2 = zeros(m,m,param_nbr);
0262 % Dg2(:, nb, :) = -gp(:,M_.lead_lag_incidence(1,nb),:);
0263 
0264 % Dg3 = -gp(:,length(yy0)+1:end,:);
0265 
0266 if kronflag==1, % kronecker products
0267     Dg0=reshape(Dg0,m^2,param_nbr);
0268     Dg1=reshape(Dg1,m^2,param_nbr);
0269     Dg2=reshape(Dg2,m^2,param_nbr);
0270     for j=1:param_nbr,
0271         Dg3(:,:,j)=Dg3(:,:,j)*M_.Sigma_e;
0272     end
0273     Dg3=reshape(Dg3,m*n,param_nbr);
0274     Om = B*M_.Sigma_e*B';
0275     Im = eye(m);
0276     Dm = duplication(m);
0277     DmPl = inv(Dm'*Dm)*Dm';
0278     Kmm = commutation(m,m);
0279     Kmn = commutation(m,n);
0280 
0281 
0282     Da = [eye(m^2),zeros(m^2,m*(m+1)/2)];
0283     Dom = [zeros(m*(m+1)/2,m^2),eye(m*(m+1)/2)];
0284 
0285 
0286     Df1Dtau = ( kron(Im,GAM0) - kron(A',GAM1) - kron(Im,GAM1*A) )*Da;
0287 
0288     Df1Dthet = kron(A',Im)*Dg0 - kron( (A')^2,Im)*Dg1 - Dg2;
0289 
0290     Df2Dtau = DmPl*( kron(GAM0,GAM0) - kron(GAM0,GAM1*A) - kron(GAM1*A,GAM0) + kron(GAM1*A,GAM1*A) )*Dm*Dom - ...
0291               DmPl*( kron(GAM0*Om,GAM1) + kron(GAM1,GAM0*Om)*Kmm - kron(GAM1*A*Om,GAM1) - kron(GAM1,GAM1*A*Om)*Kmm )*Da;
0292 
0293 
0294     Df2Dthet = DmPl*( kron(GAM0*Om,Im) + kron(Im,GAM0*Om)*Kmm - kron(Im,GAM1*A*Om)*Kmm - kron(GAM1*A*Om,Im) )*Dg0 - ...
0295         DmPl*( kron(GAM0*Om*A',Im) + kron(Im,GAM0*Om*A')*Kmm - kron(Im,GAM1*A*Om*A')*Kmm - kron(GAM1*A*Om*A',Im) )*Dg1 -...
0296         DmPl*( kron(GAM3,Im) + kron(Im,GAM3)*Kmn )*Dg3;
0297 
0298 
0299     DfDtau  = [Df1Dtau;Df2Dtau];
0300     DfDthet = [Df1Dthet;Df2Dthet];
0301 
0302     H = -DfDtau\DfDthet;
0303     x = reshape(H(1:m*m,:),m,m,param_nbr);
0304     y = reshape(Dm*H(m*m+1:end,:),m,m,param_nbr);
0305     if nauxe,
0306         x = x(nauxe+1:end,nauxe+1:end,:);
0307         y = y(nauxe+1:end,nauxe+1:end,:);
0308         dA = x;
0309         dOm = y;
0310         m = size(y,1);
0311         x = reshape(x,m*m,param_nbr);
0312         Dm = duplication(m);
0313         DmPl = inv(Dm'*Dm)*Dm';
0314         y = DmPl*reshape(y,m*m,param_nbr);
0315         H = [x;y];
0316     else
0317         dA = x;
0318         dOm = y;
0319     end
0320 
0321     Hx = [];
0322     if ~isempty(indexo),
0323         dSig = zeros(M_.exo_nbr,M_.exo_nbr);
0324         dOm = cat(3,zeros(size(dOm,1),size(dOm,1),length(indexo)),dOm);
0325         for j=1:length(indexo)
0326             dSig(indexo(j),indexo(j)) = 2*sqrt(M_.Sigma_e(indexo(j),indexo(j)));
0327             y = B*dSig*B';
0328             y = y(nauxe+1:end,nauxe+1:end);
0329             Hx(:,j) = [zeros((m-nauxe)^2,1); dyn_vech(y)];
0330             if nargout>1,
0331                 dOm(:,:,j) = y;
0332             end
0333             dSig(indexo(j),indexo(j)) = 0;
0334         end
0335     end
0336     H = [ [zeros(M_.endo_nbr,length(indexo)) Hss]; [Hx H]];
0337 
0338 elseif kronflag==-1, % perturbation
0339     fun = 'thet2tau';
0340     params0 = M_.params;
0341     H = fjaco(fun,[sqrt(diag(M_.Sigma_e(indexo,indexo))); M_.params(indx)],indx,indexo);
0342     assignin('base','M_', M_);
0343     assignin('base','oo_', oo_);
0344 
0345 else % generalized sylvester equation
0346 
0347     % solves a*x+b*x*c=d
0348     a = (GAM0-GAM1*A);
0349     inva = inv(a);
0350     b = -GAM1;
0351     c = A;
0352     elem = zeros(m,m,param_nbr);
0353     d = elem;
0354     for j=1:param_nbr,
0355         elem(:,:,j) = (Dg0(:,:,j)-Dg1(:,:,j)*A);
0356         d(:,:,j) = Dg2(:,:,j)-elem(:,:,j)*A;
0357     end
0358     xx=sylvester3(a,b,c,d);
0359     flag=1;
0360     icount=0;
0361     while flag && icount<4,
0362         [xx, flag]=sylvester3a(xx,a,b,c,d);
0363         icount=icount+1;
0364     end
0365     H=zeros(m*m+m*(m+1)/2,param_nbr+length(indexo));
0366     if nargout>1,
0367         dOm = zeros(m,m,param_nbr+length(indexo));
0368         dA=zeros(m,m,param_nbr+length(indexo));
0369         dB=zeros(m,n,param_nbr);
0370     end
0371     if ~isempty(indexo),
0372         dSig = zeros(M_.exo_nbr,M_.exo_nbr,length(indexo));
0373         for j=1:length(indexo)
0374             dSig(indexo(j),indexo(j),j) = 2*sqrt(M_.Sigma_e(indexo(j),indexo(j)));
0375             y = B*dSig(:,:,j)*B';
0376 %             y = y(nauxe+1:end,nauxe+1:end);
0377 %             H(:,j) = [zeros((m-nauxe)^2,1); dyn_vech(y)];
0378             H(:,j) = [zeros(m^2,1); dyn_vech(y)];
0379             if nargout>1,
0380                 dOm(:,:,j) = y;
0381             end
0382 %             dSig(indexo(j),indexo(j)) = 0;
0383         end
0384     end
0385     for j=1:param_nbr,
0386         x = xx(:,:,j);
0387         y = inva * (Dg3(:,:,j)-(elem(:,:,j)-GAM1*x)*B);
0388         if nargout>1,
0389             dB(:,:,j) = y;
0390         end
0391         y = y*M_.Sigma_e*B'+B*M_.Sigma_e*y';
0392 %         x = x(nauxe+1:end,nauxe+1:end);
0393 %         y = y(nauxe+1:end,nauxe+1:end);
0394         if nargout>1,
0395             dA(:,:,j+length(indexo)) = x;
0396             dOm(:,:,j+length(indexo)) = y;
0397         end
0398         H(:,j+length(indexo)) = [x(:); dyn_vech(y)];
0399     end
0400     %     for j=1:param_nbr,
0401     %         disp(['Derivatives w.r.t. ',M_.param_names(indx(j),:),', ',int2str(j),'/',int2str(param_nbr)])
0402     %         elem = (Dg0(:,:,j)-Dg1(:,:,j)*A);
0403     %         d = Dg2(:,:,j)-elem*A;
0404     %         x=sylvester3(a,b,c,d);
0405     % %         x=sylvester3a(x,a,b,c,d);
0406     %         y = inva * (Dg3(:,:,j)-(elem-GAM1*x)*B);
0407     %         y = y*B'+B*y';
0408     %         x = x(nauxe+1:end,nauxe+1:end);
0409     %         y = y(nauxe+1:end,nauxe+1:end);
0410     %         H(:,j) = [x(:); dyn_vech(y)];
0411     %     end
0412     H = [[zeros(M_.endo_nbr,length(indexo)) Hss]; H];
0413 
0414 end
0415 
0416 if nargout > 5,
0417     tot_param_nbr = param_nbr + length(indexo);
0418     tot_param_nbr_2 = tot_param_nbr*(tot_param_nbr+1)/2;
0419     d = zeros(m,m,param_nbr_2);
0420     d2A = zeros(m,m,tot_param_nbr,tot_param_nbr);
0421     d2Om = zeros(m,m,tot_param_nbr,tot_param_nbr);
0422     d2B = zeros(m,n,tot_param_nbr,tot_param_nbr);
0423     cc=triu(ones(param_nbr,param_nbr));
0424     [i2,j2]=find(cc);
0425     cc = blkdiag( zeros(length(indexo),length(indexo)), cc);
0426     [ipar2]=find(cc);
0427     ctot=triu(ones(tot_param_nbr,tot_param_nbr));
0428     ctot(1:length(indexo),1:length(indexo))=eye(length(indexo));
0429     [itot2, jtot2]=find(ctot);
0430     jcount=0;
0431     for j=1:param_nbr,
0432         for i=j:param_nbr,
0433         elem1 = (D2g0(:,:,j,i)-D2g1(:,:,j,i)*A);
0434         elem1 = D2g2(:,:,j,i)-elem1*A;
0435         elemj0 = Dg0(:,:,j)-Dg1(:,:,j)*A;
0436         elemi0 = Dg0(:,:,i)-Dg1(:,:,i)*A;
0437         elem2 = -elemj0*xx(:,:,i)-elemi0*xx(:,:,j);
0438         elem2 = elem2 + ( Dg1(:,:,j)*xx(:,:,i) + Dg1(:,:,i)*xx(:,:,j) )*A;
0439         elem2 = elem2 + GAM1*( xx(:,:,i)*xx(:,:,j) + xx(:,:,j)*xx(:,:,i));
0440         jcount=jcount+1;
0441         d(:,:,jcount) = elem1+elem2;
0442         end
0443     end
0444     xx2=sylvester3(a,b,c,d);
0445     flag=1;
0446     icount=0;
0447     while flag && icount<4,
0448         [xx2, flag]=sylvester3a(xx2,a,b,c,d);
0449         icount = icount + 1;
0450     end
0451     jcount = 0;
0452     for j=1:param_nbr,
0453         for i=j:param_nbr,
0454         jcount=jcount+1;
0455         x = xx2(:,:,jcount);
0456         elem1 = (D2g0(:,:,j,i)-D2g1(:,:,j,i)*A);
0457         elem1 = elem1 -( Dg1(:,:,j)*xx(:,:,i) + Dg1(:,:,i)*xx(:,:,j) );
0458         elemj0 = Dg0(:,:,j)-Dg1(:,:,j)*A-GAM1*xx(:,:,j);
0459         elemi0 = Dg0(:,:,i)-Dg1(:,:,i)*A-GAM1*xx(:,:,i);
0460         elem0 = elemj0*dB(:,:,i)+elemi0*dB(:,:,j);
0461         y = inva * (D2g3(:,:,j,i)-elem0-(elem1-GAM1*x)*B);
0462         d2B(:,:,j+length(indexo),i+length(indexo)) = y;
0463         d2B(:,:,i+length(indexo),j+length(indexo)) = y;
0464         y = y*M_.Sigma_e*B'+B*M_.Sigma_e*y'+ ...
0465             dB(:,:,j)*M_.Sigma_e*dB(:,:,i)'+dB(:,:,i)*M_.Sigma_e*dB(:,:,j)';
0466 %         x = x(nauxe+1:end,nauxe+1:end);
0467 %         y = y(nauxe+1:end,nauxe+1:end);
0468         d2A(:,:,j+length(indexo),i+length(indexo)) = x;
0469         d2A(:,:,i+length(indexo),j+length(indexo)) = x;
0470         d2Om(:,:,j+length(indexo),i+length(indexo)) = y;
0471         d2Om(:,:,i+length(indexo),j+length(indexo)) = y;
0472         end
0473     end    
0474     if ~isempty(indexo),
0475         d2Sig = zeros(M_.exo_nbr,M_.exo_nbr,length(indexo));
0476         for j=1:length(indexo)
0477             d2Sig(indexo(j),indexo(j),j) = 2;
0478             y = B*d2Sig(:,:,j)*B';
0479             d2Om(:,:,j,j) = y;
0480 %             y = y(nauxe+1:end,nauxe+1:end);
0481 %             H(:,j) = [zeros((m-nauxe)^2,1); dyn_vech(y)];
0482 %             H(:,j) = [zeros(m^2,1); dyn_vech(y)];
0483 %             dOm(:,:,j) = y;
0484             for i=1:param_nbr, 
0485                 y = dB(:,:,i)*dSig(:,:,j)*B'+B*dSig(:,:,j)*dB(:,:,i)';
0486                 d2Om(:,:,j,i+length(indexo)) = y;
0487                 d2Om(:,:,i+length(indexo),j) = y;
0488             end
0489         end
0490     end
0491 end
0492 
0493 return
0494 
0495 function g22 = get_2nd_deriv(gpp,m,n,i,j),
0496 
0497 g22=zeros(m,n);
0498 is=find(gpp(:,3)==i);
0499 is=is(find(gpp(is,4)==j));
0500 
0501 if ~isempty(is),
0502     g22(gpp(is,1),gpp(is,2))=gpp(is,5);
0503 end
0504 return
0505 
0506 function g22 = get_all_2nd_derivs(gpp,m,n,npar),
0507 
0508 g22=zeros(m,n,npar,npar);
0509 % c=ones(npar,npar);
0510 % c=triu(c);
0511 % ic=find(c);
0512 
0513 for is=1:length(gpp),
0514 %     d=zeros(npar,npar);
0515 %     d(gpp(is,3),gpp(is,4))=1;
0516 %     indx = find(ic==find(d));
0517     g22(gpp(is,1),gpp(is,2),gpp(is,3),gpp(is,4))=gpp(is,5);
0518     g22(gpp(is,1),gpp(is,2),gpp(is,4),gpp(is,3))=gpp(is,5);
0519 end
0520 
0521 return
0522 
0523 function r22 = get_all_resid_2nd_derivs(rpp,m,npar),
0524 
0525 r22=zeros(m,npar,npar);
0526 % c=ones(npar,npar);
0527 % c=triu(c);
0528 % ic=find(c);
0529 
0530 for is=1:length(rpp),
0531 %     d=zeros(npar,npar);
0532 %     d(rpp(is,2),rpp(is,3))=1;
0533 %     indx = find(ic==find(d));
0534     r22(rpp(is,1),rpp(is,2),rpp(is,3))=rpp(is,4);
0535     r22(rpp(is,1),rpp(is,3),rpp(is,2))=rpp(is,4);
0536 end
0537 
0538 return
getH

PURPOSE

SYNOPSIS

DESCRIPTION

CROSS-REFERENCE INFORMATION

SUBFUNCTIONS

SOURCE CODE
