Merge remote-tracking branch 'ratto/master'

matlab/dynare_identification.m

@@ -61,7 +61,7 @@ options_ident = set_default_option(options_ident,'advanced',0);
 options_ident = set_default_option(options_ident,'normalize_jacobians',1);
 options_ident = set_default_option(options_ident,'lik_init',1);
 if options_ident.gsa_sample_file,
-    GSAFolder = checkpath('gsa');
+    GSAFolder = checkpath('gsa',M_.dname);
     if options_ident.gsa_sample_file==1,
         load([GSAFolder,filesep,fname_,'_prior'],'lpmat','lpmat0','istable');
     elseif options_ident.gsa_sample_file==2,
@@ -105,11 +105,13 @@ nlags = options_ident.ar;
 periods = options_ident.periods;
 replic = options_ident.replic;
 useautocorr = options_ident.useautocorr;
+options_.order=1;
 options_.ar=nlags;
 options_.prior_mc = options_ident.prior_mc;
 options_.options_ident = options_ident;
 options_.Schur_vec_tol = 1.e-8;
 options_.nomoments=0;
+options_.analytic_derivation=1;
 
 options_ = set_default_option(options_,'datafile',[]);
 options_.mode_compute = 0;
@@ -152,7 +154,7 @@ end
 
 SampleSize = options_ident.prior_mc;
 
-if ~(exist('sylvester3mr','file')==2),
+if ~(exist('sylvester3','file')==2),
 
     dynareroot = strrep(which('dynare'),'dynare.m','');
     addpath([dynareroot 'gensylv'])

matlab/gensylv/sylvester3.m

 function x=sylvester3(a,b,c,d)
-% solves a*x+b*x*c=d
+% solves a*x+b*x*c=d where d is [n x m x p]
 
 % Copyright (C) 2005-2009 Dynare Team
 %
@@ -20,53 +20,72 @@ function x=sylvester3(a,b,c,d)
 
 n = size(a,1);
 m = size(c,1);
+p = size(d,3);
+x=zeros(n,m,p);
 if n == 1
-    x=d./(a*ones(1,m)+b*c);
+    for j=1:p,
+        x(:,:,j)=d(:,:,j)./(a*ones(1,m)+b*c);
+    end
     return
 end
 if m == 1
-    x = (a+c*b)\d;
+    for j=1:p,
+        x(:,:,j) = (a+c*b)\d(:,:,j);
+    end
     return;
 end
-x=zeros(n,m);
 [u,t]=schur(c);
 if exist('OCTAVE_VERSION')
     [aa,bb,qq,zz]=qz(full(a),full(b));
-    d=qq'*d*u;
+    for j=1:p,
+        d(:,:,j)=qq'*d(:,:,j)*u;
+    end
 else
     [aa,bb,qq,zz]=qz(full(a),full(b),'real'); % available in Matlab version 6.0
-    d=qq*d*u;
+    for j=1:p,
+        d(:,:,j)=qq*d(:,:,j)*u;
+    end
 end
 i = 1;
+c = zeros(n,1,p);
 while i < m
     if t(i+1,i) == 0
         if i == 1
-            c = zeros(n,1);
+            c = zeros(n,1,p);
         else
-            c = bb*(x(:,1:i-1)*t(1:i-1,i));
+            for j=1:p,
+                c(:,:,j) = bb*(x(:,1:i-1,j)*t(1:i-1,i));
+            end
         end
-        x(:,i)=(aa+bb*t(i,i))\(d(:,i)-c);
+        x(:,i,:)=(aa+bb*t(i,i))\squeeze(d(:,i,:)-c);
         i = i+1;
     else
         if i == n
-            c = zeros(n,1);
-            c1 = zeros(n,1);
+            c = zeros(n,1,p);
+            c1 = zeros(n,1,p);
         else
-            c = bb*(x(:,1:i-1)*t(1:i-1,i));
-            c1 = bb*(x(:,1:i-1)*t(1:i-1,i+1));
+            for j=1:p,
+                c(:,:,j) = bb*(x(:,1:i-1,j)*t(1:i-1,i));
+                c1(:,:,j) = bb*(x(:,1:i-1,j)*t(1:i-1,i+1));
+            end
         end
