getJJ

 Copyright (C) 2010-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 <http://www.gnu.org/licenses/>.

0001 function [JJ, H, gam, gp, dA, dOm, dYss] = getJJ(A, B, M_,oo_,options_,kronflag,indx,indexo,mf,nlags,useautocorr)
0002 
0003 % Copyright (C) 2010-2011 Dynare Team
0004 %
0005 % This file is part of Dynare.
0006 %
0007 % Dynare is free software: you can redistribute it and/or modify
0008 % it under the terms of the GNU General Public License as published by
0009 % the Free Software Foundation, either version 3 of the License, or
0010 % (at your option) any later version.
0011 %
0012 % Dynare is distributed in the hope that it will be useful,
0013 % but WITHOUT ANY WARRANTY; without even the implied warranty of
0014 % MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
0015 % GNU General Public License for more details.
0016 %
0017 % You should have received a copy of the GNU General Public License
0018 % along with Dynare.  If not, see <http://www.gnu.org/licenses/>.
0019 
0020 if nargin<7 || isempty(indx), indx = [1:M_.param_nbr];, end,
0021 if nargin<8 || isempty(indexo), indexo = [];, end,
0022 if nargin<10 || isempty(nlags), nlags=3; end,
0023 if nargin<11 || isempty(useautocorr), useautocorr=0; end,
0024 
0025 %   if useautocorr,
0026 warning('off','MATLAB:divideByZero')
0027 %   end
0028 if kronflag == -1,
0029     fun = 'thet2tau';
0030     params0 = M_.params;
0031     JJ = fjaco(fun,[sqrt(diag(M_.Sigma_e(indexo,indexo))); M_.params(indx)],indx,indexo,1,mf,nlags,useautocorr);
0032     M_.params = params0;
0033     params0 = M_.params;
0034     H = fjaco(fun,[sqrt(diag(M_.Sigma_e(indexo,indexo))); M_.params(indx)],indx,indexo,0,mf,nlags,useautocorr);
0035     M_.params = params0;
0036     assignin('base','M_', M_);
0037     assignin('base','oo_', oo_);
0038 else
0039     [H, dA, dOm, dYss, gp] = getH(A, B, M_,oo_,kronflag,indx,indexo);
0040     gp = reshape(gp,size(gp,1)*size(gp,2),size(gp,3));
0041     gp = [dYss; gp];
0042     %   if isempty(H),
0043     %     JJ = [];
0044     %     GAM = [];
0045     %     return
0046     %   end
0047     m = length(A);
0048     GAM =  lyapunov_symm(A,B*M_.Sigma_e*B',options_.qz_criterium,options_.lyapunov_complex_threshold,1);
0049     k = find(abs(GAM) < 1e-12);
0050     GAM(k) = 0;
0051     %   if useautocorr,
0052     sdy = sqrt(diag(GAM));
0053     sy = sdy*sdy';
0054     %   end
0055     
0056     %   BB = dOm*0;
0057     %   for j=1:length(indx),
0058     %     BB(:,:,j)= dA(:,:,j)*GAM*A'+A*GAM*dA(:,:,j)'+dOm(:,:,j);
0059     %   end
0060     %   XX =  lyapunov_symm_mr(A,BB,options_.qz_criterium,options_.lyapunov_complex_threshold,0);
0061     for j=1:length(indexo),
0062         dum =  lyapunov_symm(A,dOm(:,:,j),options_.qz_criterium,options_.lyapunov_complex_threshold,2);
0063         %     dum =  XX(:,:,j);
0064         k = find(abs(dum) < 1e-12);
0065         dum(k) = 0;
0066         if useautocorr
0067             dsy = 1/2./sdy.*diag(dum);
0068             dsy = dsy*sdy'+sdy*dsy';
0069             dum1=dum;
0070             dum1 = (dum1.*sy-dsy.*GAM)./(sy.*sy);
0071             dum1 = dum1-diag(diag(dum1))+diag(diag(dum));
0072             dumm = dyn_vech(dum1(mf,mf));
0073         else
0074             dumm = dyn_vech(dum(mf,mf));
0075         end
0076         for i=1:nlags,
0077             dum1 = A^i*dum;
0078             if useautocorr
0079                 dum1 = (dum1.*sy-dsy.*(A^i*GAM))./(sy.*sy);
0080             end
0081             dumm = [dumm; vec(dum1(mf,mf))];
0082         end
0083         JJ(:,j) = dumm;
0084     end
0085     nexo = length(indexo);
0086     for j=1:length(indx),
0087         dum =  lyapunov_symm(A,dA(:,:,j+nexo)*GAM*A'+A*GAM*dA(:,:,j+nexo)'+dOm(:,:,j+nexo),options_.qz_criterium,options_.lyapunov_complex_threshold,2);
0088         %     dum =  XX(:,:,j);
0089         k = find(abs(dum) < 1e-12);
0090         dum(k) = 0;
0091         if useautocorr
0092             dsy = 1/2./sdy.*diag(dum);
0093             dsy = dsy*sdy'+sdy*dsy';
0094             dum1=dum;
0095             dum1 = (dum1.*sy-dsy.*GAM)./(sy.*sy);
0096             dum1 = dum1-diag(diag(dum1))+diag(diag(dum));
0097             dumm = dyn_vech(dum1(mf,mf));
0098         else
0099             dumm = dyn_vech(dum(mf,mf));
0100         end
0101         for i=1:nlags,
0102             dum1 = A^i*dum;
0103             for ii=1:i,
0104                 dum1 = dum1 + A^(ii-1)*dA(:,:,j+nexo)*A^(i-ii)*GAM;
0105             end
0106             if useautocorr
0107                 dum1 = (dum1.*sy-dsy.*(A^i*GAM))./(sy.*sy);
0108             end
0109             dumm = [dumm; vec(dum1(mf,mf))];
0110         end
0111         JJ(:,j+nexo) = dumm;
0112     end
0113     
0114     JJ = [ [zeros(length(mf),nexo) dYss(mf,:)]; JJ];
0115     
0116 end
0117 if nargout >2,
0118     %     sy=sy(mf,mf);
0119     options_.ar=nlags;
0120     [GAM,stationary_vars] = th_autocovariances(oo_.dr,oo_.dr.order_var(mf),M_,options_);
0121     sy=sqrt(diag(GAM{1}));
0122     sy=sy*sy';
0123     if useautocorr,
0124         sy=sy-diag(diag(sy))+eye(length(mf));
0125         GAM{1}=GAM{1}./sy;
0126     else
0127         for j=1:nlags,
0128             GAM{j+1}=GAM{j+1}.*sy;
0129         end
0130     end
0131     gam = dyn_vech(GAM{1});
0132     for j=1:nlags,
0133         gam = [gam; vec(GAM{j+1})];
0134     end
0135 end
0136 gam = [oo_.dr.ys(oo_.dr.order_var(mf)); gam];
0137 
0138 %   if useautocorr,
0139 warning('on','MATLAB:divideByZero')
0140 %   end