%% Filter Series in Real Time
% Consider nowcasting the model in docid:econ_ug#bt_vkl6.
% Generate a random series of 100 observations from x_t.

T = 100;
A = 0.5;
B = 1;
C = 1;
D = 0.75;
Mdl = ssm(A,B,C,D);
rng(1); % For reproducibility
y = simulate(Mdl,T);

% Suppose the final 10 observations are in the forecast horizon.

fh = 10;
yf = y((end-fh+1):end); % Holdout sample responses
y = y(1:end-fh);        % In-sample responses

% Filter the observations through the model to obtain filtered states for each period.

[xhat,logL,output] = filter(Mdl,y);
xhatvar = output.FilteredStatesCov;

% xhat and xhatvar are 90-by-1 vectors of in-sample filtered states and corresponding variances, respectively. xhat(t) is the estimate of E\left(x_t |y_1 ,\ldotp \ldotp \ldotp ,y_t \right), and xhatvar is the estimate of \textrm{Var}\left(x_t |y_1 ,\ldotp \ldotp \ldotp ,y_t \right).
% Call filter again, but specify the real-time update option.

[xhatRT,logLRT,outputRT] = filter(Mdl,y,'RealTimeUpdate',true);
xhatRTvar = outputRT.FilteredStatesCov;

% xhatRT and xhatRTvar are scalars representing the estimate of E\left(x_{90} |y_1 ,\ldotp \ldotp \ldotp ,y_{90} \right) and its corresponding variance, respectively.
% Compare the filtered states and variances of period 90.

tol = 1e-10;
areMeansEqual = (xhat(end) - xhatRT) < tol
areVarsEqual = (xhatvar(end) - xhatRTvar) < tol
areLogLsEqual = (logL - logLRT) < tol;

%



% In the last period, the filtered states, their variances, and the loglikelihoods are equal.
% Nowcast the model into the forecast horizon by performing this procedure for each successive period:
% Set the initial state and its variance to their current filter estimates. This action changes the state-space model.
% As an observation becomes available, filter it through the model in real time.

% Initialize state and variance
state0 = xhatRT;        
var0 = xhatRTvar;

% Preallocate
xhatRTF = zeros(fh,1);
xhatRTvarF = zeros(fh,1);

for j = 1:fh
    Mdl.Mean0 = state0;
    Mdl.Cov0 = var0;
    [xhatRTF(j),~,outputRT] = filter(Mdl,yf(j),'RealTimeUpdate',true); % Alternatively, use update
    xhatRTvarF(j) = outputRT.FilteredStatesCov;
    state0 = xhatRTF(j);
    var0 = xhatRTvarF(j);
end

% Plot the data and nowcasts.

figure
plot((T-fh-20):T,[y(end-20:end); yf],'b-',(T-fh+1):T,xhatRTF,'r*-')
legend(["Data" "Nowcasts"],'Location',"best")

% Copyright 2021 The MathWorks, Inc.