New identification routines.
git-svn-id: https://www.dynare.org/svn/dynare/trunk@2963 ac1d8469-bf42-47a9-8791-bf33cf982152time-shift
ratto
parent
e916b1380f
commit
df921ce829
|
|
@ -0,0 +1,20 @@
|
|||
function co = cosn(H);
|
||||
|
||||
% function co = cosn(H);
|
||||
% computes the cosine of the angle between the H(:,1) and its
|
||||
% projection onto the span of H(:,2:end)
|
||||
|
||||
% Not the same as multiple correlation coefficient since the means are not
|
||||
% zero
|
||||
|
||||
y = H(:,1);
|
||||
X = H(:,2:end);
|
||||
|
||||
% y = H(:,1);
|
||||
% X = H(:,2:end);
|
||||
|
||||
yhat = X*(X\y);
|
||||
co = y'*yhat/sqrt((y'*y)*(yhat'*yhat));
|
||||
|
||||
|
||||
|
||||
|
|
@ -0,0 +1,77 @@
|
|||
function disp_identification(pdraws, idemodel, idemoments)
|
||||
|
||||
global bayestopt_
|
||||
|
||||
[SampleSize, npar] = size(pdraws);
|
||||
jok = 0;
|
||||
jokP = 0;
|
||||
jokJ = 0;
|
||||
jokPJ = 0;
|
||||
for j=1:npar,
|
||||
if any(idemodel.ind(j,:)==0),
|
||||
pno = 100*length(find(idemodel.ind(j,:)==0))/SampleSize;
|
||||
disp(['Parameter ',bayestopt_.name{j},' is not identified in the model for ',num2str(pno),'% of MC runs!' ])
|
||||
disp(' ')
|
||||
end
|
||||
if any(idemoments.ind(j,:)==0),
|
||||
pno = 100*length(find(idemoments.ind(j,:)==0))/SampleSize;
|
||||
disp(['Parameter ',bayestopt_.name{j},' is not identified by J moments for ',num2str(pno),'% of MC runs!' ])
|
||||
disp(' ')
|
||||
end
|
||||
if any(idemodel.ind(j,:)==1),
|
||||
iok = find(idemodel.ind(j,:)==1);
|
||||
jok = jok+1;
|
||||
kok(jok) = j;
|
||||
mmin(jok,1) = min(idemodel.Mco(j,iok));
|
||||
mmean(jok,1) = mean(idemodel.Mco(j,iok));
|
||||
mmax(jok,1) = max(idemodel.Mco(j,iok));
|
||||
[ipmax, jpmax] = find(abs(squeeze(idemodel.Pco(j,[1:j-1,j+1:end],iok)))>0.95);
|
||||
if ~isempty(ipmax)
|
||||
jokP = jokP+1;
|
||||
kokP(jokP) = j;
|
||||
ipmax(find(ipmax>=j))=ipmax(find(ipmax>=j))+1;
|
||||
[N,X]=hist(ipmax,[1:npar]);
|
||||
jpM(jokP)={find(N)};
|
||||
NPM(jokP)={N(find(N))./SampleSize.*100};
|
||||
pmeanM(jokP)={mean(squeeze(idemodel.Pco(j,find(N),iok))')};
|
||||
pminM(jokP)={min(squeeze(idemodel.Pco(j,find(N),iok))')};
|
||||
pmaxM(jokP)={max(squeeze(idemodel.Pco(j,find(N),iok))')};
|
||||
end
|
||||
end
|
||||
if any(idemoments.ind(j,:)==1),
|
||||
iok = find(idemoments.ind(j,:)==1);
|
||||
jokJ = jokJ+1;
|
||||
kokJ(jokJ) = j;
|
||||
mminJ(jokJ,1) = min(idemoments.Mco(j,iok));
|
||||
mmeanJ(jokJ,1) = mean(idemoments.Mco(j,iok));
|
||||
mmaxJ(jokJ,1) = max(idemoments.Mco(j,iok));
|
||||
[ipmax, jpmax] = find(abs(squeeze(idemoments.Pco(j,[1:j-1,j+1:end],iok)))>0.95);
|
||||
if ~isempty(ipmax)
|
||||
jokPJ = jokPJ+1;
|
||||
kokPJ(jokPJ) = j;
|
||||
ipmax(find(ipmax>=j))=ipmax(find(ipmax>=j))+1;
|
||||
[N,X]=hist(ipmax,[1:npar]);
|
||||
jpJ(jokPJ)={find(N)};
|
||||
NPJ(jokPJ)={N(find(N))./SampleSize.*100};
|
||||
pmeanJ(jokPJ)={mean(squeeze(idemoments.Pco(j,find(N),iok))')};
|
||||
pminJ(jokPJ)={min(squeeze(idemoments.Pco(j,find(N),iok))')};
|
||||
pmaxJ(jokPJ)={max(squeeze(idemoments.Pco(j,find(N),iok))')};
|
||||
end
|
||||
end
|
||||
end
|
||||
|
||||
dyntable('Multi collinearity in the model:',strvcat('param','min','mean','max'), ...
|
||||
strvcat(bayestopt_.name(kok)),[mmin, mmean, mmax],10,10,6);
|
||||
|
||||
dyntable('Multi collinearity for moments in J:',strvcat('param','min','mean','max'), ...
|
||||
strvcat(bayestopt_.name(kokJ)),[mminJ, mmeanJ, mmaxJ],10,10,6);
|
||||
|
||||
for j=1:length(kokP),
|
||||
dyntable([bayestopt_.name{kokP(j)},' pairwise correlations in the model'],strvcat(' ','min','mean','max'), ...
|
||||
strvcat(bayestopt_.name{jpM{j}}),[pminM{j}' pmeanM{j}' pmaxM{j}'],10,10,3);
|
||||
end
|
||||
|
||||
for j=1:length(kokPJ),
|
||||
dyntable([bayestopt_.name{kokPJ(j)},' pairwise correlations in J moments'],strvcat(' ','min','mean','max'), ...
|
||||
strvcat(bayestopt_.name{jpJ{j}}),[pminJ{j}' pmeanJ{j}' pmaxJ{j}'],10,10,3);
|
||||
end
|
||||
|
|
@ -0,0 +1,65 @@
|
|||
function [pdraws, idemodel, idemoments] = dynare_identification()
|
||||
|
||||
% main
|
||||
|
||||
global M_ options_ oo_ bayestopt_ estim_params_
|
||||
|
||||
|
||||
options_ = set_default_option(options_,'datafile',[]);
|
||||
options_.mode_compute = 0;
|
||||
[data,rawdata]=dynare_estimation_init([],1);
|
||||
% computes a first linear solution to set up various variables
|
||||
dynare_resolve;
|
||||
|
||||
options_.prior_mc=2000;
|
||||
|
||||
SampleSize = options_.prior_mc;
|
||||
|
||||
% results = prior_sampler(0,M_,bayestopt_,options_,oo_);
|
||||
|
||||
prior_draw(1,bayestopt_);
|
||||
IdentifDirectoryName = CheckPath('identification');
|
||||
|
||||
indx = estim_params_.param_vals(:,1);
|
||||
indexo=[];
|
||||
if ~isempty(estim_params_.var_exo)
|
||||
indexo = estim_params_.var_exo(:,1);
|
||||
end
|
||||
useautocorr = 0;
|
||||
nlags = 3;
|
||||
|
||||
iteration = 0;
|
||||
loop_indx = 0;
|
||||
|
||||
h = waitbar(0,'Monte Carlo identification checks ...');
|
||||
|
||||
while iteration < SampleSize,
|
||||
loop_indx = loop_indx+1;
|
||||
params = prior_draw();
|
||||
set_all_parameters(params);
|
||||
|
||||
[JJ, H] = getJJ(M_,oo_,options_,0,indx,indexo,bayestopt_.mf2,nlags,useautocorr);
|
||||
|
||||
if ~isempty(JJ),
|
||||
iteration = iteration + 1;
|
||||
pdraws(iteration,:) = params';
|
||||
[idemodel.Mco(:,iteration), idemoments.Mco(:,iteration), ...
|
||||
idemodel.Pco(:,:,iteration), idemoments.Pco(:,:,iteration), ...
|
||||
idemodel.cond(iteration), idemoments.cond(iteration), ...
|
||||
idemodel.ee(:,iteration), idemoments.ee(:,iteration), ...
|
||||
idemodel.ind(:,iteration), idemoments.ind(:,iteration), ...
|
||||
idemodel.indno{iteration}, idemoments.indno{iteration}] = ...
|
||||
identification_checks(H,JJ, bayestopt_);
|
||||
|
||||
waitbar(iteration/SampleSize,h)
|
||||
end
|
||||
end
|
||||
|
||||
close(h)
|
||||
|
||||
|
||||
save([IdentifDirectoryName '/' M_.fname '_identif'], 'pdraws', 'idemodel', 'idemoments')
|
||||
|
||||
|
||||
disp_identification(pdraws, idemodel, idemoments)
|
||||
|
||||
|
|
@ -0,0 +1,298 @@
|
|||
function [H, A0, B0, dA, dOm, info] = getH(M_,oo_,kronflag,indx,indexo)
|
||||
% computes derivative of reduced form linear model w.r.t. deep params
|
||||
|
||||
if nargin<3 | isempty(kronflag), kronflag = 0; end
|
||||
if nargin<4 | isempty(indx), indx = [1:M_.param_nbr];, end,
|
||||
if nargin<5 | isempty(indexo), indexo = [];, end,
|
||||
|
||||
|
||||
[I,J]=find(M_.lead_lag_incidence');
|
||||
yy0=oo_.steady_state(I);
|
||||
% yy0=[];
|
||||
% for j=1:size(M_.lead_lag_incidence,1);
|
||||
% yy0 = [ yy0; oo_.steady_state(find(M_.lead_lag_incidence(j,:)))];
|
||||
% end
|
||||
[df, gp] = feval([M_.fname,'_params_derivs'],yy0, oo_.exo_steady_state', M_.params, 1);
|
||||
[residual, g1 ] = feval([M_.fname,'_dynamic'],yy0, oo_.exo_steady_state', M_.params,1);
|
||||
|
||||
[residual, g1, g2 ] = feval([M_.fname,'_dynamic'],yy0, oo_.exo_steady_state', M_.params,1);
|
||||
[residual, gg1] = feval([M_.fname,'_static'],oo_.steady_state, oo_.exo_steady_state', M_.params);
|
||||
% df = feval([M_.fname,'_model_derivs'],yy0, oo_.exo_steady_state', M_.params, 1);
|
||||
dyssdtheta = -gg1\df;
|
||||
dyssdtheta = dyssdtheta(I,:);
|
||||
[nr, nc]=size(g2);
|
||||
nc = sqrt(nc);
|
||||
ns = max(max(M_.lead_lag_incidence));
|
||||
gp2 = gp*0;
|
||||
for j=1:nr,
|
||||
[I J]=ind2sub([nc nc],find(g2(j,:)));
|
||||
for i=1:nc,
|
||||
is = find(I==i);
|
||||
is = is(find(J(is)<=ns));
|
||||
if ~isempty(is),
|
||||
g20=full(g2(j,find(g2(j,:))));
|
||||
gp2(j,i,:)=g20(is)*dyssdtheta(J(is),:);
|
||||
end
|
||||
end
|
||||
end
|
||||
|
||||
|
||||
gp = gp+gp2;
|
||||
gp = gp(:,:,indx);
|
||||
param_nbr = length(indx);
|
||||
|
||||
% order_var = [oo_.dr.order_var; ...
|
||||
% [size(oo_dr.ghx,2)+1:size(oo_dr.ghx,2)+size(oo_.dr.transition_auxiliary_variables,1)]' ];
|
||||
% [A(order_var,order_var),B(order_var,:)]=dynare_resolve;
|
||||
[A,B,ys,info]=dynare_resolve;
|
||||
if info(1) > 0
|
||||
H = [];
|
||||
A0 = [];
|
||||
B0 = [];
|
||||
dA = [];
|
||||
dOm = [];
|
||||
return
|
||||
end
|
||||
|
||||
m = size(A,1);
|
||||
n = size(B,2);
|
||||
|
||||
|
||||
|
||||
klen = M_.maximum_endo_lag + M_.maximum_endo_lead + 1;
|
||||
k1 = M_.lead_lag_incidence(find([1:klen] ~= M_.maximum_endo_lag+1),:);
|
||||
a = g1(:,nonzeros(k1'));
|
||||
da = gp(:,nonzeros(k1'),:);
|
||||
kstate = oo_.dr.kstate;
|
||||
|
||||
GAM1 = zeros(M_.endo_nbr,M_.endo_nbr);
|
||||
Dg1 = zeros(M_.endo_nbr,M_.endo_nbr,param_nbr);
|
||||
% nf = find(M_.lead_lag_incidence(M_.maximum_endo_lag+2,:));
|
||||
% GAM1(:, nf) = -g1(:,M_.lead_lag_incidence(M_.maximum_endo_lag+2,nf));
|
||||
|
||||
k = find(kstate(:,2) == M_.maximum_endo_lag+2 & kstate(:,3));
|
||||
GAM1(:, kstate(k,1)) = -a(:,kstate(k,3));
|
||||
Dg1(:, kstate(k,1), :) = -da(:,kstate(k,3),:);
|
||||
k = find(kstate(:,2) > M_.maximum_endo_lag+2 & kstate(:,3));
|
||||
kk = find(kstate(:,2) > M_.maximum_endo_lag+2 & ~kstate(:,3));
|
||||
nauxe = 0;
|
||||
if ~isempty(k),
|
||||
ax(:,k)= -a(:,kstate(k,3));
|
||||
ax(:,kk)= 0;
|
||||
dax(:,k,:)= -da(:,kstate(k,3),:);
|
||||
dax(:,kk,:)= 0;
|
||||
nauxe=size(ax,2);
|
||||
GAM1 = [ax GAM1];
|
||||
Dg1 = cat(2,dax,Dg1);
|
||||
clear ax
|
||||
end
|
||||
|
||||
|
||||
[junk,cols_b,cols_j] = find(M_.lead_lag_incidence(M_.maximum_endo_lag+1, ...
|
||||
oo_.dr.order_var));
|
||||
GAM0 = zeros(M_.endo_nbr,M_.endo_nbr);
|
||||
Dg0 = zeros(M_.endo_nbr,M_.endo_nbr,param_nbr);
|
||||
GAM0(:,cols_b) = g1(:,cols_j);
|
||||
Dg0(:,cols_b,:) = gp(:,cols_j,:);
|
||||
% GAM0 = g1(:,M_.lead_lag_incidence(M_.maximum_endo_lag+1,:));
|
||||
|
||||
|
||||
k = find(kstate(:,2) == M_.maximum_endo_lag+1 & kstate(:,4));
|
||||
GAM2 = zeros(M_.endo_nbr,M_.endo_nbr);
|
||||
Dg2 = zeros(M_.endo_nbr,M_.endo_nbr,param_nbr);
|
||||
% nb = find(M_.lead_lag_incidence(1,:));
|
||||
% GAM2(:, nb) = -g1(:,M_.lead_lag_incidence(1,nb));
|
||||
GAM2(:, kstate(k,1)) = -a(:,kstate(k,4));
|
||||
Dg2(:, kstate(k,1), :) = -da(:,kstate(k,4),:);
|
||||
naux = 0;
|
||||
k = find(kstate(:,2) < M_.maximum_endo_lag+1 & kstate(:,4));
|
||||
kk = find(kstate(:,2) < M_.maximum_endo_lag+1 );
|
||||
if ~isempty(k),
|
||||
ax(:,k)= -a(:,kstate(k,4));
|
||||
ax = ax(:,kk);
|
||||
dax(:,k,:)= -da(:,kstate(k,4),:);
|
||||
dax = dax(:,kk,:);
|
||||
naux = size(ax,2);
|
||||
GAM2 = [GAM2 ax];
|
||||
Dg2 = cat(2,Dg2,dax);
|
||||
end
|
||||
|
||||
GAM0 = blkdiag(GAM0,eye(naux));
|
||||
Dg0 = cat(1,Dg0,zeros(naux+nauxe,M_.endo_nbr,param_nbr));
|
||||
Dg0 = cat(2,Dg0,zeros(m+nauxe,naux,param_nbr));
|
||||
Dg0 = cat(2,zeros(m+nauxe,nauxe,param_nbr),Dg0);
|
||||
|
||||
GAM2 = [GAM2 ; A(M_.endo_nbr+1:end,:)];
|
||||
Dg2 = cat(1,Dg2,zeros(naux+nauxe,m,param_nbr));
|
||||
Dg2 = cat(2,zeros(m+nauxe,nauxe,param_nbr),Dg2);
|
||||
GAM2 = [zeros(m+nauxe,nauxe) [GAM2; zeros(nauxe,m)]];
|
||||
GAM0 = [[zeros(m,nauxe);(eye(nauxe))] [GAM0; zeros(nauxe,m)]];
|
||||
|
||||
GAM3 = -g1(:,length(yy0)+1:end);
|
||||
% GAM3 = -g1(oo_.dr.order_var,length(yy0)+1:end);
|
||||
GAM3 = [GAM3; zeros(naux+nauxe,size(GAM3,2))];
|
||||
% Dg3 = -gp(oo_.dr.order_var,length(yy0)+1:end,:);
|
||||
Dg3 = -gp(:,length(yy0)+1:end,:);
|
||||
Dg3 = cat(1,Dg3,zeros(naux+nauxe,size(GAM3,2),param_nbr));
|
||||
|
||||
auxe = zeros(0,1);
|
||||
k0 = kstate(find(kstate(:,2) >= M_.maximum_endo_lag+2),:);;
|
||||
i0 = find(k0(:,2) == M_.maximum_endo_lag+2);
|
||||
for i=M_.maximum_endo_lag+3:M_.maximum_endo_lag+2+M_.maximum_endo_lead,
|
||||
i1 = find(k0(:,2) == i);
|
||||
n1 = size(i1,1);
|
||||
j = zeros(n1,1);
|
||||
for j1 = 1:n1
|
||||
j(j1) = find(k0(i0,1)==k0(i1(j1),1));
|
||||
end
|
||||
auxe = [auxe; i0(j)];
|
||||
i0 = i1;
|
||||
end
|
||||
auxe = [(1:size(auxe,1))' auxe(end:-1:1)];
|
||||
n_ir1 = size(auxe,1);
|
||||
nr = m + n_ir1;
|
||||
GAM1 = [[GAM1 zeros(size(GAM1,1),naux)]; zeros(naux+nauxe,m+nauxe)];
|
||||
GAM1(m+1:end,:) = sparse(auxe(:,1),auxe(:,2),ones(n_ir1,1),n_ir1,nr);
|
||||
Dg1 = cat(2,Dg1,zeros(M_.endo_nbr,naux,param_nbr));
|
||||
Dg1 = cat(1,Dg1,zeros(naux+nauxe,m+nauxe,param_nbr));
|
||||
|
||||
Ax = A;
|
||||
A1 = A;
|
||||
Bx = B;
|
||||
B1 = B;
|
||||
for j=1:M_.maximum_endo_lead-1,
|
||||
A1 = A1*A;
|
||||
B1 = A*B1;
|
||||
k = find(kstate(:,2) == M_.maximum_endo_lag+2+j );
|
||||
Ax = [A1(kstate(k,1),:); Ax];
|
||||
Bx = [B1(kstate(k,1),:); Bx];
|
||||
end
|
||||
Ax = [zeros(m+nauxe,nauxe) Ax];
|
||||
A0 = A;
|
||||
A=Ax; clear Ax A1;
|
||||
B0=B;
|
||||
B = Bx; clear Bx B1;
|
||||
|
||||
m = size(A,1);
|
||||
n = size(B,2);
|
||||
|
||||
% Dg1 = zeros(m,m,param_nbr);
|
||||
% Dg1(:, nf, :) = -gp(:,M_.lead_lag_incidence(3,nf),:);
|
||||
|
||||
% Dg0 = gp(:,M_.lead_lag_incidence(2,:),:);
|
||||
|
||||
% Dg2 = zeros(m,m,param_nbr);
|
||||
% Dg2(:, nb, :) = -gp(:,M_.lead_lag_incidence(1,nb),:);
|
||||
|
||||
% Dg3 = -gp(:,length(yy0)+1:end,:);
|
||||
|
||||
if kronflag==1, % kronecker products
|
||||
Dg0=reshape(Dg0,m^2,param_nbr);
|
||||
Dg1=reshape(Dg1,m^2,param_nbr);
|
||||
Dg2=reshape(Dg2,m^2,param_nbr);
|
||||
Dg3=reshape(Dg3,m*n,param_nbr);
|
||||
Om = B*B';
|
||||
Im = eye(m);
|
||||
Dm = duplication(m);
|
||||
DmPl = inv(Dm'*Dm)*Dm';
|
||||
Kmm = commutation(m,m);
|
||||
Kmn = commutation(m,n);
|
||||
|
||||
|
||||
Da = [eye(m^2),zeros(m^2,m*(m+1)/2)];
|
||||
Dom = [zeros(m*(m+1)/2,m^2),eye(m*(m+1)/2)];
|
||||
|
||||
|
||||
Df1Dtau = ( kron(Im,GAM0) - kron(A',GAM1) - kron(Im,GAM1*A) )*Da;
|
||||
|
||||
Df1Dthet = kron(A',Im)*Dg0 - kron( (A')^2,Im)*Dg1 - Dg2;
|
||||
|
||||
Df2Dtau = DmPl*( kron(GAM0,GAM0) - kron(GAM0,GAM1*A) - kron(GAM1*A,GAM0) + kron(GAM1*A,GAM1*A) )*Dm*Dom - ...
|
||||
DmPl*( kron(GAM0*Om,GAM1) + kron(GAM1,GAM0*Om)*Kmm - kron(GAM1*A*Om,GAM1) - kron(GAM1,GAM1*A*Om)*Kmm )*Da;
|
||||
|
||||
|
||||
Df2Dthet = DmPl*( kron(GAM0*Om,Im) + kron(Im,GAM0*Om)*Kmm - kron(Im,GAM1*A*Om)*Kmm - kron(GAM1*A*Om,Im) )*Dg0 - ...
|
||||
DmPl*( kron(GAM0*Om*A',Im) + kron(Im,GAM0*Om*A')*Kmm - kron(Im,GAM1*A*Om*A')*Kmm - kron(GAM1*A*Om*A',Im) )*Dg1 -...
|
||||
DmPl*( kron(GAM3,Im) + kron(Im,GAM3)*Kmn )*Dg3;
|
||||
|
||||
|
||||
DfDtau = [Df1Dtau;Df2Dtau];
|
||||
DfDthet = [Df1Dthet;Df2Dthet];
|
||||
|
||||
H = -DfDtau\DfDthet;
|
||||
x = reshape(H(1:m*m,:),m,m,param_nbr);
|
||||
y = reshape(Dm*H(m*m+1:end,:),m,m,param_nbr);
|
||||
x = x(nauxe+1:end,nauxe+1:end,:);
|
||||
y = y(nauxe+1:end,nauxe+1:end,:);
|
||||
m = size(y,1);
|
||||
x = reshape(x,m*m,param_nbr);
|
||||
Dm = duplication(m);
|
||||
DmPl = inv(Dm'*Dm)*Dm';
|
||||
y = DmPl*reshape(y,m*m,param_nbr);
|
||||
H = [x;y];
|
||||
|
||||
elseif kronflag==-1, % perturbation
|
||||
fun = 'thet2tau';
|
||||
params0 = M_.params;
|
||||
H = fdjac(fun,[sqrt(diag(M_.Sigma_e(indexo,indexo))); M_.params(indx)],indx,indexo);
|
||||
assignin('base','M_', M_);
|
||||
assignin('base','oo_', oo_);
|
||||
|
||||
else % generalized sylvester equation
|
||||
|
||||
% solves a*x+b*x*c=d
|
||||
a = (GAM0-GAM1*A);
|
||||
inva = inv(a);
|
||||
b = -GAM1;
|
||||
c = A;
|
||||
elem = zeros(m,m,param_nbr);
|
||||
d = elem;
|
||||
for j=1:param_nbr,
|
||||
elem(:,:,j) = (Dg0(:,:,j)-Dg1(:,:,j)*A);
|
||||
d(:,:,j) = Dg2(:,:,j)-elem(:,:,j)*A;
|
||||
end
|
||||
xx=sylvester3mr(a,b,c,d);
|
||||
if ~isempty(indexo),
|
||||
dSig = zeros(M_.exo_nbr,M_.exo_nbr);
|
||||
for j=1:length(indexo)
|
||||
dSig(indexo(j),indexo(j)) = 2*sqrt(M_.Sigma_e(indexo(j),indexo(j)));
|
||||
y = B*dSig*B';
|
||||
y = y(nauxe+1:end,nauxe+1:end);
|
||||
H(:,j) = [zeros((m-nauxe)^2,1); vech(y)];
|
||||
if nargout>3,
|
||||
dOm(:,:,j) = y;
|
||||
end
|
||||
dSig(indexo(j),indexo(j)) = 0;
|
||||
end
|
||||
end
|
||||
for j=1:param_nbr,
|
||||
x = xx(:,:,j);
|
||||
y = inva * (Dg3(:,:,j)-(elem(:,:,j)-GAM1*x)*B);
|
||||
y = y*M_.Sigma_e*B'+B*M_.Sigma_e*y';
|
||||
x = x(nauxe+1:end,nauxe+1:end);
|
||||
y = y(nauxe+1:end,nauxe+1:end);
|
||||
if nargout>3,
|
||||
dA(:,:,j+length(indexo)) = x;
|
||||
dOm(:,:,j+length(indexo)) = y;
|
||||
end
|
||||
H(:,j+length(indexo)) = [x(:); vech(y)];
|
||||
end
|
||||
% for j=1:param_nbr,
|
||||
% disp(['Derivatives w.r.t. ',M_.param_names(indx(j),:),', ',int2str(j),'/',int2str(param_nbr)])
|
||||
% elem = (Dg0(:,:,j)-Dg1(:,:,j)*A);
|
||||
% d = Dg2(:,:,j)-elem*A;
|
||||
% x=sylvester3(a,b,c,d);
|
||||
% % x=sylvester3a(x,a,b,c,d);
|
||||
% y = inva * (Dg3(:,:,j)-(elem-GAM1*x)*B);
|
||||
% y = y*B'+B*y';
|
||||
% x = x(nauxe+1:end,nauxe+1:end);
|
||||
% y = y(nauxe+1:end,nauxe+1:end);
|
||||
% H(:,j) = [x(:); vech(y)];
|
||||
% end
|
||||
|
||||
end
|
||||
|
||||
|
||||
return
|
||||
|
||||
|
||||
|
|
@ -0,0 +1,100 @@
|
|||
function [JJ, H, A, B, GAM] = getJJ(M_,oo_,options_,kronflag,indx,indexo,mf,nlags,useautocorr)
|
||||
|
||||
if nargin<5 | isempty(indx), indx = [1:M_.param_nbr];, end,
|
||||
if nargin<6 | isempty(indexo), indexo = [];, end,
|
||||
if nargin<8 | isempty(nlags), nlags=3; end,
|
||||
if nargin<9 | isempty(useautocorr), useautocorr=0; end,
|
||||
|
||||
if useautocorr,
|
||||
warning('off','MATLAB:divideByZero')
|
||||
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;
|
||||
assignin('base','M_', M_);
|
||||
assignin('base','oo_', oo_);
|
||||
else
|
||||
[H, A, B, dA, dOm, info] = getH(M_,oo_,kronflag,indx,indexo);
|
||||
if info(1) > 0
|
||||
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
|
||||
JJ(:,j+nexo) = dumm;
|
||||
end
|
||||
|
||||
if nargout >4,
|
||||
[GAM,stationary_vars] = th_autocovariances(oo_.dr,oo_.dr.order_var(mf),M_,options_);
|
||||
end
|
||||
end
|
||||
|
||||
if useautocorr,
|
||||
warning('on','MATLAB:divideByZero')
|
||||
end
|
||||
|
|
@ -0,0 +1,115 @@
|
|||
function [McoH, McoJ, PcoH, PcoJ, condH, condJ, eH, eJ, ind01, ind02, indnoH, indnoJ] = identification_checks(H,JJ, bayestopt_)
|
||||
|
||||
|
||||
% My suggestion is to have the following steps for identification check in
|
||||
% dynare:
|
||||
|
||||
% 1. check rank of H at theta
|
||||
npar = size(H,2);
|
||||
indno = {};
|
||||
ind1 = find(vnorm(H)~=0);
|
||||
H1 = H(:,ind1);
|
||||
[e1,e2] = eig(H1'*H1);
|
||||
eH = NaN(npar,1);
|
||||
eH(ind1) = e1(:,1);
|
||||
condH = cond(H1'*H1);
|
||||
ind2 = find(vnorm(JJ)~=0);
|
||||
JJ1 = JJ(:,ind2);
|
||||
[ee1,ee2] = eig(JJ1'*JJ1);
|
||||
eJ = NaN(npar,1);
|
||||
eJ(ind2) = ee1(:,1);
|
||||
condJ = cond(JJ1'*JJ1);
|
||||
|
||||
if rank(H)<npar
|
||||
ixno = 0;
|
||||
% - find out which parameters are involved,
|
||||
% using something like the vnorm and the eigenvalue decomposition of H;
|
||||
% disp('Some parameters are NOT identified in the model: H rank deficient')
|
||||
% disp(' ')
|
||||
if length(ind1)<npar,
|
||||
ixno = ixno + 1;
|
||||
indnoH(ixno) = {find(~ismember([1:npar],ind1))};
|
||||
% disp('Not identified params')
|
||||
% disp(bayestopt_.name(indnoH{1}))
|
||||
% disp(' ')
|
||||
end
|
||||
e0 = find(abs(diag(e2))<eps);
|
||||
for j=1:length(e0),
|
||||
ixno = ixno + 1;
|
||||
indnoH(ixno) = {ind1(find(e1(:,e0(j))))};
|
||||
% disp('Perfectly collinear parameters')
|
||||
% disp(bayestopt_.name(indnoH{ixno}))
|
||||
% disp(' ')
|
||||
end
|
||||
else % rank(H)==length(theta), go to 2
|
||||
% 2. check rank of J
|
||||
% disp('All parameters are identified at theta in the model (rank of H)')
|
||||
% disp(' ')
|
||||
end
|
||||
|
||||
if rank(JJ)<npar
|
||||
ixno = 0;
|
||||
% - find out which parameters are involved
|
||||
% disp('Some parameters are NOT identified by the moments included in J')
|
||||
% disp(' ')
|
||||
if length(ind2)<npar,
|
||||
ixno = ixno + 1;
|
||||
indnoJ(ixno) = {find(~ismember([1:npar],ind2))};
|
||||
end
|
||||
ee0 = find(abs(diag(ee2))<eps);
|
||||
for j=1:length(ee0),
|
||||
ixno = ixno + 1;
|
||||
indnoJ(ixno) = {ind2(find(ee1(:,ee0(j))))};
|
||||
% disp('Perfectly collinear parameters in moments J')
|
||||
% disp(bayestopt_.name(indnoJ{ixno}))
|
||||
% disp(' ')
|
||||
end
|
||||
else %rank(J)==length(theta) =>
|
||||
% disp('All parameters are identified at theta by the moments included in J')
|
||||
end
|
||||
|
||||
|
||||
% rank(H1)==size(H1,2)
|
||||
% rank(JJ1)==size(JJ1,2)
|
||||
|
||||
% to find near linear dependence problems I use
|
||||
|
||||
McoH = NaN(npar,1);
|
||||
McoJ = NaN(npar,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
|
||||
|
||||
% format long % some are nearly 1
|
||||
% McoJ
|
||||
|
||||
|
||||
% here there is no exact linear dependence, but there are several
|
||||
% near-dependencies, mostly due to strong pairwise colliniearities, which can
|
||||
% be checked using
|
||||
|
||||
PcoH = NaN(npar,npar);
|
||||
PcoJ = NaN(npar,npar);
|
||||
for ii = 1:size(H1,2);
|
||||
for jj = 1:size(H1,2);
|
||||
PcoH(ind1(ii),ind1(jj)) = [cosn([H1(:,ii),H1(:,jj)])];
|
||||
end
|
||||
end
|
||||
|
||||
for ii = 1:size(JJ1,2);
|
||||
for jj = 1:size(JJ1,2);
|
||||
PcoJ(ind2(ii),ind2(jj)) = [cosn([JJ1(:,ii),JJ1(:,jj)])];
|
||||
end
|
||||
end
|
||||
|
||||
ind01 = zeros(npar,1);
|
||||
ind02 = zeros(npar,1);
|
||||
ind01(ind1) = 1;
|
||||
ind02(ind2) = 1;
|
||||
|
||||
|
||||
|
||||
|
||||
|
|
@ -0,0 +1,70 @@
|
|||
function y = vnorm(A,varargin)
|
||||
% VNORM - Return the vector norm along specified dimension of A
|
||||
%
|
||||
% VNORM(A) returns the 2-norm along the first non-singleton
|
||||
% dimension of A
|
||||
% VNORM(A,dim) return the 2-norm along the dimension 'dim'
|
||||
% VNORM(A,dim,normtype) returns the norm specified by normtype
|
||||
% along the dimension 'dim'
|
||||
% VNORM(A,[],normtype) returns the norm specified by normtype along
|
||||
% the first non-singleton dimension of A
|
||||
%
|
||||
% normtype may be one of {inf,-inf,positive integer}.
|
||||
% For a given vector, v, these norms are defined as
|
||||
% inf: max(abs(v))
|
||||
% -inf: min(abs(v))
|
||||
% p (where p is a positive integer): sum(abs(v).^p)^(1/p)
|
||||
%
|
||||
% Examples:
|
||||
% A = [8 1 6; 3 5 7; 4 -9 2];
|
||||
%
|
||||
% %Columnwise 2-norm (Euclidean norm)
|
||||
% vnorm(A,1) = [9.4340 10.3441 9.4340];
|
||||
% vnorm(A,[],2) % Same as above (since first non-singleton dimensions
|
||||
% % is columnwise and default norm is 2-norm.
|
||||
% vnorm(A,[],[])% Again, same as above
|
||||
%
|
||||
% % Row-wise maximum of absolute values
|
||||
% vnorm(A,2,inf) = [8 7 9]';
|
||||
%
|
||||
% % Columnwise minimum of absolute values
|
||||
% vnorm(A,[],-inf) = [3 1 2];
|
||||
%
|
||||
% % Error: Use the inf type and not the string 'inf'
|
||||
% vnorm(A,[],'inf') % Wrong
|
||||
% vnorm(A,[],inf) % Correct
|
||||
|
||||
dim = [];
|
||||
ntype = [];
|
||||
|
||||
if nargin>1
|
||||
dim = varargin{1};
|
||||
if isempty(dim)
|
||||
idx = find(size(A)~=1);
|
||||
dim = idx(1);
|
||||
elseif dim~=floor(dim) || dim<1
|
||||
error('Dimension must be positive integer');
|
||||
end
|
||||
if nargin>2
|
||||
ntype = varargin{2};
|
||||
end
|
||||
end
|
||||
|
||||
|
||||
if isempty(ntype)
|
||||
y = sqrt(sum( abs(A).^2 , dim) );
|
||||
elseif ntype==1
|
||||
y = sum( abs(A) , dim );
|
||||
elseif isinf(ntype)
|
||||
if ntype > 0
|
||||
y=max(abs(A), [], dim);
|
||||
else
|
||||
y=min(abs(A), [], dim);
|
||||
end
|
||||
elseif ntype~=floor(ntype) || ntype<1
|
||||
error(['Norm type must be one of inf,-inf or a positive ' ...
|
||||
'integer']);
|
||||
else
|
||||
y = (sum( abs(A).^ntype , dim) ).^(1/ntype);
|
||||
end
|
||||
|
||||
Loading…
Reference in New Issue