Skip to content
GitLab
Commit 94d215d7 authored by ratto's avatar ratto
Browse files

Modifications linked to previous commit: 

Added first order moments 
Added LRE analysis for trivial no-identification

git-svn-id: https://www.dynare.org/svn/dynare/trunk@3361 ac1d8469-bf42-47a9-8791-bf33cf982152
parent cda0f571
Hide whitespace changes
Inline Side-by-side
matlab/getH.m
View file @ 94d215d7
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
%
......
matlab/getJJ.m
View file @ 94d215d7
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
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;
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
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;
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
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
gam = vech(GAM{1});
for j=1:nlags,
gam = [gam; vec(GAM{j+1})];
end
gam = [oo_.dr.ys(oo_.dr.order_var(mf)); gam];
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
matlab/identification_checks.m
View file @ 94d215d7
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)<npar
ixno = 0;
% - find out which parameters are involved,
......@@ -106,12 +121,16 @@ end
McoH = NaN(npar,1);
McoJ = NaN(npar,1);
McoGP = NaN(npar0,1);
for ii = 1:size(H1,2);
McoH(ind1(ii),:) = [cosn([H1(:,ii),H1(:,find([1:1:size(H1,2)]~=ii))])];
end
for ii = 1:size(JJ1,2);
McoJ(ind2(ii),:) = [cosn([JJ1(:,ii),JJ1(:,find([1:1:size(JJ1,2)]~=ii))])];
end
for ii = 1:size(gp1,2);
McoGP(ind3(ii),:) = [cosn([gp1(:,ii),gp1(:,find([1:1:size(gp1,2)]~=ii))])];
end
% format long % some are nearly 1
% McoJ
......@@ -123,6 +142,7 @@ end
PcoH = NaN(npar,npar);
PcoJ = NaN(npar,npar);
PcoGP = NaN(npar0,npar0);
for ii = 1:size(H1,2);
for jj = 1:size(H1,2);
PcoH(ind1(ii),ind1(jj)) = [cosn([H1(:,ii),H1(:,jj)])];
......@@ -135,6 +155,13 @@ for ii = 1:size(JJ1,2);
end
end
for ii = 1:size(gp1,2);
for jj = 1:size(gp1,2);
PcoGP(ind3(ii),ind3(jj)) = [cosn([gp1(:,ii),gp1(:,jj)])];
end
end
ind01 = zeros(npar,1);
ind02 = zeros(npar,1);
ind01(ind1) = 1;
......
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment