Merge commit '34fb9f1c62c63e7e4ab41656a7690167449e596f'

doc/dynare.texi

@@ -2158,7 +2158,7 @@ Specifies the correlation of two variables.
 In an estimation context, it is also possible to specify variances and
 covariances on endogenous variables: in that case, these values are
 interpreted as the calibration of the measurement errors on these
-variables.
+variables. This requires the @code{var_obs}-command to be specified before the @code{shocks}-block.
 
 Here is an example:
 
@@ -3747,14 +3747,14 @@ alternatives:
 
 @item stderr @var{VARIABLE_NAME}
 Indicates that the standard error of the exogenous variable
-@var{VARIABLE_NAME}, or of the observation error associated with
+@var{VARIABLE_NAME}, or of the observation error/measurement errors associated with
 endogenous observed variable @var{VARIABLE_NAME}, is to be estimated
 
 @item corr @var{VARIABLE_NAME1}, @var{VARIABLE_NAME2}
 Indicates that the correlation between the exogenous variables
 @var{VARIABLE_NAME1} and @var{VARIABLE_NAME2}, or the correlation of
-the observation errors associated with endogenous observed variables
+the observation errors/measurement errors associated with endogenous observed variables
-@var{VARIABLE_NAME1} and @var{VARIABLE_NAME2}, is to be estimated
+@var{VARIABLE_NAME1} and @var{VARIABLE_NAME2}, is to be estimated. Note that correlations set by previous @code{shocks}-blocks or @code{estimation}-commands are kept at their value set prior to estimation if they are not estimated again subsequently. Thus, the treatment is the same as in the case of deep parameters set during model calibration and not estimated.
 
 @item @var{PARAMETER_NAME}
 The name of a model parameter to be estimated

matlab/GetPosteriorParametersStatistics.m

