Skip to content
GitLab
Projects
Groups
Snippets
/
Help
Help
Support
Community forum
Keyboard shortcuts
?
Submit feedback
Contribute to GitLab
Sign in / Register
Toggle navigation
Menu
Open sidebar
MichelJuillard
dynare
Commits
222d8d18
Commit
222d8d18
authored
Oct 05, 2015
by
MichelJuillard
Browse files
More fixing related to objective_function_penalty_base
parent
dc3c01f2
Changes
15
Expand all
Hide whitespace changes
Inline
Side-by-side
matlab/TaRB_optimizer_wrapper.m
View file @
222d8d18
function
[
fval
,
DLIK
,
Hess
,
exit_flag
]
=
TaRB_optimizer_wrapper
(
optpar
,
par_vector
,
parameterindices
,
TargetFun
,
varargin
)
function
[
fval
,
info
,
exit_flag
,
DLIK
,
Hess
,
SteadyState
,
trend_coeff
]
=
TaRB_optimizer_wrapper
(
optpar
,
par_vector
,
parameterindices
,
TargetFun
,
varargin
)
% function [fval,DLIK,Hess,exit_flag] = TaRB_optimizer_wrapper(optpar,par_vector,parameterindices,TargetFun,varargin)
% Wrapper function for target function used in TaRB algorithm; reassembles
% full parameter vector before calling target function
...
...
@@ -36,5 +36,5 @@ function [fval,DLIK,Hess,exit_flag] = TaRB_optimizer_wrapper(optpar,par_vector,p
% along with Dynare. If not, see <http://www.gnu.org/licenses/>.
par_vector
(
parameterindices
,:)
=
optpar
;
%reassemble parameter
[
fval
,
DLIK
,
Hess
,
exit_flag
]
=
feval
(
TargetFun
,
par_vector
,
varargin
{:});
%call target function
[
fval
,
info
,
exit_flag
,
DLIK
,
Hess
,
SteadyState
,
trend_coeff
]
=
feval
(
TargetFun
,
par_vector
,
varargin
{:});
%call target function
matlab/dsge_likelihood.m
View file @
222d8d18
This diff is collapsed.
Click to expand it.
matlab/dsge_likelihood_1.m
0 → 100644
View file @
222d8d18
This diff is collapsed.
Click to expand it.
matlab/dsge_var_likelihood.m
View file @
222d8d18
function
[
fval
,
grad
,
hess
,
exit_flag
,
SteadyState
,
trend_coeff
,
info
,
PHI
,
SIGMAu
,
iXX
,
prior
]
=
dsge_var_likelihood
(
xparam1
,
DynareDataset
,
DynareInfo
,
DynareOptions
,
Model
,
EstimatedParameters
,
BayesInfo
,
BoundsInfo
,
DynareResults
)
% Evaluates the posterior kernel of the bvar-dsge model.
function
[
fval
,
grad
,
hess
,
exit_flag
,
info
,
PHI
,
SIGMAu
,
iXX
,
prior
]
=
dsge_var_likelihood
(
xparam1
,
DynareDataset
,
DynareOptions
,
Model
,
EstimatedParameters
,
BayesInfo
,
DynareResults
)
% Evaluates the posterior kernel of the bvar-dsge model.
Deprecated interface.
%
% INPUTS
% o xparam1 [double] Vector of model's parameters.
...
...
@@ -8,10 +8,6 @@ function [fval,grad,hess,exit_flag,SteadyState,trend_coeff,info,PHI,SIGMAu,iXX,p
% OUTPUTS
% o fval [double] Value of the posterior kernel at xparam1.
% o cost_flag [integer] Zero if the function returns a penalty, one otherwise.
% o SteadyState [double] Steady state vector possibly recomputed
% by call to dynare_results()
% o trend_coeff [double] place holder for trend coefficients,
% currently not supported by dsge_var
% o info [integer] Vector of informations about the penalty.
% 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).
...
...
@@ -38,275 +34,6 @@ function [fval,grad,hess,exit_flag,SteadyState,trend_coeff,info,PHI,SIGMAu,iXX,p
% You should have received a copy of the GNU General Public License
% along with Dynare. If not, see <http://www.gnu.org/licenses/>.
persistent
dsge_prior_weight_idx
grad
=
[];
hess
=
[];
exit_flag
=
[];
info
=
0
;
PHI
=
[];
SIGMAu
=
[];
iXX
=
[];
prior
=
[];
SteadyState
=
[];
trend_coeff
=
[];
% Initialization of of the index for parameter dsge_prior_weight in Model.params.
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.
nx
=
EstimatedParameters
.
nvx
+
EstimatedParameters
.
np
;
% Get the number of observed variables in the VAR model.
NumberOfObservedVariables
=
DynareDataset
.
vobs
;
% Get the number of observations.
NumberOfObservations
=
DynareDataset
.
nobs
;
% 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
;
if
~
DynareOptions
.
noconstant
NumberOfParameters
=
NumberOfParameters
+
1
;
end
% Get empirical second order moments for the observed variables.
mYY
=
evalin
(
'base'
,
'mYY'
);
mYX
=
evalin
(
'base'
,
'mYX'
);
mXY
=
evalin
(
'base'
,
'mXY'
);
mXX
=
evalin
(
'base'
,
'mXX'
);
% Initialize some of the output arguments.
fval
=
[];
exit_flag
=
1
;
% Return, with endogenous penalty, if some dsge-parameters are smaller than the lower bound of the prior domain.
if
DynareOptions
.
mode_compute
~=
1
&&
any
(
xparam1
<
BoundsInfo
.
lb
)
fval
=
Inf
;
exit_flag
=
0
;
info
(
1
)
=
41
;
k
=
find
(
xparam1
<
BoundsInfo
.
lb
);
info
(
2
)
=
sum
((
BoundsInfo
.
lb
(
k
)
-
xparam1
(
k
))
.^
2
);
return
;
end
% Return, with endogenous penalty, if some dsge-parameters are greater than the upper bound of the prior domain.
if
DynareOptions
.
mode_compute
~=
1
&&
any
(
xparam1
>
BoundsInfo
.
ub
)
fval
=
Inf
;
exit_flag
=
0
;
info
(
1
)
=
42
;
k
=
find
(
xparam1
>
BoundsInfo
.
ub
);
info
(
2
)
=
sum
((
xparam1
(
k
)
-
BoundsInfo
.
ub
(
k
))
.^
2
);
return
;
end
% Get the variance of each structural innovation.
Q
=
Model
.
Sigma_e
;
for
i
=
1
:
EstimatedParameters
.
nvx
k
=
EstimatedParameters
.
var_exo
(
i
,
1
);
Q
(
k
,
k
)
=
xparam1
(
i
)
*
xparam1
(
i
);
end
offset
=
EstimatedParameters
.
nvx
;
% Update Model.params and Model.Sigma_e.
Model
.
params
(
EstimatedParameters
.
param_vals
(:,
1
))
=
xparam1
(
offset
+
1
:
end
);
Model
.
Sigma_e
=
Q
;
% Get the weight of the dsge prior.
dsge_prior_weight
=
Model
.
params
(
dsge_prior_weight_idx
);
% Is the dsge prior proper?
if
dsge_prior_weight
<
(
NumberOfParameters
+
NumberOfObservedVariables
)/
...
NumberOfObservations
;
fval
=
Inf
;
exit_flag
=
0
;
info
(
1
)
=
51
;
info
(
2
)
=
abs
(
NumberOfObservations
*
dsge_prior_weight
-
(
NumberOfParameters
+
NumberOfObservedVariables
));
% info(2)=dsge_prior_weight;
% info(3)=(NumberOfParameters+NumberOfObservedVariables)/DynareDataset.nobs;
return
end
%------------------------------------------------------------------------------
% 2. call model setup & reduction program
%------------------------------------------------------------------------------
% 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
||
info
(
1
)
==
7
||
info
(
1
)
==
8
||
...
info
(
1
)
==
22
||
info
(
1
)
==
24
||
info
(
1
)
==
25
||
info
(
1
)
==
10
fval
=
Inf
;
info
(
2
)
=
0.1
;
exit_flag
=
0
;
return
elseif
info
(
1
)
==
3
||
info
(
1
)
==
4
||
info
(
1
)
==
19
||
info
(
1
)
==
20
||
info
(
1
)
==
21
fval
=
Inf
;
info
(
2
)
=
0.1
;
exit_flag
=
0
;
return
end
% Define the mean/steady state vector.
if
~
DynareOptions
.
noconstant
if
DynareOptions
.
loglinear
constant
=
transpose
(
log
(
SteadyState
(
BayesInfo
.
mfys
)));
else
constant
=
transpose
(
SteadyState
(
BayesInfo
.
mfys
));
end
else
constant
=
zeros
(
1
,
NumberOfObservedVariables
);
end
%------------------------------------------------------------------------------
% 3. theoretical moments (second order)
%------------------------------------------------------------------------------
% Compute the theoretical second order moments
tmp0
=
lyapunov_symm
(
T
,
R
*
Q
*
R
'
,
DynareOptions
.
lyapunov_fixed_point_tol
,
DynareOptions
.
qz_criterium
,
DynareOptions
.
lyapunov_complex_threshold
,
[],
[],
DynareOptions
.
debug
);
mf
=
BayesInfo
.
mf1
;
% Get the non centered second order moments
TheoreticalAutoCovarianceOfTheObservedVariables
=
zeros
(
NumberOfObservedVariables
,
NumberOfObservedVariables
,
NumberOfLags
+
1
);
TheoreticalAutoCovarianceOfTheObservedVariables
(:,:,
1
)
=
tmp0
(
mf
,
mf
)
+
constant
'*
constant
;
for
lag
=
1
:
NumberOfLags
tmp0
=
T
*
tmp0
;
TheoreticalAutoCovarianceOfTheObservedVariables
(:,:,
lag
+
1
)
=
tmp0
(
mf
,
mf
)
+
constant
'*
constant
;
end
% Build the theoretical "covariance" between Y and X
GYX
=
zeros
(
NumberOfObservedVariables
,
NumberOfParameters
);
for
i
=
1
:
NumberOfLags
GYX
(:,(
i
-
1
)
*
NumberOfObservedVariables
+
1
:
i
*
NumberOfObservedVariables
)
=
TheoreticalAutoCovarianceOfTheObservedVariables
(:,:,
i
+
1
);
end
if
~
DynareOptions
.
noconstant
GYX
(:,
end
)
=
constant
'
;
end
% Build the theoretical "covariance" between X and X
GXX
=
kron
(
eye
(
NumberOfLags
),
TheoreticalAutoCovarianceOfTheObservedVariables
(:,:,
1
));
for
i
=
1
:
NumberOfLags
-
1
tmp1
=
diag
(
ones
(
NumberOfLags
-
i
,
1
),
i
);
tmp2
=
diag
(
ones
(
NumberOfLags
-
i
,
1
),
-
i
);
GXX
=
GXX
+
kron
(
tmp1
,
TheoreticalAutoCovarianceOfTheObservedVariables
(:,:,
i
+
1
));
GXX
=
GXX
+
kron
(
tmp2
,
TheoreticalAutoCovarianceOfTheObservedVariables
(:,:,
i
+
1
)
'
);
end
if
~
DynareOptions
.
noconstant
% Add one row and one column to GXX
GXX
=
[
GXX
,
kron
(
ones
(
NumberOfLags
,
1
),
constant
'
)
;
[
kron
(
ones
(
1
,
NumberOfLags
),
constant
)
,
1
]
];
end
GYY
=
TheoreticalAutoCovarianceOfTheObservedVariables
(:,:,
1
);
assignin
(
'base'
,
'GYY'
,
GYY
);
assignin
(
'base'
,
'GXX'
,
GXX
);
assignin
(
'base'
,
'GYX'
,
GYX
);
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
*
NumberOfObservations
*
TheoreticalAutoCovarianceOfTheObservedVariables
(:,:,
1
)
+
mYY
;
tmp1
=
dsge_prior_weight
*
NumberOfObservations
*
GYX
+
mYX
;
tmp2
=
inv
(
dsge_prior_weight
*
NumberOfObservations
*
GXX
+
mXX
);
SIGMAu
=
tmp0
-
tmp1
*
tmp2
*
tmp1
'; clear('
tmp0
'
);
[
SIGMAu_is_positive_definite
,
penalty
]
=
ispd
(
SIGMAu
);
if
~
SIGMAu_is_positive_definite
fval
=
Inf
;
info
(
1
)
=
52
;
info
(
2
)
=
penalty
;
exit_flag
=
0
;
return
;
end
SIGMAu
=
SIGMAu
/
(
NumberOfObservations
*
(
1
+
dsge_prior_weight
));
PHI
=
tmp2
*
tmp1
'; clear('
tmp1
'
);
prodlng1
=
sum
(
gammaln
(
.
5
*
((
1
+
dsge_prior_weight
)
*
NumberOfObservations
-
...
NumberOfObservedVariables
*
NumberOfLags
...
+
1
-
(
1
:
NumberOfObservedVariables
)
'
)));
prodlng2
=
sum
(
gammaln
(
.
5
*
(
dsge_prior_weight
*
NumberOfObservations
-
...
NumberOfObservedVariables
*
NumberOfLags
...
+
1
-
(
1
:
NumberOfObservedVariables
)
'
)));
lik
=
.
5
*
NumberOfObservedVariables
*
log
(
det
(
dsge_prior_weight
*
NumberOfObservations
*
GXX
+
mXX
))
...
+
.
5
*
((
dsge_prior_weight
+
1
)
*
NumberOfObservations
-
NumberOfParameters
)
*
log
(
det
((
dsge_prior_weight
+
1
)
*
NumberOfObservations
*
SIGMAu
))
...
-
.
5
*
NumberOfObservedVariables
*
log
(
det
(
dsge_prior_weight
*
NumberOfObservations
*
GXX
))
...
-
.
5
*
(
dsge_prior_weight
*
NumberOfObservations
-
NumberOfParameters
)
*
log
(
det
(
dsge_prior_weight
*
NumberOfObservations
*
(
GYY
-
GYX
*
inv
(
GXX
)
*
GYX
'
)))
...
+
.
5
*
NumberOfObservedVariables
*
NumberOfObservations
*
log
(
2
*
pi
)
...
-
.
5
*
log
(
2
)
*
NumberOfObservedVariables
*
((
dsge_prior_weight
+
1
)
*
NumberOfObservations
-
NumberOfParameters
)
...
+
.
5
*
log
(
2
)
*
NumberOfObservedVariables
*
(
dsge_prior_weight
*
NumberOfObservations
-
NumberOfParameters
)
...
-
prodlng1
+
prodlng2
;
else
% Evaluation of the likelihood of the dsge-var model when the dsge prior weight is infinite.
iGXX
=
inv
(
GXX
);
SIGMAu
=
GYY
-
GYX
*
iGXX
*
transpose
(
GYX
);
PHI
=
iGXX
*
transpose
(
GYX
);
lik
=
NumberOfObservations
*
(
log
(
det
(
SIGMAu
))
+
NumberOfObservedVariables
*
log
(
2
*
pi
)
+
...
trace
(
inv
(
SIGMAu
)
*
(
mYY
-
transpose
(
mYX
*
PHI
)
-
mYX
*
PHI
+
transpose
(
PHI
)
*
mXX
*
PHI
)/
NumberOfObservations
));
lik
=
.
5
*
lik
;
% Minus likelihood
end
if
isnan
(
lik
)
info
(
1
)
=
45
;
info
(
2
)
=
0.1
;
fval
=
Inf
;
exit_flag
=
0
;
return
end
if
imag
(
lik
)
~=
0
info
(
1
)
=
46
;
info
(
2
)
=
0.1
;
fval
=
Inf
;
exit_flag
=
0
;
return
end
% 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
);
if
isnan
(
fval
)
info
(
1
)
=
47
;
info
(
2
)
=
0.1
;
fval
=
Inf
;
exit_flag
=
0
;
return
end
if
imag
(
fval
)
~=
0
info
(
1
)
=
48
;
info
(
2
)
=
0.1
;
fval
=
Inf
;
exit_flag
=
0
;
return
end
if
(
nargout
==
10
)
if
isinf
(
dsge_prior_weight
)
iXX
=
iGXX
;
else
iXX
=
tmp2
;
end
end
if
(
nargout
==
11
)
if
isinf
(
dsge_prior_weight
)
iXX
=
iGXX
;
else
iXX
=
tmp2
;
end
iGXX
=
inv
(
GXX
);
prior
.
SIGMAstar
=
GYY
-
GYX
*
iGXX
*
GYX
'
;
prior
.
PHIstar
=
iGXX
*
transpose
(
GYX
);
prior
.
ArtificialSampleSize
=
fix
(
dsge_prior_weight
*
NumberOfObservations
);
prior
.
DF
=
prior
.
ArtificialSampleSize
-
NumberOfParameters
-
NumberOfObservedVariables
;
prior
.
iGXX
=
iGXX
;
end
if
fval
==
Inf
pause
end
\ No newline at end of file
[
fval
,
info
,
exit_flag
,
grad
,
hess
,
SteadyState
,
trend_coeff
,
PHI
,
SIGMAu
,
iXX
,
prior
]
=
...
dsge_var_likelihood_1
(
xparam1
,
DynareDataset
,
DynareInfo
,
DynareOptions
,
Model
,
...
EstimatedParameters
,
BayesInfo
,
BoundsInfo
,
DynareResults
);
\ No newline at end of file
matlab/dsge_var_likelihood_1.m
0 → 100644
View file @
222d8d18
function
[
fval
,
info
,
exit_flag
,
grad
,
hess
,
SteadyState
,
trend_coeff
,
PHI
,
SIGMAu
,
iXX
,
prior
]
=
dsge_var_likelihood_1
(
xparam1
,
DynareDataset
,
DynareInfo
,
DynareOptions
,
Model
,
EstimatedParameters
,
BayesInfo
,
BoundsInfo
,
DynareResults
)
% Evaluates the posterior kernel of the bvar-dsge model.
%
% INPUTS
% o xparam1 [double] Vector of model's parameters.
% o gend [integer] Number of observations (without conditionning observations for the lags).
%
% OUTPUTS
% o fval [double] Value of the posterior kernel at xparam1.
% o cost_flag [integer] Zero if the function returns a penalty, one otherwise.
% o SteadyState [double] Steady state vector possibly recomputed
% by call to dynare_results()
% o trend_coeff [double] place holder for trend coefficients,
% currently not supported by dsge_var
% o info [integer] Vector of informations about the penalty.
% 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 iXX [double] inv(X'X).
% o prior [double] a matlab structure describing the dsge-var prior.
%
% SPECIAL REQUIREMENTS
% None.
% Copyright (C) 2006-2012 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/>.
persistent
dsge_prior_weight_idx
grad
=
[];
hess
=
[];
exit_flag
=
[];
info
=
0
;
PHI
=
[];
SIGMAu
=
[];
iXX
=
[];
prior
=
[];
SteadyState
=
[];
trend_coeff
=
[];
% Initialization of of the index for parameter dsge_prior_weight in Model.params.
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.
nx
=
EstimatedParameters
.
nvx
+
EstimatedParameters
.
np
;
% Get the number of observed variables in the VAR model.
NumberOfObservedVariables
=
DynareDataset
.
vobs
;
% Get the number of observations.
NumberOfObservations
=
DynareDataset
.
nobs
;
% 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
;
if
~
DynareOptions
.
noconstant
NumberOfParameters
=
NumberOfParameters
+
1
;
end
% Get empirical second order moments for the observed variables.
mYY
=
evalin
(
'base'
,
'mYY'
);
mYX
=
evalin
(
'base'
,
'mYX'
);
mXY
=
evalin
(
'base'
,
'mXY'
);
mXX
=
evalin
(
'base'
,
'mXX'
);
% Initialize some of the output arguments.
fval
=
[];
exit_flag
=
1
;
% Return, with endogenous penalty, if some dsge-parameters are smaller than the lower bound of the prior domain.
if
DynareOptions
.
mode_compute
~=
1
&&
any
(
xparam1
<
BoundsInfo
.
lb
)
fval
=
Inf
;
exit_flag
=
0
;
info
(
1
)
=
41
;
k
=
find
(
xparam1
<
BoundsInfo
.
lb
);
info
(
2
)
=
sum
((
BoundsInfo
.
lb
(
k
)
-
xparam1
(
k
))
.^
2
);
return
;
end
% Return, with endogenous penalty, if some dsge-parameters are greater than the upper bound of the prior domain.
if
DynareOptions
.
mode_compute
~=
1
&&
any
(
xparam1
>
BoundsInfo
.
ub
)
fval
=
Inf
;
exit_flag
=
0
;
info
(
1
)
=
42
;
k
=
find
(
xparam1
>
BoundsInfo
.
ub
);
info
(
2
)
=
sum
((
xparam1
(
k
)
-
BoundsInfo
.
ub
(
k
))
.^
2
);
return
;
end
% Get the variance of each structural innovation.
Q
=
Model
.
Sigma_e
;
for
i
=
1
:
EstimatedParameters
.
nvx
k
=
EstimatedParameters
.
var_exo
(
i
,
1
);
Q
(
k
,
k
)
=
xparam1
(
i
)
*
xparam1
(
i
);
end
offset
=
EstimatedParameters
.
nvx
;
% Update Model.params and Model.Sigma_e.
Model
.
params
(
EstimatedParameters
.
param_vals
(:,
1
))
=
xparam1
(
offset
+
1
:
end
);
Model
.
Sigma_e
=
Q
;
% Get the weight of the dsge prior.
dsge_prior_weight
=
Model
.
params
(
dsge_prior_weight_idx
);
% Is the dsge prior proper?
if
dsge_prior_weight
<
(
NumberOfParameters
+
NumberOfObservedVariables
)/
...
NumberOfObservations
;
fval
=
Inf
;
exit_flag
=
0
;
info
(
1
)
=
51
;
info
(
2
)
=
abs
(
NumberOfObservations
*
dsge_prior_weight
-
(
NumberOfParameters
+
NumberOfObservedVariables
));
% info(2)=dsge_prior_weight;
% info(3)=(NumberOfParameters+NumberOfObservedVariables)/DynareDataset.nobs;
return
end
%------------------------------------------------------------------------------
% 2. call model setup & reduction program
%------------------------------------------------------------------------------
% 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
||
info
(
1
)
==
7
||
info
(
1
)
==
8
||
...
info
(
1
)
==
22
||
info
(
1
)
==
24
||
info
(
1
)
==
25
||
info
(
1
)
==
10
fval
=
Inf
;
info
(
2
)
=
0.1
;
exit_flag
=
0
;
return
elseif
info
(
1
)
==
3
||
info
(
1
)
==
4
||
info
(
1
)
==
19
||
info
(
1
)
==
20
||
info
(
1
)
==
21
fval
=
Inf
;
info
(
2
)
=
0.1
;
exit_flag
=
0
;
return
end
% Define the mean/steady state vector.
if
~
DynareOptions
.
noconstant
if
DynareOptions
.
loglinear
constant
=
transpose
(
log
(
SteadyState
(
BayesInfo
.
mfys
)));
else
constant
=
transpose
(
SteadyState
(
BayesInfo
.
mfys
));
end
else
constant
=
zeros
(
1
,
NumberOfObservedVariables
);
end
%------------------------------------------------------------------------------
% 3. theoretical moments (second order)
%------------------------------------------------------------------------------
% Compute the theoretical second order moments
tmp0
=
lyapunov_symm
(
T
,
R
*
Q
*
R
'
,
DynareOptions
.
lyapunov_fixed_point_tol
,
DynareOptions
.
qz_criterium
,
DynareOptions
.
lyapunov_complex_threshold
,
[],
[],
DynareOptions
.
debug
);
mf
=
BayesInfo
.
mf1
;
% Get the non centered second order moments
TheoreticalAutoCovarianceOfTheObservedVariables
=
zeros
(
NumberOfObservedVariables
,
NumberOfObservedVariables
,
NumberOfLags
+
1
);
TheoreticalAutoCovarianceOfTheObservedVariables
(:,:,
1
)
=
tmp0
(
mf
,
mf
)
+
constant
'*
constant
;
for
lag
=
1
:
NumberOfLags
tmp0
=
T
*
tmp0
;
TheoreticalAutoCovarianceOfTheObservedVariables
(:,:,
lag
+
1
)
=
tmp0
(
mf
,
mf
)
+
constant
'*
constant
;
end
% Build the theoretical "covariance" between Y and X
GYX
=
zeros
(
NumberOfObservedVariables
,
NumberOfParameters
);
for
i
=
1
:
NumberOfLags
GYX
(:,(
i
-
1
)
*
NumberOfObservedVariables
+
1
:
i
*
NumberOfObservedVariables
)
=
TheoreticalAutoCovarianceOfTheObservedVariables
(:,:,
i
+
1
);
end
if
~
DynareOptions
.
noconstant
GYX
(:,
end
)
=
constant
'
;
end
% Build the theoretical "covariance" between X and X
GXX
=
kron
(
eye
(
NumberOfLags
),
TheoreticalAutoCovarianceOfTheObservedVariables
(:,:,
1
));
for
i
=
1
:
NumberOfLags
-
1
tmp1
=
diag
(
ones
(
NumberOfLags
-
i
,
1
),
i
);
tmp2
=
diag
(
ones
(
NumberOfLags
-
i
,
1
),
-
i
);
GXX
=
GXX
+
kron
(
tmp1
,
TheoreticalAutoCovarianceOfTheObservedVariables
(:,:,
i
+
1
));
GXX
=
GXX
+
kron
(
tmp2
,
TheoreticalAutoCovarianceOfTheObservedVariables
(:,:,
i
+
1
)
'
);
end
if
~
DynareOptions
.
noconstant
% Add one row and one column to GXX
GXX
=
[
GXX
,
kron
(
ones
(
NumberOfLags
,
1
),
constant
'
)
;
[
kron
(
ones
(
1
,
NumberOfLags
),
constant
)
,
1
]
];
end
GYY
=
TheoreticalAutoCovarianceOfTheObservedVariables
(:,:,
1
);
assignin
(
'base'
,
'GYY'
,
GYY
);
assignin
(
'base'
,
'GXX'
,
GXX
);
assignin
(
'base'
,
'GYX'
,
GYX
);
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
*
NumberOfObservations
*
TheoreticalAutoCovarianceOfTheObservedVariables
(:,:,
1
)
+
mYY
;
tmp1
=
dsge_prior_weight
*
NumberOfObservations
*
GYX
+
mYX
;
tmp2
=
inv
(
dsge_prior_weight
*
NumberOfObservations
*
GXX
+
mXX
);
SIGMAu
=
tmp0
-
tmp1
*
tmp2
*
tmp1
'; clear('
tmp0
'
);
[
SIGMAu_is_positive_definite
,
penalty
]
=
ispd
(
SIGMAu
);
if
~
SIGMAu_is_positive_definite
fval
=
Inf
;
info
(
1
)
=
52
;
info
(
2
)
=
penalty
;
exit_flag
=
0
;
return
;
end
SIGMAu
=
SIGMAu
/
(
NumberOfObservations
*
(
1
+
dsge_prior_weight
));
PHI
=
tmp2
*
tmp1
'; clear('
tmp1
'
);
prodlng1
=
sum
(
gammaln
(
.
5
*
((
1
+
dsge_prior_weight
)
*
NumberOfObservations
-
...
NumberOfObservedVariables
*
NumberOfLags
...
+
1
-
(
1
:
NumberOfObservedVariables
)
'
)));
prodlng2
=
sum
(
gammaln
(
.
5
*
(
dsge_prior_weight
*
NumberOfObservations
-
...
NumberOfObservedVariables
*
NumberOfLags
...
+
1
-
(
1
:
NumberOfObservedVariables
)
'
)));
lik
=
.
5
*
NumberOfObservedVariables
*
log
(
det
(
dsge_prior_weight
*
NumberOfObservations
*
GXX
+
mXX
))
...
+
.
5
*
((
dsge_prior_weight
+
1
)
*
NumberOfObservations
-
NumberOfParameters
)
*
log
(
det
((
dsge_prior_weight
+
1
)
*
NumberOfObservations
*
SIGMAu
))
...
-
.
5
*
NumberOfObservedVariables
*
log
(
det
(
dsge_prior_weight
*
NumberOfObservations
*
GXX
))
...
-
.
5
*
(
dsge_prior_weight
*
NumberOfObservations
-
NumberOfParameters
)
*
log
(
det
(
dsge_prior_weight
*
NumberOfObservations
*
(
GYY
-
GYX
*
inv
(
GXX
)
*
GYX
'
)))
...
+
.
5
*
NumberOfObservedVariables
*
NumberOfObservations
*
log
(
2
*
pi
)
...
-
.
5
*
log
(
2
)
*
NumberOfObservedVariables
*
((
dsge_prior_weight
+
1
)
*
NumberOfObservations
-
NumberOfParameters
)
...
+
.
5
*
log
(
2
)
*
NumberOfObservedVariables
*
(
dsge_prior_weight
*
NumberOfObservations
-
NumberOfParameters
)
...
-
prodlng1
+
prodlng2
;
else
% Evaluation of the likelihood of the dsge-var model when the dsge prior weight is infinite.
iGXX
=
inv
(
GXX
);
SIGMAu
=
GYY
-
GYX
*
iGXX
*
transpose
(
GYX
);
PHI
=
iGXX
*
transpose
(
GYX
);
lik
=
NumberOfObservations
*
(
log
(
det
(
SIGMAu
))
+
NumberOfObservedVariables
*
log
(
2
*
pi
)
+
...
trace
(
inv
(
SIGMAu
)
*
(
mYY
-
transpose
(
mYX
*
PHI
)
-
mYX
*
PHI
+
transpose
(
PHI
)
*
mXX
*
PHI
)/
NumberOfObservations
));
lik
=
.
5
*
lik
;
% Minus likelihood
end
if
isnan
(
lik
)
info
(
1
)
=
45
;
info
(
2
)
=
0.1
;
fval
=
Inf
;
exit_flag
=
0
;
return
end
if
imag
(
lik
)
~=
0
info
(
1
)
=
46
;
info
(
2
)
=
0.1
;
fval
=
Inf
;
exit_flag
=
0
;
return
end
% 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
);
if
isnan
(
fval
)
info
(
1
)
=
47
;
info
(
2
)
=
0.1
;
fval
=
Inf
;
exit_flag
=
0
;