v4.1:: Changes related to lyapunov_symm.

+ Fixed a bug introduced in the previous commit.
+ Added a threshold parameter for the complex blocks in the upper
triangular matrix T in lyapunov_symm.m (new field in options_)
+ Cosmetic changes.


git-svn-id: https://www.dynare.org/svn/dynare/trunk@2414 ac1d8469-bf42-47a9-8791-bf33cf982152

matlab/DsgeLikelihood.m

@@ -170,7 +170,7 @@ function [fval,cost_flag,ys,trend_coeff,info] = DsgeLikelihood(xparam1,gend,data
       if kalman_algo ~= 2
           kalman_algo = 1;
       end
-      Pstar = lyapunov_symm(T,R*Q*R',options_.qz_criterium);
+      Pstar = lyapunov_symm(T,R*Q*R',options_.qz_criterium,options_.lyapunov_complex_threshold);
       Pinf	= [];
   elseif options_.lik_init == 2	% Old Diffuse Kalman filter
       if kalman_algo ~= 2

matlab/DsgeLikelihood_hh.m

@@ -171,7 +171,7 @@ function [fval,llik,cost_flag,ys,trend_coeff,info] = DsgeLikelihood_hh(xparam1,g
       if kalman_algo ~= 2
           kalman_algo = 1;
       end
-      Pstar = lyapunov_symm(T,R*Q*R',options_.qz_criterium);
+      Pstar = lyapunov_symm(T,R*Q*R',options_.qz_criterium,options_.lyapunov_complex_threshold);
       Pinf	= [];
   elseif options_.lik_init == 2	% Old Diffuse Kalman filter
       if kalman_algo ~= 2

matlab/DsgeSmoother.m

@@ -113,7 +113,7 @@ function [alphahat,etahat,epsilonhat,ahat,SteadyState,trend_coeff,aK,T,R,P,PK,d,
       if kalman_algo ~= 2
           kalman_algo = 1;
       end
-      Pstar = lyapunov_symm(T,R*Q*transpose(R),options_.qz_criterium);
+      Pstar = lyapunov_symm(T,R*Q*transpose(R),options_.qz_criterium,options_.lyapunov_complex_threshold);
       Pinf	= [];
   elseif options_.lik_init == 2 % Old Diffuse Kalman filter
       if kalman_algo ~= 2

matlab/DsgeVarLikelihood.m

@@ -137,7 +137,7 @@ end
 %------------------------------------------------------------------------------
 % 3. theoretical moments (second order)
 %------------------------------------------------------------------------------
-tmp0 = lyapunov_symm(T,R*Q*R',options_.qz_criterium);% I compute the variance-covariance matrix
+tmp0 = lyapunov_symm(T,R*Q*R',options_.qz_criterium,options_.lyapunov_complex_threshold);% I compute the variance-covariance matrix
 mf  = bayestopt_.mf1;          % of the restricted state vector.
 
 % Get the non centered second order moments

matlab/UnivariateSpectralDensity.m

@@ -96,7 +96,7 @@ for i=M_.maximum_lag:-1:2
 end
 Gamma = zeros(nvar,cutoff+1);
 [A,B] = kalman_transition_matrix(dr,ikx',1:nx,dr.transition_auxiliary_variables,M_.exo_nbr);
-[vx, u] =  lyapunov_symm(A,B*M_.Sigma_e*B',options_.qz_criterium);
+[vx, u] =  lyapunov_symm(A,B*M_.Sigma_e*B',options_.qz_criterium,options_.lyapunov_complex_threshold);
 iky = iv(ivar);
 if ~isempty(u)
     iky = iky(find(any(abs(ghx(iky,:)*u) < options_.Schur_vec_tol,2)));

matlab/calib_obj.m

@@ -28,11 +28,11 @@ function f=calib_obj(M_.Sigma_e,A,ghu1,ghx,ghu,targets,var_weights,iy,nar)
   b=ghu1*M_.Sigma_e*ghu1';
   vx = [];
   if isempty(vx)
-    vx = lyapunov_symm(A,b,options_.qz_criterium);
+    vx = lyapunov_symm(A,b,options_.qz_criterium,options_.lyapunov_complex_threshold);
   else
     [vx,status] = bicgstab_(@f_var,b(:),vx(:),1e-8,50,A,nx);
     if status
-      vx = lyapunov_symm(A,b,options_.qz_criterium);
+      vx = lyapunov_symm(A,b,options_.qz_criterium,options_.lyapunov_complex_threshold);
     else
       vx=reshape(vx,nx,nx);
     end

matlab/calib_obj2.m

@@ -26,7 +26,7 @@ function objective=calib_obj2(M_.Sigma_e,A,ghu1,ghx,ghu,targets,var_weights,iy,n
   M_.Sigma_e=diag(M_.Sigma_e);
   nx = size(ghx,2);
   b=ghu1*M_.Sigma_e*ghu1';
-  vx = lyapunov_symm(A,b,options_.qz_criterium);
+  vx = lyapunov_symm(A,b,options_.qz_criterium,options_.lyapunov_complex_threshold);
   oo_.gamma_y{1} = ghx*vx*ghx'+ ghu*M_.Sigma_e*ghu';
   if ~isempty(targets{1})
     objective{1} = sqrt(oo_.gamma_y{1}(iy{1}));

matlab/dsge_posterior_kernel.m

@@ -170,7 +170,7 @@ function [fval,cost_flag,ys,trend_coeff,info] = DsgeLikelihood(xparam1,gend,data
       if kalman_algo ~= 2
           kalman_algo = 1;
       end
-      Pstar = lyapunov_symm(T,R*Q*R',options_.qz_criterium);
+      Pstar = lyapunov_symm(T,R*Q*R',options_.qz_criterium,options_.lyapunov_complex_threshold);
       Pinf	= [];
   elseif options_.lik_init == 2	% Old Diffuse Kalman filter
       if kalman_algo ~= 2

matlab/get_variance_of_endogenous_variables.m

@@ -45,7 +45,7 @@ function [vx1,i_ns] = get_variance_of_endogenous_variables(dr,i_var)
   
   [A,B] = kalman_transition_matrix(dr,nstatic+(1:npred),1:nc,dr.transition_auxiliary_variables,M_.exo_nbr);
   
-  [vx,u] = lyapunov_symm(A,B*Sigma_e*B',options_.qz_criterium);
+  [vx,u] = lyapunov_symm(A,B*Sigma_e*B',options_.qz_criterium,options_.lyapunov_complex_threshold);
   
   if size(u,2) > 0
     i_stat = find(any(abs(ghx*u) < options_.Schur_vec_tol,2));

matlab/global_initialization.m

@@ -44,6 +44,7 @@ function global_initialization()
   options_.Schur_vec_tol = 1e-8; % used to find nonstationary variables
                                  % in Schur decomposition of the
                                  % transition matrix
+  options_.lyapunov_complex_threshold = 1e-15;
   options_.solve_tolf = eps^(1/3);
   options_.solve_tolx = eps^(2/3);
   options_.solve_maxit = 500;
@@ -112,7 +113,7 @@ function global_initialization()
   options_.drop = 100;
   options_.simul_algo = 0;
   options_.model_mode = 0;
-  % if mjdgges.dll (or .mexw32 or ....) doesn't exit, matlab/qz is added to the path. 
+  % if mjdgges.dll (or .mexw32 or ....) doesn't exist, matlab/qz is added to the path. 
   % There exists now qz/mjdgges.m that contains the calls to the old Sims code 
   % Hence, if mjdgges.m is visible exist(...)==2, 
   % this means that the DLL isn't avaiable and use_qzdiv is set to 1
@@ -151,6 +152,7 @@ function global_initialization()
   options_.loglinear = 0;
   options_.markowitz = 0.5;
   options_.mh_conf_sig = 0.90;
+  options_.prior_interval = 0.90;
   options_.mh_drop = 0.5;
   options_.mh_jscale = 0.2;
   options_.mh_init_scale = 2*options_.mh_jscale;

matlab/lyapunov_symm.m

-function [x,u] = lyapunov_symm(a,b,qz_criterium,method)
+function [x,u] = lyapunov_symm(a,b,qz_criterium,lyapunov_complex_threshold,method)
 % Solves the Lyapunov equation x-a*x*a' = b, for b and x symmetric matrices.
 % If a has some unit roots, the function computes only the solution of the stable subsystem.
 %  
 % INPUTS:
-%   a        [double]    n*n matrix.  
-%   b        [double]    n*n matrix
-%   method   [integer]   Scalar, if method=0 [default] then U, T, n and k are not persistent.  
-%                                   method=1 then U, T, n and k are declared as persistent variables and the schur decomposition is triggered.    
-%                                   method=2 then U, T, n and k are declared as persistent variables and the schur decomposition is not performed.
+%   a                           [double]    n*n matrix.
+%   b                           [double]    n*n matrix.
+%   qz_criterium                [double]    scalar, unit root threshold for eigenvalues in a.
+%   lyapunov_complex_threshold  [double]    scalar, complex block threshold for the upper triangular matrix T.
+%   method                      [integer]   Scalar, if method=0 [default] then U, T, n and k are not persistent.  
+%                                                      method=1 then U, T, n and k are declared as persistent 
+%                                                               variables and the schur decomposition is triggered.    
+%                                                      method=2 then U, T, n and k are declared as persistent 
+%                                                               variables and the schur decomposition is not performed.
 % OUTPUTS
 %   x:      [double]    m*m solution matrix of the lyapunov equation, where m is the dimension of the stable subsystem.
 %   u:      [double]    Schur vectors associated with unit roots  
@@ -34,7 +38,7 @@ function [x,u] = lyapunov_symm(a,b,qz_criterium,method)
 %
 % You should have received a copy of the GNU General Public License
 % along with Dynare.  If not, see <http://www.gnu.org/licenses/>.
-    if nargin<4
+    if nargin<5
         method = 0;
     end
     if method
@@ -64,25 +68,21 @@ function [x,u] = lyapunov_symm(a,b,qz_criterium,method)
         if exist('ordschur','builtin')
             % Selects stable roots
             [U,T] = ordschur(U,T,e1);
-            % B = U(:,k+1:end)'*b*U(:,k+1:end);
             T = T(k+1:end,k+1:end);
         elseif k > 0
             % Problem for Matlab version that don't have ordschur
             error(['lyapunov_sym: you need a Matlab version > 6.5 to handle models' ...
                    ' with unit roots'])
-            %else
-            % No unit root
-            %n = length(e1);
-            %B=u'*b*u;
         end
     end
-    B = U(:,k+1:end)'*b*U(:,k+1:end);    
+    
+    B = U(:,k+1:end)'*b*U(:,k+1:end);
     
     x = zeros(n,n);
     i = n;
-
+    
     while i >= 2
-        if T(i,i-1) == 0
+        if abs(T(i,i-1))<lyapunov_complex_threshold
             if i == n
                 c = zeros(n,1);
             else
@@ -105,8 +105,8 @@ function [x,u] = lyapunov_symm(a,b,qz_criterium,method)
                      T(i-1,i-1)*T(1:i,i+1:end)*x(i+1:end,i-1) + ...
                      T(i-1,i)*T(1:i,i+1:end)*x(i+1:end,i);
             end
-            q = [eye(i)-T(1:i,1:i)*T(i,i) - T(1:i,1:i)*T(i,i-1) ; ...
-                 -T(1:i,1:i)*T(i-1,i) eye(i) - T(1:i,1:i)*T(i-1,i-1) ];
+            q = [  eye(i)-T(1:i,1:i)*T(i,i) ,  -T(1:i,1:i)*T(i,i-1) ; ...
+                   -T(1:i,1:i)*T(i-1,i)     ,   eye(i)-T(1:i,1:i)*T(i-1,i-1) ];
             z =  q\[ B(1:i,i)+c ; B(1:i,i-1) + c1 ];
             x(1:i,i) = z(1:i);
             x(1:i,i-1) = z(i+1:end);

matlab/th_autocovariances.m

@@ -93,7 +93,7 @@ function [Gamma_y,ivar] = th_autocovariances(dr,ivar,M_,options_,nodecomposition
     % state space representation for state variables only
     [A,B] = kalman_transition_matrix(dr,ipred,1:nx,dr.transition_auxiliary_variables,M_.exo_nbr);
     if options_.order == 2 | options_.hp_filter == 0
-        [vx, u] =  lyapunov_symm(A,B*M_.Sigma_e*B',options_.qz_criterium);
+        [vx, u] =  lyapunov_symm(A,B*M_.Sigma_e*B',options_.qz_criterium,options_.lyapunov_complex_threshold);
         iky = iv(ivar);
         if ~isempty(u)
             iky = iky(find(all(abs(ghx(iky,:)*u) < options_.Schur_vec_tol,2)));
@@ -132,10 +132,10 @@ function [Gamma_y,ivar] = th_autocovariances(dr,ivar,M_,options_,nodecomposition
             b1 = b1*cs;
             b2(:,exo_names_orig_ord) = ghu(iky,:);
             b2 = b2*cs;
-            vx  = lyapunov_symm(A,b1*b1',options_.qz_criterium,1);
+            vx  = lyapunov_symm(A,b1*b1',options_.qz_criterium,options_.lyapunov_complex_threshold,1);
             vv = diag(aa*vx*aa'+b2*b2');
             for i=1:M_.exo_nbr
-                vx1 = lyapunov_symm(A,b1(:,i)*b1(:,i)',options_.qz_criterium,2);
+                vx1 = lyapunov_symm(A,b1(:,i)*b1(:,i)',options_.qz_criterium,options_.lyapunov_complex_threshold,2);
                 Gamma_y{nar+2}(:,i) = abs(diag(aa*vx1*aa'+b2(:,i)*b2(:,i)'))./vv;
             end
         end