diff --git a/matlab/cosn.m b/matlab/cosn.m
new file mode 100644
index 000000000..f7564c47e
--- /dev/null
+++ b/matlab/cosn.m
@@ -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));
+
+
+
diff --git a/matlab/disp_identification.m b/matlab/disp_identification.m
new file mode 100644
index 000000000..c00285d28
--- /dev/null
+++ b/matlab/disp_identification.m
@@ -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
diff --git a/matlab/dynare_identification.m b/matlab/dynare_identification.m
new file mode 100644
index 000000000..d5cd57e92
--- /dev/null
+++ b/matlab/dynare_identification.m
@@ -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)
+
diff --git a/matlab/getH.m b/matlab/getH.m
new file mode 100644
index 000000000..7b730be81
--- /dev/null
+++ b/matlab/getH.m
@@ -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
+
+
diff --git a/matlab/getJJ.m b/matlab/getJJ.m
new file mode 100644
index 000000000..ffa4c7dcf
--- /dev/null
+++ b/matlab/getJJ.m
@@ -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
\ No newline at end of file
diff --git a/matlab/identification_checks.m b/matlab/identification_checks.m
new file mode 100644
index 000000000..6eacf0ae1
--- /dev/null
+++ b/matlab/identification_checks.m
@@ -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;
+
+
+
+
diff --git a/matlab/vnorm.m b/matlab/vnorm.m
new file mode 100644
index 000000000..dbf5790f4
--- /dev/null
+++ b/matlab/vnorm.m
@@ -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
+