Merge pull request #983 from JohannesPfeifer/lyap_fix
Two fixes regarding solution of lyapunov equationstime-shift
commit
13cb9930f4
|
@ -128,13 +128,13 @@ if options_.lik_init == 1 % Kalman filter
|
|||
kalman_algo = 1;
|
||||
end
|
||||
if options_.lyapunov_fp == 1
|
||||
Pstar = lyapunov_symm(T,R*Q*R',options_.lyapunov_fixed_point_tol,options_.lyapunov_complex_threshold, 3, [], options_.debug);
|
||||
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_.lyapunov_complex_threshold, 4, R, options_.debug);
|
||||
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_.qz_criterium,options_.lyapunov_complex_threshold, [], [], options_.debug);
|
||||
Pstar = lyapunov_symm(T,R*Q*R',options_.lyapunov_fixed_point_tol,options_.qz_criterium,options_.lyapunov_complex_threshold, [], [], options_.debug);
|
||||
end;
|
||||
Pinf = [];
|
||||
elseif options_.lik_init == 2 % Old Diffuse Kalman filter
|
||||
|
@ -171,13 +171,13 @@ elseif options_.lik_init == 5 % Old diffuse Kalman filter only for th
|
|||
R_tmp = R(stable, :);
|
||||
T_tmp = T(stable,stable);
|
||||
if options_.lyapunov_fp == 1
|
||||
Pstar_tmp = lyapunov_symm(T_tmp,R_tmp*Q*R_tmp',options_.lyapunov_fixed_point_tol,options_.lyapunov_complex_threshold, 3, [], options_.debug);
|
||||
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_.lyapunov_complex_threshold, 4, R_tmp, options_.debug);
|
||||
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_.qz_criterium,options_.lyapunov_complex_threshold, [], [], options_.debug);
|
||||
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;
|
||||
Pinf = [];
|
||||
|
|
|
@ -93,7 +93,7 @@ for i=M_.maximum_lag:-1:2
|
|||
end
|
||||
Gamma = zeros(nvar,cutoff+1);
|
||||
[A,B] = kalman_transition_matrix(dr,ikx',1:nx,M_.exo_nbr);
|
||||
[vx, u] = lyapunov_symm(A,B*M_.Sigma_e*B',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);
|
||||
if ~isempty(u)
|
||||
iky = iky(find(any(abs(ghx(iky,:)*u) < options_.Schur_vec_tol,2)));
|
||||
|
|
|
@ -370,13 +370,13 @@ switch DynareOptions.lik_init
|
|||
kalman_algo = 1;
|
||||
end
|
||||
if DynareOptions.lyapunov_fp == 1
|
||||
Pstar = lyapunov_symm(T,R*Q'*R',DynareOptions.lyapunov_fixed_point_tol,DynareOptions.lyapunov_complex_threshold, 3, [], DynareOptions.debug);
|
||||
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.lyapunov_complex_threshold, 4, R, DynareOptions.debug);
|
||||
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.qz_criterium,DynareOptions.lyapunov_complex_threshold, [], [], DynareOptions.debug);
|
||||
Pstar = lyapunov_symm(T,R*Q*R',DynareOptions.lyapunov_fixed_point_tol,DynareOptions.qz_criterium,DynareOptions.lyapunov_complex_threshold, [], [], DynareOptions.debug);
|
||||
end;
|
||||
Pinf = [];
|
||||
a = zeros(mm,1);
|
||||
|
@ -472,7 +472,7 @@ switch DynareOptions.lik_init
|
|||
if err
|
||||
disp(['dsge_likelihood:: I am not able to solve the Riccati equation, so I switch to lik_init=1!']);
|
||||
DynareOptions.lik_init = 1;
|
||||
Pstar = lyapunov_symm(T,R*Q*R',DynareOptions.qz_criterium,DynareOptions.lyapunov_complex_threshold, [], [], DynareOptions.debug);
|
||||
Pstar = lyapunov_symm(T,R*Q*R',DynareOptions.lyapunov_fixed_point_tol,DynareOptions.qz_criterium,DynareOptions.lyapunov_complex_threshold, [], [], DynareOptions.debug);
|
||||
end
|
||||
Pinf = [];
|
||||
a = zeros(mm,1);
|
||||
|
@ -493,13 +493,13 @@ switch DynareOptions.lik_init
|
|||
R_tmp = R(stable, :);
|
||||
T_tmp = T(stable,stable);
|
||||
if DynareOptions.lyapunov_fp == 1
|
||||
Pstar_tmp = lyapunov_symm(T_tmp,R_tmp*Q*R_tmp',DynareOptions.lyapunov_fixed_point_tol,DynareOptions.lyapunov_complex_threshold, 3, [], DynareOptions.debug);
|
||||
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.lyapunov_complex_threshold, 4, R_tmp, DynareOptions.debug);
|
||||
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.lyapunov_complex_threshold, [], [], DynareOptions.debug);
|
||||
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;
|
||||
Pinf = [];
|
||||
|
@ -579,14 +579,14 @@ if analytic_derivation,
|
|||
for i=1:EstimatedParameters.nvx
|
||||
k =EstimatedParameters.var_exo(i,1);
|
||||
DQ(k,k,i) = 2*sqrt(Q(k,k));
|
||||
dum = lyapunov_symm(T,DOm(:,:,i),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);
|
||||
% dum(kk) = 0;
|
||||
DP(:,:,i)=dum;
|
||||
if full_Hess
|
||||
for j=1:i,
|
||||
jcount=jcount+1;
|
||||
dum = lyapunov_symm(T,dyn_unvech(D2Om(:,jcount)),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);
|
||||
% dum(kk) = 0;
|
||||
D2P(:,jcount)=dyn_vech(dum);
|
||||
|
@ -606,7 +606,7 @@ if analytic_derivation,
|
|||
offset = offset + EstimatedParameters.nvn;
|
||||
if DynareOptions.lik_init==1,
|
||||
for j=1:EstimatedParameters.np
|
||||
dum = lyapunov_symm(T,DT(:,:,j+offset)*Pstar*T'+T*Pstar*DT(:,:,j+offset)'+DOm(:,:,j+offset),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);
|
||||
% dum(kk) = 0;
|
||||
DP(:,:,j+offset)=dum;
|
||||
|
@ -620,7 +620,7 @@ if analytic_derivation,
|
|||
D2Tij = reshape(D2T(:,jcount),size(T));
|
||||
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;
|
||||
dum = lyapunov_symm(T,tmp,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;
|
||||
D2P(:,jcount) = dyn_vech(dum);
|
||||
% D2P(:,:,j+offset,i) = dum;
|
||||
|
|
|
@ -162,7 +162,7 @@ end
|
|||
%------------------------------------------------------------------------------
|
||||
|
||||
% Compute the theoretical second order moments
|
||||
tmp0 = lyapunov_symm(T,R*Q*R',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;
|
||||
|
||||
% Get the non centered second order moments
|
||||
|
|
|
@ -64,7 +64,7 @@ else
|
|||
% return
|
||||
% end
|
||||
m = length(A);
|
||||
GAM = lyapunov_symm(A,B*M_.Sigma_e*B',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);
|
||||
GAM(k) = 0;
|
||||
% if useautocorr,
|
||||
|
@ -78,7 +78,7 @@ else
|
|||
% 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,[],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);
|
||||
k = find(abs(dum) < 1e-12);
|
||||
dum(k) = 0;
|
||||
|
@ -103,7 +103,7 @@ else
|
|||
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,[],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);
|
||||
k = find(abs(dum) < 1e-12);
|
||||
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);
|
||||
|
||||
[vx,u] = lyapunov_symm(A,B*Sigma_e*B',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
|
||||
i_stat = find(any(abs(ghx*u) < options_.Schur_vec_tol,2)); %only set those variances of objective function for which variance is finite
|
||||
|
|
|
@ -1,14 +1,12 @@
|
|||
function [x,u] = lyapunov_symm(a,b,third_argument,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, R, debug) % --*-- Unitary tests --*--
|
||||
% 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.
|
||||
% third_argument [double] scalar, if method <= 2 :
|
||||
% qz_criterium = third_argument unit root threshold for eigenvalues in a,
|
||||
% elseif method = 3 :
|
||||
% tol =third_argument the convergence criteria for fixed_point algorithm.
|
||||
% qz_criterium [double] unit root threshold for eigenvalues
|
||||
% lyapunov_fixed_point_tol [double] convergence criteria for fixed_point algorithm.
|
||||
% 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
|
||||
|
@ -44,11 +42,11 @@ function [x,u] = lyapunov_symm(a,b,third_argument,lyapunov_complex_threshold,met
|
|||
% 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<5 || isempty(method)
|
||||
if nargin<6 || isempty(method)
|
||||
method = 0;
|
||||
end
|
||||
|
||||
if nargin<7
|
||||
if nargin<8
|
||||
debug = 0;
|
||||
end
|
||||
|
||||
|
@ -57,7 +55,6 @@ if method == 3
|
|||
if ~isempty(method1)
|
||||
method = method1;
|
||||
end;
|
||||
tol = third_argument;
|
||||
if debug
|
||||
fprintf('lyapunov_symm:: [method=%d] \n',method);
|
||||
end
|
||||
|
@ -73,7 +70,7 @@ if method == 3
|
|||
end;
|
||||
at = a';
|
||||
%fixed point iterations
|
||||
while evol > tol && it_fp < max_it_fp;
|
||||
while evol > lyapunov_fixed_point_tol && it_fp < max_it_fp;
|
||||
X_old = X;
|
||||
X = a * X * at + b;
|
||||
evol = max(sum(abs(X - X_old))); %norm_1
|
||||
|
@ -110,7 +107,6 @@ elseif method == 4
|
|||
return;
|
||||
end;
|
||||
|
||||
qz_criterium = third_argument;
|
||||
if method
|
||||
persistent U T k n
|
||||
else
|
||||
|
@ -210,8 +206,8 @@ u = U(:,1:k);
|
|||
%$
|
||||
%$ % DynareOptions.lyapunov_fp == 1
|
||||
%$ try
|
||||
%$ Pstar1_small = lyapunov_symm(T_small,R_small*Q_small*R_small',lyapunov_fixed_point_tol,lyapunov_fixed_point_tol,3, [], 0);
|
||||
%$ Pstar1_large = lyapunov_symm(T_large,R_large*Q_large*R_large',lyapunov_fixed_point_tol,lyapunov_fixed_point_tol,3, [], 0);
|
||||
%$ 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;
|
||||
|
@ -229,8 +225,8 @@ u = U(:,1:k);
|
|||
%$ % 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,lyapunov_complex_threshold,4,R_small,0);
|
||||
%$ Pstar3_large = lyapunov_symm(T_large,Q_large,lyapunov_fixed_point_tol,lyapunov_complex_threshold,4,R_large,0);
|
||||
%$ 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;
|
||||
|
@ -241,8 +237,8 @@ u = U(:,1:k);
|
|||
%$
|
||||
%$ % Standard
|
||||
%$ try
|
||||
%$ Pstar4_small = lyapunov_symm(T_small,R_small*Q_small*R_small',qz_criterium, lyapunov_complex_threshold, [], [], 0);
|
||||
%$ Pstar4_large = lyapunov_symm(T_large,R_large*Q_large*R_large',qz_criterium, lyapunov_complex_threshold, [], [], 0);
|
||||
%$ 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;
|
||||
|
|
|
@ -271,7 +271,7 @@ ReducedForm.mf1 = mf1;
|
|||
switch DynareOptions.particle.initialization
|
||||
case 1% Initial state vector covariance is the ergodic variance associated to the first order Taylor-approximation of the model.
|
||||
StateVectorMean = ReducedForm.constant(mf0);
|
||||
StateVectorVariance = lyapunov_symm(ReducedForm.ghx(mf0,:),ReducedForm.ghu(mf0,:)*ReducedForm.Q*ReducedForm.ghu(mf0,:)',1e-12,1e-12,[],[],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).
|
||||
StateVectorMean = ReducedForm.constant(mf0);
|
||||
old_DynareOptionsperiods = DynareOptions.periods;
|
||||
|
|
|
@ -1,6 +1,14 @@
|
|||
function X=disclyap_fast(G,V,tol,ch)
|
||||
function [X,exitflag]=disclyap_fast(G,V,tol,ch)
|
||||
% function X=disclyap_fast(G,V,ch)
|
||||
%
|
||||
% Inputs:
|
||||
% - G [double] first input matrix
|
||||
% - V [double] second input matrix
|
||||
% - tol [scalar] tolerance criterion
|
||||
% - ch empty of non-empty if non-empty: check positive-definiteness
|
||||
% Outputs:
|
||||
% - X [double] solution matrix
|
||||
% - exitflag [scalar] 0 if solution is found, 1 otherwise
|
||||
%
|
||||
% Solve the discrete Lyapunov Equation
|
||||
% X=G*X*G'+V
|
||||
% Using the Doubling Algorithm
|
||||
|
@ -11,7 +19,7 @@ function X=disclyap_fast(G,V,tol,ch)
|
|||
% Joe Pearlman and Alejandro Justiniano
|
||||
% 3/5/2005
|
||||
|
||||
% Copyright (C) 2010-2012 Dynare Team
|
||||
% Copyright (C) 2010-2015 Dynare Team
|
||||
%
|
||||
% This file is part of Dynare.
|
||||
%
|
||||
|
@ -34,6 +42,7 @@ else
|
|||
flag_ch = 1;
|
||||
end
|
||||
s=size(G,1);
|
||||
exitflag=0;
|
||||
|
||||
%tol = 1e-16;
|
||||
|
||||
|
@ -41,13 +50,20 @@ P0=V;
|
|||
A0=G;
|
||||
|
||||
matd=1;
|
||||
while matd > tol
|
||||
iter=1;
|
||||
while matd > tol && iter< 2000
|
||||
P1=P0+A0*P0*A0';
|
||||
A1=A0*A0;
|
||||
matd=max( max( abs( P1 - P0 ) ) );
|
||||
P0=P1;
|
||||
A0=A1;
|
||||
iter=iter+1;
|
||||
end
|
||||
if iter==5000
|
||||
X=NaN(P0);
|
||||
exitflag=1;
|
||||
return
|
||||
end
|
||||
clear A0 A1 P1;
|
||||
|
||||
X=(P0+P0')/2;
|
||||
|
@ -56,6 +72,7 @@ X=(P0+P0')/2;
|
|||
if flag_ch==1
|
||||
[C,p]=chol(X);
|
||||
if p ~= 0
|
||||
exitflag=1;
|
||||
error('X is not positive definite')
|
||||
end
|
||||
end
|
|
@ -129,7 +129,7 @@ end;
|
|||
% and compute 2nd order mean correction on stationary variables (in case of
|
||||
% HP filtering, this mean correction is computed *before* filtering)
|
||||
if options_.order == 2 || options_.hp_filter == 0
|
||||
[vx, u] = lyapunov_symm(A,B*M_.Sigma_e*B',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
|
||||
iky = inv_order_var(ivar);
|
||||
else
|
||||
|
@ -188,11 +188,11 @@ if options_.hp_filter == 0
|
|||
b1 = b1*cs;
|
||||
b2(:,exo_names_orig_ord) = ghu(iky,:);
|
||||
b2 = b2*cs;
|
||||
vx = lyapunov_symm(A,b1*b1',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');
|
||||
vv2 = 0;
|
||||
for i=1:M_.exo_nbr
|
||||
vx1 = lyapunov_symm(A,b1(:,i)*b1(:,i)',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)'));
|
||||
Gamma_y{nar+2}(stationary_vars,i) = vx2;
|
||||
vv2 = vv2 +vx2;
|
||||
|
|
|
@ -51,7 +51,7 @@ elseif flagmoments==-1
|
|||
tau=[oo_.dr.ys(oo_.dr.order_var); g1(:)];
|
||||
|
||||
else
|
||||
GAM = lyapunov_symm(A,B*M_.Sigma_e*B',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);
|
||||
GAM(k) = 0;
|
||||
if useautocorr,
|
||||
|
|
|
@ -60,7 +60,7 @@ V0 = oo_.var(k4,k4);
|
|||
Atest = [dr2.M1(k3,:) dr2.M2(k3,:) dr2.M4(k3,:); eye(6) zeros(6,10);zeros(4,16)];
|
||||
Btest = [dr2.M3(k3,:); zeros(6,4); eye(4)];
|
||||
|
||||
V1=lyapunov_symm(Atest,Btest*M_.Sigma_e*Btest',1+1e-6,options_.lyapunov_complex_threshold);
|
||||
V1=lyapunov_symm(Atest,Btest*M_.Sigma_e*Btest',options_.lyapunov_fixed_point_tol,1+1e-6,options_.lyapunov_complex_threshold);
|
||||
|
||||
if max(max(abs(V1(1:6,1:6)-V0)))
|
||||
disp('Test OK')
|
||||
|
|
Loading…
Reference in New Issue