port to occbin smoother the same computational improvements done for standard one under smother_redux option. This also require to provide occbin reduced state-space matrices as output argument of missing_DiffuseKalmanSmootherH3_Z.m
parent
12c4e03d7b
commit
9304b535f4
|
@ -319,7 +319,7 @@ if kalman_algo == 2 || kalman_algo == 4
|
||||||
a_initial = zeros(np,1);
|
a_initial = zeros(np,1);
|
||||||
a_initial=set_Kalman_smoother_starting_values(a_initial,M_,oo_,options_);
|
a_initial=set_Kalman_smoother_starting_values(a_initial,M_,oo_,options_);
|
||||||
a_initial=ST*a_initial; %set state prediction for first Kalman step;
|
a_initial=ST*a_initial; %set state prediction for first Kalman step;
|
||||||
[alphahat,epsilonhat,etahat,ahat,P,aK,PK,decomp,state_uncertainty, aahat, eehat, d, regimes_,TT,RR,CC] = missing_DiffuseKalmanSmootherH3_Z(a_initial,ST, ...
|
[alphahat,epsilonhat,etahat,ahat,P,aK,PK,decomp,state_uncertainty, aahat, eehat, d, regimes_,TT,RR,CC,TTx,RRx,CCx] = missing_DiffuseKalmanSmootherH3_Z(a_initial,ST, ...
|
||||||
Z,R1,Q,diag(H), ...
|
Z,R1,Q,diag(H), ...
|
||||||
Pinf,Pstar,data1,vobs,np,smpl,data_index, ...
|
Pinf,Pstar,data1,vobs,np,smpl,data_index, ...
|
||||||
options_.nk,kalman_tol,diffuse_kalman_tol, ...
|
options_.nk,kalman_tol,diffuse_kalman_tol, ...
|
||||||
|
@ -386,9 +386,11 @@ else
|
||||||
else
|
else
|
||||||
opts_simul = options_.occbin.simul;
|
opts_simul = options_.occbin.simul;
|
||||||
end
|
end
|
||||||
|
% reconstruct smoothed variables
|
||||||
aaa=zeros(M_.endo_nbr,gend);
|
aaa=zeros(M_.endo_nbr,gend);
|
||||||
aaa(oo_.dr.restrict_var_list,:)=alphahat;
|
aaa(oo_.dr.restrict_var_list,:)=alphahat;
|
||||||
for k=2:gend
|
iTx = zeros(size(TTx));
|
||||||
|
for k=1:gend
|
||||||
if isoccbin
|
if isoccbin
|
||||||
A = TT(:,:,k);
|
A = TT(:,:,k);
|
||||||
B = RR(:,:,k);
|
B = RR(:,:,k);
|
||||||
|
@ -396,20 +398,53 @@ else
|
||||||
else
|
else
|
||||||
C=0;
|
C=0;
|
||||||
end
|
end
|
||||||
aaa(:,k) = C+A*aaa(:,k-1)+B*etahat(:,k);
|
iT = pinv(TTx(:,:,k));
|
||||||
|
% store pinv
|
||||||
|
iTx(:,:,k) = iT;
|
||||||
|
Tstar = A(~ismember(1:M_.endo_nbr,oo_.dr.restrict_var_list),oo_.dr.restrict_var_list);
|
||||||
|
Rstar = B(~ismember(1:M_.endo_nbr,oo_.dr.restrict_var_list),:);
|
||||||
|
Cstar = C(~ismember(1:M_.endo_nbr,oo_.dr.restrict_var_list));
|
||||||
|
AS = Tstar*iT;
|
||||||
|
BS = Rstar-AS*RRx(:,:,k);
|
||||||
|
CS = Cstar-AS*CCx(:,k);
|
||||||
|
static_var_list = ~ismember(1:M_.endo_nbr,oo_.dr.restrict_var_list);
|
||||||
|
ilagged = any(abs(AS*TTx(:,:,k)-Tstar)'>1.e-12);
|
||||||
|
static_var_list0 = static_var_list;
|
||||||
|
static_var_list0(static_var_list) = ilagged;
|
||||||
|
static_var_list(static_var_list) = ~ilagged;
|
||||||
|
aaa(static_var_list,k) = AS(~ilagged,:)*alphahat(:,k)+BS(~ilagged,:)*etahat(:,k)+CS(~ilagged);
|
||||||
|
if any(ilagged) && k>1
|
||||||
|
aaa(static_var_list0,k) = Tstar(ilagged,:)*alphahat(:,k-1)+Rstar(ilagged,:)*etahat(:,k)+Cstar(ilagged);
|
||||||
|
end
|
||||||
|
|
||||||
end
|
end
|
||||||
alphahat=aaa;
|
alphahat=aaa;
|
||||||
aaa=zeros(M_.endo_nbr,gend);
|
|
||||||
|
% reconstruct updated variables
|
||||||
bbb=zeros(M_.endo_nbr,gend);
|
bbb=zeros(M_.endo_nbr,gend);
|
||||||
bbb(oo_.dr.restrict_var_list,:)=ahat;
|
bbb(oo_.dr.restrict_var_list,:)=ahat; % this is t|t
|
||||||
aaa(oo_.dr.restrict_var_list,:)=aahat;
|
for k=1:gend
|
||||||
for k=d+2:gend
|
|
||||||
if isoccbin
|
if isoccbin
|
||||||
A = TT(:,:,k);
|
A = TT(:,:,k);
|
||||||
B = RR(:,:,k);
|
B = RR(:,:,k);
|
||||||
C = CC(:,k);
|
C = CC(:,k);
|
||||||
bbb(:,k) = C+A*aaa(:,k-1)+B*eehat(:,k);
|
iT = iTx(:,:,k);
|
||||||
else
|
Tstar = A(~ismember(1:M_.endo_nbr,oo_.dr.restrict_var_list),oo_.dr.restrict_var_list);
|
||||||
|
Rstar = B(~ismember(1:M_.endo_nbr,oo_.dr.restrict_var_list),:);
|
||||||
|
Cstar = C(~ismember(1:M_.endo_nbr,oo_.dr.restrict_var_list));
|
||||||
|
AS = Tstar*iT;
|
||||||
|
BS = Rstar-AS*RRx(:,:,k);
|
||||||
|
CS = Cstar-AS*CCx(:,k);
|
||||||
|
static_var_list = ~ismember(1:M_.endo_nbr,oo_.dr.restrict_var_list);
|
||||||
|
ilagged = any(abs(AS*TTx(:,:,k)-Tstar)'>1.e-12);
|
||||||
|
static_var_list0 = static_var_list;
|
||||||
|
static_var_list0(static_var_list) = ilagged;
|
||||||
|
static_var_list(static_var_list) = ~ilagged;
|
||||||
|
bbb(static_var_list,k) = AS(~ilagged,:)*ahat(:,k)+BS(~ilagged,:)*eehat(:,k)+CS(~ilagged);
|
||||||
|
if any(ilagged) && k>d+1
|
||||||
|
bbb(static_var_list0,k) = Tstar(ilagged,:)*aahat(:,k-1)+Rstar(ilagged,:)*eehat(:,k)+Cstar(ilagged);
|
||||||
|
end
|
||||||
|
elseif k>d+1
|
||||||
opts_simul.curb_retrench = options_.occbin.smoother.curb_retrench;
|
opts_simul.curb_retrench = options_.occbin.smoother.curb_retrench;
|
||||||
opts_simul.waitbar = options_.occbin.smoother.waitbar;
|
opts_simul.waitbar = options_.occbin.smoother.waitbar;
|
||||||
opts_simul.maxit = options_.occbin.smoother.maxit;
|
opts_simul.maxit = options_.occbin.smoother.maxit;
|
||||||
|
@ -433,8 +468,6 @@ else
|
||||||
bbb(oo_.dr.inv_order_var,k) = out.piecewise(1,:) - out.ys';
|
bbb(oo_.dr.inv_order_var,k) = out.piecewise(1,:) - out.ys';
|
||||||
end
|
end
|
||||||
end
|
end
|
||||||
% do not overwrite accurate computations using reduced st. space
|
|
||||||
bbb(oo_.dr.restrict_var_list,:)=ahat;
|
|
||||||
ahat0=ahat;
|
ahat0=ahat;
|
||||||
ahat=bbb;
|
ahat=bbb;
|
||||||
if ~isempty(P)
|
if ~isempty(P)
|
||||||
|
|
|
@ -486,6 +486,9 @@ varargout{1} = regimes_;
|
||||||
varargout{2} = TTT;
|
varargout{2} = TTT;
|
||||||
varargout{3} = RRR;
|
varargout{3} = RRR;
|
||||||
varargout{4} = CCC;
|
varargout{4} = CCC;
|
||||||
|
varargout{5} = TT;
|
||||||
|
varargout{6} = RR;
|
||||||
|
varargout{7} = CC;
|
||||||
% $$$ P_s=tril(P(:,:,t))+tril(P(:,:,t),-1)';
|
% $$$ P_s=tril(P(:,:,t))+tril(P(:,:,t),-1)';
|
||||||
% $$$ P1_s=tril(P1(:,:,t))+tril(P1(:,:,t),-1)';
|
% $$$ P1_s=tril(P1(:,:,t))+tril(P1(:,:,t),-1)';
|
||||||
% $$$ Fi_s = Fi(:,t);
|
% $$$ Fi_s = Fi(:,t);
|
||||||
|
|
Loading…
Reference in New Issue