diff --git a/doc/dynare.texi b/doc/dynare.texi
index 09c08cbd44e0a542b95899fb2d5deb802343249c..19d1c4e1de01c40e822b28f06f0fc4fa128a7822 100644
--- a/doc/dynare.texi
+++ b/doc/dynare.texi
@@ -5358,6 +5358,7 @@ decomposition of the above k-step ahead filtered values. Stores results in @code
@item diffuse_filter
+@anchor{diffuse_filter}
Uses the diffuse Kalman filter (as described in
@cite{Durbin and Koopman (2012)} and @cite{Koopman and Durbin
(2003)}) to estimate models with non-stationary observed variables.
@@ -7238,8 +7239,8 @@ for identification analysis. Default: @code{0}
@item ar = @var{INTEGER}
Maximum number of lags for moments in identification analysis. Default: @code{1}
-@item lik_init = @var{INTEGER}
-@xref{lik_init}.
+@item diffuse_filter = @var{INTEGER}
+@xref{diffuse_filter}.
@end table
diff --git a/matlab/dynare_identification.m b/matlab/dynare_identification.m
index 33f2c151e7901c05cb76793853f90ced902a7f08..65f24154f091ba2567c821c22e32d8c1cf9b6e0a 100644
--- a/matlab/dynare_identification.m
+++ b/matlab/dynare_identification.m
@@ -59,8 +59,24 @@ options_ident = set_default_option(options_ident,'periods',300);
options_ident = set_default_option(options_ident,'replic',100);
options_ident = set_default_option(options_ident,'advanced',0);
options_ident = set_default_option(options_ident,'normalize_jacobians',1);
+
+%Deal with non-stationary cases
+if isfield(options_ident,'diffuse_filter') %set lik_init and options_
+ options_ident.lik_init=3; %overwrites any other lik_init set
+ options_.diffuse_filter=options_ident.diffuse_filter; %make options_ inherit diffuse filter; will overwrite any conflicting lik_init in dynare_estimation_init
+else
+ if options_.diffuse_filter==1 %warning if estimation with diffuse filter was done, but not passed
+ warning('IDENTIFICATION:: Previously the diffuse_filter option was used, but it was not passed to the identification command. This may result in problems if your model contains unit roots.')
+ end
+ if isfield(options_ident,'lik_init')
+ if options_ident.lik_init==3 %user specified diffuse filter using the lik_init option
+ options_ident.analytic_derivation=0; %diffuse filter not compatible with analytic derivation
+ end
+ end
+end
options_ident = set_default_option(options_ident,'lik_init',1);
options_ident = set_default_option(options_ident,'analytic_derivation',1);
+
if isfield(options_ident,'nograph'),
options_.nograph=options_ident.nograph;
end
@@ -70,6 +86,9 @@ end
if isfield(options_ident,'graph_format'),
options_.graph_format=options_ident.graph_format;
end
+if isfield(options_ident,'prior_trunc'),
+ options_.prior_trunc=options_ident.prior_trunc;
+end
if options_ident.gsa_sample_file,
GSAFolder = checkpath('gsa',M_.dname);
@@ -122,8 +141,7 @@ options_.prior_mc = options_ident.prior_mc;
options_.options_ident = options_ident;
options_.Schur_vec_tol = 1.e-8;
options_.nomoments=0;
-options_.analytic_derivation=1;
-
+options_ = set_default_option(options_,'analytic_derivation',1);
options_ = set_default_option(options_,'datafile','');
options_.mode_compute = 0;
options_.plot_priors = 0;
diff --git a/matlab/gensylv/sylvester3.m b/matlab/gensylv/sylvester3.m
index aabdc6650c447c4a18422d03435b336d48000006..be4d8c574f39ad28bcd1e3b94cb4e9d210398ff4 100644
--- a/matlab/gensylv/sylvester3.m
+++ b/matlab/gensylv/sylvester3.m
@@ -52,6 +52,7 @@ else
end
i = 1;
c = zeros(n,1,p);
+c1 = zeros(n,1,p);
while i < m
if t(i+1,i) == 0
if i == 1
diff --git a/matlab/getH.m b/matlab/getH.m
index 17ed4ba449794aa2164de20375642a31703f3695..afdb8e1698684b097a971fd3dc8410fd9d8bea7f 100644
--- a/matlab/getH.m
+++ b/matlab/getH.m
@@ -19,10 +19,18 @@ function [H, dA, dOm, Hss, gp, d2A, d2Om, H2ss] = getH(A, B, M_,oo_,options_,kro
% 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<6 || isempty(kronflag), kronflag = 0; end
-if nargin<7 || isempty(indx), indx = [1:M_.param_nbr]; end,
-if nargin<8 || isempty(indexo), indexo = []; end,
-if nargin<9 || isempty(iv), iv = (1:length(A))'; end,
+if nargin<6 || isempty(kronflag)
+ kronflag = 0;
+end
+if nargin<7 || isempty(indx)
+ indx = [];
+end
+if nargin<8 || isempty(indexo)
+ indexo = [];
+end
+if nargin<9 || isempty(iv)
+ iv = (1:length(A))';
+end
[I,J]=find(M_.lead_lag_incidence');
yy0=oo_.dr.ys(I);
diff --git a/matlab/getJJ.m b/matlab/getJJ.m
index e27ed22b8d08355dcdc005725223766a6bbc9141..de13063ebf90fc102a7cf7442f600d9cf567528c 100644
--- a/matlab/getJJ.m
+++ b/matlab/getJJ.m
@@ -17,10 +17,18 @@ function [JJ, H, gam, gp, dA, dOm, dYss] = getJJ(A, B, M_,oo_,options_,kronflag,
% 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<7 || isempty(indx), indx = [1:M_.param_nbr];, end,
-if nargin<8 || isempty(indexo), indexo = [];, end,
-if nargin<10 || isempty(nlags), nlags=3; end,
-if nargin<11 || isempty(useautocorr), useautocorr=0; end,
+if nargin<7 || isempty(indx)
+% indx = [1:M_.param_nbr];
+end,
+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,
warning('off','MATLAB:divideByZero')
diff --git a/matlab/gsa/map_ident_.m b/matlab/gsa/map_ident_.m
index 867dfcfd63e133bf6319fe6e34f85c1c966627a0..61db1cafdc73ab2057b1cebfc59c1fe5e297f73c 100644
--- a/matlab/gsa/map_ident_.m
+++ b/matlab/gsa/map_ident_.m
@@ -64,7 +64,7 @@ if opt_gsa.load_ident_files==0,
if opt_gsa.morris==2,
pdraws = dynare_identification(options_.options_ident,[lpmatx lpmat(istable,:)]);
% [pdraws, TAU, GAM] = dynare_identification(options_.options_ident,[lpmatx lpmat(istable,:)]);
- if max(max(abs(pdraws-[lpmatx lpmat(istable,:)])))==0,
+ if ~isempty(pdraws) && max(max(abs(pdraws-[lpmatx lpmat(istable,:)])))==0,
disp(['Sample check OK ', num2str(max(max(abs(pdraws-[lpmatx lpmat(istable,:)]))))]),
clear pdraws;
end
diff --git a/matlab/gsa/mc_moments.m b/matlab/gsa/mc_moments.m
index 5c2cd034abc202e3d72218457c88391e069b2982..b2ea935385a7fd79363496d8b5f6f6b9c7a15f56 100644
--- a/matlab/gsa/mc_moments.m
+++ b/matlab/gsa/mc_moments.m
@@ -20,7 +20,7 @@ function [vdec, cc, ac] = mc_moments(mm, ss, dr)
global options_ M_ estim_params_ oo_
[nr1, nc1, nsam] = size(mm);
- nobs=size(options_.varobs,1);
+ nobs=size(options_.varobs,2);
disp('Computing theoretical moments ...')
h = dyn_waitbar(0,'Theoretical moments ...');
vdec = zeros(nobs,M_.exo_nbr,nsam);
@@ -37,7 +37,7 @@ global options_ M_ estim_params_ oo_
cc(:,:,j)=triu(corr);
dum=[];
for i=1:options_.ar
- dum=[dum, autocorr{i}];
+ dum=[dum, autocorr{i}];
end
ac(:,:,j)=dum;
dyn_waitbar(j/nsam,h)
diff --git a/matlab/gsa/prior_draw_gsa.m b/matlab/gsa/prior_draw_gsa.m
index 9c6b2cb16797e4237aa6e2e1da6198fcde9dd176..5dc30b79ff3e24ee54ae80a7738c5b0506c3e928 100644
--- a/matlab/gsa/prior_draw_gsa.m
+++ b/matlab/gsa/prior_draw_gsa.m
@@ -42,7 +42,7 @@ function pdraw = prior_draw_gsa(init,rdraw)
% You should have received a copy of the GNU General Public License
% along with Dynare. If not, see <http://www.gnu.org/licenses/>.
-global bayestopt_ options_
+global bayestopt_ options_ estim_params_ M_
persistent npar pshape p6 p7 p3 p4 lbcum ubcum
if init
@@ -55,7 +55,18 @@ if init
pdraw = zeros(npar,1);
lbcum = zeros(npar,1);
ubcum = ones(npar,1);
- bounds = prior_bounds(bayestopt_,options_);
+ [junk1,junk2,junk3,lb,ub,junk4] = set_prior(estim_params_,M_,options_); %Prepare bounds
+ if ~isempty(bayestopt_) && any(bayestopt_.pshape > 0)
+ % Set prior bounds
+ bounds = prior_bounds(bayestopt_,options_);
+ bounds.lb = max(bounds.lb,lb);
+ bounds.ub = min(bounds.ub,ub);
+ else % estimated parameters but no declared priors
+ % No priors are declared so Dynare will estimate the model by
+ % maximum likelihood with inequality constraints for the parameters.
+ bounds.lb = lb;
+ bounds.ub = ub;
+ end
% set bounds for cumulative probabilities
for i = 1:npar
switch pshape(i)
@@ -63,8 +74,8 @@ if init
p4(i) = min(p4(i),bounds.ub(i));
p3(i) = max(p3(i),bounds.lb(i));
case 3% Gaussian prior.
- lbcum(i) = 0.5 * erfc(-(bounds.lb(i)-p6(i))/p7(i) ./ sqrt(2));;
- ubcum(i) = 0.5 * erfc(-(bounds.ub(i)-p6(i))/p7(i) ./ sqrt(2));;
+ lbcum(i) = 0.5 * erfc(-(bounds.lb(i)-p6(i))/p7(i) ./ sqrt(2));
+ ubcum(i) = 0.5 * erfc(-(bounds.ub(i)-p6(i))/p7(i) ./ sqrt(2));
case 2% Gamma prior.
lbcum(i) = gamcdf(bounds.lb(i)-p3(i),p6(i),p7(i));
ubcum(i) = gamcdf(bounds.ub(i)-p3(i),p6(i),p7(i));
diff --git a/matlab/gsa/set_shocks_param.m b/matlab/gsa/set_shocks_param.m
index f0f7542527a3fb163ef4bfae9b08d1365ccf9eb5..8902dfb102cd06f01a5830788cd4a25d0ecad8e4 100644
--- a/matlab/gsa/set_shocks_param.m
+++ b/matlab/gsa/set_shocks_param.m
@@ -1,6 +1,8 @@
function set_shocks_param(xparam1)
+% function set_shocks_param(xparam1)
+% Set the structural and measurement error variances and covariances
-% Copyright (C) 2012 Dynare Team
+% Copyright (C) 2012-2015 Dynare Team
%
% This file is part of Dynare.
%
@@ -18,31 +20,86 @@ function set_shocks_param(xparam1)
% along with Dynare. If not, see <http://www.gnu.org/licenses/>.
global estim_params_ M_
-
- nvx = estim_params_.nvx;
- ncx = estim_params_.ncx;
- np = estim_params_.np;
- Sigma_e = M_.Sigma_e;
- offset = 0;
- if nvx
- offset = offset + nvx;
+
+nvx = estim_params_.nvx;
+ncx = estim_params_.ncx;
+nvn = estim_params_.nvn;
+ncn = estim_params_.ncn;
+
+Sigma_e = M_.Sigma_e;
+Correlation_matrix = M_.Correlation_matrix;
+
+H = M_.H;
+Correlation_matrix_ME = M_.Correlation_matrix_ME;
+% setting shocks variance on the diagonal of Covariance matrix; used later
+% for updating covariances
+if nvx
var_exo = estim_params_.var_exo;
for i=1:nvx
- k = var_exo(i,1);
- Sigma_e(k,k) = xparam1(i)^2;
+ k =var_exo(i,1);
+ Sigma_e(k,k) = xparam1(i)^2;
end
- end
-
- if ncx
- offset = offset + estim_params_.nvn;
+end
+% update offset
+offset = nvx;
+
+% setting measument error variance; on the diagonal of Covariance matrix; used later
+% for updating covariances
+if nvn
+ for i=1:nvn
+ k = estim_params_.nvn_observable_correspondence(i,1);
+ H(k,k) = xparam1(i+offset)^2;
+ end
+end
+
+% update offset
+offset = nvx+nvn;
+
+% setting shocks covariances
+if ncx
corrx = estim_params_.corrx;
for i=1:ncx
- k1 = corrx(i,1);
- k2 = corrx(i,2);
- Sigma_e(k1,k2) = xparam1(i+offset)*sqrt(Sigma_e_(k1,k1)*Sigma_e_(k2,k2));
- Sigma_e(k2,k1) = Sigma_e_(k1,k2);
+ k1 = corrx(i,1);
+ k2 = corrx(i,2);
+ Correlation_matrix(k1,k2) = xparam1(i+offset);
+ Correlation_matrix(k2,k1) = Correlation_matrix(k1,k2);
end
- end
-
+end
+%build covariance matrix from correlation matrix and variances already on
+%diagonal
+Sigma_e = diag(sqrt(diag(Sigma_e)))*Correlation_matrix*diag(sqrt(diag(Sigma_e)));
+%if calibrated covariances, set them now to their stored value
+if isfield(estim_params_,'calibrated_covariances')
+ Sigma_e(estim_params_.calibrated_covariances.position)=estim_params_.calibrated_covariances.cov_value;
+end
+% update offset
+offset = nvx+nvn+ncx;
+
+% setting measurement error covariances
+if ncn
+ corrn_observable_correspondence = estim_params_.corrn_observable_correspondence;
+ for i=1:ncn
+ k1 = corrn_observable_correspondence(i,1);
+ k2 = corrn_observable_correspondence(i,2);
+ Correlation_matrix_ME(k1,k2) = xparam1(i+offset);
+ Correlation_matrix_ME(k2,k1) = Correlation_matrix_ME(k1,k2);
+ end
+end
+%build covariance matrix from correlation matrix and variances already on
+%diagonal
+H = diag(sqrt(diag(H)))*Correlation_matrix_ME*diag(sqrt(diag(H)));
+%if calibrated covariances, set them now to their stored value
+if isfield(estim_params_,'calibrated_covariances_ME')
+ H(estim_params_.calibrated_covariances_ME.position)=estim_params_.calibrated_covariances_ME.cov_value;
+end
+
- M_.Sigma_e = Sigma_e;
\ No newline at end of file
+% updating matrices in M
+if nvx || ncx
+ M_.Sigma_e = Sigma_e;
+ M_.Correlation_matrix=Correlation_matrix;
+end
+if nvn || ncn
+ M_.H = H;
+ M_.Correlation_matrix_ME=Correlation_matrix_ME;
+end
\ No newline at end of file
diff --git a/matlab/gsa/stab_map_.m b/matlab/gsa/stab_map_.m
index 876a75e8a1b3e11000b1284366177755975e37fe..85ef3676127e1440a629ad1b2af64dd8efe06563 100644
--- a/matlab/gsa/stab_map_.m
+++ b/matlab/gsa/stab_map_.m
@@ -87,7 +87,7 @@ nshock = estim_params_.nvx;
nshock = nshock + estim_params_.nvn;
nshock = nshock + estim_params_.ncx;
nshock = nshock + estim_params_.ncn;
-lpmat0=[];
+lpmat0=zeros(Nsam,0);
xparam1=[];
pshape = bayestopt_.pshape(nshock+1:end);
@@ -96,7 +96,22 @@ p2 = bayestopt_.p2(nshock+1:end);
p3 = bayestopt_.p3(nshock+1:end);
p4 = bayestopt_.p4(nshock+1:end);
-bounds = prior_bounds(bayestopt_,options_);
+[junk1,junk2,junk3,lb,ub,junk4] = set_prior(estim_params_,M_,options_); %Prepare bounds
+if ~isempty(bayestopt_) && any(bayestopt_.pshape > 0)
+ % Set prior bounds
+ bounds = prior_bounds(bayestopt_,options_);
+ bounds.lb = max(bounds.lb,lb);
+ bounds.ub = min(bounds.ub,ub);
+else % estimated parameters but no declared priors
+ % No priors are declared so Dynare will estimate the model by
+ % maximum likelihood with inequality constraints for the parameters.
+ bounds.lb = lb;
+ bounds.ub = ub;
+ if opt_gsa.prior_range==0
+ warning('GSA:: When using ML, sampling from the prior is not possible. Setting prior_range=1')
+ opt_gsa.prior_range=1;
+ end
+end
if nargin==0,
OutputDirectoryName='';
@@ -143,7 +158,7 @@ if fload==0,
end
end
else %ilptau==0
- %[lpmat] = rand(Nsam,np);
+ [lpmat] = NaN(Nsam,np);
for j=1:np,
lpmat(:,j) = randperm(Nsam)'./(Nsam+1); %latin hypercube
end
@@ -151,7 +166,7 @@ if fload==0,
end
end
% try
- dummy=prior_draw_gsa(1);
+ dummy=prior_draw_gsa(1); %initialize persistent variables
% catch
% if pprior,
% if opt_gsa.prior_range==0;
@@ -242,6 +257,7 @@ if fload==0,
else
d = chol(inv(hh));
lp=randn(Nsam*2,nshock+np)*d+kron(ones(Nsam*2,1),xparam1');
+ lnprior=zeros(1,Nsam*2);
for j=1:Nsam*2,
lnprior(j) = any(lp(j,:)'<=bounds.lb | lp(j,:)'>=bounds.ub);
end
@@ -263,6 +279,7 @@ if fload==0,
iwrong=zeros(1,Nsam);
inorestriction=zeros(1,Nsam);
irestriction=zeros(1,Nsam);
+ infox=zeros(1,Nsam);
for j=1:Nsam,
M_ = set_all_parameters([lpmat0(j,:) lpmat(j,:)]',estim_params_,M_);
%try stoch_simul([]);
@@ -273,7 +290,7 @@ if fload==0,
[Tt,Rr,SteadyState,info,M_,options_,oo_] = dynare_resolve(M_,options_,oo_,'restrict');
end
infox(j,1)=info(1);
- if infox(j,1)==0 && ~exist('T'),
+ if infox(j,1)==0 && ~exist('T','var'),
dr_=oo_.dr;
if prepSA,
try
@@ -321,7 +338,7 @@ if fload==0,
% bayestopt_.restrict_columns, ...
% bayestopt_.restrict_aux);
end
- if ~exist('nspred'),
+ if ~exist('nspred','var'),
nspred = dr_.nspred; %size(dr_.ghx,2);
nboth = dr_.nboth;
nfwrd = dr_.nfwrd;
@@ -339,7 +356,7 @@ if fload==0,
istable(j)=0;
if isfield(dr_,'eigval')
egg(:,j) = sort(dr_.eigval);
- if exist('nspred')
+ if exist('nspred','var')
if any(isnan(egg(1:nspred,j)))
iwrong(j)=j;
else
@@ -349,7 +366,7 @@ if fload==0,
end
end
else
- if exist('egg'),
+ if exist('egg','var'),
egg(:,j)=ones(size(egg,1),1).*NaN;
end
iwrong(j)=j;
@@ -457,6 +474,7 @@ else
stoch_simul([]);
%T=zeros(size(dr_.ghx,1),size(dr_.ghx,2)+size(dr_.ghu,2),length(istable));
ntrans=length(istable);
+ yys=NaN(length(ys_),ntrans);
for j=1:ntrans,
M_.params(estim_params_.param_vals(:,1)) = lpmat(istable(j),:)';
%stoch_simul([]);
@@ -465,12 +483,12 @@ else
%[Tt,Rr,SteadyState,info] = dynare_resolve(bayestopt_.restrict_var_list,...
% bayestopt_.restrict_columns,...
% bayestopt_.restrict_aux);
- if ~exist('T')
+ if ~exist('T','var')
T=zeros(size(dr_.ghx,1),size(dr_.ghx,2)+size(dr_.ghu,2),ntrans);
end
dr_ = oo_.dr;
T(:,:,j) = [dr_.ghx dr_.ghu];
- if ~exist('nspred')
+ if ~exist('nspred','var')
nspred = dr_.nspred; %size(dr_.ghx,2);
nboth = dr_.nboth;
nfwrd = dr_.nfwrd;
diff --git a/matlab/gsa/th_moments.m b/matlab/gsa/th_moments.m
index 9dbe7f8f68a4d6f632566403a472dce8a41bb50d..519d53bed7e98225e9b37e8e706af19b5d8a1127 100644
--- a/matlab/gsa/th_moments.m
+++ b/matlab/gsa/th_moments.m
@@ -1,6 +1,7 @@
function [vdec, corr, autocorr, z, zz] = th_moments(dr,var_list)
+% [vdec, corr, autocorr, z, zz] = th_moments(dr,var_list)
-% Copyright (C) 2012 Dynare Team
+% Copyright (C) 2012-2015 Dynare Team
%
% This file is part of Dynare.
%
@@ -17,38 +18,38 @@ function [vdec, corr, autocorr, z, zz] = th_moments(dr,var_list)
% You should have received a copy of the GNU General Public License
% along with Dynare. If not, see <http://www.gnu.org/licenses/>.
- global M_ oo_ options_
-
- nvar = size(var_list,1);
- if nvar == 0
+global M_ oo_ options_
+
+nvar = size(var_list,2);
+if nvar == 0
nvar = length(dr.order_var);
ivar = [1:nvar]';
- else
+else
ivar=zeros(nvar,1);
for i=1:nvar
- i_tmp = strmatch(var_list(i,:),M_.endo_names,'exact');
- if isempty(i_tmp)
- error (['One of the variable specified does not exist']) ;
- else
- ivar(i) = i_tmp;
- end
+ i_tmp = strmatch(var_list(:,i),M_.endo_names,'exact');
+ if isempty(i_tmp)
+ error(['One of the variables specified does not exist']) ;
+ else
+ ivar(i) = i_tmp;
+ end
end
- end
-
- [gamma_y,stationary_vars] = th_autocovariances(dr,ivar,M_, options_);
- m = dr.ys(ivar(stationary_vars));
+end
+
+[gamma_y,stationary_vars] = th_autocovariances(dr,ivar,M_, options_);
+m = dr.ys(ivar(stationary_vars));
+
-
% i1 = find(abs(diag(gamma_y{1})) > 1e-12);
- i1 = [1:length(ivar)];
- s2 = diag(gamma_y{1});
- sd = sqrt(s2);
+i1 = [1:length(ivar)];
+s2 = diag(gamma_y{1});
+sd = sqrt(s2);
+
+
+z = [ m sd s2 ];
+mean = m;
+var = gamma_y{1};
-
- z = [ m sd s2 ];
- mean = m;
- var = gamma_y{1};
-
%'THEORETICAL MOMENTS';
%'MEAN','STD. DEV.','VARIANCE');
@@ -56,28 +57,27 @@ z;
%'VARIANCE DECOMPOSITION (in percent)';
if M_.exo_nbr>1,
-vdec = 100*gamma_y{options_.ar+2}(i1,:);
+ vdec = 100*gamma_y{options_.ar+2}(i1,:);
else
-vdec = 100*ones(size(gamma_y{1}(i1,1)));
-end
+ vdec = 100*ones(size(gamma_y{1}(i1,1)));
+end
%'MATRIX OF CORRELATIONS';
if options_.opt_gsa.useautocorr,
corr = gamma_y{1}(i1,i1)./(sd(i1)*sd(i1)');
corr = corr-diag(diag(corr))+diag(diag(gamma_y{1}(i1,i1)));
else
- corr = gamma_y{1}(i1,i1);
+ corr = gamma_y{1}(i1,i1);
end
- if options_.ar > 0
-%'COEFFICIENTS OF AUTOCORRELATION';
+if options_.ar > 0
+ %'COEFFICIENTS OF AUTOCORRELATION';
for i=1:options_.ar
- if options_.opt_gsa.useautocorr,
- autocorr{i} = gamma_y{i+1}(i1,i1);
- else
- autocorr{i} = gamma_y{i+1}(i1,i1).*(sd(i1)*sd(i1)');
- end
- zz(:,i) = diag(gamma_y{i+1}(i1,i1));
+ if options_.opt_gsa.useautocorr,
+ autocorr{i} = gamma_y{i+1}(i1,i1);
+ else
+ autocorr{i} = gamma_y{i+1}(i1,i1).*(sd(i1)*sd(i1)');
+ end
+ zz(:,i) = diag(gamma_y{i+1}(i1,i1));
end
- end
-
+end
+
-
\ No newline at end of file
diff --git a/matlab/identification_analysis.m b/matlab/identification_analysis.m
index eb209885c3e263ae79ef3f3d412f57e699ee2442..6e5e9521c01c5bdcf0d83d6da525ffef1ce2810a 100644
--- a/matlab/identification_analysis.m
+++ b/matlab/identification_analysis.m
@@ -86,7 +86,7 @@ if info(1)==0,
if init,
indJJ = (find(max(abs(JJ'),[],1)>1.e-8));
if isempty(indJJ) && any(any(isnan(JJ)))
- error('There are NaN in the JJ matrix. Please check whether your model has units roots and you forgot to set lik_init~=1.' )
+ error('There are NaN in the JJ matrix. Please check whether your model has units roots and you forgot to set diffuse_filter=1.' )
end
while length(indJJ)<nparam && nlags<10,
disp('The number of moments with non-zero derivative is smaller than the number of parameters')
@@ -118,14 +118,14 @@ if info(1)==0,
ide_strength_J_prior=NaN(1,nparam);
if init, %~isempty(indok),
normaliz = abs(params);
- if prior_exist,
+ if prior_exist,
if ~isempty(estim_params_.var_exo),
- normaliz1 = estim_params_.var_exo(:,7); % normalize with prior standard deviation
+ normaliz1 = estim_params_.var_exo(:,7)'; % normalize with prior standard deviation
else
normaliz1=[];
end
if ~isempty(estim_params_.param_vals),
- normaliz1 = [normaliz1; estim_params_.param_vals(:,7)]'; % normalize with prior standard deviation
+ normaliz1 = [normaliz1 estim_params_.param_vals(:,7)']; % normalize with prior standard deviation
end
% normaliz = max([normaliz; normaliz1]);
normaliz1(isinf(normaliz1)) = 1;
diff --git a/matlab/identification_checks.m b/matlab/identification_checks.m
index 78b12e1744b7564f8bc06bb9a6cda8eac19711f2..a4dc5553292f43966c4579907fea543ae3051e33 100644
--- a/matlab/identification_checks.m
+++ b/matlab/identification_checks.m
@@ -82,7 +82,7 @@ end
ixnoJ = 0;
-if rankJ<npar || rankJJ<npar || min(1-McoJ)<1.e-10
+if npar>0 && (rankJ<npar || rankJJ<npar || min(1-McoJ)<1.e-10)
% - find out which parameters are involved
% disp('Some parameters are NOT identified by the moments included in J')
% disp(' ')
diff --git a/matlab/plot_identification.m b/matlab/plot_identification.m
index 549224ef0dfae6b783857acdd8e702353d916111..16e4a83d116b79958c053f866fa36d51cb19cb0f 100644
--- a/matlab/plot_identification.m
+++ b/matlab/plot_identification.m
@@ -56,14 +56,22 @@ if SampleSize == 1,
subplot(211)
mmm = (idehess.ide_strength_J);
[ss, is] = sort(mmm);
- bar(log([idehess.ide_strength_J(:,is)' idehess.ide_strength_J_prior(:,is)']))
+ if ~all(isnan(idehess.ide_strength_J_prior))
+ bar(log([idehess.ide_strength_J(:,is)' idehess.ide_strength_J_prior(:,is)']))
+ else
+ bar(log([idehess.ide_strength_J(:,is)' ]))
+ end
set(gca,'xlim',[0 nparam+1])
set(gca,'xticklabel','')
dy = get(gca,'ylim');
for ip=1:nparam,
text(ip,dy(1),name{is(ip)},'rotation',90,'HorizontalAlignment','right','interpreter','none')
end
- legend('relative to param value','relative to prior std','Location','Best')
+ if ~all(isnan(idehess.ide_strength_J_prior))
+ legend('relative to param value','relative to prior std','Location','Best')
+ else
+ legend('relative to param value','Location','Best')
+ end
if idehess.flag_score,
title('Identification strength with asymptotic Information matrix (log-scale)')
else
@@ -71,14 +79,23 @@ if SampleSize == 1,
end
subplot(212)
- bar(log([idehess.deltaM(is) idehess.deltaM_prior(is)]))
+ if ~all(isnan(idehess.deltaM_prior))
+ bar(log([idehess.deltaM(is) idehess.deltaM_prior(is)]))
+ else
+ bar(log([idehess.deltaM(is)]))
+ end
+
set(gca,'xlim',[0 nparam+1])
set(gca,'xticklabel','')
dy = get(gca,'ylim');
for ip=1:nparam,
text(ip,dy(1),name{is(ip)},'rotation',90,'HorizontalAlignment','right','interpreter','none')
end
- legend('relative to param value','relative to prior std','Location','Best')
+ if ~all(isnan(idehess.deltaM_prior))
+ legend('relative to param value','relative to prior std','Location','Best')
+ else
+ legend('relative to param value','Location','Best')
+ end
if idehess.flag_score,
title('Sensitivity component with asymptotic Information matrix (log-scale)')
else
@@ -91,27 +108,30 @@ if SampleSize == 1,
skipline()
disp('Press ENTER to plot advanced diagnostics'), pause(5),
end
- hh = dyn_figure(options_,'Name',[tittxt, ' - Sensitivity plot']);
- subplot(211)
- mmm = (siJnorm)'./max(siJnorm);
- mmm1 = (siHnorm)'./max(siHnorm);
- mmm=[mmm mmm1];
- mmm1 = (siLREnorm)'./max(siLREnorm);
- offset=length(siHnorm)-length(siLREnorm);
- mmm1 = [NaN(offset,1); mmm1];
- mmm=[mmm mmm1];
-
- bar(log(mmm(is,:).*100))
- set(gca,'xlim',[0 nparam+1])
- set(gca,'xticklabel','')
- dy = get(gca,'ylim');
- for ip=1:nparam,
- text(ip,dy(1),name{is(ip)},'rotation',90,'HorizontalAlignment','right','interpreter','none')
+ if all(isnan([siJnorm';siHnorm';siLREnorm']))
+ fprintf('\nIDENTIFICATION: Skipping sensitivity plot, because standard deviation of parameters is NaN, possibly due to the use of ML.\n')
+ else
+ hh = dyn_figure(options_,'Name',[tittxt, ' - Sensitivity plot']);
+ subplot(211)
+ mmm = (siJnorm)'./max(siJnorm);
+ mmm1 = (siHnorm)'./max(siHnorm);
+ mmm=[mmm mmm1];
+ mmm1 = (siLREnorm)'./max(siLREnorm);
+ offset=length(siHnorm)-length(siLREnorm);
+ mmm1 = [NaN(offset,1); mmm1];
+ mmm=[mmm mmm1];
+
+ bar(log(mmm(is,:).*100))
+ set(gca,'xlim',[0 nparam+1])
+ set(gca,'xticklabel','')
+ dy = get(gca,'ylim');
+ for ip=1:nparam,
+ text(ip,dy(1),name{is(ip)},'rotation',90,'HorizontalAlignment','right','interpreter','none')
+ end
+ legend('Moments','Model','LRE model','Location','Best')
+ title('Sensitivity bars using derivatives (log-scale)')
+ dyn_saveas(hh,[IdentifDirectoryName '/' M_.fname '_sensitivity_' tittxt1 ],options_);
end
- legend('Moments','Model','LRE model','Location','Best')
- title('Sensitivity bars using derivatives (log-scale)')
- dyn_saveas(hh,[IdentifDirectoryName '/' M_.fname '_sensitivity_' tittxt1 ],options_);
-
% identificaton patterns
for j=1:size(idemoments.cosnJ,2),
pax=NaN(nparam,nparam);
diff --git a/matlab/set_all_parameters.m b/matlab/set_all_parameters.m
index 687259fc5d09abd4351304572d39f91b7d45bee9..5c1fbc8f0dc309eb13ea59589f2b789dd11dc50c 100644
--- a/matlab/set_all_parameters.m
+++ b/matlab/set_all_parameters.m
@@ -32,8 +32,12 @@ function M = set_all_parameters(xparam1,estim_params,M)
%! @sp 2
%! @end deftypefn
%@eod:
+%
+% Remarks: Changes to this file also need to be ported to
+% gsa/set_shocks_param.m
+
-% Copyright (C) 2003-2013 Dynare Team
+% Copyright (C) 2003-2015 Dynare Team
%
% This file is part of Dynare.
%
diff --git a/tests/Makefile.am b/tests/Makefile.am
index 706ff4074dd2f961e09f10f8938c9256180fe054..1eefe3098c30b9578abec5a08684c77fd7a66af1 100644
--- a/tests/Makefile.am
+++ b/tests/Makefile.am
@@ -6,6 +6,7 @@ MODFILES = \
estimation/t_proposal/fs2000_student.mod \
gsa/ls2003.mod \
gsa/ls2003a.mod \
+ gsa/cod_ML_morris/cod_ML_morris.mod \
ramst.mod \
ramst_a.mod \
ramst_static_tag.mod \
@@ -124,6 +125,9 @@ MODFILES = \
seeds.mod \
identification/kim/kim2.mod \
identification/as2007/as2007.mod \
+ identification/ident_unit_root/ident_unit_root.mod \
+ identification/rbc_ident/rbc_ident_std_as_structural_par.mod \
+ identification/rbc_ident/rbc_ident_varexo_only.mod \
simul/example1.mod \
simul/Solow_no_varexo.mod \
simul/simul_ZLB_purely_forward.mod \
@@ -304,6 +308,8 @@ deterministic_simulations/rbc_det_exo_lag_2c.o.trs: deterministic_simulations/rb
initval_file/ramst_initval_file.m.trs: initval_file/ramst_initval_file_data.m.tls
initval_file/ramst_initval_file.o.trs: initval_file/ramst_initval_file_data.o.tls
+identification/rbc_ident/rbc_ident_varexo_only.m.trs: identification/rbc_ident/rbc_ident_std_as_structural_par.m.trs
+identification/rbc_ident/rbc_ident_varexo_only.o.trs: identification/rbc_ident/rbc_ident_std_as_structural_par.o.trs
# Matlab TRS Files
diff --git a/tests/gsa/cod_ML_morris/cod_ML_morris.mod b/tests/gsa/cod_ML_morris/cod_ML_morris.mod
new file mode 100644
index 0000000000000000000000000000000000000000..6e38d896f81a7880649339cf82ba7e904dfa2534
--- /dev/null
+++ b/tests/gsa/cod_ML_morris/cod_ML_morris.mod
@@ -0,0 +1,135 @@
+var c, y_ro, psi, s, tt, pi_h, mc_ro, i_ro, pi_ro, pi_eu, pi_f, y_eu, mc_eu, i_eu, g, a, cp, bnr, gs, as ;
+varexo eps_a, eps_g, eps_cp, eps_s, eps_m, eps_as, eps_gs, eps_ms;
+parameters alfa, niu, delt, teta_h, teta_f, bet, fi, sigm, h, ro, psi_pi, tetas, rhos, psi_pis, psi_y, psi_deltay, rho_g, rho_a, rho_cp, rho_s, rho_gs, rho_as ;
+
+alfa = 0.425 ;
+bet = 0.99 ;
+
+niu = 0.6 ;
+delt = 0.3 ;
+teta_h = 0.7 ;
+teta_f = 0.7 ;
+fi = 2 ;
+sigm = 0.5 ;
+h = 0.8 ;
+ro = 0.5 ;
+psi_pi = 1.5 ;
+tetas = 0.7 ;
+rhos = 0.5 ;
+psi_pis = 1.5 ;
+psi_y = 0.25 ;
+psi_deltay = 0.25 ;
+rho_g = 0.5 ;
+rho_a = 0.5 ;
+rho_cp = 0.5 ;
+rho_s = 0.5 ;
+rho_gs = 0.5 ;
+rho_as = 0.5 ;
+
+model (linear) ;
+
+// (1) mk clearing condition
+(1-alfa) * c = y_ro - alfa * niu* (2-alfa)* tt - alfa * niu * psi - alfa * y_eu ;
+// (2) terms of trade
+tt - tt(-1) = pi_f - pi_h ;
+// (3) terms of trade, real exchange rate, law of one price gap
+s = psi + (1 - alfa) * tt;
+// (4) domestic producers price setting
+pi_h - delt * pi_h(-1) = (1 - teta_h) * (1- bet * teta_h) / teta_h * mc_ro + bet * ( pi_h(+1) - delt * pi_h) ;
+// (5) domestic real mg costs
+mc_ro = fi * y_ro - (1 + fi) * a + alfa * tt + sigm / (1-h) * (c - h* c(-1)) ;
+// (6) importers price setting
+pi_f - delt * pi_f(-1) = (1 - teta_f) * (1 - bet * teta_f) / teta_f * psi + bet * (pi_f (+1) - delt * pi_f) + cp;
+// (7) complete mk assumption
+c - h * c(-1) = y_eu - h * y_eu(-1) + (1 - h) / sigm * (psi + (1 - alfa) * tt) + g ;
+// (8) domestic inflation
+pi_ro = pi_h + alfa * (tt - tt(-1)) ;
+// (9) UIP
+i_ro - pi_ro(+1) - i_eu + pi_eu(+1) = s(+1) - s + bnr ;
+// (10) domestic monetary policy rule
+i_ro = ro * i_ro(-1) + (1 - ro) * (psi_pi * pi_ro + psi_y * y_ro + psi_deltay * (y_ro - y_ro(-1)) ) + eps_m ;
+// (11) foreign euler equation
+y_eu - h * y_eu(-1) = y_eu(+1) - h * y_eu - (1 - h) / sigm * (i_eu - pi_eu(+1)) + gs - gs(+1) ;
+// (12) foreign producers price setting
+pi_eu = (1 - tetas) * (1 - bet * tetas) / tetas * mc_eu + bet * pi_eu(+1) ;
+// (13) foreign real mg costs
+mc_eu = fi * y_eu - (1 + fi) * as + alfa * tt + sigm / (1-h) * (y_eu - h * y_eu(-1)) ;
+// (14) foreign monetary policy rule
+i_eu = rhos * i_eu(-1) + (1 - rhos) * (psi_pis * pi_eu + psi_y * y_eu) + eps_ms ;
+
+// exogenous processes
+// (15) domestic shock in preferences (demand sh)
+g = rho_g * g(-1) + eps_g ;
+// (16) domestic technological shock (supply sh)
+a = rho_a * a(-1) + eps_a ;
+// (17) exchange rate shock
+bnr = rho_s * bnr(-1) + eps_s ;
+// (18) foreign shock in preferences
+gs = rho_gs * gs(-1) + eps_gs ;
+// (19) foreign technological shock
+as = rho_as * as(-1) + eps_as ;
+// (20) cost-push shock
+cp = rho_cp * cp(-1) + eps_cp ;
+
+end ;
+
+shocks;
+var eps_g; stderr 1 ;
+var eps_a; stderr 1 ;
+var eps_cp; stderr 1 ;
+var eps_s; stderr 1 ;
+var eps_m ; stderr 1 ;
+var eps_gs; stderr 1 ;
+var eps_as; stderr 1 ;
+var eps_ms ; stderr 1 ;
+end ;
+
+/* steady ;
+*/
+
+check ;
+
+estimated_params ;
+niu, 0.6, 0.01, 0.999 ;
+delt, 0.3, 0.01, 0.999 ;
+teta_h, 0.7, 0.01, 0.999 ;
+teta_f, 0.7, 0.01, 0.999 ;
+fi, 2, 0.01, 5 ;
+sigm, 0.5, 0.01, 0.999 ;
+h, 0.8, 0.01, 0.999 ;
+ro, 0.5, 0.01, 0.999 ;
+psi_pi, 1.4, 0.01, 8 ;
+tetas, 0.6, 0.01, 0.999 ;
+rhos, 0.5, 0.01, 0.999 ;
+psi_pis, 1.4, 0.01, 8 ;
+psi_y, 0.25, 0.01, 0.999 ;
+psi_deltay, 0.25, 0.01, 0.999 ;
+rho_g, 0.3, 0.01, 0.999 ;
+rho_a, 0.3, 0.01, 0.999 ;
+rho_cp, 0.3, 0.01, 0.999 ;
+rho_s, 0.3, 0.01, 0.999 ;
+rho_gs, 0.3, 0.01, 0.999 ;
+rho_as, 0.3, 0.01, 0.999 ;
+
+stderr eps_g, 1, 0.01, 10 ;
+stderr eps_a, 1, 0.01, 10 ;
+stderr eps_cp, 1, 0.01, 10 ;
+stderr eps_s, 1, 0.01, 10 ;
+stderr eps_gs, 1, 0.01, 10 ;
+stderr eps_as, 1, 0.01, 10 ;
+% corr eps_g, eps_a, 0, -1, 1;
+end;
+
+varobs y_ro, pi_ro, i_ro, s, y_eu, pi_eu, i_eu, tt ;
+
+dynare_sensitivity (identification=1, nsam = 2000, lik_only = 1, morris=2) ;
+
+stoch_simul(order=2,irf=20) y_ro pi_ro i_ro s ;
+//order 1 - impulse response functions are simply the algebraic forward iteration of the model's policy
+//order 2 - impulse response functions will be the result of actual Monte Carlo simulations of future shocks
+//specify sucient periods in IRFs to see the graphs return to steady state
+
+check ;
+
+/* shock_decomposition y_ro ;
+*/
diff --git a/tests/identification/ident_unit_root/ident_unit_root.mod b/tests/identification/ident_unit_root/ident_unit_root.mod
new file mode 100644
index 0000000000000000000000000000000000000000..ef8b424d76a0280f1e850d17bdd07278bad8990e
--- /dev/null
+++ b/tests/identification/ident_unit_root/ident_unit_root.mod
@@ -0,0 +1,45 @@
+//Tests Identification command with ML and unit roots/diffuse filter option
+
+var y delta_y x z;
+
+varexo eps_x eps_z;
+
+parameters rho sigma_z sigma_x;
+
+// set parameter values
+sigma_z=0.001;
+sigma_x=0.01;
+rho=0.9;
+
+model;
+z=rho*z(-1)+sigma_z*eps_z;
+x=x(-1)+sigma_x*eps_x;
+y=x+z;
+delta_y=y-y(-1);
+end;
+
+steady_state_model;
+x=0;
+z=0;
+y=0;
+delta_y=0;
+end;
+
+//set shock variances
+shocks;
+ var eps_z=1;
+ var eps_x=1;
+end;
+
+steady;
+check;
+varobs y delta_y;
+stoch_simul(order=1,irf=0);
+
+
+estimated_params;
+rho, 0.9;
+sigma_z, 0.01;
+sigma_x, 0.01;
+end;
+identification(diffuse_filter,advanced=1,prior_trunc=0);
\ No newline at end of file
diff --git a/tests/identification/rbc_ident/rbc_ident_std_as_structural_par.mod b/tests/identification/rbc_ident/rbc_ident_std_as_structural_par.mod
new file mode 100644
index 0000000000000000000000000000000000000000..60ec3406ff621d27c9924a6fc0af27522833037b
--- /dev/null
+++ b/tests/identification/rbc_ident/rbc_ident_std_as_structural_par.mod
@@ -0,0 +1,95 @@
+% Real Business Cycle Model
+
+close all;
+
+%----------------------------------------------------------------
+% 1. Defining variables
+%----------------------------------------------------------------
+
+var y c k l x z g ysim xsim;
+
+varexo eps_z eps_g;
+parameters gn gz betta delta psi sigma theta rho_z eps_z_sigma rho_g eps_g_sigma q12 zbar gbar;
+
+
+%----------------------------------------------------------------
+% 2. Calibration
+%----------------------------------------------------------------
+
+gn = 0.00243918275778010;
+gz = 0.00499789993972673;
+betta = 0.9722^(1/4)/(1+gz);
+delta = 1-(1-0.0464)^(1/4);
+psi = 2.24;
+sigma = 1.000001;
+theta = 0.35;
+
+zbar = 0.0023;
+gbar = -0.0382;
+
+rho_z = 0.8;
+sig_z = 0.0086;
+rho_g = 0.8;
+sig_g = 0.0248;
+
+q12 = -0.0002;
+
+
+
+%----------------------------------------------------------------
+% 3. Model
+%----------------------------------------------------------------
+
+
+model;
+% Intratemporal Optimality Condition
+(psi*c)/(1-l) = (1-theta)*(k(-1)^theta)*(l^(-theta))*(exp(z)^(1-theta));
+% Intertemporal Optimality Condition
+(c^(-sigma))*((1-l)^(psi*(1-sigma))) = betta*((c(+1))^(-sigma))*((1-l(+1))^(psi*(1-sigma)))*(theta*(k^(theta-1))*((exp(z(+1))*l(+1))^(1-theta))+(1-delta));
+% Aggregate Resource Constraint
+y = c + exp(g) + x;
+% Capital Accumulation Law
+(1+gz)*(1+gn)*k = (1-delta)*k(-1) + x;
+% Production Function
+y = (k(-1)^theta)*((exp(z)*l)^(1-theta));
+
+% Stochastic Processes
+z = zbar + rho_z*z(-1) + eps_z_sigma*eps_z;
+g = gbar + rho_g*g(-1) + eps_g_sigma*eps_g;
+
+ysim = y;
+xsim = x;
+
+end;
+
+
+%----------------------------------------------------------------
+% 4. Computation
+%----------------------------------------------------------------
+
+
+shocks;
+var eps_z ; stderr 1;
+var eps_g ; stderr 1;
+end;
+
+
+initval;
+z = 0;
+k = (betta*theta)^(1/(1-theta))*exp(zbar/(1-rho_z));
+c = (1-betta*theta)/(betta*theta)*k;
+y = exp(zbar/(1-rho_z))^(1-theta)*k^theta;
+l = k/(((1-betta*(1-delta))/(betta*theta*(exp(zbar/(1-rho_z))^(1-theta))))^(1/(theta-1)));
+end;
+
+
+estimated_params;
+eps_z_sigma, 0.0001, 0,10;
+eps_g_sigma, 0.0001, 0,10;
+% betta, 0.9722^(1/4)/(1+gz),0,1;
+%corr eps_z, eps_g, 0.0001, -1,1;
+end;
+
+varobs ysim xsim;
+
+identification(advanced=1);
diff --git a/tests/identification/rbc_ident/rbc_ident_varexo_only.mod b/tests/identification/rbc_ident/rbc_ident_varexo_only.mod
new file mode 100644
index 0000000000000000000000000000000000000000..4d3aa8c16b11afb6c715c3f2b18d7ca607100ccb
--- /dev/null
+++ b/tests/identification/rbc_ident/rbc_ident_varexo_only.mod
@@ -0,0 +1,105 @@
+% Real Business Cycle Model
+
+close all;
+
+%----------------------------------------------------------------
+% 1. Defining variables
+%----------------------------------------------------------------
+
+var y c k l x z g ysim xsim;
+
+varexo eps_z eps_g;
+parameters gn gz betta delta psi sigma theta rho_z sig_z rho_g sig_g q12 zbar gbar;
+
+
+%----------------------------------------------------------------
+% 2. Calibration
+%----------------------------------------------------------------
+
+gn = 0.00243918275778010;
+gz = 0.00499789993972673;
+betta = 0.9722^(1/4)/(1+gz);
+delta = 1-(1-0.0464)^(1/4);
+psi = 2.24;
+sigma = 1.000001;
+theta = 0.35;
+
+zbar = 0.0023;
+gbar = -0.0382;
+
+rho_z = 0.8;
+sig_z = 0.0086;
+rho_g = 0.8;
+sig_g = 0.0248;
+
+q12 = -0.0002;
+
+
+
+%----------------------------------------------------------------
+% 3. Model
+%----------------------------------------------------------------
+
+
+model;
+% Intratemporal Optimality Condition
+(psi*c)/(1-l) = (1-theta)*(k(-1)^theta)*(l^(-theta))*(exp(z)^(1-theta));
+% Intertemporal Optimality Condition
+(c^(-sigma))*((1-l)^(psi*(1-sigma))) = betta*((c(+1))^(-sigma))*((1-l(+1))^(psi*(1-sigma)))*(theta*(k^(theta-1))*((exp(z(+1))*l(+1))^(1-theta))+(1-delta));
+% Aggregate Resource Constraint
+y = c + exp(g) + x;
+% Capital Accumulation Law
+(1+gz)*(1+gn)*k = (1-delta)*k(-1) + x;
+% Production Function
+y = (k(-1)^theta)*((exp(z)*l)^(1-theta));
+
+% Stochastic Processes
+z = zbar + rho_z*z(-1) + eps_z;
+g = gbar + rho_g*g(-1) + eps_g;
+
+ysim = y;
+xsim = x;
+
+end;
+
+
+%----------------------------------------------------------------
+% 4. Computation
+%----------------------------------------------------------------
+
+
+shocks;
+var eps_z ; stderr sig_z;
+var eps_g ; stderr sig_g;
+end;
+
+
+initval;
+z = 0;
+k = (betta*theta)^(1/(1-theta))*exp(zbar/(1-rho_z));
+c = (1-betta*theta)/(betta*theta)*k;
+y = exp(zbar/(1-rho_z))^(1-theta)*k^theta;
+l = k/(((1-betta*(1-delta))/(betta*theta*(exp(zbar/(1-rho_z))^(1-theta))))^(1/(theta-1)));
+end;
+
+
+estimated_params;
+stderr eps_z, 0.0001, 0,10;
+stderr eps_g, 0.0001, 0,10;
+% betta, 0.9,0,1;
+%corr eps_z, eps_g, 0.0001, -1,1;
+end;
+
+varobs ysim xsim;
+
+identification(advanced=1);
+temp=load([M_.dname filesep 'identification' filesep M_.fname '_ML_Starting_value_identif'])
+temp_comparison=load(['rbc_ident_std_as_structural_par' filesep 'identification' filesep 'rbc_ident_std_as_structural_par' '_ML_Starting_value_identif'])
+
+if max(abs(temp.idehess_point.ide_strength_J - temp_comparison.idehess_point.ide_strength_J))>1e-8 ||...
+ max(abs(temp.idehess_point.deltaM - temp_comparison.idehess_point.deltaM))>1e-8 ||...
+ max(abs(temp.idemoments_point.GAM- temp.idemoments_point.GAM))>1e-8 ||...
+ max(abs(temp.idemoments_point.GAM- temp.idemoments_point.GAM))>1e-8
+ error('Results do not match')
+
+end
\ No newline at end of file