From 706651641f75e6ac745cfebb1dc015189967f4a9 Mon Sep 17 00:00:00 2001 From: michel Date: Mon, 1 Dec 2008 12:50:53 +0000 Subject: [PATCH] adding missing data smoother + small corrections git-svn-id: https://www.dynare.org/svn/dynare/dynare_v4@2278 ac1d8469-bf42-47a9-8791-bf33cf982152 --- matlab/DsgeSmoother.m | 16 +- matlab/describe_missing_data.m | 2 +- matlab/missing_DiffuseKalmanSmoother1.m | 207 ++++++++++++++ matlab/missing_DiffuseKalmanSmoother1_Z.m | 238 ++++++++++++++++ matlab/missing_DiffuseKalmanSmoother3.m | 284 +++++++++++++++++++ matlab/missing_DiffuseKalmanSmoother3_Z.m | 316 ++++++++++++++++++++++ 6 files changed, 1060 insertions(+), 3 deletions(-) create mode 100644 matlab/missing_DiffuseKalmanSmoother1.m create mode 100644 matlab/missing_DiffuseKalmanSmoother1_Z.m create mode 100644 matlab/missing_DiffuseKalmanSmoother3.m create mode 100644 matlab/missing_DiffuseKalmanSmoother3_Z.m diff --git a/matlab/DsgeSmoother.m b/matlab/DsgeSmoother.m index 826fcc4d6..eab3545f6 100644 --- a/matlab/DsgeSmoother.m +++ b/matlab/DsgeSmoother.m @@ -264,13 +264,25 @@ function [alphahat,etahat,epsilonhat,ahat,SteadyState,trend_coeff,aK,T,R,P,PK,d, end else if kalman_algo == 1 - [alphahat,etahat,ahat,aK] = DiffuseKalmanSmoother1(T,R,Q,Pinf,Pstar,Y,trend,nobs,np,smpl,mf); + if missing_value + [alphahat,etahat,ahat,aK] = missing_DiffuseKalmanSmoother1(T,R,Q, ... + Pinf,Pstar,Y,trend,nobs,np,smpl,mf,data_index); + else + [alphahat,etahat,ahat,aK] = DiffuseKalmanSmoother1(T,R,Q, ... + Pinf,Pstar,Y,trend,nobs,np,smpl,mf); + end if all(alphahat(:)==0) kalman_algo = 2; end end if kalman_algo == 2 - [alphahat,etahat,ahat,aK] = DiffuseKalmanSmoother3(T,R,Q,Pinf,Pstar,Y,trend,nobs,np,smpl,mf); + if missing_value + [alphahat,etahat,ahat,aK] = missing_DiffuseKalmanSmoother3(T,R,Q, ... + Pinf,Pstar,Y,trend,nobs,np,smpl,mf,data_index); + else + [alphahat,etahat,ahat,aK] = DiffuseKalmanSmoother3(T,R,Q, ... + Pinf,Pstar,Y,trend,nobs,np,smpl,mf); + end end if kalman_algo == 3 data1 = Y - trend; diff --git a/matlab/describe_missing_data.m b/matlab/describe_missing_data.m index b311acd80..97d16a9f5 100644 --- a/matlab/describe_missing_data.m +++ b/matlab/describe_missing_data.m @@ -17,6 +17,6 @@ number_of_observations = length(variable_index); if ~missing_observations_counter no_more_missing_observations = 0; else - tmp = find(missing_observations_counter>=(gend*nvarobs-number_of_observations)) + tmp = find(missing_observations_counter>=(gend*nvarobs-number_of_observations)); no_more_missing_observations = tmp(1); end \ No newline at end of file diff --git a/matlab/missing_DiffuseKalmanSmoother1.m b/matlab/missing_DiffuseKalmanSmoother1.m new file mode 100644 index 000000000..d86830f52 --- /dev/null +++ b/matlab/missing_DiffuseKalmanSmoother1.m @@ -0,0 +1,207 @@ +function [alphahat,etahat,atilde, aK] = DiffuseKalmanSmoother1(T,R,Q,Pinf1,Pstar1,Y,trend,pp,mm,smpl,mf,data_index) + +% function [alphahat,etahat,a, aK] = DiffuseKalmanSmoother1(T,R,Q,Pinf1,Pstar1,Y,trend,pp,mm,smpl,mf) +% Computes the diffuse kalman smoother without measurement error, in the case of a non-singular var-cov matrix +% +% INPUTS +% T: mm*mm matrix +% R: mm*rr matrix +% Q: rr*rr matrix +% Pinf1: mm*mm diagonal matrix with with q ones and m-q zeros +% Pstar1: mm*mm variance-covariance matrix with stationary variables +% Y: pp*1 vector +% trend +% pp: number of observed variables +% mm: number of state variables +% smpl: sample size +% mf: observed variables index in the state vector +% data_index [cell] 1*smpl cell of column vectors of indices. +% +% OUTPUTS +% alphahat: smoothed state variables (a_{t|T}) +% etahat: smoothed shocks +% atilde: matrix of updated variables (a_{t|t}) +% aK: 3D array of k step ahead filtered state variables (a_{t+k|t}) + +% SPECIAL REQUIREMENTS +% See "Filtering and Smoothing of State Vector for Diffuse State Space +% Models", S.J. Koopman and J. Durbin (2003, in Journal of Time Series +% Analysis, vol. 24(1), pp. 85-98). + +% Copyright (C) 2004-2008 Dynare Team +% +% This file is part of Dynare. +% +% Dynare is free software: you can redistribute it and/or modify +% it under the terms of the GNU General Public License as published by +% the Free Software Foundation, either version 3 of the License, or +% (at your option) any later version. +% +% Dynare is distributed in the hope that it will be useful, +% but WITHOUT ANY WARRANTY; without even the implied warranty of +% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +% GNU General Public License for more details. +% +% You should have received a copy of the GNU General Public License +% along with Dynare. If not, see . + +% modified by M. Ratto: +% new output argument aK (1-step to k-step predictions) +% new options_.nk: the max step ahed prediction in aK (default is 4) +% new crit1 value for rank of Pinf +% it is assured that P is symmetric + + +global options_ + +nk = options_.nk; +spinf = size(Pinf1); +spstar = size(Pstar1); +v = zeros(pp,smpl); +a = zeros(mm,smpl+1); +atilde = zeros(mm,smpl); +aK = zeros(nk,mm,smpl+1); +iF = zeros(pp,pp,smpl); +Fstar = zeros(pp,pp,smpl); +iFinf = zeros(pp,pp,smpl); +K = zeros(mm,pp,smpl); +L = zeros(mm,mm,smpl); +Linf = zeros(mm,mm,smpl); +Kstar = zeros(mm,pp,smpl); +P = zeros(mm,mm,smpl+1); +Pstar = zeros(spstar(1),spstar(2),smpl+1); Pstar(:,:,1) = Pstar1; +Pinf = zeros(spinf(1),spinf(2),smpl+1); Pinf(:,:,1) = Pinf1; +crit = options_.kalman_tol; +crit1 = 1.e-8; +steady = smpl; +rr = size(Q,1); +QQ = R*Q*transpose(R); +QRt = Q*transpose(R); +alphahat = zeros(mm,smpl); +etahat = zeros(rr,smpl); +r = zeros(mm,smpl+1); + +Z = zeros(pp,mm); +for i=1:pp; + Z(i,mf(i)) = 1; +end + +t = 0; +while rank(Pinf(:,:,t+1),crit1) & td+1 + t = t-1; + di = data_index{t}; + if isempty(di) + r(:,t) = L(:,:,t)'*r(:,t+1); + else + r(:,t) = Z(di,:)'*iF(di,di,t)*v(di,t) + L(:,:,t)'*r(:,t+1); + end + alphahat(:,t) = a(:,t) + P(:,:,t)*r(:,t); + etahat(:,t) = QRt*r(:,t); +end +if d + r0 = zeros(mm,d+1); + r0(:,d+1) = r(:,d+1); + r1 = zeros(mm,d+1); + for t = d:-1:1 + r0(:,t) = Linf(:,:,t)'*r0(:,t+1); + di = data_index{t}; + if isempty(di) + r1(:,t) = Linf(:,:,t)'*r1(:,t+1); + else + r1(:,t) = Z(di,:)'*(iFinf(di,di,t)*v(di,t)-Kstar(:,di,t)'*r0(:,t+1)) ... + + Linf(:,:,t)'*r1(:,t+1); + end + alphahat(:,t) = a(:,t) + Pstar(:,:,t)*r0(:,t) + Pinf(:,:,t)*r1(:,t); + etahat(:,t) = QRt*r0(:,t); + end +end diff --git a/matlab/missing_DiffuseKalmanSmoother1_Z.m b/matlab/missing_DiffuseKalmanSmoother1_Z.m new file mode 100644 index 000000000..60edcd8f0 --- /dev/null +++ b/matlab/missing_DiffuseKalmanSmoother1_Z.m @@ -0,0 +1,238 @@ +function [alphahat,etahat,atilde,P,aK,PK,d,decomp] = DiffuseKalmanSmoother1_Z(T,Z,R,Q,Pinf1,Pstar1,Y,pp,mm,smpl,data_index) + +% function [alphahat,etahat,a, aK] = DiffuseKalmanSmoother1(T,Z,R,Q,Pinf1,Pstar1,Y,pp,mm,smpl) +% Computes the diffuse kalman smoother without measurement error, in the case of a non-singular var-cov matrix +% +% INPUTS +% T: mm*mm matrix +% Z: pp*mm matrix +% R: mm*rr matrix +% Q: rr*rr matrix +% Pinf1: mm*mm diagonal matrix with with q ones and m-q zeros +% Pstar1: mm*mm variance-covariance matrix with stationary variables +% Y: pp*1 vector +% pp: number of observed variables +% mm: number of state variables +% smpl: sample size +% data_index [cell] 1*smpl cell of column vectors of indices. +% +% OUTPUTS +% alphahat: smoothed variables (a_{t|T}) +% etahat: smoothed shocks +% atilde: matrix of updated variables (a_{t|t}) +% aK: 3D array of k step ahead filtered state variables (a_{t+k|t) +% (meaningless for periods 1:d) +% P: 3D array of one-step ahead forecast error variance +% matrices +% PK: 4D array of k-step ahead forecast error variance +% matrices (meaningless for periods 1:d) +% d: number of periods where filter remains in diffuse part +% (should be equal to the order of integration of the model) +% decomp: decomposition of the effect of shocks on filtered values +% +% SPECIAL REQUIREMENTS +% See "Filtering and Smoothing of State Vector for Diffuse State Space +% Models", S.J. Koopman and J. Durbin (2003, in Journal of Time Series +% Analysis, vol. 24(1), pp. 85-98). + +% Copyright (C) 2004-2008 Dynare Team +% +% This file is part of Dynare. +% +% Dynare is free software: you can redistribute it and/or modify +% it under the terms of the GNU General Public License as published by +% the Free Software Foundation, either version 3 of the License, or +% (at your option) any later version. +% +% Dynare is distributed in the hope that it will be useful, +% but WITHOUT ANY WARRANTY; without even the implied warranty of +% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +% GNU General Public License for more details. +% +% You should have received a copy of the GNU General Public License +% along with Dynare. If not, see . + +% modified by M. Ratto: +% new output argument aK (1-step to k-step predictions) +% new options_.nk: the max step ahed prediction in aK (default is 4) +% new crit1 value for rank of Pinf +% it is assured that P is symmetric + + +global options_ + +d = 0; +decomp = []; +nk = options_.nk; +spinf = size(Pinf1); +spstar = size(Pstar1); +v = zeros(pp,smpl); +a = zeros(mm,smpl+1); +atilde = zeros(mm,smpl); +aK = zeros(nk,mm,smpl+nk); +PK = zeros(nk,mm,mm,smpl+nk); +iF = zeros(pp,pp,smpl); +Fstar = zeros(pp,pp,smpl); +iFinf = zeros(pp,pp,smpl); +K = zeros(mm,pp,smpl); +L = zeros(mm,mm,smpl); +Linf = zeros(mm,mm,smpl); +Kstar = zeros(mm,pp,smpl); +P = zeros(mm,mm,smpl+1); +Pstar = zeros(spstar(1),spstar(2),smpl+1); Pstar(:,:,1) = Pstar1; +Pinf = zeros(spinf(1),spinf(2),smpl+1); Pinf(:,:,1) = Pinf1; +crit = options_.kalman_tol; +crit1 = 1.e-8; +steady = smpl; +rr = size(Q,1); +QQ = R*Q*transpose(R); +QRt = Q*transpose(R); +alphahat = zeros(mm,smpl); +etahat = zeros(rr,smpl); +r = zeros(mm,smpl+1); + +t = 0; +while rank(Pinf(:,:,t+1),crit1) & td+1 + t = t-1; + di = data_index{t}; + if isempty(di) + r(:,t) = L(:,:,t)'*r(:,t+1); + else + ZZ = Z(di,:); + r(:,t) = ZZ'*iF(di,di,t)*v(di,t) + L(:,:,t)'*r(:,t+1); + end + alphahat(:,t) = a(:,t) + P(:,:,t)*r(:,t); + etahat(:,t) = QRt*r(:,t); +end +if d + r0 = zeros(mm,d+1); + r0(:,d+1) = r(:,d+1); + r1 = zeros(mm,d+1); + for t = d:-1:1 + r0(:,t) = Linf(:,:,t)'*r0(:,t+1); + di = data_index{t}; + if isempty(di) + r1(:,t) = Linf(:,:,t)'*r1(:,t+1); + else + r1(:,t) = Z(di,:)'*(iFinf(di,di,t)*v(di,t)-Kstar(:,di,t)'*r0(:,t+1)) ... + + Linf(:,:,t)'*r1(:,t+1); + end + alphahat(:,t) = a(:,t) + Pstar(:,:,t)*r0(:,t) + Pinf(:,:,t)*r1(:,t); + etahat(:,t) = QRt*r0(:,t); + end +end + +if nargout > 7 + decomp = zeros(nk,mm,rr,smpl+nk); + ZRQinv = inv(Z*QQ*Z'); + for t = max(d,1):smpl + ri_d = Z'*iF(:,:,t)*v(:,t); + + % calculate eta_tm1t + eta_tm1t = QRt*ri_d; + % calculate decomposition + Ttok = eye(mm,mm); + for h = 1:nk + for j=1:rr + eta=zeros(rr,1); + eta(j) = eta_tm1t(j); + decomp(h,:,j,t+h) = T^(h-1)*P(:,:,t)*Z'*ZRQinv*Z*R*eta; + end + end + end +end diff --git a/matlab/missing_DiffuseKalmanSmoother3.m b/matlab/missing_DiffuseKalmanSmoother3.m new file mode 100644 index 000000000..13efff232 --- /dev/null +++ b/matlab/missing_DiffuseKalmanSmoother3.m @@ -0,0 +1,284 @@ +function [alphahat,etahat,a, aK] = missing_DiffuseKalmanSmoother3(T,R,Q,Pinf1,Pstar1,Y,trend,pp,mm,smpl,mf,data_index) +% function [alphahat,etahat,a1, aK] = missing_DiffuseKalmanSmoother3(T,R,Q,Pinf1,Pstar1,Y,trend,pp,mm,smpl,mf,data_index) +% Computes the diffuse kalman smoother without measurement error, in the case of a singular var-cov matrix. +% Univariate treatment of multivariate time series. +% +% INPUTS +% T: mm*mm matrix +% R: mm*rr matrix +% Q: rr*rr matrix +% Pinf1: mm*mm diagonal matrix with with q ones and m-q zeros +% Pstar1: mm*mm variance-covariance matrix with stationary variables +% Y: pp*1 vector +% trend +% pp: number of observed variables +% mm: number of state variables +% smpl: sample size +% mf: observed variables index in the state vector +% data_index [cell] 1*smpl cell of column vectors of indices. +% +% OUTPUTS +% alphahat: smoothed state variables (a_{t|T}) +% etahat: smoothed shocks +% a: matrix of updated variables (a_{t|t}) +% aK: 3D array of k step ahead filtered state variables (a_{t+k|t}) +% +% SPECIAL REQUIREMENTS +% See "Filtering and Smoothing of State Vector for Diffuse State Space +% Models", S.J. Koopman and J. Durbin (2003, in Journal of Time Series +% Analysis, vol. 24(1), pp. 85-98). + +% Copyright (C) 2004-2008 Dynare Team +% +% This file is part of Dynare. +% +% Dynare is free software: you can redistribute it and/or modify +% it under the terms of the GNU General Public License as published by +% the Free Software Foundation, either version 3 of the License, or +% (at your option) any later version. +% +% Dynare is distributed in the hope that it will be useful, +% but WITHOUT ANY WARRANTY; without even the implied warranty of +% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +% GNU General Public License for more details. +% +% You should have received a copy of the GNU General Public License +% along with Dynare. If not, see . + +% Modified by M. Ratto +% New output argument aK: 1-step to nk-stpe ahed predictions) +% New input argument nk: max order of predictions in aK +% New option options_.diffuse_d where the DKF stops (common with +% diffuselikelihood3) +% New icc variable to count number of iterations for Finf steps +% Pstar % Pinf simmetric +% New termination of DKF iterations based on options_.diffuse_d +% Li also stored during DKF iterations !! +% some bugs corrected in the DKF part of the smoother (Z matrix and +% alphahat) + +global options_ + +nk = options_.nk; +spinf = size(Pinf1); +spstar = size(Pstar1); +v = zeros(pp,smpl); +a = zeros(mm,smpl); +a1 = zeros(mm,smpl+1); +aK = zeros(nk,mm,smpl+nk); + +if isempty(options_.diffuse_d), + smpl_diff = 1; +else + smpl_diff=rank(Pinf1); +end + +Fstar = zeros(pp,smpl_diff); +Finf = zeros(pp,smpl_diff); +Ki = zeros(mm,pp,smpl); +Li = zeros(mm,mm,pp,smpl); +Linf = zeros(mm,mm,pp,smpl_diff); +L0 = zeros(mm,mm,pp,smpl_diff); +Kstar = zeros(mm,pp,smpl_diff); +P = zeros(mm,mm,smpl+1); +P1 = P; +Pstar = zeros(spstar(1),spstar(2),smpl_diff+1); Pstar(:,:,1) = Pstar1; +Pinf = zeros(spinf(1),spinf(2),smpl_diff+1); Pinf(:,:,1) = Pinf1; +Pstar1 = Pstar; +Pinf1 = Pinf; +crit = options_.kalman_tol; +crit1 = 1.e-6; +steady = smpl; +rr = size(Q,1); +QQ = R*Q*transpose(R); +QRt = Q*transpose(R); +alphahat = zeros(mm,smpl); +etahat = zeros(rr,smpl); +r = zeros(mm,smpl+1); + +Z = zeros(pp,mm); +for i=1:pp; + Z(i,mf(i)) = 1; +end + +t = 0; +icc=0; +newRank = rank(Pinf(:,:,1),crit1); +while newRank & t < smpl + t = t+1; + a(:,t) = a1(:,t); + Pstar(:,:,t)=tril(Pstar(:,:,t))+transpose(tril(Pstar(:,:,t),-1)); + Pinf(:,:,t)=tril(Pinf(:,:,t))+transpose(tril(Pinf(:,:,t),-1)); + Pstar1(:,:,t) = Pstar(:,:,t); + Pinf1(:,:,t) = Pinf(:,:,t); + di = data_index{t}'; + for i=di + v(i,t) = Y(i,t)-a(mf(i),t)-trend(i,t); + Fstar(i,t) = Pstar(mf(i),mf(i),t); + Finf(i,t) = Pinf(mf(i),mf(i),t); + Kstar(:,i,t) = Pstar(:,mf(i),t); + if Finf(i,t) > crit & newRank, % original MJ: if Finf(i,t) > crit + icc=icc+1; + Kinf(:,i,t) = Pinf(:,mf(i),t); + Linf(:,:,i,t) = eye(mm) - Kinf(:,i,t)*Z(i,:)/Finf(i,t); + L0(:,:,i,t) = (Kinf(:,i,t)*Fstar(i,t)/Finf(i,t) - Kstar(:,i,t))*Z(i,:)/Finf(i,t); + a(:,t) = a(:,t) + Kinf(:,i,t)*v(i,t)/Finf(i,t); + Pstar(:,:,t) = Pstar(:,:,t) + ... + Kinf(:,i,t)*transpose(Kinf(:,i,t))*Fstar(i,t)/(Finf(i,t)*Finf(i,t)) - ... + (Kstar(:,i,t)*transpose(Kinf(:,i,t)) +... + Kinf(:,i,t)*transpose(Kstar(:,i,t)))/Finf(i,t); + Pinf(:,:,t) = Pinf(:,:,t) - Kinf(:,i,t)*transpose(Kinf(:,i,t))/Finf(i,t); + Pstar(:,:,t)=tril(Pstar(:,:,t))+transpose(tril(Pstar(:,:,t),-1)); + Pinf(:,:,t)=tril(Pinf(:,:,t))+transpose(tril(Pinf(:,:,t),-1)); + % new terminiation criteria by M. Ratto + P0=Pinf(:,:,t); + % newRank = any(diag(P0(mf,mf))>crit); + % if newRank==0, id = i; end, + if ~isempty(options_.diffuse_d), + newRank = (icccrit)==0; + if newRank & (any(diag(P0(mf,mf))>crit)==0 & rank(P0,crit1)==0); + disp('WARNING!! Change in OPTIONS_.DIFFUSE_D in univariate DKF') + options_.diffuse_d = icc; + newRank=0; + end + else + %newRank = any(diag(P0(mf,mf))>crit); + newRank = (any(diag(P0(mf,mf))>crit) | rank(P0,crit1)); + if newRank==0, + options_.diffuse_d = icc; + end + end, + % if newRank==0, + % options_.diffuse_d=i; %this is buggy + % end + % end new terminiation criteria by M. Ratto + elseif Fstar(i,t) > crit + %% Note that : (1) rank(Pinf)=0 implies that Finf = 0, (2) outside this loop (when for some i and t the condition + %% rank(Pinf)=0 is satisfied we have P = Pstar and F = Fstar and (3) Finf = 0 does not imply that + %% rank(Pinf)=0. [stéphane,11-03-2004]. + Li(:,:,i,t) = eye(mm)-Kstar(:,i,t)*Z(i,:)/Fstar(i,t); % we need to store Li for DKF smoother + a(:,t) = a(:,t) + Kstar(:,i,t)*v(i,t)/Fstar(i,t); + Pstar(:,:,t) = Pstar(:,:,t) - Kstar(:,i,t)*transpose(Kstar(:,i,t))/Fstar(i,t); + Pstar(:,:,t)=tril(Pstar(:,:,t))+transpose(tril(Pstar(:,:,t),-1)); + end + end + a1(:,t+1) = T*a(:,t); + for jnk=1:nk, + aK(jnk,:,t+jnk) = T^jnk*a(:,t); + end + Pstar(:,:,t+1) = T*Pstar(:,:,t)*transpose(T)+ QQ; + Pinf(:,:,t+1) = T*Pinf(:,:,t)*transpose(T); + P0=Pinf(:,:,t+1); + if newRank, + %newRank = any(diag(P0(mf,mf))>crit); + newRank = rank(P0,crit1); + end +end + + +d = t; +P(:,:,d+1) = Pstar(:,:,d+1); +Linf = Linf(:,:,:,1:d); +L0 = L0(:,:,:,1:d); +Fstar = Fstar(:,1:d); +Finf = Finf(:,1:d); +Kstar = Kstar(:,:,1:d); +Pstar = Pstar(:,:,1:d); +Pinf = Pinf(:,:,1:d); +Pstar1 = Pstar1(:,:,1:d); +Pinf1 = Pinf1(:,:,1:d); +notsteady = 1; +while notsteady & t crit + Li(:,:,i,t) = eye(mm)-Ki(:,i,t)*Z(i,:)/Fi(i,t); + a(:,t) = a(:,t) + Ki(:,i,t)*v(i,t)/Fi(i,t); + P(:,:,t) = P(:,:,t) - Ki(:,i,t)*transpose(Ki(:,i,t))/Fi(i,t); + P(:,:,t)=tril(P(:,:,t))+transpose(tril(P(:,:,t),-1)); + end + end + a1(:,t+1) = T*a(:,t); + for jnk=1:nk, + aK(jnk,:,t+jnk) = T^jnk*a(:,t); + end + P(:,:,t+1) = T*P(:,:,t)*transpose(T) + QQ; + notsteady = ~(max(max(abs(P(:,:,t+1)-P(:,:,t)))) crit + a(:,t) = a(:,t) + Ki_s(:,i)*v(i,t)/Fi_s(i); + end + end + a1(:,t+1) = T*a(:,t); + for jnk=1:nk, + aK(jnk,:,t+jnk) = T^jnk*a(:,t); + end +end +ri=zeros(mm,1); +t = smpl+1; +while t>d+1 + t = t-1; + di = flipud(data_index{t})'; + for i = di + if Fi(i,t) > crit + ri = Z(i,:)'/Fi(i,t)*v(i,t)+Li(:,:,i,t)'*ri; + end + end + r(:,t) = ri; + alphahat(:,t) = a1(:,t) + P1(:,:,t)*r(:,t); + etahat(:,t) = QRt*r(:,t); + ri = T'*ri; +end +if d + r0 = zeros(mm,d); + r0(:,d) = ri; + r1 = zeros(mm,d); + for t = d:-1:1 + di = flipud(data_index{t})'; + for i = di + if Finf(i,t) > crit & ~(t==d & i>options_.diffuse_d), + % use of options_.diffuse_d to be sure of DKF termination + %r1(:,t) = transpose(Z)*v(:,t)/Finf(i,t) + ... BUG HERE in transpose(Z) + r1(:,t) = Z(i,:)'*v(i,t)/Finf(i,t) + ... + L0(:,:,i,t)'*r0(:,t) + Linf(:,:,i,t)'*r1(:,t); + r0(:,t) = Linf(:,:,i,t)'*r0(:,t); + elseif Fstar(i,t) > crit % step needed whe Finf == 0 + r0(:,t) = Z(i,:)'/Fstar(i,t)*v(i,t)+Li(:,:,i,t)'*r0(:,t); + end + end + alphahat(:,t) = a1(:,t) + Pstar1(:,:,t)*r0(:,t) + Pinf1(:,:,t)*r1(:,t); + r(:,t) = r0(:,t); + etahat(:,t) = QRt*r(:,t); + if t > 1 + r0(:,t-1) = T'*r0(:,t); + r1(:,t-1) = T'*r1(:,t); + end + end +end + + diff --git a/matlab/missing_DiffuseKalmanSmoother3_Z.m b/matlab/missing_DiffuseKalmanSmoother3_Z.m new file mode 100644 index 000000000..e7c74a7a4 --- /dev/null +++ b/matlab/missing_DiffuseKalmanSmoother3_Z.m @@ -0,0 +1,316 @@ +function [alphahat,etahat,a,P,aK,PK,d,decomp] = missing_DiffuseKalmanSmoother3_Z(T,Z,R,Q,Pinf1,Pstar1,Y,pp,mm,smpl,data_index) +% function [alphahat,etahat,a1,P,aK,PK,d,decomp_filt] = missing_DiffuseKalmanSmoother3_Z(T,Z,R,Q,Pinf1,Pstar1,Y,pp,mm,smpl) +% Computes the diffuse kalman smoother without measurement error, in the case of a singular var-cov matrix. +% Univariate treatment of multivariate time series. +% +% INPUTS +% T: mm*mm matrix +% Z: pp*mm matrix +% R: mm*rr matrix +% Q: rr*rr matrix +% Pinf1: mm*mm diagonal matrix with with q ones and m-q zeros +% Pstar1: mm*mm variance-covariance matrix with stationary variables +% Y: pp*1 vector +% pp: number of observed variables +% mm: number of state variables +% smpl: sample size +% data_index [cell] 1*smpl cell of column vectors of indices. +% +% OUTPUTS +% alphahat: smoothed state variables (a_{t|T}) +% etahat: smoothed shocks +% a: matrix of updated variables (a_{t|t}) +% aK: 3D array of k step ahead filtered state variables (a_{t+k|t}) +% (meaningless for periods 1:d) +% P: 3D array of one-step ahead forecast error variance +% matrices +% PK: 4D array of k-step ahead forecast error variance +% matrices (meaningless for periods 1:d) +% d: number of periods where filter remains in diffuse part +% (should be equal to the order of integration of the model) +% decomp: decomposition of the effect of shocks on filtered values +% +% SPECIAL REQUIREMENTS +% See "Filtering and Smoothing of State Vector for Diffuse State Space +% Models", S.J. Koopman and J. Durbin (2003, in Journal of Time Series +% Analysis, vol. 24(1), pp. 85-98). + +% Copyright (C) 2004-2008 Dynare Team +% +% This file is part of Dynare. +% +% Dynare is free software: you can redistribute it and/or modify +% it under the terms of the GNU General Public License as published by +% the Free Software Foundation, either version 3 of the License, or +% (at your option) any later version. +% +% Dynare is distributed in the hope that it will be useful, +% but WITHOUT ANY WARRANTY; without even the implied warranty of +% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +% GNU General Public License for more details. +% +% You should have received a copy of the GNU General Public License +% along with Dynare. If not, see . + +% Modified by M. Ratto +% New output argument aK: 1-step to nk-stpe ahed predictions) +% New input argument nk: max order of predictions in aK +% New option options_.diffuse_d where the DKF stops (common with +% diffuselikelihood3) +% New icc variable to count number of iterations for Finf steps +% Pstar % Pinf simmetric +% New termination of DKF iterations based on options_.diffuse_d +% Li also stored during DKF iterations !! +% some bugs corrected in the DKF part of the smoother (Z matrix and +% alphahat) + +global options_ + +d = 0; +decomp = []; +nk = options_.nk; +spinf = size(Pinf1); +spstar = size(Pstar1); +v = zeros(pp,smpl); +a = zeros(mm,smpl); +a1 = zeros(mm,smpl+1); +aK = zeros(nk,mm,smpl+nk); + +if isempty(options_.diffuse_d), + smpl_diff = 1; +else + smpl_diff=rank(Pinf1); +end + +Fstar = zeros(pp,smpl_diff); +Finf = zeros(pp,smpl_diff); +Ki = zeros(mm,pp,smpl); +Li = zeros(mm,mm,pp,smpl); +Linf = zeros(mm,mm,pp,smpl_diff); +L0 = zeros(mm,mm,pp,smpl_diff); +Kstar = zeros(mm,pp,smpl_diff); +P = zeros(mm,mm,smpl+1); +P1 = P; +aK = zeros(nk,mm,smpl+nk); +PK = zeros(nk,mm,mm,smpl+nk); +Pstar = zeros(spstar(1),spstar(2),smpl_diff+1); Pstar(:,:,1) = Pstar1; +Pinf = zeros(spinf(1),spinf(2),smpl_diff+1); Pinf(:,:,1) = Pinf1; +Pstar1 = Pstar; +Pinf1 = Pinf; +crit = options_.kalman_tol; +crit1 = 1.e-6; +steady = smpl; +rr = size(Q,1); % number of structural shocks +QQ = R*Q*transpose(R); +QRt = Q*transpose(R); +alphahat = zeros(mm,smpl); +etahat = zeros(rr,smpl); +r = zeros(mm,smpl); + +t = 0; +icc=0; +newRank = rank(Pinf(:,:,1),crit1); +while newRank & t < smpl + t = t+1; + a(:,t) = a1(:,t); + Pstar(:,:,t)=tril(Pstar(:,:,t))+tril(Pstar(:,:,t),-1)'; + Pinf(:,:,t)=tril(Pinf(:,:,t))+tril(Pinf(:,:,t),-1)'; + Pstar1(:,:,t) = Pstar(:,:,t); + Pinf1(:,:,t) = Pinf(:,:,t); + di = data_index{t}'; + for i=di + Zi = Z(i,:); + v(i,t) = Y(i,t)-Zi*a(:,t); + Fstar(i,t) = Zi*Pstar(:,:,t)*Zi'; + Finf(i,t) = Zi*Pinf(:,:,t)*Zi'; + Kstar(:,i,t) = Pstar(:,:,t)*Zi'; + if Finf(i,t) > crit & newRank + icc=icc+1; + Kinf(:,i,t) = Pinf(:,:,t)*Zi'; + Linf(:,:,i,t) = eye(mm) - Kinf(:,i,t)*Z(i,:)/Finf(i,t); + L0(:,:,i,t) = (Kinf(:,i,t)*Fstar(i,t)/Finf(i,t) - Kstar(:,i,t))*Zi/Finf(i,t); + a(:,t) = a(:,t) + Kinf(:,i,t)*v(i,t)/Finf(i,t); + Pstar(:,:,t) = Pstar(:,:,t) + ... + Kinf(:,i,t)*Kinf(:,i,t)'*Fstar(i,t)/(Finf(i,t)*Finf(i,t)) - ... + (Kstar(:,i,t)*Kinf(:,i,t)' +... + Kinf(:,i,t)*Kstar(:,i,t)')/Finf(i,t); + Pinf(:,:,t) = Pinf(:,:,t) - Kinf(:,i,t)*Kinf(:,i,t)'/Finf(i,t); + Pstar(:,:,t)=tril(Pstar(:,:,t))+tril(Pstar(:,:,t),-1)'; + Pinf(:,:,t)=tril(Pinf(:,:,t))+tril(Pinf(:,:,t),-1)'; + % new terminiation criteria by M. Ratto + P0=Pinf(:,:,t); + if ~isempty(options_.diffuse_d), + newRank = (icccrit)==0 & rank(P0,crit1)==0); + disp('WARNING!! Change in OPTIONS_.DIFFUSE_D in univariate DKF') + options_.diffuse_d = icc; + newRank=0; + end + else + newRank = (any(diag(Z*P0*Z')>crit) | rank(P0,crit1)); + if newRank==0, + options_.diffuse_d = icc; + end + end, + % end new terminiation criteria by M. Ratto + elseif Fstar(i,t) > crit + %% Note that : (1) rank(Pinf)=0 implies that Finf = 0, (2) outside this loop (when for some i and t the condition + %% rank(Pinf)=0 is satisfied we have P = Pstar and F = Fstar and (3) Finf = 0 does not imply that + %% rank(Pinf)=0. [stéphane,11-03-2004]. + Li(:,:,i,t) = eye(mm)-Kstar(:,i,t)*Z(i,:)/Fstar(i,t); % we need to store Li for DKF smoother + a(:,t) = a(:,t) + Kstar(:,i,t)*v(i,t)/Fstar(i,t); + Pstar(:,:,t) = Pstar(:,:,t) - Kstar(:,i,t)*Kstar(:,i,t)'/Fstar(i,t); + Pstar(:,:,t)=tril(Pstar(:,:,t))+tril(Pstar(:,:,t),-1)'; + end + end + a1(:,t+1) = T*a(:,t); + for jnk=1:nk, + aK(jnk,:,t+jnk) = T^jnk*a(:,t); + end + Pstar(:,:,t+1) = T*Pstar(:,:,t)*T'+ QQ; + Pinf(:,:,t+1) = T*Pinf(:,:,t)*T'; + P0=Pinf(:,:,t+1); + if newRank, + newRank = rank(P0,crit1); + end +end + + +d = t; +P(:,:,d+1) = Pstar(:,:,d+1); +Linf = Linf(:,:,:,1:d); +L0 = L0(:,:,:,1:d); +Fstar = Fstar(:,1:d); +Finf = Finf(:,1:d); +Kstar = Kstar(:,:,1:d); +Pstar = Pstar(:,:,1:d); +Pinf = Pinf(:,:,1:d); +Pstar1 = Pstar1(:,:,1:d); +Pinf1 = Pinf1(:,:,1:d); +notsteady = 1; +while notsteady & t crit + Li(:,:,i,t) = eye(mm)-Ki(:,i,t)*Z(i,:)/Fi(i,t); + a(:,t) = a(:,t) + Ki(:,i,t)*v(i,t)/Fi(i,t); + P(:,:,t) = P(:,:,t) - Ki(:,i,t)*Ki(:,i,t)'/Fi(i,t); + P(:,:,t)=tril(P(:,:,t))+tril(P(:,:,t),-1)'; + end + end + a1(:,t+1) = T*a(:,t); + Pf = P(:,:,t); + for jnk=1:nk, + Pf = T*Pf*T' + QQ; + aK(jnk,:,t+jnk) = T^jnk*a(:,t); + PK(jnk,:,:,t+jnk) = Pf; + end + P(:,:,t+1) = T*P(:,:,t)*T' + QQ; + notsteady = ~(max(max(abs(P(:,:,t+1)-P(:,:,t)))) crit + a(:,t) = a(:,t) + Ki_s(:,i)*v(i,t)/Fi_s(i); + end + end + a1(:,t+1) = T*a(:,t); + Pf = P(:,:,t); + for jnk=1:nk, + Pf = T*Pf*T' + QQ; + aK(jnk,:,t+jnk) = T^jnk*a(:,t); + PK(jnk,:,:,t+jnk) = Pf; + end +end +ri=zeros(mm,1); +t = smpl+1; +while t > d+1 + t = t-1; + di = flipud(data_index{t})'; + for i = di + if Fi(i,t) > crit + ri = Z(i,:)'/Fi(i,t)*v(i,t)+Li(:,:,i,t)'*ri; + end + end + r(:,t) = ri; + alphahat(:,t) = a1(:,t) + P1(:,:,t)*r(:,t); + etahat(:,t) = QRt*r(:,t); + ri = T'*ri; +end +if d + r0 = zeros(mm,d); + r0(:,d) = ri; + r1 = zeros(mm,d); + for t = d:-1:1 + di = flipud(data_index{t})'; + for i = di + % if Finf(i,t) > crit & ~(t==d & i>options_.diffuse_d), % use of options_.diffuse_d to be sure of DKF termination + if Finf(i,t) > crit + r1(:,t) = Z(i,:)'*v(i,t)/Finf(i,t) + ... + L0(:,:,i,t)'*r0(:,t) + Linf(:,:,i,t)'*r1(:,t); + r0(:,t) = Linf(:,:,i,t)'*r0(:,t); + elseif Fstar(i,t) > crit % step needed whe Finf == 0 + r0(:,t) = Z(i,:)'/Fstar(i,t)*v(i,t)+Li(:,:,i,t)'*r0(:,t); + end + end + alphahat(:,t) = a1(:,t) + Pstar1(:,:,t)*r0(:,t) + Pinf1(:,:,t)*r1(:,t); + r(:,t) = r0(:,t); + etahat(:,t) = QRt*r(:,t); + if t > 1 + r0(:,t-1) = T'*r0(:,t); + r1(:,t-1) = T'*r1(:,t); + end + end +end + +if nargout > 7 + decomp = zeros(nk,mm,rr,smpl+nk); + ZRQinv = inv(Z*QQ*Z'); + for t = max(d,1):smpl + ri_d = zeros(mm,1); + di = flipud(data_index{t})'; + for i = di + if Fi(i,t) > crit + ri_d = Z(i,:)'/Fi(i,t)*v(i,t)+Li(:,:,i,t)'*ri_d; + end + end + + % calculate eta_tm1t + eta_tm1t = QRt*ri_d; + % calculate decomposition + Ttok = eye(mm,mm); + for h = 1:nk + for j=1:rr + eta=zeros(rr,1); + eta(j) = eta_tm1t(j); + decomp(h,:,j,t+h) = Ttok*P1(:,:,t)*Z'*ZRQinv*Z*R*eta; + end + Ttok = T*Ttok; + end + end +end