@@ -15,7 +15,7 @@ function oo_ = GetPosteriorParametersStatistics(estim_params_, M_, options_, bay
 % SPECIAL REQUIREMENTS
 %   None.
 
-% Copyright (C) 2006-2012 Dynare Team
+% Copyright (C) 2006-2013 Dynare Team
 %
 % This file is part of Dynare.
 %
@@ -163,7 +163,7 @@ end
 if nvn
     type = 'measurement_errors_std';
     if TeX
-        fid = TeXBegin(OutputDirectoryName,M_.fname,3,'standard deviation of measurement errors')
+        fid = TeXBegin(OutputDirectoryName,M_.fname,3,'standard deviation of measurement errors');
     end
     disp(' ')
     disp('standard deviation of measurement errors')
@@ -174,17 +174,17 @@ if nvn
             Draws = GetAllPosteriorDraws(ip,FirstMhFile,FirstLine,TotalNumberOfMhFiles,NumberOfDraws);
             [post_mean, post_median, post_var, hpd_interval, post_deciles, density] = ...
                 posterior_moments(Draws,1,options_.mh_conf_sig);
-            name = deblank(options_.varobs(estim_params_.var_endo(i,1),:));
+            name = deblank(options_.varobs(estim_params_.nvn_observable_correspondence(i,1),:));
             oo_ = Filloo(oo_,name,type,post_mean,hpd_interval,post_median,post_var,post_deciles,density);
         else
             try
-                name = deblank(options_.varobs(estim_params_.var_endo(i,1),:));
+                name = deblank(options_.varobs(estim_params_.nvn_observable_correspondence(i,1),:));
                 [post_mean,hpd_interval,post_var] = Extractoo(oo_,name,type);
             catch
                 Draws = GetAllPosteriorDraws(ip,FirstMhFile,FirstLine,TotalNumberOfMhFiles,NumberOfDraws);
                 [post_mean, post_median, post_var, hpd_interval, post_deciles, density] = ...
                     posterior_moments(Draws,1,options_.mh_conf_sig);
-                name = deblank(options_.varobs(estim_params_.var_endo(i,1),:));
+                name = deblank(options_.varobs(estim_params_.nvn_observable_correspondence(i,1),:));
                 oo_ = Filloo(oo_,name,type,post_mean,hpd_interval,post_median,post_var,post_deciles,density);
             end
         end

matlab/PlotPosteriorDistributions.m

@@ -16,7 +16,7 @@ function oo_ = PlotPosteriorDistributions(estim_params_, M_, options_, bayestopt
 % SPECIAL REQUIREMENTS
 %    none
 
-% Copyright (C) 2005-2012 Dynare Team
+% Copyright (C) 2005-2013 Dynare Team
 %
 % This file is part of Dynare.
 %
@@ -89,7 +89,7 @@ for i=1:npar
             eval(['pmod = oo_.posterior_mode.shocks_std.' name ';'])
         end
     elseif i <= nvx+nvn
-        name = deblank(options_.varobs(estim_params_.var_endo(i-nvx,1),:));
+        name = deblank(options_.varobs(estim_params_.nvn_observable_correspondence(i-nvx,1),:));
         eval(['x1 = oo_.posterior_density.measurement_errors_std.' name '(:,1);'])
         eval(['f1 = oo_.posterior_density.measurement_errors_std.' name '(:,2);'])    
         eval(['oo_.prior_density.mearsurement_errors_std.' name '(:,1) = x2;'])

matlab/dyn_forecast.m

@@ -87,8 +87,10 @@ switch task
                 i_var_obs = [ i_var_obs; tmp];
                 trend_coeffs = [trend_coeffs; oo_.Smoother.TrendCoeffs(i)];
             end
-        end         
+        end
-        trend = trend_coeffs*(gend+(1-M_.maximum_lag:horizon));
+        if ~isempty(trend_coeffs) 
+          trend = trend_coeffs*(gend+(1-M_.maximum_lag:horizon)); 
+        end
     end
     global bayestopt_
     if isfield(bayestopt_,'mean_varobs')

matlab/dynare_estimation_1.m

@@ -77,6 +77,24 @@ if ~isequal(estim_params_.ncx,nnz(tril(M_.Sigma_e,-1)))
     end
 end
 
+M_.H_is_diagonal = 1;
+if estim_params_.ncn || ~isequal(nnz(M_.H),length(M_.H))
+    M_.H_is_diagonal = 0;
+end
+
+% Set the correlation matrix of measurement errors if necessary.
+if ~isequal(estim_params_.ncn,nnz(tril(M_.H,-1)))
+    M_.Correlation_matrix_ME = diag(1./sqrt(diag(M_.H)))*M_.H*diag(1./sqrt(diag(M_.H)));
+    % Remove NaNs appearing because of variances calibrated to zero.
+    if any(isnan(M_.Correlation_matrix_ME))
+        zero_variance_idx = find(~diag(M_.H));
+        for i=1:length(zero_variance_idx)
+            M_.Correlation_matrix_ME(zero_variance_idx(i),:) = 0;
+            M_.Correlation_matrix_ME(:,zero_variance_idx(i)) = 0;
+        end
+    end
+end
+
 data = dataset_.data;
 rawdata = dataset_.rawdata;
 data_index = dataset_.missing.aindex;
@@ -554,7 +572,7 @@ if any(bayestopt_.pshape > 0) && ~options_.mh_posterior_mode_estimation
         disp(tit1)
         ip = nvx+1;
         for i=1:nvn
-            name = deblank(options_.varobs(estim_params_.var_endo(i,1),:));
+            name = deblank(options_.varobs(estim_params_.nvn_observable_correspondence(i,1),:));
             disp(sprintf('%-*s %7.3f %8.4f %7.4f %7.4f %4s %6.4f', ...
                          header_width,name,bayestopt_.p1(ip), ...
                          xparam1(ip),stdh(ip),tstath(ip), ...
@@ -653,7 +671,7 @@ elseif ~any(bayestopt_.pshape > 0) && ~options_.mh_posterior_mode_estimation
         disp(tit1)
         ip = nvx+1;
         for i=1:nvn
-            name = deblank(options_.varobs(estim_params_.var_endo(i,1),:));
+            name = deblank(options_.varobs(estim_params_.nvn_observable_correspondence(i,1),:));
             disp(sprintf('%-*s %8.4f %7.4f %7.4f',header_width,name,xparam1(ip),stdh(ip),tstath(ip)))
             eval(['oo_.mle_mode.measurement_errors_std.' name ' = xparam1(ip);']);
             eval(['oo_.mle_std.measurement_errors_std.' name ' = stdh(ip);']);
@@ -799,7 +817,7 @@ if any(bayestopt_.pshape > 0) && options_.TeX %% Bayesian estimation (posterior
         fprintf(fidTeX,'\\hline \\hline \\endlastfoot \n');
         ip = nvx+1;
         for i=1:nvn
-            idx = strmatch(options_.varobs(estim_params_.var_endo(i,1),:),M_.endo_names);
+            idx = strmatch(options_.varobs(estim_params_.nvn_observable_correspondence(i,1),:),M_.endo_names);
             fprintf(fidTeX,'$%s$ & %4s & %7.3f & %6.4f & %8.4f & %7.4f \\\\ \n',...
                     deblank(M_.endo_names_tex(idx,:)), ...
                     deblank(pnames(bayestopt_.pshape(ip)+1,:)), ...
@@ -1100,7 +1118,8 @@ if (~((any(bayestopt_.pshape > 0) && options_.mh_replic) || (any(bayestopt_.psha
                 hh = dyn_figure(options_,'Name','Smoothed observation errors');
                 NAMES = [];
                 if options_.TeX, TeXNAMES = []; end
-                for i=1:min(nstar,number_of_plots_to_draw-(nbplt-1)*nstar)
+                nstar0=min(nstar,number_of_plots_to_draw-(nbplt-1)*nstar);
+                for i=1:nstar0
                     k = (plt-1)*nstar+i;
                     subplot(nr,nc,i);
                     plot([1 gend],[0 0],'-r','linewidth',.5)
@@ -1131,7 +1150,7 @@ if (~((any(bayestopt_.pshape > 0) && options_.mh_replic) || (any(bayestopt_.psha
                 dyn_saveas(hh,[M_.fname '_SmoothedObservationErrors' int2str(plt)],options_);
                 if options_.TeX
                     fprintf(fidTeX,'\\begin{figure}[H]\n');
-                    for jj = 1:nstar
+                    for jj = 1:nstar0
                         fprintf(fidTeX,'\\psfrag{%s}[1][][0.5][0]{%s}\n',deblank(NAMES(jj,:)),deblank(TeXNAMES(jj,:)));
                     end
                     fprintf(fidTeX,'\\centering \n');
@@ -1168,7 +1187,7 @@ if (~((any(bayestopt_.pshape > 0) && options_.mh_replic) || (any(bayestopt_.psha
         for i=1:nstar0,
             k = (plt-1)*nstar+i;
             subplot(nr,nc,i);
-            plot(1:gend,yf(k,:),'--r','linewidth',1)
+            plot(1:gend,yf(k,:),'-r','linewidth',1)
             hold on
             plot(1:gend,rawdata(:,k),'--k','linewidth',1)
             hold off

matlab/get_posterior_parameters.m

@@ -12,7 +12,7 @@ function xparam = get_posterior_parameters(type)
 % SPECIAL REQUIREMENTS
 %   None.
 
-% Copyright (C) 2006-2009 Dynare Team
+% Copyright (C) 2006-2013 Dynare Team
 %
 % This file is part of Dynare.
 %
@@ -49,7 +49,7 @@ for i=1:nvx
 end
 
 for i=1:nvn
-    k1 = estim_params_.var_endo(i,1);
+    k1 = estim_params_.nvn_observable_correspondence(i,1);
     name1 = deblank(options_.varobs(k1,:));
     xparam(m) = eval(['oo_.posterior_' type '.measurement_errors_std.' name1]);
     m = m+1;
@@ -67,8 +67,8 @@ for i=1:ncx
 end
 
 for i=1:ncn
-    k1 = estim_params_.corrn(i,1);
+    k1 = estim_params_.corrn_observable_correspondence(i,1);
-    k2 = estim_params_.corrn(i,2);
+    k2 = estim_params_.corrn_observable_correspondence(i,2);
     name1 = deblank(options_.varobs(k1,:));
     name2 = deblank(options_.varobs(k2,:));
     xparam(m) = eval(['oo_.posterior_' type '.measurement_errors_corr.' name1 '_' name2]);

matlab/get_the_name.m

@@ -40,7 +40,7 @@ function [nam,texnam] = get_the_name(k,TeX,M_,estim_params_,options_)
 %! @end deftypefn
 %@eod:
 
-% Copyright (C) 2004-2011 Dynare Team
+% Copyright (C) 2004-2013 Dynare Team
 %
 % This file is part of Dynare.
 %
@@ -73,11 +73,11 @@ if k <= nvx
         texnam = ['$ SE_{' tname '} $'];
     end
 elseif  k <= (nvx+nvn)
-    vname = deblank(options_.varobs(estim_params_.var_endo(k-estim_params_.nvx,1),:));
+    vname = deblank(options_.varobs(estim_params_.nvn_observable_correspondence(k-estim_params_.nvx,1),:));
     nam=['SE_EOBS_',vname];
     if TeX
         tname  = deblank(M_.endo_names_tex(estim_params_.var_endo(k-estim_params_.nvx,1),:));
-        texnam = ['$ SE_{' tname '} $'];
+        texnam = ['$ EOBS SE_{' tname '} $'];
     end
 elseif  k <= (nvx+nvn+ncx)
     jj = k - (nvx+nvn);
@@ -97,7 +97,7 @@ elseif  k <= (nvx+nvn+ncx+ncn)
     nam=['CC_EOBS_' vname];
     if TeX
         tname  = [deblank(M_.endo_names_tex(k1,:)) ',' deblank(M_.endo_names_tex(k2,:))];
-        texnam =['$ CC_{' tname '} $'];
+        texnam =['$ EOBS CC_{' tname '} $'];
     end
 else
     jj = k - (nvx+nvn+ncx+ncn);

matlab/global_initialization.m

@@ -422,6 +422,7 @@ oo_.exo_det_simul = [];
 M_.params = [];
 M_.endo_histval = [];
 M_.Correlation_matrix = [];
+M_.Correlation_matrix_ME = [];
 
 % homotopy
 options_.homotopy_mode = 0;

matlab/prior_posterior_statistics.m

@@ -19,7 +19,7 @@ function prior_posterior_statistics(type,dataset)
 % See the comments random_walk_metropolis_hastings.m funtion.
 
 
-% Copyright (C) 2005-2012 Dynare Team
+% Copyright (C) 2005-2013 Dynare Team
 %
 % This file is part of Dynare.
 %
@@ -298,20 +298,20 @@ end
 
 if options_.filtered_vars
     pm3(endo_nbr,gend,ifil(1),B,'Updated Variables',...
-        '',varlist,'tit_tex',M_.endo_names,...
+        '',varlist,M_.endo_names_tex,M_.endo_names,...
         varlist,'UpdatedVariables',DirectoryName, ...
         '_update');
     pm3(endo_nbr,gend+1,ifil(4),B,'One step ahead forecast',...
-        '',varlist,'tit_tex',M_.endo_names,...
+        '',varlist,M_.endo_names_tex,M_.endo_names,...
         varlist,'FilteredVariables',DirectoryName,'_filter_step_ahead');
 end
 
 if options_.forecast
     pm3(endo_nbr,horizon+maxlag,ifil(6),B,'Forecasted variables (mean)',...
-        '',varlist,'tit_tex',M_.endo_names,...
+        '',varlist,M_.endo_names_tex,M_.endo_names,...
         varlist,'MeanForecast',DirectoryName,'_forc_mean');
     pm3(endo_nbr,horizon+maxlag,ifil(6),B,'Forecasted variables (point)',...
-        '',varlist,'tit_tex',M_.endo_names,...
+        '',varlist,M_.endo_names_tex,M_.endo_names,...
         varlist,'PointForecast',DirectoryName,'_forc_point');
 end
 

matlab/prior_posterior_statistics_core.m

@@ -140,7 +140,7 @@ if run_smoother
   end
 end
 if nvn
-  stock_error = NaN(endo_nbr,gend,MAX_nerro);
+  stock_error = NaN(size(varobs,1),gend,MAX_nerro);
 end
 if naK
     stock_filter_step_ahead =NaN(length(options_.filter_step_ahead),endo_nbr,gend+max(options_.filter_step_ahead),MAX_naK);

matlab/row_header_width.m

@@ -13,7 +13,7 @@ function w=row_header_width(M_,estim_params_,bayestopt_)
 % SPECIAL REQUIREMENTS
 %   None.
 
-% Copyright (C) 2006-2009 Dynare Team
+% Copyright (C) 2006-2013 Dynare Team
 %
 % This file is part of Dynare.
 %
@@ -64,7 +64,7 @@ if ncx
     end
 end
 if ncn
-    for i=1:nvn
+    for i=1:ncn
         k1 = estim_params_.corrn(i,1);
         k2 = estim_params_.corrn(i,2);
         w = max(w,length(deblank(M_.endo_names(k1,:)))...

matlab/set_all_parameters.m

@@ -71,9 +71,8 @@ offset = nvx;
 
 % setting measument error variance
 if nvn
-    var_endo = estim_params.var_endo;
     for i=1:nvn
-        k = var_endo(i,1);
+        k = estim_params.nvn_observable_correspondence(i,1);
         H(k,k) = xparam1(i+offset)^2;
     end
 end
@@ -83,7 +82,7 @@ offset = nvx+nvn;
 
 % setting shocks covariances
 if ~isempty(M.Correlation_matrix)
-    Sigma_e = diag(sqrt(diag(Sigma_e)))*M.Correlation_matrix*diag(sqrt(diag(Sigma_e)));
+    Sigma_e = diag(sqrt(diag(Sigma_e)))*M.Correlation_matrix*diag(sqrt(diag(Sigma_e))); % use of old correlation matrix is correct due to the diagonal structure and later only using the hence correctly updated diagonal entries of Sigma_e
 end
 if ncx
     corrx = estim_params.corrx;
@@ -100,11 +99,16 @@ end
 % update offset
 offset = nvx+nvn+ncx;
 % setting measurement error covariances
+if ~isempty(M.Correlation_matrix_ME)
+    H = diag(sqrt(diag(H)))*M.Correlation_matrix_ME*diag(sqrt(diag(H)));
+end
 if ncn
-    corrn = estim_params.corrn;
+    corrn_observable_correspondence = estim_params.corrn_observable_correspondence;
     for i=1:ncn
-        k1 = corr(i,1);
+        k1 = corrn_observable_correspondence(i,1);
-        k2 = corr(i,2);
+        k2 = corrn_observable_correspondence(i,2);
+        M.Correlation_matrix_ME(k1,k2) = xparam1(i+offset);
+        M.Correlation_matrix_ME(k2,k1) = M.Correlation_matrix_ME(k1,k2);
         H(k1,k2) = xparam1(i+offset)*sqrt(H(k1,k1)*H(k2,k2));
         H(k2,k1) = H(k1,k2);
     end
@@ -122,6 +126,6 @@ end
 if nvx || ncx
     M.Sigma_e = Sigma_e;
 end
-if nvn
+if nvn || ncn
     M.H = H;
 end

matlab/set_prior.m

@@ -18,7 +18,7 @@ function [xparam1, estim_params_, bayestopt_, lb, ub, M_]=set_prior(estim_params
 % SPECIAL REQUIREMENTS
 %    None
 
-% Copyright (C) 2003-2011 Dynare Team
+% Copyright (C) 2003-2013 Dynare Team
 %
 % This file is part of Dynare.
 %
@@ -41,11 +41,11 @@ ncx = size(estim_params_.corrx,1);
 ncn = size(estim_params_.corrn,1);
 np = size(estim_params_.param_vals,1);
 
-estim_params_.nvx = nvx;
+estim_params_.nvx = nvx; %exogenous shock variances
-estim_params_.nvn = nvn;
+estim_params_.nvn = nvn; %endogenous variances, i.e. measurement error
-estim_params_.ncx = ncx;
+estim_params_.ncx = ncx; %exogenous shock correlations
-estim_params_.ncn = ncn;
+estim_params_.ncn = ncn; % correlation between endogenous variables, i.e. measurement error.
-estim_params_.np = np;
+estim_params_.np = np;   % other parameters of the model
 
 xparam1 = [];
 ub = [];
@@ -74,6 +74,7 @@ if nvx
     bayestopt_.name = cellstr(M_.exo_names(estim_params_.var_exo(:,1),:));
 end
 if nvn
+    estim_params_.nvn_observable_correspondence=NaN(nvn,1); % stores number of corresponding observable
     if isequal(M_.H,0)
         nvarobs = size(options_.varobs,1);
         M_.H = zeros(nvarobs,nvarobs);
@@ -83,7 +84,7 @@ if nvn
         if isempty(obsi_)
             error(['The variable ' deblank(M_.endo_names(estim_params_.var_endo(i,1),:)) ' has to be declared as observable since you assume a measurement error on it.'])
         end
-        estim_params_.var_endo(i,1) = obsi_;
+        estim_params_.nvn_observable_correspondence(i,1)=obsi_;
     end
     xparam1 = [xparam1; estim_params_.var_endo(:,2)];
     ub = [ub; estim_params_.var_endo(:,4)]; 
@@ -94,7 +95,7 @@ if nvn
     bayestopt_.p3 = [ bayestopt_.p3; estim_params_.var_endo(:,8)];
     bayestopt_.p4 = [ bayestopt_.p4; estim_params_.var_endo(:,9)];
     bayestopt_.jscale = [ bayestopt_.jscale; estim_params_.var_endo(:,10)];
-    bayestopt_.name = [ bayestopt_.name; cellstr(options_.varobs(estim_params_.var_endo(:,1),:))];
+    bayestopt_.name = [ bayestopt_.name; cellstr(options_.varobs(estim_params_.nvn_observable_correspondence,:))];
 end
 if ncx
     xparam1 = [xparam1; estim_params_.corrx(:,3)];
@@ -111,6 +112,7 @@ if ncx
                         repmat(', ',ncx,1) , deblank(M_.exo_names(estim_params_.corrx(:,2),:))])];
 end
 if ncn
+    estim_params_.corrn_observable_correspondence=NaN(ncn,2);
     if isequal(M_.H,0)
         nvarobs = size(options_.varobs,1);
         M_.H = zeros(nvarobs,nvarobs);
@@ -125,8 +127,15 @@ if ncn
     bayestopt_.p4 = [ bayestopt_.p4; estim_params_.corrn(:,10)];
     bayestopt_.jscale = [ bayestopt_.jscale; estim_params_.corrn(:,11)];
     bayestopt_.name = [bayestopt_.name; cellstr([repmat('corr ',ncn,1) ...
-                        deblank(M_.exo_names(estim_params_.corrn(:,1),:)) ...
+                        deblank(M_.endo_names(estim_params_.corrn(:,1),:)) ...
-                        repmat(', ',ncn,1) , deblank(M_.exo_names(estim_params_.corrn(:,2),:))])];
+                        repmat(', ',ncn,1) , deblank(M_.endo_names(estim_params_.corrn(:,2),:))])];
+    for i=1:ncn
+        k1 = estim_params_.corrn(i,1);
+        k2 = estim_params_.corrn(i,2);
+        obsi1 = strmatch(deblank(M_.endo_names(k1,:)),deblank(options_.varobs),'exact'); %find correspondence to varobs to construct H in set_all_paramters
+        obsi2 = strmatch(deblank(M_.endo_names(k2,:)),deblank(options_.varobs),'exact');
+        estim_params_.corrn_observable_correspondence(i,:)=[obsi1,obsi2]; %save correspondence
+    end
 end
 if np
     xparam1 = [xparam1; estim_params_.param_vals(:,2)];

tests/measurement_errors/fs2000_corr_me_ml_mcmc/fs2000_corr_ME.mod

@@ -0,0 +1,125 @@
+/*
+ * This file is based on the cash in advance model described
+ * Frank Schorfheide (2000): "Loss function-based evaluation of DSGE models",
+ * Journal of Applied Econometrics, 15(6), 645-670.
+ *
+ * The equations are taken from J. Nason and T. Cogley (1994): "Testing the
+ * implications of long-run neutrality for monetary business cycle models",
+ * Journal of Applied Econometrics, 9, S37-S70.
+ * Note that there is an initial minus sign missing in equation (A1), p. S63.
+ *
+ * This implementation was written by Michel Juillard. Please note that the
+ * following copyright notice only applies to this Dynare implementation of the
+ * model.
+ */
+
+/*
+ * Copyright (C) 2004-2013 Dynare Team
+ *
+ * This file is part of Dynare.
+ *
+ * Dynare is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * Dynare is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with Dynare.  If not, see <http://www.gnu.org/licenses/>.
+ */
+
+var m P c e W R k d n l gy_obs gp_obs y dA;
+varexo e_a e_m;
+
+parameters alp bet gam mst rho psi del theta;
+
+alp = 0.33;
+bet = 0.99;
+gam = 0.003;
+mst = 1.011;
+rho = 0.7;
+psi = 0.787;
+del = 0.02;
+theta=0;
+
+model;
+dA = exp(gam+e_a);
+log(m) = (1-rho)*log(mst) + rho*log(m(-1))+e_m;
+-P/(c(+1)*P(+1)*m)+bet*P(+1)*(alp*exp(-alp*(gam+log(e(+1))))*k^(alp-1)*n(+1)^(1-alp)+(1-del)*exp(-(gam+log(e(+1)))))/(c(+2)*P(+2)*m(+1))=0;
+W = l/n;
+-(psi/(1-psi))*(c*P/(1-n))+l/n = 0;
+R = P*(1-alp)*exp(-alp*(gam+e_a))*k(-1)^alp*n^(-alp)/W;
+1/(c*P)-bet*P*(1-alp)*exp(-alp*(gam+e_a))*k(-1)^alp*n^(1-alp)/(m*l*c(+1)*P(+1)) = 0;
+c+k = exp(-alp*(gam+e_a))*k(-1)^alp*n^(1-alp)+(1-del)*exp(-(gam+e_a))*k(-1);
+P*c = m;
+m-1+d = l;
+e = exp(e_a);
+y = k(-1)^alp*n^(1-alp)*exp(-alp*(gam+e_a));
+gy_obs = dA*y/y(-1);
+gp_obs = (P/P(-1))*m(-1)/dA;
+end;
+
+initval;
+k = 6;
+m = mst;
+P = 2.25;
+c = 0.45;
+e = 1;
+W = 4;
+R = 1.02;
+d = 0.85;
+n = 0.19;
+l = 0.86;
+y = 0.6;
+gy_obs = exp(gam);
+gp_obs = exp(-gam);
+dA = exp(gam);
+end;
+
+varobs gp_obs gy_obs;
+
+shocks;
+var e_a; stderr 0.014;
+var e_m; stderr 0.005;
+corr gy_obs,gp_obs = 0.5;
+end;
+
+steady;
+
+
+estimated_params;
+alp, 0.356;
+gam,  0.0085;
+del, 0.01;
+stderr e_a, 0.035449;
+stderr e_m, 0.008862;
+corr e_m, e_a, 0;
+stderr gp_obs, 1;
+stderr gy_obs, 1;
+corr gp_obs, gy_obs,0;
+end;
+
+options_.TeX=1;
+estimation(mode_compute=9,order=1,datafile=fsdat_simul,mode_check,smoother,filter_decomposition,forecast = 8,filtered_vars,filter_step_ahead=[1,3],irf=20) m P c e W R k d y gy_obs;
+
+
+
+estimated_params;
+//alp, beta_pdf, 0.356, 0.02;
+gam, normal_pdf, 0.0085, 0.003;
+//del, beta_pdf, 0.01, 0.005;
+stderr e_a, inv_gamma_pdf, 0.035449, inf;
+stderr e_m, inv_gamma_pdf, 0.008862, inf;
+corr e_m, e_a, normal_pdf, 0, 0.2;
+stderr gp_obs, inv_gamma_pdf, 0.001, inf;
+//stderr gy_obs, inv_gamma_pdf, 0.001, inf;
+//corr gp_obs, gy_obs,normal_pdf, 0, 0.2;
+end;
+
+estimation(mode_compute=0,mode_file=fs2000_corr_ME_mh_mode,order=1,datafile=fsdat_simul,mode_check,smoother,filter_decomposition,mh_replic=2000, mh_nblocks=2, mh_jscale=0.8,forecast = 8,bayesian_irf,filtered_vars,filter_step_ahead=[1,3],irf=20) m P c e W R k d y;
+shock_decomposition y W R;
+//identification(advanced=1,max_dim_cova_group=3,prior_mc=250);