Merge branch 'pr#1354'
commit
c893c0620c
|
@ -156,15 +156,7 @@ if options_.lik_init == 1 % Kalman filter
|
||||||
if kalman_algo ~= 2
|
if kalman_algo ~= 2
|
||||||
kalman_algo = 1;
|
kalman_algo = 1;
|
||||||
end
|
end
|
||||||
if options_.lyapunov_fp == 1
|
Pstar=lyapunov_solver(T,R,Q,options_);
|
||||||
Pstar = lyapunov_symm(T,R*Q*R',options_.lyapunov_fixed_point_tol,options_.qz_criterium,options_.lyapunov_complex_threshold, 3, [], options_.debug);
|
|
||||||
elseif options_.lyapunov_db == 1
|
|
||||||
Pstar = disclyap_fast(T,R*Q*R',options_.lyapunov_doubling_tol);
|
|
||||||
elseif options_.lyapunov_srs == 1
|
|
||||||
Pstar = lyapunov_symm(T,Q,options_.lyapunov_fixed_point_tol,options_.qz_criterium,options_.lyapunov_complex_threshold, 4, R, options_.debug);
|
|
||||||
else
|
|
||||||
Pstar = lyapunov_symm(T,R*Q*R',options_.lyapunov_fixed_point_tol,options_.qz_criterium,options_.lyapunov_complex_threshold, [], [], options_.debug);
|
|
||||||
end;
|
|
||||||
Pinf = [];
|
Pinf = [];
|
||||||
elseif options_.lik_init == 2 % Old Diffuse Kalman filter
|
elseif options_.lik_init == 2 % Old Diffuse Kalman filter
|
||||||
if kalman_algo ~= 2
|
if kalman_algo ~= 2
|
||||||
|
@ -210,15 +202,7 @@ elseif options_.lik_init == 5 % Old diffuse Kalman filter only for th
|
||||||
end
|
end
|
||||||
R_tmp = R(stable, :);
|
R_tmp = R(stable, :);
|
||||||
T_tmp = T(stable,stable);
|
T_tmp = T(stable,stable);
|
||||||
if options_.lyapunov_fp == 1
|
Pstar_tmp=lyapunov_solver(T_tmp,R_tmp,Q,DynareOptions);
|
||||||
Pstar_tmp = lyapunov_symm(T_tmp,R_tmp*Q*R_tmp',options_.lyapunov_fixed_point_tol,options_.qz_criterium,options_.lyapunov_complex_threshold, 3, [], options_.debug);
|
|
||||||
elseif options_.lyapunov_db == 1
|
|
||||||
Pstar_tmp = disclyap_fast(T_tmp,R_tmp*Q*R_tmp',options_.lyapunov_doubling_tol);
|
|
||||||
elseif options_.lyapunov_srs == 1
|
|
||||||
Pstar_tmp = lyapunov_symm(T_tmp,Q,options_.lyapunov_fixed_point_tol,options_.qz_criterium,options_.lyapunov_complex_threshold, 4, R_tmp, options_.debug);
|
|
||||||
else
|
|
||||||
Pstar_tmp = lyapunov_symm(T_tmp,R_tmp*Q*R_tmp',options_.lyapunov_fixed_point_tol,options_.qz_criterium,options_.lyapunov_complex_threshold, [], [], options_.debug);
|
|
||||||
end
|
|
||||||
Pstar(stable, stable) = Pstar_tmp;
|
Pstar(stable, stable) = Pstar_tmp;
|
||||||
Pinf = [];
|
Pinf = [];
|
||||||
end
|
end
|
||||||
|
|
|
@ -95,7 +95,7 @@ for i=M_.maximum_lag:-1:2
|
||||||
end
|
end
|
||||||
|
|
||||||
[A,B] = kalman_transition_matrix(oo_.dr,ikx',1:nx,M_.exo_nbr);
|
[A,B] = kalman_transition_matrix(oo_.dr,ikx',1:nx,M_.exo_nbr);
|
||||||
[vx, u] = lyapunov_symm(A,B*M_.Sigma_e*B',options_.lyapunov_fixed_point_tol,options_.qz_criterium,options_.lyapunov_complex_threshold,[],[],options_.debug);
|
[vx, u] = lyapunov_symm(A,B*M_.Sigma_e*B',options_.lyapunov_fixed_point_tol,options_.qz_criterium,options_.lyapunov_complex_threshold,[],options_.debug);
|
||||||
iky = iv(ivar);
|
iky = iv(ivar);
|
||||||
if ~isempty(u)
|
if ~isempty(u)
|
||||||
iky = iky(find(any(abs(ghx(iky,:)*u) < options_.Schur_vec_tol,2)));
|
iky = iky(find(any(abs(ghx(iky,:)*u) < options_.Schur_vec_tol,2)));
|
||||||
|
|
|
@ -111,7 +111,7 @@ function [fval,info,exit_flag,DLIK,Hess,SteadyState,trend_coeff,Model,DynareOpti
|
||||||
%! @sp 2
|
%! @sp 2
|
||||||
%! @strong{This function calls:}
|
%! @strong{This function calls:}
|
||||||
%! @sp 1
|
%! @sp 1
|
||||||
%! @ref{dynare_resolve}, @ref{lyapunov_symm}, @ref{compute_Pinf_Pstar}, @ref{kalman_filter_d}, @ref{missing_observations_kalman_filter_d}, @ref{univariate_kalman_filter_d}, @ref{kalman_steady_state}, @ref{getH}, @ref{kalman_filter}, @ref{score}, @ref{AHessian}, @ref{missing_observations_kalman_filter}, @ref{univariate_kalman_filter}, @ref{priordens}
|
%! @ref{dynare_resolve}, @ref{lyapunov_symm}, @ref{lyapunov_solver}, @ref{compute_Pinf_Pstar}, @ref{kalman_filter_d}, @ref{missing_observations_kalman_filter_d}, @ref{univariate_kalman_filter_d}, @ref{kalman_steady_state}, @ref{getH}, @ref{kalman_filter}, @ref{score}, @ref{AHessian}, @ref{missing_observations_kalman_filter}, @ref{univariate_kalman_filter}, @ref{priordens}
|
||||||
%! @end deftypefn
|
%! @end deftypefn
|
||||||
%@eod:
|
%@eod:
|
||||||
|
|
||||||
|
@ -351,15 +351,7 @@ switch DynareOptions.lik_init
|
||||||
% Use standard kalman filter except if the univariate filter is explicitely choosen.
|
% Use standard kalman filter except if the univariate filter is explicitely choosen.
|
||||||
kalman_algo = 1;
|
kalman_algo = 1;
|
||||||
end
|
end
|
||||||
if DynareOptions.lyapunov_fp == 1
|
Pstar=lyapunov_solver(T,R,Q,DynareOptions);
|
||||||
Pstar = lyapunov_symm(T,R*Q'*R',DynareOptions.lyapunov_fixed_point_tol,DynareOptions.qz_criterium,DynareOptions.lyapunov_complex_threshold, 3, [], DynareOptions.debug);
|
|
||||||
elseif DynareOptions.lyapunov_db == 1
|
|
||||||
Pstar = disclyap_fast(T,R*Q*R',DynareOptions.lyapunov_doubling_tol);
|
|
||||||
elseif DynareOptions.lyapunov_srs == 1
|
|
||||||
Pstar = lyapunov_symm(T,Q,DynareOptions.lyapunov_fixed_point_tol,DynareOptions.qz_criterium,DynareOptions.lyapunov_complex_threshold, 4, R, DynareOptions.debug);
|
|
||||||
else
|
|
||||||
Pstar = lyapunov_symm(T,R*Q*R',DynareOptions.lyapunov_fixed_point_tol,DynareOptions.qz_criterium,DynareOptions.lyapunov_complex_threshold, [], [], DynareOptions.debug);
|
|
||||||
end;
|
|
||||||
Pinf = [];
|
Pinf = [];
|
||||||
a = zeros(mm,1);
|
a = zeros(mm,1);
|
||||||
Zflag = 0;
|
Zflag = 0;
|
||||||
|
@ -471,7 +463,7 @@ switch DynareOptions.lik_init
|
||||||
if err
|
if err
|
||||||
disp(['dsge_likelihood:: I am not able to solve the Riccati equation, so I switch to lik_init=1!']);
|
disp(['dsge_likelihood:: I am not able to solve the Riccati equation, so I switch to lik_init=1!']);
|
||||||
DynareOptions.lik_init = 1;
|
DynareOptions.lik_init = 1;
|
||||||
Pstar = lyapunov_symm(T,R*Q*R',DynareOptions.lyapunov_fixed_point_tol,DynareOptions.qz_criterium,DynareOptions.lyapunov_complex_threshold, [], [], DynareOptions.debug);
|
Pstar=lyapunov_solver(T,R,Q,DynareOptions);
|
||||||
end
|
end
|
||||||
Pinf = [];
|
Pinf = [];
|
||||||
a = zeros(mm,1);
|
a = zeros(mm,1);
|
||||||
|
@ -491,15 +483,7 @@ switch DynareOptions.lik_init
|
||||||
end
|
end
|
||||||
R_tmp = R(stable, :);
|
R_tmp = R(stable, :);
|
||||||
T_tmp = T(stable,stable);
|
T_tmp = T(stable,stable);
|
||||||
if DynareOptions.lyapunov_fp == 1
|
Pstar_tmp=lyapunov_solver(T_tmp,R_tmp,Q,DynareOptions);
|
||||||
Pstar_tmp = lyapunov_symm(T_tmp,R_tmp*Q*R_tmp',DynareOptions.lyapunov_fixed_point_tol,DynareOptions.qz_criterium,DynareOptions.lyapunov_complex_threshold, 3, [], DynareOptions.debug);
|
|
||||||
elseif DynareOptions.lyapunov_db == 1
|
|
||||||
Pstar_tmp = disclyap_fast(T_tmp,R_tmp*Q*R_tmp',DynareOptions.lyapunov_doubling_tol);
|
|
||||||
elseif DynareOptions.lyapunov_srs == 1
|
|
||||||
Pstar_tmp = lyapunov_symm(T_tmp,Q,DynareOptions.lyapunov_fixed_point_tol,DynareOptions.qz_criterium,DynareOptions.lyapunov_complex_threshold, 4, R_tmp, DynareOptions.debug);
|
|
||||||
else
|
|
||||||
Pstar_tmp = lyapunov_symm(T_tmp,R_tmp*Q*R_tmp',DynareOptions.qz_criterium,DynareOptions.qz_criterium,DynareOptions.lyapunov_complex_threshold, [], [], DynareOptions.debug);
|
|
||||||
end
|
|
||||||
Pstar(stable, stable) = Pstar_tmp;
|
Pstar(stable, stable) = Pstar_tmp;
|
||||||
Pinf = [];
|
Pinf = [];
|
||||||
a = zeros(mm,1);
|
a = zeros(mm,1);
|
||||||
|
@ -580,14 +564,14 @@ if analytic_derivation,
|
||||||
for i=1:EstimatedParameters.nvx
|
for i=1:EstimatedParameters.nvx
|
||||||
k =EstimatedParameters.var_exo(i,1);
|
k =EstimatedParameters.var_exo(i,1);
|
||||||
DQ(k,k,i) = 2*sqrt(Q(k,k));
|
DQ(k,k,i) = 2*sqrt(Q(k,k));
|
||||||
dum = lyapunov_symm(T,DOm(:,:,i),DynareOptions.lyapunov_fixed_point_tol,DynareOptions.qz_criterium,DynareOptions.lyapunov_complex_threshold,[],[],DynareOptions.debug);
|
dum = lyapunov_symm(T,DOm(:,:,i),DynareOptions.lyapunov_fixed_point_tol,DynareOptions.qz_criterium,DynareOptions.lyapunov_complex_threshold,[],DynareOptions.debug);
|
||||||
% kk = find(abs(dum) < 1e-12);
|
% kk = find(abs(dum) < 1e-12);
|
||||||
% dum(kk) = 0;
|
% dum(kk) = 0;
|
||||||
DP(:,:,i)=dum;
|
DP(:,:,i)=dum;
|
||||||
if full_Hess
|
if full_Hess
|
||||||
for j=1:i,
|
for j=1:i,
|
||||||
jcount=jcount+1;
|
jcount=jcount+1;
|
||||||
dum = lyapunov_symm(T,dyn_unvech(D2Om(:,jcount)),DynareOptions.lyapunov_fixed_point_tol,DynareOptions.qz_criterium,DynareOptions.lyapunov_complex_threshold,[],[],DynareOptions.debug);
|
dum = lyapunov_symm(T,dyn_unvech(D2Om(:,jcount)),DynareOptions.lyapunov_fixed_point_tol,DynareOptions.qz_criterium,DynareOptions.lyapunov_complex_threshold,[],DynareOptions.debug);
|
||||||
% kk = (abs(dum) < 1e-12);
|
% kk = (abs(dum) < 1e-12);
|
||||||
% dum(kk) = 0;
|
% dum(kk) = 0;
|
||||||
D2P(:,jcount)=dyn_vech(dum);
|
D2P(:,jcount)=dyn_vech(dum);
|
||||||
|
@ -607,7 +591,7 @@ if analytic_derivation,
|
||||||
offset = offset + EstimatedParameters.nvn;
|
offset = offset + EstimatedParameters.nvn;
|
||||||
if DynareOptions.lik_init==1,
|
if DynareOptions.lik_init==1,
|
||||||
for j=1:EstimatedParameters.np
|
for j=1:EstimatedParameters.np
|
||||||
dum = lyapunov_symm(T,DT(:,:,j+offset)*Pstar*T'+T*Pstar*DT(:,:,j+offset)'+DOm(:,:,j+offset),DynareOptions.lyapunov_fixed_point_tol,DynareOptions.qz_criterium,DynareOptions.lyapunov_complex_threshold,[],[],DynareOptions.debug);
|
dum = lyapunov_symm(T,DT(:,:,j+offset)*Pstar*T'+T*Pstar*DT(:,:,j+offset)'+DOm(:,:,j+offset),DynareOptions.lyapunov_fixed_point_tol,DynareOptions.qz_criterium,DynareOptions.lyapunov_complex_threshold,[],DynareOptions.debug);
|
||||||
% kk = find(abs(dum) < 1e-12);
|
% kk = find(abs(dum) < 1e-12);
|
||||||
% dum(kk) = 0;
|
% dum(kk) = 0;
|
||||||
DP(:,:,j+offset)=dum;
|
DP(:,:,j+offset)=dum;
|
||||||
|
@ -621,7 +605,7 @@ if analytic_derivation,
|
||||||
D2Tij = reshape(D2T(:,jcount),size(T));
|
D2Tij = reshape(D2T(:,jcount),size(T));
|
||||||
D2Omij = dyn_unvech(D2Om(:,jcount));
|
D2Omij = dyn_unvech(D2Om(:,jcount));
|
||||||
tmp = D2Tij*Pstar*T' + T*Pstar*D2Tij' + DTi*DPj*T' + DTj*DPi*T' + T*DPj*DTi' + T*DPi*DTj' + DTi*Pstar*DTj' + DTj*Pstar*DTi' + D2Omij;
|
tmp = D2Tij*Pstar*T' + T*Pstar*D2Tij' + DTi*DPj*T' + DTj*DPi*T' + T*DPj*DTi' + T*DPi*DTj' + DTi*Pstar*DTj' + DTj*Pstar*DTi' + D2Omij;
|
||||||
dum = lyapunov_symm(T,tmp,DynareOptions.lyapunov_fixed_point_tol,DynareOptions.qz_criterium,DynareOptions.lyapunov_complex_threshold,[],[],DynareOptions.debug);
|
dum = lyapunov_symm(T,tmp,DynareOptions.lyapunov_fixed_point_tol,DynareOptions.qz_criterium,DynareOptions.lyapunov_complex_threshold,[],DynareOptions.debug);
|
||||||
% dum(abs(dum)<1.e-12) = 0;
|
% dum(abs(dum)<1.e-12) = 0;
|
||||||
D2P(:,jcount) = dyn_vech(dum);
|
D2P(:,jcount) = dyn_vech(dum);
|
||||||
% D2P(:,:,j+offset,i) = dum;
|
% D2P(:,:,j+offset,i) = dum;
|
||||||
|
|
|
@ -187,7 +187,7 @@ end
|
||||||
%------------------------------------------------------------------------------
|
%------------------------------------------------------------------------------
|
||||||
|
|
||||||
% Compute the theoretical second order moments
|
% Compute the theoretical second order moments
|
||||||
tmp0 = lyapunov_symm(T,R*Q*R',DynareOptions.lyapunov_fixed_point_tol,DynareOptions.qz_criterium,DynareOptions.lyapunov_complex_threshold, [], [], DynareOptions.debug);
|
tmp0 = lyapunov_symm(T,R*Q*R',DynareOptions.lyapunov_fixed_point_tol,DynareOptions.qz_criterium,DynareOptions.lyapunov_complex_threshold, [], DynareOptions.debug);
|
||||||
mf = BayesInfo.mf1;
|
mf = BayesInfo.mf1;
|
||||||
|
|
||||||
% Get the non centered second order moments
|
% Get the non centered second order moments
|
||||||
|
|
|
@ -102,7 +102,7 @@ else
|
||||||
% return
|
% return
|
||||||
% end
|
% end
|
||||||
m = length(A);
|
m = length(A);
|
||||||
GAM = lyapunov_symm(A,B*M_.Sigma_e*B',options_.lyapunov_fixed_point_tol,options_.qz_criterium,options_.lyapunov_complex_threshold,1,[],options_.debug);
|
GAM = lyapunov_symm(A,B*M_.Sigma_e*B',options_.lyapunov_fixed_point_tol,options_.qz_criterium,options_.lyapunov_complex_threshold,1,options_.debug);
|
||||||
k = find(abs(GAM) < 1e-12);
|
k = find(abs(GAM) < 1e-12);
|
||||||
GAM(k) = 0;
|
GAM(k) = 0;
|
||||||
% if useautocorr,
|
% if useautocorr,
|
||||||
|
@ -116,7 +116,7 @@ else
|
||||||
% end
|
% end
|
||||||
% XX = lyapunov_symm_mr(A,BB,options_.qz_criterium,options_.lyapunov_complex_threshold,0);
|
% XX = lyapunov_symm_mr(A,BB,options_.qz_criterium,options_.lyapunov_complex_threshold,0);
|
||||||
for j=1:length(indexo),
|
for j=1:length(indexo),
|
||||||
dum = lyapunov_symm(A,dOm(:,:,j),options_.lyapunov_fixed_point_tol,options_.qz_criterium,options_.lyapunov_complex_threshold,2,[],options_.debug);
|
dum = lyapunov_symm(A,dOm(:,:,j),options_.lyapunov_fixed_point_tol,options_.qz_criterium,options_.lyapunov_complex_threshold,2,options_.debug);
|
||||||
% dum = XX(:,:,j);
|
% dum = XX(:,:,j);
|
||||||
k = find(abs(dum) < 1e-12);
|
k = find(abs(dum) < 1e-12);
|
||||||
dum(k) = 0;
|
dum(k) = 0;
|
||||||
|
@ -141,7 +141,7 @@ else
|
||||||
end
|
end
|
||||||
nexo = length(indexo);
|
nexo = length(indexo);
|
||||||
for j=1:length(indx),
|
for j=1:length(indx),
|
||||||
dum = lyapunov_symm(A,dA(:,:,j+nexo)*GAM*A'+A*GAM*dA(:,:,j+nexo)'+dOm(:,:,j+nexo),options_.lyapunov_fixed_point_tol,options_.qz_criterium,options_.lyapunov_complex_threshold,2,[],options_.debug);
|
dum = lyapunov_symm(A,dA(:,:,j+nexo)*GAM*A'+A*GAM*dA(:,:,j+nexo)'+dOm(:,:,j+nexo),options_.lyapunov_fixed_point_tol,options_.qz_criterium,options_.lyapunov_complex_threshold,2,options_.debug);
|
||||||
% dum = XX(:,:,j);
|
% dum = XX(:,:,j);
|
||||||
k = find(abs(dum) < 1e-12);
|
k = find(abs(dum) < 1e-12);
|
||||||
dum(k) = 0;
|
dum(k) = 0;
|
||||||
|
|
|
@ -45,7 +45,7 @@ n = length(i_var);
|
||||||
|
|
||||||
[A,B] = kalman_transition_matrix(dr,nstatic+(1:nspred),1:nc,M_.exo_nbr);
|
[A,B] = kalman_transition_matrix(dr,nstatic+(1:nspred),1:nc,M_.exo_nbr);
|
||||||
|
|
||||||
[vx,u] = lyapunov_symm(A,B*Sigma_e*B',options_.lyapunov_fixed_point_tol,options_.qz_criterium,options_.lyapunov_complex_threshold, [], [], options_.debug);
|
[vx,u] = lyapunov_symm(A,B*Sigma_e*B',options_.lyapunov_fixed_point_tol,options_.qz_criterium,options_.lyapunov_complex_threshold, [], options_.debug);
|
||||||
|
|
||||||
if size(u,2) > 0
|
if size(u,2) > 0
|
||||||
i_stat = find(any(abs(ghx*u) < options_.Schur_vec_tol,2)); %only set those variances of objective function for which variance is finite
|
i_stat = find(any(abs(ghx*u) < options_.Schur_vec_tol,2)); %only set those variances of objective function for which variance is finite
|
||||||
|
|
|
@ -0,0 +1,186 @@
|
||||||
|
function P=lyapunov_solver(T,R,Q,DynareOptions) % --*-- Unitary tests --*--
|
||||||
|
% function P=lyapunov_solver(T,R,Q,DynareOptions)
|
||||||
|
% Solves the Lyapunov equation P-T*P*T' = R*Q*R' arising in a state-space
|
||||||
|
% system, where P is the variance of the states
|
||||||
|
%
|
||||||
|
% Inputs
|
||||||
|
% T [double] n*n matrix.
|
||||||
|
% R [double] n*m matrix.
|
||||||
|
% Q [double] m*m matrix.
|
||||||
|
% DynareOptions [structure] Dynare options
|
||||||
|
%
|
||||||
|
% Outputs
|
||||||
|
% P [double] n*n matrix.
|
||||||
|
%
|
||||||
|
% Algorithms
|
||||||
|
% Default, if none of the other algorithms is selected:
|
||||||
|
% Reordered Schur decomposition (Bartels-Stewart algorithm)
|
||||||
|
% DynareOptions.lyapunov_fp == 1
|
||||||
|
% iteration-based fixed point algorithm
|
||||||
|
% DynareOptions.lyapunov_db == 1
|
||||||
|
% doubling algorithm
|
||||||
|
% DynareOptions.lyapunov_srs == 1
|
||||||
|
% Square-root solver for discrete-time Lyapunov equations (requires Matlab System Control toolbox
|
||||||
|
% or Octave control package)
|
||||||
|
|
||||||
|
% Copyright (C) 2016 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/>.
|
||||||
|
|
||||||
|
if DynareOptions.lyapunov_fp == 1
|
||||||
|
P = lyapunov_symm(T,R*Q'*R',DynareOptions.lyapunov_fixed_point_tol,DynareOptions.qz_criterium,DynareOptions.lyapunov_complex_threshold, 3, DynareOptions.debug);
|
||||||
|
elseif DynareOptions.lyapunov_db == 1
|
||||||
|
[P, errorflag] = disclyap_fast(T,R*Q*R',DynareOptions.lyapunov_doubling_tol);
|
||||||
|
if errorflag %use Schur-based method
|
||||||
|
P = lyapunov_symm(T,R*Q*R',DynareOptions.lyapunov_fixed_point_tol,DynareOptions.qz_criterium,DynareOptions.lyapunov_complex_threshold, [], DynareOptions.debug);
|
||||||
|
end
|
||||||
|
elseif DynareOptions.lyapunov_srs == 1
|
||||||
|
% works only with Matlab System Control toolbox or Octave control package,
|
||||||
|
if isoctave
|
||||||
|
if ~user_has_octave_forge_package('control')
|
||||||
|
error('lyapunov=square_root_solver not available; you must install the control package from Octave Forge')
|
||||||
|
end
|
||||||
|
else
|
||||||
|
if ~user_has_matlab_license('control_toolbox')
|
||||||
|
error('lyapunov=square_root_solver not available; you must install the control system toolbox')
|
||||||
|
end
|
||||||
|
end
|
||||||
|
chol_Q = R*chol(Q,'lower');
|
||||||
|
R_P = dlyapchol(T,chol_Q);
|
||||||
|
P = R_P' * R_P;
|
||||||
|
else
|
||||||
|
P = lyapunov_symm(T,R*Q*R',DynareOptions.lyapunov_fixed_point_tol,DynareOptions.qz_criterium,DynareOptions.lyapunov_complex_threshold, [], DynareOptions.debug);
|
||||||
|
end;
|
||||||
|
|
||||||
|
%@test:1
|
||||||
|
%$ t = NaN(10,1);
|
||||||
|
%$ options_.lyapunov_complex_threshold = 1e-15;
|
||||||
|
%$ options_.qz_zero_threshold = 1e-6;
|
||||||
|
%$ options_.qz_criterium=1-options_.qz_zero_threshold;
|
||||||
|
%$ options_.lyapunov_fixed_point_tol = 1e-10;
|
||||||
|
%$ options_.lyapunov_doubling_tol = 1e-16;
|
||||||
|
%$ options_.debug=0;
|
||||||
|
%$
|
||||||
|
%$ n_small=8;
|
||||||
|
%$ m_small=10;
|
||||||
|
%$ T_small=randn(n_small,n_small);
|
||||||
|
%$ T_small=0.99*T_small/max(abs(eigs(T_small)));
|
||||||
|
%$ tmp2=randn(m_small,m_small);
|
||||||
|
%$ Q_small=tmp2*tmp2';
|
||||||
|
%$ R_small=randn(n_small,m_small);
|
||||||
|
%$
|
||||||
|
%$ n_large=9;
|
||||||
|
%$ m_large=11;
|
||||||
|
%$ T_large=randn(n_large,n_large);
|
||||||
|
%$ T_large=0.99*T_large/max(abs(eigs(T_large)));
|
||||||
|
%$ tmp2=randn(m_large,m_large);
|
||||||
|
%$ Q_large=tmp2*tmp2';
|
||||||
|
%$ R_large=randn(n_large,m_large);
|
||||||
|
%$
|
||||||
|
%$ % DynareOptions.lyapunov_fp == 1
|
||||||
|
%$ options_.lyapunov_fp = 1;
|
||||||
|
%$ try
|
||||||
|
%$ Pstar1_small = lyapunov_solver(T_small,R_small,Q_small,options_);
|
||||||
|
%$ Pstar1_large = lyapunov_solver(T_large,R_large,Q_large,options_);
|
||||||
|
%$ t(1) = 1;
|
||||||
|
%$ catch
|
||||||
|
%$ t(1) = 0;
|
||||||
|
%$ end
|
||||||
|
%$
|
||||||
|
%$ % Dynareoptions.lyapunov_db == 1
|
||||||
|
%$ options_.lyapunov_fp = 0;
|
||||||
|
%$ options_.lyapunov_db = 1;
|
||||||
|
%$ try
|
||||||
|
%$ Pstar2_small = lyapunov_solver(T_small,R_small,Q_small,options_);
|
||||||
|
%$ Pstar2_large = lyapunov_solver(T_large,R_large,Q_large,options_);
|
||||||
|
%$ t(2) = 1;
|
||||||
|
%$ catch
|
||||||
|
%$ t(2) = 0;
|
||||||
|
%$ end
|
||||||
|
%$
|
||||||
|
%$ % Dynareoptions.lyapunov_srs == 1
|
||||||
|
%$ if (isoctave && user_has_octave_forge_package('control')) || (~isoctave && user_has_matlab_license('control_toolbox'))
|
||||||
|
%$ options_.lyapunov_db = 0;
|
||||||
|
%$ options_.lyapunov_srs = 1;
|
||||||
|
%$ try
|
||||||
|
%$ Pstar3_small = lyapunov_solver(T_small,R_small,Q_small,options_);
|
||||||
|
%$ Pstar3_large = lyapunov_solver(T_large,R_large,Q_large,options_);
|
||||||
|
%$ t(3) = 1;
|
||||||
|
%$ catch
|
||||||
|
%$ t(3) = 0;
|
||||||
|
%$ end
|
||||||
|
%$ else
|
||||||
|
%$ t(3) = 1;
|
||||||
|
%$ end
|
||||||
|
%$
|
||||||
|
%$ % Standard
|
||||||
|
%$ options_.lyapunov_srs = 0;
|
||||||
|
%$ try
|
||||||
|
%$ Pstar4_small = lyapunov_solver(T_small,R_small,Q_small,options_);
|
||||||
|
%$ Pstar4_large = lyapunov_solver(T_large,R_large,Q_large,options_);
|
||||||
|
%$ t(4) = 1;
|
||||||
|
%$ catch
|
||||||
|
%$ t(4) = 0;
|
||||||
|
%$ end
|
||||||
|
%$
|
||||||
|
%$ % Test the results.
|
||||||
|
%$
|
||||||
|
%$ if max(max(abs(Pstar1_small-Pstar2_small)))>1e-8
|
||||||
|
%$ t(5) = 0;
|
||||||
|
%$ else
|
||||||
|
%$ t(5) = 1;
|
||||||
|
%$ end
|
||||||
|
%$
|
||||||
|
%$ if (isoctave && user_has_octave_forge_package('control')) || (~isoctave && user_has_matlab_license('control_toolbox'))
|
||||||
|
%$ if max(max(abs(Pstar1_small-Pstar3_small)))>1e-8
|
||||||
|
%$ t(6) = 0;
|
||||||
|
%$ else
|
||||||
|
%$ t(6) = 1;
|
||||||
|
%$ end
|
||||||
|
%$ else
|
||||||
|
%$ t(6) = 1;
|
||||||
|
%$ end
|
||||||
|
%$
|
||||||
|
%$ if max(max(abs(Pstar1_small-Pstar4_small)))>1e-8
|
||||||
|
%$ t(7) = 0;
|
||||||
|
%$ else
|
||||||
|
%$ t(7) = 1;
|
||||||
|
%$ end
|
||||||
|
%$
|
||||||
|
%$ if max(max(abs(Pstar1_large-Pstar2_large)))>1e-8
|
||||||
|
%$ t(8) = 0;
|
||||||
|
%$ else
|
||||||
|
%$ t(8) = 1;
|
||||||
|
%$ end
|
||||||
|
%$
|
||||||
|
%$ if (isoctave && user_has_octave_forge_package('control')) || (~isoctave && user_has_matlab_license('control_toolbox'))
|
||||||
|
%$ if max(max(abs(Pstar1_large-Pstar3_large)))>1e-8
|
||||||
|
%$ t(9) = 0;
|
||||||
|
%$ else
|
||||||
|
%$ t(9) = 1;
|
||||||
|
%$ end
|
||||||
|
%$ else
|
||||||
|
%$ t(9) = 1;
|
||||||
|
%$ end
|
||||||
|
%$
|
||||||
|
%$ if max(max(abs(Pstar1_large-Pstar4_large)))>1e-8
|
||||||
|
%$ t(10) = 0;
|
||||||
|
%$ else
|
||||||
|
%$ t(10) = 1;
|
||||||
|
%$ end
|
||||||
|
%$
|
||||||
|
%$ T = all(t);
|
||||||
|
%@eof:1
|
|
@ -1,4 +1,4 @@
|
||||||
function [x,u] = lyapunov_symm(a,b,lyapunov_fixed_point_tol,qz_criterium,lyapunov_complex_threshold,method, R, debug) % --*-- Unitary tests --*--
|
function [x,u] = lyapunov_symm(a,b,lyapunov_fixed_point_tol,qz_criterium,lyapunov_complex_threshold,method,debug) % --*-- Unitary tests --*--
|
||||||
% Solves the Lyapunov equation x-a*x*a' = b, for b and x symmetric matrices.
|
% 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.
|
% If a has some unit roots, the function computes only the solution of the stable subsystem.
|
||||||
%
|
%
|
||||||
|
@ -20,12 +20,12 @@ function [x,u] = lyapunov_symm(a,b,lyapunov_fixed_point_tol,qz_criterium,lyapuno
|
||||||
%
|
%
|
||||||
% ALGORITHM
|
% ALGORITHM
|
||||||
% Uses reordered Schur decomposition (Bartels-Stewart algorithm)
|
% Uses reordered Schur decomposition (Bartels-Stewart algorithm)
|
||||||
% [method<3] or a fixed point algorithm (method==4)
|
% [method<3] or a fixed point algorithm (method==3)
|
||||||
%
|
%
|
||||||
% SPECIAL REQUIREMENTS
|
% SPECIAL REQUIREMENTS
|
||||||
% None
|
% None
|
||||||
|
|
||||||
% Copyright (C) 2006-2014 Dynare Team
|
% Copyright (C) 2006-2016 Dynare Team
|
||||||
%
|
%
|
||||||
% This file is part of Dynare.
|
% This file is part of Dynare.
|
||||||
%
|
%
|
||||||
|
@ -46,7 +46,7 @@ if nargin<6 || isempty(method)
|
||||||
method = 0;
|
method = 0;
|
||||||
end
|
end
|
||||||
|
|
||||||
if nargin<8
|
if nargin<7
|
||||||
debug = 0;
|
debug = 0;
|
||||||
end
|
end
|
||||||
|
|
||||||
|
@ -90,22 +90,7 @@ if method == 3
|
||||||
return;
|
return;
|
||||||
end;
|
end;
|
||||||
end;
|
end;
|
||||||
elseif method == 4
|
end
|
||||||
% works only with Matlab System Control toolbox or octave the control package,
|
|
||||||
if isoctave
|
|
||||||
if ~user_has_octave_forge_package('control')
|
|
||||||
error('lyapunov=square_root_solver not available; you must install the control package from Octave Forge')
|
|
||||||
end
|
|
||||||
else
|
|
||||||
if ~user_has_matlab_license('control_toolbox')
|
|
||||||
error('lyapunov=square_root_solver not available; you must install the control system toolbox')
|
|
||||||
end
|
|
||||||
end
|
|
||||||
chol_b = R*chol(b,'lower');
|
|
||||||
Rx = dlyapchol(a,chol_b);
|
|
||||||
x = Rx' * Rx;
|
|
||||||
return;
|
|
||||||
end;
|
|
||||||
|
|
||||||
if method
|
if method
|
||||||
persistent U T k n
|
persistent U T k n
|
||||||
|
@ -179,116 +164,3 @@ if i == 1
|
||||||
end
|
end
|
||||||
x = U(:,k+1:end)*x*U(:,k+1:end)';
|
x = U(:,k+1:end)*x*U(:,k+1:end)';
|
||||||
u = U(:,1:k);
|
u = U(:,1:k);
|
||||||
|
|
||||||
%@test:1
|
|
||||||
%$ t = NaN(10,1);
|
|
||||||
%$ lyapunov_complex_threshold = 1e-15;
|
|
||||||
%$ qz_zero_threshold = 1e-6;
|
|
||||||
%$ qz_criterium=1-qz_zero_threshold;
|
|
||||||
%$ lyapunov_fixed_point_tol = 1e-10;
|
|
||||||
%$ lyapunov_doubling_tol = 1e-16;
|
|
||||||
%$
|
|
||||||
%$ n_small=8;
|
|
||||||
%$ m_small=10;
|
|
||||||
%$ T_small=randn(n_small,n_small);
|
|
||||||
%$ T_small=0.99*T_small/max(abs(eigs(T_small)));
|
|
||||||
%$ tmp2=randn(m_small,m_small);
|
|
||||||
%$ Q_small=tmp2*tmp2';
|
|
||||||
%$ R_small=randn(n_small,m_small);
|
|
||||||
%$
|
|
||||||
%$ n_large=9;
|
|
||||||
%$ m_large=11;
|
|
||||||
%$ T_large=randn(n_large,n_large);
|
|
||||||
%$ T_large=0.99*T_large/max(abs(eigs(T_large)));
|
|
||||||
%$ tmp2=randn(m_large,m_large);
|
|
||||||
%$ Q_large=tmp2*tmp2';
|
|
||||||
%$ R_large=randn(n_large,m_large);
|
|
||||||
%$
|
|
||||||
%$ % DynareOptions.lyapunov_fp == 1
|
|
||||||
%$ try
|
|
||||||
%$ Pstar1_small = lyapunov_symm(T_small,R_small*Q_small*R_small',lyapunov_fixed_point_tol,qz_criterium,lyapunov_complex_threshold,3, [], 0);
|
|
||||||
%$ Pstar1_large = lyapunov_symm(T_large,R_large*Q_large*R_large',lyapunov_fixed_point_tol,qz_criterium,lyapunov_complex_threshold,3, [], 0);
|
|
||||||
%$ t(1) = 1;
|
|
||||||
%$ catch
|
|
||||||
%$ t(1) = 0;
|
|
||||||
%$ end
|
|
||||||
%$
|
|
||||||
%$ % Dynareoptions.lyapunov_db == 1
|
|
||||||
%$ try
|
|
||||||
%$ Pstar2_small = disclyap_fast(T_small,R_small*Q_small*R_small',lyapunov_doubling_tol);
|
|
||||||
%$ Pstar2_large = disclyap_fast(T_large,R_large*Q_large*R_large',lyapunov_doubling_tol);
|
|
||||||
%$ t(2) = 1;
|
|
||||||
%$ catch
|
|
||||||
%$ t(2) = 0;
|
|
||||||
%$ end
|
|
||||||
%$
|
|
||||||
%$ % Dynareoptions.lyapunov_srs == 1
|
|
||||||
%$ if (isoctave && user_has_octave_forge_package('control')) || (~isoctave && user_has_matlab_license('control_toolbox'))
|
|
||||||
%$ try
|
|
||||||
%$ Pstar3_small = lyapunov_symm(T_small,Q_small,lyapunov_fixed_point_tol,qz_criterium,lyapunov_complex_threshold,4,R_small,0);
|
|
||||||
%$ Pstar3_large = lyapunov_symm(T_large,Q_large,lyapunov_fixed_point_tol,qz_criterium,lyapunov_complex_threshold,4,R_large,0);
|
|
||||||
%$ t(3) = 1;
|
|
||||||
%$ catch
|
|
||||||
%$ t(3) = 0;
|
|
||||||
%$ end
|
|
||||||
%$ else
|
|
||||||
%$ t(3) = 1;
|
|
||||||
%$ end
|
|
||||||
%$
|
|
||||||
%$ % Standard
|
|
||||||
%$ try
|
|
||||||
%$ Pstar4_small = lyapunov_symm(T_small,R_small*Q_small*R_small',lyapunov_fixed_point_tol,qz_criterium, lyapunov_complex_threshold, [], [], 0);
|
|
||||||
%$ Pstar4_large = lyapunov_symm(T_large,R_large*Q_large*R_large',lyapunov_fixed_point_tol,qz_criterium, lyapunov_complex_threshold, [], [], 0);
|
|
||||||
%$ t(4) = 1;
|
|
||||||
%$ catch
|
|
||||||
%$ t(4) = 0;
|
|
||||||
%$ end
|
|
||||||
%$
|
|
||||||
%$ % Test the results.
|
|
||||||
%$
|
|
||||||
%$ if max(max(abs(Pstar1_small-Pstar2_small)))>1e-8
|
|
||||||
%$ t(5) = 0;
|
|
||||||
%$ else
|
|
||||||
%$ t(5) = 1;
|
|
||||||
%$ end
|
|
||||||
%$
|
|
||||||
%$ if (isoctave && user_has_octave_forge_package('control')) || (~isoctave && user_has_matlab_license('control_toolbox'))
|
|
||||||
%$ if max(max(abs(Pstar1_small-Pstar3_small)))>1e-8
|
|
||||||
%$ t(6) = 0;
|
|
||||||
%$ else
|
|
||||||
%$ t(6) = 1;
|
|
||||||
%$ end
|
|
||||||
%$ else
|
|
||||||
%$ t(6) = 1;
|
|
||||||
%$ end
|
|
||||||
%$
|
|
||||||
%$ if max(max(abs(Pstar1_small-Pstar4_small)))>1e-8
|
|
||||||
%$ t(7) = 0;
|
|
||||||
%$ else
|
|
||||||
%$ t(7) = 1;
|
|
||||||
%$ end
|
|
||||||
%$
|
|
||||||
%$ if max(max(abs(Pstar1_large-Pstar2_large)))>1e-8
|
|
||||||
%$ t(8) = 0;
|
|
||||||
%$ else
|
|
||||||
%$ t(8) = 1;
|
|
||||||
%$ end
|
|
||||||
%$
|
|
||||||
%$ if (isoctave && user_has_octave_forge_package('control')) || (~isoctave && user_has_matlab_license('control_toolbox'))
|
|
||||||
%$ if max(max(abs(Pstar1_large-Pstar3_large)))>1e-8
|
|
||||||
%$ t(9) = 0;
|
|
||||||
%$ else
|
|
||||||
%$ t(9) = 1;
|
|
||||||
%$ end
|
|
||||||
%$ else
|
|
||||||
%$ t(9) = 1;
|
|
||||||
%$ end
|
|
||||||
%$
|
|
||||||
%$ if max(max(abs(Pstar1_large-Pstar4_large)))>1e-8
|
|
||||||
%$ t(10) = 0;
|
|
||||||
%$ else
|
|
||||||
%$ t(10) = 1;
|
|
||||||
%$ end
|
|
||||||
%$
|
|
||||||
%$ T = all(t);
|
|
||||||
%@eof:1
|
|
|
@ -306,7 +306,7 @@ ReducedForm.mf1 = mf1;
|
||||||
switch DynareOptions.particle.initialization
|
switch DynareOptions.particle.initialization
|
||||||
case 1% Initial state vector covariance is the ergodic variance associated to the first order Taylor-approximation of the model.
|
case 1% Initial state vector covariance is the ergodic variance associated to the first order Taylor-approximation of the model.
|
||||||
StateVectorMean = ReducedForm.constant(mf0);
|
StateVectorMean = ReducedForm.constant(mf0);
|
||||||
StateVectorVariance = lyapunov_symm(ReducedForm.ghx(mf0,:),ReducedForm.ghu(mf0,:)*ReducedForm.Q*ReducedForm.ghu(mf0,:)',DynareOptions.lyapunov_fixed_point_tol,DynareOptions.qz_criterium,DynareOptions.lyapunov_complex_threshold,[],[],DynareOptions.debug);
|
StateVectorVariance = lyapunov_symm(ReducedForm.ghx(mf0,:),ReducedForm.ghu(mf0,:)*ReducedForm.Q*ReducedForm.ghu(mf0,:)',DynareOptions.lyapunov_fixed_point_tol,DynareOptions.qz_criterium,DynareOptions.lyapunov_complex_threshold,[],DynareOptions.debug);
|
||||||
case 2% Initial state vector covariance is a monte-carlo based estimate of the ergodic variance (consistent with a k-order Taylor-approximation of the model).
|
case 2% Initial state vector covariance is a monte-carlo based estimate of the ergodic variance (consistent with a k-order Taylor-approximation of the model).
|
||||||
StateVectorMean = ReducedForm.constant(mf0);
|
StateVectorMean = ReducedForm.constant(mf0);
|
||||||
old_DynareOptionsperiods = DynareOptions.periods;
|
old_DynareOptionsperiods = DynareOptions.periods;
|
||||||
|
|
|
@ -1,10 +1,10 @@
|
||||||
function [X,exitflag]=disclyap_fast(G,V,tol,ch)
|
function [X,exitflag]=disclyap_fast(G,V,tol,check_flag)
|
||||||
% function X=disclyap_fast(G,V,ch)
|
% function X=disclyap_fast(G,V,ch)
|
||||||
% Inputs:
|
% Inputs:
|
||||||
% - G [double] first input matrix
|
% - G [double] first input matrix
|
||||||
% - V [double] second input matrix
|
% - V [double] second input matrix
|
||||||
% - tol [scalar] tolerance criterion
|
% - tol [scalar] tolerance criterion
|
||||||
% - ch empty of non-empty if non-empty: check positive-definiteness
|
% - check_flag empty of non-empty if non-empty: check positive-definiteness
|
||||||
% Outputs:
|
% Outputs:
|
||||||
% - X [double] solution matrix
|
% - X [double] solution matrix
|
||||||
% - exitflag [scalar] 0 if solution is found, 1 otherwise
|
% - exitflag [scalar] 0 if solution is found, 1 otherwise
|
||||||
|
@ -13,13 +13,13 @@ function [X,exitflag]=disclyap_fast(G,V,tol,ch)
|
||||||
% X=G*X*G'+V
|
% X=G*X*G'+V
|
||||||
% Using the Doubling Algorithm
|
% Using the Doubling Algorithm
|
||||||
%
|
%
|
||||||
% If ch is defined then the code will check if the resulting X
|
% If check_flag is defined then the code will check if the resulting X
|
||||||
% is positive definite and generate an error message if it is not
|
% is positive definite and generate an error message if it is not
|
||||||
%
|
%
|
||||||
% Joe Pearlman and Alejandro Justiniano
|
% Joe Pearlman and Alejandro Justiniano
|
||||||
% 3/5/2005
|
% 3/5/2005
|
||||||
|
|
||||||
% Copyright (C) 2010-2015 Dynare Team
|
% Copyright (C) 2010-2016 Dynare Team
|
||||||
%
|
%
|
||||||
% This file is part of Dynare.
|
% This file is part of Dynare.
|
||||||
%
|
%
|
||||||
|
@ -36,16 +36,13 @@ function [X,exitflag]=disclyap_fast(G,V,tol,ch)
|
||||||
% You should have received a copy of the GNU General Public License
|
% You should have received a copy of the GNU General Public License
|
||||||
% along with Dynare. If not, see <http://www.gnu.org/licenses/>.
|
% along with Dynare. If not, see <http://www.gnu.org/licenses/>.
|
||||||
|
|
||||||
if nargin <= 3 || isempty( ch ) == 1
|
if nargin <= 3 || isempty( check_flag ) == 1
|
||||||
flag_ch = 0;
|
flag_ch = 0;
|
||||||
else
|
else
|
||||||
flag_ch = 1;
|
flag_ch = 1;
|
||||||
end
|
end
|
||||||
s=size(G,1);
|
|
||||||
exitflag=0;
|
exitflag=0;
|
||||||
|
|
||||||
%tol = 1e-16;
|
|
||||||
|
|
||||||
P0=V;
|
P0=V;
|
||||||
A0=G;
|
A0=G;
|
||||||
|
|
||||||
|
|
|
@ -129,7 +129,7 @@ end;
|
||||||
% and compute 2nd order mean correction on stationary variables (in case of
|
% and compute 2nd order mean correction on stationary variables (in case of
|
||||||
% HP filtering, this mean correction is computed *before* filtering)
|
% HP filtering, this mean correction is computed *before* filtering)
|
||||||
if options_.order == 2 || options_.hp_filter == 0
|
if options_.order == 2 || options_.hp_filter == 0
|
||||||
[vx, u] = lyapunov_symm(A,B*M_.Sigma_e*B',options_.lyapunov_fixed_point_tol,options_.qz_criterium,options_.lyapunov_complex_threshold,[],[],options_.debug);
|
[vx, u] = lyapunov_symm(A,B*M_.Sigma_e*B',options_.lyapunov_fixed_point_tol,options_.qz_criterium,options_.lyapunov_complex_threshold,[],options_.debug);
|
||||||
if options_.block == 0
|
if options_.block == 0
|
||||||
iky = inv_order_var(ivar);
|
iky = inv_order_var(ivar);
|
||||||
else
|
else
|
||||||
|
@ -188,11 +188,11 @@ if options_.hp_filter == 0 && ~options_.bandpass.indicator
|
||||||
b1 = b1*cs;
|
b1 = b1*cs;
|
||||||
b2(:,exo_names_orig_ord) = ghu(iky,:);
|
b2(:,exo_names_orig_ord) = ghu(iky,:);
|
||||||
b2 = b2*cs;
|
b2 = b2*cs;
|
||||||
vx = lyapunov_symm(A,b1*b1',options_.lyapunov_fixed_point_tol,options_.qz_criterium,options_.lyapunov_complex_threshold,1,[],options_.debug);
|
vx = lyapunov_symm(A,b1*b1',options_.lyapunov_fixed_point_tol,options_.qz_criterium,options_.lyapunov_complex_threshold,1,options_.debug);
|
||||||
vv = diag(aa*vx*aa'+b2*b2');
|
vv = diag(aa*vx*aa'+b2*b2');
|
||||||
vv2 = 0;
|
vv2 = 0;
|
||||||
for i=1:M_.exo_nbr
|
for i=1:M_.exo_nbr
|
||||||
vx1 = lyapunov_symm(A,b1(:,i)*b1(:,i)',options_.lyapunov_fixed_point_tol,options_.qz_criterium,options_.lyapunov_complex_threshold,2,[],options_.debug);
|
vx1 = lyapunov_symm(A,b1(:,i)*b1(:,i)',options_.lyapunov_fixed_point_tol,options_.qz_criterium,options_.lyapunov_complex_threshold,2,options_.debug);
|
||||||
vx2 = abs(diag(aa*vx1*aa'+b2(:,i)*b2(:,i)'));
|
vx2 = abs(diag(aa*vx1*aa'+b2(:,i)*b2(:,i)'));
|
||||||
Gamma_y{nar+2}(stationary_vars,i) = vx2;
|
Gamma_y{nar+2}(stationary_vars,i) = vx2;
|
||||||
vv2 = vv2 +vx2;
|
vv2 = vv2 +vx2;
|
||||||
|
|
|
@ -55,7 +55,7 @@ elseif flagmoments==-1
|
||||||
tau=[oo_.dr.ys(oo_.dr.order_var); g1(:)];
|
tau=[oo_.dr.ys(oo_.dr.order_var); g1(:)];
|
||||||
|
|
||||||
else
|
else
|
||||||
GAM = lyapunov_symm(A,B*M_.Sigma_e*B',options_.lyapunov_fixed_point_tol,options_.qz_criterium,options_.lyapunov_complex_threshold,[],[],options_.debug);
|
GAM = lyapunov_symm(A,B*M_.Sigma_e*B',options_.lyapunov_fixed_point_tol,options_.qz_criterium,options_.lyapunov_complex_threshold,[],options_.debug);
|
||||||
k = find(abs(GAM) < 1e-12);
|
k = find(abs(GAM) < 1e-12);
|
||||||
GAM(k) = 0;
|
GAM(k) = 0;
|
||||||
if useautocorr,
|
if useautocorr,
|
||||||
|
|
Loading…
Reference in New Issue