Description of DsgeSmoother

0001 function [alphahat,etahat,epsilonhat,ahat,SteadyState,trend_coeff,aK,T,R,P,PK,decomp] = DsgeSmoother(xparam1,gend,Y,data_index,missing_value)
0002 % Estimation of the smoothed variables and innovations.
0003 %
0004 % INPUTS
0005 %   o xparam1       [double]   (p*1) vector of (estimated) parameters.
0006 %   o gend          [integer]  scalar specifying the number of observations ==> varargin{1}.
0007 %   o data          [double]   (T*n) matrix of data.
0008 %   o data_index    [cell]      1*smpl cell of column vectors of indices.
0009 %   o missing_value 1 if missing values, 0 otherwise
0010 %
0011 % OUTPUTS
0012 %   o alphahat      [double]  (m*T) matrix, smoothed endogenous variables.
0013 %   o etahat        [double]  (r*T) matrix, smoothed structural shocks (r>n is the umber of shocks).
0014 %   o epsilonhat    [double]  (n*T) matrix, smoothed measurement errors.
0015 %   o ahat          [double]  (m*T) matrix, one step ahead filtered (endogenous) variables.
0016 %   o SteadyState   [double]  (m*1) vector specifying the steady state level of each endogenous variable.
0017 %   o trend_coeff   [double]  (n*1) vector, parameters specifying the slope of the trend associated to each observed variable.
0018 %   o aK            [double]  (K,n,T+K) array, k (k=1,...,K) steps ahead filtered (endogenous) variables.
0019 %   o T and R       [double]  Matrices defining the state equation (T is the (m*m) transition matrix).
0020 %    P:             3D array of one-step ahead forecast error variance
0021 %                   matrices
0022 %    PK:            4D array of k-step ahead forecast error variance
0023 %                   matrices (meaningless for periods 1:d)
0024 %
0025 % ALGORITHM
0026 %   Diffuse Kalman filter (Durbin and Koopman)
0027 %
0028 % SPECIAL REQUIREMENTS
0029 %   None
0030 
0031 % Copyright (C) 2006-2011 Dynare Team
0032 %
0033 % This file is part of Dynare.
0034 %
0035 % Dynare is free software: you can redistribute it and/or modify
0036 % it under the terms of the GNU General Public License as published by
0037 % the Free Software Foundation, either version 3 of the License, or
0038 % (at your option) any later version.
0039 %
0040 % Dynare is distributed in the hope that it will be useful,
0041 % but WITHOUT ANY WARRANTY; without even the implied warranty of
0042 % MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
0043 % GNU General Public License for more details.
0044 %
0045 % You should have received a copy of the GNU General Public License
0046 % along with Dynare.  If not, see <http://www.gnu.org/licenses/>.
0047 
0048 global bayestopt_ M_ oo_ estim_params_ options_
0049 
0050 alphahat        = [];
0051 etahat  = [];
0052 epsilonhat      = [];
0053 ahat          = [];
0054 SteadyState   = [];
0055 trend_coeff   = [];
0056 aK            = [];
0057 T             = [];
0058 R             = [];
0059 P             = [];
0060 PK            = [];
0061 decomp        = [];
0062 nobs            = size(options_.varobs,1);
0063 smpl          = size(Y,2);
0064 
0065 set_all_parameters(xparam1);
0066 
0067 %------------------------------------------------------------------------------
0068 % 2. call model setup & reduction program
0069 %------------------------------------------------------------------------------
0070 oo_.dr.restrict_var_list = bayestopt_.smoother_var_list;
0071 oo_.dr.restrict_columns = bayestopt_.smoother_restrict_columns;
0072 [T,R,SteadyState,info,M_,options_,oo_] = dynare_resolve(M_,options_,oo_);
0073 bayestopt_.mf = bayestopt_.smoother_mf;
0074 if options_.noconstant
0075     constant = zeros(nobs,1);
0076 else
0077     if options_.loglinear == 1
0078         constant = log(SteadyState(bayestopt_.mfys));
0079     else
0080         constant = SteadyState(bayestopt_.mfys);
0081     end
0082 end
0083 trend_coeff = zeros(nobs,1);
0084 if bayestopt_.with_trend == 1
0085     trend_coeff = zeros(nobs,1);
0086     t = options_.trend_coeffs;
0087     for i=1:length(t)
0088         if ~isempty(t{i})
0089             trend_coeff(i) = evalin('base',t{i});
0090         end
0091     end
0092     trend = constant*ones(1,gend)+trend_coeff*(1:gend);
0093 else
0094     trend = constant*ones(1,gend);
0095 end
0096 start = options_.presample+1;
0097 np    = size(T,1);
0098 mf    = bayestopt_.smoother_mf;
0099 % ------------------------------------------------------------------------------
0100 %  3. Initial condition of the Kalman filter
0101 % ------------------------------------------------------------------------------
0102 %
0103 %  C'est ici qu'il faut d�erminer Pinf et Pstar. Si le mod�e est stationnaire,
0104 %  alors il suffit de poser Pstar comme la solution de l'�ation de Lyapounov et
0105 %  Pinf=[].
0106 %
0107 Q = M_.Sigma_e;
0108 H = M_.H;
0109 
0110 if isequal(H,0)
0111     H = zeros(nobs,nobs);
0112 end
0113 
0114 kalman_algo = options_.kalman_algo;
0115 if options_.lik_init == 1               % Kalman filter
0116     if kalman_algo ~= 2
0117         kalman_algo = 1;
0118     end
0119     Pstar = lyapunov_symm(T,R*Q*transpose(R),options_.qz_criterium,options_.lyapunov_complex_threshold);
0120     Pinf        = [];
0121 elseif options_.lik_init == 2 % Old Diffuse Kalman filter
0122     if kalman_algo ~= 2
0123         kalman_algo = 1;
0124     end
0125     Pstar = options_.Harvey_scale_factor*eye(np);
0126     Pinf        = [];
0127 elseif options_.lik_init == 3 % Diffuse Kalman filter
0128     if kalman_algo ~= 4
0129         kalman_algo = 3;
0130     end
0131     [Z,ST,R1,QT,Pstar,Pinf] = schur_statespace_transformation(mf,T,R,Q,options_.qz_criterium);
0132 elseif options_.lik_init == 4
0133     % Start from the solution of the Riccati equation.
0134     [err, Pstar] = kalman_steady_state(transpose(T),R*Q*transpose(R),transpose(build_selection_matrix(mf,np,nobs)),H);
0135     mexErrCheck('kalman_steady_state',err);
0136     Pinf  = [];
0137     if kalman_algo~=2
0138         kalman_algo = 1;
0139     end
0140 end
0141 kalman_tol = options_.kalman_tol;
0142 riccati_tol = options_.riccati_tol;
0143 data1 = Y-trend;
0144 % -----------------------------------------------------------------------------
0145 %  4. Kalman smoother
0146 % -----------------------------------------------------------------------------
0147 
0148 if ~missing_value
0149     for i=1:smpl
0150         data_index{i}=(1:nobs)';
0151     end
0152 end
0153 
0154 if kalman_algo == 1 || kalman_algo == 2
0155     ST = T;
0156     R1 = R;
0157     Z = zeros(nobs,size(T,2));
0158     for i=1:nobs
0159         Z(i,mf(i)) = 1;
0160     end
0161 end
0162 
0163 if kalman_algo == 1 || kalman_algo == 3
0164     [alphahat,epsilonhat,etahat,ahat,P,aK,PK,decomp] = missing_DiffuseKalmanSmootherH1_Z(ST, ...
0165                                                       Z,R1,Q,H,Pinf,Pstar, ...
0166                                                       data1,nobs,np,smpl,data_index, ...
0167                                                       options_.nk,kalman_tol,options_.filter_decomposition);
0168     if isinf(alphahat)
0169         if kalman_algo == 1
0170             kalman_algo = 2;
0171         elseif kalman_algo == 3
0172             kalman_algo = 4;
0173         else
0174             error('This case shouldn''t happen')
0175         end
0176     end
0177 end
0178 
0179 if kalman_algo == 2 || kalman_algo == 4
0180     if estim_params_.ncn
0181         ST = [ zeros(nobs,nobs) Z; zeros(np,nobs) T];
0182         ns = size(Q,1);
0183         R1 = [ eye(nobs) zeros(nobs, ns); zeros(np,nobs) R];
0184         Q = [H zeros(nobs,ns); zeros(ns,nobs) Q]; 
0185         Z = [eye(nobs) zeros(nobs, np)];
0186         if kalman_algo == 4
0187             [Z,ST,R1,QT,Pstar,Pinf] = schur_statespace_transformation((1:nobs)',ST,R1,Q,options_.qz_criterium);
0188         end
0189         
0190     end
0191     [alphahat,epsilonhat,etahat,ahat,P,aK,PK,decomp] = missing_DiffuseKalmanSmootherH3_Z(ST, ...
0192                                                       Z,R1,Q,diag(H), ...
0193                                                       Pinf,Pstar,data1,nobs,np,smpl,data_index, ...
0194                                                       options_.nk,kalman_tol,...
0195                                                       options_.filter_decomposition);
0196 end
0197 
0198 if kalman_algo == 3 || kalman_algo == 4
0199     alphahat = QT*alphahat;
0200     ahat = QT*ahat;
0201     nk = options_.nk;
0202     for jnk=1:nk
0203         aK(jnk,:,:) = QT*dynare_squeeze(aK(jnk,:,:));
0204         for i=1:size(PK,4)
0205             PK(jnk,:,:,i) = QT*dynare_squeeze(PK(jnk,:,:,i))*QT';
0206         end
0207         if options_.filter_decomposition
0208             for i=1:size(decomp,4)
0209                 decomp(jnk,:,:,i) = QT*dynare_squeeze(decomp(jnk,:,:,i));
0210             end
0211         end
0212     end
0213     for i=1:size(P,4)
0214         P(:,:,i) = QT*dynare_squeeze(P(:,:,i))*QT';
0215     end
0216 end
0217 
0218 if estim_params_.ncn && (kalman_algo == 2 || kalman_algo == 4)
0219     % extracting measurement errors
0220     % removing observed variables from the state vector
0221     k = nobs+(1:np);
0222     alphahat = alphahat(k,:);
0223     ahat = ahat(k,:);
0224     aK = aK(:,k,:,:);
0225     if ~isempty(PK)
0226         PK = PK(:,k,k,:);
0227     end
0228     if ~isempty(decomp)
0229         decomp = decomp(:,k,:,:);
0230     end
0231     if ~isempty(P)
0232         P = P(k,k,:);
0233     end
0234 end
DsgeSmoother

PURPOSE

SYNOPSIS

DESCRIPTION

CROSS-REFERENCE INFORMATION

SOURCE CODE