-        z = [aa+bb*t(i,i) bb*t(i+1,i); bb*t(i,i+1) aa+bb*t(i+1,i+1)]...
-            \[d(:,i)-c;d(:,i+1)-c1];
-        x(:,i) = z(1:n);
-        x(:,i+1) = z(n+1:end);
+        bigmat = ([aa+bb*t(i,i) bb*t(i+1,i); bb*t(i,i+1) aa+bb*t(i+1,i+1)]);
+        z = bigmat\squeeze([d(:,i,:)-c;d(:,i+1,:)-c1]);
+        x(:,i,:) = z(1:n,:);
+        x(:,i+1,:) = z(n+1:end,:);
         i = i + 2;
     end
 end
 if i == m
-    c = bb*(x(:,1:m-1)*t(1:m-1,m));
-    x(:,m)=(aa+bb*t(m,m))\(d(:,m)-c);
+    for j=1:p,
+        c(:,:,j) = bb*(x(:,1:m-1,j)*t(1:m-1,m));
+    end
+    aabbt = (aa+bb*t(m,m));
+    x(:,m,:)=aabbt\squeeze(d(:,m,:)-c);
+end
+for j=1:p,
+    x(:,:,j)=zz*x(:,:,j)*u';
 end
-x=zz*x*u';
 
 % 01/25/03 MJ corrected bug for i==m (sign of c in x determination)
 % 01/31/03 MJ added 'real' to qz call

matlab/gensylv/sylvester3a.m

-function x=sylvester3a(x0,a,b,c,d)
+function [x0, flag]=sylvester3a(x0,a,b,c,dd)
 % solves iteratively ax+bxc=d
 
-% Copyright (C) 2005-2009 Dynare Team
+% Copyright (C) 2005-2012 Dynare Team
 %
 % This file is part of Dynare.
 %
@@ -20,15 +20,19 @@ function x=sylvester3a(x0,a,b,c,d)
 
 a_1 = inv(a);
 b = a_1*b;
-d = a_1*d;
-e = 1;
-iter = 1;
-while e > 1e-8 & iter < 500
-    x = d-b*x0*c;
-    e = max(max(abs(x-x0)));
-    x0 = x;
-    iter = iter + 1;
-end
-if iter == 500
-    warning('sylvester3a : Only accuracy of %10.8f is achieved after 500 iterations') 
+flag=0;
+for j=1:size(dd,3),
+    d = a_1*dd(:,:,j);
+    e = 1;
+    iter = 1;
+    while e > 1e-8 && iter < 500
+        x = d-b*x0(:,:,j)*c;
+        e = max(max(abs(x-x0(:,:,j))));
+        x0(:,:,j) = x;
+        iter = iter + 1;
+    end
+    if iter == 500
+        sprintf('sylvester3a : Only accuracy of %10.8f is achieved after 500 iterations',e);
+        flag=1;
+    end
 end
\ No newline at end of file

matlab/gensylv/sylvester3mr.m

-function x=sylvester3mr(a,b,c,d)
-% solves a*x+b*x*c=d where d is [n x m x p]
-
-% Copyright (C) 2005-2009 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/>.
-
-n = size(a,1);
-m = size(c,1);
-if length(size(d))==2,
-    x=sylvester3(a,b,c,d);
-    return
-end
-p = size(d,3);
-if n == 1
-    for j=1:p,
-        x(:,:,j)=d(:,:,j)./(a*ones(1,m)+b*c);
-    end
-    return
-end
-if m == 1
-    invacb = inv(a+c*b);
-    x = invacb*d;
-    return;
-end
-x=zeros(n,m,p);
-[u,t]=schur(c);
-if exist('OCTAVE_VERSION')
-    [aa,bb,qq,zz]=qz(full(a),full(b));
-    for j=1:p,
-        d(:,:,j)=qq'*d(:,:,j)*u;
-    end
-else
-    [aa,bb,qq,zz]=qz(full(a),full(b),'real'); % available in Matlab version 6.0
-    for j=1:p,
-        d(:,:,j)=qq*d(:,:,j)*u;
-    end
-end
-i = 1;
-c = zeros(n,1,p);
-while i < m
-    if t(i+1,i) == 0
-        if i == 1
-            c = zeros(n,1,p);
-        else
-            for j=1:p,
-                c(:,:,j) = bb*(x(:,1:i-1,j)*t(1:i-1,i));
-            end
-        end
-        %       aabbtinv = inv(aa+bb*t(i,i));
-        %       x(:,i,:)=aabbtinv*squeeze(d(:,i,:)-c);
-        x(:,i,:)=(aa+bb*t(i,i))\squeeze(d(:,i,:)-c);
-        i = i+1;
-    else
-        if i == n
-            c = zeros(n,1,p);
-            c1 = zeros(n,1,p);
-        else
-            for j=1:p,
-                c(:,:,j) = bb*(x(:,1:i-1,j)*t(1:i-1,i));
-                c1(:,:,j) = bb*(x(:,1:i-1,j)*t(1:i-1,i+1));
-            end
-        end
-        %       bigmatinv = inv([aa+bb*t(i,i) bb*t(i+1,i); bb*t(i,i+1) aa+bb*t(i+1,i+1)]);
-        %       z = bigmatinv * squeeze([d(:,i,:)-c;d(:,i+1,:)-c1]);
-        bigmat = ([aa+bb*t(i,i) bb*t(i+1,i); bb*t(i,i+1) aa+bb*t(i+1,i+1)]);
-        z = bigmat\squeeze([d(:,i,:)-c;d(:,i+1,:)-c1]);
-        x(:,i,:) = z(1:n,:);
-        x(:,i+1,:) = z(n+1:end,:);
-        i = i + 2;
-    end
-end
-if i == m
-    for j=1:p,
-        c(:,:,j) = bb*(x(:,1:m-1,j)*t(1:m-1,m));
-    end
-    %       aabbtinv = inv(aa+bb*t(m,m));
-    %     x(:,m,:)=aabbtinv*squeeze(d(:,m,:)-c);
-    aabbt = (aa+bb*t(m,m));
-    x(:,m,:)=aabbt\squeeze(d(:,m,:)-c);
-end
-for j=1:p,
-    x(:,:,j)=zz*x(:,:,j)*u';
-end
-
-% 01/25/03 MJ corrected bug for i==m (sign of c in x determination)
-% 01/31/03 MJ added 'real' to qz call

matlab/getH.m

@@ -355,7 +355,13 @@ else % generalized sylvester equation
         elem(:,:,j) = (Dg0(:,:,j)-Dg1(:,:,j)*A);
         d(:,:,j) = Dg2(:,:,j)-elem(:,:,j)*A;
     end
-    xx=sylvester3mr(a,b,c,d);
+    xx=sylvester3(a,b,c,d);
+    flag=1;
+    icount=0;
+    while flag && icount<4,
+        [xx, flag]=sylvester3a(xx,a,b,c,d);
+        icount=icount+1;
+    end
     H=zeros(m*m+m*(m+1)/2,param_nbr+length(indexo));
     if nargout>1,
         dOm = zeros(m,m,param_nbr+length(indexo));
@@ -435,7 +441,13 @@ if nargout > 5,
         d(:,:,jcount) = elem1+elem2;
         end
     end
-    xx2=sylvester3mr(a,b,c,d);
+    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
     jcount = 0;
     for j=1:param_nbr,
         for i=j:param_nbr,

matlab/identification_analysis.m

