Commit eba6164c authored by michel's avatar michel
Browse files

fixing bugs around missing observations and starting to add smoother with missing observations

git-svn-id: https://www.dynare.org/svn/dynare/dynare_v4@2274 ac1d8469-bf42-47a9-8791-bf33cf982152
parent dd73be5b
......@@ -249,16 +249,6 @@ function [fval,cost_flag,ys,trend_coeff,info] = DsgeLikelihood(xparam1,gend,data
end
end
if kalman_algo == 2
no_correlation_flag = 1;
if length(H)==1
H = zeros(nobs,1);
else
if all(all(abs(H-diag(diag(H)))<1e-14))% ie, the covariance matrix is diagonal...
H = diag(H);
else
no_correlation_flag = 1;
end
end
end
kalman_tol = options_.kalman_tol;
riccati_tol = options_.riccati_tol;
......@@ -280,6 +270,16 @@ function [fval,cost_flag,ys,trend_coeff,info] = DsgeLikelihood(xparam1,gend,data
end
end
if (kalman_algo==2)% Univariate Kalman Filter
no_correlation_flag = 1;
if length(H)==1 & H == 0
H = zeros(nobs,1);
else
if all(all(abs(H-diag(diag(H)))<1e-14))% ie, the covariance matrix is diagonal...
H = diag(H);
else
no_correlation_flag = 0;
end
end
if no_correlation_flag
LIK = univariate_kalman_filter(T,R,Q,H,Pstar,Y,start,mf,kalman_tol,riccati_tol,data_index,number_of_observations,no_more_missing_observations);
else
......@@ -288,29 +288,38 @@ function [fval,cost_flag,ys,trend_coeff,info] = DsgeLikelihood(xparam1,gend,data
end
if (kalman_algo==3)% Multivariate Diffuse Kalman Filter
if no_missing_data_flag
data1 = data - trend;
if any(any(H ~= 0))
LIK = DiffuseLikelihoodH1_Z(ST,Z,R1,Q,H,Pinf,Pstar,data1,start);
else
LIK = DiffuseLikelihood1_Z(ST,Z,R1,Q,Pinf,Pstar,data1,start);
end
if isinf(LIK)
kalman_algo = 4;
end
LIK = diffuse_kalman_filter(ST,R1,Q,H,Pinf,Pstar,Y,start,Z,kalman_tol, ...
riccati_tol);
else
error(['The diffuse filter is not yet implemented for models with missing observations'])
LIK = missing_observations_diffuse_kalman_filter(ST,R1,Q,H,Pinf, ...
Pstar,Y,start,Z,kalman_tol,riccati_tol,...
data_index,number_of_observations,...
no_more_missing_observations);
end
if isinf(LIK)
kalman_algo = 4;
end
end
if (kalman_algo==4)% Univariate Diffuse Kalman Filter
data1 = data - trend;
if any(any(H ~= 0))
if ~estim_params_.ncn
LIK = DiffuseLikelihoodH3_Z(ST,Z,R1,Q,H,Pinf,Pstar,data1,trend,start);
no_correlation_flag = 1;
if length(H)==1 & H == 0
H = zeros(nobs,1);
else
if all(all(abs(H-diag(diag(H)))<1e-14))% ie, the covariance matrix is diagonal...
H = diag(H);
else
LIK = DiffuseLikelihoodH3corr_Z(ST,Z,R1,Q,H,Pinf,Pstar,data1,trend,start);
no_correlation_flag = 0;
end
end
if no_correlation_flag
LIK = univariate_diffuse_kalman_filter(ST,R1,Q,H,Pinf,Pstar,Y, ...
start,Z,kalman_tol,riccati_tol,data_index,...
number_of_observations,no_more_missing_observations);
else
LIK = DiffuseLikelihood3_Z(ST,Z,R1,Q,Pinf,Pstar,data1,start);
LIK = univariate_diffuse_kalman_filter_corr(ST,R1,Q,H,Pinf,Pstar, ...
Y,start,Z,kalman_tol,riccati_tol,...
data_index,number_of_observations,...
no_more_missing_observations);
end
end
if imag(LIK) ~= 0
......
function [alphahat,etahat,epsilonhat,ahat,SteadyState,trend_coeff,aK,T,R,P,PK,d,decomp] = DsgeSmoother(xparam1,gend,Y)
function [alphahat,etahat,epsilonhat,ahat,SteadyState,trend_coeff,aK,T,R,P,PK,d,decomp] = DsgeSmoother(xparam1,gend,Y,data_index,missing_value)
% Estimation of the smoothed variables and innovations.
%
% INPUTS
% o xparam1 [double] (p*1) vector of (estimated) parameters.
% o gend [integer] scalar specifying the number of observations ==> varargin{1}.
% o data [double] (T*n) matrix of data.
% o data_index [cell] 1*smpl cell of column vectors of indices.
% o missing_value 1 if missing values, 0 otherwise
%
% OUTPUTS
% o alphahat [double] (m*T) matrix, smoothed endogenous variables.
......@@ -265,26 +267,37 @@ function [alphahat,etahat,epsilonhat,ahat,SteadyState,trend_coeff,aK,T,R,P,PK,d,
[alphahat,etahat,ahat,aK] = DiffuseKalmanSmoother1(T,R,Q,Pinf,Pstar,Y,trend,nobs,np,smpl,mf);
if all(alphahat(:)==0)
kalman_algo = 2;
[alphahat,etahat,ahat,aK] = DiffuseKalmanSmoother3(T,R,Q,Pinf,Pstar,Y,trend,nobs,np,smpl,mf);
end
elseif kalman_algo == 2
end
if kalman_algo == 2
[alphahat,etahat,ahat,aK] = DiffuseKalmanSmoother3(T,R,Q,Pinf,Pstar,Y,trend,nobs,np,smpl,mf);
elseif kalman_algo == 3 | kalman_algo == 4
end
if kalman_algo == 3
data1 = Y - trend;
if options_.kalman_algo == 3
if missing_value
[alphahat,etahat,ahat,P,aK,PK,d,decomp] = missing_DiffuseKalmanSmoother1_Z(ST, ...
Z,R1,Q,Pinf,Pstar,data1,nobs,np,smpl,data_index);
else
[alphahat,etahat,ahat,P,aK,PK,d,decomp] = DiffuseKalmanSmoother1_Z(ST, ...
Z,R1,Q,Pinf,Pstar,data1,nobs,np,smpl);
if all(alphahat(:)==0)
options_.kalman_algo = 4;
[alphahat,etahat,ahat,P,aK,PK,d,decomp] = DiffuseKalmanSmoother3_Z(ST, ...
Z,R1,Q,Pinf,Pstar, ...
data1,nobs,np,smpl);
end
% oo_.NewFilteredVariables = QT*filtered_values(data1,ahat,P,T,Z);
end
if all(alphahat(:)==0)
options_.kalman_algo = 4;
end
end
if kalman_algo == 4
data1 = Y - trend;
if missing_value
[alphahat,etahat,ahat,P,aK,PK,d,decomp] = missing_DiffuseKalmanSmoother3_Z(ST, ...
Z,R1,Q,Pinf,Pstar,data1,nobs,np,smpl,data_index);
else
[alphahat,etahat,ahat,P,aK,PK,d,decomp] = DiffuseKalmanSmoother3_Z(ST, ...
Z,R1,Q,Pinf,Pstar,data1,nobs,np,smpl);
Z,R1,Q,Pinf,Pstar, ...
data1,nobs,np,smpl);
end
end
if kalman_algo == 3 | kalman_algo == 4
alphahat = QT*alphahat;
ahat = QT*ahat;
nk = options_.nk;
......
function PosteriorFilterSmootherAndForecast(Y,gend, type)
function PosteriorFilterSmootherAndForecast(Y,gend, type,data_index)
% function PosteriorFilterSmootherAndForecast(Y,gend, type)
% Computes posterior filter smoother and forecasts
......@@ -138,7 +138,7 @@ for b=1:B
set_all_parameters(deep);
dr = resol(oo_.steady_state,0);
[alphahat,etahat,epsilonhat,ahat,SteadyState,trend_coeff,aK] = ...
DsgeSmoother(deep,gend,Y);
DsgeSmoother(deep,gend,Y,data_index);
if options_.loglinear
stock_smooth(dr.order_var,:,irun1) = alphahat(1:endo_nbr,:)+ ...
......
......@@ -315,12 +315,13 @@ if options_.bvar_dsge
end
[data_index,number_of_observations,no_more_missing_observations] = describe_missing_data(data,gend,n_varobs);
missing_value = ~(number_of_observations == gend*n_varobs);
initial_estimation_checks(xparam1,gend,data,data_index,number_of_observations,no_more_missing_observations);
if options_.mode_compute == 0 & length(options_.mode_file) == 0
if options_.smoother == 1
[atT,innov,measurement_error,updated_variables,ys,trend_coeff,aK,T,R,P,PK,d,decomp] = DsgeSmoother(xparam1,gend,data);
[atT,innov,measurement_error,updated_variables,ys,trend_coeff,aK,T,R,P,PK,d,decomp] = DsgeSmoother(xparam1,gend,data,data_index,missing_value);
oo_.Smoother.SteadyState = ys;
oo_.Smoother.TrendCoeffs = trend_coeff;
oo_.Smoother.integration_order = d;
......@@ -949,7 +950,7 @@ if (any(bayestopt_.pshape >0 ) & options_.mh_replic) | ...
oo_ = compute_moments_varendo(options_,M_,oo_,var_list_);
end
if options_.smoother | ~isempty(options_.filter_step_ahead) | options_.forecast
prior_posterior_statistics('posterior',data,gend);
prior_posterior_statistics('posterior',data,gend,data_index,missing_value);
end
xparam = get_posterior_parameters('mean');
set_all_parameters(xparam);
......@@ -959,7 +960,7 @@ if (~((any(bayestopt_.pshape > 0) & options_.mh_replic) | (any(bayestopt_.pshape
> 0) & options_.load_mh_file)) ...
| ~options_.smoother ) & M_.endo_nbr^2*gend < 1e7 % to be fixed
%% ML estimation, or posterior mode without metropolis-hastings or metropolis without bayesian smooth variable
[atT,innov,measurement_error,updated_variables,ys,trend_coeff,aK,T,R,P,PK,d,decomp] = DsgeSmoother(xparam1,gend,data);
[atT,innov,measurement_error,updated_variables,ys,trend_coeff,aK,T,R,P,PK,d,decomp] = DsgeSmoother(xparam1,gend,data,data_index,missing_value);
oo_.Smoother.SteadyState = ys;
oo_.Smoother.TrendCoeffs = trend_coeff;
oo_.Smoother.integration_order = d;
......
......@@ -47,7 +47,8 @@ function initial_estimation_checks(xparam1,gend,data,data_index,number_of_observ
'observed variables'])
end
if (number_of_observations==gend*nv)% No missing observations...
r = rank(data);
k = find(all(~isnan(data),2));
r = rank(data(unique(k),:));
if r < nv
error(['Estimation can''t take place because the data are perfectly' ...
' correlated']);
......
......@@ -43,7 +43,7 @@ function [LIK, lik] = missing_observations_diffuse_kalman_filter(T,R,Q,H,Pinf,Ps
%
% You should have received a copy of the GNU General Public License
% along with Dynare. If not, see <http://www.gnu.org/licenses/>.
[pp,smpl] = size(Y);
mm = size(T,2);
a = zeros(mm,1);
......@@ -59,7 +59,7 @@ function [LIK, lik] = missing_observations_diffuse_kalman_filter(T,R,Q,H,Pinf,Ps
while rank(Pinf,kalman_tol) && (t<smpl)
t = t+1;
if isempty(data_index{t},:)
if isempty(data_index{t})
a = T*a;
Pstar = T*Pstar*transpose(T)+QQ;
Pinf = T*Pinf*transpose(T);
......
......@@ -44,8 +44,9 @@ function [LIK, lik] = univariate_diffuse_kalman_filter(T,R,Q,H,Pinf,Pstar,Y,star
%
% You should have received a copy of the GNU General Public License
% along with Dynare. If not, see <http://www.gnu.org/licenses/>.
global options_
[pp ,smpl] = size(Y,1);
[pp ,smpl] = size(Y);
mm = size(T,1);
a = zeros(mm,1);
QQ = R*Q*transpose(R);
......@@ -135,7 +136,7 @@ while notsteady && (t<smpl)
a = T*a;
Pstar = T*Pstar*T' + QQ;
if t>no_more_missing_observations
notsteady = max(max(abs(P-oldP)))>riccati_tol;
notsteady = max(max(abs(Pstar-oldP)))>riccati_tol;
end
end
......
function prior_posterior_statistics(type,Y,gend)
function prior_posterior_statistics(type,Y,gend,data_index,missing_value)
% function PosteriorFilterSmootherAndForecast(Y,gend, type)
% Computes posterior filter smoother and forecasts
%
% INPUTS
% type: posterior
% prior
% gsa
% Y: data
% gend: number of observations
% type: posterior
% prior
% gsa
% Y: data
% gend: number of observations
% data_index [cell] 1*smpl cell of column vectors of indices.
% missing_value 1 if missing values, 0 otherwise
%
% OUTPUTS
% none
......@@ -152,7 +154,7 @@ for b=1:B
if run_smoother
[alphahat,etahat,epsilonhat,alphatilde,SteadyState,trend_coeff,aK] = ...
DsgeSmoother(deep,gend,Y);
DsgeSmoother(deep,gend,Y,data_index,missing_value);
if options_.loglinear
stock_smooth(dr.order_var,:,irun(1)) = alphahat(1:endo_nbr,:)+ ...
repmat(log(dr.ys(dr.order_var)),1,gend);
......
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment