diff --git a/matlab/+mom/run.m b/matlab/+mom/run.m
index 93f5e2335f432b92331d6772e5961719f49adaf0..01aa83c182fcd0116cf05d074b858a95fb2506fd 100644
--- a/matlab/+mom/run.m
+++ b/matlab/+mom/run.m
@@ -62,7 +62,7 @@ function [oo_, options_mom_, M_] = run(bayestopt_, options_, oo_, estim_params_,
% o check_mode_file
% o check_posterior_sampler_options
% o check_prior_bounds
-% o check_prior_stderr_corr
+% o check_prior_stderr_corr_skew
% o check_steady_state_changes_parameters
% o check_varobs_are_endo_and_declared_once
% o check_hessian_at_the_mode
@@ -404,7 +404,7 @@ if do_bayesian_estimation_mcmc
end
% warning if prior allows that stderr parameters are negative or corr parameters are outside the unit circle
if do_bayesian_estimation
- check_prior_stderr_corr(estim_params_,bayestopt_);
+ check_prior_stderr_corr_skew(estim_params_,bayestopt_);
% check value of prior density
[~,~,~,info] = priordens(xparam0,bayestopt_.pshape,bayestopt_.p6,bayestopt_.p7,bayestopt_.p3,bayestopt_.p4);
if any(info)
diff --git a/matlab/check_model.m b/matlab/check_model.m
index 7dc297109915cfc12379554d479e6df690885276..2e940c2a41507016f86453bf0382b795ca6631a9 100644
--- a/matlab/check_model.m
+++ b/matlab/check_model.m
@@ -36,4 +36,9 @@ if ~isequal(M_.H,0)
if ~check_consistency_covariances(M_.H)
error('The specified covariances for the measurement errors are not consistent with the variances as they imply a correlation larger than +-1')
end
+end
+
+idxSkewOutOfBounds = (abs(M_.Skew_e) > 0.995);
+if any(idxSkewOutOfBounds)
+ error('Skewness parameters for %s are larger than theoretical bound of ±0.995.',strjoin(M_.exo_names(idxSkewOutOfBounds),','))
end
\ No newline at end of file
diff --git a/matlab/default_option_values.m b/matlab/default_option_values.m
index 300b4c6114fdac8e970ca01769904bd75ffda4bd..20587b6794f17fddee05e72da096b0e8f84b250f 100644
--- a/matlab/default_option_values.m
+++ b/matlab/default_option_values.m
@@ -394,6 +394,8 @@ options_.kalman_algo = 0;
options_.kalman_filter_mex = false;
options_.fast_kalman_filter = false;
options_.kalman_tol = 1e-10;
+options_.kalman.pskf.prune_tol = 0.01;
+options_.kalman.pskf.use_in_gaussian_case = false;
options_.kalman.keep_kalman_algo_if_singularity_is_detected = false;
options_.diffuse_kalman_tol = 1e-6;
options_.use_univariate_filters_if_singularity_is_detected = 1;
@@ -555,6 +557,7 @@ options_.homotopy_force_continue = false;
% numerical hessian
hessian.use_penalized_objective = false;
+hessian.use_onesided = false;
% Robust prediction error covariance (kalman filter)
options_.rescale_prediction_error_covariance = false;
diff --git a/matlab/distributions/csnMean.m b/matlab/distributions/csnMean.m
new file mode 100644
index 0000000000000000000000000000000000000000..afef403ebb6788100c247a451713a13d82458127
--- /dev/null
+++ b/matlab/distributions/csnMean.m
@@ -0,0 +1,126 @@
+function [E_csn] = csnMean(mu, Sigma, Gamma, nu, Delta, cdfmvna_fct)
+% function [E_csn] = csnMean(mu, Sigma, Gamma, nu, Delta, cdfmvna_fct)
+% -------------------------------------------------------------------------
+% Evaluates the unconditional expectation vector of a CSN(mu,Sigma,Gamma,nu,Delta)
+% distributed random variable based on expressions derived in
+% - Dominguez-Molina, Gonzalez-Farias, Gupta (2003, sec. 2.4) - The multivariate closed skew normal distribution
+% - Rodenburger (2015) - On Approximations of Normal Integrals in the Context of Probit based Discrete Choice Modeling (2015)
+% -------------------------------------------------------------------------
+% INPUTS
+% - mu [p by 1] 1st parameter of the CSN distribution (location)
+% - Sigma [p by p] 2nd parameter of the CSN distribution (scale)
+% - Gamma [q by p] 3rd parameter of the CSN distribution (regulates skewness from Gaussian (Gamma=0) to half normal)
+% - nu [q by 1] 4th parameter of the CSN distribution (enables closure of CSN distribution under conditioning)
+% - Delta [q by q] 5th parameter of the CSN distribution (enables closure of CSN distribution und marginalization)
+% - cdfmvna_fct [string] name of function to compute log Gaussian cdf, possible values: 'logmvncdf_ME', 'mvncdf', 'qsilatmvnv', 'qsimvnv'
+% -------------------------------------------------------------------------
+% OUTPUTS
+% - E_csn [p by 1] expectation vector of the CSN(mu,Sigma,Gamma,nu,Delta) distribution
+% =========================================================================
+% Copyright (C) 2022-2023 Gaygysyz Guljanov
+%
+% This 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.
+%
+% This 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.
+% -------------------------------------------------------------------------
+% This file is part of the replication files for the paper "Pruned Skewed
+% Kalman Filter and Smoother: Pruned Skewed Kalman Filter and Smoother:
+% With Applications to the Yield Curve and Asymmetric Monetary Policy Shocks"
+% by Gaygysyz Guljanov, Willi Mutschler, Mark Trede
+% =========================================================================
+
+if nargin < 6
+ cdfmvna_fct = 'logmvncdf_ME';
+end
+q = size(Delta, 1); % skewness dimension
+derVect = nan(q, 1); % initialize gradient vector
+Delta1 = Delta + Gamma * Sigma * Gamma'; % var-cov matrix inside the expression "psi" of equation (1) of Dominguez-Molina, Gonzalez-Farias, Gupta (2003, p. 10)
+
+% differentiate the normal cdf, steps are given in Rodenburger (2015, p.28-29)
+for jj = 1:q
+ % permutation of var-cov matrix
+ logi = true(q, 1); logi(jj, 1)=false;
+ nu2 = nan(q, 1); nu2(1:q-1)=nu(logi); nu2(q)=nu(jj);
+ Delta2 = nan(size(Delta1));
+ Delta2(1:q-1, :) = Delta1(logi, :);
+ Delta2(q, :) = Delta1(jj, :);
+ Delta2 = [Delta2(:, logi), Delta2(:, jj)];
+
+ % conditional mean and var-cov matrix
+ nu2 = -nu2;
+ condMean = (Delta2(1:q-1, q)*nu2(q))/Delta2(q, q);
+ condVar = Delta2(1:q-1, 1:q-1) - (Delta2(1:q-1, q)/Delta2(q, q))*Delta2(q, 1:q-1);
+
+ % evaluate the gradient vector
+ term1 = normpdf(nu2(q), 0, sqrt(Delta2(q, q)));
+ if isempty(condVar)
+ term2 = 1;
+ else
+ % different functions to evaluate log Gaussian cdf
+ if cdfmvna_fct == "logmvncdf_ME"
+ % requires zero mean and correlation matrix as inputs
+ normalization2 = diag(1./sqrt(diag(condVar)));
+ eval_point2 = normalization2*(nu2(1:q-1) - condMean);
+ Corr_mat2 = normalization2*condVar*normalization2;
+ Corr_mat2 = 0.5*(Corr_mat2 + Corr_mat2');
+ term2 = exp(logmvncdf_ME(eval_point2, Corr_mat2));
+ elseif cdfmvna_fct == "mvncdf"
+ term2 = mvncdf(nu2(1:q-1)', condMean', condVar);
+ elseif cdfmvna_fct == "qsilatmvnv"
+ if (q-1) > 1
+ term2 = qsilatmvnv( 1000*(q-1), condVar, repmat(-Inf,q-1,1)-condMean, nu2(1:q-1)-condMean );
+ else
+ term2 = normcdf(nu2, condMean, condVar);
+ end
+ elseif cdfmvna_fct == "qsimvnv"
+ if (q-1) > 1
+ term2 = qsimvnv( 1000*(q-1), condVar, repmat(-Inf,q-1,1)-condMean, nu2(1:q-1)-condMean );
+ else
+ term2 = normcdf(nu2, condMean, condVar);
+ end
+ end
+ end
+ derVect(jj) = term1*term2;
+end
+
+% evaluate end result using equation (1) of Dominguez-Molina, Gonzalez-Farias, Gupta (2003, p. 10)
+Delta1 = 0.5*(Delta1+Delta1');
+if isempty(Gamma) || isempty(nu) || isempty(Delta1)
+ if isempty(Gamma) && isempty(nu) && isempty(Delta1)
+ term3 = 1;
+ else
+ error("Problem with Gamma, nu, Delta being empty / not empty")
+ end
+else
+ if cdfmvna_fct == "logmvncdf_ME"
+ % requires zero mean and correlation matrix as inputs
+ normalization3 = diag(1./sqrt(diag(Delta1)));
+ eval_point3 = normalization3*(zeros(q, 1) - nu);
+ Corr_mat3 = normalization3*Delta1*normalization3;
+ Corr_mat3 = 0.5*(Corr_mat3 + Corr_mat3');
+ term3 = exp(logmvncdf_ME(eval_point3, Corr_mat3));
+ elseif cdfmvna_fct == "mvncdf"
+ term3 = mvncdf(zeros(1, q), nu', Delta1);
+ elseif cdfmvna_fct == "qsilatmvnv"
+ if q > 1
+ term3 = qsilatmvnv( 1000*q, Delta1, repmat(-Inf,q,1)-nu, zeros(q, 1)-nu );
+ else
+ term3 = normcdf(0, nu, Delta1);
+ end
+ elseif cdfmvna_fct == "qsimvnv"
+ if q > 1
+ term3 = qsimvnv( 1000*q, Delta1, repmat(-Inf,q,1)-nu, zeros(q, 1)-nu );
+ else
+ term3 = normcdf(0, nu, Delta1);
+ end
+ end
+end
+E_csn = mu + Sigma * Gamma' * (derVect/term3); % equation (1) of Dominguez-Molina, Gonzalez-Farias, Gupta (2003, p. 10)
+
+end
\ No newline at end of file
diff --git a/matlab/distributions/csnQuantile.m b/matlab/distributions/csnQuantile.m
new file mode 100644
index 0000000000000000000000000000000000000000..330b42df4107440ff00f4a9f2968d5af72b884a6
--- /dev/null
+++ b/matlab/distributions/csnQuantile.m
@@ -0,0 +1,126 @@
+function [q_alph] = csnQuantile(alph, mu, Sigma, Gamma, nu, Delta, cdfmvna_fct, optim_fct)
+% function [q_alph] = csnQuantile(alph, mu, Sigma, Gamma, nu, Delta, cdfmvna_fct, optim_fct)
+% -------------------------------------------------------------------------
+% Evaluates the quantile of a CSN(mu,Sigma,Gamma,nu,Delta) distributed random variable
+% for the cumulative probability alph in the interval [0,1]
+% Idea:
+% - q_alph = min{x: F(x) >= alph}
+% - Lemma 2.2.1 of chapter 2 of Genton (2004) - "Skew-elliptical distributions
+% and their applications: a journey beyond normality" gives the cdf, F_X(x)
+% of a CSN distributed random variable X
+% - We need to solve F_X(q_alpha) = alpha, to find alpha quantile q_alpha
+% - Minimization problem: q_alpha = argmin(log(F_X(q_alpha)) - log(alpha))
+% -------------------------------------------------------------------------
+% INPUTS
+% - alph [0<alph<1] cumulative probability alph in the interval [0,1]
+% - mu [p by 1] 1st parameter of the CSN distribution (location)
+% - Sigma [p by p] 2nd parameter of the CSN distribution (scale)
+% - Gamma [q by p] 3rd parameter of the CSN distribution (regulates skewness from Gaussian (Gamma=0) to half normal)
+% - nu [q by 1] 4th parameter of the CSN distribution (enables closure of CSN distribution under conditioning)
+% - Delta [q by q] 5th parameter of the CSN distribution (enables closure of CSN distribution und marginalization)
+% - cdfmvna_fct [string] name of function to compute log Gaussian cdf, possible values: 'logmvncdf_ME', 'mvncdf', 'qsilatmvnv', 'qsimvnv'
+% - optim_fct [string] name of minimization function, possible values: 'fminunc', 'lsqnonlin'
+% -------------------------------------------------------------------------
+% OUTPUTS
+% - q_alph [p by 1] alpha quantile vector of the CSN(mu,Sigma,Gamma,nu,Delta) distribution
+% =========================================================================
+% Copyright (C) 2022-2023 Gaygysyz Guljanov, Willi Mutschler
+%
+% This 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.
+%
+% This 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.
+% -------------------------------------------------------------------------
+% This file is part of the replication files for the paper "Pruned Skewed
+% Kalman Filter and Smoother: Pruned Skewed Kalman Filter and Smoother:
+% With Applications to the Yield Curve and Asymmetric Monetary Policy Shocks"
+% by Gaygysyz Guljanov, Willi Mutschler, Mark Trede
+% =========================================================================
+
+if nargin < 7 || (nargin >= 7 && isempty(cdfmvna_fct))
+ cdfmvna_fct = 'logmvncdf_ME';
+end
+if strcmp(cdfmvna_fct,'qsilatmvnv')
+ warning('csnQuantile: qsilatmvnv not yet working, switching to ''logmvncdf_ME'' to compute the quantile');
+ cdfmvna_fct = 'logmvncdf_ME';
+end
+if strcmp(cdfmvna_fct,'qsimvnv')
+ warning('csnQuantile: qsimvnv not yet working, switching to ''logmvncdf_ME'' to compute the quantile');
+ cdfmvna_fct = 'logmvncdf_ME';
+end
+if nargin < 8
+ optim_fct = 'lsqnonlin';
+end
+
+q = size(Delta, 1);
+Cov_mat = -Sigma * Gamma';
+
+% evaluate the first log-CDF
+Var_Cov1 = Delta - Gamma * Cov_mat;
+if cdfmvna_fct == "logmvncdf_ME"
+ normalization1 = diag(1 ./ sqrt(diag(Var_Cov1)));
+ eval_point1 = normalization1 * (zeros(q, 1) - nu);
+ Corr_mat1 = normalization1 * Var_Cov1 * normalization1;
+ Corr_mat1 = 0.5*(Corr_mat1 + Corr_mat1');
+ cdf1 = logmvncdf_ME(eval_point1, Corr_mat1);
+elseif cdfmvna_fct == "mvncdf"
+ cdf1 = log(mvncdf(zeros(1, q), nu', Var_Cov1));
+elseif cdfmvna_fct == "qsilatmvnv"
+ if q > 1
+ cdf1 = log(qsilatmvnv( 1000*q, Var_Cov1, repmat(-Inf,q,1)-nu, zeros(q, 1)-nu ));
+ else
+ cdf1 = log(normcdf(0, nu, Var_Cov1));
+ end
+elseif cdfmvna_fct == "qsimvnv"
+ if q > 1
+ cdf1 = log(qsimvnv( 1000*q, Var_Cov1, repmat(-Inf,q,1)-nu, zeros(q, 1)-nu ));
+ else
+ cdf1 = log(normcdf(0, nu, Var_Cov1));
+ end
+end
+
+% prepare helper function to get quantile
+Var_Cov2 = [Sigma, Cov_mat; Cov_mat', Var_Cov1];
+
+if cdfmvna_fct == "logmvncdf_ME"
+ % requires zero mean and correlation matrix as inputs
+ normalization2 = diag(1 ./ sqrt(diag(Var_Cov2)));
+ Corr_mat2 = normalization2 * Var_Cov2 * normalization2;
+ Corr_mat2 = 0.5 * (Corr_mat2 + Corr_mat2');
+ if strcmp(optim_fct,'lsqnonlin')
+ fun = @(x) logmvncdf_ME( (normalization2 * ([x; zeros(q, 1)] - [mu; nu])), Corr_mat2 ) - cdf1 - log(alph);
+ elseif strcmp(optim_fct,'fminunc')
+ fun = @(x) ( logmvncdf_ME( (normalization2 * ([x; zeros(q, 1)] - [mu; nu])), Corr_mat2 ) - cdf1 - log(alph))^2;
+ end
+elseif cdfmvna_fct == "mvncdf"
+ if strcmp(optim_fct,'lsqnonlin')
+ fun = @(x) log(mvncdf([x', zeros(1, q)], [mu', nu'], Var_Cov2)) - cdf1 - log(alph);
+ elseif strcmp(optim_fct,'fminunc')
+ fun = @(x) (log(mvncdf([x', zeros(1, q)], [mu', nu'], Var_Cov2)) - cdf1 - log(alph))^2;
+ end
+elseif cdfmvna_fct == "qsilatmvnv"
+ error('csnQuantile: qsilatmvnv not yet working')
+ fun = @(pp_quant) (log(qsilatmvnv(1000*(size(pp_quant,1)+q),Var_Cov2,repmat(-Inf,size(pp_quant,1)+q,1)-[mu', nu']', [pp_quant', zeros(1, q)]' - [mu', nu']')) ...
+ - cdf1 ...
+ - log(alph))^2;
+elseif cdfmvna_fct == "qsimvnv"
+ error('csnQuantile: qsimvnv not yet working')
+ fun = @(pp_quant) (log(qsimvnv(1000*(size(pp_quant,1)+q),Var_Cov2,repmat(-Inf,size(pp_quant,1)+q,1)-[pp_quant', zeros(1, q)],[pp_quant', zeros(1, q)] - [mu', nu'])) ...
+ - cdf1 ...
+ - log(alph))^2;
+end
+
+% run minimization
+optim_opt = optimset('MaxIter', 10000, 'MaxFunEvals', 10000, 'Display', 'off');
+if strcmp(optim_fct,'lsqnonlin')
+ q_alph = lsqnonlin(fun, mu, [], [], optim_opt);
+elseif strcmp(optim_fct,'fminunc')
+ q_alph = fminunc(fun, mu, optim_opt);
+end
+
+end % csnQuantile
\ No newline at end of file
diff --git a/matlab/distributions/csnRandDirect.m b/matlab/distributions/csnRandDirect.m
new file mode 100644
index 0000000000000000000000000000000000000000..44e0761329bdd6910bd27925661a343c052dcba4
--- /dev/null
+++ b/matlab/distributions/csnRandDirect.m
@@ -0,0 +1,55 @@
+function [csn_draws] = csnRandDirect(n, mu, Sigma, Gamma, nu, Delta)
+% function [res] = csnRandDirect(n, mu, Sigma, D, nu, Delta)
+% -------------------------------------------------------------------------
+% Draws random numbers from CSN distribution using accept-reject method
+% adapted from the "rcsn" command in the "csn" package of the R language
+% -------------------------------------------------------------------------
+% INPUTS
+% - n [scalar] number of desired random draws
+% - mu [p by 1] 1st parameter of the CSN distribution (location)
+% - Sigma [p by p] 2nd parameter of the CSN distribution (scale)
+% - Gamma [q by p] 3rd parameter of the CSN distribution (regulates skewness from Gaussian (Gamma=0) to half normal)
+% - nu [q by 1] 4th parameter of the CSN distribution (enables closure of CSN distribution under conditioning)
+% - Delta [q by q] 5th parameter of the CSN distribution (enables closure of CSN distribution und marginalization)
+% -------------------------------------------------------------------------
+% OUTPUTS
+% - csn_draws [p by n] random draws of the CSN(mu,Sigma,Gamma,nu,Delta) distribution
+% =========================================================================
+% Copyright (C) 2015 GPL v2 Dmitry Pavlyuk, Eugene Girtcius (rcsn.R function in the "csn" R package)
+% Copyright (C) 2022-2023 Gaygysyz Guljanov
+%
+% This 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.
+%
+% This 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.
+% -------------------------------------------------------------------------
+% This file is part of the replication files for the paper "Pruned Skewed
+% Kalman Filter and Smoother: Pruned Skewed Kalman Filter and Smoother:
+% With Applications to the Yield Curve and Asymmetric Monetary Policy Shocks"
+% by Gaygysyz Guljanov, Willi Mutschler, Mark Trede
+% =========================================================================
+
+q = length(nu);
+v = nan(q, n);
+Sigma2 = Delta + Gamma * Sigma * Gamma';
+L = chol(Sigma2, 'lower');
+for i = 1:n
+ draws = -1;
+ while sum(draws >= 0) ~= q
+ draws = -nu + L * randn(q, 1);
+ end
+ v(:, i) = draws;
+end
+
+tmp = Sigma * Gamma' / Sigma2;
+Exp = mu + tmp * (v + nu);
+Var = Sigma - tmp * Gamma * Sigma;
+L = chol(Var, 'lower');
+csn_draws = Exp + L * randn(length(mu), n);
+
+end % csnRandDirect
\ No newline at end of file
diff --git a/matlab/distributions/csnSkewness_univariate.m b/matlab/distributions/csnSkewness_univariate.m
new file mode 100644
index 0000000000000000000000000000000000000000..30d412abdbbfc3a9fa22ae06876ca5743f7d6315
--- /dev/null
+++ b/matlab/distributions/csnSkewness_univariate.m
@@ -0,0 +1,38 @@
+function Skew = csnSkewness_univariate(Sigma, Gamma)
+% function Skew = csnSkewness_univariate(Sigma, Gamma)
+% -------------------------------------------------------------------------
+% computes the skewness coefficient of a univariate CSN distributed random
+% variable X ~ CSN(mu,Sigma,Gamma,nu,Delta) with nu=0 and Delta=1
+% according to equation 3.6 of Grabek, Klos, and Koloch (2011) - Skew-Normal
+% shocks in the linear state space form DSGE model
+% -------------------------------------------------------------------------
+% INPUTS
+% - Sigma [double] scale parameter of CSN distribution
+% - Gamma [double] skewness parameter of CSN distribution
+% -------------------------------------------------------------------------
+% OUTPUTS
+% - Skew [double] theoretical skewness coefficient of X
+% =========================================================================
+% Copyright (C) 2023 Gaygysyz Guljanov
+%
+% This 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.
+%
+% This 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.
+% -------------------------------------------------------------------------
+% This file is part of the replication files for the paper "Pruned Skewed
+% Kalman Filter and Smoother: Pruned Skewed Kalman Filter and Smoother:
+% With Applications to the Yield Curve and Asymmetric Monetary Policy Shocks"
+% by Gaygysyz Guljanov, Willi Mutschler, Mark Trede
+% =========================================================================
+term1 = (4 - pi) / 2;
+term2 = (sqrt(2/pi) * Gamma * Sigma / sqrt(1 + Gamma^2 * Sigma))^3;
+term3 = (Sigma - (2/pi) * (Gamma^2 * Sigma^2) / (1 + Gamma^2 * Sigma))^(3/2);
+Skew = term1 * term2 / term3;
+
+end
\ No newline at end of file
diff --git a/matlab/distributions/csnVar.m b/matlab/distributions/csnVar.m
new file mode 100644
index 0000000000000000000000000000000000000000..2e0949b0492b271f37785c3c7cb4a246ab0f8244
--- /dev/null
+++ b/matlab/distributions/csnVar.m
@@ -0,0 +1,99 @@
+function V_csn = csnVar(Sigma, Gamma, nu, Delta, cdfmvna_fct)
+% function V_csn = csnVar(Sigma, Gamma, nu, Delta, cdfmvna_fct)
+% -------------------------------------------------------------------------
+% Evaluates the unconditional covariance matrix of a CSN(mu,Sigma,Gamma,nu,Delta)
+% distributed random variable based on expressions derived in
+% - Dominguez-Molina, Gonzalez-Farias, Gupta (2003, sec. 2.4) - The multivariate closed skew normal distribution
+% -------------------------------------------------------------------------
+% INPUTS
+% - Sigma [p by p] 2nd parameter of the CSN distribution (scale)
+% - Gamma [q by p] 3rd parameter of the CSN distribution (regulates skewness from Gaussian (Gamma=0) to half normal)
+% - nu [q by 1] 4th parameter of the CSN distribution (enables closure of CSN distribution under conditioning)
+% - Delta [q by q] 5th parameter of the CSN distribution (enables closure of CSN distribution und marginalization)
+% - cdfmvna_fct [string] name of function to compute log Gaussian cdf, possible values: 'logmvncdf_ME', 'mvncdf', 'qsilatmvnv', 'qsimvnv'
+% -------------------------------------------------------------------------
+% OUTPUTS
+% - V_csn [p by p] covariance matrix of the CSN(mu,Sigma,Gamma,nu,Delta) distribution
+% =========================================================================
+% Copyright (C) 2022-2023 Gaygysyz Guljanov
+%
+% This 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.
+%
+% This 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.
+% -------------------------------------------------------------------------
+% This file is part of the replication files for the paper "Pruned Skewed
+% Kalman Filter and Smoother: Pruned Skewed Kalman Filter and Smoother:
+% With Applications to the Yield Curve and Asymmetric Monetary Policy Shocks"
+% by Gaygysyz Guljanov, Willi Mutschler, Mark Trede
+% =========================================================================
+
+if nargin < 5
+ cdfmvna_fct = 'logmvncdf_ME';
+end
+q = size(Delta, 1); % skewness dimension
+
+% Eq. (1) and (3) of Dominguez-Molina, Gonzalez-Farias, Gupta (2003, p. 10-11)
+Delta2 = Delta + Gamma * Sigma * Gamma';
+Delta2 = 0.5 * (Delta2 + Delta2'); %ensure symmetry
+% Evaluate first and second derivatives
+first_der = zeros(q, 1); % first derivatives
+second_der = zeros(q, q); % second derivatives
+for ii = 1:q
+ [term1, term2, evalp_t1, covar_t1, evalp_t2, covar_t2, mult_matr] = derivative_gaussian_cdf(nu, Delta2, ii, cdfmvna_fct);
+ first_der(ii) = term1 * term2;
+ for jj = 1:q
+ if ii ~= jj
+ which_element = jj - (ii < jj);
+ [expr1, expr2] = derivative_gaussian_cdf(-evalp_t2, covar_t2, which_element, cdfmvna_fct);
+ second_der(ii, jj) = term1 * expr1 * expr2;
+ else
+ expr3 = -evalp_t1 / covar_t1 * term1 * term2;
+ expr4 = zeros(1, q-1);
+ for kk = 1:q-1
+ [first_term, second_term] = derivative_gaussian_cdf(-evalp_t2, covar_t2, kk, cdfmvna_fct);
+ expr4(kk) = first_term * second_term;
+ end
+ expr5 = term1 * expr4 * (-mult_matr);
+ second_der(jj, jj) = expr3 + expr5;
+ end
+ end
+end
+
+% denominator in psi and Lambda of eq. (1) and (3) of Dominguez-Molina, Gonzalez-Farias, Gupta (2003, p. 10-11)
+if cdfmvna_fct == "logmvncdf_ME"
+ % requires zero mean and correlation matrix as inputs
+ normalization3 = diag(1 ./ sqrt(diag(Delta2)));
+ eval_point3 = normalization3 * (zeros(q, 1) - nu);
+ Corr_mat3 = normalization3 * Delta2 * normalization3;
+ Corr_mat3 = 0.5 * (Corr_mat3 + Corr_mat3');
+ term3 = exp(logmvncdf_ME(eval_point3, Corr_mat3));
+elseif cdfmvna_fct == "mvncdf"
+ term3 = mvncdf(zeros(1, q), nu', Delta2);
+elseif cdfmvna_fct == "qsilatmvnv"
+ if q > 1
+ term3 = qsilatmvnv( 1000*q, Delta2, repmat(-Inf,q,1)-nu, zeros(q,1)-nu );
+ else
+ term3 = normcdf(0, nu, Delta2);
+ end
+elseif cdfmvna_fct == "qsimvnv"
+ if q > 1
+ term3 = qsimvnv( 1000*q, Delta2, repmat(-Inf,q,1)-nu, zeros(q,1)-nu );
+ else
+ term3 = normcdf(0, nu, Delta2);
+ end
+end
+
+% definition of psi in eq (1) of Dominguez-Molina, Gonzalez-Farias, Gupta (2003, p. 10)
+psi = first_der / term3;
+% definition of Lambda in eq (3) of Dominguez-Molina, Gonzalez-Farias, Gupta (2003, p. 11)
+Lambda = second_der / term3; Lambda = 0.5 * (Lambda + Lambda');
+% eq (3) of Dominguez-Molina, Gonzalez-Farias, Gupta (2003, p. 11)
+V_csn = Sigma + Sigma * Gamma' * Lambda * Gamma * Sigma - Sigma * Gamma' * (psi * psi') * Gamma * Sigma;
+
+end %csnVar
\ No newline at end of file
diff --git a/matlab/distributions/csnVarSkew_To_SigmaGamma_univariate.m b/matlab/distributions/csnVarSkew_To_SigmaGamma_univariate.m
new file mode 100644
index 0000000000000000000000000000000000000000..ee53f1f04eaf694b52e68174e92944eec27ac1c5
--- /dev/null
+++ b/matlab/distributions/csnVarSkew_To_SigmaGamma_univariate.m
@@ -0,0 +1,55 @@
+function [Sigma, Gamma, error_indicator] = csnVarSkew_To_SigmaGamma_univariate(Var, Skew, verbose)
+% function [Sigma, Gamma, error_indicator] = csnVarSkew_To_SigmaGamma_univariate(Var, Skew, verbose)
+% -------------------------------------------------------------------------
+% recovers the Sigma and Gamma parameters from the Variance and Skewness of
+% a univariate CSN distributed random variable X ~ CSN(mu,Sigma,Gamma,nu,Delta)
+% where Var=V[X], Skew=skewness[X], nu=0, and Delta=1
+% this basically inverts the formulas provided in equation 3.6 of Grabek,
+% Klos, and Koloch (2011) - Skew-Normal shocks in the linear state space form DSGE model
+% -------------------------------------------------------------------------
+% INPUTS
+% - Var [double] unconditional variance of X
+% - Skew [double] theoretical skewness coefficient of X
+% - verbose [boolean] display details if transformation did not work
+% -------------------------------------------------------------------------
+% OUTPUTS
+% - Sigma [double] scale parameter of CSN distribution
+% - Gamma [double] skewness parameter of CSN distribution
+% - error_indicator [boolean] 1: if transformation failed (mostly if Gamma is extremely large and/or Sigma is extremely small)
+% =========================================================================
+% Copyright (C) 2023 Willi Mutschler
+%
+% This 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.
+%
+% This 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.
+% -------------------------------------------------------------------------
+% This file is part of the replication files for the paper "Pruned Skewed
+% Kalman Filter and Smoother: Pruned Skewed Kalman Filter and Smoother:
+% With Applications to the Yield Curve and Asymmetric Monetary Policy Shocks"
+% by Gaygysyz Guljanov, Willi Mutschler, Mark Trede
+% =========================================================================
+if nargin < 3
+ verbose = false;
+end
+error_indicator = 0;
+Sigma = (-2^(2/3) * (4-pi)^(1/3) * abs(Skew)^(2/3) + pi - 4) * Var / (pi - 4);
+Gamma = sign(Skew) * ( 2^(1/3) * (4-pi)^(2/3)*sqrt(pi)*sqrt( abs(Skew)^(2/3) / ( ( -2 * 2^(1/3) * (4-pi)^(2/3) * (pi-2) * abs(Skew)^(4/3) + 2^(2/3) * (4-pi)^(7/3) * abs(Skew)^(2/3) + 2 * (pi - 4)^2 ) * Var ) ) );
+
+if (abs(csnVar(Sigma,Gamma,0,1) - Var)>1e-13) || (abs(csnSkewness_univariate(Sigma,Gamma) - Skew) > 1e-7)
+ if verbose
+ fprintf('Sigma: %f\n',Sigma);
+ fprintf('Gamma: %f\n',Gamma);
+ fprintf('Var: %f\n',csnVar(Sigma,Gamma,0,1));
+ fprintf('Skew: %f\n',csnSkewness_univariate(Sigma,Gamma));
+ fprintf('abs(csnVar(Sigma,Gamma,0,1) - Var): %f\n',abs(csnVar(Sigma,Gamma,0,1) - Var));
+ fprintf('abs(csnSkewness_univariate(Sigma,Gamma) - Skew): %f\n',abs(csnSkewness_univariate(Sigma,Gamma) - Skew));
+ warning('could not transform the parameters')
+ end
+ error_indicator = 1;
+end
\ No newline at end of file
diff --git a/matlab/distributions/derivative_gaussian_cdf.m b/matlab/distributions/derivative_gaussian_cdf.m
new file mode 100644
index 0000000000000000000000000000000000000000..6df52464e069ab83db14a7199678488ed0390f26
--- /dev/null
+++ b/matlab/distributions/derivative_gaussian_cdf.m
@@ -0,0 +1,104 @@
+function [term1, term2, evalp_t1, covar_t1, evalp_t2, covar_t2, mult_matr] = derivative_gaussian_cdf(mu, Sigma, which_element, cdfmvna_fct)
+% function [term1, term2, evalp_t1, covar_t1, evalp_t2, covar_t2, mult_matr] = derivative_gaussian_cdf(mu, Sigma, which_element, cdfmvna_fct)
+% -------------------------------------------------------------------------
+% Differentiates the multivariate normal cdf according to steps given in
+% Rodenburger (2015, p.28-29) - On Approximations of Normal Integrals in the Context of Probit based Discrete Choice Modeling (2015)
+% - First derivative is equal to term1 * term2
+% - Output variables "evalp_t1, covar_t1, evalp_t2, covar_t2, mult_matr"
+% are used for differentiating the first derivative for the second time,
+% for details look at the double loop in the code of csnVar() function above
+% -------------------------------------------------------------------------
+% INPUTS
+% - mu [p by 1] mean of normal distribution
+% - Sigma [p by p] covariance matrix of normal distribution
+% - which_element [scalar] indicator (1,...,p) for which variable in multivariate normal cdf the derivative is computed for
+% - cdfmvna_fct [string] name of function to compute log Gaussian cdf, possible values: 'logmvncdf_ME', 'mvncdf', 'qsilatmvnv', 'qsimvnv'
+% -------------------------------------------------------------------------
+% OUTPUTS
+% - term1 [double] evaluated normal PDF, see the first part of the derivative
+% - term2 [double] evaluated normal CDF, see the second part of the derivative
+% - evalp_t1 [1 by 1] evaluation point of normal PDF in term1
+% - covar_t1 [1 by 1] variance term of normal PDF in term1
+% - evalp_t2 [(p-1) by 1] evaluation point of normal CDF in term2
+% - covar_t2 [(p-1) by (p-1)] covariance matrix of normal CDF in term2
+% - mult_matr [(p-1) by 1] auxiliary expression Sigma_12/Sigma_22, this expression shows up when evaluating conditional mean and conditional Variance out of partitions of mean and Variance of joint distribution
+% =========================================================================
+% Copyright (C) 2022-2023 Gaygysyz Guljanov
+%
+% This 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.
+%
+% This 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.
+% -------------------------------------------------------------------------
+% This file is part of the replication files for the paper "Pruned Skewed
+% Kalman Filter and Smoother: Pruned Skewed Kalman Filter and Smoother:
+% With Applications to the Yield Curve and Asymmetric Monetary Policy Shocks"
+% by Gaygysyz Guljanov, Willi Mutschler, Mark Trede
+% =========================================================================
+
+if nargin < 4
+ cdfmvna_fct = 'logmvncdf_ME';
+end
+
+p = length(mu);
+logic_vec = true(p, 1);
+logic_vec(which_element, 1) = false;
+
+% permute the matrices
+mu2 = zeros(p, 1);
+mu2(1:p-1) = mu(logic_vec);
+mu2(p) = mu(which_element);
+
+Sigma2 = zeros(size(Sigma));
+Sigma2(1:p-1, :) = Sigma(logic_vec, :);
+Sigma2(p, :) = Sigma(which_element, :);
+Sigma2 = [Sigma2(:, logic_vec), Sigma2(:, which_element)];
+
+% Conditional mean and conditional var-cov matrix
+eval_point = -mu2;
+mult_matr = Sigma2(1:p-1, p) / Sigma2(p, p);
+cond_mean = mult_matr * eval_point(p);
+cond_var = Sigma2(1:p-1, 1:p-1) - mult_matr * Sigma2(p, 1:p-1);
+
+evalp_t1 = eval_point(p);
+covar_t1 = Sigma2(p, p);
+
+evalp_t2 = eval_point(1:p-1) - cond_mean;
+covar_t2 = cond_var;
+
+% Evaluate the first term in the derivative
+term1 = normpdf(evalp_t1, 0, sqrt(covar_t1));
+
+% Evaluate the second term in the derivative
+if isempty(cond_var)
+ term2 = 1;
+else
+ if cdfmvna_fct == "logmvncdf_ME"
+ stdnrd = diag(1./sqrt(diag(cond_var)));
+ haspl = stdnrd * evalp_t2;
+ Uytgesme = stdnrd * cond_var * stdnrd;
+ Uytgesme = 0.5 * (Uytgesme + Uytgesme');
+ term2 = exp(logmvncdf_ME(haspl, Uytgesme));
+ elseif cdfmvna_fct == "mvncdf"
+ term2 = mvncdf(eval_point(1:p-1)', cond_mean', cond_var);
+ elseif cdfmvna_fct == "qsilatmvnv"
+ if (p-1) > 1
+ term2 = qsilatmvnv( 1000*(p-1), cond_var, repmat(-Inf,p-1,1)-cond_mean, eval_point(1:p-1)-cond_mean );
+ else
+ term2 = normcdf(eval_point(1:p-1), cond_mean, cond_var);
+ end
+ elseif cdfmvna_fct == "qsimvnv"
+ if (p-1) > 1
+ term2 = qsimvnv( 1000*(p-1), cond_var, repmat(-Inf,p-1,1)-cond_mean, eval_point(1:p-1)-cond_mean );
+ else
+ term2 = normcdf(eval_point(1:p-1), cond_mean, cond_var);
+ end
+ end
+end
+
+end % derivative_gaussian_cdf
\ No newline at end of file
diff --git a/matlab/distributions/logmvncdf_ME.m b/matlab/distributions/logmvncdf_ME.m
new file mode 100644
index 0000000000000000000000000000000000000000..86550738205309b84e1c266fb38452ba19202333
--- /dev/null
+++ b/matlab/distributions/logmvncdf_ME.m
@@ -0,0 +1,93 @@
+function log_cdf = logmvncdf_ME(Zj, R, cutoff)
+% log_cdf = logmvncdf_ME(Zj, R, cutoff)
+% -------------------------------------------------------------------------
+% Approximates Gaussian log(CDF) function according to Mendell and Elston (1974)
+% -------------------------------------------------------------------------
+% INPUTS
+% - Zj [n by 1] column vector of points where Gaussian CDF is evaluated at
+% - R [n by n] correlation matrix
+% - cutoff [2 by 1] optional threshold points at which values in Zj are too low/high to be evaluated, defaults to [6 38]
+% -------------------------------------------------------------------------
+% OUTPUTS
+% log_cdf [double] approximate value of Gaussian log(CDF)
+% =========================================================================
+% Copyright (C) 2015 Dietmar Bauer (original implementation)
+% Copyright (C) 2022-2023 Gaygysyz Guljanov
+%
+% This 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.
+%
+% This 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.
+% -------------------------------------------------------------------------
+% This file is part of the replication files for the paper "Pruned Skewed
+% Kalman Filter and Smoother: Pruned Skewed Kalman Filter and Smoother:
+% With Applications to the Yield Curve and Asymmetric Monetary Policy Shocks"
+% by Gaygysyz Guljanov, Willi Mutschler, Mark Trede
+% =========================================================================
+
+if nargin < 3
+ cutoff = [6 38]; % default values
+end
+% cutoff tails as these are too low to be evaluated
+Zj(Zj>cutoff(1)) = cutoff(1);
+Zj(Zj<-cutoff(1)) = -cutoff(1);
+n = length(Zj); % remaining dimension of Z
+
+% first element
+cdf_val = phid(Zj(1));
+pdf_val = phip(Zj(1));
+log_cdf = log(cdf_val); % perform all calculations in logs
+
+for jj = 1:(n-1)
+ ajjm1 = pdf_val / cdf_val;
+
+ % Update Zj and Rij
+ tZ = Zj + ajjm1 * R(:, 1); % update Zj
+ R_jj = R(:, 1) * R(1, :);
+ tRij = R - R_jj * ( ajjm1 + Zj(1) ) * ajjm1; % update Rij
+
+ % Convert Rij (i.e. Covariance matrix) to Correlation matrix
+ COV_TO_CORR = sqrt( diag(tRij) );
+ Zj = tZ ./ COV_TO_CORR;
+ R = tRij ./ (COV_TO_CORR * COV_TO_CORR');
+
+ % Cutoff those dimensions if they are too low to be evaluated
+ Zj(Zj>cutoff(2)) = cutoff(2); % use larger cutoff
+ Zj(Zj<-cutoff(2)) = -cutoff(2); % use larger cutoff
+
+ % Evaluate jj's probability
+ cdf_val = phid(Zj(2));
+ pdf_val = phip(Zj(2));
+
+ % Delete unnecessary parts of updated Zj and Corr_mat
+ last_el = n - jj + 1;
+ Zj = Zj(2:last_el);
+ R = R(2:last_el, 2:last_el);
+
+ % Overall probability
+ log_cdf = log_cdf + log(cdf_val);
+
+end
+
+end % logmvncdf_ME
+
+
+%%%%%%%%%%%%%%
+% Normal PDF %
+%%%%%%%%%%%%%%
+function p = phip(z)
+ p = exp(-z.^2/2)/sqrt(2*pi);
+end
+
+
+%%%%%%%%%%%%%%
+% Normal CDF %
+%%%%%%%%%%%%%%
+function p = phid(z)
+ p = erfc( -z/sqrt(2) )/2;
+end
diff --git a/matlab/distributions/singular_skew_normal_2_closed_skew_normal.m b/matlab/distributions/singular_skew_normal_2_closed_skew_normal.m
new file mode 100644
index 0000000000000000000000000000000000000000..0573dea017e5822107882e4cd8a9a86b17ef8a05
--- /dev/null
+++ b/matlab/distributions/singular_skew_normal_2_closed_skew_normal.m
@@ -0,0 +1,70 @@
+function [S_hat, mu_hat, Lambda_hat, Gamma_hat] = singular_skew_normal_2_closed_skew_normal(mu_bar, Sigma_bar, Gamma_bar, tol, hmdel)
+% function [S_hat, mu_hat, Lambda_hat, Gamma_hat] = singular_skew_normal_2_closed_skew_normal(mu_bar, Sigma_bar, Gamma_bar, tol, hmdel)
+% -------------------------------------------------------------------------
+% write Singular Skew Normal (SSN) distribution as a linear transformation
+% of proper Closed Skew Normal (CSN) distribution
+% -----------------------------------------------------------------------
+% INPUTS
+% mu_bar: first parameter of the SSN distribution
+% Sigma_bar: second parameter of the SSN distribution
+% Gamma_bar: third parameter of the SSN distribution
+% tol: tolerance level showing which values should be assumed to be numerically zero
+% hmdel: minimum how many dimensions should be deleted, if known already
+% -----------------------------------------------------------------------
+% OUTPUT
+% parameters of transformed SSN distribution that is CSN
+% S_hat: multiplying matrix, which has more rows than columns
+% mu_hat : first parameter of the proper CSN
+% Lambda_hat: second parameter of the proper CSN
+% Gamma_hat : third parameter of the proper CSN
+% =========================================================================
+% Copyright (C) 2023 Gaygysyz Guljanov
+%
+% This 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.
+%
+% This 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.
+% -------------------------------------------------------------------------
+% This file is part of the replication files for the paper "Pruned Skewed
+% Kalman Filter and Smoother: Pruned Skewed Kalman Filter and Smoother:
+% With Applications to the Yield Curve and Asymmetric Monetary Policy Shocks"
+% by Gaygysyz Guljanov, Willi Mutschler, Mark Trede
+% =========================================================================
+if sum(mu_bar ~= 0)
+ error("mu_bar should be a zero vector")
+end
+[S_mat, Lambda_mat] = schur(Sigma_bar); % get Schur decomposition
+% Create a hold_vec, which tells which dimensions to hold
+if hmdel > 0
+ % Reorder the diagonal elements of Lambda, in an increasing order
+ len = length(Lambda_mat);
+ diag_el = diag(Lambda_mat);
+ [~, min_indices] = mink(diag_el, len);
+ % Should we delete more than hmdel
+ hmdel_tmp = hmdel;
+ ii = 1;
+ while diag_el(min_indices(hmdel_tmp + ii)) / diag_el(min_indices(len)) <= tol
+ ii = ii + 1;
+ end
+ hmdel_tmp = hmdel_tmp + ii - 1;
+ if hmdel_tmp >= len
+ error("Sigma matrix should at least have rank one");
+ end
+ % Which vector to hold, i.e. to not cut
+ hold_vec = true(len, 1);
+ hold_vec(min_indices(1:hmdel_tmp)) = false;
+else
+ hold_vec = diag(Lambda_mat) > tol;
+end
+% Hold the relevant dimensions, i.e. delete the unnecessary ones
+S_hat = S_mat(:, hold_vec);
+mu_hat = zeros(sum(hold_vec), 1);
+Lambda_hat = Lambda_mat(hold_vec, hold_vec);
+Gamma_hat = Gamma_bar * S_hat;
+
+end % singular_skew_normal_2_closed_skew_normal
diff --git a/matlab/estimation/GetPosteriorParametersStatistics.m b/matlab/estimation/GetPosteriorParametersStatistics.m
index ebaf93666523353a3a2df5c5d54f5f594743395d..a5b3102ff9a01420ca605407879bce44c96f27ff 100644
--- a/matlab/estimation/GetPosteriorParametersStatistics.m
+++ b/matlab/estimation/GetPosteriorParametersStatistics.m
@@ -39,6 +39,7 @@ nvx = estim_params_.nvx;
nvn = estim_params_.nvn;
ncx = estim_params_.ncx;
ncn = estim_params_.ncn;
+nsx = estim_params_.nsx;
np = estim_params_.np ;
latexFolder = CheckPath('latex',M_.dname);
@@ -87,7 +88,7 @@ if np
skipline()
disp(type)
disp(tit2)
- ip = nvx+nvn+ncx+ncn+1;
+ ip = nvx+nvn+ncx+ncn+nsx+1;
for i=1:np
if options_.mh_replic || (options_.load_mh_file && ~options_.load_results_after_load_mh) || ishssmc(options_) || isonline(options_)
draws = getalldraws(ip);
@@ -302,6 +303,48 @@ if ncn
end
end
+if nsx
+ type = 'shocks_skew';
+ if TeX
+ fid = TeXBegin(latexFolder, FileName,6, 'skewness of structural shocks');
+ end
+ skipline()
+ disp('skewness of shocks')
+ disp(tit2)
+ ip = nvx+nvn+ncx+ncn+1;
+ for i=1:nsx
+ if options_.mh_replic || (options_.load_mh_file && ~options_.load_results_after_load_mh)
+ draws = getalldraws(ip);
+ [post_mean, post_median, post_var, hpd_interval, post_deciles, density] = posterior_moments(draws, options_.mh_conf_sig);
+ k = estim_params_.skew_exo(i,1);
+ name = M_.exo_names{k};
+ oo_ = Filloo(oo_, name, type, post_mean, hpd_interval, post_median, post_var, post_deciles, density);
+ M_.Skew_e(k) = post_mean;
+ else
+ try
+ k = estim_params_.skew_exo(i,1);
+ name = M_.exo_names{k};
+ [post_mean, hpd_interval, post_var] = Extractoo(oo_, name, type);
+ catch
+ draws = getalldraws(ip);
+ [post_mean, post_median, post_var, hpd_interval, post_deciles, density] = posterior_moments(draws, options_.mh_conf_sig);
+ k = estim_params_.var_exo(i,1);
+ name = M_.exo_names{k};
+ oo_ = Filloo(oo_, name, type, post_mean, hpd_interval, post_median, post_var, post_deciles, density);
+ M_.Skew_e(k) = post_mean;
+ end
+ end
+ dprintf(pformat, header_width, name, bayestopt_.p1(ip), post_mean, hpd_interval, pnames{bayestopt_.pshape(ip)+1}, bayestopt_.p2(ip));
+ if TeX
+ name = M_.exo_names_tex{k};
+ TeXCore(fid,name, pnames{bayestopt_.pshape(ip)+1}, bayestopt_.p1(ip), bayestopt_.p2(ip), post_mean, sqrt(post_var), hpd_interval);
+ end
+ ip = ip+1;
+ end
+ if TeX
+ TeXEnd(fid, 6, 'skewness of structural shocks');
+ end
+end
%
%% subfunctions:
diff --git a/matlab/estimation/PlotPosteriorDistributions.m b/matlab/estimation/PlotPosteriorDistributions.m
index 633ab1cf56a8363a581349259e6b7f9e891f070c..3d2d5bc211aba518109d7b34838371f3a6d77687 100644
--- a/matlab/estimation/PlotPosteriorDistributions.m
+++ b/matlab/estimation/PlotPosteriorDistributions.m
@@ -40,8 +40,9 @@ nvx = estim_params_.nvx;
nvn = estim_params_.nvn;
ncx = estim_params_.ncx;
ncn = estim_params_.ncn;
+nsx = estim_params_.nsx;
np = estim_params_.np ;
-npar = nvx+nvn+ncx+ncn+np;
+npar = nvx+nvn+ncx+ncn+nsx+np;
MaxNumberOfPlotPerFigure = 9;% The square root must be an integer!
nn = sqrt(MaxNumberOfPlotPerFigure);
@@ -109,8 +110,17 @@ for i=1:npar
if ~issmc(options_) && ~options_.mh_posterior_mode_estimation
pmod = oo_.posterior_mode.measurement_errors_corr.(name);
end
+ elseif i <= nvx+nvn+ncx+ncn+nsx
+ name = M_.exo_names{estim_params_.skew_exo(i-(nvx+nvn+ncx+ncn),1)};
+ x1 = oo_.posterior_density.shocks_skew.(name)(:,1);
+ f1 = oo_.posterior_density.shocks_skew.(name)(:,2);
+ oo_.prior_density.shocks_skew.(name)(:,1) = x2;
+ oo_.prior_density.shocks_skew.(name)(:,2) = f2;
+ if ~issmc(options_) && ~options_.mh_posterior_mode_estimation
+ pmod = oo_.posterior_mode.shocks_skew.(name);
+ end
else
- j = i - (nvx+nvn+ncx+ncn);
+ j = i - (nvx+nvn+ncx+ncn+nsx);
name = M_.param_names{estim_params_.param_vals(j,1)};
x1 = oo_.posterior_density.parameters.(name)(:,1);
f1 = oo_.posterior_density.parameters.(name)(:,2);
diff --git a/matlab/estimation/PosteriorIRF.m b/matlab/estimation/PosteriorIRF.m
index 00366f0900626d97f03458c858bb42d911dd0f5e..a50ac92d5f1b700cec0624c5f37a9e477b940457 100644
--- a/matlab/estimation/PosteriorIRF.m
+++ b/matlab/estimation/PosteriorIRF.m
@@ -67,7 +67,7 @@ end
irf_shocks_indx = getIrfShocksIndx(M_, options_);
% Set various parameters & Check or create directories
-npar = estim_params_.nvx+estim_params_.nvn+estim_params_.ncx+estim_params_.ncn+estim_params_.np ;
+npar = estim_params_.nvx+estim_params_.nvn+estim_params_.ncx+estim_params_.ncn+estim_params_.nsx+estim_params_.np ;
nvobs = dataset_.vobs;
gend = dataset_.nobs;
diff --git a/matlab/estimation/check_bounds_and_definiteness_estimation.m b/matlab/estimation/check_bounds_and_definiteness_estimation.m
index 54410b1d4a1b26c23b3ae88935f9c7a5aaa03c70..db9f397ef769a440f2c89d44f61598349cd3f405 100644
--- a/matlab/estimation/check_bounds_and_definiteness_estimation.m
+++ b/matlab/estimation/check_bounds_and_definiteness_estimation.m
@@ -105,3 +105,30 @@ if ~issquare(H) || estim_params_.ncn || isfield(estim_params_,'calibrated_covari
end
end
end
+
+% skewness transformation failed
+if any(M_.csn.error_indicator)
+ warning('Skewness transformation failed, check if M_.csn.Gamma_e is extremely large and/or M_.csn.Sigma_e is extremely small');
+ penalty = sum(abs(M_.Skew_e));
+ fval = Inf;
+ exit_flag = 0;
+ info(1) = 56;
+ info(4) = penalty;
+ return
+end
+
+% theoretical bound on skewness coefficient
+if any(diag(M_.csn.Gamma_e))
+ skewnessBound = abs((sqrt(2)*(pi-4))/(pi-2)^(3/2));
+ penalty = 0;
+ for jp = 1:M_.exo_nbr
+ penalty = penalty + max(0,abs(csnSkewness_univariate(M_.csn.Sigma_e(jp,jp),M_.csn.Gamma_e(jp,jp))) - skewnessBound);
+ end
+ if penalty > 0
+ fval = Inf;
+ exit_flag = 0;
+ info(1) = 57;
+ info(4) = penalty;
+ return
+ end
+end
\ No newline at end of file
diff --git a/matlab/estimation/display_estimation_results_table.m b/matlab/estimation/display_estimation_results_table.m
index f24729cb30b9cf76ddfa255e254cab0c9611c995..c5d139778d2578e70c5919fb13b9853b0b6cea77 100644
--- a/matlab/estimation/display_estimation_results_table.m
+++ b/matlab/estimation/display_estimation_results_table.m
@@ -41,8 +41,9 @@ nvx = estim_params_.nvx; % Variance of the structural innovations (number of pa
nvn = estim_params_.nvn; % Variance of the measurement innovations (number of parameters).
ncx = estim_params_.ncx; % Covariance of the structural innovations (number of parameters).
ncn = estim_params_.ncn; % Covariance of the measurement innovations (number of parameters).
+nsx = estim_params_.nsx; % Skewness of the structural innovations (number of parameters).
np = estim_params_.np ; % Number of deep parameters.
-nx = nvx+nvn+ncx+ncn+np; % Total number of parameters to be estimated.
+nx = nvx+nvn+ncx+ncn+nsx+np; % Total number of parameters to be estimated.
skipline()
disp(['RESULTS FROM ' upper(table_title) ' ESTIMATION'])
@@ -57,7 +58,7 @@ else
tit1 = sprintf('%-*s %10s %7s %6s\n', header_width, ' ', 'Estimate', 's.d.', 't-stat');
end
if np
- ip = nvx+nvn+ncx+ncn+1;
+ ip = nvx+nvn+ncx+ncn+nsx+1;
disp('parameters')
disp(tit1)
for i=1:np
@@ -173,6 +174,29 @@ if ncn
skipline()
end
+if nsx
+ disp('skewness of shocks')
+ disp(tit1)
+ ip = nvx+nvn+ncx+ncn+1;
+ for i=1:nsx
+ k = estim_params_.skew_exo(i,1);
+ name = sprintf('%s', M_.exo_names{k});
+ NAME = sprintf('%s', M_.exo_names{k});
+ if contains(field_name,'posterior')
+ fprintf('%-*s %10.4f %8.4f %7.4f %6s %6.4f \n', ...
+ header_width, name, bayestopt_.p1(ip), xparam1(ip), stdh(ip), ...
+ pnames{bayestopt_.pshape(ip)+1}, bayestopt_.p2(ip));
+ else
+ fprintf('%-*s %10.4f %7.4f %7.4f \n',header_width, name, xparam1(ip), ...
+ stdh(ip), tstath(ip));
+ end
+ oo_.(sprintf('%s_mode', field_name)).shocks_skew.(NAME) = xparam1(ip);
+ oo_.(sprintf('%s_std_at_mode', field_name)).shocks_skew.(NAME) = stdh(ip);
+ ip = ip+1;
+ end
+ skipline()
+end
+
if any(xparam1(1:nvx+nvn)<0)
warning(sprintf('Some estimated standard deviations are negative.\n Dynare internally works with variances so that the sign does not matter.\n Nevertheless, it is recommended to impose either prior restrictions (Bayesian Estimation)\n or a lower bound (ML) to assure positive values.'))
end
@@ -184,7 +208,7 @@ if any(bayestopt_.pshape > 0) && options_.TeX %% Bayesian estimation (posterior
filename = [latexDirectoryName '/' M_.fname '_Posterior_Mode_1.tex'];
fidTeX = fopen(filename,'w');
TeXBegin_Bayesian(fidTeX,1,'parameters')
- ip = nvx+nvn+ncx+ncn+1;
+ ip = nvx+nvn+ncx+ncn+nsx+1;
for i=1:np
fprintf(fidTeX,'$%s$ & %s & %7.3f & %6.4f & %8.4f & %7.4f \\\\ \n',...
M_.param_names_tex{estim_params_.param_vals(i,1)}, ...
@@ -271,12 +295,30 @@ if any(bayestopt_.pshape > 0) && options_.TeX %% Bayesian estimation (posterior
end
TeXEnd(fidTeX)
end
+ if nsx
+ TeXfile = [latexDirectoryName '/' M_.fname '_Posterior_Mode_6.tex'];
+ fidTeX = fopen(TeXfile,'w');
+ TeXBegin_Bayesian(fidTeX,6,'skeweness of structural shocks')
+ ip = nvx+nvn+ncx+ncn+1;
+ for i=1:nsx
+ k = estim_params_.skew_exo(i,1);
+ fprintf(fidTeX,'$%s$ & %4s & %7.3f & %6.4f & %8.4f & %7.4f \\\\ \n',...
+ M_.exo_names_tex{k}, ...
+ pnames{bayestopt_.pshape(ip)+1}, ...
+ bayestopt_.p1(ip), ...
+ bayestopt_.p2(ip), ...
+ xparam1(ip), ...
+ stdh(ip));
+ ip = ip+1;
+ end
+ TeXEnd(fidTeX)
+ end
elseif all(bayestopt_.pshape == 0) && options_.TeX %% MLE and GMM Latex output
if np
filename = [latexDirectoryName '/' M_.fname '_' LaTeXtitle '_Mode_1.tex'];
fidTeX = fopen(filename, 'w');
TeXBegin_ML(fidTeX, 1, 'parameters', table_title, LaTeXtitle)
- ip = nvx+nvn+ncx+ncn+1;
+ ip = nvx+nvn+ncx+ncn+nsx+1;
for i=1:np
fprintf(fidTeX,'$%s$ & %8.4f & %7.4f & %7.4f\\\\ \n',...
M_.param_names_tex{estim_params_.param_vals(i,1)}, ...
@@ -353,6 +395,22 @@ elseif all(bayestopt_.pshape == 0) && options_.TeX %% MLE and GMM Latex output
end
TeXEnd(fidTeX)
end
+ if nsx
+ filename = [latexDirectoryName '/' M_.fname '_' LaTeXtitle '_Mode_6.tex'];
+ fidTeX = fopen(filename, 'w');
+ TeXBegin_ML(fidTeX, 6, 'skewness of structural shocks', table_title, LaTeXtitle)
+ ip = nvx+nvn+ncx+ncn+1;
+ for i=1:nsx
+ k = estim_params_.skew_exo(i,1);
+ fprintf(fidTeX, '$%s$ & %8.4f & %7.4f & %7.4f \\\\ \n', ...
+ M_.exo_names_tex{k}, ...
+ xparam1(ip), ...
+ stdh(ip), ...
+ tstath(ip));
+ ip = ip+1;
+ end
+ TeXEnd(fidTeX)
+ end
end
diff --git a/matlab/estimation/do_parameter_initialization.m b/matlab/estimation/do_parameter_initialization.m
index ac6c0a0b2e0c230a874631f1c2c6ecd4ef4ec8d3..f16e5cdd6559467939ea352488c6a6f3388936d0 100644
--- a/matlab/estimation/do_parameter_initialization.m
+++ b/matlab/estimation/do_parameter_initialization.m
@@ -44,16 +44,18 @@ nvx = size(estim_params_.var_exo,1);
nvn = size(estim_params_.var_endo,1);
ncx = size(estim_params_.corrx,1);
ncn = size(estim_params_.corrn,1);
+nsx = size(estim_params_.skew_exo,1);
np = size(estim_params_.param_vals,1);
estim_params_.nvx = nvx; %exogenous shock variances
estim_params_.nvn = nvn; %endogenous variances, i.e. measurement error
estim_params_.ncx = ncx; %exogenous shock correlations
estim_params_.ncn = ncn; % correlation between endogenous variables, i.e. measurement error.
+estim_params_.nsx = nsx; % skewness parameters of independent shocks
estim_params_.np = np; % other parameters of the model
-xparam1_explicitly_initialized = NaN(nvx+nvn+ncx+ncn+np,1);
-xparam1_properly_calibrated = NaN(nvx+nvn+ncx+ncn+np,1);
+xparam1_explicitly_initialized = NaN(nvx+nvn+ncx+ncn+nsx+np,1);
+xparam1_properly_calibrated = NaN(nvx+nvn+ncx+ncn+nsx+np,1);
offset=0;
if nvx
@@ -125,6 +127,19 @@ if ncn
end
end
offset=offset+ncn;
+if nsx
+ initialized_par_index=find(~isnan(estim_params_.skew_exo(:,2)));
+ calibrated_par_index=find(isnan(estim_params_.skew_exo(:,2)) & ~isnan(xparam1_calib(offset+1:offset+nsx,1)));
+ uninitialized_par_index=find(isnan(estim_params_.skew_exo(:,2)) & isnan(xparam1_calib(offset+1:offset+nsx,1)));
+ xparam1_explicitly_initialized(offset+initialized_par_index,1) = estim_params_.skew_exo(initialized_par_index,2);
+ estim_params_.skew_exo(calibrated_par_index,2)=xparam1_calib(offset+calibrated_par_index,1);
+ xparam1_properly_calibrated(offset+calibrated_par_index,1) = xparam1_calib(offset+calibrated_par_index,1);
+ if uninitialized_par_index
+ fprintf('PARAMETER INITIALIZATION: Warning, some estimated skewness coefficients of shocks are not\n')
+ fprintf('PARAMETER INITIALIZATION: initialized. They will be initialized with the prior mean.\n')
+ end
+end
+offset=offset+nsx;
if np
initialized_par_index=find(~isnan(estim_params_.param_vals(:,2)));
calibrated_par_index=find(isnan(estim_params_.param_vals(:,2)) & ~isnan(xparam1_calib(offset+1:offset+np,1)));
diff --git a/matlab/estimation/dsge_likelihood.m b/matlab/estimation/dsge_likelihood.m
index 6a72e9e38a4b97c9e3325e36d49582a01e56713a..c85dea0fc7cf6fba19706e8052647469f86f56d8 100644
--- a/matlab/estimation/dsge_likelihood.m
+++ b/matlab/estimation/dsge_likelihood.m
@@ -77,6 +77,11 @@ if ~isempty(xparam1)
xparam1 = xparam1(:);
end
+% transform parameters back from unbounded to bounded (only during optimization)
+if estim_params_.transformed && any(estim_params_.transform)
+ xparam1(estim_params_.transform==1) = transform_to_bounded(xparam1(estim_params_.transform==1),BoundsInfo.lb(estim_params_.transform==1),BoundsInfo.ub(estim_params_.transform==1));
+end
+
% Set flag related to analytical derivatives.
analytic_derivation = options_.analytic_derivation;
@@ -135,7 +140,12 @@ if options_.occbin.likelihood.status
occbin_.info= {options_, dr,endo_steady_state,exo_steady_state,exo_det_steady_state, M_, occbin_options, TTx, RRx, CCx,T0,R0};
else
% Linearize the model around the deterministic steady state and extract the matrices of the state equation (T and R).
- [T,R,SteadyState,info,dr, M_.params] = dynare_resolve(M_,options_,dr, endo_steady_state, exo_steady_state, exo_det_steady_state,'restrict');
+ if isfield(M_,'csn')
+ is_restrict_state_space = false;
+ [T,R,SteadyState,info,dr, M_.params] = dynare_resolve(M_,options_,dr, endo_steady_state, exo_steady_state, exo_det_steady_state);
+ else
+ [T,R,SteadyState,info,dr, M_.params] = dynare_resolve(M_,options_,dr, endo_steady_state, exo_steady_state, exo_det_steady_state,'restrict');
+ end
occbin_.status = false;
end
@@ -626,13 +636,32 @@ if ((kalman_algo==1) || (kalman_algo==3))% Multivariate Kalman Filter
T,Q,R,Z,Zflag,H,diffuse_periods, ...
options_.presample);
else
- [LIK,lik] = kalman_filter(Y,diffuse_periods+1,size(Y,2), ...
- a_0_given_tm1,Pstar, ...
- kalman_tol, riccati_tol, ...
- options_.rescale_prediction_error_covariance, ...
- options_.presample, ...
- T,Q,R,H,Z,mm,pp,rr,Zflag,diffuse_periods, ...
- analytic_deriv_info{:});
+ if all(M_.Skew_e==0) && ~options_.kalman.pskf.use_in_gaussian_case
+ [LIK,lik] = kalman_filter(Y,diffuse_periods+1,size(Y,2), ...
+ a_0_given_tm1,Pstar, ...
+ kalman_tol, riccati_tol, ...
+ options_.rescale_prediction_error_covariance, ...
+ options_.presample, ...
+ T,Q,R,H,Z,mm,pp,rr,Zflag,diffuse_periods, ...
+ analytic_deriv_info{:});
+ else
+ Z = zeros(length(bayestopt_.mf),size(T,1));
+ for i = 1:length(bayestopt_.mf)
+ Z(i,bayestopt_.mf(i)) = 1;
+ end
+ Zflag = 1;
+ Gamma_0 = zeros(size(Pstar)); % initialize at Gaussian distribution
+ nu_0 = zeros(size(a)); % normalization
+ Delta_0 = eye(size(a,1)); % normalization
+ [LIK, lik] = kalman_filter_pruned_skewed( ...
+ Y,diffuse_periods+1,size(Y,2), ...
+ a,Pstar,Gamma_0,nu_0,Delta_0, ...
+ kalman_tol,riccati_tol,options_.rescale_prediction_error_covariance,...
+ options_.presample,diffuse_periods, ...
+ T,R,Z,Zflag, ...
+ M_.csn.mu_e,M_.csn.Sigma_e,M_.csn.Gamma_e,M_.csn.nu_e,M_.csn.Delta_e,H, ...
+ options_.kalman.pskf.prune_tol);
+ end
end
end
else
diff --git a/matlab/estimation/dynare_estimation_1.m b/matlab/estimation/dynare_estimation_1.m
index e82e4389ddcb1e9669c851d9068a8734f33bb00d..8f174bc2d47d1eed09bac2768975dc3d8189d49b 100644
--- a/matlab/estimation/dynare_estimation_1.m
+++ b/matlab/estimation/dynare_estimation_1.m
@@ -154,8 +154,9 @@ nvx = estim_params_.nvx; % Variance of the structural innovations (number of pa
nvn = estim_params_.nvn; % Variance of the measurement innovations (number of parameters).
ncx = estim_params_.ncx; % Covariance of the structural innovations (number of parameters).
ncn = estim_params_.ncn; % Covariance of the measurement innovations (number of parameters).
+nsx = estim_params_.nsx; % Skewness of the structural innovations (number of parameters).
np = estim_params_.np ; % Number of deep parameters.
-nx = nvx+nvn+ncx+ncn+np; % Total number of parameters to be estimated.
+nx = nvx+nvn+ncx+ncn+nsx+np; % Total number of parameters to be estimated.
if ~isempty(estim_params_)
M_ = set_all_parameters(xparam1,estim_params_,M_);
@@ -245,7 +246,15 @@ if ~isequal(options_.mode_compute,0) && ~options_.mh_posterior_mode_estimation &
end
end
end
- [xparam1, fval, ~, hh, options_, Scale, new_rat_hess_info] = dynare_minimize_objective(objective_function,xparam1,current_optimizer,options_,[bounds.lb bounds.ub],bayestopt_.name,bayestopt_,hh,dataset_,dataset_info,options_,M_,estim_params_,bayestopt_,bounds,oo_.dr, oo_.steady_state,oo_.exo_steady_state,oo_.exo_det_steady_state);
+ % transform parameters from bounded to unbounded domain
+ bounds0 = bounds;
+ if any(estim_params_.transform)
+ estim_params_.transformed = true;
+ xparam1(estim_params_.transform==1) = transform_to_unbounded(xparam1(estim_params_.transform==1),bounds0.lb(estim_params_.transform==1),bounds0.ub(estim_params_.transform==1));
+ bounds.lb(estim_params_.transform==1) = -Inf;
+ bounds.ub(estim_params_.transform==1) = Inf;
+ end
+ [xparam1, fval, ~, hh, options_, Scale, new_rat_hess_info] = dynare_minimize_objective(objective_function,xparam1,current_optimizer,options_,[bounds.lb bounds.ub],bayestopt_.name,bayestopt_,hh,dataset_,dataset_info,options_,M_,estim_params_,bayestopt_,bounds0,oo_.dr, oo_.steady_state,oo_.exo_steady_state,oo_.exo_det_steady_state);
fprintf('\nFinal value of minus the log posterior (or likelihood):%f \n', fval);
if isnumeric(current_optimizer)
@@ -261,6 +270,14 @@ if ~isequal(options_.mode_compute,0) && ~options_.mh_posterior_mode_estimation &
bayestopt_.jscale(:) = options_.mh_jscale;
end
end
+
+ % transform parameters back from unbounded to bounded
+ if any(estim_params_.transform)
+ estim_params_.transformed = false;
+ bounds = bounds0;
+ xparam1(estim_params_.transform==1) = transform_to_bounded(xparam1(estim_params_.transform==1),bounds0.lb(estim_params_.transform==1),bounds0.ub(estim_params_.transform==1));
+ end
+
if ~isnumeric(current_optimizer) || ~isequal(current_optimizer,6) %always already computes covariance matrix
if options_.cova_compute == 1 %user did not request covariance not to be computed
if options_.analytic_derivation && strcmp(func2str(objective_function),'dsge_likelihood')
@@ -275,9 +292,17 @@ if ~isequal(options_.mode_compute,0) && ~options_.mh_posterior_mode_estimation &
% the hessian from outer product gradient of optimizer 5 below
if options_.hessian.use_penalized_objective
penalized_objective_function = str2func('penalty_objective_function');
- hh = hessian(penalized_objective_function, xparam1, options_.gstep, objective_function, fval, dataset_, dataset_info, options_, M_, estim_params_, bayestopt_, bounds,oo_.dr,oo_.steady_state,oo_.exo_steady_state,oo_.exo_det_steady_state);
+ if isfield(options_.hessian,'use_onesided') && options_.hessian.use_onesided == 1
+ hh = hessian_with_bounds(penalized_objective_function, xparam1, options_.gstep, bounds, objective_function, fval, dataset_, dataset_info, options_, M_, estim_params_, bayestopt_, bounds,oo_.dr,oo_.steady_state,oo_.exo_steady_state,oo_.exo_det_steady_state);
+ else
+ hh = hessian(penalized_objective_function, xparam1, options_.gstep, objective_function, fval, dataset_, dataset_info, options_, M_, estim_params_, bayestopt_, bounds,oo_.dr,oo_.steady_state,oo_.exo_steady_state,oo_.exo_det_steady_state);
+ end
else
- hh = hessian(objective_function, xparam1, options_.gstep, dataset_, dataset_info, options_, M_, estim_params_, bayestopt_, bounds,oo_.dr,oo_.steady_state,oo_.exo_steady_state,oo_.exo_det_steady_state);
+ if isfield(options_.hessian,'use_onesided') && options_.hessian.use_onesided == 1
+ hh = hessian_with_bounds(objective_function, xparam1, options_.gstep, bounds, dataset_, dataset_info, options_, M_, estim_params_, bayestopt_, bounds,oo_.dr,oo_.steady_state,oo_.exo_steady_state,oo_.exo_det_steady_state);
+ else
+ hh = hessian(objective_function, xparam1, options_.gstep, dataset_, dataset_info, options_, M_, estim_params_, bayestopt_, bounds,oo_.dr,oo_.steady_state,oo_.exo_steady_state,oo_.exo_det_steady_state);
+ end
end
hh = reshape(hh, nx, nx);
elseif isnumeric(current_optimizer) && isequal(current_optimizer,5)
@@ -587,7 +612,9 @@ elseif options_.partial_information ||...
% smoothing not yet supported
else
%% ML estimation, or posterior mode without Metropolis-Hastings or Metropolis without Bayesian smoothed variables
- oo_=save_display_classical_smoother_results(xparam1,M_,oo_,options_,bayestopt_,dataset_,dataset_info,estim_params_);
+ if ~isfield(M_,'csn')
+ oo_=save_display_classical_smoother_results(xparam1,M_,oo_,options_,bayestopt_,dataset_,dataset_info,estim_params_);
+ end
end
if options_.forecast == 0 || options_.mh_replic > 0 || options_.load_mh_file
diff --git a/matlab/estimation/dynare_estimation_init.m b/matlab/estimation/dynare_estimation_init.m
index ffa982c0c29a434a867a655288db6e00f0d217b1..d9dec60e96708eaf2a0a3f66b346ea227aa7e303 100644
--- a/matlab/estimation/dynare_estimation_init.m
+++ b/matlab/estimation/dynare_estimation_init.m
@@ -138,7 +138,7 @@ else
end
% Set priors over the estimated parameters.
-if ~isempty(estim_params_) && ~(isfield(estim_params_,'nvx') && (size(estim_params_.var_exo,1)+size(estim_params_.var_endo,1)+size(estim_params_.corrx,1)+size(estim_params_.corrn,1)+size(estim_params_.param_vals,1))==0)
+if ~isempty(estim_params_) && ~(isfield(estim_params_,'nvx') && (size(estim_params_.skew_exo,1)+size(estim_params_.var_exo,1)+size(estim_params_.var_endo,1)+size(estim_params_.corrx,1)+size(estim_params_.corrn,1)+size(estim_params_.param_vals,1))==0)
[xparam1,estim_params_,bayestopt_,lb,ub,M_] = set_prior(estim_params_,M_,options_);
end
@@ -152,12 +152,12 @@ if exist([M_.fname '_prior_restrictions.m'],'file')
end
% Check that the provided mode_file is compatible with the current estimation settings.
-if ~isempty(estim_params_) && ~(isfield(estim_params_,'nvx') && sum(estim_params_.nvx+estim_params_.nvn+estim_params_.ncx+estim_params_.ncn+estim_params_.np)==0) && ~isempty(options_.mode_file) && ~options_.mh_posterior_mode_estimation
+if ~isempty(estim_params_) && ~(isfield(estim_params_,'nvx') && sum(estim_params_.nvx+estim_params_.nvn+estim_params_.ncx+estim_params_.ncn+estim_params_.nsx+estim_params_.np)==0) && ~isempty(options_.mode_file) && ~options_.mh_posterior_mode_estimation
[xparam1, hh] = check_mode_file(xparam1, hh, options_, bayestopt_);
end
%check for calibrated covariances before updating parameters
-if ~isempty(estim_params_) && ~(isfield(estim_params_,'nvx') && sum(estim_params_.nvx+estim_params_.nvn+estim_params_.ncx+estim_params_.ncn+estim_params_.np)==0)
+if ~isempty(estim_params_) && ~(isfield(estim_params_,'nvx') && sum(estim_params_.nvx+estim_params_.nvn+estim_params_.ncx+estim_params_.ncn+estim_params_.nsx+estim_params_.np)==0)
estim_params_=check_for_calibrated_covariances(estim_params_,M_);
end
@@ -180,7 +180,7 @@ if ~isempty(bayestopt_) && all(bayestopt_.pshape==0) && any(isnan(xparam1))
error('ML estimation requires all estimated parameters to be initialized, either in an estimated_params or estimated_params_init-block ')
end
-if ~isempty(estim_params_) && ~(all(strcmp(fieldnames(estim_params_),'full_calibration_detected')) || (isfield(estim_params_,'nvx') && sum(estim_params_.nvx+estim_params_.nvn+estim_params_.ncx+estim_params_.ncn+estim_params_.np)==0))
+if ~isempty(estim_params_) && ~(all(strcmp(fieldnames(estim_params_),'full_calibration_detected')) || (isfield(estim_params_,'nvx') && sum(estim_params_.nvx+estim_params_.nvn+estim_params_.ncx+estim_params_.ncn+estim_params_.nsx+estim_params_.np)==0))
if ~isempty(bayestopt_) && any(bayestopt_.pshape > 0)
% Plot prior densities.
if ~options_.nograph && options_.plot_priors
@@ -210,7 +210,7 @@ if ~isempty(estim_params_) && ~(all(strcmp(fieldnames(estim_params_),'full_calib
end
end
-if isempty(estim_params_) || all(strcmp(fieldnames(estim_params_),'full_calibration_detected')) || (isfield(estim_params_,'nvx') && sum(estim_params_.nvx+estim_params_.nvn+estim_params_.ncx+estim_params_.ncn+estim_params_.np)==0) % If estim_params_ is empty (e.g. when running the smoother on a calibrated model)
+if isempty(estim_params_) || all(strcmp(fieldnames(estim_params_),'full_calibration_detected')) || (isfield(estim_params_,'nvx') && sum(estim_params_.nvx+estim_params_.nvn+estim_params_.ncx+estim_params_.ncn+estim_params_.nsx+estim_params_.np)==0) % If estim_params_ is empty (e.g. when running the smoother on a calibrated model)
if ~options_.smoother
error('Estimation: the ''estimated_params'' block is mandatory (unless you are running a smoother)')
end
@@ -229,11 +229,13 @@ if isempty(estim_params_) || all(strcmp(fieldnames(estim_params_),'full_calibrat
estim_params_.var_endo=[];
estim_params_.corrx=[];
estim_params_.corrn=[];
+ estim_params_.skew_exo=[];
estim_params_.param_vals=[];
estim_params_.nvx = 0;
estim_params_.nvn = 0;
estim_params_.ncx = 0;
estim_params_.ncn = 0;
+ estim_params_.nsx = 0;
estim_params_.np = 0;
bounds.lb = [];
bounds.ub = [];
@@ -333,7 +335,7 @@ if options_.analytic_derivation
if isfield(options_,'identification_check_endogenous_params_with_no_prior')
M_local.params = M_local.params*1.01; %vary parameters
else
- M_local.params(estim_params_.param_vals(:,1)) = xparam1(estim_params_.nvx+estim_params_.ncx+estim_params_.nvn+estim_params_.ncn+1:end); %set parameters
+ M_local.params(estim_params_.param_vals(:,1)) = xparam1(estim_params_.nvx+estim_params_.ncx+estim_params_.nvn+estim_params_.nsx+estim_params_.ncn+1:end); %set parameters
M_local.params(estim_params_.param_vals(:,1)) = M_local.params(estim_params_.param_vals(:,1))*1.01; %vary parameters
end
if options_.diffuse_filter || options_.steadystate.nocheck
@@ -416,8 +418,9 @@ nvx = estim_params_.nvx;
ncx = estim_params_.ncx;
nvn = estim_params_.nvn;
ncn = estim_params_.ncn;
+nsx = estim_params_.nsx;
if estim_params_.np
- M_local.params(estim_params_.param_vals(:,1)) = xparam1(nvx+ncx+nvn+ncn+1:end);
+ M_local.params(estim_params_.param_vals(:,1)) = xparam1(nvx+ncx+nvn+ncn+nsx+1:end);
end
[oo_.steady_state, params,info] = evaluate_steady_state(oo_.steady_state,[oo_.exo_steady_state; oo_.exo_det_steady_state],M_local,options_,steadystate_check_flag);
diff --git a/matlab/estimation/execute_prior_posterior_function.m b/matlab/estimation/execute_prior_posterior_function.m
index 898cf63005952007353f40404ec49f83ee749613..d39c1f15b88f715b54436656f6f73884d7312cd6 100644
--- a/matlab/estimation/execute_prior_posterior_function.m
+++ b/matlab/estimation/execute_prior_posterior_function.m
@@ -62,7 +62,7 @@ if strcmpi(type,'posterior')
elseif strcmpi(type,'prior')
% Get informations about the prior distribution.
if isempty(bayestopt_)
- if ~isempty(estim_params_) && ~(isfield(estim_params_,'nvx') && (size(estim_params_.var_exo,1)+size(estim_params_.var_endo,1)+size(estim_params_.corrx,1)+size(estim_params_.corrn,1)+size(estim_params_.param_vals,1))==0)
+ if ~isempty(estim_params_) && ~(isfield(estim_params_,'nvx') && (size(estim_params_.var_exo,1)+size(estim_params_.var_endo,1)+size(estim_params_.corrx,1)+size(estim_params_.corrn,1)+size(estim_params_.skew_exo,1)+size(estim_params_.param_vals,1))==0)
[~,estim_params_,bayestopt_,~,~,M_] = set_prior(estim_params_,M_,options_);
else
error('The prior distributions are not properly set up.')
diff --git a/matlab/estimation/fill_mh_mode.m b/matlab/estimation/fill_mh_mode.m
index 0d2581ce108a843c92e368fe0f84750114b733f1..7cd90d5c7bf2f95e45d4dffc40f2b577f717ffc1 100644
--- a/matlab/estimation/fill_mh_mode.m
+++ b/matlab/estimation/fill_mh_mode.m
@@ -37,10 +37,11 @@ nvx = estim_params_.nvx; % Variance of the structural innovations (number of pa
nvn = estim_params_.nvn; % Variance of the measurement innovations (number of parameters).
ncx = estim_params_.ncx; % Covariance of the structural innovations (number of parameters).
ncn = estim_params_.ncn; % Covariance of the measurement innovations (number of parameters).
+nsx = estim_params_.nsx; % Skewness of the structural innovations (numberr of parameters).
np = estim_params_.np ; % Number of deep parameters.
if np
- ip = nvx+nvn+ncx+ncn+1;
+ ip = nvx+nvn+ncx+ncn+nsx+1;
for i=1:np
k = estim_params_.param_vals(i,1);
name = M_.param_names{k};
@@ -68,7 +69,6 @@ if nvn
ip = ip+1;
end
end
-
if ncx
ip = nvx+nvn+1;
for i=1:ncx
@@ -80,7 +80,6 @@ if ncx
ip = ip+1;
end
end
-
if ncn
ip = nvx+nvn+ncx+1;
for i=1:ncn
@@ -92,3 +91,13 @@ if ncn
ip = ip+1;
end
end
+if nsx
+ ip = nvx+nvn+ncx+ncn+1;
+ for i=1:nsx
+ k = estim_params_.skew_exo(i,1);
+ name = M_.exo_names{k};
+ oo_.([field_name '_mode']).shocks_skew.(name)= xparam1(ip);
+ oo_.([field_name '_std_at_mode']).shocks_skew.(name) = stdh(ip);
+ ip = ip+1;
+ end
+end
\ No newline at end of file
diff --git a/matlab/estimation/get_posterior_parameters.m b/matlab/estimation/get_posterior_parameters.m
index ffec01826c5485dbe5c17b7591bfbf418cb09ddf..679d570d8a41bc477d7bee437376439b3c5438cd 100644
--- a/matlab/estimation/get_posterior_parameters.m
+++ b/matlab/estimation/get_posterior_parameters.m
@@ -39,9 +39,10 @@ nvx = estim_params_.nvx;
nvn = estim_params_.nvn;
ncx = estim_params_.ncx;
ncn = estim_params_.ncn;
+nsx = estim_params_.nsx;
np = estim_params_.np;
-xparam = zeros(nvx+nvn+ncx+ncn+np,1);
+xparam = zeros(nvx+nvn+ncx+ncn+nsx+np,1);
m = 1;
for i=1:nvx
@@ -76,6 +77,13 @@ for i=1:ncn
m = m+1;
end
+for i=1:nsx
+ k1 = estim_params_.skew_exo(i,1);
+ name1 = M_.exo_names{k1};
+ xparam(m) = oo_.([field1 type]).shocks_skew.(name1);
+ m = m+1;
+end
+
FirstDeep = m;
for i=1:np
diff --git a/matlab/estimation/initial_estimation_checks.m b/matlab/estimation/initial_estimation_checks.m
index 03596d646e0ae361b0770b4244e8a600e8e9f253..99dc35ab7edb3c1f7aa69f103d895a76c5ca3387 100644
--- a/matlab/estimation/initial_estimation_checks.m
+++ b/matlab/estimation/initial_estimation_checks.m
@@ -166,7 +166,7 @@ end
% check and display warning if negative values of stderr or corr params are outside unit circle for Bayesian estimation
if any(bayestopt_.pshape)
- check_prior_stderr_corr(estim_params_,bayestopt_);
+ check_prior_stderr_corr_skew(estim_params_,bayestopt_);
end
% display warning if some parameters are still NaN
diff --git a/matlab/estimation/prior_posterior_statistics.m b/matlab/estimation/prior_posterior_statistics.m
index a0fb6128cc73a382992d562163771867004122db..9ecc415942419b54bf7e1da1176091d9a8a62d32 100644
--- a/matlab/estimation/prior_posterior_statistics.m
+++ b/matlab/estimation/prior_posterior_statistics.m
@@ -48,7 +48,7 @@ Y = transpose(dataset_.data);
gend = dataset_.nobs;
nvn = estim_params_.nvn;
-npar = estim_params_.nvx+nvn+estim_params_.ncx+estim_params_.ncn+estim_params_.np;
+npar = estim_params_.nvx+nvn+estim_params_.ncx+estim_params_.ncn+nvn+estim_params_.nsx+estim_params_.np;
naK = length(options_.filter_step_ahead);
MaxNumberOfBytes=options_.MaxNumberOfBytes;
diff --git a/matlab/estimation/selec_posterior_draws.m b/matlab/estimation/selec_posterior_draws.m
index 161f8e8e39f284442928c94723c02036e7ebd6c6..cb94d33b4563e14f992fcca3a6c015cae650f93d 100644
--- a/matlab/estimation/selec_posterior_draws.m
+++ b/matlab/estimation/selec_posterior_draws.m
@@ -47,6 +47,7 @@ npar = estim_params_.nvx;
npar = npar + estim_params_.nvn;
npar = npar + estim_params_.ncx;
npar = npar + estim_params_.ncn;
+npar = npar + estim_params_.nsx;
npar = npar + estim_params_.np;
% Select one task:
diff --git a/matlab/estimation/set_prior.m b/matlab/estimation/set_prior.m
index 3069a8a09598ace5870d106aee51f254b80af1d9..5d16342ac3da64054dcc0602c66fd2332957c593 100644
--- a/matlab/estimation/set_prior.m
+++ b/matlab/estimation/set_prior.m
@@ -39,13 +39,17 @@ nvx = size(estim_params_.var_exo,1);
nvn = size(estim_params_.var_endo,1);
ncx = size(estim_params_.corrx,1);
ncn = size(estim_params_.corrn,1);
+nsx = size(estim_params_.skew_exo,1);
np = size(estim_params_.param_vals,1);
estim_params_.nvx = nvx; %exogenous shock variances
estim_params_.nvn = nvn; %endogenous variances, i.e. measurement error
estim_params_.ncx = ncx; %exogenous shock correlations
estim_params_.ncn = ncn; % correlation between endogenous variables, i.e. measurement error.
+estim_params_.nsx = nsx; %exogenous shock skewness
estim_params_.np = np; % other parameters of the model
+estim_params_.transform = []; % index for transformed parameters
+estim_params_.transformed = false; % indicator whether current parameters are transformed (only during optimization)
xparam1 = [];
ub = []; % Upper bound imposed for optimization.
@@ -55,23 +59,24 @@ bayestopt_.p1 = []; % prior mean
bayestopt_.p2 = []; % prior standard deviation
bayestopt_.p3 = []; % lower bound of the distribution, only considering whether a generalized distribution is used, not when the prior is truncated
bayestopt_.p4 = []; % upper bound of the distribution, only considering whether a generalized distribution is used, not when the prior is truncated
-bayestopt_.p5 = zeros(nvx+nvn+ncx+ncn+np,1); % prior mode
+bayestopt_.p5 = zeros(nvx+nvn+ncx+ncn+nsx+np,1); % prior mode
bayestopt_.p6 = []; % first hyper-parameter (\alpha for the BETA and GAMMA distributions, s for the INVERSE GAMMAs, expectation for the GAUSSIAN distribution, lower bound for the UNIFORM distribution).
bayestopt_.p7 = []; % second hyper-parameter (\beta for the BETA and GAMMA distributions, \nu for the INVERSE GAMMAs, standard deviation for the GAUSSIAN distribution, upper bound for the UNIFORM distribution).
bayestopt_.jscale = []; %jscale might subsequently be overwritten by mode_compute=6 or check_posterior_sampler_options
bayestopt_.name = {};
if nvx
- xparam1 = estim_params_.var_exo(:,2);
- ub = estim_params_.var_exo(:,4);
- lb = estim_params_.var_exo(:,3);
- bayestopt_.pshape = estim_params_.var_exo(:,5);
- bayestopt_.p1 = estim_params_.var_exo(:,6);
- bayestopt_.p2 = estim_params_.var_exo(:,7);
- bayestopt_.p3 = estim_params_.var_exo(:,8); %take generalized distribution into account
- bayestopt_.p4 = estim_params_.var_exo(:,9); %take generalized distribution into account
- bayestopt_.jscale = estim_params_.var_exo(:,10);
- bayestopt_.name = M_.exo_names(estim_params_.var_exo(:,1));
+ xparam1 = [xparam1; estim_params_.var_exo(:,2)];
+ ub = [ub; estim_params_.var_exo(:,4)];
+ lb = [lb; estim_params_.var_exo(:,3)];
+ estim_params_.transform = [estim_params_.transform; estim_params_.var_exo(:,11)];
+ bayestopt_.pshape = [ bayestopt_.pshape; estim_params_.var_exo(:,5)];
+ bayestopt_.p1 = [ bayestopt_.p1; estim_params_.var_exo(:,6) ];
+ bayestopt_.p2 = [ bayestopt_.p2; estim_params_.var_exo(:,7) ];
+ bayestopt_.p3 = [ bayestopt_.p3; estim_params_.var_exo(:,8) ]; %take generalized distribution into account
+ bayestopt_.p4 = [ bayestopt_.p4; estim_params_.var_exo(:,9) ]; %take generalized distribution into account
+ bayestopt_.jscale = [ bayestopt_.jscale; estim_params_.var_exo(:,10) ];
+ bayestopt_.name = [ bayestopt_.name; M_.exo_names(estim_params_.var_exo(:,1)) ];
end
if nvn
estim_params_.nvn_observable_correspondence=NaN(nvn,1); % stores number of corresponding observable
@@ -90,6 +95,7 @@ if nvn
xparam1 = [xparam1; estim_params_.var_endo(:,2)];
ub = [ub; estim_params_.var_endo(:,4)];
lb = [lb; estim_params_.var_endo(:,3)];
+ estim_params_.transform = [estim_params_.transform; estim_params_.var_endo(:,11)];
bayestopt_.pshape = [ bayestopt_.pshape; estim_params_.var_endo(:,5)];
bayestopt_.p1 = [ bayestopt_.p1; estim_params_.var_endo(:,6)];
bayestopt_.p2 = [ bayestopt_.p2; estim_params_.var_endo(:,7)];
@@ -102,6 +108,7 @@ if ncx
xparam1 = [xparam1; estim_params_.corrx(:,3)];
ub = [ub; max(min(estim_params_.corrx(:,5),1),-1)];
lb = [lb; min(max(estim_params_.corrx(:,4),-1),1)];
+ estim_params_.transform = [estim_params_.transform; estim_params_.corrx(:,12)];
bayestopt_.pshape = [ bayestopt_.pshape; estim_params_.corrx(:,6)];
bayestopt_.p1 = [ bayestopt_.p1; estim_params_.corrx(:,7)];
bayestopt_.p2 = [ bayestopt_.p2; estim_params_.corrx(:,8)];
@@ -122,10 +129,12 @@ if ncn
nvarobs = length(options_.varobs);
M_.H = zeros(nvarobs,nvarobs);
M_.Correlation_matrix_ME = eye(nvarobs);
+ M_.Sigma_eps = M_.H;
end
xparam1 = [xparam1; estim_params_.corrn(:,3)];
ub = [ub; max(min(estim_params_.corrn(:,5),1),-1)];
lb = [lb; min(max(estim_params_.corrn(:,4),-1),1)];
+ estim_params_.transform = [estim_params_.transform; estim_params_.corrn(:,12)];
bayestopt_.pshape = [ bayestopt_.pshape; estim_params_.corrn(:,6)];
bayestopt_.p1 = [ bayestopt_.p1; estim_params_.corrn(:,7)];
bayestopt_.p2 = [ bayestopt_.p2; estim_params_.corrn(:,8)];
@@ -145,10 +154,28 @@ if ncn
estim_params_.corrn_observable_correspondence(i,:)=[obsi1, obsi2];
end
end
+if nsx
+ xparam1 = [xparam1; estim_params_.skew_exo(:,2)];
+ ub = [ub; estim_params_.skew_exo(:,4)];
+ lb = [lb; estim_params_.skew_exo(:,3)];
+ estim_params_.transform = [estim_params_.transform; estim_params_.skew_exo(:,11)];
+ bayestopt_.pshape = [ bayestopt_.pshape; estim_params_.skew_exo(:,5)];
+ bayestopt_.p1 = [ bayestopt_.p1; estim_params_.skew_exo(:,6) ];
+ bayestopt_.p2 = [ bayestopt_.p2; estim_params_.skew_exo(:,7) ];
+ bayestopt_.p3 = [ bayestopt_.p3; estim_params_.skew_exo(:,8) ]; %take generalized distribution into account
+ bayestopt_.p4 = [ bayestopt_.p4; estim_params_.skew_exo(:,9) ]; %take generalized distribution into account
+ bayestopt_.jscale = [ bayestopt_.jscale; estim_params_.skew_exo(:,10) ];
+ baseid = length(bayestopt_.name);
+ bayestopt_.name = [bayestopt_.name; cell(nsx, 1)];
+ for i = 1:nsx
+ bayestopt_.name(baseid+i) = {sprintf('skew %s', M_.exo_names{estim_params_.skew_exo(i,1)})};
+ end
+end
if np
xparam1 = [xparam1; estim_params_.param_vals(:,2)];
ub = [ub; estim_params_.param_vals(:,4)];
lb = [lb; estim_params_.param_vals(:,3)];
+ estim_params_.transform = [estim_params_.transform; estim_params_.param_vals(:,11)];
bayestopt_.pshape = [ bayestopt_.pshape; estim_params_.param_vals(:,5)];
bayestopt_.p1 = [ bayestopt_.p1; estim_params_.param_vals(:,6)];
bayestopt_.p2 = [ bayestopt_.p2; estim_params_.param_vals(:,7)];
diff --git a/matlab/get_all_parameters.m b/matlab/get_all_parameters.m
index 8c7614a81d1b9ebf63fe9f158093f90baa740bf1..5a95e19e5d6ba675308d64195b77dc2d027098d7 100644
--- a/matlab/get_all_parameters.m
+++ b/matlab/get_all_parameters.m
@@ -37,20 +37,24 @@ if ~isempty(estim_params_)
ncx = estim_params_.ncx;
nvn = estim_params_.nvn;
ncn = estim_params_.ncn;
+ nsx = estim_params_.nsx;
np = estim_params_.np;
else
nvx = 0;
ncx = 0;
nvn = 0;
ncn = 0;
+ nsx = 0;
np = 0;
end
Sigma_e = M_.Sigma_e;
Correlation_matrix = M_.Correlation_matrix;
H = M_.H;
Correlation_matrix_ME = M_.Correlation_matrix_ME;
+Skew_e = M_.Skew_e;
+
+xparam1=NaN(nvx+nvn+ncx+ncn+nsx+np,1);
-xparam1=NaN(nvx+ncx+nvn+ncn+np,1);
% stderrs of the exogenous shocks
if nvx
var_exo = estim_params_.var_exo;
@@ -59,21 +63,18 @@ if nvx
xparam1(i)=sqrt(Sigma_e(k,k));
end
end
-% update offset
-offset = nvx;
+offset = nvx; % update offset
-% setting measument error variance
+% stderr of measurement errors
if nvn
for i=1:nvn
k = estim_params_.nvn_observable_correspondence(i,1);
xparam1(offset+i)=sqrt(H(k,k));
end
end
+offset = nvx+nvn; % update offset
-% update offset
-offset = nvx+nvn;
-
-% correlations among shocks (ncx)
+% correlations among shocks
if ncx
corrx = estim_params_.corrx;
for i=1:ncx
@@ -82,9 +83,9 @@ if ncx
xparam1(i+offset)=Correlation_matrix(k1,k2);
end
end
-% update offset
-offset = nvx+nvn+ncx;
+offset = nvx+nvn+ncx; % update offset
+% correlations among measurement errors
if ncn
corrn_observable_correspondence = estim_params_.corrn_observable_correspondence;
for i=1:ncn
@@ -93,10 +94,17 @@ if ncn
xparam1(i+offset)=Correlation_matrix_ME(k1,k2);
end
end
+offset = nvx+nvn+ncx+ncn; % update offset
+if nsx
+ skew_exo = estim_params_.skew_exo;
+ for i=1:nsx
+ k = skew_exo(i,1);
+ xparam1(i+offset)=Skew_e(k);
+ end
+end
% update offset
-offset = nvx+ncx+nvn+ncn;
-
+offset = nvx+nvn+ncx+ncn+nsx;
% structural parameters
if np
diff --git a/matlab/get_error_message.m b/matlab/get_error_message.m
index 7156e3c253c4e7d2dc1b60b35d92491e4d016227..4c5b97a9818d5fa4e9193fd011376e9c8cbee107 100644
--- a/matlab/get_error_message.m
+++ b/matlab/get_error_message.m
@@ -113,6 +113,10 @@ switch info(1)
message = 'You are estimating a DSGE-VAR model, but the implied covariance matrix of the VAR''s innovations, based on the artificial sample, is not positive definite!';
case 55
message = 'Fast Kalman filter only works with stationary models [lik_init=1] or stationary observables for non-stationary models [lik_init=3]';
+ case 56
+ message = 'The skewness transformation failed.';
+ case 57
+ message = 'The skewness coefficient of structural shocks is beyond its theoretical bound.';
case 61 %Discretionary policy
message = 'Discretionary policy: maximum number of iterations has been reached. Procedure failed.';
case 62
diff --git a/matlab/get_the_name.m b/matlab/get_the_name.m
index a373699a7f8ee96cc5fca7fd54d91df102af4a70..b8826e653c03a8ea40972b1f7e23353bd88575ad 100644
--- a/matlab/get_the_name.m
+++ b/matlab/get_the_name.m
@@ -43,6 +43,7 @@ nvx = estim_params_.nvx;
nvn = estim_params_.nvn;
ncx = estim_params_.ncx;
ncn = estim_params_.ncn;
+nsx = estim_params_.nsx;
if k <= nvx
vname = M_.exo_names{estim_params_.var_exo(k,1)};
@@ -78,8 +79,16 @@ elseif k <= (nvx+nvn+ncx+ncn)
tname = sprintf('%s,%s', M_.endo_names_tex{k1}, M_.endo_names_tex{k2});
texnam = sprintf('$ \\rho^{ME}_{%s} $', tname);
end
-else
+elseif k <= (nvx+nvn+ncx+ncn+nsx)
jj = k - (nvx+nvn+ncx+ncn);
+ vname = M_.exo_names{estim_params_.skew_exo(jj,1)};
+ nam = sprintf('SKEW_%s', vname);
+ if TeX
+ tname = M_.exo_names_tex{estim_params_.skew_exo(jj,1)};
+ texnam = sprintf('$ skew_{%s} $', tname);
+ end
+else
+ jj = k - (nvx+nvn+ncx+ncn+nsx);
jj1 = estim_params_.param_vals(jj,1);
nam = M_.param_names{jj1};
if TeX
diff --git a/matlab/hessian_with_bounds.m b/matlab/hessian_with_bounds.m
new file mode 100644
index 0000000000000000000000000000000000000000..5ae817046163804fa93c79d7285f771e7d57fabd
--- /dev/null
+++ b/matlab/hessian_with_bounds.m
@@ -0,0 +1,144 @@
+% This file is part of the replication files for the paper "Pruned Skewed
+% Kalman Filter and Smoother: Pruned Skewed Kalman Filter and Smoother:
+% With Applications to the Yield Curve and Asymmetric Monetary Policy Shocks"
+% by Gaygysyz Guljanov, Willi Mutschler, Mark Trede
+
+function hessian_mat = hessian_with_bounds(func, x, gstep, bounds, varargin)
+% Computes second order partial derivatives using two-sided finite
+% difference method.
+% based on get_hessian.m from Dynare with the following changes:
+% - if a step up or down of a parameter passes the lower or upper bounds,
+% then the partial derivatives are computed using one-sided finite
+% differences.
+%
+% INPUTS
+% func [string] name of the function
+% x [double] vector, the Hessian of "func" is evaluated at x.
+% gstep [double] scalar, size of epsilon.
+% varargin [void] list of additional arguments for "func".
+%
+% OUTPUTS
+% hessian_mat [double] Hessian matrix
+%
+% ALGORITHM
+% Uses Abramowitz and Stegun (1965) formulas 25.3.23
+% \[
+% \frac{\partial^2 f_{0,0}}{\partial {x^2}} = \frac{1}{h^2}\left( f_{1,0} - 2f_{0,0} + f_{ - 1,0} \right)
+% \]
+% and 25.3.27 p. 884
+%
+% \[
+% \frac{\partial ^2f_{0,0}}{\partial x\partial y} = \frac{-1}{2h^2}\left(f_{1,0} + f_{-1,0} + f_{0,1} + f_{0,-1} - 2f_{0,0} - f_{1,1} - f_{-1,-1} \right)
+% \]
+%
+% SPECIAL REQUIREMENTS
+% none
+%
+
+% Copyright (C) 2001-2017 Dynare Team
+% Copyright (C) 2023 Willi Mutschler
+%
+% This 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.
+%
+% It 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.
+verbose = false;
+
+if ~isa(func, 'function_handle')
+ func = str2func(func);
+end
+ub = bounds.ub;
+lb = bounds.lb;
+n = size(x,1);
+h1 = max(abs(x), sqrt(gstep(1))*ones(n, 1))*eps^(1/6)*gstep(2);
+h_1 = h1;
+xh1 = x+h1;
+h1 = xh1-x;
+xh1 = x-h_1;
+h_1 = x-xh1;
+xh1 = x;
+f0 = feval(func, x, varargin{:});
+f1 = zeros(size(f0, 1), n);
+f_1 = f1;
+
+for i=1:n
+ %do step up
+ xh1(i) = x(i)+h1(i);
+ if xh1(i) <= ub(i)
+ f1(:,i) = feval(func, xh1, varargin{:});
+ else
+ if verbose
+ fprintf(' - don''t do step up for parameter %d\n',i);
+ end
+ f1(:,i) = f0;
+ end
+ %do step down
+ xh1(i) = x(i)-h_1(i);
+ if xh1(i) >= lb(i)
+ f_1(:,i) = feval(func, xh1, varargin{:});
+ else
+ if verbose
+ fprintf(' - don''t do step down for parameter %d\n',i);
+ end
+ f_1(:,i) = f0;
+ end
+ %reset parameter
+ xh1(i) = x(i);
+end
+
+xh_1 = xh1;
+temp = f1+f_1-f0*ones(1, n); %term f_(1,0)+f_(-1,0)-f_(0,0) used later
+
+hessian_mat = zeros(size(f0,1), n*n);
+
+for i=1:n
+ if i > 1
+ %fill symmetric part of Hessian based on previously computed results
+ k = [i:n:n*(i-1)];
+ hessian_mat(:,(i-1)*n+1:(i-1)*n+i-1) = hessian_mat(:,k);
+ end
+ hessian_mat(:,(i-1)*n+i) = (f1(:,i)+f_1(:,i)-2*f0)./(h1(i)*h_1(i)); %formula 25.3.23
+ for j=i+1:n
+ %step in up direction
+ xh1(i) = x(i)+h1(i);
+ if xh1(i) > ub(i)
+ xh1(i) = x(i);
+ if verbose
+ fprintf(' - don''t do cross step up for parameter %d\n',i);
+ end
+ end
+ xh1(j) = x(j)+h_1(j);
+ if xh1(j) > ub(j)
+ xh1(j) = x(j);
+ if verbose
+ fprintf(' - don''t do cross step up for parameter %d\n',j);
+ end
+ end
+ %step in down direction
+ xh_1(i) = x(i)-h1(i);
+ if xh_1(i) < lb(i)
+ xh_1(i) = x(i);
+ if verbose
+ fprintf(' - don''t do cross step down for parameter %d\n',i);
+ end
+ end
+ xh_1(j) = x(j)-h_1(j);
+ if xh_1(j) < lb(j)
+ xh_1(j) = x(j);
+ if verbose
+ fprintf(' - don''t do cross step down for parameter %d\n',j);
+ end
+ end
+ hessian_mat(:,(i-1)*n+j) =-(-feval(func, xh1, varargin{:})-feval(func, xh_1, varargin{:})+temp(:,i)+temp(:,j))./(2*h1(i)*h_1(j)); %formula 25.3.27
+ %reset grid points
+ xh1(i) = x(i);
+ xh1(j) = x(j);
+ xh_1(i) = x(i);
+ xh_1(j) = x(j);
+ end
+end
\ No newline at end of file
diff --git a/matlab/kalman/csnPruneParams.m b/matlab/kalman/csnPruneParams.m
new file mode 100644
index 0000000000000000000000000000000000000000..c0a48176335bc8cca365baa33758f77c78cac928
--- /dev/null
+++ b/matlab/kalman/csnPruneParams.m
@@ -0,0 +1,91 @@
+function [Sigma, Gamma, nu, Delta] = csnPruneParams(Sigma, Gamma, nu, Delta, tol, debug_dim)
+% [Sigma, Gamma, nu, Delta] = csnPruneParams(Sigma, Gamma, nu, Delta, tol, debug_dim)
+% -------------------------------------------------------------------------
+% Reduces the dimension of a CSN distributed random variable
+% Idea:
+% 1) Representation of CSN distribution:
+% X ~ CSN(mu,Sigma,Gamma,nu,Delta) is equivalent to X = W|Z>=0 where
+% [W;Z] ~ N([mu;-nu], [Sigma, Sigma*Gamma'; Gamma*Sigma, Delta+Gamma*Sigma*Gamma']
+% 2) Correlation introduces skewness:
+% Skewness in CSN is based on the correlation between W and Z:
+% - CSN and N are the same distribution if Gamma=0
+% - if correlation is small, then CSN very close to N
+% 3) Prune dimensions if correlation in absolute value is below tol.
+% This reduces the overall skewness dimension to qq < q.
+% -------------------------------------------------------------------------
+% INPUTS
+% - Sigma [p by p] 2nd parameter of the CSN distribution (scale)
+% - Gamma [q by p] 3rd parameter of the CSN distribution (regulates skewness from Gaussian (Gamma=0) to half normal)
+% - nu [q by 1] 4th parameter of the CSN distribution (enables closure of CSN distribution under conditioning)
+% - Delta [q by q] 5th parameter of the CSN distribution (enables closure of CSN distribution und marginalization)
+% - tol [double] threshold value for correlation (in absolute value) below which correlated variables are pruned
+% - debug_dim [boolean] optional notification if there are more than debug_dim dimensions left unpruned
+% where p is the normal dimension and q the skewness dimension
+% -------------------------------------------------------------------------
+% OUTPUTS
+% - Sigma [p by p] 2nd parameter of the pruned CSN distribution (scale)
+% - Gamma [qq by p] 3rd parameter of the pruned CSN distribution (regulates skewness from Gaussian (Gamma=0) to half normal)
+% - nu [qq by 1] 4th parameter of the pruned CSN distribution (enables closure of CSN distribution under conditioning)
+% - Delta [qq by q] 5th parameter of the pruned CSN distribution (enables closure of CSN distribution und marginalization)
+% where qq < q is the skewness dimension of the pruned distribution
+% =========================================================================
+% Copyright (C) 2022-2023 Gaygysyz Guljanov, Willi Mutschler
+%
+% This 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.
+%
+% This 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.
+% -------------------------------------------------------------------------
+% This file is part of the replication files for the paper "Pruned Skewed
+% Kalman Filter and Smoother: Pruned Skewed Kalman Filter and Smoother:
+% With Applications to the Yield Curve and Asymmetric Monetary Policy Shocks"
+% by Gaygysyz Guljanov, Willi Mutschler, Mark Trede
+% =========================================================================
+if nargin < 6
+ debug_dim = [];
+end
+
+p = length(Sigma); % normal dimension
+
+% create correlation matrix of conditional definition of CSN distributed variable
+P = [Sigma, Sigma*Gamma'; Gamma*Sigma, Delta+Gamma*Sigma*Gamma']; P = 0.5*(P + P'); % covariance matrix
+n = length(P);
+try
+ % create correlation matrix from covariance matrix
+ normalization = diag(1./sqrt(diag(P)));
+ R = abs(normalization*P*normalization);
+catch
+ % error handling
+ Sigma = 0;
+ Gamma = 0;
+ nu = 0;
+ Delta = 0;
+ return
+end
+
+% prune dimensions in R if they are in absolute value lower than tol
+R = R-diag(repelem(Inf, n));
+logi2 = max(R(p+1:end, 1:p), [], 2);
+logi2 = (logi2 < tol);
+
+% prune dimensions in CSN parameters
+Gamma(logi2, :) = [];
+nu(logi2) = [];
+Delta(logi2, :) = [];
+Delta(:, logi2) = [];
+
+% Notify if there are more than ... dimensions left unpruned
+if ~isempty(debug_dim)
+ if sum(logi2 > tol) > debug_dim
+ disp('----');
+ disp(max(logi2)); disp(sum(logi2 > tol));
+ disp('----');
+ end
+end
+
+end % csnPruneParams
\ No newline at end of file
diff --git a/matlab/kalman/csn_statespace_linear_transform.m b/matlab/kalman/csn_statespace_linear_transform.m
new file mode 100644
index 0000000000000000000000000000000000000000..2726d7421c34e6f3a543839257f29c308eb8f589
--- /dev/null
+++ b/matlab/kalman/csn_statespace_linear_transform.m
@@ -0,0 +1,66 @@
+function [mu_tr, Sigma_tr, Gamma_tr, nu_tr, Delta_tr] = csn_statespace_linear_transform(G_bar, mu_bar, Sigma_bar, Gamma_bar, nu_bar, Delta_bar)
+% function [mu_tr, Sigma_tr, Gamma_tr, nu_tr, Delta_tr] = csn_statespace_linear_transform(G_bar, mu_bar, Sigma_bar, Gamma_bar, nu_bar, Delta_bar)
+% -------------------------------------------------------------------------
+% Compute CSN parameters for joint distribution of [x_t', eta_t']' from
+% state transition equation: x(t) = G*x(t-1) + R*eta(t)
+% 1) make proper csn random variable first (in case of singularity)
+% 2) do either rank deficient or full rank linear transformation
+% -----------------------------------------------------------------------
+% INPUTS
+% - G_bar = [G, R] multiplying the joint distribution of x_t and eta_t
+% - mu_bar = [mu_t_t; mu_eta]
+% - nu_bar = [nu_t_t; nu_eta];
+% - Sigma_bar = blkdiag_two(Sigma_t_t, Sigma_eta);
+% - Gamma_bar = blkdiag_two(Gamma_t_t, Gamma_eta);
+% - Delta_bar = blkdiag_two(Delta_t_t, Delta_eta);
+% -----------------------------------------------------------------------
+% OUTPUT
+% transformed parameters of CSN joint distribution of x_t and eta_t:
+% mu_tr, Sigma_tr, Gamma_tr, nu_tr, Delta_tr
+% =========================================================================
+% Copyright (C) 2023 Gaygysyz Guljanov
+%
+% This 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.
+%
+% This 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.
+% -------------------------------------------------------------------------
+% This file is part of the replication files for the paper "Pruned Skewed
+% Kalman Filter and Smoother: Pruned Skewed Kalman Filter and Smoother:
+% With Applications to the Yield Curve and Asymmetric Monetary Policy Shocks"
+% by Gaygysyz Guljanov, Willi Mutschler, Mark Trede
+% =========================================================================
+tol = 1e-8;
+
+mu_tr = G_bar * mu_bar;
+Sigma_tr = G_bar * Sigma_bar * G_bar';
+mu_bar = zeros(size(mu_bar)); % make mu_bar a zero vector
+
+%% singular skew-normal to closed skew-normal
+if rcond(Sigma_bar) < 1e-10
+ % Sigma_bar is rank deficient
+ [S_hat, ~, Sigma_bar, Gamma_bar] = singular_skew_normal_2_closed_skew_normal(mu_bar, Sigma_bar, Gamma_bar, tol, 0);
+ G_bar = G_bar * S_hat;
+end
+
+%% Linear Transformation
+% Check if A_mat is rank deficient
+dimensions = size(G_bar);
+if dimensions(1) > dimensions(2)
+ is_singular = rcond(G_bar' * G_bar) < 1e-8;
+else
+ is_singular = rcond(Sigma_tr) < 1e-8;
+end
+% If not full rank, use rank deficient linear transformation
+if is_singular
+ [Gamma_tr, nu_tr, Delta_tr] = csn_statespace_linear_transform_rank_deficient(G_bar, Sigma_bar, Gamma_bar, nu_bar, Delta_bar, tol);
+ return
+end
+[Gamma_tr, nu_tr, Delta_tr] = csn_statespace_linear_transform_full_rank(G_bar, Sigma_bar, Sigma_tr, Gamma_bar, nu_bar, Delta_bar);
+
+end % csn_statespace_linear_transform
\ No newline at end of file
diff --git a/matlab/kalman/csn_statespace_linear_transform_full_rank.m b/matlab/kalman/csn_statespace_linear_transform_full_rank.m
new file mode 100644
index 0000000000000000000000000000000000000000..0fda18103b668a7a9fc844a450502521dc157e92
--- /dev/null
+++ b/matlab/kalman/csn_statespace_linear_transform_full_rank.m
@@ -0,0 +1,57 @@
+function [Gamma_tr, nu_tr, Delta_tr] = csn_statespace_linear_transform_full_rank(G_bar, Sigma_bar, Sigma_tr, Gamma_bar, nu_bar, Delta_bar)
+% function [Gamma_tr, nu_tr, Delta_tr] = csn_statespace_linear_transform_full_rank(G_bar, Sigma_bar, Sigma_tr, Gamma_bar, nu_bar, Delta_bar)
+% -------------------------------------------------------------------------
+% Compute CSN parameters for joint distribution of [x_t', eta_t']' from
+% state transition equation: x(t) = G*x(t-1) + R*eta(t)
+% helper function for full rank linear transformation of CSN distributed random variable
+% -----------------------------------------------------------------------
+% INPUTS
+% - G_bar = [G, R] multiplying the joint distribution of x_t and eta_t
+% - Sigma_bar = blkdiag_two(Sigma_t_t, Sigma_eta);
+% - Sigma_tr : transformed (rank deficient) parameter matrix
+% - Gamma_bar = blkdiag_two(Gamma_t_t, Gamma_eta);
+% - nu_bar = [nu_t_t; nu_eta];
+% - Delta_bar = blkdiag_two(Delta_t_t, Delta_eta);
+% -----------------------------------------------------------------------
+% OUTPUT
+% transformed parameters of rank-deficient CSN joint distribution of x_t and eta_t:
+% Gamma_tr, nu_tr, Delta_tr (returns only skewness parameters, as other parameters are easy to obtain)
+% =========================================================================
+% Copyright (C) 2023 Gaygysyz Guljanov
+%
+% This 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.
+%
+% This 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.
+% -------------------------------------------------------------------------
+% This file is part of the replication files for the paper "Pruned Skewed
+% Kalman Filter and Smoother: Pruned Skewed Kalman Filter and Smoother:
+% With Applications to the Yield Curve and Asymmetric Monetary Policy Shocks"
+% by Gaygysyz Guljanov, Willi Mutschler, Mark Trede
+% =========================================================================
+[nrow, ncol] = size(G_bar);
+
+%% more rows
+if nrow >= ncol
+ nu_tr = nu_bar;
+ Delta_tr = Delta_bar;
+ if nrow == ncol
+ Gamma_tr = Gamma_bar / G_bar;
+ return
+ end
+ Gamma_tr = Gamma_bar / (G_bar' * G_bar) * G_bar';
+ return
+end
+
+%% more columns
+Tmp_mat = Sigma_bar * G_bar' / Sigma_tr;
+Gamma_tr = Gamma_bar * Tmp_mat;
+nu_tr = nu_bar;
+Delta_tr = Delta_bar + Gamma_bar * Sigma_bar * Gamma_bar' - Gamma_bar * Tmp_mat * G_bar * Sigma_bar * Gamma_bar';
+
+end
\ No newline at end of file
diff --git a/matlab/kalman/csn_statespace_linear_transform_rank_deficient.m b/matlab/kalman/csn_statespace_linear_transform_rank_deficient.m
new file mode 100644
index 0000000000000000000000000000000000000000..6fda2a197417d1a33accab38225d36f7e8669f2e
--- /dev/null
+++ b/matlab/kalman/csn_statespace_linear_transform_rank_deficient.m
@@ -0,0 +1,60 @@
+function [Gamma_tr, nu_tr, Delta_tr] = csn_statespace_linear_transform_rank_deficient(G_bar, Sigma_bar, Gamma_bar, nu_bar, Delta_bar, tol)
+% function [Gamma_tr, nu_tr, Delta_tr] = csn_statespace_linear_transform(G_bar, mu_bar, Sigma_bar, Gamma_bar, nu_bar, Delta_bar)
+% -------------------------------------------------------------------------
+% Compute CSN parameters for joint distribution of [x_t', eta_t']' from
+% state transition equation: x(t) = G*x(t-1) + R*eta(t)
+% helper function for rank deficient linear transformation of CSN distributed random variable
+% -----------------------------------------------------------------------
+% INPUTS
+% - G_bar = [G, R] multiplying the joint distribution of x_t and eta_t
+% - Sigma_bar = blkdiag_two(Sigma_t_t, Sigma_eta);
+% - Gamma_bar = blkdiag_two(Gamma_t_t, Gamma_eta);
+% - nu_bar = [nu_t_t; nu_eta];
+% - Delta_bar = blkdiag_two(Delta_t_t, Delta_eta);
+% - tol: tolerance level showing which values should be assumed to be numerically zero
+% -----------------------------------------------------------------------
+% OUTPUT
+% transformed parameters of rank-deficient CSN joint distribution of x_t and eta_t:
+% Gamma_tr, nu_tr, Delta_tr (returns only skewness parameters, as other parameters are easy to obtain)
+% =========================================================================
+% Copyright (C) 2023 Gaygysyz Guljanov
+%
+% This 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.
+%
+% This 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.
+% -------------------------------------------------------------------------
+% This file is part of the replication files for the paper "Pruned Skewed
+% Kalman Filter and Smoother: Pruned Skewed Kalman Filter and Smoother:
+% With Applications to the Yield Curve and Asymmetric Monetary Policy Shocks"
+% by Gaygysyz Guljanov, Willi Mutschler, Mark Trede
+% =========================================================================
+[nrow, ncol] = size(G_bar);
+difference = abs(nrow - ncol);
+[S_mat, Lambda_mat, T_mat] = svd(G_bar); % apply singular value decomposition
+if nrow >= ncol % more rows
+ S_mat = S_mat(:, 1:end-difference);
+ Lambda_mat = Lambda_mat(1:end-difference, :);
+else
+ T_mat = T_mat(:, 1:end-difference);
+ Lambda_mat = Lambda_mat(:, 1:end-difference);
+end
+hold_vec = abs(diag(Lambda_mat)) > tol; % which dimensions to delete
+S_mat = S_mat(:, hold_vec); % delete respective columns
+Lambda_mat = Lambda_mat(hold_vec, hold_vec); % delete respective rows and columns
+T_mat = T_mat(:, hold_vec); % delete respective columns
+
+% Apply the final formulas of theorem
+Tmp_mat = T_mat' * Sigma_bar * T_mat;
+Tmp_mat2 = Gamma_bar * Sigma_bar;
+Tmp_mat3 = Tmp_mat2 * T_mat / Tmp_mat;
+Gamma_tr = (Tmp_mat3 / Lambda_mat) / (S_mat' * S_mat) * S_mat';
+nu_tr = nu_bar;
+Delta_tr = Delta_bar + Tmp_mat2 * Gamma_bar' - Tmp_mat3 * T_mat' * Sigma_bar * Gamma_bar';
+
+end % csn_statespace_linear_transform_rank_deficient
\ No newline at end of file
diff --git a/matlab/kalman/likelihood/kalman_filter_pruned_skewed.m b/matlab/kalman/likelihood/kalman_filter_pruned_skewed.m
new file mode 100644
index 0000000000000000000000000000000000000000..cb9d909fd21f1880eddb6aa5bbe3d8f6760e8938
--- /dev/null
+++ b/matlab/kalman/likelihood/kalman_filter_pruned_skewed.m
@@ -0,0 +1,457 @@
+function [LIK, LIKK] = kalman_filter_pruned_skewed(Y,start,last,mu_tm1_tm1,Sigma_tm1_tm1,Gamma_tm1_tm1,nu_tm1_tm1,Delta_tm1_tm1,kalman_tol,riccati_tol,rescale_prediction_error_covariance,presample,diffuse_periods,G,R,F,Fflag,mu_eta,Sigma_eta,Gamma_eta,nu_eta,Delta_eta,Sigma_eps,prune_tol,cdfmvna_fct)
+% -------------------------------------------------------------------------
+% Evaluate (1) log-likelihood value and (2) filtered states
+% of linear state space model with csn distributed innovations and normally distributed noise:
+% x(t) = G*x(t-1) + R*eta(t) [state transition equation]
+% y(t) = F*x(t) + eps(t) [observation equation]
+% eta(t) ~ CSN(mu_eta,Sigma_eta,Gamma_eta,nu_eta,Delta_eta) [innovations, shocks]
+% eps(t) ~ N(0,Sigma_eps) [noise, measurement error]
+% Dimensions:
+% x(t) is (x_nbr by 1) state vector
+% y(t) is (y_nbr by 1) control vector, i.e. observable variables
+% eta(t) is (eta_nbr by 1) vector of innovations
+% eps(t) is (y_nbr by 1) vector of noise (measurement errors)
+% Note that mu_eta is such that E[eta(t)] = 0
+% -------------------------------------------------------------------------
+% INPUTS
+% - Y [y_nbr by obs_nbr] matrix with data
+% - start [integer scalar] first period
+% - last [integer scalar] last period, (last-first) has to be inferior to obs_nbr
+% - mu_tm1_tm1 [x_nbr by 1] initial filtered value of location parameter of CSN distributed states x (does not equal expectation vector unless Gamma_tm1_tm1=0), i.e. mu_0_0
+% - Sigma_tm1_tm1 [x_nbr by x_nbr] initial filtered value of scale parameter of CSN distributed states x (does not equal covariance matrix unless Gamma_tm1_tm1=0), i.e. Sigma_0_0
+% - Gamma_tm1_tm1 [skewx_dim by x_nbr] initial filtered value of skewness size parameter of CSN distributed states x (regulates skewness from Gaussian (Gamma=0) to half normal), i.e. Gamma_0_0
+% - nu_tm1_tm1 [skewx_dim by x_nbr] initial filtered value of skewness conditioning parameter of CSN distributed states x (enables closure of CSN distribution under conditioning), i.e. nu_0_0
+% - Delta_tm1_tm1 [skewx_dim by skewx_dim] initial filtered value of skewness marginalization parameter of CSN distributed states x (enables closure of CSN distribution under marginalization), i.e. Delta_0_0
+% - kalman_tol [double scalar] tolerance parameter for rcond and invertibility of covariance matrix
+% - riccati_tol [double scalar] tolerance parameter for iteration over the Riccati equation
+% - presample [integer scalar] number of initial iterations to be discarded when evaluating the likelihood
+% - diffuse_periods [integer scalar] number of diffuse filter periods in the initialization step
+% - G [x_nbr by x_nbr] state transition matrix mapping previous states to current states
+% - R [x_nbr by eta_nbr] state transition matrix mapping current innovations to current states
+% - F [y_nbr by x_nbr] observation equation matrix mapping current states into current observables (if Fflag=1)
+% [y_nbr by 1] indices for picking observables in transition matrix (if Fflag=0)
+% - Fflag [boolean] 0: if F is a vector of indices targeting observed variables in state vector
+% 1: if F is matrix in observation equation
+% - mu_eta [eta_nbr by 1] value of location parameter of CSN distributed innovations eta (does not equal expectation vector unless Gamma_eta=0), i.e. mu_eta
+% - Sigma_eta [eta_nbr by eta_nbr] value of scale parameter of CSN distributed innovations eta (does not equal covariance matrix unless Gamma_eta=0), i.e. Sigma_eta
+% - Gamma_eta [skeweta_dim by eta_nbr] value of skewness size parameter of CSN distributed innovations eta (regulates skewness from Gaussian (Gamma=0) to half normal), i.e. Gamma_eta
+% - nu_eta [skeweta_dim by eta_nbr] value of skewness conditioning parameter of CSN distributed innovations eta (enables closure of CSN distribution under conditioning), i.e. nu_eta
+% - Delta_eta [skeweta_dim by skeweta_dim] value of skewness marginalization parameter of CSN distributed innovations eta (enables closure of CSN distribution under marginalization), i.e. Delta_eta
+% - Sigma_eps [y_nbr by y_nbr] scale parameter of normally distributed measurement errors eps (equals covariance matrix), if no measurement errors this is a zero scalar
+% - G_singular [boolean] indicator if prediction step is done on joint distribution, [x_t', eta_t']', useful for models where G is singular
+% - prune_tol [double] correlation threshold to prune redundant skewness dimensions, if set to 0 no pruning will be done
+% - cdfmvna_fct [string] name of function to compute log Gaussian cdf, possible values: 'logmvncdf_ME', 'mvncdf', 'qsilatmvnv', 'qsimvnv'
+% -------------------------------------------------------------------------
+% OUTPUTS
+% - LIK [double scalar] value of negative log likelihood
+% - LIKK [(last-start+1) by 1] vector of densities for each observation
+% =========================================================================
+% Copyright (C) 2024 Gaygysyz Guljanov, Willi Mutschler
+%
+% This 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.
+%
+% This 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.
+% -------------------------------------------------------------------------
+% This file is part of the replication files for the paper "Pruned Skewed
+% Kalman Filter and Smoother: Pruned Skewed Kalman Filter and Smoother:
+% With Applications to the Yield Curve and Asymmetric Monetary Policy Shocks"
+% by Gaygysyz Guljanov, Willi Mutschler, Mark Trede
+% =========================================================================
+
+% Set defaults.
+verbose = false; % for development use
+
+if nargin < 25
+ cdfmvna_fct = 'logmvncdf_ME';
+end
+% initializations
+smpl = last-start+1; % sample size
+t = start; % time index
+likk = zeros(smpl,1); % vector gathering the densities
+LIK = Inf; % default value of the log likelihood
+oldK_Gauss = Inf; % old Gaussian Kalman gain for switching to steady-state Kalman filter (not yet implemented)
+oldK_Skewed = Inf; % old Skewed Kalman gain for switching to steady-state Kalman filter (not yet implemented)
+old_mu = Inf;
+old_Sigma = Inf;
+old_Gamma = Inf;
+old_nu = Inf;
+old_Delta = Inf;
+notsteady = [1 1]; % initialize indicator whether to break the loop and call the steady-state Kalman filter (not yet implemented, so just print convergence period)
+notsteady_Skew = [1 1]; % initialize indicator whether to break the loop and call the steady-state Kalman filter (not yet implemented, so just print convergence period)
+notsteady_mu = [1 1];
+notsteady_Sigma = [1 1];
+notsteady_Gamma = [1 1];
+notsteady_nu = [1 1];
+notsteady_Delta = [1 1];
+Omega_singular = true; % indicator whether Omega is singular
+LIKK = []; % output vector for individal density contributions
+rescale_prediction_error_covariance0 = rescale_prediction_error_covariance; % store option
+const2pi = -0.5*size(Y,1)*log(2*pi);
+
+if rcond(G) < kalman_tol
+ G_singular = true;
+ G_bar = [G, R]; % multiplying matrix of joint distribution
+else
+ if verbose
+ warning('G is non-singular');
+ end
+ G_singular = false;
+ mu_eta = R*mu_eta;
+ Sigma_eta = R*Sigma_eta*R';
+ Gamma_eta = Gamma_eta/(R'*R)*R';
+ Gamma_eta_X_Sigma_eta = Gamma_eta*Sigma_eta;
+ Delta22_common = Delta_eta + Gamma_eta_X_Sigma_eta*Gamma_eta';
+end
+
+while notsteady(1) && t<=last
+ s = t-start+1;
+
+ %%%%%%%%%%%%%%
+ % PREDICTION %
+ %%%%%%%%%%%%%%
+ if G_singular
+ % Parameters of joint distribution, [x_t', eta_t']'
+ mu_bar_tm1_tm1 = [mu_tm1_tm1; mu_eta];
+ nu_bar_tm1_tm1 = [nu_tm1_tm1; nu_eta];
+ Sigma_bar_tm1_tm1 = blkdiag_two(Sigma_tm1_tm1, Sigma_eta);
+ Gamma_bar_tm1_tm1 = blkdiag_two(Gamma_tm1_tm1, Gamma_eta);
+ Delta_bar_tm1_tm1 = blkdiag_two(Delta_tm1_tm1, Delta_eta);
+ % Linear transformation of the joint distribution
+ [mu_t_tm1, Sigma_t_tm1, Gamma_t_tm1, nu_t_tm1, Delta_t_tm1] = csn_statespace_linear_transform(G_bar, mu_bar_tm1_tm1, Sigma_bar_tm1_tm1, Gamma_bar_tm1_tm1, nu_bar_tm1_tm1, Delta_bar_tm1_tm1);
+ else
+ % auxiliary matrices
+ Gamma_tm1_tm1_X_Sigma_tm1_tm1 = Gamma_tm1_tm1*Sigma_tm1_tm1;
+ Gamma_tm1_tm1_X_Sigma_tm1_tm1_X_GT = Gamma_tm1_tm1_X_Sigma_tm1_tm1*G';
+ mu_t_tm1 = G*mu_tm1_tm1 + mu_eta;
+ Sigma_t_tm1 = G*Sigma_tm1_tm1*G' + Sigma_eta;
+ Sigma_t_tm1 = 0.5*(Sigma_t_tm1 + Sigma_t_tm1'); % ensure symmetry
+ invSigma_t_tm1 = inv(Sigma_t_tm1);
+ Gamma_t_tm1 = [Gamma_tm1_tm1_X_Sigma_tm1_tm1_X_GT; Gamma_eta_X_Sigma_eta]*invSigma_t_tm1;
+ nu_t_tm1 = [nu_tm1_tm1; nu_eta];
+ Delta11_t_tm1 = Delta_tm1_tm1 + Gamma_tm1_tm1_X_Sigma_tm1_tm1*Gamma_tm1_tm1' - Gamma_tm1_tm1_X_Sigma_tm1_tm1_X_GT*invSigma_t_tm1*Gamma_tm1_tm1_X_Sigma_tm1_tm1_X_GT';
+ Delta22_t_tm1 = Delta22_common - Gamma_eta_X_Sigma_eta*invSigma_t_tm1*Gamma_eta_X_Sigma_eta';
+ Delta12_t_tm1 = -Gamma_tm1_tm1_X_Sigma_tm1_tm1_X_GT*invSigma_t_tm1*Gamma_eta_X_Sigma_eta';
+ Delta_t_tm1 = [Delta11_t_tm1 , Delta12_t_tm1; Delta12_t_tm1' , Delta22_t_tm1];
+ Delta_t_tm1 = 0.5*(Delta_t_tm1 + Delta_t_tm1'); % ensure symmetry
+ end
+ % pruning redundant skewness dimensions to speed up filtering
+ if prune_tol > 0
+ [Sigma_t_tm1, Gamma_t_tm1, nu_t_tm1, Delta_t_tm1] = csnPruneParams(Sigma_t_tm1,Gamma_t_tm1,nu_t_tm1,Delta_t_tm1,prune_tol);
+ end
+ % Omega is Gaussian prediction error covariance
+ if Fflag
+ prediction_error = Y(:,t)-F*mu_t_tm1;
+ Omega = F*Sigma_t_tm1*F' + Sigma_eps;
+ else
+ prediction_error = Y(:,t)-mu_t_tm1(F);
+ Omega = Sigma_t_tm1(F,F) + Sigma_eps;
+ end
+ badly_conditioned_Omega = false;
+ if rescale_prediction_error_covariance
+ sigOmega = sqrt(diag(Omega)); % standard deviations
+ if any(diag(Omega)<kalman_tol) || rcond(Omega./(sigOmega*sigOmega'))<kalman_tol % Omega./(sigOmega*sigOmega') effectively converts Omega to a correlation matrix
+ badly_conditioned_Omega = true;
+ end
+ else
+ if rcond(Omega)<kalman_tol
+ sigOmega = sqrt(diag(Omega)); % standard deviations
+ if any(diag(Omega)<kalman_tol) || rcond(Omega./(sigOmega*sigOmega'))<kalman_tol % Omega./(sigOmega*sigOmega') effectively converts Omega to a correlation matrix
+ badly_conditioned_Omega = true;
+ else
+ rescale_prediction_error_covariance = 1;
+ end
+ end
+ end
+ if badly_conditioned_Omega
+ % if ~all(abs(Omega(:))<kalman_tol), then use univariate filter (will remove
+ % observations with zero variance prediction error), otherwise this is a
+ % pathological case and the draw is discarded
+ if verbose
+ warning('Omega is badly conditioned, discard draw as univariate filter will not be used (this overwrites Dynare''s default behavior');
+ end
+ return
+ else
+ Omega_singular = false;
+ if rescale_prediction_error_covariance % Omega needs to be rescaled to avoid numerical instability
+ % compute the log-determinant of covariance matrix Omega in a numerically stable way:
+ % let: corrOmega = Omega./(sigOmega*sigOmega') be the correlation matrix corresponding to Omega
+ % then: Omega = (sigOmega*sigOmega') .* corrOmega = diag(sigOmega) * corrOmega * diag(sigOmega)
+ % determinant: det(Omega) = det(diag(sigOmega)) * det(corrOmega) * det(diag(sigOmega)) = det(corrOmega) * prod(sigOmega)^2
+ % taking logs: log(det(Omega)) = log(det(corrOmega)) + 2*sum(log(sigOmega))
+ log_detOmega = log(det(Omega./(sigOmega*sigOmega'))) + 2*sum(log(sigOmega));
+ % compute inv(Omega) in a numerically stable way
+ % let: corrOmega = Omega./(sigOmega*sigOmega') be the correlation matrix corresponding to Omega
+ % then: Omega = (sigOmega*sigOmega') .* corrOmega = diag(sigOmega) * corrOmega * diag(sigOmega)
+ % inv(Omega) = diag(1./sigOmega) * inv(corrOmega) * diag(1./sigOmega) = inv(corrOmega) ./ (sigOmega*sigOmega')
+ invOmega = inv(Omega./(sigOmega*sigOmega'))./(sigOmega*sigOmega');
+ rescale_prediction_error_covariance = rescale_prediction_error_covariance0; % reset option as it might have been updated above
+ else
+ log_detOmega = log(det(Omega)); % compute the log-determinant of Omega directly
+ invOmega = inv(Omega); % compute the inverse of Omega directly
+ end
+
+ % Gaussian Kalman Gain
+ if Fflag
+ K_Gauss = Sigma_t_tm1*F'*invOmega;
+ else
+ K_Gauss = Sigma_t_tm1(:,F)*invOmega;
+ end
+ K_Skewed = Gamma_t_tm1*K_Gauss; % Skewed Kalman Gain
+
+ %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+ % LOG-LIKELIHOOD CONTRIBUTIONS %
+ %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+ % The conditional distribution of y(t) given y(t-1) is:
+ % (y(t)|y(t-1)) ~ CSN(mu_y,Sigma_y,Gamma_y,nu_y,Delta_y)
+ % = mvncdf(Gamma_y*(y(t)-mu_y),nu_y,Delta_y) / mvncdf(0,nu_y,Delta_y+Gamma_y*Sigma_y*Gamma_y') * mvnpdf(y(t),mu_y,Sigma_y)
+ % where:
+ % mu_y = F*mu_t_tm1 + mu_eps = y_predicted
+ % Sigma_y = F*Sigma_t_tm1*F' + Sigma_eps = Omega
+ % Gamma_y = Gamma_t_tm1*Sigma_t_tm1*F'*inv(F*Sigma_t_tm1*F' + Sigma_eps) = K_Skewed
+ % nu_y = nu_t_tm1
+ % Delta_y = Delta_t_tm1 + Gamma_t_tm1*Sigma_t_tm1*Gamma_t_tm1'...
+ % - Gamma_t_tm1*Sigma_t_tm1*F'*inv(F*Sigma_t_tm1*F')*F*Sigma_t_tm1*Gamma_t_tm1'...
+ % + (Gamma_t_tm1*Sigma_t_tm1*F'*inv(F*Sigma_t_tm1*F') - Gamma_t_tm1*Sigma_t_tm1*F'*inv(F*Sigma_t_tm1*F' + Sigma_eps))*F*Sigma_t_tm1*Gamma_t_tm1';
+ % = Delta_t_tm1 + (Gamma_t_tm1-K_Skewed*F)*Sigma_t_tm1*Gamma_t_tm1'
+ tmp = Gamma_t_tm1*Sigma_t_tm1;
+ Delta_y = Delta_t_tm1 + tmp*Gamma_t_tm1' - K_Skewed*F*tmp';
+ Delta_y = 0.5 * (Delta_y + Delta_y'); %ensure symmetry
+
+ % evaluate Gaussian cdfs, i.e.
+ % - bottom one: mvncdf(0,nu_y,Delta_y + Gamma_y*Sigma_y*Gamma_y')
+ % - top one: mvncdf(Gamma_y*(y(t)-mu_y),nu_y,Delta_y)
+
+ cdf_bottom_cov = Delta_y + K_Skewed*Omega*K_Skewed';
+ cdf_bottom_cov = 0.5*(cdf_bottom_cov + cdf_bottom_cov'); % ensure symmetry
+ if strcmp(cdfmvna_fct,'logmvncdf_ME')
+ % requires zero mean and correlation matrix as inputs
+ normalization_Delta_y = diag(1./sqrt(diag(cdf_bottom_cov)));
+ cdf_bottom_cov = normalization_Delta_y*cdf_bottom_cov*normalization_Delta_y; % this is now a correlation matrix!
+ cdf_bottom_cov = 0.5*(cdf_bottom_cov + cdf_bottom_cov'); % ensure symmetry
+ if ~isempty(cdf_bottom_cov)
+ try
+ log_gaussian_cdf_bottom = logmvncdf_ME(-normalization_Delta_y*nu_t_tm1, cdf_bottom_cov);
+ catch
+ if verbose
+ warning('log_gaussian_cdf_bottom something wrong');
+ end
+ return
+ end
+ else
+ log_gaussian_cdf_bottom = 0;
+ end
+
+ normalization_Delta_y = diag(1./sqrt(diag(Delta_y)));
+ Delta_y = normalization_Delta_y*Delta_y*normalization_Delta_y; % this is now a correlation matrix!
+ Delta_y = 0.5*(Delta_y + Delta_y'); % ensure symmetry
+ if ~isempty(Delta_y)
+ try
+ log_gaussian_cdf_top = logmvncdf_ME(normalization_Delta_y*(K_Skewed*prediction_error - nu_t_tm1), Delta_y);
+ catch
+ if verbose
+ warning('log_gaussian_cdf_top something wrong')
+ end
+ return
+ end
+ else
+ log_gaussian_cdf_top = 0;
+ end
+ elseif strcmp(cdfmvna_fct,'mvncdf')
+ try
+ log_gaussian_cdf_bottom = log(mvncdf(zeros(size(nu_t_tm1,1),1), nu_t_tm1, cdf_bottom_cov));
+ catch
+ if verbose
+ warning('log_gaussian_cdf_bottom something wrong')
+ end
+ return
+ end
+ try
+ log_gaussian_cdf_top = log(mvncdf(K_Skewed*prediction_error, nu_t_tm1, Delta_y));
+ catch
+ if verbose
+ warning('log_gaussian_cdf_top something wrong')
+ end
+ return
+ end
+ elseif strcmp(cdfmvna_fct,'qsilatmvnv')
+ k = size(nu_t_tm1,1);
+ if k > 1
+ try
+ log_gaussian_cdf_bottom = log(qsilatmvnv( 1000*k, cdf_bottom_cov, repmat(-Inf,k,1)-nu_t_tm1, zeros(k,1)-nu_t_tm1 ));
+ catch
+ if verbose
+ warning('log_gaussian_cdf_bottom something wrong');
+ end
+ return
+ end
+ try
+ log_gaussian_cdf_top = log(qsilatmvnv( 1000*k, Delta_y, repmat(-Inf,k,1)-nu_t_tm1, K_Skewed*prediction_error-nu_t_tm1 ));
+ catch
+ if verbose
+ warning('log_gaussian_cdf_top something wrong');
+ end
+ return
+ end
+ else
+ try
+ log_gaussian_cdf_bottom = log(normcdf(0, nu_t_tm1, cdf_bottom_cov));
+ catch
+ if verbose
+ warning('log_gaussian_cdf_bottom something wrong');
+ end
+ return
+ end
+ try
+ log_gaussian_cdf_top = log(normcdf(K_Skewed*prediction_error, nu_t_tm1, Delta_y));
+ catch
+ if verbose
+ warning('log_gaussian_cdf_top something wrong');
+ end
+ return
+ end
+ end
+ elseif strcmp(cdfmvna_fct,'qsimvnv')
+ k = size(nu_t_tm1,1);
+ if k > 1
+ try
+ log_gaussian_cdf_bottom = log(qsimvnv( 1000*k, cdf_bottom_cov, repmat(-Inf,k,1)-nu_t_tm1, zeros(k,1)-nu_t_tm1 ));
+ catch
+ if verbose
+ warning('log_gaussian_cdf_bottom something wrong');
+ end
+ return
+ end
+ try
+ log_gaussian_cdf_top = log(qsimvnv( 1000*k, Delta_y, repmat(-Inf,k,1)-nu_t_tm1, K_Skewed*prediction_error-nu_t_tm1 ));
+ catch
+ if verbose
+ warning('log_gaussian_cdf_top something wrong');
+ end
+ return
+ end
+ else
+ try
+ log_gaussian_cdf_bottom = log(normcdf(0, nu_t_tm1, cdf_bottom_cov));
+ catch
+ if verbose
+ warning('log_gaussian_cdf_bottom something wrong');
+ end
+ return
+ end
+ try
+ log_gaussian_cdf_top = log(normcdf(K_Skewed*prediction_error, nu_t_tm1, Delta_y));
+ catch
+ if verbose
+ warning('log_gaussian_cdf top something wrong');
+ end
+ return
+ end
+ end
+ end
+
+ % evaluate Gaussian pdf
+ % log_gaussian_pdf = log(mvnpdf(Y(:,t), y_predicted, Omega));
+ try
+ log_gaussian_pdf = const2pi - 0.5*log_detOmega - 0.5*transpose(prediction_error)*invOmega*prediction_error;
+ catch
+ if verbose
+ warning('log_gaussian_pdf something wrong');
+ end
+ return
+ end
+
+ % likelihood contribution
+ likk(s) = -log_gaussian_cdf_bottom + log_gaussian_pdf + log_gaussian_cdf_top;
+
+ %%%%%%%%%%%%%
+ % FILTERING %
+ %%%%%%%%%%%%%
+ mu_t_t = mu_t_tm1 + K_Gauss*prediction_error;
+ if Fflag
+ Sigma_t_t = Sigma_t_tm1 - K_Gauss*F*Sigma_t_tm1;
+ else
+ Sigma_t_t = Sigma_t_tm1 - K_Gauss*Sigma_t_tm1(F,:);
+ end
+ Gamma_t_t = Gamma_t_tm1;
+ nu_t_t = nu_t_tm1 - K_Skewed*prediction_error;
+ Delta_t_t = Delta_t_tm1;
+ Sigma_t_t = 0.5*(Sigma_t_t + Sigma_t_t'); % ensure symmetry
+ Delta_t_t = 0.5*(Delta_t_t + Delta_t_t'); % ensure symmetry
+
+ % assign for next time step
+ mu_tm1_tm1 = mu_t_t;
+ Sigma_tm1_tm1 = Sigma_t_t;
+ Gamma_tm1_tm1 = Gamma_t_t;
+ nu_tm1_tm1 = nu_t_t;
+ Delta_tm1_tm1 = Delta_t_t;
+
+ % check for steady-state Kalman filter (not yet implemented)
+ if verbose
+ if max(abs(K_Gauss(:)-oldK_Gauss)) < riccati_tol
+ if notsteady(2) == 1
+ fprintf('Gaussian Kalman gain converged after %d/%d periods (but steady-state PSKF filter is not implemented yet)\n',s,smpl);
+ notsteady(2) = 0;
+ end
+ else
+ if notsteady(2) == 0
+ error('Gaussian Kalman gain diverged after %d/%d periods\n',s,smpl);
+ notsteady(2) = 1;
+ end
+ end
+ if numel(K_Skewed) == numel(oldK_Skewed) && (max(abs(K_Skewed(:)-oldK_Skewed)) < riccati_tol)
+ if notsteady_Skew(2) == 1
+ fprintf('Skewed Kalman gain converged after %d/%d periods with skewness dimension: %d\n',s,smpl,size(K_Skewed,1));
+ notsteady_Skew(2) = 0;
+ end
+ else
+ if notsteady_Skew(2) == 0
+ error('Skewed Kalman gain diverged after %d/%d periods\n',s,smpl);
+ notsteady_Skew(2) = 1;
+ end
+ end
+ end
+ oldK_Gauss = K_Gauss(:);
+ oldK_Skewed = K_Skewed(:);
+ old_mu = mu_t_t(:);
+ old_Sigma = Sigma_t_t(:);
+ old_Gamma = Gamma_t_t(:);
+ old_nu = nu_t_t(:);
+ old_Delta = Delta_t_t(:);
+ end
+ t = t+1;
+end
+
+if verbose
+ fprintf('Skewness Dimension: %d\n',size(Gamma_t_t,1));
+end
+
+if Omega_singular
+ error('The variance of the forecast error remains singular until the end of the sample')
+end
+
+% Compute minus the log-likelihood.
+if presample>diffuse_periods
+ LIK = -1*sum(likk(1+(presample-diffuse_periods):end));
+else
+ LIK = -1*sum(likk);
+end
+LIKK = -1*likk;
+
+
+%% auxiliary functions
+function res_mat = blkdiag_two(mat1, mat2)
+ % Makes a block diagonal matrix out of two matrices
+ [nrow_mat1, ncol_mat1] = size(mat1); [nrow_mat2, ncol_mat2] = size(mat2);
+ upper_mat = zeros(nrow_mat1, ncol_mat2);
+ lower_mat = zeros(nrow_mat2, ncol_mat1);
+ res_mat = [mat1, upper_mat; lower_mat, mat2];
+end % blkdiag_two
+
+
+end % kalman_filter_pruned_skewed
\ No newline at end of file
diff --git a/matlab/name2index.m b/matlab/name2index.m
index 78b121690918eb5e95f04eaef07a724e440e2977..ce35b765df27233269406a70f73309520216f649 100644
--- a/matlab/name2index.m
+++ b/matlab/name2index.m
@@ -1,7 +1,8 @@
function i = name2index(options_, M_, estim_params_, type, name1, name2 )
% Returns the index associated to an estimated object (deep parameter,
% variance of a structural shock or measurement error, covariance between
-% two structural shocks, covariance between two measurement errors).
+% two structural shocks, covariance between two measurement errors, or
+% skewness of structural shocks).
%
% INPUTS:
% options_ [structure] Dynare structure.
@@ -16,7 +17,7 @@ function i = name2index(options_, M_, estim_params_, type, name1, name2 )
% SPECIAL REQUIREMENTS
% none
-% Copyright © 2008-2018 Dynare Team
+% Copyright © 2008-2024 Dynare Team
%
% This file is part of Dynare.
%
@@ -37,13 +38,14 @@ nvx = estim_params_.nvx;
nvn = estim_params_.nvn;
ncx = estim_params_.ncx;
ncn = estim_params_.ncn;
+nsx = estim_params_.nsx;
npa = estim_params_.np ;
-nnn = nvx+nvn+ncx+ncn+npa;
+nnn = nvx+nvn+ncx+ncn+nsx+npa;
i = [];
if strcmpi(type,'DeepParameter')
- i = nvx + nvn + ncx + ncn + ...
+ i = nvx + nvn + ncx + ncn + nsx + ...
strmatch(name1, M_.param_names(estim_params_.param_vals(:,1)), 'exact');
if nargin>5
disp('The last input argument is useless!')
@@ -114,4 +116,12 @@ if strcmpi(type,'MeasurementError')
disp('Off diagonal terms of the covariance matrix are not estimated (measurement equation)')
end
end
+end
+
+if strcmpi(type,'StructuralShockSkewness')
+ i = strmatch(name1, M_.exo_names(estim_params_.skew_exo(:,1)), 'exact');
+ if isempty(i)
+ disp(['The skewness of ' name1 ' is not an estimated parameter!'])
+ return
+ end
end
\ No newline at end of file
diff --git a/matlab/optimization/transform_to_bounded.m b/matlab/optimization/transform_to_bounded.m
new file mode 100644
index 0000000000000000000000000000000000000000..d2b087c69c0a438c0e1100ef3fc4cf3e350a1658
--- /dev/null
+++ b/matlab/optimization/transform_to_bounded.m
@@ -0,0 +1,40 @@
+function X = transform_to_bounded(Y,a,b)
+% function X = transform_to_bounded(Y,a,b)
+% -------------------------------------------------------------------------
+% transforms an unbounded parameter vector back to original bounded parameters
+% using an inverse log-odds transform as described at
+% https://mc-stan.org/docs/reference-manual/logit-transform-jacobian.html
+% -------------------------------------------------------------------------
+% INPUTS
+% - Y [double vector] unbounded parameter vector
+% - a [double vector] lower bound of original parameters
+% - b [double vector] upper bound of original parameters
+% -------------------------------------------------------------------------
+% OUTPUTS
+% - X [double vector] original bounded parameter vector
+% =========================================================================
+% Copyright (C) 2023 Willi Mutschler
+%
+% This 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.
+%
+% This 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.
+% -------------------------------------------------------------------------
+% This file is part of the replication files for the paper "Pruned Skewed
+% Kalman Filter and Smoother: Pruned Skewed Kalman Filter and Smoother:
+% With Applications to the Yield Curve and Asymmetric Monetary Policy Shocks"
+% by Gaygysyz Guljanov, Willi Mutschler, Mark Trede
+% =========================================================================
+
+X = a + (b-a).*inv_logit(Y);
+
+function u = inv_logit(v)
+ u = 1./(1+exp(-v));
+end %inv_logit
+
+end % transform_to_unbounded
\ No newline at end of file
diff --git a/matlab/optimization/transform_to_bounded_derivative.m b/matlab/optimization/transform_to_bounded_derivative.m
new file mode 100644
index 0000000000000000000000000000000000000000..2538e3d260aae9ff146f524598406d518854b4cf
--- /dev/null
+++ b/matlab/optimization/transform_to_bounded_derivative.m
@@ -0,0 +1,39 @@
+function dX_dY = transform_to_bounded_derivative(Y,a,b)
+% function dX_dY = transform_to_bounded_derivative(Y,a,b)
+% -------------------------------------------------------------------------
+% Jacobian of inverse logit transform as described at
+% https://mc-stan.org/docs/reference-manual/logit-transform-jacobian.html
+% -------------------------------------------------------------------------
+% INPUTS
+% - Y [double vector] unbounded parameter vector
+% - a [double vector] lower bound of original parameters
+% - b [double vector] upper bound of original parameters
+% -------------------------------------------------------------------------
+% OUTPUTS
+% - dX_dY [double vector] computes Jacobian of original bounded parameter vector
+% wrt to unbounded parameters
+% =========================================================================
+% Copyright (C) 2023 Willi Mutschler
+%
+% This 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.
+%
+% This 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.
+% -------------------------------------------------------------------------
+% This file is part of the replication files for the paper "Pruned Skewed
+% Kalman Filter and Smoother: Pruned Skewed Kalman Filter and Smoother:
+% With Applications to the Yield Curve and Asymmetric Monetary Policy Shocks"
+% by Gaygysyz Guljanov, Willi Mutschler, Mark Trede
+% =========================================================================
+dX_dY = (b-a).*inv_logit(Y).*(1-inv_logit(Y));
+
+function u = inv_logit(v)
+ u = 1./(1+exp(-v));
+end %inv_logit
+
+end % transform_to_unbounded
\ No newline at end of file
diff --git a/matlab/optimization/transform_to_unbounded.m b/matlab/optimization/transform_to_unbounded.m
new file mode 100644
index 0000000000000000000000000000000000000000..ed0979594ab078cb9d04785e79a40065c73939ff
--- /dev/null
+++ b/matlab/optimization/transform_to_unbounded.m
@@ -0,0 +1,39 @@
+function Y = transform_to_unbounded(X,a,b)
+% function Y = transform_to_unbounded(X,a,b)
+% -------------------------------------------------------------------------
+% transforms an original bounded parameter vector to unbounded parameters
+% using a log-odds transform as described at
+% https://mc-stan.org/docs/reference-manual/logit-transform-jacobian.html
+% -------------------------------------------------------------------------
+% INPUTS
+% - X [double vector] original bounded parameter vector
+% - a [double vector] lower bound of original parameters
+% - b [double vector] upper bound of original parameters
+% -------------------------------------------------------------------------
+% OUTPUTS
+% - Y [double vector] unbounded parameter vector
+% =========================================================================
+% Copyright (C) 2023 Willi Mutschler
+%
+% This 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.
+%
+% This 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.
+% -------------------------------------------------------------------------
+% This file is part of the replication files for the paper "Pruned Skewed
+% Kalman Filter and Smoother: Pruned Skewed Kalman Filter and Smoother:
+% With Applications to the Yield Curve and Asymmetric Monetary Policy Shocks"
+% by Gaygysyz Guljanov, Willi Mutschler, Mark Trede
+% =========================================================================
+Y = logit((X-a)./(b-a));
+
+function logit_u = logit(u)
+ logit_u = log(u./(1-u));
+end % logit
+
+end % transform_to_unbounded
\ No newline at end of file
diff --git a/matlab/row_header_width.m b/matlab/row_header_width.m
index 763be482fbba52f3b51b2e1de9dd4d6d084e2a93..d77b5fe4906ea8640d4a321881571d614486fef5 100644
--- a/matlab/row_header_width.m
+++ b/matlab/row_header_width.m
@@ -34,6 +34,7 @@ np = estim_params_.np;
nvx = estim_params_.nvx;
nvn = estim_params_.nvn;
ncx = estim_params_.ncx;
+nsx = estim_params_.nsx;
ncn = estim_params_.ncn;
w = 0;
@@ -41,7 +42,7 @@ if np
w = cellofchararraymaxlength(bayestopt_.name);
end
if nvx
- w = max(w, cellofchararraymaxlength(M_.endo_names(estim_params_.var_exo(1:nvx,1))));
+ w = max(w, cellofchararraymaxlength(M_.exo_names(estim_params_.var_exo(1:nvx,1))));
end
if nvn
w = max(w, cellofchararraymaxlength(M_.endo_names(estim_params_.var_endo(1:nvn,1))));
@@ -61,3 +62,6 @@ if ncn
w = max(w, length(M_.endo_names{k1})+length(M_.endo_names{k2}));
end
end
+if nsx
+ w = max(w, cellofchararraymaxlength(M_.exo_names(estim_params_.skew_exo(1:nsx,1))));
+end
\ No newline at end of file
diff --git a/matlab/set_all_parameters.m b/matlab/set_all_parameters.m
index ed090a1768b04cfa2f6db10d3b6d9f67743d85b8..0da0a946e657a48d9d460d360bbfa7f925948fb3 100644
--- a/matlab/set_all_parameters.m
+++ b/matlab/set_all_parameters.m
@@ -54,6 +54,7 @@ nvx = estim_params_.nvx;
ncx = estim_params_.ncx;
nvn = estim_params_.nvn;
ncn = estim_params_.ncn;
+nsx = estim_params_.nsx;
np = estim_params_.np;
if nvx || ncx
Sigma_e = M_.Sigma_e;
@@ -70,8 +71,7 @@ if nvx
Sigma_e(k,k) = xparam1(i)^2;
end
end
-% update offset
-offset = nvx;
+offset = nvx; % update offset
% setting measument error variance; on the diagonal of Covariance matrix; used later
% for updating covariances
@@ -81,9 +81,7 @@ if nvn
H(k,k) = xparam1(i+offset)^2;
end
end
-
-% update offset
-offset = nvx+nvn;
+offset = nvx+nvn; % update offset
% setting shocks covariances
if ncx
@@ -95,8 +93,7 @@ if ncx
Correlation_matrix(k2,k1) = Correlation_matrix(k1,k2);
end
end
-% update offset
-offset = nvx+nvn+ncx;
+offset = nvx+nvn+ncx; % update offset
% setting measurement error covariances
if ncn
@@ -108,11 +105,18 @@ if ncn
Correlation_matrix_ME(k2,k1) = Correlation_matrix_ME(k1,k2);
end
end
+offset = nvx+nvn+ncx+ncn; % update offset
+
+% setting skewness coefficient of structural shocks
+if nsx
+ skew_exo = estim_params_.skew_exo;
+ for i=1:nsx
+ M_.Skew_e(skew_exo(i,1)) = xparam1(i+offset);
+ end
+end
+offset = nvx+nvn+ncx+ncn+nsx; % update offset
-% update offset
-offset = nvx+ncx+nvn+ncn;
% setting structural parameters
-%
if np
M_.params(estim_params_.param_vals(:,1)) = xparam1(offset+1:end);
end
@@ -139,4 +143,26 @@ if nvn || ncn
end
M_.H = H;
M_.Correlation_matrix_ME=Correlation_matrix_ME;
+end
+% update parameters of CSN distributed shocks
+M_.csn.error_indicator = zeros(M_.exo_nbr,1);
+for jexo = 1:M_.exo_nbr
+ if M_.Skew_e(jexo,1) == 0 % Gaussian
+ M_.csn.Sigma_e(jexo,jexo) = M_.Sigma_e(jexo,jexo);
+ M_.csn.Gamma_e(jexo,jexo) = 0;
+ else % CSN
+ if abs(M_.Skew_e(jexo,1)) >= 0.995
+ %warning('Skewness parameter very close to theoretical bound')
+ end
+ [M_.csn.Sigma_e(jexo,jexo),M_.csn.Gamma_e(jexo,jexo),M_.csn.error_indicator(jexo)] = csnVarSkew_To_SigmaGamma_univariate(M_.Sigma_e(jexo,jexo),M_.Skew_e(jexo,1),false);
+ end
+end
+if all(M_.csn.error_indicator == 0)
+ if iszero(diag(M_.csn.Gamma_e))
+ M_.csn.mu_e = zeros(M_.exo_nbr,1);
+ else
+ M_.csn.mu_e = -csnMean(zeros(M_.exo_nbr,1),M_.csn.Sigma_e,M_.csn.Gamma_e,M_.csn.nu_e,M_.csn.Delta_e); % ensures E[shocks]=0
+ end
+else
+ M_.csn.mu_e = [];
end
\ No newline at end of file
diff --git a/matlab/set_parameters_locally.m b/matlab/set_parameters_locally.m
index b3fce1696a6d2fe5bd3c759879aedc444711be0e..75669f0fb177f755acd2e8f441d1aaf0acbaa2ba 100644
--- a/matlab/set_parameters_locally.m
+++ b/matlab/set_parameters_locally.m
@@ -39,6 +39,7 @@ nvx = estim_params_.nvx;
ncx = estim_params_.ncx;
nvn = estim_params_.nvn;
ncn = estim_params_.ncn;
+nsx = estim_params_.nsx;
np = estim_params_.np;
Sigma_e = M_.Sigma_e;
Correlation_matrix = M_.Correlation_matrix;
@@ -76,10 +77,42 @@ end
% and update offset
offset = offset + ncx + ncn;
+% setting skewness coefficient of structural shocks
+if nsx
+ skew_exo = estim_params_.skew_exo;
+ for i=1:nsx
+ M_.Skew_e(skew_exo(i,1)) = xparam1(i+offset);
+ end
+end
+offset = offset + nsx; % update offset
+
% structural parameters
if np
M_.params(estim_params_.param_vals(:,1)) = xparam1(offset+1:end);
end
M_.Sigma_e = Sigma_e;
-M_.Correlation_matrix=Correlation_matrix;
\ No newline at end of file
+M_.Correlation_matrix=Correlation_matrix;
+
+% update parameters of CSN distributed shocks
+M_.csn.error_indicator = zeros(M_.exo_nbr,1);
+for jexo = 1:M_.exo_nbr
+ if M_.Skew_e(jexo,1) == 0 % Gaussian
+ M_.csn.Sigma_e(jexo,jexo) = M_.Sigma_e(jexo,jexo);
+ M_.csn.Gamma_e(jexo,jexo) = 0;
+ else % CSN
+ if abs(M_.Skew_e(jexo,1)) >= 0.995
+ %warning('Skewness parameter very close to theoretical bound')
+ end
+ [M_.csn.Sigma_e(jexo,jexo),M_.csn.Gamma_e(jexo,jexo),M_.csn.error_indicator(jexo)] = csnVarSkew_To_SigmaGamma_univariate(M_.Sigma_e(jexo,jexo),M_.Skew_e(jexo,1),false);
+ end
+end
+if all(M_.csn.error_indicator == 0)
+ if iszero(diag(M_.csn.Gamma_e))
+ M_.csn.mu_e = zeros(M_.exo_nbr,1);
+ else
+ M_.csn.mu_e = -csnMean(zeros(M_.exo_nbr,1),M_.csn.Sigma_e,M_.csn.Gamma_e,M_.csn.nu_e,M_.csn.Delta_e); % ensures E[shocks]=0
+ end
+else
+ M_.csn.mu_e = [];
+end
\ No newline at end of file
diff --git a/matlab/stochastic_solver/simult.m b/matlab/stochastic_solver/simult.m
index 4787d0e6198ff69480b7aed25699dcbdc753e460..e98e07aac1aadb387cb70a5bf07f8b66f1df0818 100644
--- a/matlab/stochastic_solver/simult.m
+++ b/matlab/stochastic_solver/simult.m
@@ -77,11 +77,42 @@ nxs = length(i_exo_var);
exo_simul = zeros(options_.periods,M_.exo_nbr);
chol_S = chol(M_.Sigma_e(i_exo_var,i_exo_var));
+% update parameters of CSN distributed shocks with unit variance for simulation (as we re-scale with chol_S)
+M_.csn.error_indicator = zeros(M_.exo_nbr,1);
+for jexo = 1:M_.exo_nbr
+ if M_.Skew_e(jexo,1) == 0 % Gaussian
+ M_.csn.Sigma_e(jexo,jexo) = 1;
+ M_.csn.Gamma_e(jexo,jexo) = 0;
+ else % CSN
+ if abs(M_.Skew_e(jexo,1)) >= 0.995
+ %warning('Skewness parameter very close to theoretical bound')
+ end
+ [M_.csn.Sigma_e(jexo,jexo),M_.csn.Gamma_e(jexo,jexo),M_.csn.error_indicator(jexo)] = csnVarSkew_To_SigmaGamma_univariate(1,M_.Skew_e(jexo,1),true);
+ end
+end
+M_.csn.mu_e = nan(M_.exo_nbr,1);
+if all(M_.csn.error_indicator(i_exo_var) == 0)
+ if iszero(diag(M_.csn.Gamma_e(i_exo_var,i_exo_var)))
+ M_.csn.mu_e(i_exo_var,1) = zeros(nxs,1);
+ else
+ M_.csn.mu_e(i_exo_var,1) = -csnMean(zeros(nxs,1),M_.csn.Sigma_e(i_exo_var,i_exo_var),M_.csn.Gamma_e(i_exo_var,i_exo_var),M_.csn.nu_e(i_exo_var,1),M_.csn.Delta_e(i_exo_var,i_exo_var)); % ensures E[shocks]=0
+ end
+else
+ M_.csn.mu_e = [];
+end
+
for i=1:replic
if ~isempty(M_.Sigma_e)
% we fill the shocks row wise to have the same values
% independently of the length of the simulation
- exo_simul(:,i_exo_var) = randn(nxs,options_.periods)'*chol_S;
+ if iszero(diag(M_.csn.Gamma_e(i_exo_var,i_exo_var)))
+ exo_simul(:,i_exo_var) = transpose(randn(nxs,options_.periods))*chol_S;
+ else
+ if ~isdiag(chol_S)
+ error('Correlated CSN shocks are currently not supported.')
+ end
+ exo_simul(:,i_exo_var) = transpose(csnRandDirect(options_.periods, M_.csn.mu_e(i_exo_var,1), M_.csn.Sigma_e(i_exo_var,i_exo_var), M_.csn.Gamma_e(i_exo_var,i_exo_var), M_.csn.nu_e(i_exo_var,1), M_.csn.Delta_e(i_exo_var,i_exo_var)))*chol_S;
+ end
end
y_ = simult_(M_,options_,y0,dr,exo_simul,order);
% elimninating initial value
diff --git a/matlab/stochastic_solver/stoch_simul.m b/matlab/stochastic_solver/stoch_simul.m
index 376de52368289594b0e0a1abfa3f4275cd53a234..fba0fb77ea79e75a09c5388ff1473fa87586bea5 100644
--- a/matlab/stochastic_solver/stoch_simul.m
+++ b/matlab/stochastic_solver/stoch_simul.m
@@ -141,6 +141,19 @@ if ~options_.noprint
lh = cellofchararraymaxlength(labels)+2;
dyn_latex_table(M_, options_, my_title, 'covar_ex_shocks', headers, labels, M_.Sigma_e, lh, 10, 6);
end
+ if ~all(M_.Skew_e==0)
+ my_title='SKEWNESS OF EXOGENOUS SHOCKS';
+ labels = M_.exo_names;
+ headers = vertcat('Variables', {'Skewness'});
+ lh = cellofchararraymaxlength(labels)+2;
+ dyntable(options_, my_title, headers, labels, M_.Skew_e, lh, 10, 6);
+ if options_.TeX
+ labels = M_.exo_names_tex;
+ headers = vertcat('Variables', {'Skewness'});
+ lh = cellofchararraymaxlength(labels)+2;
+ dyn_latex_table(M_, options_, my_title, 'skew_ex_shocks', headers, labels, M_.Skew_e, lh, 10, 6);
+ end
+ end
if ~all(diag(M_.H)==0)
my_title='MATRIX OF COVARIANCE OF MEASUREMENT ERRORS';
labels = cellfun(@(x) horzcat('SE_', x), options_.varobs, 'UniformOutput', false);
diff --git a/matlab/utilities/estimation/check_prior_stderr_corr.m b/matlab/utilities/estimation/check_prior_stderr_corr_skew.m
similarity index 81%
rename from matlab/utilities/estimation/check_prior_stderr_corr.m
rename to matlab/utilities/estimation/check_prior_stderr_corr_skew.m
index 867e0ba7cd4876dc1c53f4ae5f72b01dda6305a0..23d5dd6d66c40ce0f52943704872a2a6d7bd59ef 100644
--- a/matlab/utilities/estimation/check_prior_stderr_corr.m
+++ b/matlab/utilities/estimation/check_prior_stderr_corr_skew.m
@@ -1,8 +1,9 @@
-function check_prior_stderr_corr(estim_params_,bayestopt_)
-% function check_prior_stderr_corr(estim_params_,bayestopt_)
+function check_prior_stderr_corr_skew(estim_params_,bayestopt_)
+% function check_prior_stderr_corr_skew(estim_params_,bayestopt_)
% -------------------------------------------------------------------------
-% Check if the prior allows for negative standard deviations and
-% correlations larger than +-1. If so, issue a warning.
+% Check if the prior allows for negative standard deviations,
+% correlations larger than +-1, or skewness larger than +-1.
+% If so, issue a warning.
% -------------------------------------------------------------------------
% INPUTS
% o estim_params_: [struct] information on estimated parameters
@@ -50,4 +51,9 @@ offset = nvx+nvn+ncx;
ncn = estim_params_.ncn; % number of corr parameters for measurement errors
if ncn && (any(bayestopt_.p3(1+offset:offset+ncn)<-1) || any(bayestopt_.p4(1+offset:offset+ncn)>1))
warning('Your prior allows for correlations between measurement errors larger than +-1 and will not integrate to 1 due to truncation. Please change your prior')
+end
+offset = nvx+nvn+ncx+ncn;
+nsx = estim_params_.nsx; % number of skewness parameters for structural shocks
+if nsx && (any(bayestopt_.p3(1+offset:offset+nsx)<-1) || any(bayestopt_.p4(1+offset:offset+nsx)>1))
+ warning('Your prior allows for skewness parameters larger than +-1 and will not integrate to 1 due to truncation. Please change your prior')
end
\ No newline at end of file