@@ -138,6 +138,10 @@ if info(1)==0,
             [fval,cost_flag,ys,trend_coeff,info,M_,options_,bayestopt_,oo_,DLIK,AHess] = dsge_likelihood(params',data_info,options_,M_,estim_params_,bayestopt_,oo_,derivatives_info);
 %                 fval = DsgeLikelihood(xparam1,data_info,options_,M_,estim_params_,bayestopt_,oo_);
             AHess=-AHess;
+            if min(eig(AHess))<0,
+                error('Analytic Hessian is not positive semi-definite!')
+            end
+%             chol(AHess);
             ide_hess.AHess= AHess;
             deltaM = sqrt(diag(AHess));
             iflag=any((deltaM.*deltaM)==0);
@@ -153,7 +157,9 @@ if info(1)==0,
             cmm = NaN(size(siJ,1),size(siJ,1));
             ind1=find(ide_hess.ind0);
             cmm = siJ(:,ind1)*((AHess(ind1,ind1))\siJ(:,ind1)');
-            chh = siH(:,ind1)*((AHess(ind1,ind1))\siH(:,ind1)');
+            temp1=((AHess(ind1,ind1))\siH(:,ind1)');
+            diag_chh=sum(siH(:,ind1)'.*temp1)';
+%             chh = siH(:,ind1)*((AHess(ind1,ind1))\siH(:,ind1)');
             ind1=ind1(ind1>offset);
             clre = siLRE(:,ind1-offset)*((AHess(ind1,ind1))\siLRE(:,ind1-offset)');
             rhoM=sqrt(1./diag(inv(tildaM(indok,indok))));
@@ -190,12 +196,16 @@ if info(1)==0,
 %             rhoM=sqrt(1-1./diag(inv(tildaM)));
 %             rhoM=(1-1./diag(inv(tildaM)));
             ind1=find(ide_hess.ind0);
-            chh = siH(:,ind1)*((MIM(ind1,ind1))\siH(:,ind1)');
+            temp1=((MIM(ind1,ind1))\siH(:,ind1)');
+            diag_chh=sum(siH(:,ind1)'.*temp1)';
+%             chh = siH(:,ind1)*((MIM(ind1,ind1))\siH(:,ind1)');
             ind1=ind1(ind1>offset);
             clre = siLRE(:,ind1-offset)*((MIM(ind1,ind1))\siLRE(:,ind1-offset)');
-            rhoM(indok)=sqrt(1./diag(inv(tildaM(indok,indok))));
-            normaliz(indok) = (sqrt(diag(inv(tildaM(indok,indok))))./deltaM(indok))'; %sqrt(diag(inv(MIM(indok,indok))))';
-%             deltaM = deltaM.*abs(params');
+            if ~isempty(indok),
+                rhoM(indok)=sqrt(1./diag(inv(tildaM(indok,indok))));
+                normaliz(indok) = (sqrt(diag(inv(tildaM(indok,indok))))./deltaM(indok))'; %sqrt(diag(inv(MIM(indok,indok))))';
+            end
+            %             deltaM = deltaM.*abs(params')
             flag_score=0;
         end
         ide_strength_J(indok) = (1./(normaliz(indok)'./abs(params(indok)')));
@@ -212,7 +222,8 @@ if info(1)==0,
 %         isok = find((abs(TAU)>=1.e-8));
 %         quant(isok,:) = siH(isok,:)./repmat(TAU(isok,1),1,nparam);
 %         quant(inok,:) = siH(inok,:)./repmat(mean(abs(TAU)),length(inok),nparam);
-        quant = siH./repmat(sqrt(diag(chh)),1,nparam);
+%         quant = siH./repmat(sqrt(diag(chh)),1,nparam);
+        quant = siH./repmat(sqrt(diag_chh),1,nparam);
         siHnorm = vnorm(quant).*normaliz1;
         %                 siHnorm = vnorm(siH./repmat(TAU,1,nparam)).*normaliz;
         quant=[];

matlab/simulated_moment_uncertainty.m

@@ -22,7 +22,7 @@ global options_ oo_
 
 mm=zeros(length(indx),replic);
 
-disp('Evaluting simulated moment uncertainty ... please wait')
+disp('Evaluating simulated moment uncertainty ... please wait')
 disp(['Doing ',int2str(replic),' replicas of length ',int2str(periods),' periods.'])
 noprint0 = options_.noprint;
 for j=1:replic;