function [ESR, varSR] = SRpredt(X, tau, Phi, Xbar, gamma, Omega, sigtau2, delta0, delta1, h)
% forecast of shadow rate
% Input:
%         X:       ordinary factors                      3xT
%         tau:     hidden factor                         1xT
%         Phi:     stationary lading matrix for Xtilde   3x3
%         Xbar:    constant term of X                    3x1
%         gamma    loading for tau on X                  3x1
%         Omega    covariance matrix for X               3x3
%         sigtau2  variacne of innovation for tau        scalar
%         delta0                                         scalar
%         delta1                                         3x1
%         h        forecasting horizon                   scalar
%
% Output: 
%         ESR:    expexted shadow rate      1xT
%         varSR:  variance of shadow rate   scalar
%
% Physical process of state variables 
% Z_t = [X_t', \tau_t]'
% Z_{t+1} = \mu_Z + \Phi_Z * Z_t + v_{t+1}
% Z_{t+h} = (eye(3) + \Phi_Z + \Phi_Z^2 +...+ \Phi_Z^(h-1))*\mu_Z +
%           \Phi_Z^h * Z_t +
%           v_{t+h} + \Phi_Z*v_{t+h-1} + \Phi_Z^2*v_{t+h-2} +...+ \Phi_Z^(h-1)*v_{t+1}
%
% The conditional expectation and variance of Z_{t+h}
% E_t(Z_{t+h}) = (eye(3) + \Phi_Z + \Phi_Z^2 +...+ \Phi_Z^(h-1))*\mu_Z +
%               \Phi_Z^h * Z_t
% Var_t(Z_{t+h}) = \Omega_v + \Phi_Z*\Omega_v*\Phi_Z' + \Phi_Z^2*\Omega_v*\Phi_Z^2' + ....  
%
% \mu_Z = [(eye(3) - \Phi)*Xbar; 0]
% \Phi_Z = [\Phi (eye(3) - \Phi)*gamma; 0 0 0 1]
% v_t = [u_t', \eta_t]'
% \Omega_v = [\Omega sigtau2*gamma; sigtau2*gamma' sigtau2]
%
% s_t = \delta0 + [\delta1', 0] * Z_t
% E_t(s_{t+h}) = \delta0 + [\delta1', 0] * E_t(Z_{t+h})
% Var_t(s_{t+h}) = [\delta1', 0] * Var_t(Z_{t+h}) * [\delta1', 0]'


muZ = [(eye(3) - Phi)*Xbar; 0];
PhiZ = [Phi (eye(3)-Phi)*gamma; zeros(1,3) 1];
OmegaZ = [Omega sigtau2*gamma; sigtau2*gamma' sigtau2];


Phi_exp = eye(4);
Phi_sum = eye(4);
Omega_sum = OmegaZ;

for i = 1:h-1
    Phi_exp = Phi_exp * PhiZ;
    Phi_sum = Phi_sum + Phi_exp;
    Omega_sum = Omega_sum + Phi_exp * OmegaZ * Phi_exp' ;
end

Phi_h = Phi_exp * PhiZ;
EZ = Phi_sum * muZ + Phi_h * [X; tau];   % 4xT
ESR = delta0 + [delta1' 0] * EZ;         % 1xT
varZ = Omega_sum;
varSR = [delta1' 0] * varZ * [delta1; 0];


end
