diff --git a/matlab/getH.m b/matlab/getH.m index 5d3b06efa..c7d9caef1 100644 --- a/matlab/getH.m +++ b/matlab/getH.m @@ -1,4 +1,4 @@ -function [H, dA, dOm, Hss, info] = getH(A, B, M_,oo_,kronflag,indx,indexo) +function [H, dA, dOm, Hss, gp, info] = getH(A, B, M_,oo_,kronflag,indx,indexo) % computes derivative of reduced form linear model w.r.t. deep params % diff --git a/matlab/getJJ.m b/matlab/getJJ.m index d2f9f8a75..d93316066 100644 --- a/matlab/getJJ.m +++ b/matlab/getJJ.m @@ -1,4 +1,4 @@ -function [JJ, H, gam] = getJJ(A, B, M_,oo_,options_,kronflag,indx,indexo,mf,nlags,useautocorr) +function [JJ, H, gam, gp] = getJJ(A, B, M_,oo_,options_,kronflag,indx,indexo,mf,nlags,useautocorr) % Copyright (C) 2009 Dynare Team % @@ -22,118 +22,118 @@ if nargin<8 | isempty(indexo), indexo = [];, end, if nargin<10 | isempty(nlags), nlags=3; end, if nargin<11 | isempty(useautocorr), useautocorr=0; end, -if useautocorr, + if useautocorr, warning('off','MATLAB:divideByZero') -end + end if kronflag == -1, - fun = 'thet2tau'; - params0 = M_.params; - JJ = fdjac(fun,[sqrt(diag(M_.Sigma_e(indexo,indexo))); M_.params(indx)],indx,indexo,1,mf,nlags,useautocorr); - M_.params = params0; - params0 = M_.params; - H = fdjac(fun,[sqrt(diag(M_.Sigma_e(indexo,indexo))); M_.params(indx)],indx,indexo,0,mf,nlags,useautocorr); - M_.params = params0; - assignin('base','M_', M_); - assignin('base','oo_', oo_); + fun = 'thet2tau'; + params0 = M_.params; + JJ = fdjac(fun,[sqrt(diag(M_.Sigma_e(indexo,indexo))); M_.params(indx)],indx,indexo,1,mf,nlags,useautocorr); + M_.params = params0; + params0 = M_.params; + H = fdjac(fun,[sqrt(diag(M_.Sigma_e(indexo,indexo))); M_.params(indx)],indx,indexo,0,mf,nlags,useautocorr); + M_.params = params0; + assignin('base','M_', M_); + assignin('base','oo_', oo_); else - [H, dA, dOm, dYss] = getH(A, B, M_,oo_,kronflag,indx,indexo); - % if isempty(H), - % JJ = []; - % GAM = []; - % return - % end - m = length(A); - - GAM = lyapunov_symm(A,B*M_.Sigma_e*B',options_.qz_criterium,options_.lyapunov_complex_threshold,1); - k = find(abs(GAM) < 1e-12); - GAM(k) = 0; - % if useautocorr, - sdy = sqrt(diag(GAM)); - sy = sdy*sdy'; - % end - - % BB = dOm*0; - % for j=1:length(indx), - % BB(:,:,j)= dA(:,:,j)*GAM*A'+A*GAM*dA(:,:,j)'+dOm(:,:,j); - % end - % XX = lyapunov_symm_mr(A,BB,options_.qz_criterium,options_.lyapunov_complex_threshold,0); - for j=1:length(indexo), - dum = lyapunov_symm(A,dOm(:,:,j),options_.qz_criterium,options_.lyapunov_complex_threshold,2); - % dum = XX(:,:,j); - k = find(abs(dum) < 1e-12); - dum(k) = 0; - if useautocorr - dsy = 1/2./sdy.*diag(dum); - dsy = dsy*sdy'+sdy*dsy'; - dum1=dum; - dum1 = (dum1.*sy-dsy.*GAM)./(sy.*sy); - dum1 = dum1-diag(diag(dum1))+diag(diag(dum)); - dumm = vech(dum1(mf,mf)); - else - dumm = vech(dum(mf,mf)); - end - for i=1:nlags, - dum1 = A^i*dum; - if useautocorr - dum1 = (dum1.*sy-dsy.*(A^i*GAM))./(sy.*sy); - end - dumm = [dumm; vec(dum1(mf,mf))]; - end - JJ(:,j) = dumm; + [H, dA, dOm, dYss, gp] = getH(A, B, M_,oo_,kronflag,indx,indexo); + gp = reshape(gp,size(gp,1)*size(gp,2),size(gp,3)); + % if isempty(H), + % JJ = []; + % GAM = []; + % return + % end + m = length(A); + GAM = lyapunov_symm(A,B*M_.Sigma_e*B',options_.qz_criterium,options_.lyapunov_complex_threshold,1); + k = find(abs(GAM) < 1e-12); + GAM(k) = 0; +% if useautocorr, + sdy = sqrt(diag(GAM)); + sy = sdy*sdy'; +% end + +% BB = dOm*0; +% for j=1:length(indx), +% BB(:,:,j)= dA(:,:,j)*GAM*A'+A*GAM*dA(:,:,j)'+dOm(:,:,j); +% end +% XX = lyapunov_symm_mr(A,BB,options_.qz_criterium,options_.lyapunov_complex_threshold,0); + for j=1:length(indexo), + dum = lyapunov_symm(A,dOm(:,:,j),options_.qz_criterium,options_.lyapunov_complex_threshold,2); +% dum = XX(:,:,j); + k = find(abs(dum) < 1e-12); + dum(k) = 0; + if useautocorr + dsy = 1/2./sdy.*diag(dum); + dsy = dsy*sdy'+sdy*dsy'; + dum1=dum; + dum1 = (dum1.*sy-dsy.*GAM)./(sy.*sy); + dum1 = dum1-diag(diag(dum1))+diag(diag(dum)); + dumm = vech(dum1(mf,mf)); + else + dumm = vech(dum(mf,mf)); end - nexo = length(indexo); - for j=1:length(indx), - 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); - % dum = XX(:,:,j); - k = find(abs(dum) < 1e-12); - dum(k) = 0; - if useautocorr - dsy = 1/2./sdy.*diag(dum); - dsy = dsy*sdy'+sdy*dsy'; - dum1=dum; - dum1 = (dum1.*sy-dsy.*GAM)./(sy.*sy); - dum1 = dum1-diag(diag(dum1))+diag(diag(dum)); - dumm = vech(dum1(mf,mf)); - else - dumm = vech(dum(mf,mf)); - end - for i=1:nlags, - dum1 = A^i*dum; - for ii=1:i, - dum1 = dum1 + A^(ii-1)*dA(:,:,j+nexo)*A^(i-ii)*GAM; - end - if useautocorr - dum1 = (dum1.*sy-dsy.*(A^i*GAM))./(sy.*sy); - end - dumm = [dumm; vec(dum1(mf,mf))]; - end - JJ(:,j+nexo) = dumm; + for i=1:nlags, + dum1 = A^i*dum; + if useautocorr + dum1 = (dum1.*sy-dsy.*(A^i*GAM))./(sy.*sy); + end + dumm = [dumm; vec(dum1(mf,mf))]; end - - JJ = [ [zeros(length(mf),nexo) dYss(mf,:)]; JJ]; - - if nargout >2, - % sy=sy(mf,mf); - options_.ar=nlags; - [GAM,stationary_vars] = th_autocovariances(oo_.dr,oo_.dr.order_var(mf),M_,options_); - sy=sqrt(diag(GAM{1})); - sy=sy*sy'; - if useautocorr, - sy=sy-diag(diag(sy))+eye(length(mf)); - GAM{1}=GAM{1}./sy; - else - for j=1:nlags, - GAM{j+1}=GAM{j+1}.*sy; - end - end - gam = vech(GAM{1}); - for j=1:nlags, - gam = [gam; vec(GAM{j+1})]; - end + JJ(:,j) = dumm; + end + nexo = length(indexo); + for j=1:length(indx), + 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); +% dum = XX(:,:,j); + k = find(abs(dum) < 1e-12); + dum(k) = 0; + if useautocorr + dsy = 1/2./sdy.*diag(dum); + dsy = dsy*sdy'+sdy*dsy'; + dum1=dum; + dum1 = (dum1.*sy-dsy.*GAM)./(sy.*sy); + dum1 = dum1-diag(diag(dum1))+diag(diag(dum)); + dumm = vech(dum1(mf,mf)); + else + dumm = vech(dum(mf,mf)); end - gam = [oo_.dr.ys(oo_.dr.order_var(mf)); gam]; + for i=1:nlags, + dum1 = A^i*dum; + for ii=1:i, + dum1 = dum1 + A^(ii-1)*dA(:,:,j+nexo)*A^(i-ii)*GAM; + end + if useautocorr + dum1 = (dum1.*sy-dsy.*(A^i*GAM))./(sy.*sy); + end + dumm = [dumm; vec(dum1(mf,mf))]; + end + JJ(:,j+nexo) = dumm; + end + + JJ = [ [zeros(length(mf),nexo) dYss(mf,:)]; JJ]; + + if nargout >2, +% sy=sy(mf,mf); + options_.ar=nlags; + [GAM,stationary_vars] = th_autocovariances(oo_.dr,oo_.dr.order_var(mf),M_,options_); + sy=sqrt(diag(GAM{1})); + sy=sy*sy'; + if useautocorr, + sy=sy-diag(diag(sy))+eye(length(mf)); + GAM{1}=GAM{1}./sy; + else + for j=1:nlags, + GAM{j+1}=GAM{j+1}.*sy; + end + end + gam = vech(GAM{1}); + for j=1:nlags, + gam = [gam; vec(GAM{j+1})]; + end + end + gam = [oo_.dr.ys(oo_.dr.order_var(mf)); gam]; end -if useautocorr, + if useautocorr, warning('on','MATLAB:divideByZero') -end \ No newline at end of file + end \ No newline at end of file diff --git a/matlab/identification_checks.m b/matlab/identification_checks.m index 376c3b449..7ef13f1a8 100644 --- a/matlab/identification_checks.m +++ b/matlab/identification_checks.m @@ -1,4 +1,4 @@ -function [McoH, McoJ, PcoH, PcoJ, condH, condJ, eH, eJ, ind01, ind02, indnoH, indnoJ] = identification_checks(H,JJ, bayestopt_) +function [McoH, McoJ, McoGP, PcoH, PcoJ, PcoGP, condH, condJ, condGP, eH, eJ, eGP, ind01, ind02, indnoH, indnoJ] = identification_checks(H,JJ, gp, bayestopt_) % Copyright (C) 2008 Dynare Team % @@ -22,8 +22,10 @@ function [McoH, McoJ, PcoH, PcoJ, condH, condJ, eH, eJ, ind01, ind02, indnoH, in % 1. check rank of H at theta npar = size(H,2); +npar0 = size(gp,2); indnoH = {}; indnoJ = {}; +indnoLRE = {}; ind1 = find(vnorm(H)~=0); H1 = H(:,ind1); covH = H1'*H1; @@ -50,6 +52,19 @@ eJ(find(vnorm(JJ)==0),1:length(find(vnorm(JJ)==0)))=eye(length(find(vnorm(JJ)==0 % condJ = cond(JJ1'*JJ1); condJ = cond(JJ1); +ind3 = find(vnorm(gp)~=0); +gp1 = gp(:,ind3); +covgp = gp1'*gp1; +sdgp = sqrt(diag(covgp)); +sdgp = sdgp*sdgp'; +[ex1,ex2] = eig( (gp1'*gp1)./sdgp ); +% eJ = NaN(npar,length(ind2)); +eGP = zeros(npar0,npar0); +eGP(ind3,length(find(vnorm(gp)==0))+1:end) = ex1; +eGP(find(vnorm(gp)==0),1:length(find(vnorm(gp)==0)))=eye(length(find(vnorm(gp)==0))); +% condJ = cond(JJ1'*JJ1); +condGP = cond(gp1); + if rank(H)