2019-03-26 15:36:55 +01:00
function ds = olsgibbs ( ds, eqtag, BetaPriorExpectation, BetaPriorVariance, s2, nu, ndraws, discarddraws, thin, fitted_names_dict, model_name, param_names, ds_range)
%function ds = olsgibbs(ds, eqtag, BetaPriorExpectation, BetaPriorVariance, s2, nu, ndraws, discarddraws, thin, fitted_names_dict, model_name, param_names, ds_range)
2019-01-11 11:34:48 +01:00
% Implements Gibbs Sampling for univariate linear model.
2018-09-28 16:09:50 +02:00
%
% INPUTS
2018-10-06 16:08:28 +02:00
% - ds [dseries] dataset.
% - eqtag [string] name of equation tag to estimate.
% - BetaPriorExpectation [double] vector with n elements, prior expectation of β.
% - BetaPriorVariance [double] n*n matrix, prior variance of β.
% - s2 [double] scalar, first hyperparameter for h.
% - nu [integer] scalar, second hyperparameter for h.
% - ndraws [integer] scalar, total number of draws (Gibbs sampling)
% - discarddraws [integer] scalar, number of draws to be discarded.
% - thin [integer] scalar, if thin == N, save every Nth draw (default is 1).
2018-10-24 17:13:06 +02:00
% - fitted_names_dict [cell] Nx2 or Nx3 cell array to be used in naming fitted
% values; first column is the equation tag,
% second column is the name of the
% associated fitted value, third column
% (if it exists) is the function name of
% the transformation to perform on the
% fitted value.
2019-02-28 11:49:50 +01:00
% - model_name [string] name to use in oo_ and inc file
2019-03-07 19:27:26 +01:00
% - param_names [cellstr] list of parameters to estimate (if
% empty, estimate all)
2019-03-26 15:36:55 +01:00
% - ds_range [dates] range of dates to use in estimation
2018-09-28 16:09:50 +02:00
%
% OUTPUTS
2018-10-24 17:13:06 +02:00
% - ds [dseries] dataset updated with fitted value
2018-09-28 16:09:50 +02:00
%
% SPECIAL REQUIREMENTS
2019-01-24 12:42:08 +01:00
% dynare must have been run with the option: json=compute
2018-09-28 16:09:50 +02:00
2021-02-04 16:54:09 +01:00
% Copyright (C) 2018-2021 Dynare Team
2018-09-28 16:09:50 +02:00
%
% 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/>.
global M_ oo_ options_
2019-01-11 11:34:48 +01:00
%% Check input
2019-03-26 15:36:55 +01:00
if nargin < 7 || nargin > 13
error ( ' Incorrect number of arguments' )
2018-09-28 16:09:50 +02:00
end
2018-10-06 16:08:28 +02:00
2019-01-11 11:34:48 +01:00
if isempty ( ds ) || ~ isdseries ( ds )
2018-09-28 16:09:50 +02:00
error ( ' The 1st argument must be a dseries' )
end
2018-10-06 16:08:28 +02:00
2018-09-28 16:09:50 +02:00
if ~ ischar ( eqtag )
error ( ' The 2nd argument must be a string' )
end
2018-10-06 16:08:28 +02:00
if ~ isvector ( BetaPriorExpectation )
2018-09-28 16:09:50 +02:00
error ( ' The 3rd argument must be a vector' )
else
2018-10-06 16:08:28 +02:00
if ~ isempty ( BetaPriorExpectation )
BetaPriorExpectation = transpose ( BetaPriorExpectation ( : ) ) ;
2018-09-28 16:09:50 +02:00
end
end
2018-10-06 16:08:28 +02:00
if ~ ismatrix ( BetaPriorVariance ) || length ( BetaPriorExpectation ) ~= length ( BetaPriorVariance )
error ( ' The 4th argument (BetaPriorVariance) must be a square matrix with the same dimension as the third argument (BetaPriorExpectation)' )
2018-09-28 16:09:50 +02:00
else
2018-10-06 16:08:28 +02:00
warning ( ' off' , ' MATLAB:singularMatrix' )
BetaInversePriorVariance = eye ( length ( BetaPriorVariance ) ) / BetaPriorVariance ;
warning ( ' on' , ' MATLAB:singularMatrix' )
2018-09-28 16:09:50 +02:00
end
2018-10-06 16:08:28 +02:00
2018-09-28 16:09:50 +02:00
if ~ isreal ( s2 )
2018-10-06 16:08:28 +02:00
error ( ' The 5th argument (s2) must be a double' )
2018-09-28 16:09:50 +02:00
end
2018-10-06 16:08:28 +02:00
2018-09-28 16:09:50 +02:00
if ~ isint ( nu )
2018-10-06 16:08:28 +02:00
error ( ' The 6th argument (nu) must be an integer' )
2018-09-28 16:09:50 +02:00
end
2018-10-06 16:08:28 +02:00
2018-09-28 16:09:50 +02:00
if ~ isint ( ndraws )
2018-10-06 16:08:28 +02:00
error ( ' The 7th argument (ndraws) must be an integer' )
2018-09-28 16:09:50 +02:00
end
2018-10-06 16:08:28 +02:00
2019-01-22 00:05:08 +01:00
if nargin < = 7
2018-09-28 16:09:50 +02:00
discarddraws = 0 ;
else
if ~ isint ( discarddraws )
2018-10-06 16:08:28 +02:00
error ( ' The 8th argument (discardeddraws), if provided, must be an integer' )
else
2019-01-22 00:05:08 +01:00
if discarddraws > = ndraws
2018-10-06 16:08:28 +02:00
error ( ' The 8th argument (discardeddraws) must be smaller than the 7th argument (ndraws)' )
end
2018-09-28 16:09:50 +02:00
end
end
2018-10-06 16:08:28 +02:00
2019-01-22 00:05:08 +01:00
if nargin < = 8
2018-09-28 16:09:50 +02:00
thin = 1 ;
else
if ~ isint ( thin )
2018-10-06 16:08:28 +02:00
error ( ' The 9th argument, must be an integer' )
2018-09-28 16:09:50 +02:00
end
end
2019-01-22 00:05:08 +01:00
if nargin < = 9
2018-10-24 17:13:06 +02:00
fitted_names_dict = { } ;
else
2019-01-22 00:05:08 +01:00
if ~ isempty ( fitted_names_dict ) && ...
2018-10-24 17:13:06 +02:00
( ~ iscell ( fitted_names_dict ) || ...
2019-01-22 00:05:08 +01:00
( size ( fitted_names_dict , 2 ) < 2 || size ( fitted_names_dict , 2 ) > 3 ) )
2018-10-24 17:13:06 +02:00
error ( ' The 10th argument must be an Nx2 or Nx3 cell array' ) ;
end
end
2019-02-28 11:49:50 +01:00
if nargin < = 10
model_name = eqtag ;
else
if ~ isvarname ( model_name )
error ( ' The 11th argument must be a valid string' ) ;
end
end
2019-03-07 19:27:26 +01:00
if nargin < = 11
param_names = { } ;
else
if ~ isempty ( param_names ) && ~ iscellstr ( param_names )
error ( ' The 12th argument, if provided, must be a cellstr' )
end
end
2019-03-26 15:36:55 +01:00
if nargin < = 12
ds_range = ds . dates ;
else
if isempty ( ds_range )
ds_range = ds . dates ;
else
if ds_range ( 1 ) < ds . firstdate || ds_range ( end ) > lastdate ( ds )
error ( ' There is a problem with the 13th argument: the date range does not correspond to that of the dseries' )
end
end
end
2019-01-11 11:34:48 +01:00
%% Parse equation
2019-03-26 15:36:55 +01:00
[ Y , lhssub , X , fp , lp ] = common_parsing ( ds ( ds_range ) , get_ast ( { eqtag } ) , true , param_names ) ;
2019-01-16 12:38:02 +01:00
lhsname = Y { 1 } . name ;
Y = Y { 1 } . data ;
2019-01-11 11:34:48 +01:00
X = X { 1 } ;
fp = fp { 1 } ;
lp = lp { 1 } ;
pnames = X . name ;
2019-01-16 12:38:02 +01:00
N = size ( X . data , 1 ) ;
X = X . data ;
2018-09-28 16:09:50 +02:00
2019-01-11 11:34:48 +01:00
%% Estimation (see Koop, Gary. Bayesian Econometrics. 2003. Chapter 4.2)
2018-10-06 16:08:28 +02:00
PosteriorDegreesOfFreedom = N + nu ;
n = length ( pnames ) ;
assert ( n == length ( BetaPriorExpectation ) , ' the length prior mean for beta must be the same as the number of parameters in the equation to be estimated.' ) ;
h = 1.0 / s2 * nu ; % Initialize h to the prior expectation.
2018-09-28 16:09:50 +02:00
2018-10-06 16:08:28 +02:00
periods = 1 ;
linee = 1 ;
2018-09-28 16:09:50 +02:00
% Posterior Simulation
2019-02-28 11:49:50 +01:00
oo_ . olsgibbs . ( model_name ) . draws = zeros ( floor ( ( ndraws - discarddraws ) / thin ) , n + 3 ) ;
2018-10-06 16:08:28 +02:00
for i = 1 : discarddraws
% Set conditional distribution of β
InverseConditionalPoseriorVariance = BetaInversePriorVariance + h * ( X ' * X ) ;
cholConditionalPosteriorVariance = chol ( InverseConditionalPoseriorVariance \ eye ( n ) , ' upper' ) ;
ConditionalPosteriorExpectation = ( BetaPriorExpectation * BetaInversePriorVariance + h * ( Y ' * X ) ) / InverseConditionalPoseriorVariance ;
% Draw beta | Y, h
beta = rand_multivariate_normal ( ConditionalPosteriorExpectation , cholConditionalPosteriorVariance , n ) ;
% draw h | Y, beta
resids = Y - X * transpose ( beta ) ;
s2_ = ( resids ' * resids + nu * s2 ) / PosteriorDegreesOfFreedom ;
h = gamrnd ( PosteriorDegreesOfFreedom / 2.0 , 2.0 / ( PosteriorDegreesOfFreedom * s2_ ) ) ;
end
2018-09-28 16:09:50 +02:00
2019-12-10 17:03:22 +01:00
hh = dyn_waitbar ( 0 , ' Please wait. Gibbs sampler...' ) ;
set ( hh , ' Name' , ' Olsgibbs estimation.' ) ;
2018-10-06 16:08:28 +02:00
for i = discarddraws + 1 : ndraws
2019-12-10 17:03:22 +01:00
if ~ mod ( i , 10 )
2020-04-06 09:13:36 +02:00
dyn_waitbar ( ( i - discarddraws ) / ( ndraws - discarddraws ) , hh , ' Please wait. Gibbs sampler...' ) ;
2019-12-10 17:03:22 +01:00
end
2018-10-06 16:08:28 +02:00
% Set conditional distribution of β
InverseConditionalPoseriorVariance = BetaInversePriorVariance + h * ( X ' * X ) ;
cholConditionalPosteriorVariance = chol ( InverseConditionalPoseriorVariance \ eye ( n ) , ' upper' ) ;
ConditionalPosteriorExpectation = ( BetaPriorExpectation * BetaInversePriorVariance + h * ( Y ' * X ) ) / InverseConditionalPoseriorVariance ;
% Draw beta | Y, h
beta = rand_multivariate_normal ( ConditionalPosteriorExpectation , cholConditionalPosteriorVariance , n ) ;
2018-09-28 16:09:50 +02:00
% draw h | Y, beta
2018-10-06 16:08:28 +02:00
resids = Y - X * transpose ( beta ) ;
s2_ = ( resids ' * resids + nu * s2 ) / PosteriorDegreesOfFreedom ;
h = gamrnd ( PosteriorDegreesOfFreedom / 2.0 , 2.0 / ( PosteriorDegreesOfFreedom * s2_ ) ) ;
R2 = 1 - var ( resids ) / var ( Y ) ;
if isequal ( periods , thin )
2019-02-28 11:49:50 +01:00
oo_ . olsgibbs . ( model_name ) . draws ( linee , 1 : n ) = beta ;
oo_ . olsgibbs . ( model_name ) . draws ( linee , n + 1 ) = h ;
oo_ . olsgibbs . ( model_name ) . draws ( linee , n + 2 ) = s2_ ;
oo_ . olsgibbs . ( model_name ) . draws ( linee , n + 3 ) = R2 ;
2018-10-06 16:08:28 +02:00
periods = 1 ;
linee = linee + 1 ;
else
periods = periods + 1 ;
2018-09-28 16:09:50 +02:00
end
end
2019-12-10 17:03:22 +01:00
dyn_waitbar_close ( hh ) ;
2018-09-28 16:09:50 +02:00
2019-01-11 11:34:48 +01:00
%% Save posterior moments.
2019-02-28 11:49:50 +01:00
oo_ . olsgibbs . ( model_name ) . posterior . mean . beta = mean ( oo_ . olsgibbs . ( model_name ) . draws ( : , 1 : n ) ) ' ;
oo_ . olsgibbs . ( model_name ) . posterior . mean . h = mean ( oo_ . olsgibbs . ( model_name ) . draws ( : , n + 1 ) ) ;
oo_ . olsgibbs . ( model_name ) . posterior . variance . beta = cov ( oo_ . olsgibbs . ( model_name ) . draws ( : , 1 : n ) ) ;
oo_ . olsgibbs . ( model_name ) . posterior . variance . h = var ( oo_ . olsgibbs . ( model_name ) . draws ( : , n + 1 ) ) ;
oo_ . olsgibbs . ( model_name ) . s2 = mean ( oo_ . olsgibbs . ( model_name ) . draws ( : , n + 2 ) ) ;
oo_ . olsgibbs . ( model_name ) . R2 = mean ( oo_ . olsgibbs . ( model_name ) . draws ( : , n + 3 ) ) ;
2018-09-28 16:09:50 +02:00
2018-10-24 17:13:06 +02:00
% Yhat
idx = 0 ;
yhatname = [ eqtag ' _olsgibbs_FIT' ] ;
if ~ isempty ( fitted_names_dict )
idx = strcmp ( fitted_names_dict ( : , 1 ) , eqtag ) ;
if any ( idx )
yhatname = fitted_names_dict { idx , 2 } ;
end
end
2019-12-02 16:22:11 +01:00
oo_ . olsgibbs . ( model_name ) . Yhat = dseries ( X * oo_ . olsgibbs . ( model_name ) . posterior . mean . beta , fp , yhatname ) ;
2019-12-10 14:08:50 +01:00
oo_ . olsgibbs . ( model_name ) . YhatOrig = oo_ . olsgibbs . ( model_name ) . Yhat ;
oo_ . olsgibbs . ( model_name ) . Yobs = dseries ( Y , fp , lhsname ) ;
2019-12-02 16:22:11 +01:00
% Residuals
2019-12-03 10:50:11 +01:00
oo_ . olsgibbs . ( model_name ) . resid = Y - oo_ . olsgibbs . ( model_name ) . Yhat ;
2018-10-24 17:13:06 +02:00
% Apply correcting function for Yhat if it was passed
2019-12-02 16:22:11 +01:00
oo_ . olsgibbs . ( model_name ) . Yhat = oo_ . olsgibbs . ( model_name ) . Yhat + lhssub { 1 } ;
2018-10-24 17:13:06 +02:00
if any ( idx ) ...
&& length ( fitted_names_dict ( idx , : ) ) == 3 ...
&& ~ isempty ( fitted_names_dict { idx , 3 } )
2019-02-28 11:49:50 +01:00
oo_ . olsgibbs . ( model_name ) . Yhat = ...
feval ( fitted_names_dict { idx , 3 } , oo_ . olsgibbs . ( model_name ) . Yhat ) ;
2018-10-24 17:13:06 +02:00
end
2021-02-04 16:54:09 +01:00
ds . ( oo_ . olsgibbs . ( model_name ) . Yhat . name { : } ) = oo_ . olsgibbs . ( model_name ) . Yhat ;
2018-10-24 17:13:06 +02:00
2018-10-06 16:08:28 +02:00
% Compute and save posterior densities.
for i = 1 : n
2019-02-28 11:49:50 +01:00
xx = oo_ . olsgibbs . ( model_name ) . draws ( : , i ) ;
2018-10-06 16:08:28 +02:00
nn = length ( xx ) ;
bandwidth = mh_optimal_bandwidth ( xx , nn , 0 , ' gaussian' ) ;
[ x , f ] = kernel_density_estimate ( xx , 512 , nn , bandwidth , ' gaussian' ) ;
2019-02-28 11:49:50 +01:00
oo_ . olsgibbs . ( model_name ) . posterior . density . ( pnames { i } ) = [ x , f ] ;
2018-09-28 16:09:50 +02:00
end
2018-10-06 16:08:28 +02:00
% Update model's parameters with posterior mean.
2019-02-26 14:32:19 +01:00
idxs = zeros ( length ( pnames ) , 1 ) ;
2018-10-06 16:08:28 +02:00
for j = 1 : length ( pnames )
2019-02-26 14:32:19 +01:00
idxs ( j ) = find ( strcmp ( M_ . param_names , pnames { j } ) ) ;
2019-02-28 11:49:50 +01:00
M_ . params ( idxs ( j ) ) = oo_ . olsgibbs . ( model_name ) . posterior . mean . beta ( j ) ;
2018-09-28 16:09:50 +02:00
end
2019-05-10 14:50:44 +02:00
oo_ . olsgibbs . ( model_name ) . pnames = pnames ;
2018-09-28 16:09:50 +02:00
2019-02-26 14:32:19 +01:00
% Write .inc file
2019-02-28 11:49:50 +01:00
write_param_init_inc_file ( ' olsgibbs' , model_name , idxs , oo_ . olsgibbs . ( model_name ) . posterior . mean . beta ) ;
2019-02-26 14:32:19 +01:00
2019-01-11 11:34:48 +01:00
%% Print Output
2018-10-06 16:08:28 +02:00
if ~ options_ . noprint
ttitle = [ ' Bayesian estimation (with Gibbs sampling) of equation ' ' ' eqtag ' ' ' ' ] ;
2019-01-16 12:38:02 +01:00
preamble = { [ ' Dependent Variable: ' lhsname { : } ] , ...
2018-10-06 16:08:28 +02:00
sprintf ( ' No. Independent Variables: %d' , size ( X , 2 ) ) , ...
sprintf ( ' Observations: %d from %s to %s\n' , size ( X , 1 ) , fp . char , lp . char ) } ;
2019-02-28 11:49:50 +01:00
afterward = { sprintf ( ' s^2: %f' , oo_ . olsgibbs . ( model_name ) . s2 ) , sprintf ( ' R^2: %f' , oo_ . olsgibbs . ( model_name ) . R2 ) } ;
dyn_table ( ttitle , preamble , afterward , pnames , { ' Posterior mean' , ' Posterior std.' } , 4 , [ oo_ . olsgibbs . ( model_name ) . posterior . mean . beta , sqrt ( diag ( oo_ . olsgibbs . ( model_name ) . posterior . variance . beta ) ) ] ) ;
2019-01-11 11:34:48 +01:00
end
end