stepan
/
dynare
1
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

Fixed bug in dynare_resolve (wrong calling sequence introduced in commit #013c599ec92f7d6e5fc3f351a58d9aa5ba401410). 

Removed globals from DsgeVarLikelihood and changed the calling sequence. As in DsgeLikelihood, the penalty is now a
persistent variable.

Added a global structure for the data: dataset_.

Removed globals from dsgevar_posterior_density and mode_check.

Simplification of the clode, definition of the variable objective_function at the top of dynare_estimation_1 (equal
to 'DsgeLikelihood' or 'DsgeVarLikelihood').
time-shift
Stéphane Adjemian (Charybdis) 2011-09-22 11:17:31 +02:00
parent f1ffeb29bb
commit e0fa737cee
12 changed files with 215 additions and 206 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

2
matlab/DsgeLikelihood.m
View File

@ -255,7 +255,7 @@ Model.H = H;
%------------------------------------------------------------------------------ %------------------------------------------------------------------------------
% Linearize the model around the deterministic sdteadystate and extract the matrices of the state equation (T and R). % Linearize the model around the deterministic sdteadystate and extract the matrices of the state equation (T and R).
[T,R,SteadyState,info,Model,DynareOptions,DynareResults] = dynare_resolve(Model,DynareOptions,DynareResults); [T,R,SteadyState,info,Model,DynareOptions,DynareResults] = dynare_resolve(Model,DynareOptions,DynareResults,'restrict');
% Return, with endogenous penalty when possible, if dynare_resolve issues an error code (defined in resol). % Return, with endogenous penalty when possible, if dynare_resolve issues an error code (defined in resol).
if info(1) == 1 || info(1) == 2 || info(1) == 5 || info(1) == 22 || info(1) == 24 if info(1) == 1 || info(1) == 2 || info(1) == 5 || info(1) == 22 || info(1) == 24

228
matlab/DsgeVarLikelihood.m
View File

@ -1,18 +1,18 @@
function [fval,cost_flag,info,PHI,SIGMAu,iXX,prior] = DsgeVarLikelihood(xparam1,gend) function [fval,exit_flag,info,PHI,SIGMAu,iXX,prior] = DsgeVarLikelihood(xparam1,DynareDataset,DynareOptions,Model,EstimatedParameters,BayesInfo,DynareResults)
% Evaluates the posterior kernel of the bvar-dsge model. % Evaluates the posterior kernel of the bvar-dsge model.
% %
% INPUTS % INPUTS
% o xparam1 [double] Vector of model's parameters. % o xparam1 [double] Vector of model's parameters.
% o gend [integer] Number of observations (without conditionning observations for the lags). % o gend [integer] Number of observations (without conditionning observations for the lags).
% %
% OUTPUTS % OUTPUTS
% o fval [double] Value of the posterior kernel at xparam1. % o fval [double] Value of the posterior kernel at xparam1.
% o cost_flag [integer] Zero if the function returns a penalty, one otherwise. % o cost_flag [integer] Zero if the function returns a penalty, one otherwise.
% o info [integer] Vector of informations about the penalty. % o info [integer] Vector of informations about the penalty.
% o PHI [double] Stacked BVAR-DSGE autoregressive matrices (at the mode associated to xparam1). % o PHI [double] Stacked BVAR-DSGE autoregressive matrices (at the mode associated to xparam1).
% o SIGMAu [double] Covariance matrix of the BVAR-DSGE (at the mode associated to xparam1). % o SIGMAu [double] Covariance matrix of the BVAR-DSGE (at the mode associated to xparam1).
% o iXX [double] inv(X'X). % o iXX [double] inv(X'X).
% o prior [double] a matlab structure describing the dsge-var prior. % o prior [double] a matlab structure describing the dsge-var prior.
% %
% SPECIAL REQUIREMENTS % SPECIAL REQUIREMENTS
% None. % None.
@ -34,139 +34,180 @@ function [fval,cost_flag,info,PHI,SIGMAu,iXX,prior] = DsgeVarLikelihood(xparam1,
% You should have received a copy of the GNU General Public License % You should have received a copy of the GNU General Public License
% along with Dynare. If not, see <http://www.gnu.org/licenses/>. % along with Dynare. If not, see <http://www.gnu.org/licenses/>.
global bayestopt_ estim_params_ M_ options_ oo_ % Declaration of the persistent variables.
persistent penalty dsge_prior_weight_idx
nvx = estim_params_.nvx; % Initialization of the penalty
nvn = estim_params_.nvn; if ~nargin || isempty(penalty)
ncx = estim_params_.ncx; penalty = 1e8;
ncn = estim_params_.ncn; if ~nargin, return, end
np = estim_params_.np; end
nx = nvx+nvn+ncx+ncn+np; if nargin==1
ns = nvx+nvn+ncx+ncn; penalty = xparam1;
return
end
NumberOfObservedVariables = size(options_.varobs,1); % Initialization of of the index for parameter dsge_prior_weight in Model.params.
NumberOfLags = options_.dsge_varlag; if isempty(dsge_prior_weight_idx)
dsge_prior_weight_idx = strmatch('dsge_prior_weight',Model.param_names);
end
% Get the number of estimated (dsge) parameters.
ns = EstimatedParameters.nvx + ...
EstimatedParameters.nvn + ...
EstimatedParameters.ncx + ...
EstimatedParameters.ncn;
nx = ns + EstimatedParameters.np;
% Get the number of observed variables in the VAR model.
NumberOfObservedVariables = DynareDataset.info.nvobs;
% Get the number of lags in the VAR model.
NumberOfLags = DynareOptions.dsge_varlag;
% Get the number of parameters in the VAR model.
NumberOfParameters = NumberOfObservedVariables*NumberOfLags ; NumberOfParameters = NumberOfObservedVariables*NumberOfLags ;
if ~options_.noconstant if ~DynareOptions.noconstant
NumberOfParameters = NumberOfParameters + 1; NumberOfParameters = NumberOfParameters + 1;
end end
% Get empirical second order moments for the observed variables.
mYY = evalin('base', 'mYY'); mYY = evalin('base', 'mYY');
mYX = evalin('base', 'mYX'); mYX = evalin('base', 'mYX');
mXY = evalin('base', 'mXY'); mXY = evalin('base', 'mXY');
mXX = evalin('base', 'mXX'); mXX = evalin('base', 'mXX');
% Initialize some of the output arguments.
fval = []; fval = [];
cost_flag = 1; exit_flag = 1;
if ~isequal(options_.mode_compute,1) && any(xparam1 < bayestopt_.lb) % Return, with endogenous penalty, if some dsge-parameters are smaller than the lower bound of the prior domain.
k = find(xparam1 < bayestopt_.lb); if DynareOptions.mode_compute ~= 1 && any(xparam1 < BayesInfo.lb)
fval = bayestopt_.penalty+sum((bayestopt_.lb(k)-xparam1(k)).^2); k = find(xparam1 < BayesInfo.lb);
cost_flag = 0; fval = penalty+sum((BayesInfo.lb(k)-xparam1(k)).^2);
exit_flag = 0;
info = 41; info = 41;
return; return;
end end
if ~isequal(options_.mode_compute,11) && any(xparam1 > bayestopt_.ub) % Return, with endogenous penalty, if some dsge-parameters are greater than the upper bound of the prior domain.
k = find(xparam1 > bayestopt_.ub); if DynareOptions.mode_compute ~= 1 && any(xparam1 > BayesInfo.ub)
fval = bayestopt_.penalty+sum((xparam1(k)-bayestopt_.ub(k)).^2); k = find(xparam1 > BayesInfo.ub);
cost_flag = 0; fval = penalty+sum((xparam1(k)-BayesInfo.ub(k)).^2);
exit_flag = 0;
info = 42; info = 42;
return; return;
end end
Q = M_.Sigma_e; % Get the variance of each structural innovation.
for i=1:estim_params_.nvx Q = Model.Sigma_e;
k = estim_params_.var_exo(i,1); for i=1:EstimatedParameters.nvx
k = EstimatedParameters.var_exo(i,1);
Q(k,k) = xparam1(i)*xparam1(i); Q(k,k) = xparam1(i)*xparam1(i);
end end
offset = estim_params_.nvx; offset = EstimatedParameters.nvx;
if estim_params_.nvn
% Check that the user does not estimate measurment errors.
% TODO Check that the user does not declare non estimated measurement errors...
if EstimatedParameters.nvn
disp('DsgeVarLikelihood :: Measurement errors are not implemented!') disp('DsgeVarLikelihood :: Measurement errors are not implemented!')
return return
end end
if estim_params_.ncx
% Check that the user does not estimate off diagonal elements in the covariance matrix of the structural innovation.
% TODO Check that Q is a diagonal matrix...
if EstimatedParameters.ncx
disp('DsgeVarLikelihood :: Correlated structural innovations are not implemented!') disp('DsgeVarLikelihood :: Correlated structural innovations are not implemented!')
return return
end end
M_.params(estim_params_.param_vals(:,1)) = xparam1(offset+1:end); % Update Model.params and Model.Sigma_e.
M_.Sigma_e = Q; Model.params(EstimatedParameters.param_vals(:,1)) = xparam1(offset+1:end);
Model.Sigma_e = Q;
%% Weight of the dsge prior: % Get the weight of the dsge prior.
dsge_prior_weight = M_.params(strmatch('dsge_prior_weight',M_.param_names)); dsge_prior_weight = Model.params(dsge_prior_weight_idx);
% Is the DSGE prior proper?
if dsge_prior_weight<(NumberOfParameters+NumberOfObservedVariables)/gend; % Is the dsge prior proper?
fval = bayestopt_.penalty+abs(gend*dsge_prior_weight-(NumberOfParameters+NumberOfObservedVariables)); if dsge_prior_weight<(NumberOfParameters+NumberOfObservedVariables)/DynareDataset.info.ntobs;
cost_flag = 0; fval = penalty+abs(DynareDataset.info.ntobs*dsge_prior_weight-(NumberOfParameters+NumberOfObservedVariables));
exit_flag = 0;
info = 51; info = 51;
return; return
end end
%------------------------------------------------------------------------------ %------------------------------------------------------------------------------
% 2. call model setup & reduction program % 2. call model setup & reduction program
%------------------------------------------------------------------------------ %------------------------------------------------------------------------------
[T,R,SteadyState,info,M_,options_,oo_] = dynare_resolve(M_,options_,oo_);
% Solve the Dsge model and get the matrices of the reduced form solution. T and R are the matrices of the
% state equation
[T,R,SteadyState,info,Model,DynareOptions,DynareResults] = dynare_resolve(Model,DynareOptions,DynareResults,'restrict');
% Return, with endogenous penalty when possible, if dynare_resolve issues an error code (defined in resol).
if info(1) == 1 || info(1) == 2 || info(1) == 5 if info(1) == 1 || info(1) == 2 || info(1) == 5
fval = bayestopt_.penalty+1; fval = penalty+1;
cost_flag = 0; info = info(1);
exit_flag = 0;
return return
elseif info(1) == 3 || info(1) == 4 || info(1) == 19 || info(1) == 20 || info(1) == 21 elseif info(1) == 3 || info(1) == 4 || info(1) == 19 || info(1) == 20 || info(1) == 21
fval = bayestopt_.penalty+info(2); fval = penalty+info(2);
cost_flag = 0; info = info(1);
exit_flag = 0;
return return
end end
if ~options_.noconstant % Define the mean/steady state vector.
if options_.loglinear if ~DynareOptions.noconstant
constant = transpose(log(SteadyState(bayestopt_.mfys))); if DynareOptions.loglinear
constant = transpose(log(SteadyState(BayesInfo.mfys)));
else else
constant = transpose(SteadyState(bayestopt_.mfys)); constant = transpose(SteadyState(BayesInfo.mfys));
end end
else else
constant = zeros(1,NumberOfObservedVariables); constant = zeros(1,NumberOfObservedVariables);
end end
if bayestopt_.with_trend == 1
disp('DsgeVarLikelihood :: Linear trend is not yet implemented!') % Dsge-VAR with deterministic trends is not implemented
return if BayesInfo.with_trend == 1
error('DsgeVarLikelihood :: Linear trend is not yet implemented!')
end end
%------------------------------------------------------------------------------ %------------------------------------------------------------------------------
% 3. theoretical moments (second order) % 3. theoretical moments (second order)
%------------------------------------------------------------------------------ %------------------------------------------------------------------------------
tmp0 = lyapunov_symm(T,R*Q*R',options_.qz_criterium,options_.lyapunov_complex_threshold);% I compute the variance-covariance matrix
mf = bayestopt_.mf1; % of the restricted state vector. % Compute the theoretical second order moments
tmp0 = lyapunov_symm(T,R*Q*R',DynareOptions.qz_criterium,DynareOptions.lyapunov_complex_threshold);
mf = BayesInfo.mf1;
% Get the non centered second order moments % Get the non centered second order moments
TheoreticalAutoCovarianceOfTheObservedVariables = ... TheoreticalAutoCovarianceOfTheObservedVariables = zeros(NumberOfObservedVariables,NumberOfObservedVariables,NumberOfLags+1);
zeros(NumberOfObservedVariables,NumberOfObservedVariables,NumberOfLags+1);
TheoreticalAutoCovarianceOfTheObservedVariables(:,:,1) = tmp0(mf,mf)+constant'*constant; TheoreticalAutoCovarianceOfTheObservedVariables(:,:,1) = tmp0(mf,mf)+constant'*constant;
for lag = 1:NumberOfLags for lag = 1:NumberOfLags
tmp0 = T*tmp0; tmp0 = T*tmp0;
TheoreticalAutoCovarianceOfTheObservedVariables(:,:,lag+1) = tmp0(mf,mf) ... TheoreticalAutoCovarianceOfTheObservedVariables(:,:,lag+1) = tmp0(mf,mf) + constant'*constant;
+ constant'*constant;
end end
% Build the theoretical "covariance" between Y and X % Build the theoretical "covariance" between Y and X
GYX = zeros(NumberOfObservedVariables,NumberOfParameters); GYX = zeros(NumberOfObservedVariables,NumberOfParameters);
for i=1:NumberOfLags for i=1:NumberOfLags
GYX(:,(i-1)*NumberOfObservedVariables+1:i*NumberOfObservedVariables) = ... GYX(:,(i-1)*NumberOfObservedVariables+1:i*NumberOfObservedVariables) = TheoreticalAutoCovarianceOfTheObservedVariables(:,:,i+1);
TheoreticalAutoCovarianceOfTheObservedVariables(:,:,i+1);
end end
if ~options_.noconstant if ~DynareOptions.noconstant
GYX(:,end) = constant'; GYX(:,end) = constant';
end end
% Build the theoretical "covariance" between X and X % Build the theoretical "covariance" between X and X
GXX = kron(eye(NumberOfLags), ... GXX = kron(eye(NumberOfLags), TheoreticalAutoCovarianceOfTheObservedVariables(:,:,1));
TheoreticalAutoCovarianceOfTheObservedVariables(:,:,1));
for i = 1:NumberOfLags-1 for i = 1:NumberOfLags-1
tmp1 = diag(ones(NumberOfLags-i,1),i); tmp1 = diag(ones(NumberOfLags-i,1),i);
tmp2 = diag(ones(NumberOfLags-i,1),-i); tmp2 = diag(ones(NumberOfLags-i,1),-i);
GXX = GXX + kron(tmp1,TheoreticalAutoCovarianceOfTheObservedVariables(:,:,i+1)); GXX = GXX + kron(tmp1,TheoreticalAutoCovarianceOfTheObservedVariables(:,:,i+1));
GXX = GXX + kron(tmp2,TheoreticalAutoCovarianceOfTheObservedVariables(:,:,i+1)'); GXX = GXX + kron(tmp2,TheoreticalAutoCovarianceOfTheObservedVariables(:,:,i+1)');
end end
if ~options_.noconstant if ~DynareOptions.noconstant
% Add one row and one column to GXX % Add one row and one column to GXX
GXX = [GXX , kron(ones(NumberOfLags,1),constant') ; [ kron(ones(1,NumberOfLags),constant) , 1] ]; GXX = [GXX , kron(ones(NumberOfLags,1),constant') ; [ kron(ones(1,NumberOfLags),constant) , 1] ];
end end
@ -177,45 +218,46 @@ assignin('base','GYY',GYY);
assignin('base','GXX',GXX); assignin('base','GXX',GXX);
assignin('base','GYX',GYX); assignin('base','GYX',GYX);
if ~isinf(dsge_prior_weight) if ~isinf(dsge_prior_weight)% Evaluation of the likelihood of the dsge-var model when the dsge prior weight is finite.
tmp0 = dsge_prior_weight*gend*TheoreticalAutoCovarianceOfTheObservedVariables(:,:,1) + mYY ; tmp0 = dsge_prior_weight*DynareDataset.info.ntobs*TheoreticalAutoCovarianceOfTheObservedVariables(:,:,1) + mYY ;
tmp1 = dsge_prior_weight*gend*GYX + mYX; tmp1 = dsge_prior_weight*DynareDataset.info.ntobs*GYX + mYX;
tmp2 = inv(dsge_prior_weight*gend*GXX+mXX); tmp2 = inv(dsge_prior_weight*DynareDataset.info.ntobs*GXX+mXX);
SIGMAu = tmp0 - tmp1*tmp2*tmp1'; clear('tmp0'); SIGMAu = tmp0 - tmp1*tmp2*tmp1'; clear('tmp0');
if ~ispd(SIGMAu) if ~ispd(SIGMAu)
v = diag(SIGMAu); v = diag(SIGMAu);
k = find(v<0); k = find(v<0);
fval = bayestopt_.penalty + sum(v(k).^2); fval = penalty + sum(v(k).^2);
info = 52; info = 52;
cost_flag = 0; exit_flag = 0;
return; return;
end end
SIGMAu = SIGMAu / (gend*(1+dsge_prior_weight)); SIGMAu = SIGMAu / (DynareDataset.info.ntobs*(1+dsge_prior_weight));
PHI = tmp2*tmp1'; clear('tmp1'); PHI = tmp2*tmp1'; clear('tmp1');
prodlng1 = sum(gammaln(.5*((1+dsge_prior_weight)*gend- ... prodlng1 = sum(gammaln(.5*((1+dsge_prior_weight)*DynareDataset.info.ntobs- ...
NumberOfObservedVariables*NumberOfLags ... NumberOfObservedVariables*NumberOfLags ...
+1-(1:NumberOfObservedVariables)'))); +1-(1:NumberOfObservedVariables)')));
prodlng2 = sum(gammaln(.5*(dsge_prior_weight*gend- ... prodlng2 = sum(gammaln(.5*(dsge_prior_weight*DynareDataset.info.ntobs- ...
NumberOfObservedVariables*NumberOfLags ... NumberOfObservedVariables*NumberOfLags ...
+1-(1:NumberOfObservedVariables)'))); +1-(1:NumberOfObservedVariables)')));
lik = .5*NumberOfObservedVariables*log(det(dsge_prior_weight*gend*GXX+mXX)) ... lik = .5*NumberOfObservedVariables*log(det(dsge_prior_weight*DynareDataset.info.ntobs*GXX+mXX)) ...
+ .5*((dsge_prior_weight+1)*gend-NumberOfParameters)*log(det((dsge_prior_weight+1)*gend*SIGMAu)) ... + .5*((dsge_prior_weight+1)*DynareDataset.info.ntobs-NumberOfParameters)*log(det((dsge_prior_weight+1)*DynareDataset.info.ntobs*SIGMAu)) ...
- .5*NumberOfObservedVariables*log(det(dsge_prior_weight*gend*GXX)) ... - .5*NumberOfObservedVariables*log(det(dsge_prior_weight*DynareDataset.info.ntobs*GXX)) ...
- .5*(dsge_prior_weight*gend-NumberOfParameters)*log(det(dsge_prior_weight*gend*(GYY-GYX*inv(GXX)*GYX'))) ... - .5*(dsge_prior_weight*DynareDataset.info.ntobs-NumberOfParameters)*log(det(dsge_prior_weight*DynareDataset.info.ntobs*(GYY-GYX*inv(GXX)*GYX'))) ...
+ .5*NumberOfObservedVariables*gend*log(2*pi) ... + .5*NumberOfObservedVariables*DynareDataset.info.ntobs*log(2*pi) ...
- .5*log(2)*NumberOfObservedVariables*((dsge_prior_weight+1)*gend-NumberOfParameters) ... - .5*log(2)*NumberOfObservedVariables*((dsge_prior_weight+1)*DynareDataset.info.ntobs-NumberOfParameters) ...
+ .5*log(2)*NumberOfObservedVariables*(dsge_prior_weight*gend-NumberOfParameters) ... + .5*log(2)*NumberOfObservedVariables*(dsge_prior_weight*DynareDataset.info.ntobs-NumberOfParameters) ...
- prodlng1 + prodlng2; - prodlng1 + prodlng2;
else else% Evaluation of the likelihood of the dsge-var model when the dsge prior weight is infinite.
iGXX = inv(GXX); iGXX = inv(GXX);
SIGMAu = GYY - GYX*iGXX*transpose(GYX); SIGMAu = GYY - GYX*iGXX*transpose(GYX);
PHI = iGXX*transpose(GYX); PHI = iGXX*transpose(GYX);
lik = gend * ( log(det(SIGMAu)) + NumberOfObservedVariables*log(2*pi) + ... lik = DynareDataset.info.ntobs * ( log(det(SIGMAu)) + NumberOfObservedVariables*log(2*pi) + ...
trace(inv(SIGMAu)*(mYY - transpose(mYX*PHI) - mYX*PHI + transpose(PHI)*mXX*PHI)/gend)); trace(inv(SIGMAu)*(mYY - transpose(mYX*PHI) - mYX*PHI + transpose(PHI)*mXX*PHI)/DynareDataset.info.ntobs));
lik = .5*lik;% Minus likelihood lik = .5*lik;% Minus likelihood
end end
lnprior = priordens(xparam1,bayestopt_.pshape,bayestopt_.p6,bayestopt_.p7,bayestopt_.p3,bayestopt_.p4); % Add the (logged) prior density for the dsge-parameters.
lnprior = priordens(xparam1,BayesInfo.pshape,BayesInfo.p6,BayesInfo.p7,BayesInfo.p3,BayesInfo.p4);
fval = (lik-lnprior); fval = (lik-lnprior);
if (nargout == 6) if (nargout == 6)
@ -232,10 +274,10 @@ if (nargout==7)
else else
iXX = tmp2; iXX = tmp2;
end end
iGXX = inv(GXX); iGXX = inv(GXX);
prior.SIGMAstar = GYY - GYX*iGXX*GYX'; prior.SIGMAstar = GYY - GYX*iGXX*GYX';
prior.PHIstar = iGXX*transpose(GYX); prior.PHIstar = iGXX*transpose(GYX);
prior.ArtificialSampleSize = fix(dsge_prior_weight*gend); prior.ArtificialSampleSize = fix(dsge_prior_weight*DynareDataset.info.ntobs);
prior.DF = prior.ArtificialSampleSize - NumberOfParameters - NumberOfObservedVariables; prior.DF = prior.ArtificialSampleSize - NumberOfParameters - NumberOfObservedVariables;
prior.iGXX = iGXX; prior.iGXX = iGXX;
end end

9
matlab/PosteriorIRF.m
View File

@ -34,7 +34,7 @@ function PosteriorIRF(type)
% along with Dynare. If not, see <http://www.gnu.org/licenses/>. % along with Dynare. If not, see <http://www.gnu.org/licenses/>.
global options_ estim_params_ oo_ M_ bayestopt_ global options_ estim_params_ oo_ M_ bayestopt_ dataset_
% Set the number of periods % Set the number of periods
if isempty(options_.irf) || ~options_.irf if isempty(options_.irf) || ~options_.irf
options_.irf = 40; options_.irf = 40;
@ -64,8 +64,8 @@ np = estim_params_.np ;
npar = nvx+nvn+ncx+ncn+np; npar = nvx+nvn+ncx+ncn+np;
offset = npar-np; clear('nvx','nvn','ncx','ncn','np'); offset = npar-np; clear('nvx','nvn','ncx','ncn','np');
nvobs = size(options_.varobs,1); nvobs = dataset_.info.nvobs;
gend = options_.nobs; gend = dataset_.info.ntobs;
MaxNumberOfPlotPerFigure = 9; MaxNumberOfPlotPerFigure = 9;
nn = sqrt(MaxNumberOfPlotPerFigure); nn = sqrt(MaxNumberOfPlotPerFigure);
MAX_nirfs_dsge = ceil(options_.MaxNumberOfBytes/(options_.irf*nvar*M_.exo_nbr)/8); MAX_nirfs_dsge = ceil(options_.MaxNumberOfBytes/(options_.irf*nvar*M_.exo_nbr)/8);
@ -230,7 +230,8 @@ else
'options_', options_, ... 'options_', options_, ...
'bayestopt_', bayestopt_, ... 'bayestopt_', bayestopt_, ...
'estim_params_', estim_params_, ... 'estim_params_', estim_params_, ...
'oo_', oo_); 'oo_', oo_, ...
'dataset_',dataset_);
% which files have to be copied to run remotely % which files have to be copied to run remotely
NamFileInput(1,:) = {'',[M_.fname '_static.m']}; NamFileInput(1,:) = {'',[M_.fname '_static.m']};

31
matlab/PosteriorIRF_core1.m
View File

@ -40,7 +40,7 @@ function myoutput=PosteriorIRF_core1(myinputs,fpar,npar,whoiam, ThisMatlab)
% along with Dynare. If not, see <http://www.gnu.org/licenses/>. % along with Dynare. If not, see <http://www.gnu.org/licenses/>.
global options_ estim_params_ oo_ M_ bayestopt_ global options_ estim_params_ oo_ M_ bayestopt_ dataset_
if nargin<4, if nargin<4,
whoiam=0; whoiam=0;
@ -151,6 +151,7 @@ while fpar<npar
stock_param(irun2,:) = deep; stock_param(irun2,:) = deep;
set_parameters(deep); set_parameters(deep);
[dr,info,M_,options_,oo_] = resol(0,M_,options_,oo_); [dr,info,M_,options_,oo_] = resol(0,M_,options_,oo_);
oo_.dr = dr;
if info(1) if info(1)
nosaddle = nosaddle + 1; nosaddle = nosaddle + 1;
fpar = fpar - 1; fpar = fpar - 1;
@ -188,24 +189,24 @@ while fpar<npar
end end
if MAX_nirfs_dsgevar if MAX_nirfs_dsgevar
IRUN = IRUN+1; IRUN = IRUN+1;
[fval,cost_flag,info,PHI,SIGMAu,iXX] = DsgeVarLikelihood(deep',gend); [fval,cost_flag,info,PHI,SIGMAu,iXX] = DsgeVarLikelihood(deep',dataset_,options_,M_,estim_params_,bayestopt_,oo_);
dsge_prior_weight = M_.params(strmatch('dsge_prior_weight',M_.param_names)); dsge_prior_weight = M_.params(strmatch('dsge_prior_weight',M_.param_names));
DSGE_PRIOR_WEIGHT = floor(gend*(1+dsge_prior_weight)); DSGE_PRIOR_WEIGHT = floor(dataset_.info.ntobs*(1+dsge_prior_weight));
SIGMA_inv_upper_chol = chol(inv(SIGMAu*gend*(dsge_prior_weight+1))); SIGMA_inv_upper_chol = chol(inv(SIGMAu*dataset_.info.ntobs*(dsge_prior_weight+1)));
explosive_var = 1; explosive_var = 1;
while explosive_var while explosive_var
% draw from the marginal posterior of SIGMA % draw from the marginal posterior of SIGMA
SIGMAu_draw = rand_inverse_wishart(nvobs, DSGE_PRIOR_WEIGHT-NumberOfParametersPerEquation, ... SIGMAu_draw = rand_inverse_wishart(dataset_.info.nvobs, DSGE_PRIOR_WEIGHT-NumberOfParametersPerEquation, ...
SIGMA_inv_upper_chol); SIGMA_inv_upper_chol);
% draw from the conditional posterior of PHI % draw from the conditional posterior of PHI
PHI_draw = rand_matrix_normal(NumberOfParametersPerEquation,nvobs, PHI, ... PHI_draw = rand_matrix_normal(NumberOfParametersPerEquation,dataset_.info.nvobs, PHI, ...
chol(SIGMAu_draw)', chol(iXX)'); chol(SIGMAu_draw)', chol(iXX)');
Companion_matrix(1:nvobs,:) = transpose(PHI_draw(1:NumberOfLagsTimesNvobs,:)); Companion_matrix(1:dataset_.info.nvobs,:) = transpose(PHI_draw(1:NumberOfLagsTimesNvobs,:));
% Check for stationarity % Check for stationarity
explosive_var = any(abs(eig(Companion_matrix))>1.000000001); explosive_var = any(abs(eig(Companion_matrix))>1.000000001);
end end
% Get the mean % Get the mean
mu = zeros(1,nvobs); mu = zeros(1,dataset_.info.nvobs);
% Get rotation % Get rotation
if dsge_prior_weight > 0 if dsge_prior_weight > 0
Atheta(oo_.dr.order_var,M_.exo_names_orig_ord) = oo_.dr.ghu*sqrt(M_.Sigma_e); Atheta(oo_.dr.order_var,M_.exo_names_orig_ord) = oo_.dr.ghu*sqrt(M_.Sigma_e);
@ -215,24 +216,24 @@ while fpar<npar
SIGMAu_chol = chol(SIGMAu_draw)'; SIGMAu_chol = chol(SIGMAu_draw)';
SIGMAtrOMEGA = SIGMAu_chol*OMEGAstar'; SIGMAtrOMEGA = SIGMAu_chol*OMEGAstar';
PHIpower = eye(NumberOfLagsTimesNvobs); PHIpower = eye(NumberOfLagsTimesNvobs);
irfs = zeros (options_.irf,nvobs*M_.exo_nbr); irfs = zeros (options_.irf,dataset_.info.nvobs*M_.exo_nbr);
tmp3 = PHIpower(1:nvobs,1:nvobs)*SIGMAtrOMEGA; tmp3 = PHIpower(1:dataset_.info.nvobs,1:dataset_.info.nvobs)*SIGMAtrOMEGA;
irfs(1,:) = tmp3(:)'; irfs(1,:) = tmp3(:)';
for t = 2:options_.irf for t = 2:options_.irf
PHIpower = Companion_matrix*PHIpower; PHIpower = Companion_matrix*PHIpower;
tmp3 = PHIpower(1:nvobs,1:nvobs)*SIGMAtrOMEGA; tmp3 = PHIpower(1:dataset_.info.nvobs,1:dataset_.info.nvobs)*SIGMAtrOMEGA;
irfs(t,:) = tmp3(:)'+kron(ones(1,M_.exo_nbr),mu); irfs(t,:) = tmp3(:)'+kron(ones(1,M_.exo_nbr),mu);
end end
tmp_dsgevar = kron(ones(options_.irf,1),mu); tmp_dsgevar = kron(ones(options_.irf,1),mu);
for j = 1:(nvobs*M_.exo_nbr) for j = 1:(dataset_.info.nvobs*M_.exo_nbr)
if max(irfs(:,j)) - min(irfs(:,j)) > 1e-10 if max(irfs(:,j)) - min(irfs(:,j)) > 1e-10
tmp_dsgevar(:,j) = (irfs(:,j)); tmp_dsgevar(:,j) = (irfs(:,j));
end end
end end
if IRUN < MAX_nirfs_dsgevar if IRUN < MAX_nirfs_dsgevar
stock_irf_bvardsge(:,:,:,IRUN) = reshape(tmp_dsgevar,options_.irf,nvobs,M_.exo_nbr); stock_irf_bvardsge(:,:,:,IRUN) = reshape(tmp_dsgevar,options_.irf,dataset_.info.nvobs,M_.exo_nbr);
else else
stock_irf_bvardsge(:,:,:,IRUN) = reshape(tmp_dsgevar,options_.irf,nvobs,M_.exo_nbr); stock_irf_bvardsge(:,:,:,IRUN) = reshape(tmp_dsgevar,options_.irf,dataset_.info.nvobs,M_.exo_nbr);
instr = [MhDirectoryName '/' M_.fname '_irf_bvardsge' ... instr = [MhDirectoryName '/' M_.fname '_irf_bvardsge' ...
int2str(NumberOfIRFfiles_dsgevar) '.mat stock_irf_bvardsge;'];, int2str(NumberOfIRFfiles_dsgevar) '.mat stock_irf_bvardsge;'];,
eval(['save ' instr]); eval(['save ' instr]);
@ -241,7 +242,7 @@ while fpar<npar
end end
NumberOfIRFfiles_dsgevar = NumberOfIRFfiles_dsgevar+1; NumberOfIRFfiles_dsgevar = NumberOfIRFfiles_dsgevar+1;
IRUN =0; IRUN =0;
stock_irf_dsgevar = zeros(options_.irf,nvobs,M_.exo_nbr,MAX_nirfs_dsgevar); stock_irf_dsgevar = zeros(options_.irf,dataset_.info.nvobs,M_.exo_nbr,MAX_nirfs_dsgevar);
end end
end end
if irun == MAX_nirfs_dsge || irun == npar || fpar == npar if irun == MAX_nirfs_dsge || irun == npar || fpar == npar

34
matlab/dsgevar_posterior_density.m
View File

@ -1,20 +1,20 @@
function bvar = dsgevar_posterior_density(deep) function bvar = dsgevar_posterior_density(deep,DynareDataset,DynareOptions,Model,EstimatedParameters,BayesInfo,DynareResults)
% This function characterizes the posterior distribution of a bvar with % This function characterizes the posterior distribution of a bvar with
% a dsge prior (as in Del Negro and Schorfheide 2003) for a given value % a dsge prior (as in Del Negro and Schorfheide 2003) for a given value
% of the deep parameters (structural parameters + the size of the % of the deep parameters (structural parameters + the size of the
% shocks + dsge_prior_weight). % shocks + dsge_prior_weight).
% %
% INPUTS % INPUTS
% deep: [double] a vector with the deep parameters. % deep: [double] a vector with the deep parameters.
% %
% OUTPUTS % OUTPUTS
% bvar: a matlab structure with prior and posterior densities. % bvar: a matlab structure with prior and posterior densities.
% %
% ALGORITHM % ALGORITHM
% ... % ...
% SPECIAL REQUIREMENTS % SPECIAL REQUIREMENTS
% none % none
% %
% Copyright (C) 1996-2008 Dynare Team % Copyright (C) 1996-2008 Dynare Team
% %
@ -33,8 +33,6 @@ function bvar = dsgevar_posterior_density(deep)
% You should have received a copy of the GNU General Public License % You should have received a copy of the GNU General Public License
% along with Dynare. If not, see <http://www.gnu.org/licenses/>. % along with Dynare. If not, see <http://www.gnu.org/licenses/>.
global options_ M_
gend = options_.nobs; gend = options_.nobs;
dsge_prior_weight = M_.params(strmatch('dsge_prior_weight',M_.param_names)); dsge_prior_weight = M_.params(strmatch('dsge_prior_weight',M_.param_names));
DSGE_PRIOR_WEIGHT = floor(gend*(1+dsge_prior_weight)); DSGE_PRIOR_WEIGHT = floor(gend*(1+dsge_prior_weight));
@ -42,14 +40,14 @@ DSGE_PRIOR_WEIGHT = floor(gend*(1+dsge_prior_weight));
bvar.NumberOfLags = options_.varlag; bvar.NumberOfLags = options_.varlag;
bvar.NumberOfVariables = size(options_.varobs,1); bvar.NumberOfVariables = size(options_.varobs,1);
bvar.Constant = 'no'; bvar.Constant = 'no';
bvar.NumberOfEstimatedParameters = bvar.NumberOfLags*bvar.NumberOfVariables; bvar.NumberOfEstimatedParameters = bvar.NumberOfLags*bvar.NumberOfVariables;
if ~options_.noconstant if ~options_.noconstant
bvar.Constant = 'yes'; bvar.Constant = 'yes';
bvar.NumberOfEstimatedParameters = bvar.NumberOfEstimatedParameters + ... bvar.NumberOfEstimatedParameters = bvar.NumberOfEstimatedParameters + ...
bvar.NumberOfVariables; bvar.NumberOfVariables;
end end
[fval,cost_flag,info,PHI,SIGMAu,iXX,prior] = DsgeVarLikelihood(deep',gend); [fval,cost_flag,info,PHI,SIGMAu,iXX,prior] = DsgeVarLikelihood(deep',DynareDataset,DynareOptions,Model,EstimatedParameters,BayesInfo,DynareResults);
% Conditionnal posterior density of the lagged matrices (given Sigma) -> % Conditionnal posterior density of the lagged matrices (given Sigma) ->
% Matric-variate normal distribution. % Matric-variate normal distribution.
@ -58,12 +56,12 @@ bvar.LaggedMatricesConditionalOnSigma.posterior.arg1 = PHI;
bvar.LaggedMatricesConditionalOnSigma.posterior.arg2 = 'Sigma'; bvar.LaggedMatricesConditionalOnSigma.posterior.arg2 = 'Sigma';
bvar.LaggedMatricesConditionalOnSigma.posterior.arg3 = iXX; bvar.LaggedMatricesConditionalOnSigma.posterior.arg3 = iXX;
% Marginal posterior density of the covariance matrix -> Inverted Wishart. % Marginal posterior density of the covariance matrix -> Inverted Wishart.
bvar.Sigma.posterior.density = 'inverse wishart'; bvar.Sigma.posterior.density = 'inverse wishart';
bvar.Sigma.posterior.arg1 = SIGMAu*DSGE_PRIOR_WEIGHT; bvar.Sigma.posterior.arg1 = SIGMAu*DSGE_PRIOR_WEIGHT;
bvar.Sigma.posterior.arg2 = DSGE_PRIOR_WEIGHT-bvar.NumberOfEstimatedParameters; bvar.Sigma.posterior.arg2 = DSGE_PRIOR_WEIGHT-bvar.NumberOfEstimatedParameters;
% Marginal posterior density of the lagged matrices -> Generalized % Marginal posterior density of the lagged matrices -> Generalized
% Student distribution (See appendix B.5 in Zellner (1971)). % Student distribution (See appendix B.5 in Zellner (1971)).
bvar.LaggedMatrices.posterior.density = 'matric-variate student'; bvar.LaggedMatrices.posterior.density = 'matric-variate student';
bvar.LaggedMatrices.posterior.arg1 = inv(iXX);%P bvar.LaggedMatrices.posterior.arg1 = inv(iXX);%P
@ -80,12 +78,12 @@ bvar.LaggedMatricesConditionalOnSigma.prior.arg1 = prior.PHIstar;
bvar.LaggedMatricesConditionalOnSigma.prior.arg2 = 'Sigma'; bvar.LaggedMatricesConditionalOnSigma.prior.arg2 = 'Sigma';
bvar.LaggedMatricesConditionalOnSigma.prior.arg3 = prior.iGXX; bvar.LaggedMatricesConditionalOnSigma.prior.arg3 = prior.iGXX;
% Marginal posterior density of the covariance matrix -> Inverted Wishart. % Marginal posterior density of the covariance matrix -> Inverted Wishart.
bvar.Sigma.prior.density = 'inverse wishart'; bvar.Sigma.prior.density = 'inverse wishart';
bvar.Sigma.prior.arg1 = prior.SIGMAstar*prior.ArtificialSampleSize; bvar.Sigma.prior.arg1 = prior.SIGMAstar*prior.ArtificialSampleSize;
bvar.Sigma.prior.arg2 = prior.DF; bvar.Sigma.prior.arg2 = prior.DF;
% Marginal posterior density of the lagged matrices -> Generalized % Marginal posterior density of the lagged matrices -> Generalized
% Student distribution (See appendix B.5 in Zellner (1971)). % Student distribution (See appendix B.5 in Zellner (1971)).
bvar.LaggedMatrices.prior.density = 'matric-variate student'; bvar.LaggedMatrices.prior.density = 'matric-variate student';
bvar.LaggedMatrices.prior.arg1 = inv(prior.iGXX);%P bvar.LaggedMatrices.prior.arg1 = inv(prior.iGXX);%P

33
matlab/dynare_estimation_1.m
View File

@ -29,7 +29,7 @@ function dynare_estimation_1(var_list_,dname)
% You should have received a copy of the GNU General Public License % You should have received a copy of the GNU General Public License
% along with Dynare. If not, see <http://www.gnu.org/licenses/>. % along with Dynare. If not, see <http://www.gnu.org/licenses/>.
global M_ options_ oo_ estim_params_ bayestopt_ global M_ options_ oo_ estim_params_ bayestopt_ dataset_
if ~options_.dsge_var if ~options_.dsge_var
objective_function = str2func('DsgeLikelihood'); objective_function = str2func('DsgeLikelihood');
@ -208,11 +208,6 @@ if ~isequal(options_.mode_compute,0) && ~options_.mh_posterior_mode_estimation
else else
nit=1000; nit=1000;
end end
if ~options_.dsge_var
[xparam1,hh,gg,fval,invhess] = newrat('DsgeLikelihood',xparam1,hh,gg,igg,crit,nit,flag,dataset_,options_,M_,estim_params_,bayestopt_,oo_);
else
[xparam1,hh,gg,fval,invhess] = newrat('DsgeVarLikelihood',xparam1,hh,gg,igg,crit,nit,flag,gend);
end
[xparam1,hh,gg,fval,invhess] = newrat(objective_function,xparam1,hh,gg,igg,crit,nit,flag,dataset_,options_,M_,estim_params_,bayestopt_,oo_); [xparam1,hh,gg,fval,invhess] = newrat(objective_function,xparam1,hh,gg,igg,crit,nit,flag,dataset_,options_,M_,estim_params_,bayestopt_,oo_);
parameter_names = bayestopt_.name; parameter_names = bayestopt_.name;
save([M_.fname '_mode.mat'],'xparam1','hh','gg','fval','invhess','parameter_names'); save([M_.fname '_mode.mat'],'xparam1','hh','gg','fval','invhess','parameter_names');
@ -362,7 +357,7 @@ if ~isequal(options_.mode_compute,0) && ~options_.mh_posterior_mode_estimation
end end
if options_.cova_compute == 0 if options_.cova_compute == 0
hh = NaN(length(xparam1),length(xparam1)); hh = [];%NaN(length(xparam1),length(xparam1));
end end
if ~options_.mh_posterior_mode_estimation && options_.cova_compute if ~options_.mh_posterior_mode_estimation && options_.cova_compute
@ -380,11 +375,7 @@ if ~options_.mh_posterior_mode_estimation && options_.cova_compute
end end
if options_.mode_check == 1 && ~options_.mh_posterior_mode_estimation if options_.mode_check == 1 && ~options_.mh_posterior_mode_estimation
if options_.cova_compute mode_check('objective_function',xparam1,hh,dataset_,options_,M_,estim_params_,bayestopt_,oo_);
mode_check(xparam1,0,hh,gend,data,lb,ub,data_index,number_of_observations,no_more_missing_observations);
else
mode_check(xparam1,0,[],gend,data,lb,ub,data_index,number_of_observations,no_more_missing_observations);
end
end end
if ~options_.mh_posterior_mode_estimation if ~options_.mh_posterior_mode_estimation
@ -509,16 +500,10 @@ if any(bayestopt_.pshape > 0) && ~options_.mh_posterior_mode_estimation
estim_params_nbr = size(xparam1,1); estim_params_nbr = size(xparam1,1);
scale_factor = -sum(log10(diag(invhess))); scale_factor = -sum(log10(diag(invhess)));
log_det_invhess = -estim_params_nbr*log(scale_factor)+log(det(scale_factor*invhess)); log_det_invhess = -estim_params_nbr*log(scale_factor)+log(det(scale_factor*invhess));
if ~options_.dsge_var likelihood = feval(objective_function,xparam1,dataset_,options_,M_,estim_params_,bayestopt_,oo_);
md_Laplace = .5*estim_params_nbr*log(2*pi) + .5*log_det_invhess ... oo_.MarginalDensity.LaplaceApproximation = .5*estim_params_nbr*log(2*pi) + .5*log_det_invhess - likelihood;
- DsgeLikelihood(xparam1,dataset_,options_,M_,estim_params_,bayestopt_,oo_);
else
md_Laplace = .5*estim_params_nbr*log(2*pi) + .5*log_det_invhess ...
- DsgeVarLikelihood(xparam1,gend);
end
oo_.MarginalDensity.LaplaceApproximation = md_Laplace;
disp(' ') disp(' ')
disp(sprintf('Log data density [Laplace approximation] is %f.',md_Laplace)) disp(sprintf('Log data density [Laplace approximation] is %f.',oo_.MarginalDensity.LaplaceApproximation))
disp(' ') disp(' ')
end end
elseif ~any(bayestopt_.pshape > 0) && options_.mh_posterior_mode_estimation elseif ~any(bayestopt_.pshape > 0) && options_.mh_posterior_mode_estimation
@ -807,11 +792,7 @@ if (any(bayestopt_.pshape >0 ) && options_.mh_replic) || ...
invhess = compute_mh_covariance_matrix; invhess = compute_mh_covariance_matrix;
end end
if options_.cova_compute if options_.cova_compute
if options_.dsge_var feval(options_.posterior_sampling_method,objective_function,options_.proposal_distribution,xparam1,invhess,bounds,dataset_,options_,M_,estim_params_,bayestopt_,oo_);
feval(options_.posterior_sampling_method,'DsgeVarLikelihood',options_.proposal_distribution,xparam1,invhess,bounds,gend);
else
feval(options_.posterior_sampling_method,'DsgeLikelihood',options_.proposal_distribution,xparam1,invhess,bounds,dataset_,options_,M_,estim_params_,bayestopt_,oo_);
end
else else
error('I Cannot start the MCMC because the hessian of the posterior kernel at the mode was not computed.') error('I Cannot start the MCMC because the hessian of the posterior kernel at the mode was not computed.')
end end

9
matlab/dynare_resolve.m
View File

@ -1,4 +1,4 @@
function [A,B,ys,info,Model,DynareOptions,DynareResults] = dynare_resolve(Model,DynareOptions,DynareResults) function [A,B,ys,info,Model,DynareOptions,DynareResults] = dynare_resolve(Model,DynareOptions,DynareResults,mode)
% Computes the linear approximation and the matrices A and B of the transition equation. % Computes the linear approximation and the matrices A and B of the transition equation.
%@info: %@info:
@ -80,15 +80,18 @@ if info(1) > 0
return return
end end
if nargin == 0 switch nargin
case 3
endo_nbr = Model.endo_nbr; endo_nbr = Model.endo_nbr;
nstatic = DynareResults.dr.nstatic; nstatic = DynareResults.dr.nstatic;
npred = DynareResults.dr.npred; npred = DynareResults.dr.npred;
iv = (1:endo_nbr)'; iv = (1:endo_nbr)';
ic = [ nstatic+(1:npred) endo_nbr+(1:size(DynareResults.dr.ghx,2)-npred) ]'; ic = [ nstatic+(1:npred) endo_nbr+(1:size(DynareResults.dr.ghx,2)-npred) ]';
else case 4
iv = DynareResults.dr.restrict_var_list; iv = DynareResults.dr.restrict_var_list;
ic = DynareResults.dr.restrict_columns; ic = DynareResults.dr.restrict_columns;
otherwise
error('dynare_resolve:: Error in the calling sequence!')
end end
if nargout==1 if nargout==1

4
matlab/gsa/stab_map_.m
View File

@ -247,7 +247,7 @@ if fload==0,
M_.params(estim_params_.param_vals(:,1)) = lpmat(j,:)'; M_.params(estim_params_.param_vals(:,1)) = lpmat(j,:)';
%try stoch_simul([]); %try stoch_simul([]);
try try
[Tt,Rr,SteadyState,infox{j},M_,options_,oo_] = dynare_resolve(M_,options_,oo_); [Tt,Rr,SteadyState,infox{j},M_,options_,oo_] = dynare_resolve(M_,options_,oo_,'restrict');
if infox{j}(1)==0 && ~exist('T'), if infox{j}(1)==0 && ~exist('T'),
dr_=oo_.dr; dr_=oo_.dr;
T=zeros(size(dr_.ghx,1),size(dr_.ghx,2)+size(dr_.ghu,2),Nsam); T=zeros(size(dr_.ghx,1),size(dr_.ghx,2)+size(dr_.ghu,2),Nsam);
@ -402,7 +402,7 @@ else
for j=1:ntrans, for j=1:ntrans,
M_.params(estim_params_.param_vals(:,1)) = lpmat(istable(j),:)'; M_.params(estim_params_.param_vals(:,1)) = lpmat(istable(j),:)';
%stoch_simul([]); %stoch_simul([]);
[Tt,Rr,SteadyState,info,M_,options_,oo_] = dynare_resolve(M_,options_,oo_); [Tt,Rr,SteadyState,info,M_,options_,oo_] = dynare_resolve(M_,options_,oo_,'restrict');
% This syntax is not compatible with the current version of dynare_resolve [stepan]. % This syntax is not compatible with the current version of dynare_resolve [stepan].
%[Tt,Rr,SteadyState,info] = dynare_resolve(bayestopt_.restrict_var_list,... %[Tt,Rr,SteadyState,info] = dynare_resolve(bayestopt_.restrict_var_list,...
% bayestopt_.restrict_columns,... % bayestopt_.restrict_columns,...

1
matlab/gsa/thet2tau.m
View File

@ -17,7 +17,6 @@ M_.params(indx) = params(length(indexo)+1:end);
if ~isempty(indexo) if ~isempty(indexo)
M_.Sigma_e(indexo,indexo) = diag(params(1:length(indexo)).^2); M_.Sigma_e(indexo,indexo) = diag(params(1:length(indexo)).^2);
end end
% [A(oo_.dr.order_var,oo_.dr.order_var),B(oo_.dr.order_var,:)]=dynare_resolve;
[A,B,plouf,plouf,M_,options_,oo_] = dynare_resolve(M_,options_,oo_); [A,B,plouf,plouf,M_,options_,oo_] = dynare_resolve(M_,options_,oo_);
if flagmoments==0, if flagmoments==0,
tau = [oo_.dr.ys(oo_.dr.order_var); A(:); dyn_vech(B*M_.Sigma_e*B')]; tau = [oo_.dr.ys(oo_.dr.order_var); A(:); dyn_vech(B*M_.Sigma_e*B')];

2
matlab/initial_estimation_checks.m
View File

@ -35,7 +35,7 @@ if DynareDataset.info.nvobs>Model.exo_nbr+EstimatedParameters.nvn
end end
if DynareOptions.dsge_var if DynareOptions.dsge_var
[fval,cost_flag,info] = DsgeVarLikelihood(xparam1,gend); [fval,cost_flag,info] = DsgeVarLikelihood(xparam1,DynareDataset,DynareOptions,Model,EstimatedParameters,BayesInfo,DynareResults);
else else
[fval,cost_flag,ys,trend_coeff,info] = DsgeLikelihood(xparam1,DynareDataset,DynareOptions,Model,EstimatedParameters,BayesInfo,DynareResults); [fval,cost_flag,ys,trend_coeff,info] = DsgeLikelihood(xparam1,DynareDataset,DynareOptions,Model,EstimatedParameters,BayesInfo,DynareResults);
end end

67
matlab/mode_check.m
View File

@ -1,8 +1,8 @@
function mode_check(x,fval,hessian,gend,data,lb,ub,data_index,number_of_observations,no_more_missing_observations) function mode_check(func,x,fval,hessian,DynareDataset,DynareOptions,Model,EstimatedParameters,BayesInfo,DynareResults)
% function mode_check(x,fval,hessian,gend,data,lb,ub) % function mode_check(x,fval,hessian,gend,data,lb,ub)
% Checks the maximum likelihood mode % Checks the maximum likelihood mode
% %
% INPUTS % INPUTS
% x: mode % x: mode
% fval: value at the maximum likelihood mode % fval: value at the maximum likelihood mode
@ -14,7 +14,7 @@ function mode_check(x,fval,hessian,gend,data,lb,ub,data_index,number_of_observat
% %
% OUTPUTS % OUTPUTS
% none % none
% %
% SPECIAL REQUIREMENTS % SPECIAL REQUIREMENTS
% none % none
@ -35,18 +35,12 @@ function mode_check(x,fval,hessian,gend,data,lb,ub,data_index,number_of_observat
% You should have received a copy of the GNU General Public License % You should have received a copy of the GNU General Public License
% along with Dynare. If not, see <http://www.gnu.org/licenses/>. % along with Dynare. If not, see <http://www.gnu.org/licenses/>.
global bayestopt_ M_ options_ estim_params_ TeX = DynareOptions.TeX;
TeX = options_.TeX;
if ~isempty(hessian); if ~isempty(hessian);
[ s_min, k ] = min(diag(hessian)); [ s_min, k ] = min(diag(hessian));
end end
if options_.dsge_var
fval = DsgeVarLikelihood(x,gend); fval = feval(fun,x,DynareDataset,DynareOptions,Model,EstimatedParameters,BayesInfo,DynareResults);
else
fval = DsgeLikelihood(x,gend,data,data_index,number_of_observations,no_more_missing_observations);
end
bayestopt_.penalty=fval;
if ~isempty(hessian); if ~isempty(hessian);
disp(' ') disp(' ')
@ -55,14 +49,14 @@ if ~isempty(hessian);
disp(sprintf('Fval obtained by the minimization routine: %f', fval)) disp(sprintf('Fval obtained by the minimization routine: %f', fval))
disp(' ') disp(' ')
if s_min<eps if s_min<eps
disp(sprintf('Most negative variance %f for parameter %d (%s = %f)', s_min, k , bayestopt_.name{k}, x(k))) disp(sprintf('Most negative variance %f for parameter %d (%s = %f)', s_min, k , BayesInfo.name{k}, x(k)))
end end
end end
[nbplt,nr,nc,lr,lc,nstar] = pltorg(length(x)); [nbplt,nr,nc,lr,lc,nstar] = pltorg(length(x));
if TeX if TeX
fidTeX = fopen([M_.fname '_CheckPlots.TeX'],'w'); fidTeX = fopen([Model.fname '_CheckPlots.TeX'],'w');
fprintf(fidTeX,'%% TeX eps-loader file generated by mode_check.m (Dynare).\n'); fprintf(fidTeX,'%% TeX eps-loader file generated by mode_check.m (Dynare).\n');
fprintf(fidTeX,['%% ' datestr(now,0) '\n']); fprintf(fidTeX,['%% ' datestr(now,0) '\n']);
fprintf(fidTeX,' \n'); fprintf(fidTeX,' \n');
@ -77,7 +71,7 @@ for plt = 1:nbplt,
for k=1:min(nstar,length(x)-(plt-1)*nstar) for k=1:min(nstar,length(x)-(plt-1)*nstar)
subplot(nr,nc,k) subplot(nr,nc,k)
kk = (plt-1)*nstar+k; kk = (plt-1)*nstar+k;
[name,texname] = get_the_name(kk,TeX,M_,estim_params_,options_); [name,texname] = get_the_name(kk,TeX,Model,EstimatedParameters,DynareOptions);
if TeX if TeX
if isempty(NAMES) if isempty(NAMES)
NAMES = name; NAMES = name;
@ -88,32 +82,23 @@ for plt = 1:nbplt,
end end
end end
xx = x; xx = x;
l1 = max(lb(kk),0.5*x(kk)); l1 = max(BayesInfo.lb(kk),0.5*x(kk));
l2 = min(ub(kk),1.5*x(kk)); l2 = min(BayesInfo.ub(kk),1.5*x(kk));
z = [l1:(l2-l1)/20:l2]; z = [l1:(l2-l1)/20:l2];
if options_.mode_check_nolik==0, if DynareOptions.mode_check_nolik==0,
y = zeros(length(z),2); y = zeros(length(z),2);
dy = priordens(xx,bayestopt_.pshape,bayestopt_.p6,bayestopt_.p7,bayestopt_.p3,bayestopt_.p4); dy = priordens(xx,BayesInfo.pshape,BayesInfo.p6,BayesInfo.p7,BayesInfo.p3,BayesInfo.p4);
end end
for i=1:length(z) for i=1:length(z)
xx(kk) = z(i); xx(kk) = z(i);
if options_.dsge_var [fval, exit_flag] = feval(fun,x,DynareDataset,DynareOptions,Model,EstimatedParameters,BayesInfo,DynareResults);
[fval,cost_flag] = DsgeVarLikelihood(xx,gend); if exit_flag
if cost_flag y(i,1) = fval;
y(i,1) = fval;
else
y(i,1) = NaN;
end
else else
[fval,cost_flag] = DsgeLikelihood(xx,gend,data,data_index,number_of_observations,no_more_missing_observations); y(i,1) = NaN;
if cost_flag
y(i,1) = fval;
else
y(i,1) = NaN;
end
end end
if options_.mode_check_nolik==0 if DynareOptions.mode_check_nolik==0
lnprior = priordens(xx,bayestopt_.pshape,bayestopt_.p6,bayestopt_.p7,bayestopt_.p3,bayestopt_.p4); lnprior = priordens(xx,BayesInfo.pshape,BayesInfo.p6,BayesInfo.p7,BayesInfo.p3,BayesInfo.p4);
y(i,2) = (y(i,1)+lnprior-dy); y(i,2) = (y(i,1)+lnprior-dy);
end end
end end
@ -130,7 +115,7 @@ for plt = 1:nbplt,
axis tight axis tight
drawnow drawnow
end end
if options_.mode_check_nolik==0, if DynareOptions.mode_check_nolik==0,
if exist('OCTAVE_VERSION'), if exist('OCTAVE_VERSION'),
axes('outerposition',[0.3 0.93 0.42 0.07],'box','on'), axes('outerposition',[0.3 0.93 0.42 0.07],'box','on'),
else else
@ -142,12 +127,12 @@ for plt = 1:nbplt,
text(0.25,0.5,'log-post') text(0.25,0.5,'log-post')
text(0.69,0.5,'log-lik kernel') text(0.69,0.5,'log-lik kernel')
end end
eval(['print -depsc2 ' M_.fname '_CheckPlots' int2str(plt) '.eps']); eval(['print -depsc2 ' Model.fname '_CheckPlots' int2str(plt) '.eps']);
if ~exist('OCTAVE_VERSION') if ~exist('OCTAVE_VERSION')
eval(['print -dpdf ' M_.fname '_CheckPlots' int2str(plt)]); eval(['print -dpdf ' Model.fname '_CheckPlots' int2str(plt)]);
saveas(hh,[M_.fname '_CheckPlots' int2str(plt) '.fig']); saveas(hh,[Model.fname '_CheckPlots' int2str(plt) '.fig']);
end end
if options_.nograph, close(hh), end if DynareOptions.nograph, close(hh), end
if TeX if TeX
% TeX eps loader file % TeX eps loader file
fprintf(fidTeX,'\\begin{figure}[H]\n'); fprintf(fidTeX,'\\begin{figure}[H]\n');
@ -155,7 +140,7 @@ for plt = 1:nbplt,
fprintf(fidTeX,'\\psfrag{%s}[1][][0.5][0]{%s}\n',deblank(NAMES(jj,:)),deblank(TeXNAMES(jj,:))); fprintf(fidTeX,'\\psfrag{%s}[1][][0.5][0]{%s}\n',deblank(NAMES(jj,:)),deblank(TeXNAMES(jj,:)));
end end
fprintf(fidTeX,'\\centering \n'); fprintf(fidTeX,'\\centering \n');
fprintf(fidTeX,'\\includegraphics[scale=0.5]{%s_CheckPlots%s}\n',M_.fname,int2str(plt)); fprintf(fidTeX,'\\includegraphics[scale=0.5]{%s_CheckPlots%s}\n',Model.fname,int2str(plt));
fprintf(fidTeX,'\\caption{Check plots.}'); fprintf(fidTeX,'\\caption{Check plots.}');
fprintf(fidTeX,'\\label{Fig:CheckPlots:%s}\n',int2str(plt)); fprintf(fidTeX,'\\label{Fig:CheckPlots:%s}\n',int2str(plt));
fprintf(fidTeX,'\\end{figure}\n'); fprintf(fidTeX,'\\end{figure}\n');

1
matlab/thet2tau.m
View File

@ -36,7 +36,6 @@ M_.params(indx) = params(length(indexo)+1:end);
if ~isempty(indexo) if ~isempty(indexo)
M_.Sigma_e(indexo,indexo) = diag(params(1:length(indexo)).^2); M_.Sigma_e(indexo,indexo) = diag(params(1:length(indexo)).^2);
end end
% [A(oo_.dr.order_var,oo_.dr.order_var),B(oo_.dr.order_var,:)]=dynare_resolve;
[A,B,tele,tubbies,M_,options_,oo_] = dynare_resolve(M_,options_,oo_); [A,B,tele,tubbies,M_,options_,oo_] = dynare_resolve(M_,options_,oo_);
if flagmoments==0, if flagmoments==0,
tau = [oo_.dr.ys(oo_.dr.order_var); A(:); dyn_vech(B*M_.Sigma_e*B')]; tau = [oo_.dr.ys(oo_.dr.order_var); A(:); dyn_vech(B*M_.Sigma_e*B')];